c COARE bulk flux version 2.0
c Based on Liu et al. 1979 and Liu's original computer code.
c First modified by David Rogers October 12, 1993
c Adopt blending between Kansas and Free-convection forms by Fairall,
c Rogers and Bradley October 20 1993
c Modified by Fairall and Bradley June 22 1994 to include
c 1   cool skin effect(needs radiative fluxes)
c     equations quoted are from Fairall,Bradley and Godfrey (unpub ms 1994)
c     NB if IR SST used, set Jcool=0
c 2   sensible heat flux (ocean cooling) due to rain at wet bulb temp.
c     formalism by Gosnel,Webster and Fairall (unpub ms 1994)
c 3   a simplified version of Price, Weller and Pinkel (PWP) model for solar warm layer
c     this option requires a long data file in time order, with resolution at least 1 hour
c 4   Subroutine H_ADJUST added by Meghan Cronin 6/25/94 to adjust ws,qq,Tair
c     to specified heights (e.g. 10 m). Wind speed can be adjusted to standard
c     height different than the humidity and temperature standard height. Changes
c     within the main subroutine are noted. 
c....................................................................
c   input:
c     Intime (run start time) days          SPA       
c     hUm (wind measurement height) m       MC
c     hTm (T&rh measurement height) m       MC
c     hUs (wind standard height) m       MC
c     hTs (T&rh standard height) m       MC
c     ts_depthx (depth of sst instrument) m 
c     ws (wind speed) 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) NB Code works 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
c     zi (boundary-layer depth; use 600m if data unavailable) m
c     Jcool (=1 for cool skin calculation; =0 if SST by IR radiometer)
c     Jwarm (=0 to ignore; otherwise =2 in first line of data and =1 all other lines
c
c   output:
c
c     HF W/m**2   (Sensible heat flux)
c     EF W/m**2   (Latent heat flux)
c     RF W/m**2   (Rainfall heat flux)
c     TAU m**2/s**2
c     Ustar m/s
c     Qstar kg/kg
c     Tstar C
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 - ht/L where L is the Obukhov length
c     Z0 - roughness length
c     zot - roughness length for temperature
c     zoq - roughness length for humidity
c     T0 - skin temperature C
c
      subroutine bulk_v2e(m,intime,
     & hUm,hTm,hUs,hTs,ts_depth,                       !    MC
     & ws,sst,atb,qq,ws_h,Ta_h,qq_h,          !    MC
     & 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)
      real*8 hUm(m),hTm(m),hUs(m),hTs(m),ws_h(m),Ta_h(m),qq_h(m) 
      real*8 ws(m),sst(m),atb(m),qq(m),pp(m),zi(m),rain(m)
      real*8 QH(m),QE(m),TAU(m),Ustar(m),Qstar(m),Tstar(m)
      real*8 rl(m),rs(m),RF(m),T0(m),ts_depth(m)
      real*8 CD(m),CE(m),CH(m),RR(m),RT(m),RQ(m),Zl(m),ZO(m)
      real*8 Jcool(m),Jwarm(m),zot(m),zoq(m),dt_wrm(m),dter(m)
      real*8 intime(m),jtime,qcol_ac,tau_ac
      integer jamset
      common /old/jtime,qcol_ac,tau_ac,tau_old,rf_old,hf_old,ef_old,jamset
       loc=10                     ! local offset from GMT
c
      qcol_ac=0
      tau_ac=0
      jtime=0
      jamset=0
      tau_old=0
      hf_old=0
      ef_old=0
      
      do i=1,m
      call bulk_flux(intime(i),loc,hUm(i),hTm(i),hUs(i),hTs(i),
     & ws(i),sst(i),atb(i),qq(i),ws_h(i),Ta_h(i),qq_h(i),
     & rs(i),rl(i),rain(i),pp(i),zi(i),Jcool(i),Jwarm(i),QH(i),
     & QE(i),RF(i),TAU(i),Ustar(i),Tstar(i),Qstar(i),
     & CD(i),CH(i),CE(i),RR(i),RT(i),RQ(i),ZL(i),ZO(i),
     & zot(i),zoq(i),dt_wrm(i),dter(i),
     & T0(i),ts_depth(i))
      enddo
c
      return
      end
c
      subroutine bulk_flux(intime,loc,
     & hUm,hTm,hUs,hTs,                       !    MC
     & ws,sst,atb,qq,ws_h,Ta_h,qq_h,          !    MC
     & 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)
      real*8 hUm,hTm,hUs,hTs,ws_h,Ta_h,qq_h,ts_depthx
      real*8 ws,sst,atb,qq,pp,zix,rainx
      real*8 HF,EF,TAU,Ustar,Qstar,Tstar
      real*8 rl,rs,RF,T0,intime
      real*8 CD,CE,CH,RRx,RTx,RQx,Zlx,ZOx
      real*8 Jcoolx,Jwarmx,zotx,zoqx,dt_wrmx,dterx
      real*8 jtime,qcol_ac,tau_ac
      integer jamset
      common /old/jtime,qcol_ac,tau_ac,tau_old,rf_old,hf_old,ef_old,jamset
      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,hl,rhoa,rhow,rgas,toK,
     &visa,visw,von,fdg
      common/wgust/DU_Wg,Wg     !MC
      Jcool=Jcoolx
      Jwarm=Jwarmx
c............. MC added
         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
c............ MC end
         U=ws        !wind speed m/s
         TS=sst      !surface temp. Celsius
         T=atb       !air temp. Celsius
         P=pp        !pressure mb
         zi=zix
         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
c         ts_depth=2.0           ! Franklin thermosalinograph intake
         ts_depth=ts_depthx    ! depth of sst measurement
c -------------------- correct SST with PWP model ------------
      if(jwarm.eq.0) go to 15                   ! by-pass warm layer
         newtime=mod(loc+intime*24,24)*3600           !run time in secs SPA
      if(jwarm.eq.2.) then                      ! first line of data
         jump=1                                 
         go to 16                               ! set jtime and pass thru' ASL
      end if
      if(newtime.gt.21600.and.jump.eq.1) goto 16 ! 6 am too late to start
      if(newtime.lt.jtime) then                  ! reset at midnight
         jamset=0                                ! test threshold q morning only
         fxp=.5
         tk_pwp=19
         ts_pwp=ts
         tau_ac=0
         qcol_ac=0
         dt_wrm=0.
         jump=0
         rich=.65                                  ! critical Rich. No
         ctd1=sqrt(2*rich*cpw/(al*grav*rhow))        ! Chris integrates u*^2 so
         ctd2=sqrt(2*al*grav/(rich*rhow))/(cpw**1.5) ! has /rhoa in both of these
         go to 16
      else
         dtime=newtime-jtime                       ! 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  ! threshold to start integrating
         jamset=1
         tau_ac=tau_ac+tau_old*dtime                  ! tau from previous pass
         if(qcol_ac+q_pwp*dtime.gt.0) then
           do index=1,5                           ! iterate for warm layer thickness
           fxp=1.-(0.28*0.014*(1-exp(-tk_pwp/0.014))
     &       +0.27*0.357*(1-exp(-tk_pwp/0.357))
     &       +.45*12.82*(1-exp(-tk_pwp/12.82)))/tk_pwp        ! solar absorb. prof
           qjoule=(fxp*rns-qr_out)*dtime
           if(qcol_ac+qjoule.gt.0.)             
     &     tk_pwp=min(19,ctd1*tau_ac/sqrt(qcol_ac+qjoule))
           end do
         else
           fxp=.76
           tk_pwp=19
         endif
        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               ! pwp layer deeper than sensor
        dsea=dt_wrm
      else
        dsea=dt_wrm*ts_depth/tk_pwp             ! linear temperature profile
      endif
        ts=ts+dsea        
16    jtime=newtime

c--------------------------------- end warm layer ------------------------
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 ----------------
         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
            TAU=rhoa*USR*usr*u/sqrt(u*u+wg*wg)
            HF=-cpa*rhoa*USR*TSR
            EF=-hl*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*hl/(rgas*(T+toK)**2)                     ! Clausius-Clapeyron
         alfac= 1/(1+0.622*(dqs_dt*hl*dwat)/(cpa*dtmp))     ! wet bulb factor
         RF= rain*alfac*cpw*((TS-T)+(QS-Q)*hl/cpa)/3600.
c
c -------------------------------- cool skin parameters ---------------
c
        dterx=dter
        tktx=tkt
        T0=ts-dter
        tau_old=tau
        ef_old=ef
        hf_old=hf
c-------------------------------- warm layer parameter ----------------
        dt_wrmx=dt_wrm
C
C	COMPUTE TRANSFER COEFFICIENT
C

          CD=(USR/U)**2
          CH=USR*TSR/(U*(T-TS+.0098*zt+dter))
          CE=USR*QSR/(U*(Q-QS+dqer))                                      
            Ustar=USR
            Tstar=TSR
            Qstar=QSR
            RRx=RR
            RTx=RT
            RQx=RQ
            ZLx=ZL
            ZOx=ZO
            zotx=zot
            zoqx=zoq
c........... MC added
            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    
c......... MC end
         else
c input parameters out of range
            EF=-999.
            HF=-999.
            TAU=-999.
            Ustar=-999.
            Tstar=-999.
            Qstar=-999.
            RRx=-999.
            RTx=-999.
            RQx=-999.
            ZLx=-999.
            ZOx=-999.
c......... MC added
            ws_h=-999.
            Ta_h=-999.
            qq_h=-999.
c........... MC end
         endif
      return
      end
c ---------------------------------------------------------------------
      SUBROUTINE ASL(Jcool,IER)
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,hl,rhoa,rhow,rgas,toK,
     &visa,visw,von,fdg
      common/wgust/DU_Wg,Wg      ! MC
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
c---------------------------  Factors  -------------------------------
         Beta=1.2     ! evaluated from Fairalls low windspeed turbulence data
         Von=0.4      ! von Karman's "constant"
         fdg=1.00     ! Fairall's LKB rr to von karman adjustment
         toK=273.16   ! Celsius to Kelvin
         grav=9.72    ! gravity equatorial value (ref. IGPP-SIO)
c--------------------------- Air constants ---------------------------
         Rgas=287.1                  ! J/kg/K     gas const. dry air
         hl=(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
c      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*hl*QS/(rgas*(TS+toK)**2) ! correction for dq;slope of sat. vap.
c
c---------------------------- Initialise everything  ---------------------
       IER=0
       ZL=0.                         
       Dter=0.                        ! cool skin Dt
       Dqer=0.                        ! cool skin Dq
c----------------------------  Initial guesses  ---------------------------
       US=0.                         !surface current = 0.
       Wg=0.5                        !Gustiness factor initial guess
       ZO=.0005                      !roughness initial guess
       tkt=.001                      ! guess sublayer thickness
       DU=U-US
       DU_Wg=(DU**2.+Wg**2.)**.5     !include gustiness in wind spd. difference
       DT=T-TS+.0098*zt              !potential temperature diff        
       DQ=Q-QS
       USR=.04*DU_Wg                 !
       TSR=.04*DT                    !initial guesses                        
       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
c -----------------------------  Iterate 20 times  ------------------------
      do index=1,20
       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/(ALOG(ZU/ZO)-PUZ)              !Gustiness effect incl.
       RR=ZO*USR/VISA
       CALL LKB(RR,RT,1)
       IF(RT.NE.-999.) GOTO 21
       IER = -2                                     !error - return
       RETURN
   21 CALL LKB(RR,RQ,2)
      IF(RQ.NE.-999.) GOTO 22
      IER = -2                                      !error - return
      RETURN
   22 zot=rt*visa/usr
      zoq=rq*visa/usr
      S = (ALOG(ZT/zot)-PTZ)/(von*fdg)       !coeff fdg=1.04 included following
      D = (ALOG(ZQ/zoq)-PQZ)/(von*fdg)       !Fairall observations during COARE. 
                                             !NOTE coeff changed to 1.
      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.
c --------------------------------  Cool skin  ---------------------------------
      if(Jcool.eq.0) goto 24
      hsb=-rhoa*cpa*usr*tsr
      hlb=-rhoa*hl*usr*qsr
      qout=rnl+hsb+hlb
      dels=rns*(.137+11*tkt-6.6e-5/tkt*(1-exp(tkt/8.0e-4))) ! Eq.16 Shortwave
      qcol=qout-dels
      if(qcol.gt.0.) then
        alq=Al*qcol+be*hlb*cpw/hl                     ! 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
   24 continue
c--------------------------------- End cool skin  -------------------------------
      end do
c---------------------------------  End iterations  -----------------------------
      RETURN                                         ! to bulk_flux
      END
c
c
      Subroutine humidity(T,P,Qsat)                                 
c
c     Tetens' formula for saturation vp Buck(1981) JAM 20, 1527-1532 
c
      Qsat = (1.0007+3.46e-6*P)*6.1121*exp(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     New scalar RR relation from Moana Wave data.
C
      DIMENSION A(9,2),B(9,2),RAN(9)
      DATA A/0.177,2.7e3,1.03,1.026,1.625,4.661,34.904,1667.19,5.88E5,
     10.292,3.7e3,1.4,1.393,1.956,4.994,30.709,1448.68,2.98E5/
      DATA B/0.,4.28,0,-0.599,-1.018,-1.475,-2.067,-2.907,-3.935,
     10.,4.28,0,-0.528,-0.870,-1.297,-1.845,-2.682,-3.616/
      DATA RAN/0.11,.16,1.00,3.0,10.0,30.0,100.,300.,1000./
      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
      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.*ALOG((1.+CHIK*CHIK)/2.)
        GOTO 12
   11   PSIK=2.*ALOG((1.+CHIK)/2.)+ALOG((1.+CHIK*CHIK)/2.)
     1 -2.*ATAN(CHIK)+2.*ATAN(1.)
   12   CHIC=(1.-12.87*ZL)**.333    !for very unstable conditions
        PSIC=1.5*ALOG((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
      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
      common/const/al,beta,cpa,cpw,grav,hl,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
      IF(TVSR.EQ.0.)GOTO 10
      OB=TV*USR*USR/(grav*VON*TVSR)
      ZL=Z/OB
      GOTO 99
   10 ZL=0.
   99 RETURN
      END
       
c-------------------------------------------------------------------------
c
      SUBROUTINE H_ADJUST(ZUs,ZTs,ZQs,U_hs,T_hs,Q_hs,IHUMID)
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
      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,hl,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 = (ALOG(ZTs*USR/VISA/RT)-PTZs)/(von*fdg)
      D = (ALOG(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*(ALOG(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
c Alternatively, you could add the delta correction to the top measurement.
c It shouldn't make a difference as long as log profiles are forced to go
c through both measurements.
c       CALL ZETA(T,Q,USR,TSR,QSR,ZU,ZUL)
c       CALL ZETA(T,Q,USR,TSR,QSR,ZT,ZTL)
c       CALL ZETA(T,Q,USR,TSR,QSR,ZQ,ZQL)
c       PUZ = PSI(1,ZUL)
c       PTZ = PSI(2,ZTL)
c       PQZ = PSI(2,ZQL)
c       U_wg_hs=DU_Wg + USR*(ALOG(ZUs/ZU)-(PUZs-PUZ))/0.4
c       U_hs = sqrt(U_wg_hs**2 - Wg**2)
c       T_hs=T-.0098*(ZTs-ZT)+TSR*2.2*(ALOG(ZTs/ZT)-PTZs+PTZ)*1.15
c       Q_hs=(Q + QSR*2.2*(ALOG(ZQs/ZQ) -PQZs +PQZ)*1.15)*1000
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