Index: /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F =================================================================== --- /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F (revision 335) +++ /LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F (revision 336) @@ -1176,8 +1176,177 @@ c Calculer le frottement au sol (Cdrag) c + CALL clcdrag(knon, nsrf, is_oce, zxli, u, v, t, q, zgeop, + . ts, qsurf, rugos, pcfm, pcfh, zcdn, zri) +c +c Calculer les coefficients turbulents dans l'atmosphere +c + DO i = 1, knon + itop(i) = isommet + ENDDO + + DO k = 2, isommet + DO i = 1, knon + zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 + . +(v(i,k)-v(i,k-1))**2) + zmgeom(i)=zgeop(i,k)-zgeop(i,k-1) + zdphi =zmgeom(i) / 2.0 + zt = (t(i,k)+t(i,k-1)) * 0.5 + zq = (q(i,k)+q(i,k-1)) * 0.5 +c +c calculer Qs et dQs/dT: +c + IF (thermcep) THEN + zdelta = MAX(0.,SIGN(1.,RTT-zt)) + zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) + . + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta + zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k) + zqs = MIN(0.5,zqs) + zcor = 1./(1.-RETV*zqs) + zqs = zqs*zcor + zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor) + ELSE + IF (zt .LT. t_coup) THEN + zqs = qsats(zt) / pplay(i,k) + zdqs = dqsats(zt,zqs) + ELSE + zqs = qsatl(zt) / pplay(i,k) + zdqs = dqsatl(zt,zqs) + ENDIF + ENDIF +c +c calculer la fraction nuageuse (processus humide): +c + zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq) + zfr = MAX(0.0,MIN(1.0,zfr)) + IF (.NOT.richum) zfr = 0.0 +c +c calculer le nombre de Richardson: +c + IF (tvirtu) THEN + ztvd =( t(i,k) + . + zdphi/RCPD/(1.+RVTMP2*zq) + . *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) + . )*(1.+RETV*q(i,k)) + ztvu =( t(i,k-1) + . - zdphi/RCPD/(1.+RVTMP2*zq) + . *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) + . )*(1.+RETV*q(i,k-1)) + zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu)) + zri(i) = zri(i) + . + zmgeom(i)*zmgeom(i)/RG*gamt(k) + . *(paprs(i,k)/101325.0)**RKAPPA + . /(zdu2*0.5*(ztvd+ztvu)) +c + ELSE ! calcul de Ridchardson compatible LMD5 +c + zri(i) =(RCPD*(t(i,k)-t(i,k-1)) + . -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k) + . *(pplay(i,k)-pplay(i,k-1)) + . )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k))) + zri(i) = zri(i) + + . zmgeom(i)*zmgeom(i)*gamt(k)/RG +cSB . /(paprs(i,k)/101325.0)**RKAPPA + . *(paprs(i,k)/101325.0)**RKAPPA + . /(zdu2*0.5*(t(i,k-1)+t(i,k))) + ENDIF +c +c finalement, les coefficients d'echange sont obtenus: +c + zcdn=SQRT(zdu2) / zmgeom(i) * RG +c + IF (opt_ec) THEN + z2geomf=zgeop(i,k-1)+zgeop(i,k) + zalm2=(0.5*ckap/RG*z2geomf + . /(1.+0.5*ckap/rg/clam*z2geomf))**2 + zalh2=(0.5*ckap/rg*z2geomf + . /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2 + IF (zri(i).LT.0.0) THEN ! situation instable + zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3 + . / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG) + zscf = SQRT(-zri(i)*zscf) + zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf) + zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf) + pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm) + pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh) + ELSE ! situation stable + zscf=SQRT(1.+cd*zri(i)) + pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf) + pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf) + ENDIF + ELSE + zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) + . /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2 + pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable)) + pcfm(i,k)= zl2(i)* pcfm(i,k) + pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different + ENDIF + ENDDO + ENDDO +c +c Au-dela du sommet, pas de diffusion turbulente: +c + DO i = 1, knon + IF (itop(i)+1 .LE. klev) THEN + DO k = itop(i)+1, klev + pcfh(i,k) = 0.0 + pcfm(i,k) = 0.0 + ENDDO + ENDIF + ENDDO +c + RETURN + END + + SUBROUTINE clcdrag(knon, nsrf, is_oce, zxli, + . u, v, t, q, zgeop, + . ts, qsurf, rugos, + . pcfm, pcfh, zcdn, zri) +c ================================================================= c +c Objet : calcul cdrags pour le moment et les flux chaleur sensible, latente (pcfm,pcfh) +c et du nombre de Richardson zri +c ================================================================= c + IMPLICIT NONE +#include "dimensions.h" +#include "dimphy.h" +#include "YOMCST.h" +#include "indicesol.h" +c + INTEGER knon, nsrf, is_oce + REAL ts(klon), qsurf(klon) + REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev) + REAL zgeop(klon,klev) + REAL rugos(klon), zri(klon) +c + REAL zcdn + REAL pcfm(klon,klev), pcfh(klon,klev) +c +c Quelques constantes et options: +c + REAL ckap, cb, cc, cd, cepdu2 + PARAMETER (ckap=0.35) + PARAMETER (cb=5.0) + PARAMETER (cc=5.0) + PARAMETER (cd=5.0) + PARAMETER (cepdu2 =(0.1)**2) +c +c Variables locales + INTEGER i + REAL zdu2, zdphi, ztsolv, ztvd, zscf, zucf, zcr + REAL friv, frih + REAL zcfm1(klon), zcfm2(klon) + REAL zcfh1(klon), zcfh2(klon) +c +c Fonctions thermodynamiques et fonctions d'instabilite + REAL fsta, fins, x + LOGICAL zxli + fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x)) + fins(x) = SQRT(1.0-18.0*x) +c +c Calculer le frottement au sol (Cdrag) +c DO i = 1, knon zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2) zdphi=zgeop(i,1) -c ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1) +c ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1) ztsolv = ts(i) * (1.0+RETV*qsurf(i)) ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1))) @@ -1215,125 +1384,8 @@ IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25) ENDIF - ENDDO - -c -c Calculer les coefficients turbulents dans l'atmosphere -c - DO i = 1, knon - itop(i) = isommet - ENDDO - - DO k = 2, isommet - DO i = 1, knon - zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2 - . +(v(i,k)-v(i,k-1))**2) - zmgeom(i)=zgeop(i,k)-zgeop(i,k-1) - zdphi =zmgeom(i) / 2.0 - zt = (t(i,k)+t(i,k-1)) * 0.5 - zq = (q(i,k)+q(i,k-1)) * 0.5 -c -c calculer Qs et dQs/dT: -c - IF (thermcep) THEN - zdelta = MAX(0.,SIGN(1.,RTT-zt)) - zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta) - . + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta - zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k) - zqs = MIN(0.5,zqs) - zcor = 1./(1.-RETV*zqs) - zqs = zqs*zcor - zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor) - ELSE - IF (zt .LT. t_coup) THEN - zqs = qsats(zt) / pplay(i,k) - zdqs = dqsats(zt,zqs) - ELSE - zqs = qsatl(zt) / pplay(i,k) - zdqs = dqsatl(zt,zqs) - ENDIF - ENDIF -c -c calculer la fraction nuageuse (processus humide): -c - zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq) - zfr = MAX(0.0,MIN(1.0,zfr)) - IF (.NOT.richum) zfr = 0.0 -c -c calculer le nombre de Richardson: -c - IF (tvirtu) THEN - ztvd =( t(i,k) - . + zdphi/RCPD/(1.+RVTMP2*zq) - . *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) - . )*(1.+RETV*q(i,k)) - ztvu =( t(i,k-1) - . - zdphi/RCPD/(1.+RVTMP2*zq) - . *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) ) - . )*(1.+RETV*q(i,k-1)) - zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu)) - zri(i) = zri(i) - . + zmgeom(i)*zmgeom(i)/RG*gamt(k) - . *(paprs(i,k)/101325.0)**RKAPPA - . /(zdu2*0.5*(ztvd+ztvu)) -c - ELSE ! calcul de Ridchardson compatible LMD5 -c - zri(i) =(RCPD*(t(i,k)-t(i,k-1)) - . -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k) - . *(pplay(i,k)-pplay(i,k-1)) - . )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k))) - zri(i) = zri(i) + - . zmgeom(i)*zmgeom(i)*gamt(k)/RG -cSB . /(paprs(i,k)/101325.0)**RKAPPA - . *(paprs(i,k)/101325.0)**RKAPPA - . /(zdu2*0.5*(t(i,k-1)+t(i,k))) - ENDIF -c -c finalement, les coefficients d'echange sont obtenus: -c - zcdn=SQRT(zdu2) / zmgeom(i) * RG -c - IF (opt_ec) THEN - z2geomf=zgeop(i,k-1)+zgeop(i,k) - zalm2=(0.5*ckap/RG*z2geomf - . /(1.+0.5*ckap/rg/clam*z2geomf))**2 - zalh2=(0.5*ckap/rg*z2geomf - . /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2 - IF (zri(i).LT.0.0) THEN ! situation instable - zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3 - . / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG) - zscf = SQRT(-zri(i)*zscf) - zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf) - zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf) - pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm) - pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh) - ELSE ! situation stable - zscf=SQRT(1.+cd*zri(i)) - pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf) - pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf) - ENDIF - ELSE - zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1)) - . /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2 - pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable)) - pcfm(i,k)= zl2(i)* pcfm(i,k) - pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different - ENDIF - ENDDO - ENDDO -c -c Au-dela du sommet, pas de diffusion turbulente: -c - DO i = 1, knon - IF (itop(i)+1 .LE. klev) THEN - DO k = itop(i)+1, klev - pcfh(i,k) = 0.0 - pcfm(i,k) = 0.0 - ENDDO - ENDIF - ENDDO -c + END DO RETURN END + SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, . pcfm, pcfh)