!$Author: wagner $
!$Date: 2006-08-23 21:30:40 $
!$Revision: 1.4 $
!$Source: /weru/cvs/weps/weps_ref/src/test_friction_velocity.f95,v $
!-----------------------------------------------------------------------

PROGRAM test_friction_velocity

  USE Friction_Velocity

  IMPLICIT NONE


  ! For "gopt90" getopt commandline parsing routine
  integer :: gopt90
  character(len=1024) :: optstring = 'DlvVt:Tph?'
  character(len=1024) :: optarg, argv
  integer :: argc, c
  logical :: print_flag = .false.
  integer :: result = 0

  INTEGER :: TEST_WU_conversion_functions_1
  INTEGER :: TEST_WUF_conversion_functions_1
  INTEGER :: TEST_Rrg_reduction_function_1
  INTEGER :: TEST_Rcrow_reduction_function_1
  INTEGER :: TEST_BRcd_functions_1
  INTEGER :: TEST_relative_angle_functions_1

  INTEGER :: Incomplete_tests

  REAL :: ridge_ht, ridge_sp, ridge_dir
  REAL :: surf_rr
  REAL :: crop_lai, crop_sai, crop_ht, crop_row_sp
  REAL :: res_lai, res_sai, res_ht
  REAL :: bm_lai, bm_sai, bm_ht
  REAL :: ws, ws_dir, anemht, wzo
  REAL :: br_cd


  !real :: WUF_ref
  !real :: WUF_ref_1
  !real :: WU_ref
  !real :: WU_star
!      real :: WU_star_bv (wu_star_v, brcd)
 !     real :: WUS (wu_star, wu_star_bv)
  !real :: WZo_rr
  !real :: WZo_rg
 !     real :: WZo_rgrr (wzo_rg, wzo_rr)
  !    real :: WZo_v (bht, wzo, brcd)

!      real :: BRcd (br_lai, br_sai)
 !     real :: BRsai (Rrg, Csai, bd_sai)
  !    real :: BRlai (Rcrow, Rrg, Clai, bd_lai)
   !   real :: Rcrow (crop_ht, row_sp, ridge_ht)
  !  real :: Rrg (crop_ht, ridge_ht)

  crop_lai = 0.0
  crop_sai = 0.0
  crop_ht = 0.0
  crop_row_sp = 0.0

  res_lai = 0.0
  res_sai = 0.0
  res_ht = 0.0

  bm_lai = crop_lai + res_lai
  bm_sai = crop_sai + res_lai
  bm_ht = max(crop_ht, res_ht)

  ridge_ht = 0.0
  ridge_dir = 0.0

  surf_rr = 0.0

  DO
    c = gopt90(optstring, argc, argv);  if (c == -1) EXIT
    select case (char(c))

    case ('t')
      write(*,*) 'The optarg is: ', trim(argv), ichar(argv)
      if (trim(argv) == "1") then       ! Test wind speed conversion functions
         result = TEST_WU_conversion_functions_1(print_flag)
         write(*,*) 'Test ',trim(argv), ' result is: ', result
      else if (trim(argv) == "2") then  ! Test friction velocity conversion functions
         result = TEST_WUF_conversion_functions_1(print_flag)
         write(*,*) 'Test ',trim(argv), ' result is: ', result
      else if (trim(argv) == "3") then  ! Test reduction of _effective_ leaf and stem area when crop partly sheltered in furrow
         result = TEST_Rrg_reduction_function_1(print_flag)
         write(*,*) 'Test ',trim(argv), ' result is: ', result
      else if (trim(argv) == "4") then  ! Test reduction of _effective_ leaf area when crop rows spaced > 5 crop heights
         result = TEST_Rcrow_reduction_function_1(print_flag)
         write(*,*) 'Test ',trim(argv), ' result is: ', result
      else if (trim(argv) == "5") then  ! Test BRcd functions
         result = TEST_BRcd_functions_1(print_flag)
         write(*,*) 'Test ',trim(argv), ' result is: ', result
      end if
    case ('T')   !Test code not yet completed
         result = Incomplete_tests(print_flag)
         write(*,*) 'Test ',trim(argv), ' result is: ', result
    case ('p')   !Must be selected before test(s) are specified (verbose print flag)
      print_flag = .true.
