Ignore:
Timestamp:
May 4, 2019, 12:35:08 PM (6 years ago)
Author:
jyg
Message:

New computations of the lifted parcel in Emanuel's
convective scheme, controlled by two new flags:

Integer flag icvflag_Tpa:

0 -> Ice fraction in the adiabatic ascents =

function of environment tempertaure.
Computation in two steps: the first is ice
free; the second adds ice (Arnaud Jam
scheme)

1 -> Ice fraction in the adiabatic ascents =

function of environment tempertaure.
Computation in one single step.

2 -> Ice fraction in the adiabatic ascents =

function of ascent temperature.
Computation in one single step.

Logical flag qsat_depends_on_qt:

When True, specific hunidity at saturation

is a function of the amount of condsensate.

Else it is equal to qsat.

Location:
LMDZ6/trunk/libf/phylmd
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/cv3_routines.F90

    r3435 r3492  
    125125     tlcrit=-55.0
    126126     CALL getin_p('tlcrit',tlcrit)
     127     qsat_depends_on_qt = .FALSE.
     128     CALL getin_p('qsat_depends_on_qt',qsat_depends_on_qt)
    127129
    128130    WRITE (*, *) 't_top_max=', t_top_max
    … …  
    236238! debug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
    237239      lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15)
    238       lf(i, k) = lf0 - clmci*(t(i,k)-273.15)
     240!!      lf(i, k) = lf0 - clmci*(t(i,k)-273.15)   ! erreur de signe !!
     241      lf(i, k) = lf0 + clmci*(t(i,k)-273.15)
    239242      cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k)
    240243      cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k)
    … …  
    11531156  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ep, sigp, hp
    11541157  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: buoy
    1155   REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: frac
     1158  REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: frac   ! ice fraction
    11561159
    11571160!local variables:
    11581161  INTEGER i, j, k
    1159   REAL tg, qg, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
     1162  REAL smallestreal
     1163  REAL tg, qg, dqgdT, ahg, alv, alf, s, tc, es, esi, denom, rg, tca, elacrit
    11601164  REAL als
    1161   REAL qsat_new, snew, qi(nloc, nd)
    1162   REAL by, defrac, pden, tbis
    1163   REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
    1164   LOGICAL lcape(nloc)
    1165   INTEGER iposit(nloc)
    1166   REAL fracg
    1167   REAL deltap
     1165  REAL                                               :: qsat_new, snew
     1166  REAL, DIMENSION (nloc,nd)                          :: qi
     1167  REAL, DIMENSION (nloc,nd)                          :: qhsat    ! specific humidity at saturation
     1168  REAL, DIMENSION (nloc,nd)                          :: dqhsatdT ! dqhsat/dT
     1169  REAL, DIMENSION (nloc)                             :: ah0, cape, capem, byp
     1170  LOGICAL, DIMENSION (nloc)                          :: lcape
     1171  INTEGER, DIMENSION (nloc)                          :: iposit
     1172  REAL                                               :: by, defrac, pden, tbis
     1173  REAL                                               :: fracg
     1174  REAL                                               :: deltap
     1175  REAL, SAVE                                         :: Tx, Tm
     1176  DATA Tx/263.15/, Tm/243.15/
     1177!$OMP THREADPRIVATE(Tx, Tm)
     1178  REAL                                               :: aa, bb, dd, ddelta, discr
     1179  REAL                                               :: ff, fp
     1180  REAL                                               :: coefx, coefm, Zx, Zm, Ux, U, Um
    11681181
    11691182  IF (prt_level >= 10) THEN
    1170     print *,'cv3_undilute2.0. t(1,k), q(1,k), qs(1,k) ', &
    1171                         (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
     1183    print *,'cv3_undilute2.0. icvflag_Tpa, t(1,k), q(1,k), qs(1,k) ', &
     1184                        icvflag_Tpa, (k, t(1,k), q(1,k), qs(1,k), k = 1,nl)
    11721185  ENDIF
     1186  smallestreal=tiny(smallestreal)
    11731187
    11741188! =====================================================================
    … …  
    11981212  END DO
    11991213
     1214  DO k = minorig, nl
     1215    DO i = 1,ncum
     1216      qhsat(i,k) = qs(i,k)
     1217    ENDDO
     1218  ENDDO
     1219
     1220!jyg<
     1221! =====================================================================
     1222! --- SET THE THE FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
     1223! =====================================================================
     1224  DO k = 1, nl
     1225    DO i = 1, ncum
     1226      ep(i, k) = 0.0
     1227      sigp(i, k) = spfac
     1228    END DO
     1229  END DO
     1230!>jyg
     1231!
    12001232
    12011233! ***  Find lifted parcel quantities above cloud base    ***
    12021234
    1203 
     1235!----------------------------------------------------------------------------
     1236!
     1237  IF (icvflag_Tpa == 2) THEN
     1238!
     1239!----------------------------------------------------------------------------
     1240!
     1241    DO k = minorig + 1, nl
     1242      DO i = 1,ncum
     1243        tp(i,k) = t(i,k)
     1244      ENDDO
     1245!!      alv = lv0 - clmcpv*(t(i,k)-273.15)
     1246!!      alf = lf0 + clmci*(t(i,k)-273.15)
     1247!!      als = alf + alv
     1248      DO j = 1,4
     1249        DO i = 1, ncum
     1250! ori       if(k.ge.(icb(i)+1))then
     1251          IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1252            tg = tp(i, k)
     1253            IF (tg .gt. Tx) THEN
     1254              es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
     1255              qg = eps*es/(p(i,k)-es*(1.-eps))
     1256            ELSE
     1257              esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
     1258              qg = eps*esi/(p(i,k)-esi*(1.-eps))
     1259            ENDIF
     1260! Ice fraction
     1261            ff = 0.
     1262            fp = 1./(Tx - Tm)
     1263            IF (tg < Tx) THEN
     1264              IF (tg > Tm) THEN
     1265                ff = (Tx - tg)*fp
     1266              ELSE
     1267                ff = 1.
     1268              ENDIF ! (tg > Tm)
     1269            ENDIF ! (tg < Tx)
     1270! Intermediate variables
     1271            aa = cpd + (cl-cpd)*qnk(i) + lv(i,k)*lv(i,k)*qg/(rrv*tg*tg)
     1272            ahg = (cpd + (cl-cpd)*qnk(i))*tg + lv(i,k)*qg - &
     1273                  lf(i,k)*ff*(qnk(i) - qg) + gz(i,k)
     1274            dd = lf(i,k)*lv(i,k)*qg/(rrv*tg*tg)
     1275            ddelta = lf(i,k)*(qnk(i) - qg)
     1276            bb = aa + ddelta*fp + dd*fp*(Tx-tg)
     1277! Compute Zx and Zm
     1278            coefx = aa
     1279            coefm = aa + dd
     1280            IF (tg .gt. Tx) THEN
     1281              Zx = ahg            + coefx*(Tx - tg)
     1282              Zm = ahg - ddelta   + coefm*(Tm - tg)
     1283            ELSE
     1284              IF (tg .gt. Tm) THEN
     1285                Zx = ahg          + (coefx +fp*ddelta)*(Tx - Tg)
     1286                Zm = ahg          + (coefm +fp*ddelta)*(Tm - Tg)
     1287              ELSE
     1288                Zx = ahg + ddelta + coefx*(Tx - tg)
     1289                Zm = ahg          + coefm*(Tm - tg)
     1290              ENDIF ! (tg .gt. Tm)
     1291            ENDIF ! (tg .gt. Tx)
     1292! Compute the masks Um, U, Ux
     1293            Um = (sign(1., Zm-ah0(i))+1.)/2.
     1294            Ux = (sign(1., ah0(i)-Zx)+1.)/2.
     1295            U = (1. - Um)*(1. - Ux)
     1296! Compute the updated parcell temperature Tp : 3 cases depending on tg value
     1297            IF (tg .gt. Tx) THEN
     1298              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg + ddelta*fp*(Tx-tg))
     1299              Tp(i,k) = tg + &
     1300                  Um*  (ah0(i) - ahg + ddelta)           /(aa + dd) + &
     1301                  U *2*(ah0(i) - ahg + ddelta*fp*(Tx-tg))/(bb + sqrt(discr)) + &
     1302                  Ux*  (ah0(i) - ahg)                    /aa
     1303            ELSEIF (tg .gt. Tm) THEN
     1304              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg)
     1305              Tp(i,k) = tg + &
     1306                  Um*  (ah0(i) - ahg + ddelta*fp*(tg-Tm))/(aa + dd) + &
     1307                  U *2*(ah0(i) - ahg)                    /(bb + sqrt(discr)) + &
     1308                  Ux*  (ah0(i) - ahg + ddelta*fp*(tg-Tx))/aa
     1309            ELSE
     1310              discr = bb*bb - 4*dd*fp*(ah0(i) - ahg + ddelta*fp*(Tm-tg))
     1311              Tp(i,k) = tg + &
     1312                  Um*  (ah0(i) - ahg)                    /(aa + dd) + &
     1313                  U *2*(ah0(i) - ahg + ddelta*fp*(Tm-tg))/(bb + sqrt(discr)) + &
     1314                  Ux*  (ah0(i) - ahg - ddelta)           /aa
     1315            ENDIF ! (tg .gt. Tx)
     1316!
     1317!!     print *,' j, k, Um, U, Ux, aa, bb, discr, dd, ddelta ', j, k, Um, U, Ux, aa, bb, discr, dd, ddelta
     1318!!     print *,' j, k, ah0(i), ahg, tg, qg, tp(i,k), ff ', j, k, ah0(i), ahg, tg, qg, tp(i,k), ff
     1319          END IF ! (k>=(icbs(i)+1))
     1320        END DO ! i = 1, ncum
     1321      END DO ! j = 1,4
     1322      DO i = 1, ncum
     1323        IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1324          tg = tp(i, k)
     1325          IF (tg .gt. Tx) THEN
     1326            es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
     1327            qg = eps*es/(p(i,k)-es*(1.-eps))
     1328          ELSE
     1329            esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
     1330            qg = eps*esi/(p(i,k)-esi*(1.-eps))
     1331          ENDIF
     1332          clw(i, k) = qnk(i) - qg
     1333          clw(i, k) = max(0.0, clw(i,k))
     1334          tvp(i, k) = max(0., tp(i,k)*(1.+qg/eps-qnk(i)))
     1335! print*,tvp(i,k),'tvp'
     1336          IF (clw(i,k)<1.E-11) THEN
     1337            tp(i, k) = tv(i, k)
     1338            tvp(i, k) = tv(i, k)
     1339          END IF ! (clw(i,k)<1.E-11)
     1340        END IF ! (k>=(icbs(i)+1))
     1341      END DO ! i = 1, ncum
     1342    END DO ! k = minorig + 1, nl
     1343!----------------------------------------------------------------------------
     1344!
     1345  ELSE IF (icvflag_Tpa == 1) THEN  ! (icvflag_Tpa == 2)
     1346!
     1347!----------------------------------------------------------------------------
     1348!
     1349!  Ice fraction
     1350!
     1351    DO k = minorig + 1, nl
     1352      DO i = 1, ncum
     1353          frac(i, k) = (Tx - t(i,k))/(Tx - Tm)
     1354          frac(i, k) = min(max(frac(i,k),0.0), 1.0)
     1355      END DO
     1356    END DO
     1357! Below cloud base, set ice fraction to cloud base value
     1358    DO k = 1, nl
     1359      DO i = 1, ncum
     1360        IF (k<icb(i)) THEN
     1361          frac(i,k) = frac(i,icb(i))
     1362        END IF
     1363      END DO
     1364    END DO
     1365
     1366    DO k = minorig + 1, nl
     1367      DO i = 1,ncum
     1368        tp(i,k) = t(i,k)
     1369      ENDDO
     1370!!      alv = lv0 - clmcpv*(t(i,k)-273.15)
     1371!!      alf = lf0 + clmci*(t(i,k)-273.15)
     1372!!      als = alf + alv
     1373      DO j = 1,4
     1374        DO i = 1, ncum
     1375! ori       if(k.ge.(icb(i)+1))then
     1376          IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1377            tg = tp(i, k)
     1378            IF (tg .gt. Tx) THEN
     1379              es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
     1380              qg = eps*es/(p(i,k)-es*(1.-eps))
     1381              dqgdT = lv(i,k)*qg/(rrv*tg*tg)
     1382            ELSE
     1383              esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
     1384              qg = eps*esi/(p(i,k)-esi*(1.-eps))
     1385              dqgdT = (lv(i,k)+lf(i,k))*qg/(rrv*tg*tg)
     1386            ENDIF
     1387            IF (qsat_depends_on_qt) THEN
     1388              dqgdT = dqgdT*(1.-qnk(i))/(1.-qg)**2
     1389              qg = qg*(1.-qnk(i))/(1.-qg)           
     1390            ENDIF
     1391            ahg = (cpd + (cl-cpd)*qnk(i))*tg + lv(i,k)*qg - &
     1392                  lf(i,k)*frac(i,k)*(qnk(i) - qg) + gz(i,k)
     1393            Tp(i,k) = tg + (ah0(i) - ahg)/ &
     1394                    (cpd + (cl-cpd)*qnk(i) + (lv(i,k)+frac(i,k)*lf(i,k))*dqgdT)
     1395!!   print *,'undilute2 iterations k, Tp(i,k), ah0(i), ahg ', &
     1396!!                                 k, Tp(i,k), ah0(i), ahg
     1397          END IF ! (k>=(icbs(i)+1))
     1398        END DO ! i = 1, ncum
     1399      END DO ! j = 1,4
     1400      DO i = 1, ncum
     1401        IF (k>=(icbs(i)+1)) THEN                                ! convect3
     1402          tg = tp(i, k)
     1403          IF (tg .gt. Tx) THEN
     1404            es = 6.112*exp(17.67*(tg - 273.15)/(tg + 243.5 - 273.15))
     1405            qg = eps*es/(p(i,k)-es*(1.-eps))
     1406          ELSE
     1407            esi = exp(23.33086-(6111.72784/tg)+0.15215*log(tg))
     1408            qg = eps*esi/(p(i,k)-esi*(1.-eps))
     1409          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
     1414          clw(i, k) = max(0.0, clw(i,k))
     1415          tvp(i, k) = max(0., tp(i,k)*(1.+qg/eps-qnk(i)))
     1416! print*,tvp(i,k),'tvp'
     1417          IF (clw(i,k)<1.E-11) THEN
     1418            tp(i, k) = tv(i, k)
     1419            tvp(i, k) = tv(i, k)
     1420          END IF ! (clw(i,k)<1.E-11)
     1421        END IF ! (k>=(icbs(i)+1))
     1422      END DO ! i = 1, ncum
     1423    END DO ! k = minorig + 1, nl
     1424!
     1425!----------------------------------------------------------------------------
     1426!
     1427  ELSE ! (icvflag_Tpa == 2) ELSE IF(icvflag_Tpa == 1)
     1428!
     1429!----------------------------------------------------------------------------
     1430!
    12041431  DO k = minorig + 1, nl
    12051432    DO i = 1, ncum
    … …  
    13581585  END DO
    13591586
    1360   IF (prt_level >= 10) THEN
    1361     print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
    1362                           (k, tp(1,k), tvp(1,k), k = 1,nl)
    1363   ENDIF
    1364 
     1587!----------------------------------------------------------------------------
     1588!
     1589  ENDIF ! (icvflag_Tpa == 2) ELSEIF (icvflag_Tpa == 1) ELSE
     1590!
     1591!----------------------------------------------------------------------------
     1592!
    13651593! =====================================================================
    1366 ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF
    1367 ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD
     1594! --- SET THE PRECIPITATION EFFICIENCIES
    13681595! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I)
    13691596! =====================================================================
    1370 !
    1371 !jyg<
    1372   DO k = 1, nl
    1373     DO i = 1, ncum
    1374       ep(i, k) = 0.0
    1375       sigp(i, k) = spfac
    1376     END DO
    1377   END DO
    1378 !>jyg
    13791597!
    13801598  IF (flag_epkeorig/=1) THEN
    … …  
    14131631    END DO
    14141632  END IF
     1633!
     1634!   =========================================================================
     1635  IF (prt_level >= 10) THEN
     1636    print *,'cv3_undilute2.1. tp(1,k), tvp(1,k) ', &
     1637                          (k, tp(1,k), tvp(1,k), k = 1,nl)
     1638  ENDIF
    14151639!
    14161640! =====================================================================
    … …  
    16511875      DO i = 1, ncum
    16521876        IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN
     1877!jyg< frac computation moved to beginning of cv3_undilute2.
     1878!     kept here for compatibility test with CMip6 version
    16531879          frac(i, k) = 1. - (t(i,k)-243.15)/(263.15-243.15)
    16541880          frac(i, k) = min(max(frac(i,k),0.0), 1.0)
    … …  
    16591885    END DO
    16601886! Below cloud base, set ice fraction to cloud base value
    1661     DO k = 1, nl
    1662       DO i = 1, ncum
    1663         IF (k<icb(i)) THEN
    1664           frac(i,k) = frac(i,icb(i))
    1665         END IF
    1666       END DO
    1667     END DO
    1668 !
    1669   ELSE
     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
     1894!
     1895  ELSE   ! (cvflag_ice)
    16701896!
    16711897    DO k = minorig + 1, nl

  • LMDZ6/trunk/libf/phylmd/cv3param.h

    r2905 r3492  
    1212      logical ok_convstop
    1313      logical ok_intermittent
     14      logical qsat_depends_on_qt
    1415      integer noff, minorig, nl, nlp, nlm
    1516      integer cv_flag_feed
    … …  
    4546                      ,ok_optim_yield &
    4647                      ,ok_entrain &
    47                       ,ok_homo_tend
     48                      ,ok_homo_tend &
     49                      ,qsat_depends_on_qt
    4850!$OMP THREADPRIVATE(/cv3param/)
    4951

  • LMDZ6/trunk/libf/phylmd/cv_driver.F90

    r3409 r3492  
    568568                ,cape,ep,hp,icb,inb,clw,nk,t,h,lv &
    569569                ,epmax_diag)
    570         ! on écrase ep et recalcule hp
     570        ! on écrase ep et recalcule hp
    571571    END IF
    572572
    … …  
    681681! ==================================================================
    682682SUBROUTINE cv_flag(iflag_ice_thermo)
     683
     684  USE ioipsl_getin_p_mod, ONLY : getin_p
     685
    683686  IMPLICIT NONE
    684687
    … …  
    693696  cvflag_grav = .TRUE.
    694697  cvflag_ice = iflag_ice_thermo >= 1
     698  !
     699! si icvflag_Tpa=0, alors la fraction de glace dans l'ascendance adiabatique est
     700  ! fonction de la temperature de l'environnement et la temperature de l'ascendance est
     701  ! calculee en deux itérations, une en supposant qu'il n'y a pas de glace et l'autre
     702  ! en ajoutant la glace (ancien schéma d'Arnaud Jam).
     703! si icvflag_Tpa=1, alors la fraction de glace dans l'ascendance adiabatique est
     704  ! fonction de la temperature de l'environnement et la temperature de l'ascendance est
     705  ! calculee en une seule iteration.
     706! si icvflag_Tpa=2, alors la fraction de glace dans l'ascendance adiabatique est
     707  ! fonction de la temperature de l'ascendance et la temperature de l'ascendance est
     708  ! calculee en une seule iteration.
     709  icvflag_Tpa=0 
     710  call getin_p('icvflag_Tpa', icvflag_Tpa)
    695711
    696712  RETURN

  • LMDZ6/trunk/libf/phylmd/cvflag.h

    r1992 r3492  
    44      logical cvflag_grav
    55      logical cvflag_ice
     6      integer icvflag_Tpa
    67
    7       COMMON /cvflag/ cvflag_grav, cvflag_ice 
     8      COMMON /cvflag/ icvflag_Tpa, cvflag_grav, cvflag_ice
    89!$OMP THREADPRIVATE(/cvflag/)
Note: See TracChangeset for help on using the changeset viewer.