!$Author: wagner $
!$Date: 2006-08-30 14:41:04 $
!$Revision: 1.6 $
!$Source: /weru/cvs/weps/weps_ref/src/surf_friction_velocity.f95,v $
!-----------------------------------------------------------------------

PROGRAM surf_friction_velocity

  USE Friction_Velocity

  IMPLICIT NONE

  REAL :: ws                   !! wind speed (m/s)
  REAL :: ws_dir               !! wind direction (degrees from North)
  REAL :: anemht               !! anemometer height (m)
  REAL :: wzo                  !! aerodynamic roughness @ anemometer location (mm)

  REAL :: ridge_ht             !! oriented roughness - ridge height (mm)
  REAL :: ridge_sp             !! oriented roughness - ridge spacing (mm)
  REAL :: ridge_dir            !! oriented roughness - ridge direction (degrees from North)

  REAL :: surf_rr              !! random roughness (mm)

  REAL :: crop_lai             !! crop leaf area index (m^2/m^2)
  REAL :: crop_sai             !! crop stem area index (m^2/m^2)
  REAL :: crop_ht              !! crop height (mm)
  REAL :: crop_row_sp          !! crop row spacing (broadcast == 0.0) (mm)
  INTEGER :: ridge_plant_flag  !! crop planted in ridge (1) or furrow (0)

  REAL :: res_lai              !! residue leaf area index (m^2/m^2)
  REAL :: res_sai              !! residue stem area index (m^2/m^2)
  REAL :: res_ht               !! residue height (mm)

  REAL :: rep_bio_ht           !! representative average combined height of standing crop and residue (mm)

  REAL :: wu                   !! wind speed at "_reference_" conditions (10m ht and 25mm aerodynamic roughness)
  REAL :: wuf                  !! friction velocity for wind speed at "_reference_" conditions (10m ht and 25mm aerodynamic roughness)
  REAL :: wus                  !! friction velocity  at surface (m/s)

  REAL, PARAMETER :: PI = 3.1415927    !! PI
  REAL :: bio_drag
  REAL :: rel_ridge_sp
  INTEGER :: i = 0

  DO WHILE (i >= 0)
    i = i + 1
    read(5,*, END = 99) ws, ws_dir, anemht, wzo                &
           ,surf_rr, ridge_ht, ridge_sp, ridge_dir             &
           ,res_ht, res_lai, res_sai                           &
           ,crop_ht, crop_lai, crop_sai, crop_row_sp, ridge_plant_flag
    write(6,*) i, ws, ws_dir, anemht, wzo                      &
           ,surf_rr, ridge_ht, ridge_sp, ridge_dir             &
           ,res_ht, res_lai, res_sai                           &
           ,crop_ht, crop_lai, crop_sai, crop_row_sp, ridge_plant_flag
    write(6,*)

    ! Reference conditions are: 10m anemometer height and 25mm aerodynamic roughness
    wu  = WU_ref(ws,anemht,wzo)
    ! "Reference" friction velocity is the friction velocity at the anemometer location
    wuf  = WUF_ref(ws,anemht,wzo)
    write(6,*) i, 'Wind speed at "reference" conditions: ', wu, '(m/s)'
    write(6,*) i, '"Reference" friction velocity: ', wuf, '(m/s)'
    write(6,*)

    ! Determine the relative angle between wind direction and ridge direction
    ! and from that the "relative" ridge spacing in the wind direction
    if (rel_angle(ws_dir, ridge_dir) <= 1.0) then  ! wind parallel to ridges
        rel_ridge_sp = 0.0  ! set to zero since relative ridge spacing is infinite
    else
        rel_ridge_sp = ridge_sp/sin(rel_angle(ws_dir, ridge_dir)*PI/180.0)
    endif
    write(6,*) i, 'Aerodynamic roughness due to random roughness: ', WZo_rr(surf_rr), '(mm)'
    !write(6,*) i, ws_dir, ridge_dir, rel_angle(ws_dir, ridge_dir), ridge_sp, rel_ridge_sp
    write(6,*) i, 'Aerodynamic roughness due to oriented roughness (ridges): ', WZo_rg(ridge_ht,rel_ridge_sp), '(mm)'
    write(6,*) i, 'Aerodynamic roughness due to total roughness: ', WZo_rgrr(surf_rr,ridge_ht,rel_ridge_sp), '(mm)'
    write(6,*)

    !  if (crop_ht /= 0.0 .and. res_ht /= 0.0) then  ! compute crop/residue effects on friction velocity
    rep_bio_ht = rep_biomass_ht (crop_ht, crop_sai, res_ht, res_sai)
    bio_drag = BRcd (crop_ht, ridge_ht, crop_row_sp, ridge_plant_flag, crop_lai, res_lai, crop_sai, res_sai)
    write(6,*) i, 'Biodrag coefficient', bio_drag
    write(6,*) i, 'Aerodynamic roughness above canopy: ',WZo_v (crop_ht, WZo_rgrr(surf_rr,ridge_ht,rel_ridge_sp), bio_drag), '(mm)'
    write(6,*) i, '"Effective" aerodynamic roughness at surface below canopy: ',WZo_bv (surf_rr, ridge_ht, rel_ridge_sp, crop_ht, bio_drag), '(mm)'

    write(6,*) i, 'Friction velocity at surface below canopy: ',WUS_bv (wuf, rep_bio_ht, WZo_rgrr(surf_rr,ridge_ht,rel_ridge_sp), bio_drag), '(m/s)'
    write(6,*) i, 'Friction velocity on surface with no canopy: ',WUS_b(wuf, WZo_rgrr(surf_rr,ridge_ht,rel_ridge_sp)), '(m/s)'
    write(6,*)
    wus = min(WUS_b(wuf, WZo_rgrr(surf_rr,ridge_ht,rel_ridge_sp)), WUS_bv (wuf, rep_bio_ht, WZo_rgrr(surf_rr,ridge_ht,rel_ridge_sp), bio_drag) )
    write(6,*) i, 'Surface friction velocity value to use for wind erosion computations: ', wus, '(m/s)'
    write(6,*)
   END DO
    goto 999

 99  write(6,*) 'Reached end of data file - exiting'
999  continue


STOP
END PROGRAM surf_friction_velocity