!     write(*,*) 'This is another help option'
    case ('h')
      write(*,*) 'c is : ', c, char(c)
      c = ichar('?')
      write(*,*) 'c is : ', c, char(c)
      call display_cmd_args()
      stop
    case ('?')
      write(*,*) 'option is: ', char(c)
      call display_cmd_args()
      stop
    case default
      write(*,*) 'Default option: ', c
    end select
  END DO
STOP
END PROGRAM test_friction_velocity

SUBROUTINE display_cmd_args()
  write(*,*) 'Valid command line options:'
  write(*,*) '-t1      Test wind speed (WU) functions'
  write(*,*) '-t2      Test friction velocity (WUF) functions'
  write(*,*) '-t3      Test sai/lai (Rrg) reduction function'
  write(*,*) '-t4      Test crop lai row spacing (Rcrow) reduction function'
  write(*,*) '-T       Run test code not yet completed'
  write(*,*) '-p       Enable "print" flag'
  write(*,*) '-? or -h Display this help screen'
  RETURN
END SUBROUTINE

!! Evaluates wind speed conversion functions by specifying a range of
!! wind speeds at various anemometer heights and surface conditions
!! (different aerodynamic roughness values), converting them first to a
!! wind speed at _reference_ conditions (10 meter height and 25mm
!! aerodynamic roughness), then converting back to the original height
!! and roughness conditions.  The difference in the original and twice
!! converted wind speeds are then checked for errors outside the specified
!! tolerance value.

INTEGER FUNCTION TEST_WU_conversion_functions_1(p_flag)

  USE Friction_Velocity

  IMPLICIT NONE

  LOGICAL, INTENT(IN) :: p_flag        !! print flag
  REAL, PARAMETER :: m_to_mm = 1000.0  !! m to mm conversion factor (mm/m)
  REAL :: tolerance = 2.0e-6           !! tolerance value used for pass/fail tests

  REAL :: ws = 0.0, anemht = 0.01 , wzo = 0.0
  REAL :: max_diff = 0.0


  IF (p_flag) write(*,fmt="(10(a11))") 'ws', 'anemht', 'wzo', 'WU_ref', 'WU_from_ref', 'diff', 'max_diff'

  DO WHILE (ws <= 30.0)
     anemht = 0.01
     DO WHILE (anemht <= 30.0)
        wzo = 0.0
        DO WHILE (wzo <= 100.0)
           if (anemht*m_to_mm > wzo) then  !anemht cannot be <= wzo
              max_diff = max (max_diff, abs(ws - WU_from_WUref(WU_ref(ws,anemht,wzo),anemht,wzo)))
              if (p_flag) then
                 write(6,*) ws, anemht, wzo, WU_ref(ws,anemht,wzo), &
                   WU_from_WUref(WU_ref(ws,anemht,wzo),anemht,wzo), &
                   ws - WU_from_WUref(WU_ref(ws,anemht,wzo),anemht,wzo), max_diff
              end if
           end if
           if (wzo < 1.0) then
              wzo = wzo + 0.01
           else if (wzo < 10.0) then
              wzo = wzo + 0.1
           else if (wzo < 20.0) then
              wzo = wzo + 1.0
           else
             wzo = wzo + 5
           end if
        END DO
        anemht = 0.5 + anemht
     END DO
        ws = ws + 0.5
  END DO

  ! Test the optional args for WU_ref (Not working yet - LEW)
