C     Program   flux.f  FORTRAN   August 1999  R^2
C     This is the new version of Carol Anne Clayson's
C     code for TOGA-COARE surface fluxes.
C
C     Inputs:  wind speed  in m/s = wdsp , 
C              Note - wind speeds less than 0.1 m/s will bomb !
C              sea surface temperature  in deg. C = sst ,
C              surface air temperature  in deg. C = airt ,
C              specific humidity  in gm/kg = qa ,
C              shortwave ( solar ) radiation  in W/m^2 = swr ,
C              downwelling longwave radiation  in W/m^2 = lwr .
C
C     Outputs:  momentum flux  in N/m^2 = mf ,
C               sensible heat flux  in W/m^2 = shf ,
C               latent heat flux  in W/m^2 = lhf . 
C
      subroutine flux(wdsp,sst,airt,qa,swr,lwr,mf,shf,lhf)
      implicit none
      integer*2 n,iflag
      real*4 wdsp,sst,airt,qa,swr,lwr,mf,shf,lhf
      real*4 zu,zt,zq,pres,qs,uplwr,rho,cp
      real*4 zl,du,dt,dq,ustar,qstar,thetastar,zo
      real*4 puz,ztl,zql,ptz,pqz,rr,qnet,rent,st,da
      real*4 zot,zoq,skint,zln,watt,test
C
C     Data for measurement heights in meters :
C     zu is for wind, zt is for temp. , zq is for specific humidity.
      data zu /10.0/
      data zt /10.0/
      data zq /10.0/
C     Assumed air pressure at sea level in hPa :
      data pres /1008.0/
C
C     Obtain saturation specific humidity ( qs ) for input sst :
      call getqs(sst,pres,qs)
C
C     Convert input qa to pure mass ratio :
      qa = qa/1000.0
C
C     Calculate upward going long-wave radiation in W/m^2 :
C     Assumed surface emissivity = 0.9860 .
      uplwr = 0.986*5.6705E-08*(sst + 273.15)**4
C
C     Obtain moist air density ( rho ) in kg/m^3 :
      call getrho(airt,pres,qa,rho)
C
C     The specific heat of air + water vapor at constant
C     pressure ( cp ) in J/(K kg) :
      cp = 1004.0*(1.0 + 0.9*qa)
C
C     Initialize variables for the iteration loop :
C
C     The dimensionaless height ( zeta ) :
      zl = 1.0
C     Delta wind speed ( momentum ) :
      du = wdsp
C     Delta temperature :
      dt = airt - sst
C     Delta humidity :
      dq = qa - qs
C     Friction velocity and u* scaled quantities :
      ustar = 0.04 * du
      qstar = 0.003 * dq
      thetastar = 0.003 * dt
C     Water temperature for flux calculations :
      watt = sst
C
C     Main iteration loop :
      do 333 n = 1,50
C     Call to getzo - obtain roughness length for wind :
        call getzo(ustar,wdsp,qa,thetastar,airt,qstar,watt,zu,zo)
C     Call to psi - obtain stability function for wind :
        iflag = 1
        call psi(iflag,zl,puz)
        ztl = zl*zt/zu
        zql = zl*zq/zu
C     Calls to psi - obtain stability functions for temp. & humidity :
        iflag = 2
        call psi(iflag,ztl,ptz)
        call psi(iflag,zql,pqz)
        ustar = du*0.40/(log(zu/zo)-puz)
C     The following code line ( rr ) is not used for anything at present.
C     Note, kinematic viscosity of air used is 1.4E-05 m^2/s :
        rr = zo*ustar/1.40E-05
C     Calculate fluxes for renewal time estimate :
        call getflux(rho,watt,qa,ustar,thetastar,qstar,mf,shf,lhf)
        qnet = swr + lwr - uplwr - shf - lhf
C     Calculate renewal time :
        call getrent(zo,ustar,rho,cp,airt,qnet,rent)
C     Obtain the Stanton number :
        call getst(ustar,rent,st)
C     Obtain the Dalton number :
        call getda(ustar,rent,da)
C     Call to getzot - obtain roughness length for temp. :
        call getzot(ustar,zo,st,zot)
C     Call to getzoq - obtain roughness length for humidity :
        call getzoq(ustar,zo,da,zoq)
        thetastar = dt*0.40/(log(zt/zot) - ptz)
        qstar = dq*0.40/(log(zq/zoq) - pqz)
        call getflux(rho,watt,qa,ustar,thetastar,qstar,mf,shf,lhf)
