; number of test cases:
  n = 62

  lhf = 0.0
  shf = 0.0 
  tau = [0.0, 0.0]
  u_star = [0.0, 0.0]
  t_star = 0.0
  q_star = 0.0 
  z_over_L = 0.0 
  wave_age = 0.0 
  dom_phs_spd = 0.0
  h_sig = 0.0
  ww_stab = 0.0
  zo_m = [0.0, 0.0]
  u_at_z = 0.0
  t_at_z = 0.0
  q_at_z = 0.0

  CONV_CRIT = 0.00005   ;  convergence critereon (fractional change)  []  */
  CONVECT = 1.25        ;  convective parameter  */
  warn = 1L             ;  warning are given     */
  DTR = 3.14159 / 180.0 ;  conversion from degrees to radians
  eqv_neut = 0L         ;  output winds are winds rather than equivalent neutral winds  */
  z_wanted = 15.0       ;  height to which winds, potential temp,

; open file of input data, and read first line (the header)
  OPENR, test_data, "testdata99.dat", /GET_LUN
  junk = ' '
  READF, test_data, junk  

; write header for output
  print, 'run  U   |ustar| ustar1 ustar2 tstar   qstar     ' + $
   'zref/L     cp     wa   Hsig   tau1   tau2    shf   lhf' + $
   '    u(z)   t(z)   q(z)'

  FOR i = 1, n DO BEGIN 

; read in test data
    READF, test_data, 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, $
     press_sfc, salinity, CONVECT, astab   

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

    rel_wind_ang = wind_ang - wave_ang

    dyn_in_prm = LONG( dyn_in_prm )
    ss_prm = LONG( ss_prm )
    air_moist_prm = LONG( air_moist_prm )
    sfc_moist_prm = LONG( sfc_moist_prm )
    astab = LONG( astab )
 
;    print, num, dyn_in_prm, dyn_in_val, wind_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, press_sfc, salinity, $
;     CONVECT, astab

; run the height adjustment code 
    result = ht_adj99( dyn_in_prm, dyn_in_val, rel_wind_ang, CONVECT, conv_crit, $
     press_sfc, air_moist_prm, air_moist_val, sfc_moist_prm, sfc_moist_val, $
     salinity, ss_prm, ss_val, t_air, t_skin, ref_ht_wind, ref_ht_tq, $
     astab, warn, shf, lhf, tau, u_star, t_star, q_star, z_over_L, $ 
     wave_age, dom_phs_spd, h_sig, ww_stab, zo_m, eqv_neut, z_wanted, u_at_z, $
     t_at_z, q_at_z )

;   determine zonal and meridional stresses 
    tau_x = tau(0) * cos(-wave_ang*DTR) - tau(1) * sin(-wave_ang*DTR);
    tau_y = tau(0) * sin(-wave_ang*DTR) + tau(1) * cos(-wave_ang*DTR);

; write out the results 
    IF ( result LT 0 ) THEN print, 'non-convergence '

    print, 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 )', $
     i, dyn_in_val, sqrt( u_star(0) * u_star(0) + $
     u_star(1) * u_star(1) ), u_star(0), u_star(1), t_star, $
     q_star, z_over_L, dom_phs_spd, wave_age, h_sig, tau(0), $
     tau(1), shf, lhf, u_at_z, t_at_z, q_at_z 

  ENDFOR  


END