!  ws = 15.19; anemht = 10.0; wzo = 25.0
!  write(6,*) WU_ref(ws,anemht,wzo) !, WU_ref(ws,anemht), WU_ref(ws,wzo=wzo)

  if (p_flag) write(6,*) 'max_diff: ', max_diff, 'TINY(ws): ', TINY(ws), 'test_tolerance: ', tolerance
  IF (max_diff > tolerance) THEN
     TEST_WU_conversion_functions_1 = 1    !Test fails
     if (p_flag) write(6,*) 'WU test fails'
  ELSE
     TEST_WU_conversion_functions_1 = 0    !Test succeeds
     if (p_flag) write(6,*) 'WU test succeeds'
  END IF
END FUNCTION TEST_WU_conversion_functions_1

!! Evaluates wind speed to _reference_ friction velocity functions by
!! specifying a range of wind speeds at various anemometer heights and
!! surface conditions (different aerodynamic roughness values).
!! Two methods of computing reference_ friction velocity were applied
!! and the differences compared.
!!
!! Method 1:  Converts given wind speed at specified aerodynamic roughness
!!            and anemometer height to _reference_ conditions
!!            (10 meter height and 25mm aerodynamic roughness),
!!            then computes the _reference_ friction velocity from the
!!            the just computed _reference_ wind speed.
!! Method 2:  Converts given wind speed at specified aerodynamic roughness
!!            and anemometer height directly to the _reference_ friction
!!            velocity.
!! 
!! The difference in the two methods are then checked for errors outside
!! the specified tolerance value (both methods should yield the same values).

INTEGER FUNCTION TEST_WUF_conversion_functions_1(p_flag)

  USE Friction_Velocity

  IMPLICIT NONE

  LOGICAL, INTENT(IN) :: p_flag        !! print flag
  REAL, PARAMETER :: m_to_mm = 1000.0  !! m to mm conversion factor (mm/m)
  REAL :: tolerance = 0.0              !! tolerance value used for pass/fail tests

  REAL :: ws = 0.0, anemht = 0.01 , wzo = 0.0
  REAL :: max_diff = 0.0

  IF (p_flag) write(*,fmt="(10(a11))") 'ws', 'anemht', 'wzo', 'WU_ref', 'WUF_ref_0', 'WUF_ref_1', 'diff', 'max_diff'

  DO WHILE (ws <= 30.0)
     anemht = 0.01
     DO WHILE (anemht <= 30.0)
        wzo = 0.0
        DO WHILE (wzo <= 100.0)
           if (anemht*m_to_mm > wzo) then  !anemht cannot be <= wzo
              max_diff = max (max_diff, abs(WUF_ref(WU_ref(ws,anemht,wzo)) - WUF_ref(ws,anemht,wzo)))
              if (p_flag) then
                 write(6,*) ws, anemht, wzo, WU_ref(ws,anemht,wzo), &
                   WUF_ref(WU_ref(ws,anemht,wzo)), &
                   WUF_ref(ws,anemht,wzo), &
                   WUF_ref(WU_ref(ws,anemht,wzo)) - WUF_ref(ws,anemht,wzo), max_diff
              end if
           end if
           if (wzo < 1.0) then
              wzo = wzo + 0.01
           else if (wzo < 10.0) then
              wzo = wzo + 0.1
           else if (wzo < 20.0) then
              wzo = wzo + 1.0
           else
             wzo = wzo + 5
           end if
        END DO
        anemht = 0.5 + anemht
     END DO
        ws = ws + 0.5
  END DO

  if (p_flag) write(6,*) 'max_diff: ', max_diff, 'TINY(ws): ', TINY(ws), 'test_tolerance: ', tolerance
  IF (max_diff > tolerance) THEN
     TEST_WUF_conversion_functions_1 = 1   !Test fails
     if (p_flag) write(6,*) 'WUF test fails'
  ELSE
     TEST_WUF_conversion_functions_1 = 0   !Test succeeds
     if (p_flag) write(6,*) 'WUF test succeeds'
  END IF
END FUNCTION TEST_WUF_conversion_functions_1

