pkg/cheapaml/cheapaml_coare3_flux.F

ced0783fba Jean*0001 #include "CHEAPAML_OPTIONS.h"
                0002
27b07f3033 Jean*0003 C--   File cheapaml_coare3_flux.F:
                0004 C--    Contents:
                0005 C--    o CHEAPAML_COARE3_FLUX
                0006 C--    o PSIU (Function)
                0007 C--    o PSIT (Function)
                0008
                0009 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0010
                0011 CBOP
                0012 C     !ROUTINE: CHEAPAML_COARE3_FLUX
ced0783fba Jean*0013 C     !INTERFACE:
27b07f3033 Jean*0014       SUBROUTINE CHEAPAML_COARE3_FLUX(
8fd83faf35 Jean*0015      I                    i,j,bi,bj, iceOrNot,
                0016      I                    tSurf, windSq,
b4dc6cd434 Jean*0017      O                    hf, ef, evap, Rnl, ssqt, q100, cdq, cdu,
8fd83faf35 Jean*0018      O                    dSensdTs, dEvapdTs, dLWdTs, dQAdTs,
0f79ecb0a7 Jean*0019      I                    myIter, myThid )
27b07f3033 Jean*0020
                0021 C     !DESCRIPTION:
                0022
                0023 C     !USES:
                0024       IMPLICIT NONE
ced0783fba Jean*0025 C     === Global variables ===
                0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
                0029 #include "CHEAPAML.h"
                0030
27b07f3033 Jean*0031 C     !INPUT PARAMETERS:
8fd83faf35 Jean*0032 C     i, j     :: local indices of current grid-point
                0033 C     bi, bj   :: current tile indices
                0034 C     iceOrNot :: 0=open water, 1=ice cover, 2=ice+snow
                0035 C     tSurf    :: surface temperature
                0036 C     windSq   :: relative wind (vs surface motion) speed square
                0037 C     myIter   :: Current iteration number in simulation
                0038 C     myThid   :: My Thread Id number
