c.......................................................................
c PROGRAM:  bulk_v2.5
c COARE bulk flux version 2.5b
c
c  Received from David July 1996 as refluxc3 - bug removed and gravity subroutine added
c
c  This version made from Frank Bradley's MCTEX code, built from reflux3 (from David Rogers)
c  Prepared for Flux97 workshop (NCAR, Boulder, 14-16 May 1997) for final COARE archive
c  Input set up for test data file test2_5b.dat, 4 days of Moana Wave data 26-29 Nov 1992
c  prepared from Chris Fairall's hourly data file wavhr2_5.asc 31/10/96
c  Because Chris' Tsea was at only 0.05m depth, add Ts from MSP profiler at 6m to test warm layer code
c  Profiler data from Hemantha Wijesekera
c
c  READ statement and input data stripped down to essentials
c
c 1. Based on Liu et al. 1979 and Liu's original computer code.
c
c 2. Adopt blending between Kansas and Free-convection forms by Fairall,
c    Rogers and Bradley October 20 1993 (Fairall et al 1995 JGR,101,3747-3764)
c
c 3. Modified by Fairall and Bradley June 22 1994 to include
c    a) cool skin effect (needs radiative fluxes)
c    equations quoted are from Fairall et al. 1995 JGR,101,1295-1308
c    b) sensible heat flux (ocean cooling) due to rain at wet bulb temp.
c    formalism by Gosnel, Fairall and Webster 1995 JGR,100, 18437-18442)
c    c) a simplified version of Price, Weller and Pinkel (PWP) model for solar 
c    warm layer
c    This option requires a long data file in time order, with resolution at 
c    least 1 hour, including radiative fluxes
c
c 4. Version 2f includes minor modifications and corrections to warm layer and
c    cool skin effect. (2/95)
c
c 5. Version 2.5 further code modifications (changes to Tau to output in kinematic units
c    which is required by the warm layer calculation), 6% decrease of scalar transfer 
c    coefficients by revising Kristensen separation correction. See proceedings COARE Flux
c    Group Workshop, Honolulu, Aug. 1995 (Bradley and Weller, eds., UCAR, Boulder)
c                 
c    The Webb correction to the latent heat flux, the sensible heat flux due to rainfall
c    and the rain momentum flux are now output.  Also the gustiness factor, wg.
c    BUT ef_webb has NOT been added to the latent heat flux (QE)internally.
c    Similarly RF has NOT been added to QS nor Tau_r to TAU
c
c    The form of the transfer coefficients has been changed 
c    to correspond directly with the form of the equations summarized below.  
c
c    Latitude and Longitude is now a required input.
c
c    Note the wind speed should be input relative to the sea surface.  
c
c....................................................................
c    Summary and form of the basic flux equations (Fairall et al 1995 JGR,101,3747-3764)
c
c Sensible heat flux:
c
c     QS=-rhoa*Cpa*Ustar*Tstar = rhoa*Cpa*Ch*S*(Ts-Pt)
c
c Latent heat flux:
c
c     QE=-rhoa*Le*Ustar*Qstar  = rhoa*Le*Ce*S*(qs-q)
c
c Wind stress
c
c     Tau_i=-rhoa*Ustar**2  = rhoa*Cd*S*(us_i-u_i)
c
c where 
c     Ch, Ce, and Cd are the transfer coefficients for sensible
c     heat, latent heat and stress, Pt is the potential temperature,
c     q is the water vapor mixing ratio, and u_i is one of the horizontal
c     wind components relative to the fixed earth measured at some
c     atmospheric reference height, z.  S is the average value of the wind speed
c     relative to the sea surface, which includes a gustiness factor, wg, such
c     that
c
c     S=SQR(u**2 + v**2 + wg**2).
c
c     Ts is sea surface interface temperature (skin temperature), qs is the 
c     interfacial value of the saturation mixing ratio, us_i is the surface current.
c     Note that the code requires the wind components relative to the sea surface (i.e., us_i=0).
c
c
c....................................................................
c   input:                                          
c     intime (time, either yymnddhhmmss.ssss or hh.hhhh)
c     hUm (wind measurement height) m   
c     hTm (T&rh measurement height) m  
c     hUs (wind standard height) m - use this to change your mean data to different ref. Height.   
c     hTs (T&rh standard height) m   
c     ts_depth (depth of sst instrument - a positive value) m - for cool skin/warm layer correction 
c     ws (wind speed relative to the sea surface) m/s
c     sst (sea surface temp.)  deg. C
c     atb (air temperature) deg. C
c     qq (specific humidity) g/kg or qq (RH as decimal) - Note, code works internally in kg/kg!!
c     rs (shortwave radiation) W/m2
c     rl (downwelling longwave) W/m2
c     rain (average rainrate) mm/hour
c     pp (pressure) mb; use 1008mb if data unavailable
c     zi (boundary-layer depth; use 600m if data unavailable) m
c     Jcool (=1 for cool skin calculation; =0 if SST measured by IR radiometer)
c     Jwarm (=1 for warm layer calculation; =0 if SST measured by IR radiometer)
c     intime (time variable); coaretime in yymnddhhmmss.ssss - GMT
c     glat (latitude) degrees  [latitude north +ve, latitude south -ve]
c     glon (longitude) degrees [Longitude east +ve, longitude west -ve]
c....................................................................
c
c   output:
c     QH W/m**2     (Sensible heat flux) - turbulent part only
c     QE W/m**2     (Latent heat flux) - turbulent part only
c     RF W/m**2     (Rainfall heat flux)
c     TAU N /m**2   (wind stress) - turbulent part only
c     Ustar m/s     (velocity scaling parameter - friction velocity)
c     Qstar kg/kg   (humidity scaling parameter)
c     Tstar C       (temperature scaling parameter)
c     CD - drag coefficient
c     CH - transfer coefficient for heat
c     CE - transfer coefficient for moisture
c     RR - Roughness Reynolds number
c     RT - Roughness Reynolds number for temperature
c     RQ - Roughness Reynolds number for moisture
c     ZL - height/L where L is the Obukhov length
c     Z0 - roughness length
c     zot - roughness length for temperature
c     zoq - roughness length for humidity
c     dt_wrm - total warm layer temperature difference C
c     tk_pwp - thickness of warm layer
c     dter - cool skin temperature difference C 
c     T0 - skin temperature C (T0 = sst - dter + dt_wrm*ts_depth/tk_pwp)
c     wg - gustiness factor m/s
c     Taur - momentum flux due to rain N/m**2
c     Ef_webb - Webb correction to the latent heat flux W/m**2
c
      program fluxes
      real*8 hUm,hTm,hUs,hTs,ws_h,Ta_h,qq_h 
      real*8 ws,sst,atb,qq,pp,zi,rain
      real*8 QH,QE,TAU,Ustar,Qstar,Tstar
      real*8 rl,rs,RF,T0,ts_depth,wg
      real*8 CD,CE,CH,RR,RT,RQ,Zl,ZO
      real*8 zot,zoq,dt_wrm,dter
      real*8 glat,glon,Ef_webb,Tau_r,xtime
      real*8 intime,sol_time      
      real*8 time_old,qcol_ac,tau_ac,tau_old,rf_old,hf_old,ef_old
      real*8 tk_pwp,fxp
      integer hh,mn,yy,mm,dd,jamset,jump
      character*17 chtime
      COMMON /old/time_old,qcol_ac,tau_ac,tau_old,rf_old,hf_old,
     &            ef_old,jamset,jump,fxp,tk_pwp   
