C test program for 'FORTRAN' version of height adjustment component 
C of air/sea boundary layer model 


      IMPLICIT NONE

C  feeds input into the sea state / flux subroutine   
      external ht_adj
      integer astab, air_moist_prm, count, dyn_in_prm, eqv_neut, 
     +        ht_adj, i, num, sfc_moist_prm, ss_prm, warn
      real    dom_phs_spd, h_sig, lhf, q_at_z, q_star, shf, t_at_z, 
     +        tau(2), t_star, u_at_z, 
     +        u_star(2), wave_age, ww_stab, z_over_L, zo_m(2)
      real    air_moist_val, CONVECT, CONV_CRIT, dtr, dyn_in_val, 
     +        pressure, ref_ht_tq, ref_ht_wind, rel_wind_ang, 
     +        salinity, sfc_moist_val, ss_val, t_air, t_skin, 
     +        tau_x, tau_y, wave_ang, wind_ang, z_wanted

      OPEN( unit=9, file='testdata99.dat' )

C     convergence critereon (fractional change)  []  
      CONV_CRIT = 0.00005

C     convective parameter 
      CONVECT = 1.25

C     warning level (1 warnings; 0 no warnings)
      warn = 1 

C     flag for wind output: 0 = wind speed, 1 = equivalent neutral wind speed
      eqv_neut = 0

C     Height for which winds, potential temperature, and humidity are adjusted.
      z_wanted = 15.0

C     conversion from degrees to radians
      dtr = 3.14159 / 180.0

C     begin loop for series of test_data 
C     skip the header line 
      READ( 9, 10)

   10 FORMAT( ) 
   11 FORMAT( A,A,A ) 
   12 FORMAT( i2, 1x, F5.2, 1x, F6.3, 1x, F6.3, 1x, F6.3, 1x, F7.4, 
     + 1x, F9.6, 1x, F8.5, 1x, F6.2, 1x, F6.2, 1x, F5.2, 1x, F6.3, 1x, 
     + F6.3, 1x, F6.2, 1x, F6.2, 1x, F6.2, 1x, F6.2, 1x, F6.3 )  
 
      WRITE(*,11) 'run  U   |ustar| ustar1 ustar2 tstar   qstar     ',
     + 'zref/L     cp     wa   Hsig   tau1   tau2    shf   lhf ',
     + '   u(z)   t(z)   q(z)'

      DO 100 i=1, 62

        READ(9,*) num, dyn_in_prm, dyn_in_val, wind_ang, wave_ang, 
     +   ss_prm, ss_val, air_moist_prm, air_moist_val, sfc_moist_prm, 
     +   sfc_moist_val, t_skin, t_air, ref_ht_wind, ref_ht_tq, pressure, 
     +   salinity, CONVECT, astab
C23456789012345678901234567890123456789012345678901234567890123456789012

C       convert angle to the 'mathimatical' coordinate system
C       conversion from meteorological direction convention
        wind_ang = 270.0 - wind_ang
C       conversion from oceanographic direction convention
        wave_ang = 90.0 - wave_ang

        rel_wind_ang = wind_ang - wave_ang

        count = ht_adj( %val(dyn_in_prm), %val(dyn_in_val),
     +   %val(rel_wind_ang), %val(CONVECT), %val(CONV_CRIT),
     +   %val(pressure), %val(air_moist_prm), %val(air_moist_val),
     +   %val(sfc_moist_prm), %val(sfc_moist_val), %val(salinity),
     +   %val(ss_prm), %val(ss_val), %val(t_air), %val(t_skin),
     +   %val(ref_ht_wind), %val(ref_ht_tq), %val(astab), %val(warn),
     +   shf, lhf, tau, u_star, t_star, q_star,
     +   z_over_L, wave_age, dom_phs_spd, h_sig, ww_stab, zo_m, 
     +   %val(eqv_neut), %val(z_wanted), u_at_z, t_at_z, q_at_z )

C       calculate meridional and zonal compoents of stress
        tau_x = tau(1) * cos(-wave_ang*DTR) - tau(2) *
     +   sin(-wave_ang*DTR)
        tau_y = tau(1) * sin(-wave_ang*DTR) + tau(2) *
     +   cos(-wave_ang*DTR)

C        check for convergence
         if ( count .EQ. 1 ) print*, 'non-convergence '  
         WRITE(*,12) i, dyn_in_val, sqrt( u_star(1) * u_star(1) + 
     +    u_star(2) * u_star(2) ), u_star(1), u_star(2), t_star, 
     +    q_star, z_over_L, dom_phs_spd, wave_age, h_sig, tau(1), 
     +    tau(2), shf, lhf, u_at_z, t_at_z, q_at_z 

  100 CONTINUE 

C     end of test program: ftest_ht_adj99.f
      END