C $Header: /u/gcmpack/MITgcm/pkg/cheapaml/cheapaml_coare3_flux.F,v 1.16 2013/06/17 13:45:14 jmc Exp $ C $Name: $ #include "CHEAPAML_OPTIONS.h" C-- File cheapaml_coare3_flux.F: C-- Contents: C-- o CHEAPAML_COARE3_FLUX C-- o PSIU (Function) C-- o PSIT (Function) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: CHEAPAML_COARE3_FLUX C !INTERFACE: SUBROUTINE CHEAPAML_COARE3_FLUX( I i,j,bi,bj, iceOrNot, I tSurf, windSq, O hf, ef, evap, Rnl, ssqt, q100, cdq, cdu, O dSensdTs, dEvapdTs, dLWdTs, dQAdTs, I myIter, myThid ) C !DESCRIPTION: C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "CHEAPAML.h" C !INPUT PARAMETERS: C i, j :: local indices of current grid-point C bi, bj :: current tile indices C iceOrNot :: 0=open water, 1=ice cover, 2=ice+snow C tSurf :: surface temperature C windSq :: relative wind (vs surface motion) speed square C myIter :: Current iteration number in simulation C myThid :: My Thread Id number INTEGER i,j,bi,bj INTEGER iceOrNot _RL tSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL windSq(1-OLx:sNx+OLx,1-OLy:sNy+OLy) INTEGER myIter, myThid C !OUTPUT PARAMETERS: C cdu :: surface drag coeff (for wind stress) _RL hf, ef, evap, Rnl, ssqt, q100, cdq _RL cdu C derivative vs surf. temp of Sensible, Evap, LW, q100 _RL dSensdTs, dEvapdTs, dLWdTs, dQAdTs CEOP C !LOCAL VARIABLES: INTEGER iter,nits _RL tau,L,psu,pst,Bf _RL CD,usr,tsr,qsr c _RL ttas,ttt,ttt2,pt,essqt _RL zo,zot,zoq,RR,zL _RL twoPI,cwave,lwave C various constants _RL u,q,Tas,tta,zi,es,qs,tsw _RL psiu,psit,zot10,Ct10,CC,Ribu _RL Du,Wg,Dt,Dq,u10,zo10,Cd10,Ch10 _RL xBeta,visa,Ribcu,QaR _RL Ct,zetu,L10,charn C Constants and coefficients (Stull 1988 p640). xBeta = 1.2 _d 0 !Given as 1.25 in Fairall et al.(1996) twoPI = 2. _d 0*PI visa = 1.326 _d -5 C default relative humidity QaR = 0.8 _d 0 C sea surface temperature without skin correction c tsw=theta(i,j,1,bi,bj) tsw = tSurf(i,j) Tas = Tair(i,j,bi,bj) C net upward long wave Rnl = 0.96 _d 0*(stefan*(tsw+celsius2K)**4) !Net longwave (up = +). C Teten''s return s air svp es in mb es = (1.0007 _d 0 + 3.46 _d -6*p0)*6.1121 _d 0 & *EXP( 17.502 _d 0*tsw/(240.97 _d 0+tsw) ) es = es*0.98 _d 0 !reduced for salinity Kraus 1972 p. 46 C- convert from mb to spec. humidity kg/kg qs = 0.62197 _d 0*es/(p0 -0.378 _d 0*es) tta = Tas+celsius2K c ttas=tta+gamma_blk*zt c ttt=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk c ttt2=tta-(CheapHgrid(i,j,bi,bj) - zt)*gamma_blk-celsius2K c pt = p0*(1.-gamma_blk*CheapHgrid(i,j,bi,bj)/ttas) c & **(gravity/gamma_blk/gasR) c essqt = (1.0007 _d 0 + 3.46 _d -6*pt)*6.1121 _d 0 c & *EXP( 17.502 _d 0*ttt2/(240.97 _d 0+ttt2) ) C- convert from mb to spec. humidity kg/kg c ssqt = 0.62197 _d 0*essqt/(pt -0.378 _d 0*essqt) C- LANL formulation C saturation no more at the top: ssqt=ssq0*EXP( lath*(ssq1-ssq2/tta) ) / p0 IF (useFreshWaterFlux) THEN q=qair(i,j,bi,bj) ELSE q=QaR*ssqt ENDIF C Wave parameters cwave=gravity*wavesp(i,j,bi,bj)/twoPI lwave=cwave*wavesp(i,j,bi,bj) C Initial guesses zo = 0.0001 _d 0 Wg = 0.5 _d 0 !Gustiness factor initial guess C Air-sea differences - includes warm layer in Dt and Dq c u = (uwind(i,j,bi,bj)-uVel(i,j,1,bi,bj))**2 c & + (vwind(i,j,bi,bj)-vVel(i,j,1,bi,bj))**2 u = windSq(i,j) Du= SQRT(u + Wg**2 ) !include gustiness in wind spd. difference u = SQRT(u) Dt=tsw-Tas-gamma_blk*zt !potential temperature difference. Dq=qs-q C **************** neutral coefficients ****************** u10 = Du*LOG(10. _d 0/zo)/LOG(zu/zo) usr = 0.035 _d 0*u10 zo10= 0.011 _d 0*usr*usr/gravity+0.11 _d 0*visa/usr Cd10= (xkar/LOG(10. _d 0/zo10))**2 Ch10= 0.00115 _d 0 Ct10= Ch10/SQRT(Cd10) zot10=10. _d 0/EXP(xkar/Ct10) Cd = (xkar/LOG(zu/zo10))**2 C standard coare3 boundary layer height zi=600. _d 0 C ************* Grachev and Fairall (JAM, 1997) ********** Ct=xkar/LOG(zt/zot10) ! Temperature transfer coefficient CC=xkar*Ct/Cd ! z/L vs Rib linear coefficient Ribcu=-zu/(zi*0.004 _d 0*xBeta**3) ! Saturation or plateau Rib Ribu=-gravity*zu*(Dt+0.61 _d 0*tta*Dq)/(tta*Du**2) IF (Ribu.LT.0. _d 0) THEN zetu=CC*Ribu/(1. _d 0+Ribu/Ribcu) ! Unstable G and F ELSE zetu=CC*Ribu*(1. _d 0 +27. _d 0/9. _d 0*Ribu/CC) ! Stable ENDIF L10=zu/zetu ! MO length IF (zetu.GT.50. _d 0) THEN nits=1 ELSE nits=3 ! number of iterations ENDIF C First guess M-O stability dependent C scaling params.(u*,t*,q*) to estimate zo and z/L usr= Du*xkar/(LOG(zu/zo10)-psiu(zu/L10)) tsr=-(Dt)*xkar/(LOG(zt/zot10)-psit(zt/L10)) qsr=-(Dq)*xkar/(LOG(zq/zot10)-psit(zq/L10)) charn=0.011 _d 0 !then modify Charnock for high wind speeds Chris data IF (Du.GT.10. _d 0) charn=0.011 _d 0 & + (0.018 _d 0-0.011 _d 0)*(Du-10.)/(18.-10.) IF (Du.GT.18. _d 0) charn=0.018 _d 0 C **** Iterate across u*(t*,q*),zo(zot,zoq) and z/L including cool skin **** DO iter=1,nits IF (WAVEMODEL.EQ.'Smith') THEN zo=charn*usr*usr/gravity + 0.11 _d 0*visa/usr !after Smith 1988 ELSEIF (WAVEMODEL.EQ.'Oost') THEN zo=(50./twoPI)*lwave*(usr/cwave)**4.5 _d 0 & + 0.11 _d 0*visa/usr !Oost et al. ELSEIF (WAVEMODEL.EQ.'TayYel') THEN zo=1200. _d 0*wavesh(i,j,bi,bj)*(wavesh(i,j,bi,bj)/lwave)**4.5 & + 0.11 _d 0*visa/usr !Taylor and Yelland ENDIF rr=zo*usr/visa C *** zoq and zot fitted to results from several ETL cruises ************ IF ( rr.LE.0. ) THEN WRITE(errorMessageUnit,'(A,I8,I4,A,5I4)') & 'CHEAPAML_COARE3_FLUX: myIter,iter=', myIter, iter, & ' , in: i,j,bi,bj,thid=', i, j, bi, bj, myThid WRITE(errorMessageUnit,'(A,1P4E17.9)') & ' rr,zo,usr,visa=', rr, zo, usr, visa WRITE(errorMessageUnit,'(A,1P4E17.9)') & ' L,zu,zL,zt =', L, zu, zL, zt WRITE(errorMessageUnit,'(A,1P4E16.8)') & ' ln(zu/zo),psu,diff,zL*=', LOG(zu/zo), psu, LOG(zu/zo)-psu, & ( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr ) & /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) ) WRITE(errorMessageUnit,'(A,1P4E17.9)') & ' tsr,tta,q,qsr =', tsr, tta, q, qsr CALL MDS_FLUSH( errorMessageUnit, myThid ) CALL MDS_FLUSH( standardMessageUnit, myThid ) ENDIF zoq = MIN( 1.15 _d -4, 5.5 _d -5/rr**0.6 _d 0 ) zot = zoq zL=xkar*gravity*zu*( tsr*(1.+0.61 _d 0*q)+0.61 _d 0*tta*qsr ) & /( tta*usr*usr*(1. _d 0+0.61 _d 0*q) ) L=zu/zL psu=psiu(zu/L) pst=psit(zt/L) usr=Du*xkar/(LOG(zu/zo)-psiu(zu/L)) tsr=-(Dt)*xkar/(LOG(zt/zot)-psit(zt/L)) qsr=-(Dq)*xkar/(LOG(zq/zoq)-psit(zq/L)) Bf=-gravity/tta*usr*(tsr+0.61 _d 0*tta*qsr) IF (Bf.GT.0. _d 0) THEN Wg=xBeta*(Bf*zi)**.333 _d 0 ELSE Wg=0.2 _d 0 ENDIF Du=SQRT(u**2 + Wg**2) !include gustiness in wind spd. ENDDO C compute surface fluxes and other parameters tau=rhoa*usr*usr !stress N/m2 hf=-cpair*rhoa*usr*tsr !sensible W/m2 ef=-lath*rhoa*usr*qsr !latent W/m2 evap=-rhoa*usr*qsr cdq = evap/Dq cdu = tau/Du q100=qs+qsr*(LOG(100. _d 0/zoq)-psit(100. _d 0/L)) C-- compute derivative of surface fluxes relatice to Tsurf C- dSensdTs = -cpair*rhoa*usr*(tsr/Dt) dSensdTs = cpair*rhoa*usr & *xkar/(LOG(zt/zot10)-psit(zt/L10)) C- dEvapdTs = -rhoa*usr* d/dTs(qsr) C d/dTs(qsr)= (-xkar/(LOG(zq/zoq)-psit(zq/L)) )* d/dTs(qs) C d/dTs(qs) = 0.62197 _d 0*p0/(p0 -0.378 _d 0*es)**2 * d/dTs(es) C d/dTs(es) = (0.98)* es * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2 C- this simplifies (using qs) to: dEvapdTs = rhoa*usr*( xkar/(LOG(zq/zoq)-psit(zq/L)) ) & * qs*p0/(p0 -0.378 _d 0*es) c & *0.98 _d 0 & * 17.502 _d 0 * 240.97 _d 0 / (240.97 _d 0+tsw)**2 if (iceornot.EQ.0) THEN c dLWdTs = 4. _d 0*ocean_emissivity*stefan*tsr*tsr*tsr dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr ELSEIF (iceornot.EQ.2) THEN c dLWdTs = 4. _d 0*snow_emissivity*stefan*tsr*tsr*tsr dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr ELSEIF (iceornot.EQ.1) THEN c dLWdTs = 4. _d 0*ice_emissivity*stefan*tsr*tsr*tsr dLWdTs = 4. _d 0 * 0.96 _d 0 *stefan*tsr*tsr*tsr ENDIF C-- for now, ignores derivative of q100 relative to Tsurf: dQAdTs = 0. RETURN END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP 0 C !ROUTINE: PSIU C !INTERFACE: _RL FUNCTION psiu(zL) C !DESCRIPTION: C psiu and psit evaluate stability function for wind speed and scalars C matching Kansas and free convection forms with weighting f C convective form follows Fairall et al (1996) with profile constants C from Grachev et al (2000) BLM C stable form from Beljaars and Holtslag (1991) C !USES: IMPLICIT NONE #include "EEPARAMS.h" C !INPUT PARAMETERS: _RL zL C !LOCAL VARIABLES: _RL x,y,psik,psic,f,c CEOP IF (zL.LT.0.0) THEN x = (1. - 15.*zL)**.25 !Kansas unstable psik=2.*LOG((1.+x)/2.)+LOG((1.+x*x)/2.)-2.*ATAN(x)+2.*ATAN(oneRL) y = (1. - 10.15 _d 0*zL)**.3333 _d 0 !Convective psic = 1.5*LOG((1.+y+y*y)/3.) & - SQRT(3. _d 0)*ATAN( (1.+2.*y)/SQRT(3. _d 0) ) & + 4.*ATAN(oneRL)/SQRT(3. _d 0) f = zL*zL/(1.+zL*zL) psiu = (1.-f)*psik+f*psic ELSE c = MIN( 50. _d 0, 0.35 _d 0*zL ) !Stable c psiu=-((1.+1.*zL)**1.+.6667*(zL-14.28)/EXP(c)+8.525) psiu = -( (1.+zL) + 0.6667 _d 0*(zL-14.28 _d 0)/EXP(c) & + 8.525 _d 0 ) ENDIF RETURN END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP 0 C !ROUTINE: PSIT C !INTERFACE: _RL FUNCTION psit(zL) C !DESCRIPTION: C !USES: IMPLICIT NONE #include "EEPARAMS.h" C !INPUT PARAMETERS: _RL zL C !LOCAL VARIABLES: _RL x,y,psik,psic,f,c CEOP IF (zL.LT.0.0) THEN x = (1. - 15.*zL)**.5 !Kansas unstable psik = 2.*LOG((1.+x)/2.) y = (1. - 34.15 _d 0*zL)**.3333 _d 0 !Convective psic = 1.5*LOG((1.+y+y*y)/3.) & - SQRT(3. _d 0)*ATAN( (1.+2.*y)/SQRT(3. _d 0) ) & + 4.*ATAN(oneRL)/SQRT(3. _d 0) f = zL*zL/(1.+zL*zL) psit = (1.-f)*psik+f*psic ELSE c = MIN( 50. _d 0, 0.35 _d 0*zL ) !Stable psit = -( (1.+2.*zL/3.)**1.5 & + 0.6667 _d 0*(zL-14.28 _d 0)/EXP(c) & + 8.525 _d 0 ) ENDIF RETURN END