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