C     Calculate new sea surface temperature :
        call getskin(zo,rho,ustar,sst,lwr,uplwr,shf,lhf,swr,skint)
        watt = skint
        call getqs(skint,pres,qs)
        dt = airt - skint
        dq = qa - qs
        uplwr = 0.986*5.6705E-08*(skint + 273.15)**4
        call zeta(airt,qa,ustar,thetastar,qstar,zu,zln)
C     Convergence test :
        test = abs((zl - zln)/(zl + 1.0E-08))
        if(test .lt. 1.0E-03) goto 111
        zl = zln
        if( n .eq. 50 ) then
          write(6,*) 'Main loop fails to converge'
          mf = -999.9
          shf = -999.9
          lhf = -999.9
          return
        endif
  333 continue
C
C     Compute the final fluxes :
  111  call getflux(rho,watt,qa,ustar,thetastar,qstar,mf,shf,lhf)
      return
      end
C     End of routine flux .
C
C
C
C     Program  getrho.f
C     Returns mean moist air density using temperature in C and
C     air pressure in mbars from driver.  Rho is in kg/m^3 .
C     Note - returns  rho = 1.29227  kg/m^3  for input of
C     temp = 273.15 K , qa = 0.0, and  pres = 1.01325E+05 Pa .
      subroutine getrho(temp,pres,qa,rho)
      implicit none
      real*4 temp,pres,qa,rho,atemp,apres
      real*4 airmass,avognum,k_boltz,vtemp
C     Convert to good units:
      atemp = temp + 273.15
      apres = pres * 100.0
