Changeset 3496 for LMDZ6/trunk/libf/phylmd/cv3_routines.F90
- Timestamp:
- May 10, 2019, 12:17:35 PM (5 years ago)
- File:
-
- 1 edited
-
LMDZ6/trunk/libf/phylmd/cv3_routines.F90 (modified) (49 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/cv3_routines.F90
r3492 r3496 35 35 36 36 include "cv3param.h" 37 include "cvflag.h" 37 38 include "conema3.h" 38 39 … … 125 126 tlcrit=-55.0 126 127 CALL getin_p('tlcrit',tlcrit) 128 ejectliq=0. 129 CALL getin_p('ejectliq',ejectliq) 130 ejectice=0. 131 CALL getin_p('ejectice',ejectice) 132 cvflag_prec_eject = .FALSE. 133 CALL getin_p('cvflag_prec_eject',cvflag_prec_eject) 127 134 qsat_depends_on_qt = .FALSE. 128 135 CALL getin_p('qsat_depends_on_qt',qsat_depends_on_qt) 136 adiab_ascent_mass_flux_depends_on_ejectliq = .FALSE. 137 CALL getin_p('adiab_ascent_mass_flux_depends_on_ejectliq',adiab_ascent_mass_flux_depends_on_ejectliq) 129 138 130 139 WRITE (*, *) 't_top_max=', t_top_max … … 172 181 173 182 include "cv3param.h" 183 include "cvflag.h" 174 184 175 185 !inputs: … … 292 302 USE mod_phys_lmdz_transfert_para, ONLY : bcast 293 303 USE add_phys_tend_mod, ONLY: fl_cor_ebil 304 USE print_control_mod, ONLY: prt_level 294 305 IMPLICIT NONE 295 306 … … 519 530 END DO 520 531 ENDIF 532 IF (prt_level .GE. 10) THEN 533 print *,'cv3_feed : iflag(1), pfeed(1), plcl(1), wghti(1,k) ', & 534 iflag(1), pfeed(1), plcl(1), (wghti(1,k),k=1,10) 535 ENDIF 521 536 522 537 ! ------------------------------------------------------------------- … … 1108 1123 tnk, qnk, gznk, hnk, t, q, qs, gz, & 1109 1124 p, ph, h, tv, lv, lf, pbase, buoybase, plcl, & 1110 inb, tp, tvp, clw, hp, ep, sigp, buoy, frac) 1125 inb, tp, tvp, clw, hp, ep, sigp, buoy, & 1126 frac_a, frac_s, qpreca, qta) 1111 1127 USE print_control_mod, ONLY: prt_level 1112 1128 IMPLICIT NONE … … 1156 1172 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: ep, sigp, hp 1157 1173 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: buoy 1158 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: frac ! ice fraction 1174 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: frac_a, frac_s 1175 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qpreca 1176 REAL, DIMENSION (nloc, nd), INTENT (OUT) :: qta 1159 1177 1160 1178 !local variables: … … 1162 1180 REAL smallestreal 1163 1181 REAL tg, qg, dqgdT, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit 1182 REAL :: phinu2p 1164 1183 REAL als 1165 1184 REAL :: qsat_new, snew 1166 1185 REAL, DIMENSION (nloc,nd) :: qi 1186 REAL, DIMENSION (nloc,nd) :: ha ! moist static energy of adiabatic ascents 1187 ! taking into account precip ejection 1188 REAL, DIMENSION (nloc,nd) :: hla ! liquid water static energy of adiabatic ascents 1189 ! taking into account precip ejection 1190 REAL, DIMENSION (nloc,nd) :: qcld ! specific cloud water 1167 1191 REAL, DIMENSION (nloc,nd) :: qhsat ! specific humidity at saturation 1168 1192 REAL, DIMENSION (nloc,nd) :: dqhsatdT ! dqhsat/dT 1193 REAL, DIMENSION (nloc,nd) :: frac ! ice fraction function of envt temperature 1194 REAL, DIMENSION (nloc,nd) :: qps ! specific solid precipitation 1195 REAL, DIMENSION (nloc,nd) :: qpl ! specific liquid precipitation 1169 1196 REAL, DIMENSION (nloc) :: ah0, cape, capem, byp 1170 1197 LOGICAL, DIMENSION (nloc) :: lcape 1171 1198 INTEGER, DIMENSION (nloc) :: iposit 1199 REAL :: denomm1 1172 1200 REAL :: by, defrac, pden, tbis 1173 1201 REAL :: fracg … … 1196 1224 END DO 1197 1225 1226 1198 1227 ! ===================================================================== 1199 1228 ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES … … 1211 1240 qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i) 1212 1241 END DO 1242 ! 1243 ! Ice fraction 1244 ! 1245 IF (cvflag_ice) THEN 1246 DO k = minorig, nl 1247 DO i = 1, ncum 1248 frac(i, k) = (Tx - t(i,k))/(Tx - Tm) 1249 frac(i, k) = min(max(frac(i,k),0.0), 1.0) 1250 END DO 1251 END DO 1252 ! Below cloud base, set ice fraction to cloud base value 1253 DO k = 1, nl 1254 DO i = 1, ncum 1255 IF (k<icb(i)) THEN 1256 frac(i,k) = frac(i,icb(i)) 1257 END IF 1258 END DO 1259 END DO 1260 ELSE 1261 DO k = 1, nl 1262 DO i = 1, ncum 1263 frac(i,k) = 0. 1264 END DO 1265 END DO 1266 ENDIF ! (cvflag_ice) 1267 1213 1268 1214 1269 DO k = minorig, nl 1215 1270 DO i = 1,ncum 1271 ha(i,k) = ah0(i) 1272 hla(i,k) = hnk(i) 1273 qta(i,k) = qnk(i) 1274 qpreca(i,k) = 0. 1275 frac_a(i,k) = 0. 1276 frac_s(i,k) = frac(i,k) 1277 qpl(i,k) = 0. 1278 qps(i,k) = 0. 1216 1279 qhsat(i,k) = qs(i,k) 1280 qcld(i,k) = max(qta(i,k)-qhsat(i,k),0.) 1281 IF (k <= icb(i)+1) THEN 1282 qhsat(i,k) = qnk(i)-clw(i,k) 1283 qcld(i,k) = clw(i,k) 1284 ENDIF 1217 1285 ENDDO 1218 1286 ENDDO … … 1347 1415 !---------------------------------------------------------------------------- 1348 1416 ! 1349 ! Ice fraction1350 !1351 DO k = minorig + 1, nl1352 DO i = 1, ncum1353 frac(i, k) = (Tx - t(i,k))/(Tx - Tm)1354 frac(i, k) = min(max(frac(i,k),0.0), 1.0)1355 END DO1356 END DO1357 ! Below cloud base, set ice fraction to cloud base value1358 DO k = 1, nl1359 DO i = 1, ncum1360 IF (k<icb(i)) THEN1361 frac(i,k) = frac(i,icb(i))1362 END IF1363 END DO1364 END DO1365 1366 1417 DO k = minorig + 1, nl 1367 1418 DO i = 1,ncum … … 1376 1427 IF (k>=(icbs(i)+1)) THEN ! convect3 1377 1428 tg = tp(i, k) 1378 IF (tg .gt. Tx ) THEN1429 IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN 1379 1430 es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15)) 1380 1431 qg = eps*es/(p(i,k)-es*(1.-eps)) … … 1386 1437 ENDIF 1387 1438 IF (qsat_depends_on_qt) THEN 1388 dqgdT = dqgdT*(1.-q nk(i))/(1.-qg)**21389 qg = qg*(1.-q nk(i))/(1.-qg)1439 dqgdT = dqgdT*(1.-qta(i,k-1))/(1.-qg)**2 1440 qg = qg*(1.-qta(i,k-1))/(1.-qg) 1390 1441 ENDIF 1391 ahg = (cpd + (cl-cpd)*q nk(i))*tg + lv(i,k)*qg - &1392 lf(i,k)*frac(i,k)*(q nk(i) - qg) + gz(i,k)1442 ahg = (cpd + (cl-cpd)*qta(i,k-1))*tg + lv(i,k)*qg - & 1443 lf(i,k)*frac(i,k)*(qta(i,k-1) - qg) + gz(i,k) 1393 1444 Tp(i,k) = tg + (ah0(i) - ahg)/ & 1394 (cpd + (cl-cpd)*q nk(i) + (lv(i,k)+frac(i,k)*lf(i,k))*dqgdT)1445 (cpd + (cl-cpd)*qta(i,k-1) + (lv(i,k)+frac(i,k)*lf(i,k))*dqgdT) 1395 1446 !! print *,'undilute2 iterations k, Tp(i,k), ah0(i), ahg ', & 1396 1447 !! k, Tp(i,k), ah0(i), ahg … … 1401 1452 IF (k>=(icbs(i)+1)) THEN ! convect3 1402 1453 tg = tp(i, k) 1403 IF (tg .gt. Tx ) THEN1454 IF (tg .gt. Tx .OR. .NOT.cvflag_ice) THEN 1404 1455 es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15)) 1405 1456 qg = eps*es/(p(i,k)-es*(1.-eps)) … … 1408 1459 qg = eps*esi/(p(i,k)-esi*(1.-eps)) 1409 1460 ENDIF 1410 IF (qsat_depends_on_qt) THEN 1411 qg = qg*(1.-qnk(i))/(1.-qg) 1412 ENDIF 1413 clw(i, k) = qnk(i) - qg 1461 IF (qsat_depends_on_qt) THEN 1462 qg = qg*(1.-qta(i,k-1))/(1.-qg) 1463 ENDIF 1464 qhsat(i,k) = qg 1465 END IF ! (k>=(icbs(i)+1)) 1466 END DO ! i = 1, ncum 1467 DO i = 1, ncum 1468 IF (k>=(icbs(i)+1)) THEN ! convect3 1469 clw(i, k) = qta(i,k-1) - qhsat(i,k) 1414 1470 clw(i, k) = max(0.0, clw(i,k)) 1415 tvp(i, k) = max(0., tp(i,k)*(1.+q g/eps-qnk(i)))1471 tvp(i, k) = max(0., tp(i,k)*(1.+qhsat(i,k)/eps-qta(i,k-1))) 1416 1472 ! print*,tvp(i,k),'tvp' 1417 1473 IF (clw(i,k)<1.E-11) THEN … … 1421 1477 END IF ! (k>=(icbs(i)+1)) 1422 1478 END DO ! i = 1, ncum 1479 ! 1480 IF (cvflag_prec_eject) THEN 1481 DO i = 1, ncum 1482 IF (k>=(icbs(i)+1)) THEN ! convect3 1483 ! Specific precipitation (liquid and solid) and ice content 1484 ! before ejection of precipitation !!jygprl 1485 elacrit = elcrit*min(max(1.-(tp(i,k)-T0)/Tlcrit, 0.), 1.) !!jygprl 1486 !!!! qcld(i,k) = min(clw(i,k), elacrit) !!jygprl 1487 qcld(i,k) = min(clw(i,k), elacrit*(1.-qta(i,k-1))/(1.-elacrit)) !!jygprl 1488 phinu2p = qhsat(i,k-1) + qcld(i,k-1) - (qhsat(i,k) + qcld(i,k)) !!jygprl 1489 qpl(i,k) = qpl(i,k-1) + (1.-frac(i,k))*phinu2p !!jygprl 1490 qps(i,k) = qps(i,k-1) + frac(i,k) *phinu2p !!jygprl 1491 qi(i,k) = (1.-ejectliq)*clw(i,k)*frac(i,k) + & !!jygprl 1492 ejectliq*(qps(i,k-1) + frac(i,k)*(phinu2p+qcld(i,k))) !!jygprl 1493 !! 1494 ! ===================================================================================== 1495 ! Ejection of precipitation from adiabatic ascents if requested (cvflag_prec_eject=True): 1496 ! Compute the steps of total water (qta), of moist static energy (ha), of specific 1497 ! precipitation (qpl and qps) and of specific cloud water (qcld) associated with precipitation 1498 ! ejection. 1499 ! ===================================================================================== 1500 ! 1501 ! Verif 1502 qpreca(i,k) = ejectliq*qpl(i,k) + ejectice*qps(i,k) !!jygprl 1503 frac_a(i,k) = ejectice*qps(i,k)/max(qpreca(i,k),smallestreal) !!jygprl 1504 frac_s(i,k) = (1.-ejectliq)*frac(i,k) + & !!jygprl 1505 ejectliq*(1. - (qpl(i,k)+(1.-frac(i,k))*qcld(i,k))/max(clw(i,k),smallestreal)) !!jygprl 1506 ! 1507 denomm1 = 1./(1. - qpreca(i,k)) 1508 ! 1509 qta(i,k) = qta(i,k-1) - & 1510 qpreca(i,k)*(1.-qta(i,k-1))*denomm1 1511 ha(i,k) = ha(i,k-1) + & 1512 ( qpreca(i,k)*(-(1.-qta(i,k-1))*(cl-cpd)*tp(i,k) + & 1513 lv(i,k)*qhsat(i,k) - lf(i,k)*(frac_s(i,k)*qcld(i,k)+qps(i,k))) + & 1514 lf(i,k)*ejectice*qps(i,k))*denomm1 1515 hla(i,k) = hla(i,k-1) + & 1516 ( qpreca(i,k)*(-(1.-qta(i,k-1))*(cpv-cpd)*tp(i,k) - & 1517 lv(i,k)*((1.-frac_s(i,k))*qcld(i,k)+qpl(i,k)) - & 1518 (lv(i,k)+lf(i,k))*(frac_s(i,k)*qcld(i,k)+qps(i,k))) + & 1519 lv(i,k)*ejectliq*qpl(i,k) + (lv(i,k)+lf(i,k))*ejectice*qps(i,k))*denomm1 1520 qpl(i,k) = qpl(i,k)*(1.-ejectliq)*denomm1 1521 qps(i,k) = qps(i,k)*(1.-ejectice)*denomm1 1522 qcld(i,k) = qcld(i,k)*denomm1 1523 qhsat(i,k) = qhsat(i,k)*(1.-qta(i,k))/(1.-qta(i,k-1)) 1524 END IF ! (k>=(icbs(i)+1)) 1525 END DO ! i = 1, ncum 1526 ENDIF ! (cvflag_prec_eject) 1527 ! 1423 1528 END DO ! k = minorig + 1, nl 1424 1529 ! 1425 1530 !---------------------------------------------------------------------------- 1426 1531 ! 1427 ELSE ! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)1532 ELSE IF (icvflag_Tpa == 0) THEN! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1) 1428 1533 ! 1429 1534 !---------------------------------------------------------------------------- … … 1587 1692 !---------------------------------------------------------------------------- 1588 1693 ! 1589 ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE 1694 ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE (icvflag_Tpa == 0) 1590 1695 ! 1591 1696 !---------------------------------------------------------------------------- … … 1872 1977 IF (cvflag_ice) THEN 1873 1978 ! 1979 IF (cvflag_prec_eject) THEN 1980 !! DO k = minorig + 1, nl 1981 !! DO i = 1, ncum 1982 !! IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN 1983 !! frac_s(i,k) = qi(i,k)/max(clw(i,k),smallestreal) 1984 !! frac_s(i,k) = 1. - (qpl(i,k)+(1.-frac_s(i,k))*qcld(i,k))/max(clw(i,k),smallestreal) 1985 !! END IF 1986 !! END DO 1987 !! END DO 1988 ELSE ! (cvflag_prec_eject) 1874 1989 DO k = minorig + 1, nl 1875 1990 DO i = 1, ncum … … 1877 1992 !jyg< frac computation moved to beginning of cv3_undilute2. 1878 1993 ! kept here for compatibility test with CMip6 version 1879 frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15) 1880 frac(i, k) = min(max(frac(i,k),0.0), 1.0) 1881 hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac(i,k)*lf(i,k))* & 1882 ep(i, k)*clw(i, k) 1994 frac_s(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15) 1995 frac_s(i, k) = min(max(frac_s(i,k),0.0), 1.0) 1883 1996 END IF 1884 1997 END DO 1885 1998 END DO 1886 ! Below cloud base, set ice fraction to cloud base value 1887 !! DO k = 1, nl 1888 !! DO i = 1, ncum 1889 !! IF (k<icb(i)) THEN 1890 !! frac(i,k) = frac(i,icb(i)) 1891 !! END IF 1892 !! END DO 1893 !! END DO 1999 ENDIF ! (cvflag_prec_eject) ELSE 2000 DO k = minorig + 1, nl 2001 DO i = 1, ncum 2002 IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN 2003 !! hp(i, k) = hnk(i) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac_s(i,k)*lf(i,k))* & !!jygprl 2004 !! ep(i, k)*clw(i, k) !!jygprl 2005 hp(i, k) = hla(i,k-1) + (lv(i,k)+(cpd-cpv)*t(i,k)+frac_s(i,k)*lf(i,k))* & !!jygprl 2006 ep(i, k)*clw(i, k) !!jygprl 2007 END IF 2008 END DO 2009 END DO 1894 2010 ! 1895 2011 ELSE ! (cvflag_ice) … … 2576 2692 SUBROUTINE cv3_unsat(nloc, ncum, nd, na, ntra, icb, inb, iflag, & 2577 2693 t, rr, rs, gz, u, v, tra, p, ph, & 2578 th, tv, lv, lf, cpn, ep, sigp, clw, &2694 th, tv, lv, lf, cpn, ep, sigp, clw, frac_s, qpreca, frac_a, qta , & !!jygprl 2579 2695 m, ment, elij, delt, plcl, coef_clos, & 2580 2696 mp, rp, up, vp, trap, wt, water, evap, fondue, ice, & 2581 2697 faci, b, sigd, & 2582 wdtrainA, wdtrain M) ! RomP2698 wdtrainA, wdtrainS, wdtrainM) ! RomP 2583 2699 USE print_control_mod, ONLY: prt_level, lunout 2584 2700 IMPLICIT NONE … … 2598 2714 REAL, DIMENSION (nloc, na), INTENT (IN) :: gz 2599 2715 REAL, DIMENSION (nloc, nd), INTENT (IN) :: u, v 2600 REAL tra(nloc, nd, ntra) 2601 REAL p(nloc, nd), ph(nloc, nd+1) 2602 REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, sigp, clw 2716 REAL, DIMENSION (nloc, nd, ntra), INTENT(IN) :: tra 2717 REAL, DIMENSION (nloc, nd), INTENT (IN) :: p 2718 REAL, DIMENSION (nloc, nd+1), INTENT (IN) :: ph 2719 REAL, DIMENSION (nloc, na), INTENT (IN) :: ep, sigp, clw !adiab ascent shedding 2720 REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_s !ice fraction in adiab ascent shedding !!jygprl 2721 REAL, DIMENSION (nloc, na), INTENT (IN) :: qpreca !adiab ascent precip !!jygprl 2722 REAL, DIMENSION (nloc, na), INTENT (IN) :: frac_a !ice fraction in adiab ascent precip !!jygprl 2723 REAL, DIMENSION (nloc, na), INTENT (IN) :: qta !adiab ascent specific total water !!jygprl 2603 2724 REAL, DIMENSION (nloc, na), INTENT (IN) :: th, tv, lv, cpn 2604 2725 REAL, DIMENSION (nloc, na), INTENT (IN) :: lf … … 2613 2734 REAL, DIMENSION (nloc, na), INTENT (OUT) :: mp, rp, up, vp 2614 2735 REAL, DIMENSION (nloc, na), INTENT (OUT) :: water, evap, wt 2615 REAL, DIMENSION (nloc, na), INTENT (OUT) :: ice, fondue, faci 2736 REAL, DIMENSION (nloc, na), INTENT (OUT) :: ice, fondue 2737 REAL, DIMENSION (nloc, na), INTENT (OUT) :: faci ! ice fraction in precipitation 2616 2738 REAL, DIMENSION (nloc, na, ntra), INTENT (OUT) :: trap 2617 2739 REAL, DIMENSION (nloc, na), INTENT (OUT) :: b … … 2621 2743 ! Distinction des wdtrain 2622 2744 ! Pa = wdtrainA Pm = wdtrainM 2623 REAL, DIMENSION (nloc, na), INTENT (OUT) :: wdtrainA, wdtrain M2745 REAL, DIMENSION (nloc, na), INTENT (OUT) :: wdtrainA, wdtrainS, wdtrainM 2624 2746 2625 2747 !local variables 2626 2748 INTEGER i, j, k, il, num1, ndp1 2749 REAL smallestreal 2627 2750 REAL tinv, delti, coef 2628 2751 REAL awat, afac, afac1, afac2, bfac … … 2631 2754 REAL ampmax, thaw 2632 2755 REAL tevap(nloc) 2633 REAL lvcp(nloc, na), lfcp(nloc, na) 2634 REAL h(nloc, na), hm(nloc, na) 2635 REAL frac(nloc, na) 2636 REAL fraci(nloc, na), prec(nloc, na) 2756 REAL, DIMENSION (nloc, na) :: lvcp, lfcp 2757 REAL, DIMENSION (nloc, na) :: h, hm 2758 REAL, DIMENSION (nloc, na) :: ma 2759 REAL, DIMENSION (nloc, na) :: frac ! ice fraction in precipitation source 2760 REAL, DIMENSION (nloc, na) :: fraci ! provisionnal ice fraction in precipitation 2761 REAL, DIMENSION (nloc, na) :: prec 2637 2762 REAL wdtrain(nloc) 2638 2763 LOGICAL lwork(nloc), mplus(nloc) … … 2641 2766 ! ------------------------------------------------------ 2642 2767 IF (prt_level .GE. 10) print *,' ->cv3_unsat, iflag(1) ', iflag(1) 2768 2769 smallestreal=tiny(smallestreal) 2643 2770 2644 2771 ! ============================= … … 2660 2787 !! RomP >>> 2661 2788 wdtrainA(:,:) = 0. 2789 wdtrainS(:,:) = 0. 2662 2790 wdtrainM(:,:) = 0. 2663 2791 !! RomP <<< … … 2715 2843 END DO 2716 2844 2845 ! 2846 ! Get adiabatic ascent mass flux 2847 ! 2848 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2849 IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN 2850 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2851 !!! Warning : this option leads to water conservation violation 2852 !!! Expert only 2853 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2854 DO il = 1, ncum 2855 ma(il, nlp) = 0. 2856 ma(il, 1) = 0. 2857 END DO 2858 2859 DO i = nl, 2, -1 2860 DO il = 1, ncum 2861 ma(il, i) = ma(il, i+1)*(1.-qta(il,i))/(1.-qta(il,i-1)) + m(il, i) 2862 END DO 2863 END DO 2864 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2865 ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq) 2866 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2867 DO il = 1, ncum 2868 ma(il, nlp) = 0. 2869 ma(il, 1) = 0. 2870 END DO 2871 2872 DO i = nl, 2, -1 2873 DO il = 1, ncum 2874 ma(il, i) = ma(il, i+1) + m(il, i) 2875 END DO 2876 END DO 2877 2878 ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE 2879 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2717 2880 2718 2881 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ … … 2739 2902 ! *** calculate detrained precipitation *** 2740 2903 2741 DO il = 1, ncum 2742 IF (i<=inb(il) .AND. lwork(il)) THEN 2743 IF (cvflag_grav) THEN 2744 wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i) 2745 wdtrainA(il, i) = wdtrain(il)/grav ! Pa RomP 2746 ELSE 2747 wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i) 2748 wdtrainA(il, i) = wdtrain(il)/10. ! Pa RomP 2749 END IF 2750 END IF 2751 END DO 2904 2905 DO il = 1, ncum 2906 IF (i<=inb(il) .AND. lwork(il)) THEN 2907 wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i) 2908 wdtrainS(il, i) = wdtrain(il)/grav ! Ps jyg 2909 !! wdtrainA(il, i) = wdtrain(il)/grav ! Ps RomP 2910 END IF 2911 END DO 2752 2912 2753 2913 IF (i>1) THEN … … 2757 2917 awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i) 2758 2918 awat = max(awat, 0.0) 2759 IF (cvflag_grav) THEN 2760 wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i) 2761 wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP 2762 ELSE 2763 wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i) 2764 wdtrainM(il, i) = wdtrain(il)/10. - wdtrainA(il, i) ! Pm RomP 2765 END IF 2919 wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i) 2920 wdtrainM(il, i) = wdtrain(il)/grav - wdtrainS(il, i) ! Pm jyg 2921 !! wdtrainM(il, i) = wdtrain(il)/grav - wdtrainA(il, i) ! Pm RomP 2766 2922 END IF 2767 2923 END DO … … 2769 2925 END IF 2770 2926 2927 IF (cvflag_prec_eject) THEN 2928 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2929 IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN 2930 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2931 !!! Warning : this option leads to water conservation violation 2932 !!! Expert only 2933 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2934 IF ( i > 1) THEN 2935 DO il = 1, ncum 2936 IF (i<=inb(il) .AND. lwork(il)) THEN 2937 wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i))/(1. - qta(il, i-1)) ! Pa jygprl 2938 wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i) 2939 END IF 2940 END DO 2941 ENDIF ! ( i > 1) 2942 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2943 ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq) 2944 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2945 IF ( i > 1) THEN 2946 DO il = 1, ncum 2947 IF (i<=inb(il) .AND. lwork(il)) THEN 2948 wdtrainA(il,i) = ma(il, i+1)*(qta(il, i-1)-qta(il,i)) ! Pa jygprl 2949 wdtrain(il) = wdtrain(il) + grav*wdtrainA(il,i) 2950 END IF 2951 END DO 2952 ENDIF ! ( i > 1) 2953 2954 ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE 2955 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2956 ENDIF ! (cvflag_prec_eject) 2957 2771 2958 2772 2959 ! *** find rain water and evaporation using provisional *** … … 2774 2961 2775 2962 2963 IF (cvflag_ice) THEN !!jygprl 2964 DO il = 1, ncum !!jygprl 2965 IF (i<=inb(il) .AND. lwork(il)) THEN !!jygprl 2966 frac(il, i) = (frac_a(il,i)*wdtrainA(il,i)+frac_s(il,i)*(wdtrainS(il,i)+wdtrainM(il,i))) / & !!jygprl 2967 max(wdtrainA(il,i)+wdtrainS(il,i)+wdtrainM(il,i),smallestreal) !!jygprl 2968 fraci(il, i) = frac(il, i) !!jygprl 2969 END IF !!jygprl 2970 END DO !!jygprl 2971 END IF !!jygprl 2972 2776 2973 DO il = 1, ncum 2777 2974 IF (i<=inb(il) .AND. lwork(il)) THEN … … 2779 2976 wt(il, i) = 45.0 2780 2977 2781 IF (cvflag_ice) THEN2782 frac(il, inb(il)) = 1. - (t(il,inb(il))-243.15)/(263.15-243.15)2783 frac(il, inb(il)) = min(max(frac(il,inb(il)),0.), 1.)2784 fraci(il, inb(il)) = frac(il, inb(il))2785 ELSE2786 CONTINUE2787 END IF2788 2978 2789 2979 IF (i<inb(il)) THEN … … 2802 2992 rp(il, i) = 0.5*(rp(il,i)+rr(il,i)) 2803 2993 END IF 2804 fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15)2805 fraci(il, i) = min(max(fraci(il,i),0.0), 1.0)2994 !! fraci(il, i) = 1. - (t(il,i)-243.15)/(263.15-243.15) 2995 !! fraci(il, i) = min(max(fraci(il,i),0.0), 1.0) 2806 2996 rp(il, i) = max(rp(il,i), 0.0) 2807 2997 rp(il, i) = amin1(rp(il,i), rs(il,i)) … … 3230 3420 t, rr, t_wake, rr_wake, s_wake, u, v, tra, & 3231 3421 gz, p, ph, h, hp, lv, lf, cpn, th, th_wake, & 3232 ep, clw, m, tp, mp, rp, up, vp, trap, &3422 ep, clw, qpreca, m, tp, mp, rp, up, vp, trap, & 3233 3423 wt, water, ice, evap, fondue, faci, b, sigd, & 3234 3424 ment, qent, hent, iflag_mix, uent, vent, & … … 3240 3430 !! tls, tps, ! useless . jyg 3241 3431 qcondc, wd, & 3242 ftd, fqd, q nk, qtc, sigt, tau_cld_cv, coefw_cld_cv)3432 ftd, fqd, qta, qtc, sigt, tau_cld_cv, coefw_cld_cv) 3243 3433 3244 3434 USE print_control_mod, ONLY: lunout, prt_level … … 3280 3470 REAL, DIMENSION (nloc, na, na, ntra), INTENT (IN) :: traent 3281 3471 REAL, DIMENSION (nloc, nd), INTENT (IN) :: tv, tvp, wghti 3282 REAL,INTENT(IN) :: tau_cld_cv, coefw_cld_cv 3472 REAL, DIMENSION (nloc, nd), INTENT (IN) :: qta 3473 REAL, DIMENSION (nloc, na),INTENT(IN) :: qpreca 3474 REAL, INTENT(IN) :: tau_cld_cv, coefw_cld_cv 3283 3475 ! 3284 3476 !input/output: … … 3309 3501 REAL :: ax, bx, cx, dx, ex 3310 3502 REAL :: cpinv, rdcp, dpinv 3503 REAL :: sigaq 3311 3504 REAL, DIMENSION (nloc) :: awat 3312 3505 REAL, DIMENSION (nloc, nd) :: lvcp, lfcp ! , mke ! unused . jyg … … 3326 3519 REAL, DIMENSION (nloc) :: sument 3327 3520 REAL, DIMENSION (nloc, nd) :: sigment, qtment ! cld 3328 REAL, DIMENSION (nloc) :: qnk3329 3521 REAL sumdq !jyg 3330 3522 ! … … 3437 3629 END DO 3438 3630 3631 ! - Adiabatic ascent mass flux "ma" and cloud base mass flux "cbmf" 3632 !----------------------------------------------------------------- 3633 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3634 IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN 3635 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3636 !!! Warning : this option leads to water conservation violation 3637 !!! Expert only 3638 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3639 DO il = 1, ncum 3640 ma(il, nlp) = 0. 3641 ma(il, 1) = 0. 3642 END DO 3643 DO k = nl, 2, -1 3644 DO il = 1, ncum 3645 ma(il, k) = ma(il, k+1)*(1.-qta(il, k))/(1.-qta(il, k-1)) + m(il, k) 3646 cbmf(il) = max(cbmf(il), ma(il,k)) 3647 END DO 3648 END DO 3649 DO k = 2,nl 3650 DO il = 1, ncum 3651 IF (k <icb(il)) THEN 3652 ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il) 3653 ENDIF 3654 END DO 3655 END DO 3656 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3657 ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq) 3658 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3659 !! Line kept for compatibility with earlier versions 3439 3660 DO k = 2, nl 3440 3661 DO il = 1, ncum … … 3445 3666 END DO 3446 3667 3668 DO il = 1, ncum 3669 ma(il, nlp) = 0. 3670 ma(il, 1) = 0. 3671 END DO 3672 DO k = nl, 2, -1 3673 DO il = 1, ncum 3674 ma(il, k) = ma(il, k+1) + m(il, k) 3675 END DO 3676 END DO 3677 DO k = 2,nl 3678 DO il = 1, ncum 3679 IF (k <icb(il)) THEN 3680 ma(il, k) = ma(il, k-1) + wghti(il,k-1)*cbmf(il) 3681 ENDIF 3682 END DO 3683 END DO 3684 3685 ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE 3686 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3687 ! 3447 3688 ! print*,'cv3_yield avant ft' 3448 3689 ! am is the part of cbmf taken from the first level … … 3581 3822 !*** Compute convective mass fluxes upwd and dnwd *** 3582 3823 3824 ! 3825 ! ================================================= 3826 ! upward fluxes | 3827 ! ------------------------------------------------ 3828 ! 3583 3829 upwd(:,:) = 0. 3584 3830 up_to(:,:) = 0. 3585 3831 up_from(:,:) = 0. 3586 dnwd(:,:) = 0. 3587 dn_to(:,:) = 0. 3588 dn_from(:,:) = 0. 3589 ! 3590 ! ================================================= 3591 ! upward fluxes | 3592 ! ------------------------------------------------ 3832 ! 3833 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3834 IF (adiab_ascent_mass_flux_depends_on_ejectliq) THEN 3835 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3836 !! The decrease of the adiabatic ascent mass flux due to ejection of precipitation 3837 !! is taken into account. 3838 !! WARNING : in the present version, taking into account the mass-flux decrease due to 3839 !! precipitation ejection leads to water conservation violation. 3840 ! 3841 ! - Upward mass flux of mixed draughts 3842 !--------------------------------------- 3843 DO i = 2, nl 3844 DO j = 1, i-1 3845 DO il = 1, ncum 3846 IF (i<=inb(il)) THEN 3847 up_to(il,i) = up_to(il,i) + ment(il,j,i) 3848 ENDIF 3849 ENDDO 3850 ENDDO 3851 ENDDO 3852 ! 3853 DO j = 3, nl 3854 DO i = 2, j-1 3855 DO il = 1, ncum 3856 IF (j<=inb(il)) THEN 3857 up_from(il,i) = up_from(il,i) + ment(il,i,j) 3858 ENDIF 3859 ENDDO 3860 ENDDO 3861 ENDDO 3862 ! 3863 ! The difference between upwd(il,i) and upwd(il,i-1) is due to updrafts ending in layer 3864 !(i-1) (theses drafts cross interface (i-1) but not interface(i)) and to updrafts starting 3865 !from layer (i-1) (theses drafts cross interface (i) but not interface(i-1)): 3866 ! 3867 DO i = 2, nlp 3868 DO il = 1, ncum 3869 IF (i<=inb(il)+1) THEN 3870 upwd(il,i) = max(0., upwd(il,i-1) - up_to(il,i-1) + up_from(il,i-1)) 3871 ENDIF 3872 ENDDO 3873 ENDDO 3874 ! 3875 ! - Total upward mass flux 3876 !--------------------------- 3877 DO i = 2, nlp 3878 DO il = 1, ncum 3879 IF (i<=inb(il)+1) THEN 3880 upwd(il,i) = upwd(il,i) + ma(il,i) 3881 ENDIF 3882 ENDDO 3883 ENDDO 3884 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3885 ELSE ! (adiab_ascent_mass_flux_depends_on_ejectliq) 3886 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3887 !! The decrease of the adiabatic ascent mass flux due to ejection of precipitation 3888 !! is not taken into account. 3889 ! 3890 ! - Upward mass flux 3891 !------------------- 3593 3892 DO i = 2, nl 3594 3893 DO il = 1, ncum … … 3613 3912 ENDDO 3614 3913 ENDDO 3615 !!DO i = 2, nl 3616 !! DO j = i+1, nl !! Permuter les boucles i et j 3914 ! 3617 3915 DO j = 3, nl 3618 3916 DO i = 2, j-1 … … 3636 3934 ENDDO 3637 3935 ENDDO 3936 3937 3938 ENDIF ! (adiab_ascent_mass_flux_depends_on_ejectliq) ELSE 3939 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3940 3638 3941 ! 3639 3942 ! ================================================= 3640 3943 ! downward fluxes | 3641 3944 ! ------------------------------------------------ 3945 dnwd(:,:) = 0. 3946 dn_to(:,:) = 0. 3947 dn_from(:,:) = 0. 3642 3948 DO i = 1, nl 3643 3949 DO j = i+1, nl … … 3650 3956 ENDDO 3651 3957 ! 3652 !!DO i = 2, nl3653 !! DO j = 1, i-1 !! Permuter les boucles i et j3654 3958 DO j = 1, nl 3655 3959 DO i = j+1, nl … … 4423 4727 !!!! 4424 4728 !!!! ENDDO 4729 4730 !! DO i = 1, nlp 4731 !! DO il = 1, ncum 4732 !! ma(il, i) = 0 4733 !! END DO 4734 !! END DO 4735 !! 4736 !! DO i = 1, nl 4737 !! DO j = i, nl 4738 !! DO il = 1, ncum 4739 !! ma(il, i) = ma(il, i) + m(il, j) 4740 !! END DO 4741 !! END DO 4742 !! END DO 4743 4744 !jyg< (loops stop at nl) 4745 !! DO i = nl + 1, nd 4746 !! DO il = 1, ncum 4747 !! ma(il, i) = 0. 4748 !! END DO 4749 !! END DO 4750 !>jyg 4751 4752 !! DO i = 1, nl 4753 !! DO il = 1, ncum 4754 !! IF (i<=(icb(il)-1)) THEN 4755 !! ma(il, i) = 0 4756 !! END IF 4757 !! END DO 4758 !! END DO 4759 4425 4760 !----------------------------------------------------------- 4426 4761 ENDIF !(.NOT.ok_optim_yield) !| … … 4447 4782 !>jyg 4448 4783 4449 DO i = 1, nlp4450 DO il = 1, ncum4451 ma(il, i) = 04452 END DO4453 END DO4454 4455 DO i = 1, nl4456 DO j = i, nl4457 DO il = 1, ncum4458 ma(il, i) = ma(il, i) + m(il, j)4459 END DO4460 END DO4461 END DO4462 4463 !jyg< (loops stop at nl)4464 !! DO i = nl + 1, nd4465 !! DO il = 1, ncum4466 !! ma(il, i) = 0.4467 !! END DO4468 !! END DO4469 !>jyg4470 4471 DO i = 1, nl4472 DO il = 1, ncum4473 IF (i<=(icb(il)-1)) THEN4474 ma(il, i) = 04475 END IF4476 END DO4477 END DO4478 4784 4479 4785 ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc … … 4562 4868 ! 14/01/15 AJ delta n'a rien à faire là... 4563 4869 DO il = 1, ncum ! cld 4564 IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld 4870 !! IF (wa(il,i)>0.0 .AND. iflag(il)<=1) & ! cld 4871 !! siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld 4872 !! *rrd*tvp(il, i)/p(il, i)/100. ! cld 4873 !! 4874 !! siga(il, i) = min(siga(il,i), 1.0) ! cld 4875 sigaq = 0. 4876 IF (wa(il,i)>0.0 .AND. iflag(il)<=1) THEN ! cld 4565 4877 siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) & ! cld 4566 *rrd*tvp(il, i)/p(il, i)/100. ! cld 4567 4568 siga(il, i) = min(siga(il,i), 1.0) ! cld 4878 *rrd*tvp(il, i)/p(il, i)/100. ! cld 4879 siga(il, i) = min(siga(il,i), 1.0) ! cld 4880 sigaq = siga(il,i)*qta(il,i-1) ! cld 4881 ENDIF 4569 4882 4570 4883 ! IM cf. FH … … 4578 4891 sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1)) ! cld 4579 4892 sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i)) ! cld 4580 qtc(il, i) = (siga(il,i)*qnk(il)+sigment(il,i)*qtment(il,i)) & ! cld 4893 !! qtc(il, i) = (siga(il,i)*qta(il,i-1)+sigment(il,i)*qtment(il,i)) & ! cld 4894 qtc(il, i) = (sigaq+sigment(il,i)*qtment(il,i)) & ! cld 4581 4895 /(siga(il,i)+sigment(il,i)) ! cld 4582 4896 sigt(il,i) = sigment(il, i) + siga(il, i) 4583 4897 4584 ! qtc(il, i) = siga(il,i)*q nk(il)+(1.-siga(il,i))*qtment(il,i) ! cld4898 ! qtc(il, i) = siga(il,i)*qta(il,i-1)+(1.-siga(il,i))*qtment(il,i) ! cld 4585 4899 ! print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i) 4586 4900 … … 4871 5185 do k=1,nl 4872 5186 do i=1,ncum 4873 4874 5187 hp(i,k)=h(i,k) 5188 enddo 4875 5189 enddo 4876 5190
Note: See TracChangeset
for help on using the changeset viewer.