27b07f3033 Jean*0039       INTEGER i,j,bi,bj
8fd83faf35 Jean*0040       INTEGER iceOrNot
                0041       _RL tSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0042       _RL windSq(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0f79ecb0a7 Jean*0043       INTEGER myIter, myThid
27b07f3033 Jean*0044 C     !OUTPUT PARAMETERS:
b4dc6cd434 Jean*0045 C     cdu :: surface drag coeff (for wind stress)
8fd83faf35 Jean*0046       _RL hf, ef, evap, Rnl, ssqt, q100, cdq
                0047       _RL cdu
                0048 C     derivative vs surf. temp of Sensible, Evap, LW, q100
                0049       _RL dSensdTs, dEvapdTs, dLWdTs, dQAdTs
27b07f3033 Jean*0050 CEOP
                0051
                0052 C     !LOCAL VARIABLES:
                0053       INTEGER iter,nits
b4dc6cd434 Jean*0054       _RL tau,L,psu,pst,Bf
                0055       _RL CD,usr,tsr,qsr
                0056 c     _RL ttas,ttt,ttt2,pt,essqt
                0057       _RL zo,zot,zoq,RR,zL
                0058       _RL twoPI,cwave,lwave
a897f57321 Nico*0059
27b07f3033 Jean*0060 C various constants
b4dc6cd434 Jean*0061       _RL u,q,Tas,tta,zi,es,qs,tsw
ced0783fba Jean*0062       _RL psiu,psit,zot10,Ct10,CC,Ribu
                0063       _RL Du,Wg,Dt,Dq,u10,zo10,Cd10,Ch10
                0064       _RL xBeta,visa,Ribcu,QaR
b4dc6cd434 Jean*0065       _RL Ct,zetu,L10,charn
2616d73cb2 Nico*0066
27b07f3033 Jean*0067 C Constants and coefficients (Stull 1988 p640).
                0068       xBeta = 1.2 _d 0    !Given as 1.25 in Fairall et al.(1996)
                0069       twoPI = 2. _d 0*PI
                0070       visa = 1.326 _d -5
                0071 C default relative humidity
                0072       QaR = 0.8 _d 0
ced0783fba Jean*0073
27b07f3033 Jean*0074 C sea surface temperature without skin correction
8fd83faf35 Jean*0075 c     tsw=theta(i,j,1,bi,bj)
                0076       tsw = tSurf(i,j)
                0077       Tas = Tair(i,j,bi,bj)
ced0783fba Jean*0078
27b07f3033 Jean*0079 C net upward long wave
                0080       Rnl = 0.96 _d 0*(stefan*(tsw+celsius2K)**4) !Net longwave (up = +).
4fa4901be6 Nico*0081
27b07f3033 Jean*0082 C Teten''s return s air svp es in mb
                0083       es = (1.0007 _d 0 + 3.46 _d -6*p0)*6.1121 _d 0
                0084      &     *EXP( 17.502 _d 0*tsw/(240.97 _d 0+tsw) )
                0085       es = es*0.98 _d 0             !reduced for salinity Kraus 1972 p. 46
b4dc6cd434 Jean*0086 C-    convert from mb to spec. humidity  kg/kg
27b07f3033 Jean*0087       qs = 0.62197 _d 0*es/(p0 -0.378 _d 0*es)
8fd83faf35 Jean*0088
27b07f3033 Jean*0089       tta = Tas+celsius2K
b4dc6cd434 Jean*0090 c     ttas=tta+gamma_blk*zt
                0091 c     ttt=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk
                0092 c     ttt2=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk-celsius2K
                0093 c     pt = p0*(1.-gamma_blk*CheapHgrid(i,j,bi,bj)/ttas)
                0094 c    &     **(gravity/gamma_blk/gasR)
                0095 c     essqt = (1.0007 _d 0 + 3.46 _d -6*pt)*6.1121 _d 0
                0096 c    &        *EXP( 17.502 _d 0*ttt2/(240.97 _d 0+ttt2) )
                0097 C-    convert from mb to spec. humidity  kg/kg
                0098 c     ssqt = 0.62197 _d 0*essqt/(pt -0.378 _d 0*essqt)
                0099 C-     LANL formulation
2616d73cb2 Nico*0100 C     saturation no more at the top:
27b07f3033 Jean*0101       ssqt=ssq0*EXP( lath*(ssq1-ssq2/tta) ) / p0
2616d73cb2 Nico*0102
27b07f3033 Jean*0103       IF (useFreshWaterFlux) THEN
                0104         q=qair(i,j,bi,bj)
                0105       ELSE
                0106         q=QaR*ssqt
                0107       ENDIF
ced0783fba Jean*0108
27b07f3033 Jean*0109 C Wave parameters
                0110       cwave=gravity*wavesp(i,j,bi,bj)/twoPI
                0111       lwave=cwave*wavesp(i,j,bi,bj)
ced0783fba Jean*0112
27b07f3033 Jean*0113 C Initial guesses
                0114       zo = 0.0001 _d 0
                0115       Wg = 0.5 _d 0                 !Gustiness factor initial guess
4fa4901be6 Nico*0116
27b07f3033 Jean*0117 C Air-sea differences - includes warm layer in Dt and Dq
8fd83faf35 Jean*0118 c     u = (uwind(i,j,bi,bj)-uVel(i,j,1,bi,bj))**2
                0119 c    &  + (vwind(i,j,bi,bj)-vVel(i,j,1,bi,bj))**2
                0120       u = windSq(i,j)
de622c6057 Jean*0121       Du= SQRT(u + Wg**2 )  !include gustiness in wind spd. difference
27b07f3033 Jean*0122       u = SQRT(u)
a897f57321 Nico*0123       Dt=tsw-Tas-gamma_blk*zt  !potential temperature difference.
0a4434c42b Jean*0124       Dq=qs-q
27b07f3033 Jean*0125
                0126 C **************** neutral coefficients ******************
                0127
                0128       u10 = Du*LOG(10. _d 0/zo)/LOG(zu/zo)
                0129       usr = 0.035 _d 0*u10
                0130       zo10= 0.011 _d 0*usr*usr/gravity+0.11 _d 0*visa/usr
                0131       Cd10= (xkar/LOG(10. _d 0/zo10))**2
                0132       Ch10= 0.00115 _d 0
                0133       Ct10= Ch10/SQRT(Cd10)
                0134       zot10=10. _d 0/EXP(xkar/Ct10)
                0135       Cd = (xkar/LOG(zu/zo10))**2
                0136
                0137 C standard coare3 boundary layer height
ced0783fba Jean*0138       zi=600. _d 0
                0139
27b07f3033 Jean*0140 C ************* Grachev and Fairall (JAM, 1997) **********
                0141
                0142       Ct=xkar/LOG(zt/zot10)   ! Temperature transfer coefficient
a897f57321 Nico*0143       CC=xkar*Ct/Cd            ! z/L vs Rib linear coefficient
0a4434c42b Jean*0144       Ribcu=-zu/(zi*0.004 _d 0*xBeta**3)  ! Saturation or plateau Rib
b4dc6cd434 Jean*0145       Ribu=-gravity*zu*(Dt+0.61 _d 0*tta*Dq)/(tta*Du**2)
27b07f3033 Jean*0146       IF (Ribu.LT.0. _d 0) THEN
ced0783fba Jean*0147           zetu=CC*Ribu/(1. _d 0+Ribu/Ribcu)   ! Unstable G and F
27b07f3033 Jean*0148       ELSE
ced0783fba Jean*0149           zetu=CC*Ribu*(1. _d 0 +27. _d 0/9. _d 0*Ribu/CC) ! Stable
27b07f3033 Jean*0150       ENDIF
ced0783fba Jean*0151       L10=zu/zetu                       ! MO length
27b07f3033 Jean*0152       IF (zetu.GT.50. _d 0) THEN
ced0783fba Jean*0153         nits=1
27b07f3033 Jean*0154       ELSE
ced0783fba Jean*0155         nits=3   ! number of iterations
27b07f3033 Jean*0156       ENDIF
                0157
                0158 C First guess M-O stability dependent
                0159 C scaling params.(u*,t*,q*) to estimate zo and z/L
                0160
                0161       usr= Du*xkar/(LOG(zu/zo10)-psiu(zu/L10))
                0162       tsr=-(Dt)*xkar/(LOG(zt/zot10)-psit(zt/L10))
                0163       qsr=-(Dq)*xkar/(LOG(zq/zot10)-psit(zq/L10))
                0164
0a4434c42b Jean*0165       charn=0.011 _d 0     !then modify Charnock for high wind speeds Chris data
27b07f3033 Jean*0166       IF (Du.GT.10. _d 0) charn=0.011 _d 0
                0167      &                    + (0.018 _d 0-0.011 _d 0)*(Du-10.)/(18.-10.)
                0168       IF (Du.GT.18. _d 0) charn=0.018 _d 0
                0169
                0170 C **** Iterate across u*(t*,q*),zo(zot,zoq) and z/L including cool skin ****
                0171
                0172       DO iter=1,nits
                0173        IF (WAVEMODEL.EQ.'Smith') THEN
ced0783fba Jean*0174         zo=charn*usr*usr/gravity + 0.11 _d 0*visa/usr    !after Smith 1988
27b07f3033 Jean*0175        ELSEIF (WAVEMODEL.EQ.'Oost') THEN
                0176         zo=(50./twoPI)*lwave*(usr/cwave)**4.5 _d 0
                0177      &    + 0.11 _d 0*visa/usr !Oost et al.
                0178        ELSEIF (WAVEMODEL.EQ.'TayYel') THEN
ced0783fba Jean*0179         zo=1200. _d 0*wavesh(i,j,bi,bj)*(wavesh(i,j,bi,bj)/lwave)**4.5
27b07f3033 Jean*0180      &    + 0.11 _d 0*visa/usr !Taylor and Yelland
                0181        ENDIF
                0182        rr=zo*usr/visa
                0183
                0184 C *** zoq and zot fitted to results from several ETL cruises ************
                0185
0f79ecb0a7 Jean*0186        IF ( rr.LE.0. ) THEN
                0187          WRITE(errorMessageUnit,'(A,I8,I4,A,5I4)')
                0188      &    'CHEAPAML_COARE3_FLUX: myIter,iter=', myIter, iter,
                0189      &    ' , in: i,j,bi,bj,thid=', i, j, bi, bj, myThid
                0190          WRITE(errorMessageUnit,'(A,1P4E17.9)')
                0191      &    ' rr,zo,usr,visa=', rr, zo, usr, visa
                0192          WRITE(errorMessageUnit,'(A,1P4E17.9)')
                0193      &    ' L,zu,zL,zt    =', L, zu, zL, zt
                0194          WRITE(errorMessageUnit,'(A,1P4E16.8)')
                0195      &    ' ln(zu/zo),psu,diff,zL*=', LOG(zu/zo), psu, LOG(zu/zo)-psu,
b4dc6cd434 Jean*0196      &          ( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr )
                0197      &         /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) )
