Changeset 393 for LMDZ.3.3/branches
- Timestamp:
- Jul 19, 2002, 11:51:34 AM (22 years ago)
- Location:
- LMDZ.3.3/branches/rel-LF/libf/phylmd
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ.3.3/branches/rel-LF/libf/phylmd/clmain.F
r389 r393 482 482 DO j = 1, knon 483 483 yrugm(j) = 0.018*ycoefm(j,1) * (yu1(j)**2+yv1(j)**2)/RG 484 $ + 0.000014 / sqrt(ycoefm(j,1) * (yu1(j)**2+yv1(j)**2))485 484 yrugm(j) = MAX(1.5e-05,yrugm(j)) 486 485 ENDDO … … 710 709 real emis_new(klon), z0_new(klon) 711 710 real pctsrf_new(klon,nbsrf) 711 c JLD 712 real zzpk 712 713 713 714 c … … 773 774 zx_dq(i,klev) = zx_coef(i,klev) / zx_buf1(i) 774 775 c 775 zx_buf2(i) = delp(i,klev) + zx_coef(i,klev) 776 zx_ch(i,klev) = (local_h(i,klev)*delp(i,klev) 776 zzpk=(pplay(i,klev)/psref(i))**RKAPPA 777 zx_buf2(i) = zzpk*delp(i,klev) + zx_coef(i,klev) 778 zx_ch(i,klev) = (local_h(i,klev)*zzpk*delp(i,klev) 777 779 . -zx_coef(i,klev)*z_gamah(i,klev))/zx_buf2(i) 778 780 zx_dh(i,klev) = zx_coef(i,klev) / zx_buf2(i) … … 788 790 zx_dq(i,k) = zx_coef(i,k) / zx_buf1(i) 789 791 c 790 zx_buf2(i) = delp(i,k)+zx_coef(i,k) 792 zzpk=(pplay(i,k)/psref(i))**RKAPPA 793 zx_buf2(i) = zzpk*delp(i,k)+zx_coef(i,k) 791 794 . +zx_coef(i,k+1)*(1.-zx_dh(i,k+1)) 792 zx_ch(i,k) = (local_h(i,k)* delp(i,k)795 zx_ch(i,k) = (local_h(i,k)*zzpk*delp(i,k) 793 796 . +zx_coef(i,k+1)*zx_ch(i,k+1) 794 797 . +zx_coef(i,k+1)*z_gamah(i,k+1) … … 810 813 zx_dq(i,1) = -1. * RG / zx_buf1(i) 811 814 c 812 zx_buf2(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2)) 813 zx_ch(i,1) = (local_h(i,1)*delp(i,1) 815 zzpk=(pplay(i,1)/psref(i))**RKAPPA 816 zx_buf2(i) = zzpk*delp(i,1) + zx_coef(i,2)*(1.-zx_dh(i,2)) 817 zx_ch(i,1) = (local_h(i,1)*zzpk*delp(i,1) 814 818 . +zx_coef(i,2)*(z_gamah(i,2)+zx_ch(i,2))) 815 819 . /zx_buf2(i) … … 1101 1105 REAL zcfm1(klon), zcfm2(klon) 1102 1106 REAL zcfh1(klon), zcfh2(klon) 1103 c$$$ REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn 1104 c$$$cPB differenciation coefficient de frottement drag et flux chaleur 1105 REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn,zcdh, rugh 1107 REAL zdphi, zdu2, ztvd, ztvu, ztsolv, zcdn 1106 1108 REAL zscf, zucf, zcr 1107 1109 REAL zt, zq, zdelta, zcvm5, zcor, zqs, zfr, zdqs … … 1186 1188 c Calculer le frottement au sol (Cdrag) 1187 1189 c 1188 c$$$c PB 1189 c$$$c essais d'itération pour l'océan 1190 c 1191 rugh = 1.3 e-4 1192 IF (nsrf.EQ.is_oce) THEN 1193 DO k=1,10 1194 c$$$ WRITE(*,*) 'k',k 1195 c$$$ WRITE(*,*) rugos(100) 1196 DO i = 1, knon 1197 zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2) 1198 zdphi=zgeop(i,1) 1199 ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1) 1200 c$$$ ztsolv = ts(i) * (1.0+RETV*qsurf(i)) 1201 ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1))) 1202 . *(1.+RETV*q(i,1)) 1203 zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd) 1204 zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2 1205 c PB ajout drag neutre pour flux chaleur 1206 zcdh = ckap**2/(log(1.+zgeop(i,1)/(RG*rugos(i))) 1207 $ * log(1.+zgeop(i,1)/(RG*rugh))) 1208 IF (zri(i) .ge. 0.) THEN ! situation stable 1209 IF (.NOT.zxli) THEN 1210 zscf=SQRT(1.+cd*ABS(zri(i))) 1211 FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1) 1212 ! zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/ zscf) 1213 zcfm1(i) = zcdn * FRIV 1214 FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 ) 1215 ! zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf) 1216 c PB avec drag neutre pour flux chaleur different 1217 c zcfh1(i) = zcdn * FRIH 1218 zcfh1(i) = zcdh * FRIH 1219 pcfm(i,1) = zcfm1(i) 1220 pcfh(i,1) = zcfh1(i) 1221 ELSE 1222 pcfm(i,1) = zcdn* fsta(zri(i)) 1223 pcfh(i,1) = zcdn* fsta(zri(i)) 1224 ENDIF 1225 ELSE ! situation instable 1226 IF (.NOT.zxli) THEN 1227 zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i)) 1228 . *(1.0+zgeop(i,1)/(RG*rugos(i))))) 1229 zcfm2(i) = zcdn*amax1((1.-2.0*cb*zri(i)*zucf),0.1) 1230 c PB ajout pour prendre drag neutre des flux chaleur different drag neutre vent 1231 zucf=1./(1.+3.0*cb*cc*zcdh*SQRT(ABS(zri(i)) 1232 . *(1.0+zgeop(i,1)/(RG*rugh)))) 1233 c$$$ zcfh2(i) = zcdn*amax1((1.-3.0*cb*zri(i)*zucf),0.1) 1234 zcfh2(i) = zcdh*amax1((1.-3.0*cb*zri(i)*zucf),0.1) 1235 pcfm(i,1) = zcfm2(i) 1236 pcfh(i,1) = zcfh2(i) 1237 ELSE 1238 pcfm(i,1) = zcdn* fins(zri(i)) 1239 pcfh(i,1) = zcdn* fins(zri(i)) 1240 ENDIF 1241 zcr = (0.0016/(zcdh*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) 1242 IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdh*(1.0+zcr**1.25)**(1./1.25) 1243 ENDIF 1244 C 1245 c$$$C PB test drag 1246 c$$$ pcfm(i,1)=zcdn 1247 c$$$ pcfh(i,1) = zcdn 1248 rugos(i)= 0.018*pcfm(i,1) * zdu2/RG 1249 $ +0.11*0.000014/sqrt(pcfm(i,1) * zdu2) 1250 rugh = 0.62*0.000014/sqrt(pcfm(i,1) * zdu2)+1.4e-5 1251 END DO 1252 END DO 1253 1254 ELSE 1255 DO i = 1, knon 1256 zdu2=max(cepdu2,u(i,1)**2+v(i,1)**2) 1257 zdphi=zgeop(i,1) 1258 ztsolv = ts(i) * (1.0+RETV*q(i,1)) ! qsol approx = q(i,1) 1259 c ztsolv = ts(i) * (1.0+RETV*qsurf(i)) 1260 ztvd=(t(i,1)+zdphi/RCPD/(1.+RVTMP2*q(i,1))) 1261 . *(1.+RETV*q(i,1)) 1262 zri(i)=zgeop(i,1)*(ztvd-ztsolv)/(zdu2*ztvd) 1263 zcdn = (ckap/log(1.+zgeop(i,1)/(RG*rugos(i))))**2 1264 c pb ajout pour avoir drage neutre flux chaleur differents 1265 zcdh = ckap**2/(log(1.+zgeop(i,1)/(RG*rugos(i))) 1266 $ * log(1.+zgeop(i,1)/(RG*1.3e-4))) 1267 IF (zri(i) .ge. 0.) THEN ! situation stable 1268 IF (.NOT.zxli) THEN 1269 zscf=SQRT(1.+cd*ABS(zri(i))) 1270 FRIV = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), 0.1) 1271 ! zcfm1(i) = zcdn/(1.+2.0*cb*zri(i)/ zscf) 1272 zcfm1(i) = zcdn * FRIV 1273 FRIH = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), 0.1 ) 1274 ! zcfh1(i) = zcdn/(1.+3.0*cb*zri(i)*zscf) 1275 c$$$C Modif pour drag neutre flux chaleur differents 1276 c$$$ zcfh1(i) = zcdn * FRIH 1277 zcfh1(i) = zcdh * FRIH 1278 pcfm(i,1) = zcfm1(i) 1279 pcfh(i,1) = zcfh1(i) 1280 ELSE 1281 pcfm(i,1) = zcdn* fsta(zri(i)) 1282 pcfh(i,1) = zcdn* fsta(zri(i)) 1283 ENDIF 1284 ELSE ! situation instable 1285 IF (.NOT.zxli) THEN 1286 zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i)) 1287 . *(1.0+zgeop(i,1)/(RG*rugos(i))))) 1288 zcfm2(i) = zcdn*amax1((1.-2.0*cb*zri(i)*zucf),0.1) 1289 C PB ajout pour drage neutre flux chaleur differents 1290 zucf=1./(1.+3.0*cb*cc*zcdn*SQRT(ABS(zri(i)) 1291 . *(1.0+zgeop(i,1)/(RG*rugos(i))))) 1292 c$$$ zcfh2(i) = zcdn*amax1((1.-3.0*cb*zri(i)*zucf),0.1) 1293 zcfh2(i) = zcdh*amax1((1.-3.0*cb*zri(i)*zucf),0.1) 1294 pcfm(i,1) = zcfm2(i) 1295 pcfh(i,1) = zcfh2(i) 1296 ELSE 1297 pcfm(i,1) = zcdn* fins(zri(i)) 1298 pcfh(i,1) = zcdn* fins(zri(i)) 1299 ENDIF 1300 zcr = (0.0016/(zcdn*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) 1301 IF(nsrf.EQ.is_oce)pcfh(i,1)=zcdn*(1.0+zcr**1.25)**(1./1.25) 1302 ENDIF 1303 C 1304 c$$$C PB test drag 1305 c$$$ pcfm(i,1)=zcdn 1306 c$$$ pcfh(i,1) = zcdn 1307 END DO 1308 c$$$CPB fin test iterations 1309 ENDIF 1190 CALL clcdrag(knon, nsrf, zxli, u, v, t, q, zgeop, 1191 . ts, qsurf, rugos, pcfm, pcfh, zcdn, zri) 1310 1192 c 1311 1193 c Calculer les coefficients turbulents dans l'atmosphere -
LMDZ.3.3/branches/rel-LF/libf/phylmd/fisrtilp.F
r383 r393 104 104 REAL zprec_cond(klon) 105 105 cAA 106 REAL zmair, zcpair, zcpeau 107 C Pour la conversion eau-neige 108 REAL zlh_solid(klon), zm_solid 106 109 c--------------------------------------------------------------- 107 110 c … … 209 212 ENDDO 210 213 c 214 c Calculer la varition de temp. de l'air du a la chaleur sensible 215 C transporter par la pluie. 216 C Il resterait a rajouter cet effet de la chaleur sensible sur les 217 C flux de surface, du a la diff. de temp. entre le 1er niveau et la 218 C surface. 219 C 220 DO i = 1, klon 221 zmair=(paprs(i,k)-paprs(i,k+1))/RG 222 zcpair=RCPD*(1.0+RVTMP2*zq(i)) 223 zcpeau=RCPD*RVTMP2 224 zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau 225 $ + zmair*zcpair*zt(i) ) 226 $ / (zmair*zcpair + zrfl(i)*dtime*zcpeau) 227 CC WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1)) 228 ENDDO 229 c 230 c 211 231 c Calculer l'evaporation de la precipitation 212 232 c 233 234 213 235 IF (evap_prec) THEN 214 236 DO i = 1, klon … … 363 385 DO i = 1, klon 364 386 zq(i) = zq(i) - zcond(i) 365 zt(i) = zt(i) + zcond(i) * RLVTT/RCPD 387 c zt(i) = zt(i) + zcond(i) * RLVTT/RCPD 388 zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) 366 389 ENDDO 367 390 c … … 484 507 IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN 485 508 snow(i) = zrfl(i) 509 zlh_solid(i) = RLSTT-RLVTT 486 510 ELSE 487 511 rain(i) = zrfl(i) 488 ENDIF 489 ENDDO 512 zlh_solid(i) = 0. 513 ENDIF 514 ENDDO 515 C 516 C For energy conservation : when snow is present, the solification 517 c latent heat is considered. 518 DO k = 1, klev 519 DO i = 1, klon 520 zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k))) 521 zmair=(paprs(i,k)-paprs(i,k+1))/RG 522 zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime 523 d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair) 524 END DO 525 END DO 490 526 c 491 527 RETURN -
LMDZ.3.3/branches/rel-LF/libf/phylmd/radlwsw.F
r367 r393 39 39 #include "dimphy.h" 40 40 #include "raddim.h" 41 #include "YOETHF.h" 41 42 c 42 43 real rmu0(klon), fract(klon), dist … … 88 89 REAL*8 ztopsw0(kdlon), ztoplw0(kdlon) 89 90 REAL*8 zsolsw0(kdlon), zsollw0(kdlon) 90 91 REAL*8 zznormcp 91 92 c 92 93 c------------------------------------------- … … 215 216 ENDDO 216 217 DO k = 1, kflev 218 c DO i = 1, kdlon 219 c heat(iof+i,k) = zheat(i,k) 220 c cool(iof+i,k) = zcool(i,k) 221 c heat0(iof+i,k) = zheat0(i,k) 222 c cool0(iof+i,k) = zcool0(i,k) 223 c ENDDO 217 224 DO i = 1, kdlon 218 heat(iof+i,k) = zheat(i,k) 219 cool(iof+i,k) = zcool(i,k) 220 heat0(iof+i,k) = zheat0(i,k) 221 cool0(iof+i,k) = zcool0(i,k) 225 C scale factor to take into account the difference between 226 C dry air and watter vapour scpecific heat capacity 227 zznormcp=1.0+RVTMP2*PWV(i,k) 228 heat(iof+i,k) = zheat(i,k)/zznormcp 229 cool(iof+i,k) = zcool(i,k)/zznormcp 230 heat0(iof+i,k) = zheat0(i,k)/zznormcp 231 cool0(iof+i,k) = zcool0(i,k)/zznormcp 222 232 ENDDO 223 233 ENDDO -
LMDZ.3.3/branches/rel-LF/libf/phylmd/suphec.F
r2 r393 123 123 C --------------------------------------------- 124 124 C 125 RCW= 4218.125 RCW=RCPV 126 126 WRITE(UNIT=6,FMT='('' *** Thermodynamic, liquid ***'')') 127 127 WRITE(UNIT=6,FMT='('' Cw = '',E13.7)') RCW … … 132 132 C -------------------------------------------- 133 133 C 134 RCS= 2106.134 RCS=RCPV 135 135 WRITE(UNIT=6,FMT='('' *** thermodynamic, solid ***'')') 136 136 WRITE(UNIT=6,FMT='('' Cs = '',E13.7)') RCS
Note: See TracChangeset
for help on using the changeset viewer.
