Changeset 336 for LMDZ.3.3/branches


Ignore:
Timestamp:
Feb 15, 2002, 1:16:35 PM (22 years ago)
Author:
lmdzadmin
Message:

Le calcul du cdrag est inclus dans une routine pour une plus grande
modularite, IMusat, JPolcher
LF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F

    r295 r336  
    11761176c Calculer le frottement au sol (Cdrag)
    11771177c
     1178      CALL clcdrag(knon, nsrf, is_oce, zxli, u, v, t, q, zgeop,
     1179     .             ts, qsurf, rugos, pcfm, pcfh, zcdn, zri)
     1180c
     1181c Calculer les coefficients turbulents dans l'atmosphere
     1182c
     1183      DO i = 1, knon
     1184         itop(i) = isommet
     1185      ENDDO
     1186
     1187      DO k = 2, isommet
     1188      DO i = 1, knon
     1189            zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
     1190     .                     +(v(i,k)-v(i,k-1))**2)
     1191            zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
     1192            zdphi =zmgeom(i) / 2.0
     1193            zt = (t(i,k)+t(i,k-1)) * 0.5
     1194            zq = (q(i,k)+q(i,k-1)) * 0.5
     1195c
     1196c           calculer Qs et dQs/dT:
     1197c
     1198            IF (thermcep) THEN
     1199              zdelta = MAX(0.,SIGN(1.,RTT-zt))
     1200              zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)
     1201     .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
     1202              zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
     1203              zqs = MIN(0.5,zqs)
     1204              zcor = 1./(1.-RETV*zqs)
     1205              zqs = zqs*zcor
     1206              zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
     1207            ELSE
     1208              IF (zt .LT. t_coup) THEN
     1209                 zqs = qsats(zt) / pplay(i,k)
     1210                 zdqs = dqsats(zt,zqs)
     1211              ELSE
     1212                 zqs = qsatl(zt) / pplay(i,k)
     1213                 zdqs = dqsatl(zt,zqs)
     1214              ENDIF
     1215            ENDIF
     1216c
     1217c           calculer la fraction nuageuse (processus humide):
     1218c
     1219            zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
     1220            zfr = MAX(0.0,MIN(1.0,zfr))
     1221            IF (.NOT.richum) zfr = 0.0
     1222c
     1223c           calculer le nombre de Richardson:
     1224c
     1225            IF (tvirtu) THEN
     1226            ztvd =( t(i,k)
     1227     .             + zdphi/RCPD/(1.+RVTMP2*zq)
     1228     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
     1229     .            )*(1.+RETV*q(i,k))
     1230            ztvu =( t(i,k-1)
     1231     .             - zdphi/RCPD/(1.+RVTMP2*zq)
     1232     .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
     1233     .            )*(1.+RETV*q(i,k-1))
     1234            zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
     1235            zri(i) = zri(i)
     1236     .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)
     1237     .               *(paprs(i,k)/101325.0)**RKAPPA
     1238     .               /(zdu2*0.5*(ztvd+ztvu))
     1239c
     1240            ELSE ! calcul de Ridchardson compatible LMD5
     1241c
     1242            zri(i) =(RCPD*(t(i,k)-t(i,k-1))
     1243     .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)
     1244     .               *(pplay(i,k)-pplay(i,k-1))
     1245     .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
     1246            zri(i) = zri(i) +
     1247     .             zmgeom(i)*zmgeom(i)*gamt(k)/RG
     1248cSB     .             /(paprs(i,k)/101325.0)**RKAPPA
     1249     .             *(paprs(i,k)/101325.0)**RKAPPA
     1250     .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))
     1251            ENDIF
     1252c
     1253c           finalement, les coefficients d'echange sont obtenus:
     1254c
     1255            zcdn=SQRT(zdu2) / zmgeom(i) * RG
     1256c
     1257          IF (opt_ec) THEN
     1258            z2geomf=zgeop(i,k-1)+zgeop(i,k)
     1259            zalm2=(0.5*ckap/RG*z2geomf
     1260     .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
     1261            zalh2=(0.5*ckap/rg*z2geomf
     1262     .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
     1263            IF (zri(i).LT.0.0) THEN  ! situation instable
     1264               zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
     1265     .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
     1266               zscf = SQRT(-zri(i)*zscf)
     1267               zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
     1268               zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
     1269               pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
     1270               pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
     1271            ELSE ! situation stable
     1272               zscf=SQRT(1.+cd*zri(i))
     1273               pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
     1274               pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
     1275            ENDIF
     1276          ELSE
     1277            zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
     1278     .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
     1279            pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
     1280            pcfm(i,k)= zl2(i)* pcfm(i,k)
     1281            pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
     1282          ENDIF
     1283      ENDDO
     1284      ENDDO
     1285c
     1286c Au-dela du sommet, pas de diffusion turbulente:
     1287c
     1288      DO i = 1, knon
     1289         IF (itop(i)+1 .LE. klev) THEN
     1290            DO k = itop(i)+1, klev
     1291               pcfh(i,k) = 0.0
     1292               pcfm(i,k) = 0.0
     1293            ENDDO
     1294         ENDIF
     1295      ENDDO
     1296c
     1297      RETURN
     1298      END
     1299
     1300      SUBROUTINE clcdrag(knon, nsrf, is_oce, zxli,
     1301     .                   u, v, t, q, zgeop,
     1302     .                   ts, qsurf, rugos,
     1303     .                   pcfm, pcfh, zcdn, zri)
     1304c ================================================================= c
     1305c Objet : calcul cdrags pour le moment et les flux chaleur sensible, latente (pcfm,pcfh)
     1306c         et du nombre de Richardson zri
     1307c ================================================================= c
     1308      IMPLICIT NONE
     1309#include "dimensions.h"
     1310#include "dimphy.h"
     1311#include "YOMCST.h"
     1312#include "indicesol.h"
     1313c
     1314      INTEGER knon, nsrf, is_oce
     1315      REAL ts(klon), qsurf(klon)
     1316      REAL u(klon,klev), v(klon,klev), t(klon,klev), q(klon,klev)
     1317      REAL zgeop(klon,klev)
     1318      REAL rugos(klon), zri(klon)
     1319c
     1320      REAL zcdn
     1321      REAL pcfm(klon,klev), pcfh(klon,klev)
     1322c
     1323c Quelques constantes et options:
     1324c
     1325      REAL ckap, cb, cc, cd, cepdu2
     1326      PARAMETER (ckap=0.35)
     1327      PARAMETER (cb=5.0)
     1328      PARAMETER (cc=5.0)
     1329      PARAMETER (cd=5.0)
     1330      PARAMETER (cepdu2 =(0.1)**2)
     1331c
     1332c Variables locales
     1333      INTEGER i
     1334      REAL zdu2, zdphi, ztsolv, ztvd, zscf, zucf, zcr
     1335      REAL friv, frih
     1336      REAL zcfm1(klon), zcfm2(klon)
     1337      REAL zcfh1(klon), zcfh2(klon)
     1338c
     1339c Fonctions thermodynamiques et fonctions d'instabilite
     1340      REAL fsta, fins, x
     1341      LOGICAL zxli
     1342      fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x))
     1343      fins(x) = SQRT(1.0-18.0*x)
     1344c
     1345c Calculer le frottement au sol (Cdrag)
     1346c
    11781347      DO i = 1, knon
    11791348         zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2)
    11801349         zdphi=zgeop(i,1)
    1181 c         ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
     1350c        ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1)
    11821351         ztsolv = ts(i) * (1.0+RETV*qsurf(i))
    11831352         ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1)))
     
    12151384           IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25)
    12161385         ENDIF
    1217       ENDDO
    1218 
    1219 c
    1220 c Calculer les coefficients turbulents dans l'atmosphere
    1221 c
    1222       DO i = 1, knon
    1223          itop(i) = isommet
    1224       ENDDO
    1225 
    1226       DO k = 2, isommet
    1227       DO i = 1, knon
    1228             zdu2=MAX(cepdu2,(u(i,k)-u(i,k-1))**2
    1229      .                     +(v(i,k)-v(i,k-1))**2)
    1230             zmgeom(i)=zgeop(i,k)-zgeop(i,k-1)
    1231             zdphi =zmgeom(i) / 2.0
    1232             zt = (t(i,k)+t(i,k-1)) * 0.5
    1233             zq = (q(i,k)+q(i,k-1)) * 0.5
    1234 c
    1235 c           calculer Qs et dQs/dT:
    1236 c
    1237             IF (thermcep) THEN
    1238               zdelta = MAX(0.,SIGN(1.,RTT-zt))
    1239               zcvm5 = R5LES*RLVTT/RCPD/(1.0+RVTMP2*zq)*(1.-zdelta)
    1240      .            + R5IES*RLSTT/RCPD/(1.0+RVTMP2*zq)*zdelta
    1241               zqs = R2ES * FOEEW(zt,zdelta) / pplay(i,k)
    1242               zqs = MIN(0.5,zqs)
    1243               zcor = 1./(1.-RETV*zqs)
    1244               zqs = zqs*zcor
    1245               zdqs = FOEDE(zt,zdelta,zcvm5,zqs,zcor)
    1246             ELSE
    1247               IF (zt .LT. t_coup) THEN
    1248                  zqs = qsats(zt) / pplay(i,k)
    1249                  zdqs = dqsats(zt,zqs)
    1250               ELSE
    1251                  zqs = qsatl(zt) / pplay(i,k)
    1252                  zdqs = dqsatl(zt,zqs)
    1253               ENDIF
    1254             ENDIF
    1255 c
    1256 c           calculer la fraction nuageuse (processus humide):
    1257 c
    1258             zfr = (zq+ratqs*zq-zqs) / (2.0*ratqs*zq)
    1259             zfr = MAX(0.0,MIN(1.0,zfr))
    1260             IF (.NOT.richum) zfr = 0.0
    1261 c
    1262 c           calculer le nombre de Richardson:
    1263 c
    1264             IF (tvirtu) THEN
    1265             ztvd =( t(i,k)
    1266      .             + zdphi/RCPD/(1.+RVTMP2*zq)
    1267      .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
    1268      .            )*(1.+RETV*q(i,k))
    1269             ztvu =( t(i,k-1)
    1270      .             - zdphi/RCPD/(1.+RVTMP2*zq)
    1271      .              *( (1.-zfr) + zfr*(1.+RLVTT*zqs/RD/zt)/(1.+zdqs) )
    1272      .            )*(1.+RETV*q(i,k-1))
    1273             zri(i) =zmgeom(i)*(ztvd-ztvu)/(zdu2*0.5*(ztvd+ztvu))
    1274             zri(i) = zri(i)
    1275      .             + zmgeom(i)*zmgeom(i)/RG*gamt(k)
    1276      .               *(paprs(i,k)/101325.0)**RKAPPA
    1277      .               /(zdu2*0.5*(ztvd+ztvu))
    1278 c
    1279             ELSE ! calcul de Ridchardson compatible LMD5
    1280 c
    1281             zri(i) =(RCPD*(t(i,k)-t(i,k-1))
    1282      .              -RD*0.5*(t(i,k)+t(i,k-1))/paprs(i,k)
    1283      .               *(pplay(i,k)-pplay(i,k-1))
    1284      .              )*zmgeom(i)/(zdu2*0.5*RCPD*(t(i,k-1)+t(i,k)))
    1285             zri(i) = zri(i) +
    1286      .             zmgeom(i)*zmgeom(i)*gamt(k)/RG
    1287 cSB     .             /(paprs(i,k)/101325.0)**RKAPPA
    1288      .             *(paprs(i,k)/101325.0)**RKAPPA
    1289      .             /(zdu2*0.5*(t(i,k-1)+t(i,k)))
    1290             ENDIF
    1291 c
    1292 c           finalement, les coefficients d'echange sont obtenus:
    1293 c
    1294             zcdn=SQRT(zdu2) / zmgeom(i) * RG
    1295 c
    1296           IF (opt_ec) THEN
    1297             z2geomf=zgeop(i,k-1)+zgeop(i,k)
    1298             zalm2=(0.5*ckap/RG*z2geomf
    1299      .             /(1.+0.5*ckap/rg/clam*z2geomf))**2
    1300             zalh2=(0.5*ckap/rg*z2geomf
    1301      .             /(1.+0.5*ckap/RG/(clam*SQRT(1.5*cd))*z2geomf))**2
    1302             IF (zri(i).LT.0.0) THEN  ! situation instable
    1303                zscf = ((zgeop(i,k)/zgeop(i,k-1))**(1./3.)-1.)**3
    1304      .                / (zmgeom(i)/RG)**3 / (zgeop(i,k-1)/RG)
    1305                zscf = SQRT(-zri(i)*zscf)
    1306                zscfm = 1.0 / (1.0+3.0*cb*cc*zalm2*zscf)
    1307                zscfh = 1.0 / (1.0+3.0*cb*cc*zalh2*zscf)
    1308                pcfm(i,k)=zcdn*zalm2*(1.-2.0*cb*zri(i)*zscfm)
    1309                pcfh(i,k)=zcdn*zalh2*(1.-3.0*cb*zri(i)*zscfh)
    1310             ELSE ! situation stable
    1311                zscf=SQRT(1.+cd*zri(i))
    1312                pcfm(i,k)=zcdn*zalm2/(1.+2.0*cb*zri(i)/zscf)
    1313                pcfh(i,k)=zcdn*zalh2/(1.+3.0*cb*zri(i)*zscf)
    1314             ENDIF
    1315           ELSE
    1316             zl2(i)=(mixlen*MAX(0.0,(paprs(i,k)-paprs(i,itop(i)+1))
    1317      .                          /(paprs(i,2)-paprs(i,itop(i)+1)) ))**2
    1318             pcfm(i,k)=sqrt(max(zcdn*zcdn*(ric-zri(i))/ric, kstable))
    1319             pcfm(i,k)= zl2(i)* pcfm(i,k)
    1320             pcfh(i,k) = pcfm(i,k) /prandtl ! h et m different
    1321           ENDIF
    1322       ENDDO
    1323       ENDDO
    1324 c
    1325 c Au-dela du sommet, pas de diffusion turbulente:
    1326 c
    1327       DO i = 1, knon
    1328          IF (itop(i)+1 .LE. klev) THEN
    1329             DO k = itop(i)+1, klev
    1330                pcfh(i,k) = 0.0
    1331                pcfm(i,k) = 0.0
    1332             ENDDO
    1333          ENDIF
    1334       ENDDO
    1335 c
     1386      END DO
    13361387      RETURN
    13371388      END
     1389
    13381390      SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t,
    13391391     .                  pcfm, pcfh)
Note: See TracChangeset for help on using the changeset viewer.