! $Header$ SUBROUTINE conccm(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, rain, snow, & kbascm, ktopcm) USE dimphy USE lmdz_yoethf USE lmdz_yomcst IMPLICIT NONE ! ====================================================================== ! Auteur(s): Z.X. Li (LMD/CNRS) date: le 14 mars 1996 ! Objet: Schema simple (avec flux de masse) pour la convection ! (schema standard du modele NCAR CCM2) ! ====================================================================== ! Entree: REAL dtime ! pas d'integration REAL paprs(klon, klev + 1) ! pression inter-couche (Pa) REAL pplay(klon, klev) ! pression au milieu de couche (Pa) REAL t(klon, klev) ! temperature (K) REAL q(klon, klev) ! humidite specifique (g/g) REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s) ! Sortie: REAL d_t(klon, klev) ! incrementation temperature REAL d_q(klon, klev) ! incrementation vapeur REAL rain(klon) ! pluie (mm/s) REAL snow(klon) ! neige (mm/s) INTEGER kbascm(klon) ! niveau du bas de convection INTEGER ktopcm(klon) ! niveau du haut de convection REAL pt(klon, klev) REAL pq(klon, klev) REAL pres(klon, klev) REAL dp(klon, klev) REAL zgeom(klon, klev) REAL cmfprs(klon) REAL cmfprt(klon) INTEGER ntop(klon) INTEGER nbas(klon) INTEGER i, k REAL zlvdcp, zlsdcp, zdelta, zz, za, zb LOGICAL usekuo ! utiliser convection profonde (schema Kuo) PARAMETER (usekuo = .TRUE.) REAL d_t_bis(klon, klev) REAL d_q_bis(klon, klev) REAL rain_bis(klon) REAL snow_bis(klon) INTEGER ibas_bis(klon) INTEGER itop_bis(klon) REAL d_ql_bis(klon, klev) REAL rneb_bis(klon, klev) ! initialiser les variables de sortie (pour securite) DO i = 1, klon rain(i) = 0.0 snow(i) = 0.0 kbascm(i) = 0 ktopcm(i) = 0 END DO DO k = 1, klev DO i = 1, klon d_t(i, k) = 0.0 d_q(i, k) = 0.0 END DO END DO ! preparer les variables d'entree (attention: l'ordre des niveaux ! verticaux augmente du haut vers le bas) DO k = 1, klev DO i = 1, klon pt(i, k) = t(i, klev - k + 1) pq(i, k) = q(i, klev - k + 1) pres(i, k) = pplay(i, klev - k + 1) dp(i, k) = paprs(i, klev + 1 - k) - paprs(i, klev + 1 - k + 1) END DO END DO DO i = 1, klon zgeom(i, klev) = rd * t(i, 1) / (0.5 * (paprs(i, 1) + pplay(i, & 1))) * (paprs(i, 1) - pplay(i, 1)) END DO DO k = 2, klev DO i = 1, klon zgeom(i, klev + 1 - k) = zgeom(i, klev + 1 - k + 1) + rd * 0.5 * (t(i, k - 1) + t(i, k)) / & paprs(i, k) * (pplay(i, k - 1) - pplay(i, k)) END DO END DO CALL cmfmca(dtime, pres, dp, zgeom, pt, pq, cmfprt, cmfprs, ntop, nbas) DO k = 1, klev DO i = 1, klon d_q(i, klev + 1 - k) = pq(i, k) - q(i, klev + 1 - k) d_t(i, klev + 1 - k) = pt(i, k) - t(i, klev + 1 - k) END DO END DO DO i = 1, klon rain(i) = cmfprt(i) * rhoh2o snow(i) = cmfprs(i) * rhoh2o kbascm(i) = klev + 1 - nbas(i) ktopcm(i) = klev + 1 - ntop(i) END DO IF (usekuo) THEN CALL conkuo(dtime, paprs, pplay, t, q, conv_q, d_t_bis, d_q_bis, & d_ql_bis, rneb_bis, rain_bis, snow_bis, ibas_bis, itop_bis) DO k = 1, klev DO i = 1, klon d_t(i, k) = d_t(i, k) + d_t_bis(i, k) d_q(i, k) = d_q(i, k) + d_q_bis(i, k) END DO END DO DO i = 1, klon rain(i) = rain(i) + rain_bis(i) snow(i) = snow(i) + snow_bis(i) kbascm(i) = min(kbascm(i), ibas_bis(i)) ktopcm(i) = max(ktopcm(i), itop_bis(i)) END DO DO k = 1, klev ! eau liquide convective est DO i = 1, klon ! dispersee dans l'air zlvdcp = rlvtt / rcpd / (1.0 + rvtmp2 * q(i, k)) zlsdcp = rlstt / rcpd / (1.0 + rvtmp2 * q(i, k)) zdelta = max(0., sign(1., rtt - t(i, k))) zz = d_ql_bis(i, k) ! re-evap. de l'eau liquide zb = max(0.0, zz) za = -max(0.0, zz) * (zlvdcp * (1. - zdelta) + zlsdcp * zdelta) d_t(i, k) = d_t(i, k) + za d_q(i, k) = d_q(i, k) + zb END DO END DO END IF END SUBROUTINE conccm SUBROUTINE cmfmca(deltat, p, dp, gz, tb, shb, cmfprt, cmfprs, cnt, cnb) USE dimphy USE lmdz_yoethf USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep USE lmdz_yomcst IMPLICIT NONE ! ----------------------------------------------------------------------- ! Moist convective mass flux procedure: ! If stratification is unstable to nonentraining parcel ascent, ! complete an adjustment making use of a simple cloud model ! Code generalized to allow specification of parcel ("updraft") ! properties, as well as convective transport of an arbitrary ! number of passive constituents (see cmrb array). ! ----------------------------Code History------------------------------- ! Original version: J. J. Hack, March 22, 1990 ! Standardized: J. Rosinski, June 1992 ! Reviewed: J. Hack, G. Taylor, August 1992 ! Adaptation au LMD: Z.X. Li, mars 1996 (reference: Hack 1994, ! J. Geophys. Res. vol 99, D3, 5551-5568). J'ai ! introduit les constantes et les fonctions thermo- ! dynamiques du Centre Europeen. J'ai elimine le ! re-indicage du code en esperant que cela pourra ! simplifier la lecture et la comprehension. ! ----------------------------------------------------------------------- INTEGER pcnst ! nombre de traceurs passifs PARAMETER (pcnst = 1) ! ------------------------------Arguments-------------------------------- ! Input arguments REAL deltat ! time step (seconds) REAL p(klon, klev) ! pressure REAL dp(klon, klev) ! delta-p REAL gz(klon, klev) ! geopotential (a partir du sol) REAL thtap(klon) ! PBL perturbation theta REAL shp(klon) ! PBL perturbation specific humidity REAL pblh(klon) ! PBL height (provided by PBL routine) REAL cmrp(klon, pcnst) ! constituent perturbations in PBL ! Updated arguments: REAL tb(klon, klev) ! temperature (t bar) REAL shb(klon, klev) ! specific humidity (sh bar) REAL cmrb(klon, klev, pcnst) ! constituent mixing ratios (cmr bar) ! Output arguments REAL cmfdt(klon, klev) ! dT/dt due to moist convection REAL cmfdq(klon, klev) ! dq/dt due to moist convection REAL cmfmc(klon, klev) ! moist convection cloud mass flux REAL cmfdqr(klon, klev) ! dq/dt due to convective rainout REAL cmfsl(klon, klev) ! convective lw static energy flux REAL cmflq(klon, klev) ! convective total water flux REAL cmfprt(klon) ! convective precipitation rate REAL cmfprs(klon) ! convective snowfall rate REAL qc(klon, klev) ! dq/dt due to rainout terms INTEGER cnt(klon) ! top level of convective activity INTEGER cnb(klon) ! bottom level of convective activity ! ------------------------------Parameters------------------------------- REAL c0 ! rain water autoconversion coefficient PARAMETER (c0 = 1.0E-4) REAL dzmin ! minimum convective depth for precipitation PARAMETER (dzmin = 0.0) REAL betamn ! minimum overshoot parameter PARAMETER (betamn = 0.10) REAL cmftau ! characteristic adjustment time scale PARAMETER (cmftau = 3600.) INTEGER limcnv ! top interface level limit for convection PARAMETER (limcnv = 1) REAL tpmax ! maximum acceptable t perturbation (degrees C) PARAMETER (tpmax = 1.50) REAL shpmax ! maximum acceptable q perturbation (g/g) PARAMETER (shpmax = 1.50E-3) REAL tiny ! arbitrary small num used in transport estimates PARAMETER (tiny = 1.0E-36) REAL eps ! convergence criteria (machine dependent) PARAMETER (eps = 1.0E-13) REAL tmelt ! freezing point of water(req'd for rain vs snow) PARAMETER (tmelt = 273.15) REAL ssfac ! supersaturation bound (detrained air) PARAMETER (ssfac = 1.001) ! ---------------------------Local workspace----------------------------- REAL gam(klon, klev) ! L/cp (d(qsat)/dT) REAL sb(klon, klev) ! dry static energy (s bar) REAL hb(klon, klev) ! moist static energy (h bar) REAL shbs(klon, klev) ! sat. specific humidity (sh bar star) REAL hbs(klon, klev) ! sat. moist static energy (h bar star) REAL shbh(klon, klev + 1) ! specific humidity on interfaces REAL sbh(klon, klev + 1) ! s bar on interfaces REAL hbh(klon, klev + 1) ! h bar on interfaces REAL cmrh(klon, klev + 1) ! interface constituent mixing ratio REAL prec(klon) ! instantaneous total precipitation REAL dzcld(klon) ! depth of convective layer (m) REAL beta(klon) ! overshoot parameter (fraction) REAL betamx ! local maximum on overshoot REAL eta(klon) ! convective mass flux (kg/m^2 s) REAL etagdt ! eta*grav*deltat REAL cldwtr(klon) ! cloud water (mass) REAL rnwtr(klon) ! rain water (mass) REAL sc(klon) ! dry static energy ("in-cloud") REAL shc(klon) ! specific humidity ("in-cloud") REAL hc(klon) ! moist static energy ("in-cloud") REAL cmrc(klon) ! constituent mix rat ("in-cloud") REAL dq1(klon) ! shb convective change (lower lvl) REAL dq2(klon) ! shb convective change (mid level) REAL dq3(klon) ! shb convective change (upper lvl) REAL ds1(klon) ! sb convective change (lower lvl) REAL ds2(klon) ! sb convective change (mid level) REAL ds3(klon) ! sb convective change (upper lvl) REAL dcmr1(klon) ! cmrb convective change (lower lvl) REAL dcmr2(klon) ! cmrb convective change (mid level) REAL dcmr3(klon) ! cmrb convective change (upper lvl) REAL flotab(klon) ! hc - hbs (mesure d'instabilite) LOGICAL ldcum(klon) ! .TRUE. si la convection existe LOGICAL etagt0 ! true if eta > 0.0 REAL dt ! current 2 delta-t (model time step) REAL cats ! modified characteristic adj. time REAL rdt ! 1./dt REAL qprime ! modified specific humidity pert. REAL tprime ! modified thermal perturbation REAL pblhgt ! bounded pbl height (max[pblh,1m]) REAL fac1 ! intermediate scratch variable REAL shprme ! intermediate specific humidity pert. REAL qsattp ! saturation mixing ratio for ! thermally perturbed PBL parcels REAL dz ! local layer depth REAL b1 ! bouyancy measure in detrainment lvl REAL b2 ! bouyancy measure in condensation lvl REAL g ! bounded vertical gradient of hb REAL tmass ! total mass available for convective exchange REAL denom ! intermediate scratch variable REAL qtest1 ! used in negative q test (middle lvl) REAL qtest2 ! used in negative q test (lower lvl) REAL fslkp ! flux lw static energy (bot interface) REAL fslkm ! flux lw static energy (top interface) REAL fqlkp ! flux total water (bottom interface) REAL fqlkm ! flux total water (top interface) REAL botflx ! bottom constituent mixing ratio flux REAL topflx ! top constituent mixing ratio flux REAL efac1 ! ratio cmrb to convectively induced change (bot lvl) REAL efac2 ! ratio cmrb to convectively induced change (mid lvl) REAL efac3 ! ratio cmrb to convectively induced change (top lvl) INTEGER i, k ! indices horizontal et vertical INTEGER km1 ! k-1 (index offset) INTEGER kp1 ! k+1 (index offset) INTEGER m ! constituent index INTEGER ktp ! temporary index used to track top INTEGER is ! nombre de points a ajuster REAL tmp1, tmp2, tmp3, tmp4 REAL zx_t, zx_p, zx_q, zx_qs, zx_gam REAL zcor, zdelta, zcvm5 REAL qhalf, sh1, sh2, shbs1, shbs2 qhalf(sh1, sh2, shbs1, shbs2) = min(max(sh1, sh2), & (shbs2 * sh1 + shbs1 * sh2) / (shbs1 + shbs2)) ! ----------------------------------------------------------------------- ! pas de traceur pour l'instant DO m = 1, pcnst DO k = 1, klev DO i = 1, klon cmrb(i, k, m) = 0.0 END DO END DO END DO ! Les perturbations de la couche limite sont zero pour l'instant DO m = 1, pcnst DO i = 1, klon cmrp(i, m) = 0.0 END DO END DO DO i = 1, klon thtap(i) = 0.0 shp(i) = 0.0 pblh(i) = 1.0 END DO ! Ensure that characteristic adjustment time scale (cmftau) assumed ! in estimate of eta isn't smaller than model time scale (deltat) dt = deltat cats = max(dt, cmftau) rdt = 1.0 / dt ! Compute sb,hb,shbs,hbs DO k = 1, klev DO i = 1, klon zx_t = tb(i, k) zx_p = p(i, k) zx_q = shb(i, k) zdelta = max(0., sign(1., rtt - zx_t)) zcvm5 = r5les * rlvtt * (1. - zdelta) + r5ies * rlstt * zdelta zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * zx_q) zx_qs = r2es * foeew(zx_t, zdelta) / zx_p zx_qs = min(0.5, zx_qs) zcor = 1. / (1. - retv * zx_qs) zx_qs = zx_qs * zcor zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) shbs(i, k) = zx_qs gam(i, k) = zx_gam END DO END DO DO k = limcnv, klev DO i = 1, klon sb(i, k) = rcpd * tb(i, k) + gz(i, k) hb(i, k) = sb(i, k) + rlvtt * shb(i, k) hbs(i, k) = sb(i, k) + rlvtt * shbs(i, k) END DO END DO ! Compute sbh, shbh DO k = limcnv + 1, klev km1 = k - 1 DO i = 1, klon sbh(i, k) = 0.5 * (sb(i, km1) + sb(i, k)) shbh(i, k) = qhalf(shb(i, km1), shb(i, k), shbs(i, km1), shbs(i, k)) hbh(i, k) = sbh(i, k) + rlvtt * shbh(i, k) END DO END DO ! Specify properties at top of model (not used, but filling anyway) DO i = 1, klon sbh(i, limcnv) = sb(i, limcnv) shbh(i, limcnv) = shb(i, limcnv) hbh(i, limcnv) = hb(i, limcnv) END DO ! Zero vertically independent control, tendency & diagnostic arrays DO i = 1, klon prec(i) = 0.0 dzcld(i) = 0.0 cnb(i) = 0 cnt(i) = klev + 1 END DO DO k = 1, klev DO i = 1, klon cmfdt(i, k) = 0. cmfdq(i, k) = 0. cmfdqr(i, k) = 0. cmfmc(i, k) = 0. cmfsl(i, k) = 0. cmflq(i, k) = 0. END DO END DO ! Begin moist convective mass flux adjustment procedure. ! Formalism ensures that negative cloud liquid water can never occur DO k = klev - 1, limcnv + 1, -1 km1 = k - 1 kp1 = k + 1 DO i = 1, klon eta(i) = 0.0 beta(i) = 0.0 ds1(i) = 0.0 ds2(i) = 0.0 ds3(i) = 0.0 dq1(i) = 0.0 dq2(i) = 0.0 dq3(i) = 0.0 ! Specification of "cloud base" conditions qprime = 0.0 tprime = 0.0 ! Assign tprime within the PBL to be proportional to the quantity ! thtap (which will be bounded by tpmax), passed to this routine by ! the PBL routine. Don't allow perturbation to produce a dry ! adiabatically unstable parcel. Assign qprime within the PBL to be ! an appropriately modified value of the quantity shp (which will be ! bounded by shpmax) passed to this routine by the PBL routine. The ! quantity qprime should be less than the local saturation value ! (qsattp=qsat[t+tprime,p]). In both cases, thtap and shp are ! linearly reduced toward zero as the PBL top is approached. pblhgt = max(pblh(i), 1.0) IF (gz(i, kp1) / rg<=pblhgt .AND. dzcld(i)==0.0) THEN fac1 = max(0.0, 1.0 - gz(i, kp1) / rg / pblhgt) tprime = min(thtap(i), tpmax) * fac1 qsattp = shbs(i, kp1) + rcpd / rlvtt * gam(i, kp1) * tprime shprme = min(min(shp(i), shpmax) * fac1, max(qsattp - shb(i, kp1), 0.0)) qprime = max(qprime, shprme) ELSE tprime = 0.0 qprime = 0.0 END IF ! Specify "updraft" (in-cloud) thermodynamic properties sc(i) = sb(i, kp1) + rcpd * tprime shc(i) = shb(i, kp1) + qprime hc(i) = sc(i) + rlvtt * shc(i) flotab(i) = hc(i) - hbs(i, k) dz = dp(i, k) * rd * tb(i, k) / rg / p(i, k) IF (flotab(i)>0.0) THEN dzcld(i) = dzcld(i) + dz ELSE dzcld(i) = 0.0 END IF END DO ! Check on moist convective instability is = 0 DO i = 1, klon IF (flotab(i)>0.0) THEN ldcum(i) = .TRUE. is = is + 1 ELSE ldcum(i) = .FALSE. END IF END DO IF (is==0) THEN DO i = 1, klon dzcld(i) = 0.0 END DO GO TO 70 END IF ! Current level just below top level => no overshoot IF (k<=limcnv + 1) THEN DO i = 1, klon IF (ldcum(i)) THEN cldwtr(i) = sb(i, k) - sc(i) + flotab(i) / (1.0 + gam(i, k)) cldwtr(i) = max(0.0, cldwtr(i)) beta(i) = 0.0 END IF END DO GO TO 20 END IF ! First guess at overshoot parameter using crude buoyancy closure ! 10% overshoot assumed as a minimum and 1-c0*dz maximum to start ! If pre-existing supersaturation in detrainment layer, beta=0 ! cldwtr is temporarily equal to RLVTT*l (l=> liquid water) DO i = 1, klon IF (ldcum(i)) THEN cldwtr(i) = sb(i, k) - sc(i) + flotab(i) / (1.0 + gam(i, k)) cldwtr(i) = max(0.0, cldwtr(i)) betamx = 1.0 - c0 * max(0.0, (dzcld(i) - dzmin)) b1 = (hc(i) - hbs(i, km1)) * dp(i, km1) b2 = (hc(i) - hbs(i, k)) * dp(i, k) beta(i) = max(betamn, min(betamx, 1.0 + b1 / b2)) IF (hbs(i, km1)<=hb(i, km1)) beta(i) = 0.0 END IF END DO ! Bound maximum beta to ensure physically realistic solutions ! First check constrains beta so that eta remains positive ! (assuming that eta is already positive for beta equal zero) ! La premiere contrainte de beta est que le flux eta doit etre positif. DO i = 1, klon IF (ldcum(i)) THEN tmp1 = (1.0 + gam(i, k)) * (sc(i) - sbh(i, kp1) + cldwtr(i)) - & (hbh(i, kp1) - hc(i)) * dp(i, k) / dp(i, kp1) tmp2 = (1.0 + gam(i, k)) * (sc(i) - sbh(i, k)) IF ((beta(i) * tmp2 - tmp1)>0.0) THEN betamx = 0.99 * (tmp1 / tmp2) beta(i) = max(0.0, min(betamx, beta(i))) END IF ! Second check involves supersaturation of "detrainment layer" ! small amount of supersaturation acceptable (by ssfac factor) ! La 2e contrainte est que la convection ne doit pas sursaturer ! la "detrainment layer", Neanmoins, une petite sursaturation ! est acceptee (facteur ssfac). IF (hb(i, km1)0.0) THEN betamx = ssfac * (tmp1 / tmp4) beta(i) = max(0.0, min(betamx, beta(i))) END IF ELSE beta(i) = 0.0 END IF ! Third check to avoid introducing 2 delta x thermodynamic ! noise in the vertical ... constrain adjusted h (or theta e) ! so that the adjustment doesn't contribute to "kinks" in h g = min(0.0, hb(i, k) - hb(i, km1)) tmp3 = (hb(i, k) - hb(i, km1) - g) * (cats / dt) / (hc(i) - hbs(i, k)) tmp1 = (1.0 + gam(i, k)) * (sc(i) - sbh(i, kp1) + cldwtr(i)) - & (hbh(i, kp1) - hc(i)) * dp(i, k) / dp(i, kp1) tmp1 = tmp1 / dp(i, k) tmp1 = tmp3 * tmp1 + (hc(i) - hbh(i, kp1)) / dp(i, k) tmp2 = tmp3 * (1.0 + gam(i, k)) * (sc(i) - sbh(i, k)) / dp(i, k) + & (hc(i) - hbh(i, k) - cldwtr(i)) * (1.0 / dp(i, k) + 1.0 / dp(i, kp1)) IF ((beta(i) * tmp2 - tmp1)>0.0) THEN betamx = 0.0 IF (tmp2/=0.0) betamx = tmp1 / tmp2 beta(i) = max(0.0, min(betamx, beta(i))) END IF END IF END DO ! Calculate mass flux required for stabilization. ! Ensure that the convective mass flux, eta, is positive by ! setting negative values of eta to zero.. ! Ensure that estimated mass flux cannot move more than the ! minimum of total mass contained in either layer k or layer k+1. ! Also test for other pathological cases that result in non- ! physical states and adjust eta accordingly. 20 CONTINUE DO i = 1, klon IF (ldcum(i)) THEN beta(i) = max(0.0, beta(i)) tmp1 = hc(i) - hbs(i, k) tmp2 = ((1.0 + gam(i, k)) * (sc(i) - sbh(i, kp1) + cldwtr(i)) - beta(i) * (1.0 + gam(& i, k)) * (sc(i) - sbh(i, k))) / dp(i, k) - (hbh(i, kp1) - hc(i)) / dp(i, kp1) eta(i) = tmp1 / (tmp2 * rg * cats) tmass = min(dp(i, k), dp(i, kp1)) / rg IF (eta(i)>tmass * rdt .OR. eta(i)<=0.0) eta(i) = 0.0 ! Check on negative q in top layer (bound beta) IF (shc(i) - shbh(i, k)<0.0 .AND. beta(i) * eta(i)/=0.0) THEN denom = eta(i) * rg * dt * (shc(i) - shbh(i, k)) / dp(i, km1) beta(i) = max(0.0, min(-0.999 * shb(i, km1) / denom, beta(i))) END IF ! Check on negative q in middle layer (zero eta) qtest1 = shb(i, k) + eta(i) * rg * dt * ((shc(i) - shbh(i, & kp1)) - (1.0 - beta(i)) * cldwtr(i) / rlvtt - beta(i) * (shc(i) - shbh(i, & k))) / dp(i, k) IF (qtest1<=0.0) eta(i) = 0.0 ! Check on negative q in lower layer (bound eta) fac1 = -(shbh(i, kp1) - shc(i)) / dp(i, kp1) qtest2 = shb(i, kp1) - eta(i) * rg * dt * fac1 IF (qtest2<0.0) THEN eta(i) = 0.99 * shb(i, kp1) / (rg * dt * fac1) END IF END IF END DO ! Calculate cloud water, rain water, and thermodynamic changes DO i = 1, klon IF (ldcum(i)) THEN etagdt = eta(i) * rg * dt cldwtr(i) = etagdt * cldwtr(i) / rlvtt / rg rnwtr(i) = (1.0 - beta(i)) * cldwtr(i) ds1(i) = etagdt * (sbh(i, kp1) - sc(i)) / dp(i, kp1) dq1(i) = etagdt * (shbh(i, kp1) - shc(i)) / dp(i, kp1) ds2(i) = (etagdt * (sc(i) - sbh(i, kp1)) + rlvtt * rg * cldwtr(i) - beta(i) * etagdt & * (sc(i) - sbh(i, k))) / dp(i, k) dq2(i) = (etagdt * (shc(i) - shbh(i, kp1)) - rg * rnwtr(i) - beta(i) * etagdt * (shc & (i) - shbh(i, k))) / dp(i, k) ds3(i) = beta(i) * (etagdt * (sc(i) - sbh(i, k)) - rlvtt * rg * cldwtr(i)) / dp(i, & km1) dq3(i) = beta(i) * etagdt * (shc(i) - shbh(i, k)) / dp(i, km1) ! Isolate convective fluxes for later diagnostics fslkp = eta(i) * (sc(i) - sbh(i, kp1)) fslkm = beta(i) * (eta(i) * (sc(i) - sbh(i, k)) - rlvtt * cldwtr(i) * rdt) fqlkp = eta(i) * (shc(i) - shbh(i, kp1)) fqlkm = beta(i) * eta(i) * (shc(i) - shbh(i, k)) ! Update thermodynamic profile (update sb, hb, & hbs later) tb(i, kp1) = tb(i, kp1) + ds1(i) / rcpd tb(i, k) = tb(i, k) + ds2(i) / rcpd tb(i, km1) = tb(i, km1) + ds3(i) / rcpd shb(i, kp1) = shb(i, kp1) + dq1(i) shb(i, k) = shb(i, k) + dq2(i) shb(i, km1) = shb(i, km1) + dq3(i) prec(i) = prec(i) + rnwtr(i) / rhoh2o ! Update diagnostic information for final budget ! Tracking temperature & specific humidity tendencies, ! rainout term, convective mass flux, convective liquid ! water static energy flux, and convective total water flux cmfdt(i, kp1) = cmfdt(i, kp1) + ds1(i) / rcpd * rdt cmfdt(i, k) = cmfdt(i, k) + ds2(i) / rcpd * rdt cmfdt(i, km1) = cmfdt(i, km1) + ds3(i) / rcpd * rdt cmfdq(i, kp1) = cmfdq(i, kp1) + dq1(i) * rdt cmfdq(i, k) = cmfdq(i, k) + dq2(i) * rdt cmfdq(i, km1) = cmfdq(i, km1) + dq3(i) * rdt cmfdqr(i, k) = cmfdqr(i, k) + (rg * rnwtr(i) / dp(i, k)) * rdt cmfmc(i, kp1) = cmfmc(i, kp1) + eta(i) cmfmc(i, k) = cmfmc(i, k) + beta(i) * eta(i) cmfsl(i, kp1) = cmfsl(i, kp1) + fslkp cmfsl(i, k) = cmfsl(i, k) + fslkm cmflq(i, kp1) = cmflq(i, kp1) + rlvtt * fqlkp cmflq(i, k) = cmflq(i, k) + rlvtt * fqlkm qc(i, k) = (rg * rnwtr(i) / dp(i, k)) * rdt END IF END DO ! Next, convectively modify passive constituents DO m = 1, pcnst DO i = 1, klon IF (ldcum(i)) THEN ! If any of the reported values of the constituent is negative in ! the three adjacent levels, nothing will be done to the profile IF ((cmrb(i, kp1, m)<0.0) .OR. (cmrb(i, k, m)<0.0) .OR. (cmrb(i, km1, & m)<0.0)) GO TO 40 ! Specify constituent interface values (linear interpolation) cmrh(i, k) = 0.5 * (cmrb(i, km1, m) + cmrb(i, k, m)) cmrh(i, kp1) = 0.5 * (cmrb(i, k, m) + cmrb(i, kp1, m)) ! Specify perturbation properties of constituents in PBL pblhgt = max(pblh(i), 1.0) IF (gz(i, kp1) / rg<=pblhgt .AND. dzcld(i)==0.) THEN fac1 = max(0.0, 1.0 - gz(i, kp1) / rg / pblhgt) cmrc(i) = cmrb(i, kp1, m) + cmrp(i, m) * fac1 ELSE cmrc(i) = cmrb(i, kp1, m) END IF ! Determine fluxes, flux divergence => changes due to convection ! Logic must be included to avoid producing negative values. A bit ! messy since there are no a priori assumptions about profiles. ! Tendency is modified (reduced) when pending disaster detected. etagdt = eta(i) * rg * dt botflx = etagdt * (cmrc(i) - cmrh(i, kp1)) topflx = beta(i) * etagdt * (cmrc(i) - cmrh(i, k)) dcmr1(i) = -botflx / dp(i, kp1) efac1 = 1.0 efac2 = 1.0 efac3 = 1.0 IF (cmrb(i, kp1, m) + dcmr1(i)<0.0) THEN efac1 = max(tiny, abs(cmrb(i, kp1, m) / dcmr1(i)) - eps) END IF IF (efac1==tiny .OR. efac1>1.0) efac1 = 0.0 dcmr1(i) = -efac1 * botflx / dp(i, kp1) dcmr2(i) = (efac1 * botflx - topflx) / dp(i, k) IF (cmrb(i, k, m) + dcmr2(i)<0.0) THEN efac2 = max(tiny, abs(cmrb(i, k, m) / dcmr2(i)) - eps) END IF IF (efac2==tiny .OR. efac2>1.0) efac2 = 0.0 dcmr2(i) = (efac1 * botflx - efac2 * topflx) / dp(i, k) dcmr3(i) = efac2 * topflx / dp(i, km1) IF (cmrb(i, km1, m) + dcmr3(i)<0.0) THEN efac3 = max(tiny, abs(cmrb(i, km1, m) / dcmr3(i)) - eps) END IF IF (efac3==tiny .OR. efac3>1.0) efac3 = 0.0 efac3 = min(efac2, efac3) dcmr2(i) = (efac1 * botflx - efac3 * topflx) / dp(i, k) dcmr3(i) = efac3 * topflx / dp(i, km1) cmrb(i, kp1, m) = cmrb(i, kp1, m) + dcmr1(i) cmrb(i, k, m) = cmrb(i, k, m) + dcmr2(i) cmrb(i, km1, m) = cmrb(i, km1, m) + dcmr3(i) END IF 40 END DO END DO ! end of m=1,pcnst loop IF (k==limcnv + 1) GO TO 60 ! on ne pourra plus glisser ! Dans la procedure de glissage ascendant, les variables thermo- ! dynamiques des couches k et km1 servent au calcul des couches ! superieures. Elles ont donc besoin d'une mise-a-jour. DO i = 1, klon IF (ldcum(i)) THEN zx_t = tb(i, k) zx_p = p(i, k) zx_q = shb(i, k) zdelta = max(0., sign(1., rtt - zx_t)) zcvm5 = r5les * rlvtt * (1. - zdelta) + r5ies * rlstt * zdelta zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * zx_q) zx_qs = r2es * foeew(zx_t, zdelta) / zx_p zx_qs = min(0.5, zx_qs) zcor = 1. / (1. - retv * zx_qs) zx_qs = zx_qs * zcor zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) shbs(i, k) = zx_qs gam(i, k) = zx_gam zx_t = tb(i, km1) zx_p = p(i, km1) zx_q = shb(i, km1) zdelta = max(0., sign(1., rtt - zx_t)) zcvm5 = r5les * rlvtt * (1. - zdelta) + r5ies * rlstt * zdelta zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * zx_q) zx_qs = r2es * foeew(zx_t, zdelta) / zx_p zx_qs = min(0.5, zx_qs) zcor = 1. / (1. - retv * zx_qs) zx_qs = zx_qs * zcor zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) shbs(i, km1) = zx_qs gam(i, km1) = zx_gam sb(i, k) = sb(i, k) + ds2(i) sb(i, km1) = sb(i, km1) + ds3(i) hb(i, k) = sb(i, k) + rlvtt * shb(i, k) hb(i, km1) = sb(i, km1) + rlvtt * shb(i, km1) hbs(i, k) = sb(i, k) + rlvtt * shbs(i, k) hbs(i, km1) = sb(i, km1) + rlvtt * shbs(i, km1) sbh(i, k) = 0.5 * (sb(i, k) + sb(i, km1)) shbh(i, k) = qhalf(shb(i, km1), shb(i, k), shbs(i, km1), shbs(i, k)) hbh(i, k) = sbh(i, k) + rlvtt * shbh(i, k) sbh(i, km1) = 0.5 * (sb(i, km1) + sb(i, k - 2)) shbh(i, km1) = qhalf(shb(i, k - 2), shb(i, km1), shbs(i, k - 2), & shbs(i, km1)) hbh(i, km1) = sbh(i, km1) + rlvtt * shbh(i, km1) END IF END DO ! Ensure that dzcld is reset if convective mass flux zero ! specify the current vertical extent of the convective activity ! top of convective layer determined by size of overshoot param. 60 CONTINUE DO i = 1, klon etagt0 = eta(i) > 0.0 IF (.NOT. etagt0) dzcld(i) = 0.0 IF (etagt0 .AND. beta(i)>betamn) THEN ktp = km1 ELSE ktp = k END IF IF (etagt0) THEN cnt(i) = min(cnt(i), ktp) cnb(i) = max(cnb(i), k) END IF END DO 70 END DO ! end of k loop ! determine whether precipitation, prec, is frozen (snow) or not DO i = 1, klon IF (tb(i, klev)