Changeset 1992 for LMDZ5/trunk/libf/phylmd/conccm.F90
- Timestamp:
- Mar 5, 2014, 2:19:12 PM (10 years ago)
- File:
-
- 1 moved
-
LMDZ5/trunk/libf/phylmd/conccm.F90 (moved) (moved from LMDZ5/trunk/libf/phylmd/conccm.F) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ5/trunk/libf/phylmd/conccm.F90
r1988 r1992 1 ! 1 2 2 ! $Header$ 3 ! 4 SUBROUTINE conccm (dtime,paprs,pplay,t,q,conv_q, 5 s d_t, d_q, rain, snow, kbascm, ktopcm) 6 c 7 USE dimphy 8 IMPLICIT none 9 c====================================================================== 10 c Auteur(s): Z.X. Li (LMD/CNRS) date: le 14 mars 1996 11 c Objet: Schema simple (avec flux de masse) pour la convection 12 c (schema standard du modele NCAR CCM2) 13 c====================================================================== 14 cym#include "dimensions.h" 15 cym#include "dimphy.h" 16 #include "YOMCST.h" 17 #include "YOETHF.h" 18 c 19 c Entree: 20 REAL dtime ! pas d'integration 21 REAL paprs(klon,klev+1) ! pression inter-couche (Pa) 22 REAL pplay(klon,klev) ! pression au milieu de couche (Pa) 23 REAL t(klon,klev) ! temperature (K) 24 REAL q(klon,klev) ! humidite specifique (g/g) 25 REAL conv_q(klon,klev) ! taux de convergence humidite (g/g/s) 26 c Sortie: 27 REAL d_t(klon,klev) ! incrementation temperature 28 REAL d_q(klon,klev) ! incrementation vapeur 29 REAL rain(klon) ! pluie (mm/s) 30 REAL snow(klon) ! neige (mm/s) 31 INTEGER kbascm(klon) ! niveau du bas de convection 32 INTEGER ktopcm(klon) ! niveau du haut de convection 33 c 34 REAL pt(klon,klev) 35 REAL pq(klon,klev) 36 REAL pres(klon,klev) 37 REAL dp(klon,klev) 38 REAL zgeom(klon,klev) 39 REAL cmfprs(klon) 40 REAL cmfprt(klon) 41 INTEGER ntop(klon) 42 INTEGER nbas(klon) 43 INTEGER i, k 44 REAL zlvdcp, zlsdcp, zdelta, zz, za, zb 45 c 46 LOGICAL usekuo ! utiliser convection profonde (schema Kuo) 47 PARAMETER (usekuo=.TRUE.) 48 c 49 REAL d_t_bis(klon,klev) 50 REAL d_q_bis(klon,klev) 51 REAL rain_bis(klon) 52 REAL snow_bis(klon) 53 INTEGER ibas_bis(klon) 54 INTEGER itop_bis(klon) 55 REAL d_ql_bis(klon,klev) 56 REAL rneb_bis(klon,klev) 57 c 58 c initialiser les variables de sortie (pour securite) 3 4 SUBROUTINE conccm(dtime, paprs, pplay, t, q, conv_q, d_t, d_q, rain, snow, & 5 kbascm, ktopcm) 6 7 USE dimphy 8 IMPLICIT NONE 9 ! ====================================================================== 10 ! Auteur(s): Z.X. Li (LMD/CNRS) date: le 14 mars 1996 11 ! Objet: Schema simple (avec flux de masse) pour la convection 12 ! (schema standard du modele NCAR CCM2) 13 ! ====================================================================== 14 ! ym#include "dimensions.h" 15 ! ym#include "dimphy.h" 16 include "YOMCST.h" 17 include "YOETHF.h" 18 19 ! Entree: 20 REAL dtime ! pas d'integration 21 REAL paprs(klon, klev+1) ! pression inter-couche (Pa) 22 REAL pplay(klon, klev) ! pression au milieu de couche (Pa) 23 REAL t(klon, klev) ! temperature (K) 24 REAL q(klon, klev) ! humidite specifique (g/g) 25 REAL conv_q(klon, klev) ! taux de convergence humidite (g/g/s) 26 ! Sortie: 27 REAL d_t(klon, klev) ! incrementation temperature 28 REAL d_q(klon, klev) ! incrementation vapeur 29 REAL rain(klon) ! pluie (mm/s) 30 REAL snow(klon) ! neige (mm/s) 31 INTEGER kbascm(klon) ! niveau du bas de convection 32 INTEGER ktopcm(klon) ! niveau du haut de convection 33 34 REAL pt(klon, klev) 35 REAL pq(klon, klev) 36 REAL pres(klon, klev) 37 REAL dp(klon, klev) 38 REAL zgeom(klon, klev) 39 REAL cmfprs(klon) 40 REAL cmfprt(klon) 41 INTEGER ntop(klon) 42 INTEGER nbas(klon) 43 INTEGER i, k 44 REAL zlvdcp, zlsdcp, zdelta, zz, za, zb 45 46 LOGICAL usekuo ! utiliser convection profonde (schema Kuo) 47 PARAMETER (usekuo=.TRUE.) 48 49 REAL d_t_bis(klon, klev) 50 REAL d_q_bis(klon, klev) 51 REAL rain_bis(klon) 52 REAL snow_bis(klon) 53 INTEGER ibas_bis(klon) 54 INTEGER itop_bis(klon) 55 REAL d_ql_bis(klon, klev) 56 REAL rneb_bis(klon, klev) 57 58 ! initialiser les variables de sortie (pour securite) 59 DO i = 1, klon 60 rain(i) = 0.0 61 snow(i) = 0.0 62 kbascm(i) = 0 63 ktopcm(i) = 0 64 END DO 65 DO k = 1, klev 66 DO i = 1, klon 67 d_t(i, k) = 0.0 68 d_q(i, k) = 0.0 69 END DO 70 END DO 71 72 ! preparer les variables d'entree (attention: l'ordre des niveaux 73 ! verticaux augmente du haut vers le bas) 74 DO k = 1, klev 75 DO i = 1, klon 76 pt(i, k) = t(i, klev-k+1) 77 pq(i, k) = q(i, klev-k+1) 78 pres(i, k) = pplay(i, klev-k+1) 79 dp(i, k) = paprs(i, klev+1-k) - paprs(i, klev+1-k+1) 80 END DO 81 END DO 82 DO i = 1, klon 83 zgeom(i, klev) = rd*t(i, 1)/(0.5*(paprs(i,1)+pplay(i, & 84 1)))*(paprs(i,1)-pplay(i,1)) 85 END DO 86 DO k = 2, klev 87 DO i = 1, klon 88 zgeom(i, klev+1-k) = zgeom(i, klev+1-k+1) + rd*0.5*(t(i,k-1)+t(i,k))/ & 89 paprs(i, k)*(pplay(i,k-1)-pplay(i,k)) 90 END DO 91 END DO 92 93 CALL cmfmca(dtime, pres, dp, zgeom, pt, pq, cmfprt, cmfprs, ntop, nbas) 94 95 DO k = 1, klev 96 DO i = 1, klon 97 d_q(i, klev+1-k) = pq(i, k) - q(i, klev+1-k) 98 d_t(i, klev+1-k) = pt(i, k) - t(i, klev+1-k) 99 END DO 100 END DO 101 102 DO i = 1, klon 103 rain(i) = cmfprt(i)*rhoh2o 104 snow(i) = cmfprs(i)*rhoh2o 105 kbascm(i) = klev + 1 - nbas(i) 106 ktopcm(i) = klev + 1 - ntop(i) 107 END DO 108 109 IF (usekuo) THEN 110 CALL conkuo(dtime, paprs, pplay, t, q, conv_q, d_t_bis, d_q_bis, & 111 d_ql_bis, rneb_bis, rain_bis, snow_bis, ibas_bis, itop_bis) 112 DO k = 1, klev 59 113 DO i = 1, klon 60 rain(i) = 0.0 61 snow(i) = 0.0 62 kbascm(i) = 0 63 ktopcm(i) = 0 64 ENDDO 65 DO k = 1, klev 114 d_t(i, k) = d_t(i, k) + d_t_bis(i, k) 115 d_q(i, k) = d_q(i, k) + d_q_bis(i, k) 116 END DO 117 END DO 118 DO i = 1, klon 119 rain(i) = rain(i) + rain_bis(i) 120 snow(i) = snow(i) + snow_bis(i) 121 kbascm(i) = min(kbascm(i), ibas_bis(i)) 122 ktopcm(i) = max(ktopcm(i), itop_bis(i)) 123 END DO 124 DO k = 1, klev ! eau liquide convective est 125 DO i = 1, klon ! dispersee dans l'air 126 zlvdcp = rlvtt/rcpd/(1.0+rvtmp2*q(i,k)) 127 zlsdcp = rlstt/rcpd/(1.0+rvtmp2*q(i,k)) 128 zdelta = max(0., sign(1.,rtt-t(i,k))) 129 zz = d_ql_bis(i, k) ! re-evap. de l'eau liquide 130 zb = max(0.0, zz) 131 za = -max(0.0, zz)*(zlvdcp*(1.-zdelta)+zlsdcp*zdelta) 132 d_t(i, k) = d_t(i, k) + za 133 d_q(i, k) = d_q(i, k) + zb 134 END DO 135 END DO 136 END IF 137 138 RETURN 139 END SUBROUTINE conccm 140 SUBROUTINE cmfmca(deltat, p, dp, gz, tb, shb, cmfprt, cmfprs, cnt, cnb) 141 USE dimphy 142 IMPLICIT NONE 143 ! ----------------------------------------------------------------------- 144 ! Moist convective mass flux procedure: 145 ! If stratification is unstable to nonentraining parcel ascent, 146 ! complete an adjustment making use of a simple cloud model 147 148 ! Code generalized to allow specification of parcel ("updraft") 149 ! properties, as well as convective transport of an arbitrary 150 ! number of passive constituents (see cmrb array). 151 ! ----------------------------Code History------------------------------- 152 ! Original version: J. J. Hack, March 22, 1990 153 ! Standardized: J. Rosinski, June 1992 154 ! Reviewed: J. Hack, G. Taylor, August 1992 155 ! Adaptation au LMD: Z.X. Li, mars 1996 (reference: Hack 1994, 156 ! J. Geophys. Res. vol 99, D3, 5551-5568). J'ai 157 ! introduit les constantes et les fonctions thermo- 158 ! dynamiques du Centre Europeen. J'ai elimine le 159 ! re-indicage du code en esperant que cela pourra 160 ! simplifier la lecture et la comprehension. 161 ! ----------------------------------------------------------------------- 162 ! ym#include "dimensions.h" 163 ! ym#include "dimphy.h" 164 INTEGER pcnst ! nombre de traceurs passifs 165 PARAMETER (pcnst=1) 166 ! ------------------------------Arguments-------------------------------- 167 ! Input arguments 168 169 REAL deltat ! time step (seconds) 170 REAL p(klon, klev) ! pressure 171 REAL dp(klon, klev) ! delta-p 172 REAL gz(klon, klev) ! geopotential (a partir du sol) 173 174 REAL thtap(klon) ! PBL perturbation theta 175 REAL shp(klon) ! PBL perturbation specific humidity 176 REAL pblh(klon) ! PBL height (provided by PBL routine) 177 REAL cmrp(klon, pcnst) ! constituent perturbations in PBL 178 179 ! Updated arguments: 180 181 REAL tb(klon, klev) ! temperature (t bar) 182 REAL shb(klon, klev) ! specific humidity (sh bar) 183 REAL cmrb(klon, klev, pcnst) ! constituent mixing ratios (cmr bar) 184 185 ! Output arguments 186 187 REAL cmfdt(klon, klev) ! dT/dt due to moist convection 188 REAL cmfdq(klon, klev) ! dq/dt due to moist convection 189 REAL cmfmc(klon, klev) ! moist convection cloud mass flux 190 REAL cmfdqr(klon, klev) ! dq/dt due to convective rainout 191 REAL cmfsl(klon, klev) ! convective lw static energy flux 192 REAL cmflq(klon, klev) ! convective total water flux 193 REAL cmfprt(klon) ! convective precipitation rate 194 REAL cmfprs(klon) ! convective snowfall rate 195 REAL qc(klon, klev) ! dq/dt due to rainout terms 196 INTEGER cnt(klon) ! top level of convective activity 197 INTEGER cnb(klon) ! bottom level of convective activity 198 ! ------------------------------Parameters------------------------------- 199 REAL c0 ! rain water autoconversion coefficient 200 PARAMETER (c0=1.0E-4) 201 REAL dzmin ! minimum convective depth for precipitation 202 PARAMETER (dzmin=0.0) 203 REAL betamn ! minimum overshoot parameter 204 PARAMETER (betamn=0.10) 205 REAL cmftau ! characteristic adjustment time scale 206 PARAMETER (cmftau=3600.) 207 INTEGER limcnv ! top interface level limit for convection 208 PARAMETER (limcnv=1) 209 REAL tpmax ! maximum acceptable t perturbation (degrees C) 210 PARAMETER (tpmax=1.50) 211 REAL shpmax ! maximum acceptable q perturbation (g/g) 212 PARAMETER (shpmax=1.50E-3) 213 REAL tiny ! arbitrary small num used in transport estimates 214 PARAMETER (tiny=1.0E-36) 215 REAL eps ! convergence criteria (machine dependent) 216 PARAMETER (eps=1.0E-13) 217 REAL tmelt ! freezing point of water(req'd for rain vs snow) 218 PARAMETER (tmelt=273.15) 219 REAL ssfac ! supersaturation bound (detrained air) 220 PARAMETER (ssfac=1.001) 221 222 ! ---------------------------Local workspace----------------------------- 223 REAL gam(klon, klev) ! L/cp (d(qsat)/dT) 224 REAL sb(klon, klev) ! dry static energy (s bar) 225 REAL hb(klon, klev) ! moist static energy (h bar) 226 REAL shbs(klon, klev) ! sat. specific humidity (sh bar star) 227 REAL hbs(klon, klev) ! sat. moist static energy (h bar star) 228 REAL shbh(klon, klev+1) ! specific humidity on interfaces 229 REAL sbh(klon, klev+1) ! s bar on interfaces 230 REAL hbh(klon, klev+1) ! h bar on interfaces 231 REAL cmrh(klon, klev+1) ! interface constituent mixing ratio 232 REAL prec(klon) ! instantaneous total precipitation 233 REAL dzcld(klon) ! depth of convective layer (m) 234 REAL beta(klon) ! overshoot parameter (fraction) 235 REAL betamx ! local maximum on overshoot 236 REAL eta(klon) ! convective mass flux (kg/m^2 s) 237 REAL etagdt ! eta*grav*deltat 238 REAL cldwtr(klon) ! cloud water (mass) 239 REAL rnwtr(klon) ! rain water (mass) 240 REAL sc(klon) ! dry static energy ("in-cloud") 241 REAL shc(klon) ! specific humidity ("in-cloud") 242 REAL hc(klon) ! moist static energy ("in-cloud") 243 REAL cmrc(klon) ! constituent mix rat ("in-cloud") 244 REAL dq1(klon) ! shb convective change (lower lvl) 245 REAL dq2(klon) ! shb convective change (mid level) 246 REAL dq3(klon) ! shb convective change (upper lvl) 247 REAL ds1(klon) ! sb convective change (lower lvl) 248 REAL ds2(klon) ! sb convective change (mid level) 249 REAL ds3(klon) ! sb convective change (upper lvl) 250 REAL dcmr1(klon) ! cmrb convective change (lower lvl) 251 REAL dcmr2(klon) ! cmrb convective change (mid level) 252 REAL dcmr3(klon) ! cmrb convective change (upper lvl) 253 REAL flotab(klon) ! hc - hbs (mesure d'instabilite) 254 LOGICAL ldcum(klon) ! .true. si la convection existe 255 LOGICAL etagt0 ! true if eta > 0.0 256 REAL dt ! current 2 delta-t (model time step) 257 REAL cats ! modified characteristic adj. time 258 REAL rdt ! 1./dt 259 REAL qprime ! modified specific humidity pert. 260 REAL tprime ! modified thermal perturbation 261 REAL pblhgt ! bounded pbl height (max[pblh,1m]) 262 REAL fac1 ! intermediate scratch variable 263 REAL shprme ! intermediate specific humidity pert. 264 REAL qsattp ! saturation mixing ratio for 265 ! ! thermally perturbed PBL parcels 266 REAL dz ! local layer depth 267 REAL b1 ! bouyancy measure in detrainment lvl 268 REAL b2 ! bouyancy measure in condensation lvl 269 REAL g ! bounded vertical gradient of hb 270 REAL tmass ! total mass available for convective exchange 271 REAL denom ! intermediate scratch variable 272 REAL qtest1 ! used in negative q test (middle lvl) 273 REAL qtest2 ! used in negative q test (lower lvl) 274 REAL fslkp ! flux lw static energy (bot interface) 275 REAL fslkm ! flux lw static energy (top interface) 276 REAL fqlkp ! flux total water (bottom interface) 277 REAL fqlkm ! flux total water (top interface) 278 REAL botflx ! bottom constituent mixing ratio flux 279 REAL topflx ! top constituent mixing ratio flux 280 REAL efac1 ! ratio cmrb to convectively induced change (bot lvl) 281 REAL efac2 ! ratio cmrb to convectively induced change (mid lvl) 282 REAL efac3 ! ratio cmrb to convectively induced change (top lvl) 283 284 INTEGER i, k ! indices horizontal et vertical 285 INTEGER km1 ! k-1 (index offset) 286 INTEGER kp1 ! k+1 (index offset) 287 INTEGER m ! constituent index 288 INTEGER ktp ! temporary index used to track top 289 INTEGER is ! nombre de points a ajuster 290 291 REAL tmp1, tmp2, tmp3, tmp4 292 REAL zx_t, zx_p, zx_q, zx_qs, zx_gam 293 REAL zcor, zdelta, zcvm5 294 295 REAL qhalf, sh1, sh2, shbs1, shbs2 296 include "YOMCST.h" 297 include "YOETHF.h" 298 include "FCTTRE.h" 299 qhalf(sh1, sh2, shbs1, shbs2) = min(max(sh1,sh2), & 300 (shbs2*sh1+shbs1*sh2)/(shbs1+shbs2)) 301 302 ! ----------------------------------------------------------------------- 303 ! pas de traceur pour l'instant 304 DO m = 1, pcnst 305 DO k = 1, klev 66 306 DO i = 1, klon 67 d_t(i,k) = 0.0 68 d_q(i,k) = 0.0 69 ENDDO 70 ENDDO 71 c 72 c preparer les variables d'entree (attention: l'ordre des niveaux 73 c verticaux augmente du haut vers le bas) 74 DO k = 1, klev 307 cmrb(i, k, m) = 0.0 308 END DO 309 END DO 310 END DO 311 312 ! Les perturbations de la couche limite sont zero pour l'instant 313 314 DO m = 1, pcnst 315 DO i = 1, klon 316 cmrp(i, m) = 0.0 317 END DO 318 END DO 319 DO i = 1, klon 320 thtap(i) = 0.0 321 shp(i) = 0.0 322 pblh(i) = 1.0 323 END DO 324 325 ! Ensure that characteristic adjustment time scale (cmftau) assumed 326 ! in estimate of eta isn't smaller than model time scale (deltat) 327 328 dt = deltat 329 cats = max(dt, cmftau) 330 rdt = 1.0/dt 331 332 ! Compute sb,hb,shbs,hbs 333 334 DO k = 1, klev 335 DO i = 1, klon 336 zx_t = tb(i, k) 337 zx_p = p(i, k) 338 zx_q = shb(i, k) 339 zdelta = max(0., sign(1.,rtt-zx_t)) 340 zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta 341 zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zx_q) 342 zx_qs = r2es*foeew(zx_t, zdelta)/zx_p 343 zx_qs = min(0.5, zx_qs) 344 zcor = 1./(1.-retv*zx_qs) 345 zx_qs = zx_qs*zcor 346 zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) 347 shbs(i, k) = zx_qs 348 gam(i, k) = zx_gam 349 END DO 350 END DO 351 352 DO k = limcnv, klev 353 DO i = 1, klon 354 sb(i, k) = rcpd*tb(i, k) + gz(i, k) 355 hb(i, k) = sb(i, k) + rlvtt*shb(i, k) 356 hbs(i, k) = sb(i, k) + rlvtt*shbs(i, k) 357 END DO 358 END DO 359 360 ! Compute sbh, shbh 361 362 DO k = limcnv + 1, klev 363 km1 = k - 1 364 DO i = 1, klon 365 sbh(i, k) = 0.5*(sb(i,km1)+sb(i,k)) 366 shbh(i, k) = qhalf(shb(i,km1), shb(i,k), shbs(i,km1), shbs(i,k)) 367 hbh(i, k) = sbh(i, k) + rlvtt*shbh(i, k) 368 END DO 369 END DO 370 371 ! Specify properties at top of model (not used, but filling anyway) 372 373 DO i = 1, klon 374 sbh(i, limcnv) = sb(i, limcnv) 375 shbh(i, limcnv) = shb(i, limcnv) 376 hbh(i, limcnv) = hb(i, limcnv) 377 END DO 378 379 ! Zero vertically independent control, tendency & diagnostic arrays 380 381 DO i = 1, klon 382 prec(i) = 0.0 383 dzcld(i) = 0.0 384 cnb(i) = 0 385 cnt(i) = klev + 1 386 END DO 387 388 DO k = 1, klev 389 DO i = 1, klon 390 cmfdt(i, k) = 0. 391 cmfdq(i, k) = 0. 392 cmfdqr(i, k) = 0. 393 cmfmc(i, k) = 0. 394 cmfsl(i, k) = 0. 395 cmflq(i, k) = 0. 396 END DO 397 END DO 398 399 ! Begin moist convective mass flux adjustment procedure. 400 ! Formalism ensures that negative cloud liquid water can never occur 401 402 DO k = klev - 1, limcnv + 1, -1 403 km1 = k - 1 404 kp1 = k + 1 405 DO i = 1, klon 406 eta(i) = 0.0 407 beta(i) = 0.0 408 ds1(i) = 0.0 409 ds2(i) = 0.0 410 ds3(i) = 0.0 411 dq1(i) = 0.0 412 dq2(i) = 0.0 413 dq3(i) = 0.0 414 415 ! Specification of "cloud base" conditions 416 417 qprime = 0.0 418 tprime = 0.0 419 420 ! Assign tprime within the PBL to be proportional to the quantity 421 ! thtap (which will be bounded by tpmax), passed to this routine by 422 ! the PBL routine. Don't allow perturbation to produce a dry 423 ! adiabatically unstable parcel. Assign qprime within the PBL to be 424 ! an appropriately modified value of the quantity shp (which will be 425 ! bounded by shpmax) passed to this routine by the PBL routine. The 426 ! quantity qprime should be less than the local saturation value 427 ! (qsattp=qsat[t+tprime,p]). In both cases, thtap and shp are 428 ! linearly reduced toward zero as the PBL top is approached. 429 430 pblhgt = max(pblh(i), 1.0) 431 IF (gz(i,kp1)/rg<=pblhgt .AND. dzcld(i)==0.0) THEN 432 fac1 = max(0.0, 1.0-gz(i,kp1)/rg/pblhgt) 433 tprime = min(thtap(i), tpmax)*fac1 434 qsattp = shbs(i, kp1) + rcpd/rlvtt*gam(i, kp1)*tprime 435 shprme = min(min(shp(i),shpmax)*fac1, max(qsattp-shb(i,kp1),0.0)) 436 qprime = max(qprime, shprme) 437 ELSE 438 tprime = 0.0 439 qprime = 0.0 440 END IF 441 442 ! Specify "updraft" (in-cloud) thermodynamic properties 443 444 sc(i) = sb(i, kp1) + rcpd*tprime 445 shc(i) = shb(i, kp1) + qprime 446 hc(i) = sc(i) + rlvtt*shc(i) 447 flotab(i) = hc(i) - hbs(i, k) 448 dz = dp(i, k)*rd*tb(i, k)/rg/p(i, k) 449 IF (flotab(i)>0.0) THEN 450 dzcld(i) = dzcld(i) + dz 451 ELSE 452 dzcld(i) = 0.0 453 END IF 454 END DO 455 456 ! Check on moist convective instability 457 458 is = 0 459 DO i = 1, klon 460 IF (flotab(i)>0.0) THEN 461 ldcum(i) = .TRUE. 462 is = is + 1 463 ELSE 464 ldcum(i) = .FALSE. 465 END IF 466 END DO 467 468 IF (is==0) THEN 75 469 DO i = 1, klon 76 pt(i,k) = t(i,klev-k+1) 77 pq(i,k) = q(i,klev-k+1) 78 pres(i,k) = pplay(i,klev-k+1) 79 dp(i,k) = paprs(i,klev+1-k)-paprs(i,klev+1-k+1) 80 ENDDO 81 ENDDO 470 dzcld(i) = 0.0 471 END DO 472 GO TO 70 473 END IF 474 475 ! Current level just below top level => no overshoot 476 477 IF (k<=limcnv+1) THEN 82 478 DO i = 1, klon 83 zgeom(i,klev) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) 84 . * (paprs(i,1)-pplay(i,1)) 85 ENDDO 86 DO k = 2, klev 479 IF (ldcum(i)) THEN 480 cldwtr(i) = sb(i, k) - sc(i) + flotab(i)/(1.0+gam(i,k)) 481 cldwtr(i) = max(0.0, cldwtr(i)) 482 beta(i) = 0.0 483 END IF 484 END DO 485 GO TO 20 486 END IF 487 488 ! First guess at overshoot parameter using crude buoyancy closure 489 ! 10% overshoot assumed as a minimum and 1-c0*dz maximum to start 490 ! If pre-existing supersaturation in detrainment layer, beta=0 491 ! cldwtr is temporarily equal to RLVTT*l (l=> liquid water) 492 493 DO i = 1, klon 494 IF (ldcum(i)) THEN 495 cldwtr(i) = sb(i, k) - sc(i) + flotab(i)/(1.0+gam(i,k)) 496 cldwtr(i) = max(0.0, cldwtr(i)) 497 betamx = 1.0 - c0*max(0.0, (dzcld(i)-dzmin)) 498 b1 = (hc(i)-hbs(i,km1))*dp(i, km1) 499 b2 = (hc(i)-hbs(i,k))*dp(i, k) 500 beta(i) = max(betamn, min(betamx,1.0+b1/b2)) 501 IF (hbs(i,km1)<=hb(i,km1)) beta(i) = 0.0 502 END IF 503 END DO 504 505 ! Bound maximum beta to ensure physically realistic solutions 506 507 ! First check constrains beta so that eta remains positive 508 ! (assuming that eta is already positive for beta equal zero) 509 ! La premiere contrainte de beta est que le flux eta doit etre positif. 510 511 DO i = 1, klon 512 IF (ldcum(i)) THEN 513 tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) - & 514 (hbh(i,kp1)-hc(i))*dp(i, k)/dp(i, kp1) 515 tmp2 = (1.0+gam(i,k))*(sc(i)-sbh(i,k)) 516 IF ((beta(i)*tmp2-tmp1)>0.0) THEN 517 betamx = 0.99*(tmp1/tmp2) 518 beta(i) = max(0.0, min(betamx,beta(i))) 519 END IF 520 521 ! Second check involves supersaturation of "detrainment layer" 522 ! small amount of supersaturation acceptable (by ssfac factor) 523 ! La 2e contrainte est que la convection ne doit pas sursaturer 524 ! la "detrainment layer", Neanmoins, une petite sursaturation 525 ! est acceptee (facteur ssfac). 526 527 IF (hb(i,km1)<hbs(i,km1)) THEN 528 tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) - & 529 (hbh(i,kp1)-hc(i))*dp(i, k)/dp(i, kp1) 530 tmp1 = tmp1/dp(i, k) 531 tmp2 = gam(i, km1)*(sbh(i,k)-sc(i)+cldwtr(i)) - hbh(i, k) + hc(i) - & 532 sc(i) + sbh(i, k) 533 tmp3 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i, k) 534 tmp4 = (dt/cats)*(hc(i)-hbs(i,k))*tmp2/(dp(i,km1)*(hbs(i,km1)-hb(i, & 535 km1))) + tmp3 536 IF ((beta(i)*tmp4-tmp1)>0.0) THEN 537 betamx = ssfac*(tmp1/tmp4) 538 beta(i) = max(0.0, min(betamx,beta(i))) 539 END IF 540 ELSE 541 beta(i) = 0.0 542 END IF 543 544 ! Third check to avoid introducing 2 delta x thermodynamic 545 ! noise in the vertical ... constrain adjusted h (or theta e) 546 ! so that the adjustment doesn't contribute to "kinks" in h 547 548 g = min(0.0, hb(i,k)-hb(i,km1)) 549 tmp3 = (hb(i,k)-hb(i,km1)-g)*(cats/dt)/(hc(i)-hbs(i,k)) 550 tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) - & 551 (hbh(i,kp1)-hc(i))*dp(i, k)/dp(i, kp1) 552 tmp1 = tmp1/dp(i, k) 553 tmp1 = tmp3*tmp1 + (hc(i)-hbh(i,kp1))/dp(i, k) 554 tmp2 = tmp3*(1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i, k) + & 555 (hc(i)-hbh(i,k)-cldwtr(i))*(1.0/dp(i,k)+1.0/dp(i,kp1)) 556 IF ((beta(i)*tmp2-tmp1)>0.0) THEN 557 betamx = 0.0 558 IF (tmp2/=0.0) betamx = tmp1/tmp2 559 beta(i) = max(0.0, min(betamx,beta(i))) 560 END IF 561 END IF 562 END DO 563 564 ! Calculate mass flux required for stabilization. 565 566 ! Ensure that the convective mass flux, eta, is positive by 567 ! setting negative values of eta to zero.. 568 ! Ensure that estimated mass flux cannot move more than the 569 ! minimum of total mass contained in either layer k or layer k+1. 570 ! Also test for other pathological cases that result in non- 571 ! physical states and adjust eta accordingly. 572 573 20 CONTINUE 574 DO i = 1, klon 575 IF (ldcum(i)) THEN 576 beta(i) = max(0.0, beta(i)) 577 tmp1 = hc(i) - hbs(i, k) 578 tmp2 = ((1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i))-beta(i)*(1.0+gam( & 579 i,k))*(sc(i)-sbh(i,k)))/dp(i, k) - (hbh(i,kp1)-hc(i))/dp(i, kp1) 580 eta(i) = tmp1/(tmp2*rg*cats) 581 tmass = min(dp(i,k), dp(i,kp1))/rg 582 IF (eta(i)>tmass*rdt .OR. eta(i)<=0.0) eta(i) = 0.0 583 584 ! Check on negative q in top layer (bound beta) 585 586 IF (shc(i)-shbh(i,k)<0.0 .AND. beta(i)*eta(i)/=0.0) THEN 587 denom = eta(i)*rg*dt*(shc(i)-shbh(i,k))/dp(i, km1) 588 beta(i) = max(0.0, min(-0.999*shb(i,km1)/denom,beta(i))) 589 END IF 590 591 ! Check on negative q in middle layer (zero eta) 592 593 qtest1 = shb(i, k) + eta(i)*rg*dt*((shc(i)-shbh(i, & 594 kp1))-(1.0-beta(i))*cldwtr(i)/rlvtt-beta(i)*(shc(i)-shbh(i, & 595 k)))/dp(i, k) 596 IF (qtest1<=0.0) eta(i) = 0.0 597 598 ! Check on negative q in lower layer (bound eta) 599 600 fac1 = -(shbh(i,kp1)-shc(i))/dp(i, kp1) 601 qtest2 = shb(i, kp1) - eta(i)*rg*dt*fac1 602 IF (qtest2<0.0) THEN 603 eta(i) = 0.99*shb(i, kp1)/(rg*dt*fac1) 604 END IF 605 END IF 606 END DO 607 608 609 ! Calculate cloud water, rain water, and thermodynamic changes 610 611 DO i = 1, klon 612 IF (ldcum(i)) THEN 613 etagdt = eta(i)*rg*dt 614 cldwtr(i) = etagdt*cldwtr(i)/rlvtt/rg 615 rnwtr(i) = (1.0-beta(i))*cldwtr(i) 616 ds1(i) = etagdt*(sbh(i,kp1)-sc(i))/dp(i, kp1) 617 dq1(i) = etagdt*(shbh(i,kp1)-shc(i))/dp(i, kp1) 618 ds2(i) = (etagdt*(sc(i)-sbh(i,kp1))+rlvtt*rg*cldwtr(i)-beta(i)*etagdt & 619 *(sc(i)-sbh(i,k)))/dp(i, k) 620 dq2(i) = (etagdt*(shc(i)-shbh(i,kp1))-rg*rnwtr(i)-beta(i)*etagdt*(shc & 621 (i)-shbh(i,k)))/dp(i, k) 622 ds3(i) = beta(i)*(etagdt*(sc(i)-sbh(i,k))-rlvtt*rg*cldwtr(i))/dp(i, & 623 km1) 624 dq3(i) = beta(i)*etagdt*(shc(i)-shbh(i,k))/dp(i, km1) 625 626 ! Isolate convective fluxes for later diagnostics 627 628 fslkp = eta(i)*(sc(i)-sbh(i,kp1)) 629 fslkm = beta(i)*(eta(i)*(sc(i)-sbh(i,k))-rlvtt*cldwtr(i)*rdt) 630 fqlkp = eta(i)*(shc(i)-shbh(i,kp1)) 631 fqlkm = beta(i)*eta(i)*(shc(i)-shbh(i,k)) 632 633 634 ! Update thermodynamic profile (update sb, hb, & hbs later) 635 636 tb(i, kp1) = tb(i, kp1) + ds1(i)/rcpd 637 tb(i, k) = tb(i, k) + ds2(i)/rcpd 638 tb(i, km1) = tb(i, km1) + ds3(i)/rcpd 639 shb(i, kp1) = shb(i, kp1) + dq1(i) 640 shb(i, k) = shb(i, k) + dq2(i) 641 shb(i, km1) = shb(i, km1) + dq3(i) 642 prec(i) = prec(i) + rnwtr(i)/rhoh2o 643 644 ! Update diagnostic information for final budget 645 ! Tracking temperature & specific humidity tendencies, 646 ! rainout term, convective mass flux, convective liquid 647 ! water static energy flux, and convective total water flux 648 649 cmfdt(i, kp1) = cmfdt(i, kp1) + ds1(i)/rcpd*rdt 650 cmfdt(i, k) = cmfdt(i, k) + ds2(i)/rcpd*rdt 651 cmfdt(i, km1) = cmfdt(i, km1) + ds3(i)/rcpd*rdt 652 cmfdq(i, kp1) = cmfdq(i, kp1) + dq1(i)*rdt 653 cmfdq(i, k) = cmfdq(i, k) + dq2(i)*rdt 654 cmfdq(i, km1) = cmfdq(i, km1) + dq3(i)*rdt 655 cmfdqr(i, k) = cmfdqr(i, k) + (rg*rnwtr(i)/dp(i,k))*rdt 656 cmfmc(i, kp1) = cmfmc(i, kp1) + eta(i) 657 cmfmc(i, k) = cmfmc(i, k) + beta(i)*eta(i) 658 cmfsl(i, kp1) = cmfsl(i, kp1) + fslkp 659 cmfsl(i, k) = cmfsl(i, k) + fslkm 660 cmflq(i, kp1) = cmflq(i, kp1) + rlvtt*fqlkp 661 cmflq(i, k) = cmflq(i, k) + rlvtt*fqlkm 662 qc(i, k) = (rg*rnwtr(i)/dp(i,k))*rdt 663 END IF 664 END DO 665 666 ! Next, convectively modify passive constituents 667 668 DO m = 1, pcnst 87 669 DO i = 1, klon 88 zgeom(i,klev+1-k) = zgeom(i,klev+1-k+1) 89 . + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) 90 . * (pplay(i,k-1)-pplay(i,k)) 91 ENDDO 92 ENDDO 93 c 94 CALL cmfmca(dtime, pres, dp, zgeom, pt, pq, 95 $ cmfprt, cmfprs, ntop, nbas) 96 C 97 DO k = 1, klev 98 DO i = 1, klon 99 d_q(i,klev+1-k) = pq(i,k) - q(i,klev+1-k) 100 d_t(i,klev+1-k) = pt(i,k) - t(i,klev+1-k) 101 ENDDO 102 ENDDO 103 c 104 DO i = 1, klon 105 rain(i) = cmfprt(i) * rhoh2o 106 snow(i) = cmfprs(i) * rhoh2o 107 kbascm(i) = klev+1 - nbas(i) 108 ktopcm(i) = klev+1 - ntop(i) 109 ENDDO 110 c 111 IF (usekuo) THEN 112 CALL conkuo(dtime, paprs, pplay, t, q, conv_q, 113 s d_t_bis, d_q_bis, d_ql_bis, rneb_bis, 114 s rain_bis, snow_bis, ibas_bis, itop_bis) 115 DO k = 1, klev 116 DO i = 1, klon 117 d_t(i,k) = d_t(i,k) + d_t_bis(i,k) 118 d_q(i,k) = d_q(i,k) + d_q_bis(i,k) 119 ENDDO 120 ENDDO 121 DO i = 1, klon 122 rain(i) = rain(i) + rain_bis(i) 123 snow(i) = snow(i) + snow_bis(i) 124 kbascm(i) = MIN(kbascm(i),ibas_bis(i)) 125 ktopcm(i) = MAX(ktopcm(i),itop_bis(i)) 126 ENDDO 127 DO k = 1, klev ! eau liquide convective est 128 DO i = 1, klon ! dispersee dans l'air 129 zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q(i,k)) 130 zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q(i,k)) 131 zdelta = MAX(0.,SIGN(1.,RTT-t(i,k))) 132 zz = d_ql_bis(i,k) ! re-evap. de l'eau liquide 133 zb = MAX(0.0,zz) 134 za = - MAX(0.0,zz) * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta) 135 d_t(i,k) = d_t(i,k) + za 136 d_q(i,k) = d_q(i,k) + zb 137 ENDDO 138 ENDDO 139 ENDIF 140 c 141 RETURN 142 END 143 SUBROUTINE cmfmca(deltat, p, dp, gz, 144 $ tb, shb, 145 $ cmfprt, cmfprs, cnt, cnb) 146 USE dimphy 147 IMPLICIT none 148 C----------------------------------------------------------------------- 149 C Moist convective mass flux procedure: 150 C If stratification is unstable to nonentraining parcel ascent, 151 C complete an adjustment making use of a simple cloud model 152 C 153 C Code generalized to allow specification of parcel ("updraft") 154 C properties, as well as convective transport of an arbitrary 155 C number of passive constituents (see cmrb array). 156 C----------------------------Code History------------------------------- 157 C Original version: J. J. Hack, March 22, 1990 158 C Standardized: J. Rosinski, June 1992 159 C Reviewed: J. Hack, G. Taylor, August 1992 160 c Adaptation au LMD: Z.X. Li, mars 1996 (reference: Hack 1994, 161 c J. Geophys. Res. vol 99, D3, 5551-5568). J'ai 162 c introduit les constantes et les fonctions thermo- 163 c dynamiques du Centre Europeen. J'ai elimine le 164 c re-indicage du code en esperant que cela pourra 165 c simplifier la lecture et la comprehension. 166 C----------------------------------------------------------------------- 167 cym#include "dimensions.h" 168 cym#include "dimphy.h" 169 INTEGER pcnst ! nombre de traceurs passifs 170 PARAMETER (pcnst=1) 171 C------------------------------Arguments-------------------------------- 172 C Input arguments 173 C 174 REAL deltat ! time step (seconds) 175 REAL p(klon,klev) ! pressure 176 REAL dp(klon,klev) ! delta-p 177 REAL gz(klon,klev) ! geopotential (a partir du sol) 178 c 179 REAL thtap(klon) ! PBL perturbation theta 180 REAL shp(klon) ! PBL perturbation specific humidity 181 REAL pblh(klon) ! PBL height (provided by PBL routine) 182 REAL cmrp(klon,pcnst) ! constituent perturbations in PBL 183 c 184 c Updated arguments: 185 c 186 REAL tb(klon,klev) ! temperature (t bar) 187 REAL shb(klon,klev) ! specific humidity (sh bar) 188 REAL cmrb(klon,klev,pcnst) ! constituent mixing ratios (cmr bar) 189 C 190 C Output arguments 191 C 192 REAL cmfdt(klon,klev) ! dT/dt due to moist convection 193 REAL cmfdq(klon,klev) ! dq/dt due to moist convection 194 REAL cmfmc(klon,klev ) ! moist convection cloud mass flux 195 REAL cmfdqr(klon,klev) ! dq/dt due to convective rainout 196 REAL cmfsl(klon,klev ) ! convective lw static energy flux 197 REAL cmflq(klon,klev ) ! convective total water flux 198 REAL cmfprt(klon) ! convective precipitation rate 199 REAL cmfprs(klon) ! convective snowfall rate 200 REAL qc(klon,klev) ! dq/dt due to rainout terms 201 INTEGER cnt(klon) ! top level of convective activity 202 INTEGER cnb(klon) ! bottom level of convective activity 203 C------------------------------Parameters------------------------------- 204 REAL c0 ! rain water autoconversion coefficient 205 PARAMETER (c0=1.0e-4) 206 REAL dzmin ! minimum convective depth for precipitation 207 PARAMETER (dzmin=0.0) 208 REAL betamn ! minimum overshoot parameter 209 PARAMETER (betamn=0.10) 210 REAL cmftau ! characteristic adjustment time scale 211 PARAMETER (cmftau=3600.) 212 INTEGER limcnv ! top interface level limit for convection 213 PARAMETER (limcnv=1) 214 REAL tpmax ! maximum acceptable t perturbation (degrees C) 215 PARAMETER (tpmax=1.50) 216 REAL shpmax ! maximum acceptable q perturbation (g/g) 217 PARAMETER (shpmax=1.50e-3) 218 REAL tiny ! arbitrary small num used in transport estimates 219 PARAMETER (tiny=1.0e-36) 220 REAL eps ! convergence criteria (machine dependent) 221 PARAMETER (eps=1.0e-13) 222 REAL tmelt ! freezing point of water(req'd for rain vs snow) 223 PARAMETER (tmelt=273.15) 224 REAL ssfac ! supersaturation bound (detrained air) 225 PARAMETER (ssfac=1.001) 226 C 227 C---------------------------Local workspace----------------------------- 228 REAL gam(klon,klev) ! L/cp (d(qsat)/dT) 229 REAL sb(klon,klev) ! dry static energy (s bar) 230 REAL hb(klon,klev) ! moist static energy (h bar) 231 REAL shbs(klon,klev) ! sat. specific humidity (sh bar star) 232 REAL hbs(klon,klev) ! sat. moist static energy (h bar star) 233 REAL shbh(klon,klev+1) ! specific humidity on interfaces 234 REAL sbh(klon,klev+1) ! s bar on interfaces 235 REAL hbh(klon,klev+1) ! h bar on interfaces 236 REAL cmrh(klon,klev+1) ! interface constituent mixing ratio 237 REAL prec(klon) ! instantaneous total precipitation 238 REAL dzcld(klon) ! depth of convective layer (m) 239 REAL beta(klon) ! overshoot parameter (fraction) 240 REAL betamx ! local maximum on overshoot 241 REAL eta(klon) ! convective mass flux (kg/m^2 s) 242 REAL etagdt ! eta*grav*deltat 243 REAL cldwtr(klon) ! cloud water (mass) 244 REAL rnwtr(klon) ! rain water (mass) 245 REAL sc (klon) ! dry static energy ("in-cloud") 246 REAL shc (klon) ! specific humidity ("in-cloud") 247 REAL hc (klon) ! moist static energy ("in-cloud") 248 REAL cmrc(klon) ! constituent mix rat ("in-cloud") 249 REAL dq1(klon) ! shb convective change (lower lvl) 250 REAL dq2(klon) ! shb convective change (mid level) 251 REAL dq3(klon) ! shb convective change (upper lvl) 252 REAL ds1(klon) ! sb convective change (lower lvl) 253 REAL ds2(klon) ! sb convective change (mid level) 254 REAL ds3(klon) ! sb convective change (upper lvl) 255 REAL dcmr1(klon) ! cmrb convective change (lower lvl) 256 REAL dcmr2(klon) ! cmrb convective change (mid level) 257 REAL dcmr3(klon) ! cmrb convective change (upper lvl) 258 REAL flotab(klon) ! hc - hbs (mesure d'instabilite) 259 LOGICAL ldcum(klon) ! .true. si la convection existe 260 LOGICAL etagt0 ! true if eta > 0.0 261 REAL dt ! current 2 delta-t (model time step) 262 REAL cats ! modified characteristic adj. time 263 REAL rdt ! 1./dt 264 REAL qprime ! modified specific humidity pert. 265 REAL tprime ! modified thermal perturbation 266 REAL pblhgt ! bounded pbl height (max[pblh,1m]) 267 REAL fac1 ! intermediate scratch variable 268 REAL shprme ! intermediate specific humidity pert. 269 REAL qsattp ! saturation mixing ratio for 270 C ! thermally perturbed PBL parcels 271 REAL dz ! local layer depth 272 REAL b1 ! bouyancy measure in detrainment lvl 273 REAL b2 ! bouyancy measure in condensation lvl 274 REAL g ! bounded vertical gradient of hb 275 REAL tmass ! total mass available for convective exchange 276 REAL denom ! intermediate scratch variable 277 REAL qtest1! used in negative q test (middle lvl) 278 REAL qtest2! used in negative q test (lower lvl) 279 REAL fslkp ! flux lw static energy (bot interface) 280 REAL fslkm ! flux lw static energy (top interface) 281 REAL fqlkp ! flux total water (bottom interface) 282 REAL fqlkm ! flux total water (top interface) 283 REAL botflx! bottom constituent mixing ratio flux 284 REAL topflx! top constituent mixing ratio flux 285 REAL efac1 ! ratio cmrb to convectively induced change (bot lvl) 286 REAL efac2 ! ratio cmrb to convectively induced change (mid lvl) 287 REAL efac3 ! ratio cmrb to convectively induced change (top lvl) 288 c 289 INTEGER i,k ! indices horizontal et vertical 290 INTEGER km1 ! k-1 (index offset) 291 INTEGER kp1 ! k+1 (index offset) 292 INTEGER m ! constituent index 293 INTEGER ktp ! temporary index used to track top 294 INTEGER is ! nombre de points a ajuster 295 C 296 REAL tmp1, tmp2, tmp3, tmp4 297 REAL zx_t, zx_p, zx_q, zx_qs, zx_gam 298 REAL zcor, zdelta, zcvm5 299 C 300 REAL qhalf, sh1, sh2, shbs1, shbs2 301 #include "YOMCST.h" 302 #include "YOETHF.h" 303 #include "FCTTRE.h" 304 qhalf(sh1,sh2,shbs1,shbs2) = MIN(MAX(sh1,sh2), 305 $ (shbs2*sh1 + shbs1*sh2)/(shbs1+shbs2)) 306 C 307 C----------------------------------------------------------------------- 308 c pas de traceur pour l'instant 309 DO m = 1, pcnst 310 DO k = 1, klev 311 DO i = 1, klon 312 cmrb(i,k,m) = 0.0 313 ENDDO 314 ENDDO 315 ENDDO 316 c 317 c Les perturbations de la couche limite sont zero pour l'instant 318 c 319 DO m = 1, pcnst 320 DO i = 1, klon 321 cmrp(i,m) = 0.0 322 ENDDO 323 ENDDO 324 DO i = 1, klon 325 thtap(i) = 0.0 326 shp(i) = 0.0 327 pblh(i) = 1.0 328 ENDDO 329 C 330 C Ensure that characteristic adjustment time scale (cmftau) assumed 331 C in estimate of eta isn't smaller than model time scale (deltat) 332 C 333 dt = deltat 334 cats = MAX(dt,cmftau) 335 rdt = 1.0/dt 336 C 337 C Compute sb,hb,shbs,hbs 338 C 339 DO k = 1, klev 340 DO i = 1, klon 341 zx_t = tb(i,k) 342 zx_p = p(i,k) 343 zx_q = shb(i,k) 344 zdelta=MAX(0.,SIGN(1.,RTT-zx_t)) 345 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta 346 zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q) 347 zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p 348 zx_qs=MIN(0.5,zx_qs) 349 zcor=1./(1.-retv*zx_qs) 350 zx_qs=zx_qs*zcor 351 zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor) 352 shbs(i,k) = zx_qs 353 gam(i,k) = zx_gam 354 ENDDO 355 ENDDO 356 C 357 DO k=limcnv,klev 358 DO i=1,klon 359 sb (i,k) = RCPD*tb(i,k) + gz(i,k) 360 hb (i,k) = sb(i,k) + RLVTT*shb(i,k) 361 hbs(i,k) = sb(i,k) + RLVTT*shbs(i,k) 362 ENDDO 363 ENDDO 364 C 365 C Compute sbh, shbh 366 C 367 DO k=limcnv+1,klev 368 km1 = k - 1 369 DO i=1,klon 370 sbh (i,k) =0.5*(sb(i,km1) + sb(i,k)) 371 shbh(i,k) =qhalf(shb(i,km1),shb(i,k),shbs(i,km1),shbs(i,k)) 372 hbh (i,k) =sbh(i,k) + RLVTT*shbh(i,k) 373 ENDDO 374 ENDDO 375 C 376 C Specify properties at top of model (not used, but filling anyway) 377 C 378 DO i=1,klon 379 sbh (i,limcnv) = sb(i,limcnv) 380 shbh(i,limcnv) = shb(i,limcnv) 381 hbh (i,limcnv) = hb(i,limcnv) 382 ENDDO 383 C 384 C Zero vertically independent control, tendency & diagnostic arrays 385 C 386 DO i=1,klon 387 prec(i) = 0.0 388 dzcld(i) = 0.0 389 cnb(i) = 0 390 cnt(i) = klev+1 391 ENDDO 392 393 DO k = 1, klev 394 DO i = 1,klon 395 cmfdt(i,k) = 0. 396 cmfdq(i,k) = 0. 397 cmfdqr(i,k) = 0. 398 cmfmc(i,k) = 0. 399 cmfsl(i,k) = 0. 400 cmflq(i,k) = 0. 401 ENDDO 402 ENDDO 403 C 404 C Begin moist convective mass flux adjustment procedure. 405 C Formalism ensures that negative cloud liquid water can never occur 406 C 407 DO 70 k=klev-1,limcnv+1,-1 408 km1 = k - 1 409 kp1 = k + 1 410 DO 10 i=1,klon 411 eta (i) = 0.0 412 beta (i) = 0.0 413 ds1 (i) = 0.0 414 ds2 (i) = 0.0 415 ds3 (i) = 0.0 416 dq1 (i) = 0.0 417 dq2 (i) = 0.0 418 dq3 (i) = 0.0 419 C 420 C Specification of "cloud base" conditions 421 C 422 qprime = 0.0 423 tprime = 0.0 424 C 425 C Assign tprime within the PBL to be proportional to the quantity 426 C thtap (which will be bounded by tpmax), passed to this routine by 427 C the PBL routine. Don't allow perturbation to produce a dry 428 C adiabatically unstable parcel. Assign qprime within the PBL to be 429 C an appropriately modified value of the quantity shp (which will be 430 C bounded by shpmax) passed to this routine by the PBL routine. The 431 C quantity qprime should be less than the local saturation value 432 C (qsattp=qsat[t+tprime,p]). In both cases, thtap and shp are 433 C linearly reduced toward zero as the PBL top is approached. 434 C 435 pblhgt = MAX(pblh(i),1.0) 436 IF (gz(i,kp1)/RG.LE.pblhgt .AND. dzcld(i).EQ.0.0) THEN 437 fac1 = MAX(0.0,1.0-gz(i,kp1)/RG/pblhgt) 438 tprime = MIN(thtap(i),tpmax)*fac1 439 qsattp = shbs(i,kp1) + RCPD/RLVTT*gam(i,kp1)*tprime 440 shprme = MIN(MIN(shp(i),shpmax)*fac1, 441 $ MAX(qsattp-shb(i,kp1),0.0)) 442 qprime = MAX(qprime,shprme) 443 ELSE 444 tprime = 0.0 445 qprime = 0.0 446 ENDIF 447 C 448 C Specify "updraft" (in-cloud) thermodynamic properties 449 C 450 sc (i) = sb (i,kp1) + RCPD*tprime 451 shc(i) = shb(i,kp1) + qprime 452 hc (i) = sc (i ) + RLVTT*shc(i) 453 flotab(i) = hc(i) - hbs(i,k) 454 dz = dp(i,k)*RD*tb(i,k)/RG/p(i,k) 455 IF (flotab(i).gt.0.0) THEN 456 dzcld(i) = dzcld(i) + dz 457 ELSE 458 dzcld(i) = 0.0 459 ENDIF 460 10 CONTINUE 461 C 462 C Check on moist convective instability 463 C 464 is = 0 465 DO i = 1, klon 466 IF (flotab(i).GT.0.0) THEN 467 ldcum(i) = .TRUE. 468 is = is + 1 469 ELSE 470 ldcum(i) = .FALSE. 471 ENDIF 472 ENDDO 473 C 474 IF (is.EQ.0) THEN 475 DO i=1,klon 476 dzcld(i) = 0.0 477 ENDDO 478 GOTO 70 479 ENDIF 480 C 481 C Current level just below top level => no overshoot 482 C 483 IF (k.le.limcnv+1) THEN 484 DO i=1,klon 485 IF (ldcum(i)) THEN 486 cldwtr(i) = sb(i,k)-sc(i)+flotab(i)/(1.0+gam(i,k)) 487 cldwtr(i) = MAX(0.0,cldwtr(i)) 488 beta(i) = 0.0 489 ENDIF 490 ENDDO 491 GOTO 20 492 ENDIF 493 C 494 C First guess at overshoot parameter using crude buoyancy closure 495 C 10% overshoot assumed as a minimum and 1-c0*dz maximum to start 496 C If pre-existing supersaturation in detrainment layer, beta=0 497 C cldwtr is temporarily equal to RLVTT*l (l=> liquid water) 498 C 499 DO i=1,klon 500 IF (ldcum(i)) THEN 501 cldwtr(i) = sb(i,k)-sc(i)+flotab(i)/(1.0+gam(i,k)) 502 cldwtr(i) = MAX(0.0,cldwtr(i)) 503 betamx = 1.0 - c0*MAX(0.0,(dzcld(i)-dzmin)) 504 b1 = (hc(i) - hbs(i,km1))*dp(i,km1) 505 b2 = (hc(i) - hbs(i,k ))*dp(i,k ) 506 beta(i) = MAX(betamn,MIN(betamx, 1.0+b1/b2)) 507 IF (hbs(i,km1).le.hb(i,km1)) beta(i) = 0.0 508 ENDIF 509 ENDDO 510 C 511 C Bound maximum beta to ensure physically realistic solutions 512 C 513 C First check constrains beta so that eta remains positive 514 C (assuming that eta is already positive for beta equal zero) 515 c La premiere contrainte de beta est que le flux eta doit etre positif. 516 C 517 DO i=1,klon 518 IF (ldcum(i)) THEN 519 tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i)) 520 $ - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1) 521 tmp2 = (1.0+gam(i,k))*(sc(i)-sbh(i,k)) 522 IF ((beta(i)*tmp2-tmp1).GT.0.0) THEN 523 betamx = 0.99*(tmp1/tmp2) 524 beta(i) = MAX(0.0,MIN(betamx,beta(i))) 525 ENDIF 526 C 527 C Second check involves supersaturation of "detrainment layer" 528 C small amount of supersaturation acceptable (by ssfac factor) 529 c La 2e contrainte est que la convection ne doit pas sursaturer 530 c la "detrainment layer", Neanmoins, une petite sursaturation 531 c est acceptee (facteur ssfac). 532 C 533 IF (hb(i,km1).lt.hbs(i,km1)) THEN 534 tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i)) 535 $ - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1) 536 tmp1 = tmp1/dp(i,k) 537 tmp2 = gam(i,km1)*(sbh(i,k)-sc(i) + cldwtr(i)) - 538 $ hbh(i,k) + hc(i) - sc(i) + sbh(i,k) 539 tmp3 = (1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i,k) 540 tmp4 = (dt/cats)*(hc(i)-hbs(i,k))*tmp2 541 $ / (dp(i,km1)*(hbs(i,km1)-hb(i,km1))) + tmp3 542 IF ((beta(i)*tmp4-tmp1).GT.0.0) THEN 543 betamx = ssfac*(tmp1/tmp4) 544 beta(i) = MAX(0.0,MIN(betamx,beta(i))) 545 ENDIF 546 ELSE 547 beta(i) = 0.0 548 ENDIF 549 C 550 C Third check to avoid introducing 2 delta x thermodynamic 551 C noise in the vertical ... constrain adjusted h (or theta e) 552 C so that the adjustment doesn't contribute to "kinks" in h 553 C 554 g = MIN(0.0,hb(i,k)-hb(i,km1)) 555 tmp3 = (hb(i,k)-hb(i,km1)-g)*(cats/dt) / (hc(i)-hbs(i,k)) 556 tmp1 = (1.0+gam(i,k))*(sc(i)-sbh(i,kp1) + cldwtr(i)) 557 $ - (hbh(i,kp1)-hc(i))*dp(i,k)/dp(i,kp1) 558 tmp1 = tmp1/dp(i,k) 559 tmp1 = tmp3*tmp1 + (hc(i) - hbh(i,kp1))/dp(i,k) 560 tmp2 = tmp3*(1.0+gam(i,k))*(sc(i)-sbh(i,k))/dp(i,k) 561 $ + (hc(i)-hbh(i,k)-cldwtr(i)) 562 $ *(1.0/dp(i,k)+1.0/dp(i,kp1)) 563 IF ((beta(i)*tmp2-tmp1).GT.0.0) THEN 564 betamx = 0.0 565 IF (tmp2.NE.0.0) betamx = tmp1/tmp2 566 beta(i) = MAX(0.0,MIN(betamx,beta(i))) 567 ENDIF 568 ENDIF 569 ENDDO 570 C 571 C Calculate mass flux required for stabilization. 572 C 573 C Ensure that the convective mass flux, eta, is positive by 574 C setting negative values of eta to zero.. 575 C Ensure that estimated mass flux cannot move more than the 576 C minimum of total mass contained in either layer k or layer k+1. 577 C Also test for other pathological cases that result in non- 578 C physical states and adjust eta accordingly. 579 C 580 20 CONTINUE 581 DO i=1,klon 582 IF (ldcum(i)) THEN 583 beta(i) = MAX(0.0,beta(i)) 584 tmp1 = hc(i) - hbs(i,k) 585 tmp2 = ((1.0+gam(i,k))*(sc(i)-sbh(i,kp1)+cldwtr(i)) - 586 $ beta(i)*(1.0+gam(i,k))*(sc(i)-sbh(i,k)))/dp(i,k) - 587 $ (hbh(i,kp1)-hc(i))/dp(i,kp1) 588 eta(i) = tmp1/(tmp2*RG*cats) 589 tmass = MIN(dp(i,k),dp(i,kp1))/RG 590 IF (eta(i).GT.tmass*rdt .OR. eta(i).LE.0.0) eta(i) = 0.0 591 C 592 C Check on negative q in top layer (bound beta) 593 C 594 IF(shc(i)-shbh(i,k).LT.0.0 .AND. beta(i)*eta(i).NE.0.0)THEN 595 denom = eta(i)*RG*dt*(shc(i) - shbh(i,k))/dp(i,km1) 596 beta(i) = MAX(0.0,MIN(-0.999*shb(i,km1)/denom,beta(i))) 597 ENDIF 598 C 599 C Check on negative q in middle layer (zero eta) 600 C 601 qtest1 = shb(i,k) + eta(i)*RG*dt*((shc(i) - shbh(i,kp1)) - 602 $ (1.0 - beta(i))*cldwtr(i)/RLVTT - 603 $ beta(i)*(shc(i) - shbh(i,k)))/dp(i,k) 604 IF (qtest1.le.0.0) eta(i) = 0.0 605 C 606 C Check on negative q in lower layer (bound eta) 607 C 608 fac1 = -(shbh(i,kp1) - shc(i))/dp(i,kp1) 609 qtest2 = shb(i,kp1) - eta(i)*RG*dt*fac1 610 IF (qtest2 .lt. 0.0) THEN 611 eta(i) = 0.99*shb(i,kp1)/(RG*dt*fac1) 612 ENDIF 613 ENDIF 614 ENDDO 615 C 616 C 617 C Calculate cloud water, rain water, and thermodynamic changes 618 C 619 DO 30 i=1,klon 620 IF (ldcum(i)) THEN 621 etagdt = eta(i)*RG*dt 622 cldwtr(i) = etagdt*cldwtr(i)/RLVTT/RG 623 rnwtr(i) = (1.0 - beta(i))*cldwtr(i) 624 ds1(i) = etagdt*(sbh(i,kp1) - sc(i))/dp(i,kp1) 625 dq1(i) = etagdt*(shbh(i,kp1) - shc(i))/dp(i,kp1) 626 ds2(i) = (etagdt*(sc(i) - sbh(i,kp1)) + 627 $ RLVTT*RG*cldwtr(i) - beta(i)*etagdt* 628 $ (sc(i) - sbh(i,k)))/dp(i,k) 629 dq2(i) = (etagdt*(shc(i) - shbh(i,kp1)) - 630 $ RG*rnwtr(i) - beta(i)*etagdt* 631 $ (shc(i) - shbh(i,k)))/dp(i,k) 632 ds3(i) = beta(i)*(etagdt*(sc(i) - sbh(i,k)) - 633 $ RLVTT*RG*cldwtr(i))/dp(i,km1) 634 dq3(i) = beta(i)*etagdt*(shc(i) - shbh(i,k))/dp(i,km1) 635 C 636 C Isolate convective fluxes for later diagnostics 637 C 638 fslkp = eta(i)*(sc(i) - sbh(i,kp1)) 639 fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) - 640 $ RLVTT*cldwtr(i)*rdt) 641 fqlkp = eta(i)*(shc(i) - shbh(i,kp1)) 642 fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k)) 643 C 644 C 645 C Update thermodynamic profile (update sb, hb, & hbs later) 646 C 647 tb (i,kp1) = tb(i,kp1) + ds1(i) / RCPD 648 tb (i,k ) = tb(i,k ) + ds2(i) / RCPD 649 tb (i,km1) = tb(i,km1) + ds3(i) / RCPD 650 shb(i,kp1) = shb(i,kp1) + dq1(i) 651 shb(i,k ) = shb(i,k ) + dq2(i) 652 shb(i,km1) = shb(i,km1) + dq3(i) 653 prec(i) = prec(i) + rnwtr(i)/rhoh2o 654 C 655 C Update diagnostic information for final budget 656 C Tracking temperature & specific humidity tendencies, 657 C rainout term, convective mass flux, convective liquid 658 C water static energy flux, and convective total water flux 659 C 660 cmfdt (i,kp1) = cmfdt (i,kp1) + ds1(i)/RCPD*rdt 661 cmfdt (i,k ) = cmfdt (i,k ) + ds2(i)/RCPD*rdt 662 cmfdt (i,km1) = cmfdt (i,km1) + ds3(i)/RCPD*rdt 663 cmfdq (i,kp1) = cmfdq (i,kp1) + dq1(i)*rdt 664 cmfdq (i,k ) = cmfdq (i,k ) + dq2(i)*rdt 665 cmfdq (i,km1) = cmfdq (i,km1) + dq3(i)*rdt 666 cmfdqr(i,k ) = cmfdqr(i,k ) + (RG*rnwtr(i)/dp(i,k))*rdt 667 cmfmc (i,kp1) = cmfmc (i,kp1) + eta(i) 668 cmfmc (i,k ) = cmfmc (i,k ) + beta(i)*eta(i) 669 cmfsl (i,kp1) = cmfsl (i,kp1) + fslkp 670 cmfsl (i,k ) = cmfsl (i,k ) + fslkm 671 cmflq (i,kp1) = cmflq (i,kp1) + RLVTT*fqlkp 672 cmflq (i,k ) = cmflq (i,k ) + RLVTT*fqlkm 673 qc (i,k ) = (RG*rnwtr(i)/dp(i,k))*rdt 674 ENDIF 675 30 CONTINUE 676 C 677 C Next, convectively modify passive constituents 678 C 679 DO 50 m=1,pcnst 680 DO 40 i=1,klon 681 IF (ldcum(i)) THEN 682 C 683 C If any of the reported values of the constituent is negative in 684 C the three adjacent levels, nothing will be done to the profile 685 C 686 IF ((cmrb(i,kp1,m).LT.0.0) .OR. 687 $ (cmrb(i,k,m).LT.0.0) .OR. 688 $ (cmrb(i,km1,m).LT.0.0)) GOTO 40 689 C 690 C Specify constituent interface values (linear interpolation) 691 C 692 cmrh(i,k ) = 0.5*(cmrb(i,km1,m) + cmrb(i,k ,m)) 693 cmrh(i,kp1) = 0.5*(cmrb(i,k ,m) + cmrb(i,kp1,m)) 694 C 695 C Specify perturbation properties of constituents in PBL 696 C 697 pblhgt = MAX(pblh(i),1.0) 698 IF (gz(i,kp1)/RG.LE.pblhgt .AND. dzcld(i).EQ.0.) THEN 699 fac1 = MAX(0.0,1.0-gz(i,kp1)/RG/pblhgt) 700 cmrc(i) = cmrb(i,kp1,m) + cmrp(i,m)*fac1 701 ELSE 702 cmrc(i) = cmrb(i,kp1,m) 703 ENDIF 704 C 705 C Determine fluxes, flux divergence => changes due to convection 706 C Logic must be included to avoid producing negative values. A bit 707 C messy since there are no a priori assumptions about profiles. 708 C Tendency is modified (reduced) when pending disaster detected. 709 C 710 etagdt = eta(i)*RG*dt 711 botflx = etagdt*(cmrc(i) - cmrh(i,kp1)) 712 topflx = beta(i)*etagdt*(cmrc(i)-cmrh(i,k)) 713 dcmr1(i) = -botflx/dp(i,kp1) 714 efac1 = 1.0 715 efac2 = 1.0 716 efac3 = 1.0 717 C 718 IF (cmrb(i,kp1,m)+dcmr1(i) .LT. 0.0) THEN 719 efac1 = MAX(tiny,ABS(cmrb(i,kp1,m)/dcmr1(i)) - eps) 720 ENDIF 721 C 722 IF (efac1.EQ.tiny .OR. efac1.GT.1.0) efac1 = 0.0 723 dcmr1(i) = -efac1*botflx/dp(i,kp1) 724 dcmr2(i) = (efac1*botflx - topflx)/dp(i,k) 725 C 726 IF (cmrb(i,k,m)+dcmr2(i) .LT. 0.0) THEN 727 efac2 = MAX(tiny,ABS(cmrb(i,k ,m)/dcmr2(i)) - eps) 728 ENDIF 729 C 730 IF (efac2.EQ.tiny .OR. efac2.GT.1.0) efac2 = 0.0 731 dcmr2(i) = (efac1*botflx - efac2*topflx)/dp(i,k) 732 dcmr3(i) = efac2*topflx/dp(i,km1) 733 C 734 IF (cmrb(i,km1,m)+dcmr3(i) .LT. 0.0) THEN 735 efac3 = MAX(tiny,ABS(cmrb(i,km1,m)/dcmr3(i)) - eps) 736 ENDIF 737 C 738 IF (efac3.EQ.tiny .OR. efac3.GT.1.0) efac3 = 0.0 739 efac3 = MIN(efac2,efac3) 740 dcmr2(i) = (efac1*botflx - efac3*topflx)/dp(i,k) 741 dcmr3(i) = efac3*topflx/dp(i,km1) 742 C 743 cmrb(i,kp1,m) = cmrb(i,kp1,m) + dcmr1(i) 744 cmrb(i,k ,m) = cmrb(i,k ,m) + dcmr2(i) 745 cmrb(i,km1,m) = cmrb(i,km1,m) + dcmr3(i) 746 ENDIF 747 40 CONTINUE 748 50 CONTINUE ! end of m=1,pcnst loop 749 C 750 IF (k.EQ.limcnv+1) GOTO 60 ! on ne pourra plus glisser 751 c 752 c Dans la procedure de glissage ascendant, les variables thermo- 753 c dynamiques des couches k et km1 servent au calcul des couches 754 c superieures. Elles ont donc besoin d'une mise-a-jour. 755 C 756 DO i = 1, klon 757 IF (ldcum(i)) THEN 758 zx_t = tb(i,k) 759 zx_p = p(i,k) 760 zx_q = shb(i,k) 761 zdelta=MAX(0.,SIGN(1.,RTT-zx_t)) 762 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta 763 zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q) 764 zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p 765 zx_qs=MIN(0.5,zx_qs) 766 zcor=1./(1.-retv*zx_qs) 767 zx_qs=zx_qs*zcor 768 zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor) 769 shbs(i,k) = zx_qs 770 gam(i,k) = zx_gam 771 c 772 zx_t = tb(i,km1) 773 zx_p = p(i,km1) 774 zx_q = shb(i,km1) 775 zdelta=MAX(0.,SIGN(1.,RTT-zx_t)) 776 zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta 777 zcvm5 = zcvm5 / RCPD / (1.0+RVTMP2*zx_q) 778 zx_qs= r2es * FOEEW(zx_t,zdelta)/zx_p 779 zx_qs=MIN(0.5,zx_qs) 780 zcor=1./(1.-retv*zx_qs) 781 zx_qs=zx_qs*zcor 782 zx_gam = FOEDE(zx_t,zdelta,zcvm5,zx_qs,zcor) 783 shbs(i,km1) = zx_qs 784 gam(i,km1) = zx_gam 785 C 786 sb (i,k ) = sb(i,k ) + ds2(i) 787 sb (i,km1) = sb(i,km1) + ds3(i) 788 hb (i,k ) = sb(i,k ) + RLVTT*shb(i,k) 789 hb (i,km1) = sb(i,km1) + RLVTT*shb(i,km1) 790 hbs(i,k ) = sb(i,k ) + RLVTT*shbs(i,k ) 791 hbs(i,km1) = sb(i,km1) + RLVTT*shbs(i,km1) 792 C 793 sbh (i,k) = 0.5*(sb(i,k) + sb(i,km1)) 794 shbh(i,k) = qhalf(shb(i,km1),shb(i,k) 795 $ ,shbs(i,km1),shbs(i,k)) 796 hbh (i,k) = sbh(i,k) + RLVTT*shbh(i,k) 797 sbh (i,km1) = 0.5*(sb(i,km1) + sb(i,k-2)) 798 shbh(i,km1) = qhalf(shb(i,k-2),shb(i,km1), 799 $ shbs(i,k-2),shbs(i,km1)) 800 hbh (i,km1) = sbh(i,km1) + RLVTT*shbh(i,km1) 801 ENDIF 802 ENDDO 803 C 804 C Ensure that dzcld is reset if convective mass flux zero 805 C specify the current vertical extent of the convective activity 806 C top of convective layer determined by size of overshoot param. 807 C 808 60 CONTINUE 809 DO i=1,klon 810 etagt0 = eta(i).gt.0.0 811 IF (.not.etagt0) dzcld(i) = 0.0 812 IF (etagt0 .and. beta(i).gt.betamn) THEN 813 ktp = km1 814 ELSE 815 ktp = k 816 ENDIF 817 IF (etagt0) THEN 818 cnt(i) = MIN(cnt(i),ktp) 819 cnb(i) = MAX(cnb(i),k) 820 ENDIF 821 ENDDO 822 70 CONTINUE ! end of k loop 823 C 824 C determine whether precipitation, prec, is frozen (snow) or not 825 C 826 DO i=1,klon 827 IF (tb(i,klev).LT.tmelt .AND. tb(i,klev-1).lt.tmelt) THEN 828 cmfprs(i) = prec(i)*rdt 829 ELSE 830 cmfprt(i) = prec(i)*rdt 831 ENDIF 832 ENDDO 833 C 834 RETURN ! we're all done ... return to calling procedure 835 END 670 IF (ldcum(i)) THEN 671 672 ! If any of the reported values of the constituent is negative in 673 ! the three adjacent levels, nothing will be done to the profile 674 675 IF ((cmrb(i,kp1,m)<0.0) .OR. (cmrb(i,k,m)<0.0) .OR. (cmrb(i,km1, & 676 m)<0.0)) GO TO 40 677 678 ! Specify constituent interface values (linear interpolation) 679 680 cmrh(i, k) = 0.5*(cmrb(i,km1,m)+cmrb(i,k,m)) 681 cmrh(i, kp1) = 0.5*(cmrb(i,k,m)+cmrb(i,kp1,m)) 682 683 ! Specify perturbation properties of constituents in PBL 684 685 pblhgt = max(pblh(i), 1.0) 686 IF (gz(i,kp1)/rg<=pblhgt .AND. dzcld(i)==0.) THEN 687 fac1 = max(0.0, 1.0-gz(i,kp1)/rg/pblhgt) 688 cmrc(i) = cmrb(i, kp1, m) + cmrp(i, m)*fac1 689 ELSE 690 cmrc(i) = cmrb(i, kp1, m) 691 END IF 692 693 ! Determine fluxes, flux divergence => changes due to convection 694 ! Logic must be included to avoid producing negative values. A bit 695 ! messy since there are no a priori assumptions about profiles. 696 ! Tendency is modified (reduced) when pending disaster detected. 697 698 etagdt = eta(i)*rg*dt 699 botflx = etagdt*(cmrc(i)-cmrh(i,kp1)) 700 topflx = beta(i)*etagdt*(cmrc(i)-cmrh(i,k)) 701 dcmr1(i) = -botflx/dp(i, kp1) 702 efac1 = 1.0 703 efac2 = 1.0 704 efac3 = 1.0 705 706 IF (cmrb(i,kp1,m)+dcmr1(i)<0.0) THEN 707 efac1 = max(tiny, abs(cmrb(i,kp1,m)/dcmr1(i))-eps) 708 END IF 709 710 IF (efac1==tiny .OR. efac1>1.0) efac1 = 0.0 711 dcmr1(i) = -efac1*botflx/dp(i, kp1) 712 dcmr2(i) = (efac1*botflx-topflx)/dp(i, k) 713 714 IF (cmrb(i,k,m)+dcmr2(i)<0.0) THEN 715 efac2 = max(tiny, abs(cmrb(i,k,m)/dcmr2(i))-eps) 716 END IF 717 718 IF (efac2==tiny .OR. efac2>1.0) efac2 = 0.0 719 dcmr2(i) = (efac1*botflx-efac2*topflx)/dp(i, k) 720 dcmr3(i) = efac2*topflx/dp(i, km1) 721 722 IF (cmrb(i,km1,m)+dcmr3(i)<0.0) THEN 723 efac3 = max(tiny, abs(cmrb(i,km1,m)/dcmr3(i))-eps) 724 END IF 725 726 IF (efac3==tiny .OR. efac3>1.0) efac3 = 0.0 727 efac3 = min(efac2, efac3) 728 dcmr2(i) = (efac1*botflx-efac3*topflx)/dp(i, k) 729 dcmr3(i) = efac3*topflx/dp(i, km1) 730 731 cmrb(i, kp1, m) = cmrb(i, kp1, m) + dcmr1(i) 732 cmrb(i, k, m) = cmrb(i, k, m) + dcmr2(i) 733 cmrb(i, km1, m) = cmrb(i, km1, m) + dcmr3(i) 734 END IF 735 40 END DO 736 END DO ! end of m=1,pcnst loop 737 738 IF (k==limcnv+1) GO TO 60 ! on ne pourra plus glisser 739 740 ! Dans la procedure de glissage ascendant, les variables thermo- 741 ! dynamiques des couches k et km1 servent au calcul des couches 742 ! superieures. Elles ont donc besoin d'une mise-a-jour. 743 744 DO i = 1, klon 745 IF (ldcum(i)) THEN 746 zx_t = tb(i, k) 747 zx_p = p(i, k) 748 zx_q = shb(i, k) 749 zdelta = max(0., sign(1.,rtt-zx_t)) 750 zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta 751 zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zx_q) 752 zx_qs = r2es*foeew(zx_t, zdelta)/zx_p 753 zx_qs = min(0.5, zx_qs) 754 zcor = 1./(1.-retv*zx_qs) 755 zx_qs = zx_qs*zcor 756 zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) 757 shbs(i, k) = zx_qs 758 gam(i, k) = zx_gam 759 760 zx_t = tb(i, km1) 761 zx_p = p(i, km1) 762 zx_q = shb(i, km1) 763 zdelta = max(0., sign(1.,rtt-zx_t)) 764 zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta 765 zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zx_q) 766 zx_qs = r2es*foeew(zx_t, zdelta)/zx_p 767 zx_qs = min(0.5, zx_qs) 768 zcor = 1./(1.-retv*zx_qs) 769 zx_qs = zx_qs*zcor 770 zx_gam = foede(zx_t, zdelta, zcvm5, zx_qs, zcor) 771 shbs(i, km1) = zx_qs 772 gam(i, km1) = zx_gam 773 774 sb(i, k) = sb(i, k) + ds2(i) 775 sb(i, km1) = sb(i, km1) + ds3(i) 776 hb(i, k) = sb(i, k) + rlvtt*shb(i, k) 777 hb(i, km1) = sb(i, km1) + rlvtt*shb(i, km1) 778 hbs(i, k) = sb(i, k) + rlvtt*shbs(i, k) 779 hbs(i, km1) = sb(i, km1) + rlvtt*shbs(i, km1) 780 781 sbh(i, k) = 0.5*(sb(i,k)+sb(i,km1)) 782 shbh(i, k) = qhalf(shb(i,km1), shb(i,k), shbs(i,km1), shbs(i,k)) 783 hbh(i, k) = sbh(i, k) + rlvtt*shbh(i, k) 784 sbh(i, km1) = 0.5*(sb(i,km1)+sb(i,k-2)) 785 shbh(i, km1) = qhalf(shb(i,k-2), shb(i,km1), shbs(i,k-2), & 786 shbs(i,km1)) 787 hbh(i, km1) = sbh(i, km1) + rlvtt*shbh(i, km1) 788 END IF 789 END DO 790 791 ! Ensure that dzcld is reset if convective mass flux zero 792 ! specify the current vertical extent of the convective activity 793 ! top of convective layer determined by size of overshoot param. 794 795 60 CONTINUE 796 DO i = 1, klon 797 etagt0 = eta(i) > 0.0 798 IF (.NOT. etagt0) dzcld(i) = 0.0 799 IF (etagt0 .AND. beta(i)>betamn) THEN 800 ktp = km1 801 ELSE 802 ktp = k 803 END IF 804 IF (etagt0) THEN 805 cnt(i) = min(cnt(i), ktp) 806 cnb(i) = max(cnb(i), k) 807 END IF 808 END DO 809 70 END DO ! end of k loop 810 811 ! determine whether precipitation, prec, is frozen (snow) or not 812 813 DO i = 1, klon 814 IF (tb(i,klev)<tmelt .AND. tb(i,klev-1)<tmelt) THEN 815 cmfprs(i) = prec(i)*rdt 816 ELSE 817 cmfprt(i) = prec(i)*rdt 818 END IF 819 END DO 820 821 RETURN ! we're all done ... return to calling procedure 822 END SUBROUTINE cmfmca
Note: See TracChangeset
for help on using the changeset viewer.