!! Calculate reduction of _effective_ leaf and stem area when crop partly sheltered in furrow
!! Evaluates "reduction of _effective_ leaf and stem area when crop partly sheltered in furrow" function
!! It specifically checks to see if the value is zero for crop heights less than 1/2 the ridge height
!! and the special case of when the crop height itself is zero.
!! It also checks to see that the reduction decreases (Rrg value increases) as the crop height increases
!! above 1/2 the ridge height.
INTEGER FUNCTION TEST_Rrg_reduction_function_1(p_flag)

  USE Friction_Velocity

  IMPLICIT NONE

  LOGICAL, INTENT(IN) :: p_flag        !! print flag
  REAL, PARAMETER :: m_to_mm = 1000.0  !! m to mm conversion factor (mm/m)
  REAL :: tolerance = 0.0              !! tolerance value used for pass/fail tests

  REAL :: crop_ht = 0.0, ridge_ht = 1.0
  REAL :: max_diff = 0.0, prev = 0.0
  INTEGER :: ridge_plant_flag = 0       ! no need to test when crop is planted on ridges

  if (p_flag) write(6,*) "crop_ht, ridge_ht, ridge_plant_flag, Rrg"
  DO WHILE ( crop_ht <= 50.0)
     if (p_flag) write(6,*) crop_ht, ridge_ht, ridge_plant_flag, Rrg(crop_ht, ridge_ht, ridge_plant_flag)
     if (crop_ht <= 0.5) then
       max_diff = max_diff + Rrg(crop_ht, ridge_ht, ridge_plant_flag) ! These should sum to zero
     else
       if (Rrg(crop_ht, ridge_ht, ridge_plant_flag) <= prev) then
          max_diff = max_diff + (prev - Rrg(crop_ht, ridge_ht, ridge_plant_flag))  !problem if we get here
       endif
     endif
     prev = Rrg(crop_ht, ridge_ht, ridge_plant_flag)  ! keep copy for next iteration
     if (crop_ht < 0.5) then
       crop_ht = crop_ht + 0.1
     else if (crop_ht < 0.6) then
       crop_ht = crop_ht + 0.01
     else if (crop_ht < 2.0) then
       crop_ht = crop_ht + 0.1
     else
       crop_ht = crop_ht + 1.0
     endif
  END DO

  if (p_flag) write(6,*) "crop_ht, ridge_ht, ridge_plant_flag, Rrg"
  crop_ht = 0.0; ridge_ht = 0.0
  DO WHILE (crop_ht <= 50.0)
     if (p_flag) write(6,*) crop_ht, ridge_ht, ridge_plant_flag, Rrg(crop_ht, ridge_ht, ridge_plant_flag)
     if (crop_ht == 0.0) then
       max_diff = max_diff + Rrg(crop_ht, ridge_ht, ridge_plant_flag) ! This should sum to zero
     else
       if (Rrg(crop_ht, ridge_ht, ridge_plant_flag) /= 1.0) then  !problem if we get here
          max_diff = max_diff + 1.0
       endif
     endif
     if (crop_ht < 0.5) then
       crop_ht = crop_ht + 0.1
     else if (crop_ht < 0.6) then
       crop_ht = crop_ht + 0.01
     else if (crop_ht < 2.0) then
       crop_ht = crop_ht + 0.1
     else
       crop_ht = crop_ht + 1.0
     endif
  END DO
  if (p_flag) write(6,*) 'max_diff: ', max_diff, 'test_tolerance: ', tolerance
  IF (max_diff > tolerance) THEN
     TEST_Rrg_reduction_function_1 = 1   !Test fails
     if (p_flag) write(6,*) 'Rrg test fails'
  ELSE
     TEST_Rrg_reduction_function_1 = 0   !Test succeeds
     if (p_flag) write(6,*) 'Rrg test succeeds'
  END IF
END FUNCTION TEST_Rrg_reduction_function_1