0f79ecb0a7 Jean*0198          WRITE(errorMessageUnit,'(A,1P4E17.9)')
b4dc6cd434 Jean*0199      &    ' tsr,tta,q,qsr  =', tsr, tta, q, qsr
0f79ecb0a7 Jean*0200          CALL MDS_FLUSH( errorMessageUnit, myThid )
                0201          CALL MDS_FLUSH( standardMessageUnit, myThid )
                0202        ENDIF
27b07f3033 Jean*0203        zoq = MIN( 1.15 _d -4, 5.5 _d -5/rr**0.6 _d 0 )
                0204        zot = zoq
                0205
b4dc6cd434 Jean*0206        zL=xkar*gravity*zu*( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr )
                0207      &                   /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) )
27b07f3033 Jean*0208        L=zu/zL
                0209        psu=psiu(zu/L)
                0210        pst=psit(zt/L)
                0211        usr=Du*xkar/(LOG(zu/zo)-psiu(zu/L))
                0212        tsr=-(Dt)*xkar/(LOG(zt/zot)-psit(zt/L))
                0213        qsr=-(Dq)*xkar/(LOG(zq/zoq)-psit(zq/L))
b4dc6cd434 Jean*0214        Bf=-gravity/tta*usr*(tsr+0.61 _d 0*tta*qsr)
27b07f3033 Jean*0215        IF (Bf.GT.0. _d 0) THEN
ced0783fba Jean*0216           Wg=xBeta*(Bf*zi)**.333 _d 0
27b07f3033 Jean*0217        ELSE
ced0783fba Jean*0218           Wg=0.2 _d 0
27b07f3033 Jean*0219        ENDIF
                0220        Du=SQRT(u**2 + Wg**2)        !include gustiness in wind spd.
                0221       ENDDO
                0222
                0223 C compute surface fluxes and other parameters
