Changeset 336 for LMDZ.3.3/branches
- Timestamp:
- Feb 15, 2002, 1:16:35 PM (23 years ago)
- File:
-
- 1 edited
-
LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
r295 r336 1176 1176 c Calculer le frottement au sol (Cdrag) 1177 1177 c 1178 CALL clcdrag(knon, nsrf, is_oce, zxli, u, v, t, q, zgeop, 1179 . ts, qsurf, rugos, pcfm, pcfh, zcdn, zri) 1180 c 1181 c Calculer les coefficients turbulents dans l'atmosphere 1182 c 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 1195 c 1196 c calculer Qs et dQs/dT: 1197 c 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 1216 c 1217 c calculer la fraction nuageuse (processus humide): 1218 c 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 1222 c 1223 c calculer le nombre de Richardson: 1224 c 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)) 1239 c 1240 ELSE ! calcul de Ridchardson compatible LMD5 1241 c 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 1248 cSB . /(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 1252 c 1253 c finalement, les coefficients d'echange sont obtenus: 1254 c 1255 zcdn=SQRT(zdu2) / zmgeom(i) * RG 1256 c 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 1285 c 1286 c Au-dela du sommet, pas de diffusion turbulente: 1287 c 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 1296 c 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) 1304 c ================================================================= c 1305 c Objet : calcul cdrags pour le moment et les flux chaleur sensible, latente (pcfm,pcfh) 1306 c et du nombre de Richardson zri 1307 c ================================================================= c 1308 IMPLICIT NONE 1309 #include "dimensions.h" 1310 #include "dimphy.h" 1311 #include "YOMCST.h" 1312 #include "indicesol.h" 1313 c 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) 1319 c 1320 REAL zcdn 1321 REAL pcfm(klon,klev), pcfh(klon,klev) 1322 c 1323 c Quelques constantes et options: 1324 c 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) 1331 c 1332 c 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) 1338 c 1339 c 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) 1344 c 1345 c Calculer le frottement au sol (Cdrag) 1346 c 1178 1347 DO i = 1, knon 1179 1348 zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2) 1180 1349 zdphi=zgeop(i,1) 1181 c 1350 c ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1) 1182 1351 ztsolv = ts(i) * (1.0+RETV*qsurf(i)) 1183 1352 ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1))) … … 1215 1384 IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25) 1216 1385 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 1336 1387 RETURN 1337 1388 END 1389 1338 1390 SUBROUTINE coefkz2(nsrf, knon, paprs, pplay,t, 1339 1391 . pcfm, pcfh)
Note: See TracChangeset
for help on using the changeset viewer.