!! Test reduction of _effective_ leaf area when crop rows spaced > 5 crop heights
INTEGER FUNCTION TEST_Rcrow_reduction_function_1(p_flag)

  USE Friction_Velocity

  IMPLICIT NONE

  LOGICAL, INTENT(IN) :: p_flag        !! print flag
  REAL, PARAMETER :: m_to_mm = 1000.0  !! m to mm conversion factor (mm/m)
  REAL :: tolerance = 0.0              !! tolerance value used for pass/fail tests

  REAL :: row_sp = 0.0, crop_ht = 0.0, ridge_ht = 1.0
  INTEGER :: ridge_plant_flag = 0
  REAL :: max_diff = 0.0, prev = 0.0

  DO WHILE (ridge_plant_flag <= 1)

  DO WHILE (row_sp <= 1000.0)
   if (p_flag) write(6,*) "crop_ht, ridge_ht, row_sp, ridge_flag, Rcrow"
   crop_ht = 0.0
   DO WHILE (crop_ht <= 500.0)
     if (p_flag) write(6,*) crop_ht, ridge_ht, row_sp, ridge_plant_flag, Rcrow(crop_ht, ridge_ht, row_sp, ridge_plant_flag)

     if (((ridge_plant_flag == 0) .and. (crop_ht <= 0.5)) .or. (crop_ht <= 0.0)) then
       max_diff = max_diff + Rcrow(crop_ht, ridge_ht, row_sp, ridge_plant_flag) ! These should sum to zero
     else
       if (Rcrow(crop_ht, ridge_ht, row_sp, ridge_plant_flag) < prev) then
          max_diff = max_diff + (prev - Rcrow(crop_ht, ridge_ht, row_sp, ridge_plant_flag))  !problem if we get here
       endif
     endif
     prev = Rcrow(crop_ht, ridge_ht, row_sp, ridge_plant_flag)  ! keep copy for next iteration
     if (crop_ht < 0.5) then
       crop_ht = crop_ht + 0.1
     else if (crop_ht < 0.6) then
       crop_ht = crop_ht + 0.01
     else if (crop_ht < 2.0) then
       crop_ht = crop_ht + 0.1
     else if (crop_ht < 50.0) then
       crop_ht = crop_ht + 1.0
     else
       crop_ht = crop_ht + 10.0
     endif
   END DO
   if (row_sp < 1.0) then
     row_sp = row_sp + 0.1
   else if (row_sp < 10.0) then
     row_sp = row_sp + 1.0
   else if (row_sp < 100.0) then
     row_sp = row_sp + 10.0
   else
     row_sp = row_sp + 100.0
   endif
  END DO
   ridge_plant_flag = ridge_plant_flag + 1
   prev = 0.0  !reset value
  END DO

 
  if (p_flag) write(6,*) 'max_diff: ', max_diff, 'test_tolerance: ', tolerance
  IF (max_diff > tolerance) THEN
     TEST_Rcrow_reduction_function_1 = 1   !Test fails
     if (p_flag) write(6,*) 'Rcrow test fails'
  ELSE
     TEST_Rcrow_reduction_function_1 = 0   !Test succeeds
     if (p_flag) write(6,*) 'Rcrow test succeeds'
  END IF
END FUNCTION TEST_Rcrow_reduction_function_1