6f33a64124 Jean*0224       tau=rhoa*usr*usr               !stress N/m2
b4dc6cd434 Jean*0225       hf=-cpair*rhoa*usr*tsr         !sensible W/m2
                0226       ef=-lath*rhoa*usr*qsr          !latent W/m2
27b07f3033 Jean*0227       evap=-rhoa*usr*qsr
                0228       cdq = evap/Dq
de622c6057 Jean*0229       cdu = tau/Du
8fd83faf35 Jean*0230
27b07f3033 Jean*0231       q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L))
                0232
8fd83faf35 Jean*0233 C--   compute derivative of surface fluxes relatice to Tsurf
                0234 C-    dSensdTs = -cpair*rhoa*usr*(tsr/Dt)
                0235       dSensdTs = cpair*rhoa*usr
                0236      &                *xkar/(LOG(zt/zot10)-psit(zt/L10))
                0237 C-    dEvapdTs  = -rhoa*usr* d/dTs(qsr)
                0238 C     d/dTs(qsr)= (-xkar/(LOG(zq/zoq)-psit(zq/L)) )* d/dTs(qs)
                0239 C     d/dTs(qs) = 0.62197 _d 0*p0/(p0 -0.378 _d 0*es)**2 * d/dTs(es)
                0240 C     d/dTs(es) = (0.98)* es * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2
                0241 C-    this simplifies (using qs) to:
                0242       dEvapdTs = rhoa*usr*( xkar/(LOG(zq/zoq)-psit(zq/L)) )
                0243      &         * qs*p0/(p0 -0.378 _d 0*es)
                0244 c    &         *0.98 _d 0
                0245      &         * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2
                0246
                0247       if (iceornot.EQ.0) THEN
                0248 c       dLWdTs = 4. _d 0*ocean_emissivity*stefan*tsr*tsr*tsr
                0249         dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
                0250       ELSEIF (iceornot.EQ.2) THEN
                0251 c       dLWdTs = 4. _d 0*snow_emissivity*stefan*tsr*tsr*tsr
                0252         dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
                0253       ELSEIF (iceornot.EQ.1) THEN
                0254 c       dLWdTs = 4. _d 0*ice_emissivity*stefan*tsr*tsr*tsr
                0255         dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr
                0256       ENDIF
                0257
                0258 C--   for now, ignores derivative of q100 relative to Tsurf:
                0259       dQAdTs = 0.
                0260
