| 1 | ! |
|---|
| 2 | ! $Header$ |
|---|
| 3 | ! |
|---|
| 4 | subroutine convect1(len,nd,ndp1,noff,minorig, |
|---|
| 5 | & t,q,qs,u,v, |
|---|
| 6 | & p,ph,iflag,ft, |
|---|
| 7 | & fq,fu,fv,precip,cbmf,delt,Ma) |
|---|
| 8 | C.............................START PROLOGUE............................ |
|---|
| 9 | C |
|---|
| 10 | C SCCS IDENTIFICATION: @(#)convect1.f 1.1 04/21/00 |
|---|
| 11 | C 19:40:52 /h/cm/library/nogaps4/src/sub/fcst/convect1.f_v |
|---|
| 12 | C |
|---|
| 13 | C CONFIGURATION IDENTIFICATION: None |
|---|
| 14 | C |
|---|
| 15 | C MODULE NAME: convect1 |
|---|
| 16 | C |
|---|
| 17 | C DESCRIPTION: |
|---|
| 18 | C |
|---|
| 19 | C convect1 The Emanuel Cumulus Convection Scheme |
|---|
| 20 | C |
|---|
| 21 | C CONTRACT NUMBER AND TITLE: None |
|---|
| 22 | C |
|---|
| 23 | C REFERENCES: Programmers K. Emanuel (MIT), Timothy F. Hogan, M. Peng (NRL) |
|---|
| 24 | C |
|---|
| 25 | C CLASSIFICATION: Unclassified |
|---|
| 26 | C |
|---|
| 27 | C RESTRICTIONS: None |
|---|
| 28 | C |
|---|
| 29 | C COMPILER DEPENDENCIES: FORTRAN 77, FORTRAN 90 |
|---|
| 30 | C |
|---|
| 31 | C COMPILE OPTIONS: Fortran 77: -Zu -Wf"-ei -o aggress" |
|---|
| 32 | C Fortran 90: -O vector3,scalar3,task1,aggress,overindex -ei -r 2 |
|---|
| 33 | C |
|---|
| 34 | C LIBRARIES OF RESIDENCE: /a/ops/lib/libfcst159.a |
|---|
| 35 | C |
|---|
| 36 | C USAGE: call convect1(len,nd,noff,minorig, |
|---|
| 37 | C & t,q,qs,u,v, |
|---|
| 38 | C & p,ph,iflag,ft, |
|---|
| 39 | C & fq,fu,fv,precip,cbmf,delt) |
|---|
| 40 | C |
|---|
| 41 | C PARAMETERS: |
|---|
| 42 | C Name Type Usage Description |
|---|
| 43 | C ---------- ---------- ------- ---------------------------- |
|---|
| 44 | C |
|---|
| 45 | C len Integer Input first (i) dimension |
|---|
| 46 | C nd Integer Input vertical (k) dimension |
|---|
| 47 | C ndp1 Integer Input nd + 1 |
|---|
| 48 | C noff Integer Input integer limit for convection (nd-noff) |
|---|
| 49 | C minorig Integer Input First level of convection |
|---|
| 50 | C t Real Input temperature |
|---|
| 51 | C q Real Input specific hum |
|---|
| 52 | C qs Real Input sat specific hum |
|---|
| 53 | C u Real Input u-wind |
|---|
| 54 | C v Real Input v-wind |
|---|
| 55 | C p Real Input full level pressure |
|---|
| 56 | C ph Real Input half level pressure |
|---|
| 57 | C iflag Integer Output iflag on latitude strip |
|---|
| 58 | C ft Real Output temp tend |
|---|
| 59 | C fq Real Output spec hum tend |
|---|
| 60 | C fu Real Output u-wind tend |
|---|
| 61 | C fv Real Output v-wind tend |
|---|
| 62 | C cbmf Real In/Out cumulus mass flux |
|---|
| 63 | C delt Real Input time step |
|---|
| 64 | C iflag Integer Output integer flag for Emanuel conditions |
|---|
| 65 | C |
|---|
| 66 | C COMMON BLOCKS: |
|---|
| 67 | C Block Name Type Usage Notes |
|---|
| 68 | C -------- -------- ---- ------ ------------------------ |
|---|
| 69 | C |
|---|
| 70 | C FILES: None |
|---|
| 71 | C |
|---|
| 72 | C DATA BASES: None |
|---|
| 73 | C |
|---|
| 74 | C NON-FILE INPUT/OUTPUT: None |
|---|
| 75 | C |
|---|
| 76 | C ERROR CONDITIONS: None |
|---|
| 77 | C |
|---|
| 78 | C ADDITIONAL COMMENTS: None |
|---|
| 79 | C |
|---|
| 80 | C.................MAINTENANCE SECTION................................ |
|---|
| 81 | C |
|---|
| 82 | C MODULES CALLED: |
|---|
| 83 | C Name Description |
|---|
| 84 | C convect2 Emanuel cumulus convection tendency calculations |
|---|
| 85 | C ------- ---------------------- |
|---|
| 86 | C LOCAL VARIABLES AND |
|---|
| 87 | C STRUCTURES: |
|---|
| 88 | C Name Type Description |
|---|
| 89 | C ------- ------ ----------- |
|---|
| 90 | C See Comments Below |
|---|
| 91 | C |
|---|
| 92 | C i Integer loop index |
|---|
| 93 | C k Integer loop index |
|---|
| 94 | c |
|---|
| 95 | C METHOD: |
|---|
| 96 | C |
|---|
| 97 | C See Emanuel, K. and M. Zivkovic-Rothman, 2000: Development and evaluation of a |
|---|
| 98 | C convective scheme for use in climate models. |
|---|
| 99 | C |
|---|
| 100 | C FILES: None |
|---|
| 101 | C |
|---|
| 102 | C INCLUDE FILES: None |
|---|
| 103 | C |
|---|
| 104 | C MAKEFILE: /a/ops/met/nogaps/src/sub/fcst/fcst159lib.mak |
|---|
| 105 | C |
|---|
| 106 | C..............................END PROLOGUE............................. |
|---|
| 107 | c |
|---|
| 108 | c |
|---|
| 109 | implicit none |
|---|
| 110 | c |
|---|
| 111 | #include "dimensions.h" |
|---|
| 112 | #include "dimphy.h" |
|---|
| 113 | c |
|---|
| 114 | integer len |
|---|
| 115 | integer nd |
|---|
| 116 | integer ndp1 |
|---|
| 117 | integer noff |
|---|
| 118 | real t(len,nd) |
|---|
| 119 | real q(len,nd) |
|---|
| 120 | real qs(len,nd) |
|---|
| 121 | real u(len,nd) |
|---|
| 122 | real v(len,nd) |
|---|
| 123 | real p(len,nd) |
|---|
| 124 | real ph(len,ndp1) |
|---|
| 125 | integer iflag(len) |
|---|
| 126 | real ft(len,nd) |
|---|
| 127 | real fq(len,nd) |
|---|
| 128 | real fu(len,nd) |
|---|
| 129 | real fv(len,nd) |
|---|
| 130 | real precip(len) |
|---|
| 131 | real cbmf(len) |
|---|
| 132 | real Ma(len,nd) |
|---|
| 133 | integer minorig |
|---|
| 134 | real delt,cpd,cpv,cl,rv,rd,lv0,g |
|---|
| 135 | real sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp |
|---|
| 136 | real alpha,entp,coeffs,coeffr,omtrain,cu |
|---|
| 137 | c |
|---|
| 138 | !------------------------------------------------------------------- |
|---|
| 139 | ! --- ARGUMENTS |
|---|
| 140 | !------------------------------------------------------------------- |
|---|
| 141 | ! --- On input: |
|---|
| 142 | ! |
|---|
| 143 | ! t: Array of absolute temperature (K) of dimension ND, with first |
|---|
| 144 | ! index corresponding to lowest model level. Note that this array |
|---|
| 145 | ! will be altered by the subroutine if dry convective adjustment |
|---|
| 146 | ! occurs and if IPBL is not equal to 0. |
|---|
| 147 | ! |
|---|
| 148 | ! q: Array of specific humidity (gm/gm) of dimension ND, with first |
|---|
| 149 | ! index corresponding to lowest model level. Must be defined |
|---|
| 150 | ! at same grid levels as T. Note that this array will be altered |
|---|
| 151 | ! if dry convective adjustment occurs and if IPBL is not equal to 0. |
|---|
| 152 | ! |
|---|
| 153 | ! qs: Array of saturation specific humidity of dimension ND, with first |
|---|
| 154 | ! index corresponding to lowest model level. Must be defined |
|---|
| 155 | ! at same grid levels as T. Note that this array will be altered |
|---|
| 156 | ! if dry convective adjustment occurs and if IPBL is not equal to 0. |
|---|
| 157 | ! |
|---|
| 158 | ! u: Array of zonal wind velocity (m/s) of dimension ND, witth first |
|---|
| 159 | ! index corresponding with the lowest model level. Defined at |
|---|
| 160 | ! same levels as T. Note that this array will be altered if |
|---|
| 161 | ! dry convective adjustment occurs and if IPBL is not equal to 0. |
|---|
| 162 | ! |
|---|
| 163 | ! v: Same as u but for meridional velocity. |
|---|
| 164 | ! |
|---|
| 165 | ! tra: Array of passive tracer mixing ratio, of dimensions (ND,NTRA), |
|---|
| 166 | ! where NTRA is the number of different tracers. If no |
|---|
| 167 | ! convective tracer transport is needed, define a dummy |
|---|
| 168 | ! input array of dimension (ND,1). Tracers are defined at |
|---|
| 169 | ! same vertical levels as T. Note that this array will be altered |
|---|
| 170 | ! if dry convective adjustment occurs and if IPBL is not equal to 0. |
|---|
| 171 | ! |
|---|
| 172 | ! p: Array of pressure (mb) of dimension ND, with first |
|---|
| 173 | ! index corresponding to lowest model level. Must be defined |
|---|
| 174 | ! at same grid levels as T. |
|---|
| 175 | ! |
|---|
| 176 | ! ph: Array of pressure (mb) of dimension ND+1, with first index |
|---|
| 177 | ! corresponding to lowest level. These pressures are defined at |
|---|
| 178 | ! levels intermediate between those of P, T, Q and QS. The first |
|---|
| 179 | ! value of PH should be greater than (i.e. at a lower level than) |
|---|
| 180 | ! the first value of the array P. |
|---|
| 181 | ! |
|---|
| 182 | ! nl: The maximum number of levels to which convection can penetrate, plus 1. |
|---|
| 183 | ! NL MUST be less than or equal to ND-1. |
|---|
| 184 | ! |
|---|
| 185 | ! delt: The model time step (sec) between calls to CONVECT |
|---|
| 186 | ! |
|---|
| 187 | !---------------------------------------------------------------------------- |
|---|
| 188 | ! --- On Output: |
|---|
| 189 | ! |
|---|
| 190 | ! iflag: An output integer whose value denotes the following: |
|---|
| 191 | ! VALUE INTERPRETATION |
|---|
| 192 | ! ----- -------------- |
|---|
| 193 | ! 0 Moist convection occurs. |
|---|
| 194 | ! 1 Moist convection occurs, but a CFL condition |
|---|
| 195 | ! on the subsidence warming is violated. This |
|---|
| 196 | ! does not cause the scheme to terminate. |
|---|
| 197 | ! 2 Moist convection, but no precip because ep(inb) lt 0.0001 |
|---|
| 198 | ! 3 No moist convection because new cbmf is 0 and old cbmf is 0. |
|---|
| 199 | ! 4 No moist convection; atmosphere is not |
|---|
| 200 | ! unstable |
|---|
| 201 | ! 6 No moist convection because ihmin le minorig. |
|---|
| 202 | ! 7 No moist convection because unreasonable |
|---|
| 203 | ! parcel level temperature or specific humidity. |
|---|
| 204 | ! 8 No moist convection: lifted condensation |
|---|
| 205 | ! level is above the 200 mb level. |
|---|
| 206 | ! 9 No moist convection: cloud base is higher |
|---|
| 207 | ! then the level NL-1. |
|---|
| 208 | ! |
|---|
| 209 | ! ft: Array of temperature tendency (K/s) of dimension ND, defined at same |
|---|
| 210 | ! grid levels as T, Q, QS and P. |
|---|
| 211 | ! |
|---|
| 212 | ! fq: Array of specific humidity tendencies ((gm/gm)/s) of dimension ND, |
|---|
| 213 | ! defined at same grid levels as T, Q, QS and P. |
|---|
| 214 | ! |
|---|
| 215 | ! fu: Array of forcing of zonal velocity (m/s^2) of dimension ND, |
|---|
| 216 | ! defined at same grid levels as T. |
|---|
| 217 | ! |
|---|
| 218 | ! fv: Same as FU, but for forcing of meridional velocity. |
|---|
| 219 | ! |
|---|
| 220 | ! ftra: Array of forcing of tracer content, in tracer mixing ratio per |
|---|
| 221 | ! second, defined at same levels as T. Dimensioned (ND,NTRA). |
|---|
| 222 | ! |
|---|
| 223 | ! precip: Scalar convective precipitation rate (mm/day). |
|---|
| 224 | ! |
|---|
| 225 | ! wd: A convective downdraft velocity scale. For use in surface |
|---|
| 226 | ! flux parameterizations. See convect.ps file for details. |
|---|
| 227 | ! |
|---|
| 228 | ! tprime: A convective downdraft temperature perturbation scale (K). |
|---|
| 229 | ! For use in surface flux parameterizations. See convect.ps |
|---|
| 230 | ! file for details. |
|---|
| 231 | ! |
|---|
| 232 | ! qprime: A convective downdraft specific humidity |
|---|
| 233 | ! perturbation scale (gm/gm). |
|---|
| 234 | ! For use in surface flux parameterizations. See convect.ps |
|---|
| 235 | ! file for details. |
|---|
| 236 | ! |
|---|
| 237 | ! cbmf: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST |
|---|
| 238 | ! BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT |
|---|
| 239 | ! ITS NEXT CALL. That is, the value of CBMF must be "remembered" |
|---|
| 240 | ! by the calling program between calls to CONVECT. |
|---|
| 241 | ! |
|---|
| 242 | ! det: Array of detrainment mass flux of dimension ND. |
|---|
| 243 | ! |
|---|
| 244 | !------------------------------------------------------------------- |
|---|
| 245 | c |
|---|
| 246 | c Local arrays |
|---|
| 247 | c |
|---|
| 248 | integer nl |
|---|
| 249 | integer nlp |
|---|
| 250 | integer nlm |
|---|
| 251 | integer i,k,n |
|---|
| 252 | real delti |
|---|
| 253 | real rowl |
|---|
| 254 | real clmcpv |
|---|
| 255 | real clmcpd |
|---|
| 256 | real cpdmcp |
|---|
| 257 | real cpvmcpd |
|---|
| 258 | real eps |
|---|
| 259 | real epsi |
|---|
| 260 | real epsim1 |
|---|
| 261 | real ginv |
|---|
| 262 | real hrd |
|---|
| 263 | real prccon1 |
|---|
| 264 | integer icbmax |
|---|
| 265 | real lv(klon,klev) |
|---|
| 266 | real cpn(klon,klev) |
|---|
| 267 | real cpx(klon,klev) |
|---|
| 268 | real tv(klon,klev) |
|---|
| 269 | real gz(klon,klev) |
|---|
| 270 | real hm(klon,klev) |
|---|
| 271 | real h(klon,klev) |
|---|
| 272 | real work(klon) |
|---|
| 273 | integer ihmin(klon) |
|---|
| 274 | integer nk(klon) |
|---|
| 275 | real rh(klon) |
|---|
| 276 | real chi(klon) |
|---|
| 277 | real plcl(klon) |
|---|
| 278 | integer icb(klon) |
|---|
| 279 | real tnk(klon) |
|---|
| 280 | real qnk(klon) |
|---|
| 281 | real gznk(klon) |
|---|
| 282 | real pnk(klon) |
|---|
| 283 | real qsnk(klon) |
|---|
| 284 | real ticb(klon) |
|---|
| 285 | real gzicb(klon) |
|---|
| 286 | real tp(klon,klev) |
|---|
| 287 | real tvp(klon,klev) |
|---|
| 288 | real clw(klon,klev) |
|---|
| 289 | c |
|---|
| 290 | real ah0(klon),cpp(klon) |
|---|
| 291 | real tg,qg,s,alv,tc,ahg,denom,es,rg |
|---|
| 292 | c |
|---|
| 293 | integer ncum |
|---|
| 294 | integer idcum(klon) |
|---|
| 295 | c |
|---|
| 296 | cpd=1005.7 |
|---|
| 297 | cpv=1870.0 |
|---|
| 298 | cl=4190.0 |
|---|
| 299 | rv=461.5 |
|---|
| 300 | rd=287.04 |
|---|
| 301 | lv0=2.501E6 |
|---|
| 302 | g=9.8 |
|---|
| 303 | C |
|---|
| 304 | C *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *** |
|---|
| 305 | C *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *** |
|---|
| 306 | C *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *** |
|---|
| 307 | C *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *** |
|---|
| 308 | C *** BETWEEN 0 C AND TLCRIT) *** |
|---|
| 309 | C *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *** |
|---|
| 310 | C *** FORMULATION *** |
|---|
| 311 | C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** |
|---|
| 312 | C *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** |
|---|
| 313 | C *** OF CLOUD *** |
|---|
| 314 | C *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *** |
|---|
| 315 | C *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *** |
|---|
| 316 | C *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** |
|---|
| 317 | C *** OF RAIN *** |
|---|
| 318 | C *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** |
|---|
| 319 | C *** OF SNOW *** |
|---|
| 320 | C *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *** |
|---|
| 321 | C *** TRANSPORT *** |
|---|
| 322 | C *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *** |
|---|
| 323 | C *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *** |
|---|
| 324 | C *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *** |
|---|
| 325 | C *** APPROACH TO QUASI-EQUILIBRIUM *** |
|---|
| 326 | C *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *** |
|---|
| 327 | C *** (DAMP MUST BE LESS THAN 1) *** |
|---|
| 328 | c |
|---|
| 329 | sigs=0.12 |
|---|
| 330 | sigd=0.05 |
|---|
| 331 | elcrit=0.0011 |
|---|
| 332 | tlcrit=-55.0 |
|---|
| 333 | omtsnow=5.5 |
|---|
| 334 | dtmax=0.9 |
|---|
| 335 | damp=0.1 |
|---|
| 336 | alpha=0.2 |
|---|
| 337 | entp=1.5 |
|---|
| 338 | coeffs=0.8 |
|---|
| 339 | coeffr=1.0 |
|---|
| 340 | omtrain=50.0 |
|---|
| 341 | c |
|---|
| 342 | cu=0.70 |
|---|
| 343 | damp=0.1 |
|---|
| 344 | c |
|---|
| 345 | c |
|---|
| 346 | c Define nl, nlp, nlm, and delti |
|---|
| 347 | c |
|---|
| 348 | nl=nd-noff |
|---|
| 349 | nlp=nl+1 |
|---|
| 350 | nlm=nl-1 |
|---|
| 351 | delti=1.0/delt |
|---|
| 352 | ! |
|---|
| 353 | !------------------------------------------------------------------- |
|---|
| 354 | ! --- SET CONSTANTS |
|---|
| 355 | !------------------------------------------------------------------- |
|---|
| 356 | ! |
|---|
| 357 | rowl=1000.0 |
|---|
| 358 | clmcpv=cl-cpv |
|---|
| 359 | clmcpd=cl-cpd |
|---|
| 360 | cpdmcp=cpd-cpv |
|---|
| 361 | cpvmcpd=cpv-cpd |
|---|
| 362 | eps=rd/rv |
|---|
| 363 | epsi=1.0/eps |
|---|
| 364 | epsim1=epsi-1.0 |
|---|
| 365 | ginv=1.0/g |
|---|
| 366 | hrd=0.5*rd |
|---|
| 367 | prccon1=86400.0*1000.0/(rowl*g) |
|---|
| 368 | ! |
|---|
| 369 | ! dtmax is the maximum negative temperature perturbation. |
|---|
| 370 | ! |
|---|
| 371 | !===================================================================== |
|---|
| 372 | ! --- INITIALIZE OUTPUT ARRAYS AND PARAMETERS |
|---|
| 373 | !===================================================================== |
|---|
| 374 | ! |
|---|
| 375 | do 20 k=1,nd |
|---|
| 376 | do 10 i=1,len |
|---|
| 377 | ft(i,k)=0.0 |
|---|
| 378 | fq(i,k)=0.0 |
|---|
| 379 | fu(i,k)=0.0 |
|---|
| 380 | fv(i,k)=0.0 |
|---|
| 381 | tvp(i,k)=0.0 |
|---|
| 382 | tp(i,k)=0.0 |
|---|
| 383 | clw(i,k)=0.0 |
|---|
| 384 | gz(i,k) = 0. |
|---|
| 385 | 10 continue |
|---|
| 386 | 20 continue |
|---|
| 387 | do 60 i=1,len |
|---|
| 388 | precip(i)=0.0 |
|---|
| 389 | iflag(i)=0 |
|---|
| 390 | 60 continue |
|---|
| 391 | c |
|---|
| 392 | !===================================================================== |
|---|
| 393 | ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY |
|---|
| 394 | !===================================================================== |
|---|
| 395 | do 110 k=1,nl+1 |
|---|
| 396 | do 100 i=1,len |
|---|
| 397 | lv(i,k)= lv0-clmcpv*(t(i,k)-273.15) |
|---|
| 398 | cpn(i,k)=cpd*(1.0-q(i,k))+cpv*q(i,k) |
|---|
| 399 | cpx(i,k)=cpd*(1.0-q(i,k))+cl*q(i,k) |
|---|
| 400 | tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1) |
|---|
| 401 | 100 continue |
|---|
| 402 | 110 continue |
|---|
| 403 | c |
|---|
| 404 | c gz = phi at the full levels (same as p). |
|---|
| 405 | c |
|---|
| 406 | do 120 i=1,len |
|---|
| 407 | gz(i,1)=0.0 |
|---|
| 408 | 120 continue |
|---|
| 409 | do 140 k=2,nlp |
|---|
| 410 | do 130 i=1,len |
|---|
| 411 | gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k)) |
|---|
| 412 | & *(p(i,k-1)-p(i,k))/ph(i,k) |
|---|
| 413 | 130 continue |
|---|
| 414 | 140 continue |
|---|
| 415 | c |
|---|
| 416 | c h = phi + cpT (dry static energy). |
|---|
| 417 | c hm = phi + cp(T-Tbase)+Lq |
|---|
| 418 | c |
|---|
| 419 | do 170 k=1,nlp |
|---|
| 420 | do 160 i=1,len |
|---|
| 421 | h(i,k)=gz(i,k)+cpn(i,k)*t(i,k) |
|---|
| 422 | hm(i,k)=gz(i,k)+cpx(i,k)*(t(i,k)-t(i,1))+lv(i,k)*q(i,k) |
|---|
| 423 | 160 continue |
|---|
| 424 | 170 continue |
|---|
| 425 | c |
|---|
| 426 | !------------------------------------------------------------------- |
|---|
| 427 | ! --- Find level of minimum moist static energy |
|---|
| 428 | ! --- If level of minimum moist static energy coincides with |
|---|
| 429 | ! --- or is lower than minimum allowable parcel origin level, |
|---|
| 430 | ! --- set iflag to 6. |
|---|
| 431 | !------------------------------------------------------------------- |
|---|
| 432 | do 180 i=1,len |
|---|
| 433 | work(i)=1.0e12 |
|---|
| 434 | ihmin(i)=nl |
|---|
| 435 | 180 continue |
|---|
| 436 | do 200 k=2,nlp |
|---|
| 437 | do 190 i=1,len |
|---|
| 438 | if((hm(i,k).lt.work(i)).and. |
|---|
| 439 | & (hm(i,k).lt.hm(i,k-1)))then |
|---|
| 440 | work(i)=hm(i,k) |
|---|
| 441 | ihmin(i)=k |
|---|
| 442 | endif |
|---|
| 443 | 190 continue |
|---|
| 444 | 200 continue |
|---|
| 445 | do 210 i=1,len |
|---|
| 446 | ihmin(i)=min(ihmin(i),nlm) |
|---|
| 447 | if(ihmin(i).le.minorig)then |
|---|
| 448 | iflag(i)=6 |
|---|
| 449 | endif |
|---|
| 450 | 210 continue |
|---|
| 451 | c |
|---|
| 452 | !------------------------------------------------------------------- |
|---|
| 453 | ! --- Find that model level below the level of minimum moist static |
|---|
| 454 | ! --- energy that has the maximum value of moist static energy |
|---|
| 455 | !------------------------------------------------------------------- |
|---|
| 456 | |
|---|
| 457 | do 220 i=1,len |
|---|
| 458 | work(i)=hm(i,minorig) |
|---|
| 459 | nk(i)=minorig |
|---|
| 460 | 220 continue |
|---|
| 461 | do 240 k=minorig+1,nl |
|---|
| 462 | do 230 i=1,len |
|---|
| 463 | if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then |
|---|
| 464 | work(i)=hm(i,k) |
|---|
| 465 | nk(i)=k |
|---|
| 466 | endif |
|---|
| 467 | 230 continue |
|---|
| 468 | 240 continue |
|---|
| 469 | !------------------------------------------------------------------- |
|---|
| 470 | ! --- Check whether parcel level temperature and specific humidity |
|---|
| 471 | ! --- are reasonable |
|---|
| 472 | !------------------------------------------------------------------- |
|---|
| 473 | do 250 i=1,len |
|---|
| 474 | if(((t(i,nk(i)).lt.250.0).or. |
|---|
| 475 | & (q(i,nk(i)).le.0.0).or. |
|---|
| 476 | & (p(i,ihmin(i)).lt.400.0)).and. |
|---|
| 477 | & (iflag(i).eq.0))iflag(i)=7 |
|---|
| 478 | 250 continue |
|---|
| 479 | !------------------------------------------------------------------- |
|---|
| 480 | ! --- Calculate lifted condensation level of air at parcel origin level |
|---|
| 481 | ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980) |
|---|
| 482 | !------------------------------------------------------------------- |
|---|
| 483 | do 260 i=1,len |
|---|
| 484 | tnk(i)=t(i,nk(i)) |
|---|
| 485 | qnk(i)=q(i,nk(i)) |
|---|
| 486 | gznk(i)=gz(i,nk(i)) |
|---|
| 487 | pnk(i)=p(i,nk(i)) |
|---|
| 488 | qsnk(i)=qs(i,nk(i)) |
|---|
| 489 | c |
|---|
| 490 | rh(i)=qnk(i)/qsnk(i) |
|---|
| 491 | rh(i)=min(1.0,rh(i)) |
|---|
| 492 | chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i)) |
|---|
| 493 | plcl(i)=pnk(i)*(rh(i)**chi(i)) |
|---|
| 494 | if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0)) |
|---|
| 495 | & .and.(iflag(i).eq.0))iflag(i)=8 |
|---|
| 496 | 260 continue |
|---|
| 497 | !------------------------------------------------------------------- |
|---|
| 498 | ! --- Calculate first level above lcl (=icb) |
|---|
| 499 | !------------------------------------------------------------------- |
|---|
| 500 | do 270 i=1,len |
|---|
| 501 | icb(i)=nlm |
|---|
| 502 | 270 continue |
|---|
| 503 | c |
|---|
| 504 | do 290 k=minorig,nl |
|---|
| 505 | do 280 i=1,len |
|---|
| 506 | if((k.ge.(nk(i)+1)).and.(p(i,k).lt.plcl(i))) |
|---|
| 507 | & icb(i)=min(icb(i),k) |
|---|
| 508 | 280 continue |
|---|
| 509 | 290 continue |
|---|
| 510 | c |
|---|
| 511 | do 300 i=1,len |
|---|
| 512 | if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9 |
|---|
| 513 | 300 continue |
|---|
| 514 | c |
|---|
| 515 | c Compute icbmax. |
|---|
| 516 | c |
|---|
| 517 | icbmax=2 |
|---|
| 518 | do 310 i=1,len |
|---|
| 519 | icbmax=max(icbmax,icb(i)) |
|---|
| 520 | 310 continue |
|---|
| 521 | ! |
|---|
| 522 | !------------------------------------------------------------------- |
|---|
| 523 | ! --- Calculates the lifted parcel virtual temperature at nk, |
|---|
| 524 | ! --- the actual temperature, and the adiabatic |
|---|
| 525 | ! --- liquid water content. The procedure is to solve the equation. |
|---|
| 526 | ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. |
|---|
| 527 | !------------------------------------------------------------------- |
|---|
| 528 | ! |
|---|
| 529 | do 320 i=1,len |
|---|
| 530 | tnk(i)=t(i,nk(i)) |
|---|
| 531 | qnk(i)=q(i,nk(i)) |
|---|
| 532 | gznk(i)=gz(i,nk(i)) |
|---|
| 533 | ticb(i)=t(i,icb(i)) |
|---|
| 534 | gzicb(i)=gz(i,icb(i)) |
|---|
| 535 | 320 continue |
|---|
| 536 | c |
|---|
| 537 | c *** Calculate certain parcel quantities, including static energy *** |
|---|
| 538 | c |
|---|
| 539 | do 330 i=1,len |
|---|
| 540 | ah0(i)=(cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) |
|---|
| 541 | & +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15))+gznk(i) |
|---|
| 542 | cpp(i)=cpd*(1.-qnk(i))+qnk(i)*cpv |
|---|
| 543 | 330 continue |
|---|
| 544 | c |
|---|
| 545 | c *** Calculate lifted parcel quantities below cloud base *** |
|---|
| 546 | c |
|---|
| 547 | do 350 k=minorig,icbmax-1 |
|---|
| 548 | do 340 i=1,len |
|---|
| 549 | tp(i,k)=tnk(i)-(gz(i,k)-gznk(i))/cpp(i) |
|---|
| 550 | tvp(i,k)=tp(i,k)*(1.+qnk(i)*epsi) |
|---|
| 551 | 340 continue |
|---|
| 552 | 350 continue |
|---|
| 553 | c |
|---|
| 554 | c *** Find lifted parcel quantities above cloud base *** |
|---|
| 555 | c |
|---|
| 556 | do 360 i=1,len |
|---|
| 557 | tg=ticb(i) |
|---|
| 558 | qg=qs(i,icb(i)) |
|---|
| 559 | alv=lv0-clmcpv*(ticb(i)-273.15) |
|---|
| 560 | c |
|---|
| 561 | c First iteration. |
|---|
| 562 | c |
|---|
| 563 | s=cpd+alv*alv*qg/(rv*ticb(i)*ticb(i)) |
|---|
| 564 | s=1./s |
|---|
| 565 | ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) |
|---|
| 566 | tg=tg+s*(ah0(i)-ahg) |
|---|
| 567 | tg=max(tg,35.0) |
|---|
| 568 | tc=tg-273.15 |
|---|
| 569 | denom=243.5+tc |
|---|
| 570 | if(tc.ge.0.0)then |
|---|
| 571 | es=6.112*exp(17.67*tc/denom) |
|---|
| 572 | else |
|---|
| 573 | es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
|---|
| 574 | endif |
|---|
| 575 | qg=eps*es/(p(i,icb(i))-es*(1.-eps)) |
|---|
| 576 | c |
|---|
| 577 | c Second iteration. |
|---|
| 578 | c |
|---|
| 579 | s=cpd+alv*alv*qg/(rv*ticb(i)*ticb(i)) |
|---|
| 580 | s=1./s |
|---|
| 581 | ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) |
|---|
| 582 | tg=tg+s*(ah0(i)-ahg) |
|---|
| 583 | tg=max(tg,35.0) |
|---|
| 584 | tc=tg-273.15 |
|---|
| 585 | denom=243.5+tc |
|---|
| 586 | if(tc.ge.0.0)then |
|---|
| 587 | es=6.112*exp(17.67*tc/denom) |
|---|
| 588 | else |
|---|
| 589 | es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) |
|---|
| 590 | end if |
|---|
| 591 | qg=eps*es/(p(i,icb(i))-es*(1.-eps)) |
|---|
| 592 | c |
|---|
| 593 | alv=lv0-clmcpv*(ticb(i)-273.15) |
|---|
| 594 | tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) |
|---|
| 595 | & -gz(i,icb(i))-alv*qg)/cpd |
|---|
| 596 | clw(i,icb(i))=qnk(i)-qg |
|---|
| 597 | clw(i,icb(i))=max(0.0,clw(i,icb(i))) |
|---|
| 598 | rg=qg/(1.-qnk(i)) |
|---|
| 599 | tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi) |
|---|
| 600 | 360 continue |
|---|
| 601 | c |
|---|
| 602 | do 380 k=minorig,icbmax |
|---|
| 603 | do 370 i=1,len |
|---|
| 604 | tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i) |
|---|
| 605 | 370 continue |
|---|
| 606 | 380 continue |
|---|
| 607 | c |
|---|
| 608 | !------------------------------------------------------------------- |
|---|
| 609 | ! --- Test for instability. |
|---|
| 610 | ! --- If there was no convection at last time step and parcel |
|---|
| 611 | ! --- is stable at icb, then set iflag to 4. |
|---|
| 612 | !------------------------------------------------------------------- |
|---|
| 613 | |
|---|
| 614 | do 390 i=1,len |
|---|
| 615 | if((cbmf(i).eq.0.0) .and.(iflag(i).eq.0).and. |
|---|
| 616 | & (tvp(i,icb(i)).le.(tv(i,icb(i))-dtmax)))iflag(i)=4 |
|---|
| 617 | 390 continue |
|---|
| 618 | |
|---|
| 619 | !===================================================================== |
|---|
| 620 | ! --- IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY |
|---|
| 621 | !===================================================================== |
|---|
| 622 | c |
|---|
| 623 | ncum=0 |
|---|
| 624 | do 400 i=1,len |
|---|
| 625 | if(iflag(i).eq.0)then |
|---|
| 626 | ncum=ncum+1 |
|---|
| 627 | idcum(ncum)=i |
|---|
| 628 | endif |
|---|
| 629 | 400 continue |
|---|
| 630 | c |
|---|
| 631 | c Call convect2, which compresses the points and computes the heating, |
|---|
| 632 | c moistening, velocity mixing, and precipiation. |
|---|
| 633 | c |
|---|
| 634 | c print*,'cpd avant convect2 ',cpd |
|---|
| 635 | if(ncum.gt.0)then |
|---|
| 636 | call convect2(ncum,idcum,len,nd,ndp1,nl,minorig, |
|---|
| 637 | & nk,icb, |
|---|
| 638 | & t,q,qs,u,v,gz,tv,tp,tvp,clw,h, |
|---|
| 639 | & lv,cpn,p,ph,ft,fq,fu,fv, |
|---|
| 640 | & tnk,qnk,gznk,plcl, |
|---|
| 641 | & precip,cbmf,iflag, |
|---|
| 642 | & delt,cpd,cpv,cl,rv,rd,lv0,g, |
|---|
| 643 | & sigs,sigd,elcrit,tlcrit,omtsnow,dtmax,damp, |
|---|
| 644 | & alpha,entp,coeffs,coeffr,omtrain,cu,Ma) |
|---|
| 645 | endif |
|---|
| 646 | c |
|---|
| 647 | return |
|---|
| 648 | end |
|---|