C     Avogadro constant ( #/mole ):
      avognum = 6.0221367E+23
C     Boltzmann constant ( J/K ):
      k_boltz = 1.380658E-23
C     Molecular weight of dry air in kg/mole:
      airmass = 2.8965E-02
C     The virtual temperature is :
      vtemp = atemp*(1.0 + 0.608*qa)
C     mean density :
      rho = ( airmass/avognum )*apres/( k_boltz*vtemp )
      return
      end
C
C
C     Program getqs.f
C     Calculates the saturation specific humidity for water
C     vapor in air for a given temperature in deg. C. and
C     air pressure in mbars (hPa).
C     Uses sixth-order polynomial fit from Flatau et al. (1992).
C     Note, this code gives a value of es = 40.42 hPa at T = 29 deg. C 
C     and 1 atmosphere.
      subroutine getqs(temp,pres,qs)
      implicit none
      integer*2 i
      real*4 ac(7),temp,rm,qs
      real*4 es,ws,pres
C     Initialize coefficients:
      data ac(1) /6.1117675/
      data ac(2) /0.44398606/
      data ac(3) /1.430533e-02/
      data ac(4) /2.650272e-04/
      data ac(5) /3.02247e-06/
      data ac(6) /2.038863e-08/
      data ac(7) /6.3878097e-10/
C
C     Mass ratio of H2O to AIR :
      rm = 18.015/28.965
C     Saturation vapor pressure in hPa = es
      es = ac(1)
      do i = 2,7
        es = es + ac(i)*temp**(i-1)
      enddo
C     Correction for sea water :
      es = 0.980*es
C     Water vapor mixing ratio (saturation) :
      ws = rm*es/(pres-es)
C     Specific humidity (saturation):
      qs = rm*es/(pres-(1.0-rm)*es)
      return
      end
C
C
C     Program   getzo.f
C     Calculates the surface roughness length for
C     momentum, zo , in meters.
C     Inputs are :
C                 ustar  in m/s
C                 thetastar  in deg. C
C                 qstar  no dim.
C                 wdsp  in m/s
C                 qa  no dim.
C                 airt  in deg. C
C                 sst  in deg. C
C                 zu  in m
C     Returns  zo  in m .
C     One has a choice as to the procedure to use.
C     Right now we use procedure 5 ( wave age stuff with
C     call to getzom ) -  flag = 5 .
C
      subroutine getzo(ustar,wdsp,qa,thetastar,airt,qstar,sst,zu,zo)
      implicit none
      integer*2 flag,i,n
      real*4 ustar,wdsp,qa,thetastar,airt,qstar,sst
      real*4 zu,zo,cd,u10,c,zon,coef,test
      real*4 betag,betac,betas,wa,usfc,moninv,bstar,zom
C
C     Set flag for procedure :
      flag = 5
C
C     Procedure number 1 :
C     Uses the Charnock relation for open ocean (Smith 1988).
C
      if( flag .eq. 1 ) then
        zo = 0.011*ustar*ustar/9.81
        return
      endif
C
C     Procedure number 2 :
C     Uses windspeed to determine the drag coefficient ( cd ).
C
      if( flag .eq. 2 ) then
        if( wdsp .lt. 10.0) then
          cd = 1.14E-03
        else
          cd = 4.9E-04 * 0.065*wdsp
        endif
        zo = zu * exp(-0.40/sqrt(cd))
        return
      endif
C
C     Procedure number 3 :
C     Uses windspeed to determine the drag coefficient ( cd ).
C
      if( flag .eq. 3 ) then
        if( wdsp .lt. 2.2 ) then
          cd = 1.08E-03/wdsp**0.15
        else if( wdsp .ge. 2.2 .and. wdsp .lt. 5.0 ) then
          cd = (0.771 + 0.0858*wdsp)*1.0E-03
        else if( wdsp .ge. 5.0 .and. wdsp .lt. 8.0 ) then
          cd = (0.867 + 0.0667*wdsp)*1.0E-03
        else if( wdsp .ge. 8.0 .and. wdsp .lt. 25.0 ) then
          cd = (1.2 + 0.025*wdsp)*1.0E-03
        else
          cd = 0.073*wdsp*1.0E-03
        endif
        zo = zu * exp(-0.40/sqrt(cd))
        return
      endif
C
C     Procedure number 4 :
C     Makes calls to drag.
C
      if( flag .eq. 4 ) then
        zo = 0.0
        do i = 1,51
          u10 = ustar*log(10.0/zo)/0.40
          call drag(u10,coef)
          cd = coef
          c = 1.0/sqrt(cd)
          zon = 10.0/exp(0.4*c)
          test = abs((zon - zo)/(zo + 1.0E-08))
          if(test .lt. 0.010) goto 77
          zo = zon
          if( i .eq. 50) then
            write(6,*)'No convergence of proc = 4 in getzo'
            stop
          endif
        enddo
  77    return
      endif
C
C     Procedure number 5 :
C     Makes calls to getbstar and getzom .
C     Note that  moninv = 1/L  where L is the local Monin-Obukhov length.
C
      if( flag .eq. 5 ) then
        betag = 1.0
        betac = 1.0
        betas = 0.0
        wa = 20.0
        usfc = 0.55*ustar
        moninv = 0.0
        n = 0
        call getbstar(qa,thetastar,airt,qstar,bstar)
        moninv = 0.40*bstar/(ustar*ustar)
        if( abs(moninv) .lt. 0.001 ) then
          moninv = 0.001*moninv/abs(moninv)
        endif
        if( abs(moninv) .gt. 10.0 ) then
          moninv = 10.00*moninv/abs(moninv)
        endif
        call getzom(sst,wdsp,zu,n,betag,betac,betas,ustar,wa,
     &  usfc,moninv,zom)
        zo = zom
        return
      endif
C
      end
C
C
C     Program   drag.f
C     Calculates the  drag coefficient  in some way.
C     Called by subroutine  getzo( ) - procedure number 4 .
C     Inputs  u10  and returns  coef ( cd ) .
C     Note - this routine looks incomplete as written -
C     there may be problems when using procedure 4 of getzo .
      subroutine drag(u10,coef)
      implicit none
      integer*2 id,k
      real*4 u10,coef,cd
      real*4 r(5),a(5),b(5),p(5)
C
C     Initialize array values :
      data r(1) /2.2/
      data r(2) /5.0/
      data r(3) /8.0/
      data r(4) /25.0/
      data r(5) /50.0/
C
      data a(1) /0.0/
      data a(2) /0.771/
      data a(3) /0.867/
      data a(4) /1.2/
      data a(5) /0.0/
C
      data b(1) /1.08/
      data b(2) /0.0858/
      data b(3) /0.0667/
      data b(4) /0.025/
      data b(5) /0.073/
C
      data p(1) /-0.15/
      data p(2) /1.0/
      data p(3) /1.0/
      data p(4) /1.0/
      data p(5) /1.0/
C
C     Set-up parameter for this routine :
      id = 0
      k = id - 2
C
      if(k .lt. 0) then
        if(u10 .gt. 50.0) then
          cd = 3.70E-03
        else if(u10 .lt. 0.3) then
          cd = 1.5E-03
        else
          if(u10 .gt. r(1)) then
            cd = (a(1) + b(1) * u10**p(1))/1000.0
          endif
          if(u10 .gt. r(2)) then
            cd = (a(2) + b(2) * u10**p(2))/1000.0
          endif
          if(u10 .gt. r(3)) then
            cd = (a(3) + b(3) * u10**p(3))/1000.0
          endif
          if(u10 .gt. r(4)) then
            cd = (a(4) + b(4) * u10**p(4))/1000.0
          endif
          if(u10 .gt. r(5)) then
            cd = (a(5) + b(5) * u10**p(5))/1000.0
          endif
        endif
      endif
C
      if(k .eq. 0) then
        cd = (0.61 + 0.063*u10)/1000.0
      endif
C
      if(k .gt. 0) then
        if(u10 .ge. 11.0) then
          cd = (0.49 + 0.065*u10)/1000.0
        else
          cd = 1.2E-03
        endif
      endif
C
      coef = cd
      return
      end
C
C
C     Program   bstar.f
C     Calculates  bstar  -  note that the buoyancy
C     flux is given by  - bstar * ustar .
C     Called from  getzo( ) - procedure number 5 .
C     Inputs are :
C                 the mean specific humidity  =  hmean
C                 thetastar  in deg. C (= K)  =  tempstar
C                 mean air temperature in deg. C  =  tmean
C                 qstar ( humidity ) no dim.  =  humstar
C     Output is :  bstar
      subroutine getbstar(hmean,tempstar,tmean,humstar,bstar)
      implicit none
      real*4 hmean,tempstar,tmean,humstar
      real*4 bstar,ta
C
C     Convert temperature to K :
      ta = tmean + 273.15
C
      bstar =  9.810*((1.0 + 0.608*hmean)*tempstar +
     &     0.608*humstar*ta)/(ta*(1.0 + 0.608*hmean))
C
      if( abs(bstar) .lt. 0.000001 ) then
        if(bstar .eq. 0.00) then
          bstar = 0.0000010
        else
          bstar = 0.000001*bstar/abs(bstar)
        endif
      endif
C
      return
      end
C
C
C     Program   psi.f
C     Calculates the stability functions for
C     the wind, temperature, and specific humidity
C     profiles.
C     Input parameters :   iflag  ( an integer )
C                          zeta  the dimensionless height
C     Returns the stability function value  =  stabfunc
C
      subroutine psi(iflag,zeta,stabfunc)
      implicit none
      integer*2 iflag
      real*4 zeta,stabfunc
      real*4 pi,a,b,c,d,gamma
      real*4 blendf,chic,chik,psic,psik
      real*4 psim,psih
C
      data pi /3.1415926536/
C     Coefficient values from Beljaars & Holtslag (1991) for
C     stable situation ( zeta > 0. ) :
      data a /1.0/
      data b /0.667/
      data c /5.0/
      data d /0.35/
C     Empirical constant gamma of Fairall et al. (1996) :
      data gamma /12.87/
C
C     Note -  iflag = 1 for wind profile function  and
C             iflag = 2 for temp. and humidity functions .
C
C     For zeta < 0. use the new Fairall version of convective form :
      if( zeta .lt. 0.0 ) then
        blendf = 1.0/(1.0 + zeta*zeta)
        chik = sqrt(sqrt(1.0 - 16.0*zeta))
C
C     Equations (14) and (13) of Fairall et al. (1996) :
        chic = (1.0-gamma*zeta)**(1.0/3.0)
        psic = 1.50*log((chic*chic + chic + 1.0)/3.0)
     &   - sqrt(3.0)*atan((2.0*chic + 1.0)/sqrt(3.0))
     &   + pi/sqrt(3.0)
C
        if( iflag .eq. 1 ) then
C     Equation (25) of Beljaars & Holtslag (1991) :
          psik = 2.0*log((1.0 + chik)/2.0)
     &    + log((1.0 + chik*chik)/2.0) - 2.0*atan(chik) + pi/2.0
C     Equation (26) of Fairall et al. (1996) :
          stabfunc = blendf*psik + (1.0 - blendf)*psic
          return
        endif
        if( iflag .eq. 2 ) then
C     Equation (26) of Beljaars & Holtslag (1991) :
          psik = 2.0*log((1.0 + chik*chik)/2.0)
          stabfunc = blendf*psik + (1.0 - blendf)*psic
          return
        endif
      endif
C
C     For zeta = 0. :
      if( zeta .eq. 0.00 ) then
        stabfunc = 0.0
        return
      endif
C
C     For zeta > 0. :
      if( zeta .gt. 0.0 ) then
        if( iflag .eq. 1 ) then
C     Equation (28) of Beljaars & Holtslag (1991) :
          psim = a*zeta + b*(zeta - c/d)*exp(-d*zeta) + b*c/d
          stabfunc = -1.0*psim
          return
        endif
        if( iflag .eq. 2 ) then
C     Equation (32) of Beljaars & Holtslag (1991) :
          psih = (1.0 + 2.0*a*zeta/3.0)**1.50
     &     + b*(zeta - c/d)*exp(-d*zeta) + b*c/d - 1.0
          stabfunc = -1.0*psih
          return
        endif
      endif
C
      end
C
C
C     Program   getflux.f
C     Calculates the fluxes.
C     Inputs are :
C                 moist air density in kg/m^3  =  rho
C                 water temp. in deg. C  =  watt
C                 specific humidity  =  qa
C                 ustar
C                 thetastar
C                 qstar
C     Returns :
C              momentum flux in N/m^2  =  mf
C              sensible heat flux in W/m^2  =  shf
C              latent heat flux in W/m^2  =  lhf
C
      subroutine getflux(rho,watt,qa,ustar,thetastar,qstar,
     &  mf,shf,lhf)
      implicit none
      real*4 rho,watt,qa,ustar,thetastar,qstar
      real*4 mf,shf,lhf
      real*4 cp,lv
C
C     The specific heat of air + water vapor at constant
C     pressure ( cp ) in J/(K kg) :
      cp = 1004.0*(1.0 + 0.9*qa)
C     The latent heat of vaporization for water ( lv ) in J/kg :
      lv = 4186.8*(597.31 - 0.56525*watt)
C     The momentum flux in N/m^2 :
      mf = rho*ustar*ustar
C     The sensible heat flux in W/m^2 :
      shf = -1.0*rho*cp*ustar*thetastar
C     The latent heat flux in W/m^2 :
      lhf = -1.0*rho*lv*ustar*qstar
C
      return
      end
C
C
C     Program   renew-time.f
C     Calculates the surface renewal timescale.
C     This routine may have problems with small values of ustar  -
C     the input ustar depends on the wind speed which may be small.
C     Inputs are :
C                  zo  in m
C                  ustar  in m/s
C                  rho  in kg/m^3
C                  cp  in J/(K kg)
C                  airt  in deg. C
C                  qnet  in W/m^2
C     Returns :   rent  in sec.
C
      subroutine getrent(zo,ustar,rho,cp,airt,qnet,rent)
      implicit none
      real*4 zo,ustar,rho,cp,airt,qnet,rent
      real*4 shear,conv,rfcr,csh,ccon,rfo
      real*4 alpha,visc,args,argc
C
C     Note - there is a strange sign convention for the Richardson
C     numbers here.
C     The critical value of the surface Richardson number :
      rfcr = -2.0E-04
C     Note - csh & ccon are empirically determined constants (no dim.).
C     The following density check seems unnecessary since
C     the input rho is always < 1.5 kg/m^3 :
      if( rho .lt. 50.0 ) then
        rfcr = -2.0E-01
        csh = 53.5824
        ccon = 0.80
      else
C     From Wick thesis :
        csh = 209.0
        ccon = 3.13
      endif
C
C     The kinematic viscosity in m^2/s from Gill pg. 75 :
      visc = 1.40E-05
      alpha = 1.0/(airt + 273.15)
C     Note - the qnet used here from flux( ) does not contain a 
C     contribution from solar radiation ( swr ) .
C     Obtain the Richardson number :
      rfo = 9.810*qnet*alpha*visc/(rho*cp*ustar**4)
      if( rfo .gt. 0.0 ) then
        rfo = -1.0*rfo
      endif
C
C     Obtain the time scale :
C     Watch out for small values of ustar !
C     These can occur for small wind speeds ( < 0.1 m/s ).
      args = visc*zo/ustar**3
C     Shear-driven timescale ( sec. ) :
      shear = csh*sqrt(args)
      argc = visc*rho*cp/(alpha*9.810*abs(qnet))
C     Convective timescale ( sec. ) :
      conv = ccon*sqrt(argc)
      rent = shear + (conv - shear)*exp(-rfcr/rfo)
C
      return
      end
C
C
C     Program   st.f
C     Determine the interfacial Stanton number.
C     Inputs are ustar in m/s and the surface
C     renewal time ( rent ) in sec.
C     Returns the Stanton number ( st ) no dim.
      subroutine getst(ustar,rent,st)
      implicit none
      real*4 ustar,rent,st
      real*4 tdc
C
C     The thermal diffusion viscosity for air from Gill p. 71
C     in m^2/s :
      tdc = 2.0E-05
C
      st = sqrt(tdc/(ustar*ustar*rent))
      return
      end
C
C
C     Program   da.f
C     Determine the interfacial Dalton number.
C     Inputs are ustar in m/s and the surface
C     renewal time ( rent ) in sec.
C     Returns the Dalton number ( da ) no dim.
      subroutine getda(ustar,rent,da)
      implicit none
      real*4 ustar,rent,da
      real*4 mdc
C
C     The molecular diffusion viscosity for H2O in air from Gill p. 68
C     in m^2/s :
      mdc = 2.4E-05
C
      da = sqrt(mdc/(ustar*ustar*rent))
      return
      end
C
C
C     Program   zot.f
C     Calculates the roughness length for
C     temperature ( heat ) using equation (21)
C     of Clayson et al. (1996).
C     Inputs are :
C                  ustar  in m/s ( not used at this time )
C                  zo  in m
C                  st  no dim.
C     Returns roughness length  zot  in m .
      subroutine getzot(ustar,zo,st,zot)
      implicit none
      real*4 ustar,zo,st,zot,prt
C
C     The turbulent Prandtl number ( prt ) :
      prt = 0.85
C
      zot = zo*exp(0.40*(5.0 - 1.0/(prt*st)))
      return
      end
C
C
C     Program   zoq.f
C     Calculates the roughness length for
C     humidity ( water vapor ) using equation (20)
C     of Clayson et al. (1996).
C     Inputs are :
C                  ustar  in m/s ( not used at this time )
C                  zo  in m
C                  da  no dim.
C     Returns roughness length  zoq  in m .
      subroutine getzoq(ustar,zo,da,zoq)
      implicit none
      real*4 ustar,zo,da,zoq,sc
C
C     The turbulent Schmidt number ( sc ) :
      sc = 0.60
C
      zoq = zo*exp(0.40*(5.0 - 1.0/(sc*da)))
      return
      end
C
C
C     Program   skin.f
C     Calculates the sea surface skin temperature.
C     Inputs are :
C           zo  in m
C           rho  in kg/m^3
C           ustar  in m/s
C           sst  in deg. C
C           lwr  in W/m^2
C           uplwr  in W/m^2
C           shf  in W/m^2
C           lhf  in W/m^2
C           swr  in W/m^2
C     Returns skin temperature ( skint ) in deg. C
C     Call to subroutine getwrt( ) .
C
C     Must set water type for solar absorption.
C
      subroutine getskin(zo,rho,ustar,sst,lwr,uplwr,shf,lhf,
     &  swr,skint)
      implicit none
      integer*2 watertype
      real*4 zo,rho,ustar,sst,lwr,uplwr,shf,lhf
      real*4 swr,skint
      real*4 sal,wrho,wcp,wkappa,tec,qnet
      real*4 wustar,wrt,qsolar
      real*4 rc,gamma1,gamma2,zsm,delta
C
C     Determine the water type for solar absorption.
C     Type I = 1
C     Type IA = 11
C     Type IB = 12
C     Type II = 2
C     Type III = 3
C     no solar absorption = 0
C
      watertype = 1
C
C     Some constants :
C     Salinity in practical salinity units :
      data sal /35.0/
C     Sea water density at 30 deg. C in kg/m^3 :
      data wrho /1021.7/
C     Specific heat of sea water in J/(K kg) :
      data wcp /4002.0/
C     Kappa from Gill pg. 71 :
      data wkappa /1.4E-07/
C     Thermal expansion coefficient for sea water in 1/deg. C
C     from Gill Table A3.1 :
      data tec /3.196E-04/
C
C     Determine qsolar from water type :
C
      if( watertype .ne. 0 ) then
        if( watertype .eq. 1 ) then
          rc = 0.58
          gamma1 = 0.35
          gamma2 = 23.0
        endif
        if( watertype .eq. 11 ) then
          rc = 0.62
          gamma1 = 0.60
          gamma2 = 20.0
        endif
        if( watertype .eq. 12 ) then
          rc = 0.67
          gamma1 = 1.00
          gamma2 = 17.0
        endif
        if( watertype .eq. 2 ) then
          rc = 0.77
          gamma1 = 1.50
          gamma2 = 14.0
        endif
        if( watertype .eq. 3 ) then
          rc = 0.78
          gamma1 = 1.40
          gamma2 = 7.9
        endif
        zsm = 0.0010
        qsolar = swr - swr*(rc*exp(-zsm/gamma1) + (1.0 - rc)*
     &         exp(-zsm/gamma2))
      endif
      if( watertype .eq. 0 ) then
        qsolar = 0.0
      endif
C
C     Net heat loss from ocean :
      qnet = lwr - uplwr - shf - lhf + qsolar
C
      wustar = sqrt(rho/wrho)*ustar
C     Obtain renewal time from getwrt( ) :
      call getwrt(zo,wustar,wrho,wcp,tec,qnet,wrt)
C
C     Calculate skint :
      delta = sqrt(wrt)*qnet/(wrho*wcp*sqrt(wkappa))
      skint = delta + sst
C
      return
      end
C
C
C     Program   renew-time.f
C     Calculates the surface renewal timescale.
C     Inputs are :
C                  zo  in m
C                  ustar  in m/s
C                  rho  in kg/m^3
C                  cp  in J/(K kg)
C                  airt  in deg. C
C                  qnet  in W/m^2
C     Returns :   rent  in sec.
C
      subroutine getwrt(zo,ustar,rho,cp,alpha,qnet,rent)
      implicit none
      real*4 zo,ustar,rho,cp,qnet,rent
      real*4 shear,conv,rfcr,csh,ccon,rfo
      real*4 alpha,visc,args,argc
C
C     Note - there is a strange sign convention for the Richardson
C     numbers here.
C     The critical value of the surface Richardson number :
      rfcr = -2.0E-04
C     Note - csh & ccon are empirically determined constants (no dim.).
      if( rho .lt. 50.0 ) then
        rfcr = -2.0E-01
        csh = 53.5824
        ccon = 0.80
      else
C     From Wick thesis :
        csh = 209.0
        ccon = 3.13
      endif
C
C     The kinematic viscosity in m^2/s from Gill pg. 75 :
      visc = 1.0E-06
C     Obtain the Richardson number :
      rfo = 9.810*qnet*alpha*visc/(rho*cp*ustar**4)
      if( rfo .gt. 0.0 ) then
        rfo = -1.0*rfo
      endif
C
C     Obtain the time scale :
      args = visc*zo/ustar**3
C     Shear-driven timescale ( sec. ) :
      shear = csh*sqrt(args)
      argc = visc*rho*cp/(alpha*9.810*abs(qnet))
C     Convective timescale ( sec. ) :
      conv = ccon*sqrt(argc)
      rent = shear + (conv - shear)*exp(-rfcr/rfo)
C
      return
      end
C
C
C     Program   zeta.f
C     Determines the dimensionless height zeta ( zln ).
C     Inputs are :
C                 air temp. in deg C  =  airt
C                 the specific humidity  =  qa
C                 the measurement height for wind in m  =  zu
C                 ustar
C                 thetastar
C                 qstar
C     Returns the dimensionless height in m  =  zln
C
      subroutine zeta(airt,qa,ustar,thetastar,qstar,zu,zln)
      implicit none
      real*4 airt,qa,ustar,thetastar,qstar,zu,zln
      real*4 vonk,ga,ta,vt,vtstar,ob
C
C     Some constants.
C     The von Karman constant :
      data vonk /0.40/
C     Acceleration due to gravity at earth's surface ( m/s^2 ) :
      data ga /9.810/
C
C     Convert to absolute temperature :
      ta = airt + 273.15
C     The virtual temperature :
      vt = ta*(1.0 + 0.608*qa)
C     New virtual temperature star :
      vtstar = thetastar*(1.0 + 0.608*qa) + 0.608*ta*qstar
C
      if( vtstar .lt. 0.0 .or. vtstar .gt. 0.0 ) then
C     The local Monin-Obukhov length :
        ob = vt*ustar*ustar/(ga*vonk*vtstar)
        zln = zu/ob
      else
        zln = 0.0
      endif
C
      return
      end
C
C
C     Program   zom.f
C     Calculates  zo  using a stress parameterization
C     process thats includes shear due to ocean waves.
C     Inputs are :
C                 sea surface temp. in deg. C  =  sst
C                 wind speed in m/s  =  wdsp
C                 measurement height in m  =  zu
C                 ustar
C                 moninv'
C         parameters - n, wa, usfc, betag, betac, betas .
C     Returns  zom  in m .
C     Called from subroutine  getzo( )
C     and makes calls to subroutine  getphiu( )  .
C
      subroutine getzom(sst,wdsp,zref,n,betag,betac,betas,ustar,wa,
     &   usfc,moninv,zom)
      implicit none
      integer*2 n
      real*4 sst,wdsp,zref,betag,betac,betas,ustar
      real*4 wa,usfc,moninv,zom
      real*4 cmin,dum,hp,hzzg
      real*4 wavelen,z0
      real*4 denwat,sfcten,zz,period,hrat,drat
      real*4 ww,wweql,pi,ga,phiu,xbh,xhz
C
C     Note - the kinematic viscosity of air used here is
C            nu = 1.50E-05  m^2/s .
C     Some values :
      pi = 3.14159265
      ga = 9.810
      hrat = 13.333
      drat = 13.333
      ww = 1.00
      wweql = 0.81
C     The density of sea water in kg/m^3  :
      denwat = 1027.0
      zz = 10000.0
C
C     The surface tension of the vapor-liquid interface for water in N/m :
      sfcten = (-1.55E-04 * (sst + 273.15)) + 0.118
C     The minimum phase speed for surface waves in m/s :
      cmin = sqrt( 2.0 * sqrt( ga*sfcten / denwat ))
C
C     This is a strange construction -
  77  if( betag .gt. 0.0 .and. (wa*ustar) .lt. cmin ) then
        betag = 0.0
        betac = 0.0
        betas = 1.0
      endif
      cmin = sqrt( 2.0 * sqrt( ga*sfcten / denwat ))
      z0 = betas * 1.50E-05 * 0.110/abs(ustar) +
     &    betac * 0.019 * sfcten / (ustar*ustar*denwat) +
     &    0.480 / wa * ustar * ustar / ga
C     Determine the wave characteristics :
      dum = ustar*ustar*wa*wa
      if( dum .gt. cmin*cmin) then
        wavelen = pi*(dum + (sqrt(dum*dum - cmin**4)))/ga
      else
        wavelen = pi*(2.0*dum)/ga
      endif
      period = wavelen / (wa * ustar)
C     Modified Kitaigorodskii's relation between wave height
C     and wave age.
C         hp = hart*zref/zz * exp(IV * wa * ww * wweql)
C     Toba's relation between wave height, friction velocity
C     and period.
      hp = 0.062 * sqrt(ga*ustar*period**3)
C     Limit the height of breaking waves.
      if( hp .gt. (0.137 * wavelen) ) then
        hp = 0.137 * wavelen
      endif
C     Determine the mean water surface speed.
      usfc = betag * ustar * wa * (pi*hp/wavelen)**2 +
     &       (1.0 - betag) * 0.55 * ustar
C     Determine wave age.
      xbh = betag * hp
      xhz = hp / z0
      call getphiu(xbh,moninv,xhz,phiu)
      hzzg = hp * zz / zref
      if( hzzg .lt. 0.00001 ) then
        hzzg = 0.00001
      endif
      wa = (log(hzzg + 1.0) + phiu) / (0.40 * ww * wweql)
C     Insure that wa is not horribly unreasonable.
      if( wa .gt. 150.0 ) then
        wa = 150.0
      endif
      if( wa .lt. 1.10 ) then
        wa = 1.08
      endif
C     Determine zref/z0 .
      zz = zref / ((betas * 1.50E-05 * 0.110/abs(ustar)) +
     &     betac * 0.019 * sfcten / (ustar*ustar*denwat) +
     &     betag * 0.480 / wa * ustar * ustar / ga )
C     Insure that the roughness length is not too absurd.
      if( zref .gt. ( zz * 25.0 ) ) then
        zz = zref/25.0
      endif
      if( (zref/zz) .gt. 1.0E-03 ) then
        zz = zref/1.0E-03
      endif
      if( betag .gt. 0.0 .and. (wa*ustar) .lt. cmin ) then
        goto 77
      endif
C
      zom = zref/zz
      return
      end
C
C
C     Program   phiu.f
C     Calculates stability function.
C     Called from  getzom( )
C
      subroutine getphiu(zref,moninv,zz,phiu)
      implicit none
      real*4 zref,moninv,zz,phiu
      real*4 zeta,zeta0
C
      if( moninv .lt. 0.0 ) then
        zeta = (1.0 - 15.0*zref*moninv)**0.250
        zeta0 = (1.0 - 15.0*zref/zz*moninv)**0.250
        phiu = log( (zeta0*zeta0 + 1.0)*(zeta0 + 1.0)**2 /
     &    ((zeta*zeta + 1.0)*(zeta + 1.0)**2) ) +
     &    2.0 * (atan(zeta) - atan(zeta0))
      else
        phiu = 7.0 * zref * moninv
      endif
C
      return
      end
C
C
C     Program cac-new-qs.f
CC    Not used at the present time !
C     Calculates the saturation specific humidity (qs)
C     for water vapor in air for a given temperature
C     in deg. C and air pressure in mbars ( hPa ).
C     Uses Carol Anne's ' simple version '.
C      subroutine newqs(temp,pres,qs)
C      implicit none
C      real*4 temp,pres,qs,eterm,rqs
C      eterm = 6.1121*exp(17.502*temp/(240.97 + temp))
C      rqs = (1.0007 + 3.46E-06*pres)*eterm*0.980
C      qs = 0.62197*(rqs/(pres - 0.378*rqs))
C      return
C      end
C
C     End of program.