c
c initialize variables
      qcol_ac=0.
      tau_ac=0.
      time_old=0.
      jamset=0.
      tau_old=0.
      hf_old=0.
      ef_old=0.
      rf_old=0.
      jamset=0
      jump=0
c
c
c open input/output files and set up fixed instrument levels
c
      open(unit=3,file='d:\123\wave\test2_5b.dat')  
      open(unit=4,file='d:\123\wave\test2_5b.out')
c
      write(4,500)
500   format(1x,'index     xtime       hsb   hlb    tub    ts     HF    
     1 EF     TAU    T0    Webb  RainF rain dt_cool dt_warm tk_pwp')
c
      hum=15.                 !height of wind measurement
      htm=15.                 !height of air temp. and RH
      hUs=10.                 !10m standard levels
      hTs=10.
c      ts_depth=0.05           !Chris Fairall's floating sensor
c      ts_depth=0.45           !from IMET buoy
      ts_depth=6.0            !Hemantha's 6m data
c
c  default values for pressure and mixed layer height
      pp=1008.
      zi=600.
c
c initial guesses to the warm layer parameters, these values simulate what is expected
c in early morning:  fxp=0.5 implies a shallow heating layer to start the integration;
c tk_pwp=19.0 implies the thickness is a maximum from the day before and is not meant to
c match this timestep's fxp.
      fxp=0.5       
      tk_pwp=19.0  
c if SST sensed by IR radiometer, set jwarm and jcool to zero
      jwarm=1
      jcool=1

c loop through data
c
      index=0
      READ(3,505)ADUM            ! skip header line of test data
505   FORMAT(A1) 
      do while (.not. eof(3))    !start of loop 
      index=index+1              !count data records (hours) 
c      
c read Chris Fairall's Moana Wave data from test file, 116 lines(hours)
c hsb,hlb,tub are bulk fluxes calculated independently by Chris for comparison
c use Hemantha's MSP data at 6m for Ts to demonstrate warm layer calculation
c
700   READ(3,705,end=910)xtime,u,tt,t,q,hsb,hlb,tub, 
     2 rs,rl,rain,xlat,xlon,ts                          
705   FORMAT(f12.0,10F6.0,3F9.0)
c
c tt is Chris' floating temperature sensor at 0.05m depth
      ws=u
      sst=ts
      atb=t
      qq=q
      glat=xlat
      glon=xlon 
      intime=xtime                    
c
c check warm layer, cool skin switches
c            
        if(Jwarm.gt.0) then
          if(index.eq.1) then
            Jwarm=2
          else
            Jwarm=1
          endif
        endif
c
c  convert time to decimal hours (e.g. yymnddhhmmss.ssss -> hh.hhhh)
c  skip this part if intime is in this format already
c
      if(intime .gt. 1e9) then                    !eg 921125182900
         write(chtime(1:17),'(f17.4)') intime
         read(chtime,12) yy,mn,dd,hh,mm,ss
12       format(5i2,f7.4)
         intime=(float(hh)+float(mm)/60.+ss/3600.) !eg 18.4833
      endif
c
c and convert to local solar time in seconds
c
         sol_time=mod(glon/15+intime+24,24)*3600   !eg 17580
c
c call bulk flux routine
c                    
      call bulk_flux(sol_time,glat,hUm,hTm,hUs,hTs,
     &   ws,sst,atb,qq,ws_h,Ta_h,qq_h,
     &   rs,rl,rain,pp,zi,Jcool,Jwarm,QH,
     &   QE,RF,TAU,Ustar,Tstar,Qstar,
     &   CD,CH,CE,RR,RT,RQ,ZL,ZO,
     &   zot,zoq,dt_wrm,dter,
     &   T0,ts_depth,wg,TAU_r,EF_webb)
c
c output results
c     
      write(4,200)index,xtime,hsb,hlb,tub,ts,
     & QH,QE,TAU,T0,Ef_webb,rf,rain,dter,dt_wrm,tk_pwp
200   format(i4,f16.0,2(2f6.0,f7.3,f7.2,2x),2f6.0,f6.1,1x,3f7.2)

      write(*,*) index
c     
      enddo                !return to beginning of loop
910   continue   
c
c end of loop
c
      stop
      end
c
c*****************************************************************
      subroutine bulk_flux(sol_time,glat,hUm,hTm,hUs,hTs,                  
     & ws,sst,atb,qq,ws_h,Ta_h,qq_h,rs,rl,rainx,pp,zix,Jcoolx,Jwarmx,
     & HF,EF,RF,TAU,Ustar,Tstar,Qstar,
     & CD,CH,CE,RRx,RTx,RQx,ZLx,ZOx,zotx,zoqx,
     & dt_wrmx,dterx,T0,ts_depthx,wgx,TAU_r,Hl_webb)
c
      real*8 hUm,hTm,hUs,hTs,ws_h,Ta_h,qq_h,ts_depthx
      real*8 ws,sst,atb,qq,pp,zix,rainx,glat
      real*8 HF,EF,TAU,Ustar,Qstar,Tstar,TAU_r,Hl_webb
      real*8 rl,rs,RF,T0,sol_time,wgx,ee
      real*8 CD,CE,CH,RRx,RTx,RQx,Zlx,ZOx
      real*8 zotx,zoqx,dt_wrmx,dterx,qjoule,qr_out,dtime      
      real*8 ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,ctd1,ctd2,rich
      real*8 dqs_dt,dwat,dtmp
c
      real*8 U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,QA
      real*8 USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt
      real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa
      real*8 visw,von,fdg,DU_Wg,Wg,newtime      
      real*8 time_old,qcol_ac,tau_ac,tau_old,rf_old,hf_old,ef_old
      real*8 fxp,tk_pwp,dsea
      integer ID,jump,jamset,ihumid
      COMMON /old/time_old,qcol_ac,tau_ac,tau_old,rf_old,hf_old,
     &            ef_old,jamset,jump,fxp,tk_pwp    
      COMMON/PIN/U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,ID
      COMMON/POUT/USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt
      COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,
     &             visa,visw,von,fdg
      COMMON/wgust/DU_Wg,Wg    
c
c initialize variables that appear in COMMON
c
      Jcool=Jcoolx
      Jwarm=Jwarmx
      ZU=hUm       !height of wind measurement
      ZT=hTm       !height of temperature measurement
      ZQ=hTm       !height of water vapor measurement
      ZUs=hUs      !standard height of wind measurement
      ZTs=hTs      !standard height of temperature measurement
      ZQs=hTs      !standard height of water vapor measurement
      U=ws         !wind speed m/s
      TS=sst       !surface temp. Celsius
      T=atb        !air temp. Celsius
      P=pp         !pressure mb
      zi=zix       !atmospheric boundary layer depth
c      
      toK=273.16
      Rnl= 0.97*(5.67e-8*(TS+toK)**4-rl)    !Net longwave (up = +)
      Rns=0.945*rs                          !Net shortwave (into water)
      rain=rainx                            !rainfall
      ts_depth=ts_depthx                    !depth of sst measurement
      call gravity(glat,grav)        
c
c Warm Layer
c
      if(Jwarm.ne.0) then    
        newtime=sol_time                    !run time in secs 
        if(Jwarm.eq.2) then                    
          jump=1
          goto 16                                 !set time_old and pass thru' ASL
        elseif(newtime.gt.21600.and.jump.eq.1) then
          goto 16                                 !6 am too late to start
        elseif(newtime.lt.time_old) then          !reset all var. at midnight
          jump=0                                  !test threshold q morning only
          jamset=0
          fxp=0.5
          tk_pwp=19.0
          tau_ac=0.0
          qcol_ac=0.0
          dt_wrm=0.0  
          goto 16
        else
          rich=.65                                    !critical Rich. No.
          ctd1=sqrt(2*rich*cpw/(al*grav*rhow))        !u*^2 integrated so
          ctd2=sqrt(2*al*grav/(rich*rhow))/(cpw**1.5) !has /rhow in both
          dtime=newtime-time_old                      !delta time
          qr_out=rnl+hf_old+ef_old+rf                 !flux out from previous pass
          q_pwp=fxp*rns-qr_out                        !effective net warming
          if(q_pwp.lt.50.and.jamset.eq.0) go to 16    !integration threshold
          jamset=1
          tau_ac=tau_ac+tau_old*dtime                 !tau from previous pass
c
          if(qcol_ac+q_pwp*dtime.gt.0) then
            do 10 index=1,5                           !iterate for warm layer thickness
              fxp=1.-(0.28*0.014*(1-dexp(-tk_pwp/0.014))   !
     &            +0.27*0.357*(1-dexp(-tk_pwp/0.357))         !
     &            +.45*12.82*(1-dexp(-tk_pwp/12.82)))/tk_pwp  !solar absorb. prof
              qjoule=(fxp*rns-qr_out)*dtime
              if((qcol_ac+qjoule.gt.0.0))
     &         tk_pwp=min(19,ctd1*tau_ac/sqrt(qcol_ac+qjoule))
   10       continue
          else
            fxp=.76
            tk_pwp=19
            qjoule=(fxp*rns-qr_out)*dtime  
          endif
c
          qcol_ac=qcol_ac+qjoule               !integrate heat input
          if(qcol_ac.gt.0) then
            dt_wrm=ctd2*(qcol_ac)**1.5/tau_ac  !pwp model warming
          else
            dt_wrm=0.
          endif         
        endif
        if(tk_pwp.lt.ts_depth) then            !sensor deeper than pwp layer
          dsea=dt_wrm                          !all warming must be added to ts
        else                                   !warming deeper than sensor
          dsea=dt_wrm*ts_depth/tk_pwp          !assume linear temperature profile
        endif
        ts=ts+dsea                             !add warming above sensor for new ts
   16   time_old=newtime
      endif
c
c end of warm layer
c
   15 call humidity(T,P,QA)         !Teten's formula returns sat. air in mb
      if(qq.lt.2.) then             !checks whether humidity in g/Kg or RH      
         R=qq
         ee=QA*R                    !convert from RH using vapour pressure      
         Q=.62197*(ee/(P-0.378*ee)) !Spec. humidity kg/kg
      else
         Q=qq/1000.                 !g/kg to kg/kg
      endif
      QA=.62197*(QA/(P-0.378*QA))   !convert from mb to spec. humidity  kg/kg
      call humidity(TS,P,QS)        !sea QS returned in mb      
      QS=QS*0.98                    !reduced for salinity Kraus 1972 p. 46
      QS=.62197*(QS/(P-0.378*QS))   !convert from mb to spec. humidity  kg/kg
c
c calculate atmospheric surface layer      
c
      call ASL(Jcool,IER)
      if(IER.ge.0) then
c
c compute surface stress (TAU), sensible heat flux (HF),  
c latent heat flux (EF) & other parameters
c
      S=sqrt(u*u + wg*wg)           !velocity incl. gustiness param.
      TAU=rhoa*USR*usr*u/S          !kinematic units
      HF=-cpa*rhoa*USR*TSR
      EF=-xlv*rhoa*USR*QSR
c
c compute heat flux due to rainfall
c
      dwat=2.11e-5*((T+toK)/toK)**1.94                    !water vapour diffusivity
      dtmp=(1.+3.309e-3*T-1.44e-6*T*T)*0.02411/(rhoa*cpa) !heat diffusivity
      dqs_dt=QA*xlv/(rgas*(T+toK)**2)                     !Clausius-Clapeyron
      alfac= 1/(1+0.622*(dqs_dt*xlv*dwat)/(cpa*dtmp))     !wet bulb factor
      W=-1.61*USR*QSR-(1+1.61*Q)*USR*TSR/(T+toK)
             
      rF= rain*alfac*cpw*((TS-T)+(QS-Q)*xlv/cpa)/3600.
c
c Webb correction to latent heat flux 
c
      Hl_webb=rhoa*xlv*W*Q
c
c compute momentum flux due to rainfall
c      
      TAU_r=0.85*rain/3600*u  
c
      dterx=dter               !cool skin parameters
      tktx=tkt          
      T0=ts-dter
      tau_old=tau 
      ef_old=ef
      hf_old=hf
      dt_wrmx=dt_wrm           !warm layer parameter
c      
c compute transfer coefficients
c
         CD=(USR/S)**2
         CH=USR*TSR/(S*(T-TS+.0098*zt+dter)) !revise 2e to 2f to include '+dter'
         CE=USR*QSR/(S*(Q-QS+dqer))                                      
         Ustar=USR
         Tstar=TSR
         Qstar=QSR
         RRx=RR
         RTx=RT
         RQx=RQ
         ZLx=ZL
         ZOx=ZO
         zotx=zot
         zoqx=zoq
         wgx=wg
         ihumid=0
         if(qq .lt. 2) ihumid=1
         call h_adjust(ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,ihumid)
         ws_h=U_hs
         Ta_h=T_hs
         qq_h=Q_hs    
      else                           !input parameters out of range
         EF=-999.
         HF=-999.
         TAU=-999.
         TAUr=-999.
         EF_webb=-999.
         Ustar=-999.
         Tstar=-999.
         Qstar=-999.
         RRx=-999.
         RTx=-999.
         RQx=-999.
         ZLx=-999.
         ZOx=-999.
         ws_h=-999.
         Ta_h=-999.
         qq_h=-999.   
         wg=-999.
      endif
      return      !return to main program
      end
c
c ------------------------------------------------------------------
      subroutine ASL(Jcool,IER)
c
c TO EVALUATE SURFACE FLUXES, SURFACE ROUGHNESS AND STABILITY OF
c THE ATMOSPHERIC SURFACE LAYER FROM BULK PARAMETERS BASED ON
c LIU ET AL. (79) JAS 36 1722-1735 
c    
      real*8 U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,D,S
      real*8 USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt
      real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa
      real*8 visw,von,fdg,DU_Wg,Wg
      real*8 PTZ,PQZ,PUZ,ZTl,ZLN,ZQl,wetc,bigc,cpv,dq,dqq,dt,dtt
      integer ID
      COMMON/PIN/U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,ID
      COMMON/POUT/USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt
      COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,
     &             visa,visw,von,fdg
      COMMON/wgust/DU_Wg,Wg     

c Factors
      Beta=1.2     ! evaluated from Fairall's low windspeed turbulence data
      Von=0.4      ! von Karman's "constant"
c      fdg=1.00     ! Fairall's LKB rr to von karman adjustment
      fdg=1.00     !based on results from Flux workshop August 1995
      toK=273.16   ! Celsius to Kelvin
c      grav=9.72    ! gravity equatorial value (ref. IGPP-SIO)
 
c Air constants and coefficients
      Rgas=287.1                  !J/kg/K     gas const. dry air
      xlv=(2.501-0.00237*TS)*1e+6  !J/kg  latent heat of vaporization at TS
      Cpa=1004.67                 !J/kg/K specific heat of dry air (Businger 1982)
      Cpv=Cpa*(1+0.84*Q)          !Moist air - currently not used (Businger 1982)
      rhoa=P*100./(Rgas*(T+toK)*(1.+.61*Q)) !kg/m3  Moist air density ( " )
      visa=1.326e-5*(1+6.542e-3*T+8.301e-6*T*T-4.84e-9*T*T*T)   !m2/s
          !Kinematic viscosity of dry air - Andreas (1989) CRREL Rep. 89-11
 
c Cool skin constants
      al=2.1e-5*(ts+3.2)**0.79     !water thermal expansion coefft.
      be=0.026                     !salinity expansion coefft.
      cpw=4000.                    !J/kg/K specific heat water
      rhow=1022.                   !kg/m3  density water
      visw=1.e-6                   !m2/s kinematic viscosity water
      tcw=0.6                      !W/m/K   Thermal conductivity water
      bigc=16.*grav*cpw*(rhow*visw)**3/(tcw*tcw*rhoa*rhoa)
      wetc=.622*xlv*QS/(rgas*(TS+toK)**2) !correction for dq;slope of sat. vap.
 
c Initialise everything
      IER=0
      ZL=0.                         
      Dter=0.                        !cool skin Dt
      Dqer=0.                        !cool skin Dq

c Initial guesses
c     next line originally required for surface current, input now requires relative wind
c      US=0.                        !surface current = 0.
c
      Wg=0.5                        !Gustiness factor initial guess
      ZO=.0005                      !roughness initial guess
      tkt=.001                      !guess sublayer thickness
      DU=U                          !assumes U is measured rel. to current
      DU_Wg=(DU**2.+Wg**2.)**.5     !include gustiness in wind spd. difference
                                    !equivalent to S in definition of fluxes
      DT=T-TS+.0098*zt              !potential temperature diff        
      DQ=Q-QS
      USR=.04*DU_Wg
      TSR=.04*DT 
      QSR=.04*DQ
      if(DU_Wg.ne.0.)then
         TA=T+toK
         RI=grav*ZU*(DT+0.61*TA*DQ)/(TA*DU_Wg**2)
      else
         IER=-1
         RI=-999.
      endif
      if(RI.gt.0.25) IER=-1
      do 200 index=1,20        !iterate 20 times 
         call ZETA(T,Q,USR,TSR,QSR,ZU,ZLN)
         ZL=ZLN
         PUZ=PSI(1,ZL)
         ZTL=ZL*ZT/ZU
         ZQL=ZL*ZQ/ZU
         PTZ=PSI(2,ZTL)
         PQZ=PSI(2,ZQL)
         ZO=0.011*USR*USR/grav + 0.11*visa/USR    !after Smith 1988 
         USR=DU_Wg*von/(dlog(ZU/ZO)-PUZ)          !Gustiness effect incl.
         RR=ZO*USR/VISA         
         call LKB(RR,RT,1)      
         if(RT.eq.-999.) then
           IER = -2                         !error - return
          return
         endif                    
   21    call LKB(RR,RQ,2)
         if(RQ.eq.-999.) then
           IER = -2                         !error - return
           return
         endif
   22    zot=rt*visa/usr
         zoq=rq*visa/usr
         S = (dlog(ZT/zot)-PTZ)/(von*fdg) !coeff fdg=1.04 included following
         D = (dlog(ZQ/zoq)-PQZ)/(von*fdg) !Fairall observations during COARE. 
                                          !NOTE coeff changed to 1.
                                          !This has been adjusted to 0.94
                                          !following discussion during FLux
                                          !workshop August 1995.
         dtt=(dt+dter)
         dqq=(dq+dqer)               ! or we could calculate new sat. hum.
         tsr=dtt/S                              !! modify
         qsr=dqq/D                              !! fluxes
         TVSR=TSR*(1.+0.61*Q)+(0.61*TA*QSR)
         Bf=-grav/TA*USR*TVSR
         if(Bf.gt.0) then
           Wg=Beta*(Bf*zi)**0.333  
         else
           Wg=0.
         endif
         DU_Wg=(DU**2.+Wg**2.)**.5        !include gustiness in wind spd.
         if(Jcool.ne.0) then              !Cool skin
           hsb=-rhoa*cpa*usr*tsr
           hlb=-rhoa*xlv*usr*qsr
           qout=rnl+hsb+hlb
           dels=rns*(.137+11*tkt-6.6e-5/tkt*(1-dexp(-tkt/8.0e-4))) ! Eq.16 
                                                    ! Shortwave, correction
                                                    ! from version 2e to 2f
           qcol=qout-dels
           if(qcol.gt.0.) then
             alq=Al*qcol+be*hlb*cpw/xlv                ! Eq. 7 Buoy flux water
             xlamx=6/(1+(bigc*alq/usr**4)**.75)**.333  ! Eq 13 Saunders coeff.
             tkt=xlamx*visw/(sqrt(rhoa/rhow)*usr)      ! Eq.11 Sublayer thickness
             dter=qcol*tkt/tcw                         ! Eq.12 Cool skin
           else
             dter=0.
           endif    
           dqer=wetc*dter
         endif             !end cool skin
  200 continue             !end iterations                 
      return               !to main subroutine, bulk_flux
      end
c
c------------------------------------------------------------------
      subroutine humidity(T,P,Qsat)                                 
c
c Tetens' formula for saturation vp Buck(1981) JAM 20, 1527-1532 
c     
      real*8 T,P,Qsat
c     
      Qsat = (1.0007+3.46e-6*P)*6.1121*dexp(17.502*T/(240.97+T))
      return
      end
c
c------------------------------------------------------------------
      subroutine LKB(RR,RT,IFLAG)
c
c TO DETERMINE THE LOWER BOUNDARY VALUE RT OF THE 
c LOGARITHMIC PROFILES OF TEMPERATURE (IFLAG=1) 
c OR HUMIDITY (IFLAG=2) IN THE ATMOSPHERE FROM ROUGHNESS 
c REYNOLD NUMBER RR BETWEEN 0 AND 1000.  OUT OF RANGE
c RR INDICATED BY RT=-999. BASED ON LIU ET AL.(1979)
c JAS 36 1722-1723 
c-------------------------------------------------------------------
c Scalar RR relation from Moana Wave data.
c
c      real*8 A(9,2),B(9,2),RAN(9),RR,RT
c      integer iflag
c      DATA A/0.177,2.7e3,1.03,1.026,1.625,4.661,34.904,1667.19,5.88E5,
c     &       0.292,3.7e3,1.4,1.393,1.956,4.994,30.709,1448.68,2.98E5/
c      DATA B/0.,4.28,0,-0.599,-1.018,-1.475,-2.067,-2.907,-3.935,
c     &       0.,4.28,0,-0.528,-0.870,-1.297,-1.845,-2.682,-3.616/
c      DATA RAN/0.11,.16,1.00,3.0,10.0,30.0,100.,300.,1000./        
c-------------------------------------------------------------------
c
c Scalar RR relation from Liu et al. 
c
      real*8 a(8,2),b(8,2),ran(8),RR,RT
      integer iflag                                            
c
      data a/0.177,1.376,1.026,1.625,4.661,34.904,1667.19,5.88e5,               
     & 0.292,1.808,1.393,1.956,4.994,30.709,1448.68,2.98e5/                     
      data b/0.,0.929,-0.599,-1.018,-1.475,-2.067,-2.907,-3.935,                
     & 0.,0.826,-0.528,-0.870,-1.297,-1.845,-2.682,-3.616/                      
      data ran/0.11,0.825,3.0,10.0,30.0,100.,300.,1000./                        
c    
      I=1
      if((RR.le.0.).or.(RR.ge.1000.)) goto 90
   10 continue
      if(RR.le.RAN(I)) goto 20
      I=I+1
      goto 10
   20 RT=A(I,IFLAG)*RR**B(I,IFLAG)
      goto 99
   90 RT=-999. 
   99 return
      end
c
c-----------------------------------------------------------------
      function PSI(ID,ZL)
C
C TO EVALUATE THE STABILITY FUNCTION PSI FOR WIND SPEED (IFLAG=1)
C OR FOR TEMPERATURE AND HUMIDITY PROFILES FROM STABILITY 
C PARAMETER ZL. SEE LIU ET AL (1979).
c Modified to include convective form following Fairall (Unpublished)
C     
      real*8 Zl,chik,psik,f,chic,psic
      integer ID
c      
      if(ZL)10,20,30                                                 
   10 F=1./(1+zl*zl) 
      CHIK=(1.-16.*ZL)**0.25
      if(ID.eq.1) goto 11
      PSIK=2.*dlog((1.+CHIK*CHIK)/2.)
      goto 12
   11 PSIK=2.*dlog((1.+CHIK)/2.)+dlog((1.+CHIK*CHIK)/2.)
     & -2.*ATAN(CHIK)+2.*ATAN(1.)
   12 CHIC=(1.-12.87*ZL)**.333    !for very unstable conditions
      PSIC=1.5*dlog((CHIC*CHIC+CHIC+1.)/3.)
     &      -(3.**.5)*ATAN((2*CHIC+1.)/(3.**.5))
     &      +4.*ATAN(1.)/(3.**0.5) 
c                     
c match Kansas and free-conv. forms with weighting F
c
      PSI= F*PSIK+(1-F)*PSIC                                        
      goto 99
   20 PSI=0.
      goto 99
   30 continue
      PSI=-4.7*ZL
   99 return
      end
c
c------------------------------------------------------------
      subroutine ZETA(T,Q,USR,TSR,QSR,Z,ZL)
C
C TO EVALUATE OBUKHOVS STABILITY PARAMETER Z/L FROM AVERAGE
C TEMP T IN DEG C, AVERAGE HUMIDITY Q IN GM/GM, HEIGHT IN M,
C AND FRICTIONAL VEL,TEMP.,HUM. IN MKS UNITS
C SEE LIU ET AL. (1979)
C     
      real*8 T,Q,OB,TVSR,TV,TA,sgn
      real*8 USR,TSR,QSR,Z,ZL
      real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa
      real*8 visw,von,fdg
      COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,
     &             visa,visw,von,fdg
c
      TA=T+toK
      TV=TA*(1.+0.61*Q)
      TVSR=TSR*(1.+0.61*Q)+0.61*TA*QSR    
      sgn=sign(1.,tvsr)               !added this to avoid program
      if(abs(tvsr) .lt. 1.e-3) then   !failure when TVSR is very small
         tvsr=sgn*tvsr
      endif
      OB=TV*USR*USR/(grav*VON*TVSR) 
      ZL=Z/OB                     
c      if(ZL .gt. 1000) ZL=1000.
      goto 99
   10 ZL=0. 
   99 return
      end
       
c-------------------------------------------------------------------------
      subroutine H_ADJUST(ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,IHUMID)
c
c This subroutine adjusts the U,T,Q variables to the specified
c standard height (ZUs,ZTs,ZQs) using the loglayer profiles.
c The DELTA correction (adjustment) is relative to the surface
c measurement.             Cronin 4/13/94
c     
      real*8 ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,ZUsL,ZTsL,ZQsL
      real*8 U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,PUZs,PTZs,PQZs
      real*8 U_wg_hs,Rho_hs,Rho_avg,QA,Rho,P_hs,ee
      real*8 USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt
      real*8 al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,visa
      real*8 visw,von,fdg,DU_Wg,Wg,S,D
      integer ID,ihumid 
c      
      COMMON/PIN/U,T,Q,TS,QS,rns,rnl,ZU,ZT,ZQ,zi,P,ID
      COMMON/POUT/USR,TSR,QSR,ZO,zot,zoq,ZL,RR,RT,RQ,RI,dter,dqer,tkt
      COMMON/const/al,beta,cpa,cpw,grav,xlv,rhoa,rhow,rgas,toK,
     &             visa,visw,von,fdg
      COMMON/wgust/DU_Wg,Wg

      call ZETA(T,Q,USR,TSR,QSR,ZUs,ZUsL)
      call ZETA(T,Q,USR,TSR,QSR,ZTs,ZTsL)
      call ZETA(T,Q,USR,TSR,QSR,ZQs,ZQsL)
      PUZs= PSI(1,ZUsL)
      PTZs= PSI(2,ZTsL)
      PQZs= PSI(2,ZQsL)
 
      S = (dlog(ZTs*USR/VISA/RT)-PTZs)/(von*fdg)
      D = (dlog(ZQs*USR/VISA/RQ)-PQZs)/(von*fdg)
      T_hs =TSR*S +TS - dter -.0098*ZTs
      Q_hs =(QSR*D + QS - dqer)*1000
      U_wg_hs = USR*(dlog(ZUs/ZO) - PUZs)/0.4
      if(U_wg_hs.ge.Wg) then
         U_hs = SQRT(U_wg_hs**2 - Wg**2)
      else
         U_hs = U_wg_hs
      endif
c
      if(IHUMID.eq.1) then    ! then need to convert sp hum into rh
         Q_hs = Q_hs/1000     ! sh kg/kg
         RHO=1./(287.*(T+273.16)*(1.+.61*Q))*P*100.
         P_hs = P - (RHO*grav*(ZTs - ZT))/100 !Approx hydrost.Pressure mb
         RHO_hs=1./(287.*(T_hs+273.16)*(1.+.61*Q_hs))*P_hs*100
         RHO_avg = (RHO + RHO_hs)/2
         P_hs = P -(RHO_avg*grav*(ZTs - ZT))/100 !hydrostatic Pressure
         call humidity(T_hs,P_hs,QA)         !Teten's formula for Pvap,sat
         ee=Q_hs*P_hs/(.62197 + .378*Q_hs)   !to get vapor pressure
         Q_hs = ee/QA                        !to get relative humidity
      endif
      return
      end
      
      Subroutine gravity(lat,g)
c       calculates g as a funciton of latitude using the 1980 IUGG formula
c         
c       Bulletin Geodesique, Vol 62, No 3, 1988 (Geodesist's Handbook)
c       p 356, 1980 Gravity Formula (IUGG, H. Moritz)
c       units are in m/sec^2 and have a relative precision of 1 part
c       in 10^10 (0.1 microGal)
c       code by M. Zumberge.
c
c       check values are:
c
c        g = 9.780326772 at latitude  0.0
c        g = 9.806199203 at latitude 45.0
c        g = 9.832186368 at latitude 90.0
c
      real*8 gamma, c1, c2, c3, c4, phi, lat, g
 
      gamma = 9.7803267715
      c1 = 0.0052790414
      c2 = 0.0000232718
      c3 = 0.0000001262
      c4 = 0.0000000007
      phi = lat * 3.14159265358979 / 180.0
      g = gamma * (1.0 
     $ + c1 * ((sin(phi))**2)
     $ + c2 * ((sin(phi))**4)
     $ + c3 * ((sin(phi))**6)
     $ + c4 * ((sin(phi))**8))
c
      return
      end
 
 
 
 