!! Test reduction of _effective_ leaf area when crop rows spaced > 5 crop heights
INTEGER FUNCTION TEST_BRcd_functions_1(p_flag)

  USE Friction_Velocity

  IMPLICIT NONE

  LOGICAL, INTENT(IN) :: p_flag        !! print flag
  REAL, PARAMETER :: m_to_mm = 1000.0  !! m to mm conversion factor (mm/m)
  REAL :: tolerance = 0.0              !! tolerance value used for pass/fail tests

  REAL :: crop_sai = 0.0, crop_lai = 0.0
  REAL :: res_sai = 0.0, res_lai = 0.0
  REAL :: tot_sai = 0.0, tot_lai = 0.0

  REAL :: row_sp = 0.0, crop_ht = 0.0, ridge_ht = 1.0
  INTEGER :: ridge_plant_flag = 0
  REAL :: max_diff = 0.0, prev = 0.0

   if (p_flag) write(6,*) "tot_lai, tot_sai, BRcd"
   tot_sai = 0.0
   DO WHILE (tot_sai <= 50.0)
     tot_lai = 0.0
     DO WHILE (tot_lai <= 50.0)
       if (p_flag) write(6,*) tot_lai, tot_sai, BRcd(tot_lai,tot_sai)
       if (tot_lai < 0.5) then
         tot_lai = tot_lai + 0.1
       else if (tot_lai < 0.6) then
         tot_lai = tot_lai + 0.01
       else if (tot_lai < 2.0) then
         tot_lai = tot_lai + 0.1
       else if (tot_lai < 50.0) then
         tot_lai = tot_lai + 1.0
       else
         tot_lai = tot_lai + 10.0
       endif
     END DO
     if (tot_sai < 0.5) then
       tot_sai = tot_sai + 0.1
     else if (tot_sai < 0.6) then
       tot_sai = tot_sai + 0.01
     else if (tot_sai < 2.0) then
       tot_sai = tot_sai + 0.1
     else if (tot_sai < 50.0) then
       tot_sai = tot_sai + 1.0
     else
       tot_sai = tot_sai + 10.0
     endif
  END DO

 
  if (p_flag) write(6,*) 'max_diff: ', max_diff, 'test_tolerance: ', tolerance
  IF (max_diff > tolerance) THEN
     TEST_BRcd_functions_1 = 1   !Test fails
     if (p_flag) write(6,*) 'BRcd test fails'
  ELSE
     TEST_BRcd_functions_1 = 0   !Test succeeds
     if (p_flag) write(6,*) 'BRcd test succeeds'
  END IF
END FUNCTION TEST_BRcd_functions_1


!! Test relative angle functions (not done yet - LEW)
INTEGER FUNCTION TEST_relative_angle_functions_1(p_flag)

  USE Friction_Velocity

  IMPLICIT NONE

  LOGICAL, INTENT(IN) :: p_flag        !! print flag
  REAL, PARAMETER :: m_to_mm = 1000.0  !! m to mm conversion factor (mm/m)
  REAL :: tolerance = 0.0              !! tolerance value used for pass/fail tests

  REAL :: ws_dir
  REAL :: ridge_dir

     ws_dir = -770.0
     DO WHILE (ws_dir < 780.0)
       write(6,*) ws_dir, mod(ws_dir, 360.0), mod(mod(ws_dir, 360.0)+360.0, 360.0), mod(mod(ws_dir, 360.0)+360.0, 180), n180_angle(ws_dir)
       ws_dir = ws_dir + 10.0
     END DO

     ws_dir = -770.0
     DO WHILE (ws_dir < 780.0)
     ridge_dir = 0.0
     DO WHILE (ridge_dir < 360.0)
       write(6,*) ws_dir, ridge_dir, n180_angle(ws_dir), n180_angle(ridge_dir), &
                  rel_angle(ws_dir, ridge_dir)
  !     if (mod(ws_dir, 360.0) < 0.0) ws_dir = mod(ws_dir, 360.0) + 360.0 
!       write(6,*) ws_dir, mod(ws_dir, 90.0), sin(ws_dir/3.14159/180.0), sin(mod(ws_dir,90)/3.14159/180.0)
       ridge_dir = ridge_dir + 10.0
     END DO
       ws_dir = ws_dir + 10.0
     END DO
   TEST_relative_angle_functions_1 = 0            !Haven't written test yet
END FUNCTION TEST_relative_angle_functions_1