27b07f3033 Jean*0261       RETURN
                0262       END
                0263
                0264 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0265
                0266 CBOP 0
                0267 C     !ROUTINE: PSIU
                0268
                0269 C     !INTERFACE:
                0270       _RL FUNCTION psiu(zL)
                0271
                0272 C     !DESCRIPTION:
                0273 C psiu and psit evaluate stability function for wind speed and scalars
                0274 C matching Kansas and free convection forms with weighting f
                0275 C convective form follows Fairall et al (1996) with profile constants
                0276 C from Grachev et al (2000) BLM
                0277 C stable form from Beljaars and Holtslag (1991)
                0278
                0279 C     !USES:
                0280       IMPLICIT NONE
                0281 #include "EEPARAMS.h"
                0282
                0283 C     !INPUT PARAMETERS:
                0284       _RL zL
                0285 C     !LOCAL VARIABLES:
                0286       _RL x,y,psik,psic,f,c
                0287 CEOP
                0288
                0289       IF (zL.LT.0.0) THEN
                0290        x = (1. - 15.*zL)**.25                   !Kansas unstable
                0291        psik=2.*LOG((1.+x)/2.)+LOG((1.+x*x)/2.)-2.*ATAN(x)+2.*ATAN(oneRL)
                0292        y = (1. - 10.15 _d 0*zL)**.3333 _d 0     !Convective
                0293        psic = 1.5*LOG((1.+y+y*y)/3.)
                0294      &      - SQRT(3. _d 0)*ATAN( (1.+2.*y)/SQRT(3. _d 0) )
                0295      &      + 4.*ATAN(oneRL)/SQRT(3. _d 0)
                0296        f = zL*zL/(1.+zL*zL)
                0297        psiu = (1.-f)*psik+f*psic
                0298       ELSE
                0299        c = MIN( 50. _d 0, 0.35 _d 0*zL )        !Stable
                0300 c      psiu=-((1.+1.*zL)**1.+.6667*(zL-14.28)/EXP(c)+8.525)
                0301        psiu = -( (1.+zL) + 0.6667 _d 0*(zL-14.28 _d 0)/EXP(c)
                0302      &          + 8.525 _d 0 )
                0303       ENDIF
                0304       RETURN
                0305       END
                0306
                0307 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0308
                0309 CBOP 0
                0310 C     !ROUTINE: PSIT
                0311
                0312 C     !INTERFACE:
                0313       _RL FUNCTION psit(zL)
                0314
                0315 C     !DESCRIPTION:
                0316
                0317 C     !USES:
                0318       IMPLICIT NONE
                0319 #include "EEPARAMS.h"
                0320
                0321 C     !INPUT PARAMETERS:
                0322       _RL zL
                0323 C     !LOCAL VARIABLES:
                0324       _RL x,y,psik,psic,f,c
                0325 CEOP
                0326
                0327       IF (zL.LT.0.0) THEN
                0328        x = (1. - 15.*zL)**.5                    !Kansas unstable
                0329        psik = 2.*LOG((1.+x)/2.)
                0330        y = (1. - 34.15 _d 0*zL)**.3333 _d 0     !Convective
                0331        psic = 1.5*LOG((1.+y+y*y)/3.)
                0332      &      - SQRT(3. _d 0)*ATAN( (1.+2.*y)/SQRT(3. _d 0) )
                0333      &      + 4.*ATAN(oneRL)/SQRT(3. _d 0)
                0334        f = zL*zL/(1.+zL*zL)
                0335        psit = (1.-f)*psik+f*psic
                0336       ELSE
                0337        c = MIN( 50. _d 0, 0.35 _d 0*zL )        !Stable
                0338        psit = -( (1.+2.*zL/3.)**1.5
                0339      &          + 0.6667 _d 0*(zL-14.28 _d 0)/EXP(c)
                0340      &          + 8.525 _d 0 )
                0341       ENDIF
                0342
                0343       RETURN
                0344       END