INTEGER FUNCTION Incomplete_tests(p_flag)

  USE Friction_Velocity

  IMPLICIT NONE

  LOGICAL, INTENT(IN) :: p_flag        !! print flag
  REAL, PARAMETER :: m_to_mm = 1000.0  !! m to mm conversion factor (mm/m)
  REAL :: tolerance = 0.0              !! tolerance value used for pass/fail tests

  REAL :: ridge_ht, ridge_sp, ridge_dir
  REAL :: surf_rr
  REAL :: crop_lai, crop_sai, crop_ht, crop_row_sp
  REAL :: res_lai, res_sai, res_ht
  REAL :: bm_lai, bm_sai, bm_ht
  REAL :: ws, anemht, wzo
  REAL :: br_cd

!     ws = 0.0; anemht=10.0; wzo=25.0
!     write(6,fmt="(4(a10))") 'ws', 'WUF_ref', 'WU_ref'
!     DO WHILE (ws <= 30.0)
!        write(6,*) ws, WUF_ref(ws), WU_star(WUF_ref(ws,wzo))
!        ws = 0.1 + ws
!     END DO

! surf_rr = 0.0
! write(6,fmt="(4(a10))") 'RR', 'aero_RR'
! DO WHILE (surf_rr <= 100.0)
!    write(6,*) surf_rr, WZo_rr(surf_rr)
!    surf_rr = 1.0 + surf_rr
! END DO

  surf_rr = 0.0
  write(6,fmt="(10(a12))") 'RR', 'aero_RR', 'WUF_ref(10)', 'WUS_b(10)', 'WZo_rep'
  DO WHILE (surf_rr <= 100.0)
     write(6,*) surf_rr, WZo_rr(surf_rr), WUF_ref(10.0), &
     WUS_b(WUF_ref(10.0),Wzo_rr(surf_rr)), WZo_rep(WUS_b(WUF_ref(10.0),(Wzo_rr(surf_rr))), WUF_ref(10.0))
     surf_rr = 2.0 + surf_rr
  END DO

! ridge_ht = 0.0; ridge_sp = 0.0
! write(6,fmt="(4(a10))") 'ridge_ht', 'ridge_sp', 'aero_ridge'
! DO WHILE (ridge_sp <= 500.0)
!    ridge_ht = 0.0
!    DO WHILE (ridge_ht <= 500.0)
!       write(6,*) ridge_ht, ridge_sp, Wzo_rg(ridge_ht, ridge_sp)
!       ridge_ht = 5.0 + ridge_ht
!    END DO
!    ridge_sp = 5.0 + ridge_sp
! END DO

  surf_rr = 3.0; ridge_ht = 10.0; ridge_sp = 1000.0; bm_ht = 300.0
  bm_lai = 0.3; bm_sai = 0.2

  br_cd = BRcd(bm_lai, bm_sai)
  ws = 0.0
  write(6,fmt="(10(a11))") 'ws', 'br_cd', 'WUS_rr', 'WUS_rg', 'WUS_rgrr', 'WUS_v', 'WUS_bv'
  DO WHILE (ws <= 15.0)
    write(6,*) ws, br_cd, &
      WUS_b(WUF_ref(ws), WZo_rr(surf_rr)),  &
      WUS_b(WUF_ref(ws), WZo_rg(ridge_ht,ridge_sp)),        &
      WUS_b(WUF_ref(ws), WZo_rgrr(surf_rr,ridge_ht,ridge_sp)), &
      WUS_b(WUF_ref(ws), WZo_v(bm_ht,WZo_rgrr(surf_rr,ridge_ht,ridge_sp),br_cd)), &
      WUS_bv(WUF_ref(ws), bm_ht, WZo_rgrr(surf_rr,ridge_ht,ridge_sp), br_cd)
    ws = ws + 1.0
  END DO
  ws = 0.0
  write(6,fmt="(10(a11))") 'ws', 'br_cd', 'WZo_rr', 'WZo_rg', 'WZo_rgrr', 'WZo_v', 'WZo_bv', 'WZo_rep'
  DO WHILE (ws <= 15.0)
    write(6,*) ws, br_cd, &
      WZo_rr(surf_rr),  &
      WZo_rg(ridge_ht,ridge_sp),        &
      WZo_rgrr(surf_rr,ridge_ht,ridge_sp), &
      WZo_v(bm_ht,WZo_rgrr(surf_rr,ridge_ht,ridge_sp),br_cd), &
      WZo_bv(surf_rr, ridge_ht, ridge_sp, bm_ht, br_cd), &
      WZo_rep(WUS_bv(WUF_ref(ws), bm_ht, WZo_rgrr(surf_rr,ridge_ht,ridge_sp), br_cd), WUF_ref(ws))
    ws = ws + 1.0
  END DO

  Incomplete_tests = 0   !Test succeeds (for now) :-)
END FUNCTION Incomplete_tests

! ----------------------------------------------------------------------

integer function gopt90(olist,nopt,oarg)

!-----------------------------------------------------------------------
!
! A "getopt" function that doesn't use the C function call.
! Returns -1 on exit with nothing to do
! Returns ascii code of option character.
! Allow case of space, or no space between option and arg.
!
! David Noone <dcn@caltech.edu> - Mon May  6 14:32:44 PDT 2002
!
!-----------------------------------------------------------------------
  implicit none
!------------------------- Input Arguments -----------------------------
  character(len=*), intent(in)	:: olist	! list of allowable options
!---------------------- Input/Output Arguments ------------------------- 
  integer, intent(inout)	:: nopt		! index of next option
!------------------------ Output Arguments -----------------------------
  character(len=*), intent(out)	:: oarg	! argument 
!------------------------- Local Variables -----------------------------
  character(len=1)	        :: copt		! option character
  integer			:: istr		! index of option in list
  integer			:: istr2	! another index in list
  integer 			:: llst		! length of option list
  integer			:: narg		! total number of arguments
  integer 			:: iargc	! arg counter
!-----------------------------------------------------------------------
  gopt90 = -1				! assume failure
!
! Check inputs
!
  llst = len_trim(olist)
  if (llst == 0) then
    write(*,*) '(GOPT90) No allowable options: olist'
    return
  end if
!
! Check the number of arguments
!
  narg = iargc()
  if (nopt <    0) return
  if (nopt > narg) return
  if (nopt == 0) nopt = 1		! first pass
!
! Read the next argument (option)
!
  call getarg(nopt,oarg)
  if (oarg(1:1) /= '-') return		! not an option, just return
  nopt = nopt + 1
  if (len_trim(oarg) <= 1) return	! no more options flagged
!
  copt = oarg(2:2)
  if (len_trim(oarg) == 2) then
    oarg = ''				! has space
  else
    oarg = oarg(3:len_trim(oarg))	! no space, assume argument
  end if
!
! Check that it is allowable
!
  istr = index(olist,copt)
  if (istr == 0) then
     write(*,*) '(GOPT90) Illegal option: -'//copt
     return
  end if
!
! Check that the olist has only one instance of copt (input check)
!
  istr2 = index(olist(istr+1:llst),copt)
  if (istr2 /= 0) then
     write(*,*) '(GOPT90) Multiple instances of allowable option: -'//copt
     return
  end if
!
! Check if this option needs an option, then get it
!
  if (istr /= llst) then		! if last in string, then cant have a :
     if (olist(istr+1:istr+1) == ':') then
        if (nopt > narg) then
           write(*,*) '(GOPT90) Missing argument: -'//copt
           return
        end if
!
! May already have the arg, if there was no space
!
        if (len_trim(oarg) < 1) then
          call getarg(nopt,oarg)
          nopt = nopt + 1
        end if
!
     else				! if no arg, don't allow if no-space
        if (len_trim(oarg) > 0) then
          write(*,*) '(GOPT90) Argument given for non-argument option: -'//copt
          return
        end if
     end if
  end if
! 
! Success, assign to output.
!
  gopt90 = ichar(copt)
  return
end function gopt90