Line data Source code
1 :
2 : module hk_conv
3 : !
4 : ! Moist convection. Primarily data used by both Zhang-McFarlane convection
5 : ! and Hack shallow convective schemes.
6 : !
7 : ! $Id$
8 : !
9 : use shr_kind_mod, only: r8 => shr_kind_r8
10 : use cam_logfile, only: iulog
11 : use spmd_utils, only: masterproc
12 : use cam_abortutils, only: endrun
13 : implicit none
14 :
15 : private
16 : save
17 : !
18 : ! Public interfaces
19 : !
20 : public mfinti ! Initialization of data for moist convection
21 : public cmfmca ! Hack shallow convection
22 : public hkconv_readnl ! read hkconv_nl namelist
23 :
24 : !
25 : ! Private data used for Hack shallow convection
26 : !
27 : real(r8), parameter :: unset_r8 = huge(1.0_r8)
28 :
29 : ! Namelist variables
30 : real(r8) :: hkconv_c0 = unset_r8
31 : real(r8) :: hkconv_cmftau = unset_r8
32 :
33 : real(r8) :: hlat ! latent heat of vaporization
34 : real(r8) :: c0 ! rain water autoconversion coefficient set from namelist input hkconv_c0
35 : real(r8) :: betamn ! minimum overshoot parameter
36 : real(r8) :: rhlat ! reciprocal of hlat
37 : real(r8) :: rcp ! reciprocal of cp
38 : real(r8) :: cmftau ! characteristic adjustment time scale set from namelist input hkconv_cmftau
39 : real(r8) :: rhoh2o ! density of liquid water (STP)
40 : real(r8) :: dzmin ! minimum convective depth for precipitation
41 : real(r8) :: tiny ! arbitrary small num used in transport estimates
42 : real(r8) :: eps ! convergence criteria (machine dependent)
43 : real(r8) :: tpmax ! maximum acceptable t perturbation (degrees C)
44 : real(r8) :: shpmax ! maximum acceptable q perturbation (g/g)
45 :
46 : integer :: iloc ! longitude location for diagnostics
47 : integer :: jloc ! latitude location for diagnostics
48 : integer :: nsloc ! nstep for which to produce diagnostics
49 : !
50 : logical :: rlxclm ! logical to relax column versus cloud triplet
51 :
52 : real(r8) cp ! specific heat of dry air
53 : real(r8) grav ! gravitational constant
54 : real(r8) rgrav ! reciprocal of grav
55 : real(r8) rgas ! gas constant for dry air
56 : integer limcnv ! top interface level limit for convection
57 :
58 :
59 :
60 :
61 : contains
62 1536 : subroutine hkconv_readnl(nlfile)
63 :
64 : use namelist_utils, only: find_group_name
65 : use units, only: getunit, freeunit
66 : use mpishorthand
67 :
68 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
69 :
70 : ! Local variables
71 : integer :: unitn, ierr
72 : character(len=*), parameter :: subname = 'hkconv_readnl'
73 :
74 : namelist /hkconv_nl/ hkconv_cmftau, hkconv_c0
75 : !-----------------------------------------------------------------------------
76 :
77 1536 : if (masterproc) then
78 2 : unitn = getunit()
79 2 : open( unitn, file=trim(nlfile), status='old' )
80 2 : call find_group_name(unitn, 'hkconv_nl', status=ierr)
81 2 : if (ierr == 0) then
82 2 : read(unitn, hkconv_nl, iostat=ierr)
83 2 : if (ierr /= 0) then
84 0 : call endrun(subname // ':: ERROR reading namelist')
85 : end if
86 : end if
87 2 : close(unitn)
88 2 : call freeunit(unitn)
89 :
90 : ! set local variables
91 2 : cmftau = hkconv_cmftau
92 2 : c0 = hkconv_c0
93 :
94 : end if
95 :
96 : #ifdef SPMD
97 : ! Broadcast namelist variables
98 1536 : call mpibcast(cmftau, 1, mpir8, 0, mpicom)
99 1536 : call mpibcast(c0, 1, mpir8, 0, mpicom)
100 : #endif
101 :
102 1536 : end subroutine hkconv_readnl
103 :
104 : !================================================================================================
105 :
106 1536 : subroutine mfinti (rair ,cpair ,gravit ,latvap ,rhowtr,limcnv_in )
107 : !-----------------------------------------------------------------------
108 : !
109 : ! Purpose:
110 : ! Initialize moist convective mass flux procedure common block, cmfmca
111 : !
112 : ! Method:
113 : ! <Describe the algorithm(s) used in the routine.>
114 : ! <Also include any applicable external references.>
115 : !
116 : ! Author: J. Hack
117 : !
118 : !-----------------------------------------------------------------------
119 : use spmd_utils, only: masterproc
120 : !------------------------------Arguments--------------------------------
121 : !
122 : ! Input arguments
123 : !
124 : real(r8), intent(in) :: rair ! gas constant for dry air
125 : real(r8), intent(in) :: cpair ! specific heat of dry air
126 : real(r8), intent(in) :: gravit ! acceleration due to gravity
127 : real(r8), intent(in) :: latvap ! latent heat of vaporization
128 : real(r8), intent(in) :: rhowtr ! density of liquid water (STP)
129 : integer, intent(in) :: limcnv_in ! top interface level limit for convection
130 :
131 : ! local variables
132 : character(len=32) :: hgrid ! horizontal grid specifier
133 : !
134 : !-----------------------------------------------------------------------
135 : !
136 : ! Initialize physical constants for moist convective mass flux procedure
137 : !
138 1536 : cp = cpair ! specific heat of dry air
139 1536 : hlat = latvap ! latent heat of vaporization
140 1536 : grav = gravit ! gravitational constant
141 1536 : rgas = rair ! gas constant for dry air
142 1536 : rhoh2o = rhowtr ! density of liquid water (STP)
143 :
144 1536 : limcnv = limcnv_in
145 :
146 : ! Initialize free parameters for moist convective mass flux procedure
147 : ! cmftau - characteristic adjustment time scale
148 : ! c0 - rain water autoconversion coeff (1/m)
149 :
150 1536 : if (masterproc) then
151 2 : write(iulog,*) 'tuning parameters hk_conv: cmftau',cmftau
152 2 : write(iulog,*) 'tuning parameters hk_conv: c0',c0
153 : endif
154 1536 : dzmin = 0.0_r8 ! minimum cloud depth to precipitate (m)
155 1536 : betamn = 0.10_r8 ! minimum overshoot parameter
156 :
157 :
158 1536 : tpmax = 1.50_r8 ! maximum acceptable t perturbation (deg C)
159 1536 : shpmax = 1.50e-3_r8 ! maximum acceptable q perturbation (g/g)
160 1536 : rlxclm = .true. ! logical variable to specify that relaxation
161 : ! time scale should applied to column as
162 : ! opposed to triplets individually
163 : !
164 : ! Initialize miscellaneous (frequently used) constants
165 : !
166 1536 : rhlat = 1.0_r8/hlat ! reciprocal latent heat of vaporization
167 1536 : rcp = 1.0_r8/cp ! reciprocal specific heat of dry air
168 1536 : rgrav = 1.0_r8/grav ! reciprocal gravitational constant
169 : !
170 : ! Initialize diagnostic location information for moist convection scheme
171 : !
172 1536 : iloc = 1 ! longitude point for diagnostic info
173 1536 : jloc = 1 ! latitude point for diagnostic info
174 1536 : nsloc = 1 ! nstep value at which to begin diagnostics
175 : !
176 : ! Initialize other miscellaneous parameters
177 : !
178 1536 : tiny = 1.0e-36_r8 ! arbitrary small number (scalar transport)
179 1536 : eps = 1.0e-13_r8 ! convergence criteria (machine dependent)
180 : !
181 1536 : return
182 : end subroutine mfinti
183 :
184 1495368 : subroutine cmfmca(lchnk ,ncol , &
185 : nstep ,ztodt ,pmid ,pdel , &
186 : rpdel ,zm ,tpert ,qpert ,phis , &
187 : pblh ,t ,q ,cmfdt ,dq , &
188 : cmfmc ,cmfdqr ,cmfsl ,cmflq ,precc , &
189 : qc ,cnt ,cnb ,icwmr ,rliq , &
190 : pmiddry ,pdeldry ,rpdeldry)
191 : !-----------------------------------------------------------------------
192 : !
193 : ! Purpose:
194 : ! Moist convective mass flux procedure:
195 : !
196 : ! Method:
197 : ! If stratification is unstable to nonentraining parcel ascent,
198 : ! complete an adjustment making successive use of a simple cloud model
199 : ! consisting of three layers (sometimes referred to as a triplet)
200 : !
201 : ! Code generalized to allow specification of parcel ("updraft")
202 : ! properties, as well as convective transport of an arbitrary
203 : ! number of passive constituents (see q array). The code
204 : ! is written so the water vapor field is passed independently
205 : ! in the calling list from the block of other transported
206 : ! constituents, even though as currently designed, it is the
207 : ! first component in the constituents field.
208 : !
209 : ! Author: J. Hack
210 : !
211 : ! BAB: changed code to report tendencies in cmfdt and dq, instead of
212 : ! updating profiles. Cmfdq contains water only, made it a local variable
213 : ! made dq (all constituents) the argument.
214 : !
215 : !-----------------------------------------------------------------------
216 :
217 : !#######################################################################
218 : !# #
219 : !# Debugging blocks are marked this way for easy identification #
220 : !# #
221 : !#######################################################################
222 : use constituents, only: pcnst
223 : use constituents, only: cnst_get_type_byind
224 : use ppgrid, only: pcols, pver, pverp
225 : #if ( defined DIAGNS )
226 : use phys_grid, only: get_lat_all_p, get_lon_all_p
227 : #endif
228 : use wv_saturation, only: qsat
229 :
230 : real(r8) ssfac ! supersaturation bound (detrained air)
231 : parameter (ssfac = 1.001_r8)
232 :
233 : !------------------------------Arguments--------------------------------
234 : !
235 : ! Input arguments
236 : !
237 : integer, intent(in) :: lchnk ! chunk identifier
238 : integer, intent(in) :: ncol ! number of atmospheric columns
239 : integer, intent(in) :: nstep ! current time step index
240 :
241 : real(r8), intent(in) :: ztodt ! 2 delta-t (seconds)
242 : real(r8), intent(in) :: pmid(pcols,pver) ! pressure
243 : real(r8), intent(in) :: pdel(pcols,pver) ! delta-p
244 : real(r8), intent(in) :: pmiddry(pcols,pver) ! pressure
245 : real(r8), intent(in) :: pdeldry(pcols,pver) ! delta-p
246 : real(r8), intent(in) :: rpdel(pcols,pver) ! 1./pdel
247 : real(r8), intent(in) :: rpdeldry(pcols,pver) ! 1./pdel
248 : real(r8), intent(in) :: zm(pcols,pver) ! height abv sfc at midpoints
249 : real(r8), intent(in) :: tpert(pcols) ! PBL perturbation theta
250 : real(r8), intent(in) :: qpert(pcols,pcnst) ! PBL perturbation specific humidity
251 : real(r8), intent(in) :: phis(pcols) ! surface geopotential
252 : real(r8), intent(in) :: pblh(pcols) ! PBL height (provided by PBL routine)
253 : real(r8), intent(in) :: t(pcols,pver) ! temperature (t bar)
254 : real(r8), intent(in) :: q(pcols,pver,pcnst) ! specific humidity (sh bar)
255 : !
256 : ! Output arguments
257 : !
258 : real(r8), intent(out) :: cmfdt(pcols,pver) ! dt/dt due to moist convection
259 : real(r8), intent(out) :: cmfmc(pcols,pverp) ! moist convection cloud mass flux
260 : real(r8), intent(out) :: cmfdqr(pcols,pver) ! dq/dt due to convective rainout
261 : real(r8), intent(out) :: cmfsl(pcols,pver ) ! convective lw static energy flux
262 : real(r8), intent(out) :: cmflq(pcols,pver ) ! convective total water flux
263 : real(r8), intent(out) :: precc(pcols) ! convective precipitation rate
264 : ! JJH mod to explicitly export cloud water
265 : real(r8), intent(out) :: qc(pcols,pver) ! dq/dt due to export of cloud water
266 : real(r8), intent(out) :: cnt(pcols) ! top level of convective activity
267 : real(r8), intent(out) :: cnb(pcols) ! bottom level of convective activity
268 : real(r8), intent(out) :: dq(pcols,pver,pcnst) ! constituent tendencies
269 : real(r8), intent(out) :: icwmr(pcols,pver)
270 : real(r8), intent(out) :: rliq(pcols)
271 : !
272 : !---------------------------Local workspace-----------------------------
273 : !
274 : real(r8) pm(pcols,pver) ! pressure
275 : real(r8) pd(pcols,pver) ! delta-p
276 : real(r8) rpd(pcols,pver) ! 1./pdel
277 :
278 : real(r8) cmfdq(pcols,pver) ! dq/dt due to moist convection
279 : real(r8) gam(pcols,pver) ! 1/cp (d(qsat)/dT)
280 : real(r8) sb(pcols,pver) ! dry static energy (s bar)
281 : real(r8) hb(pcols,pver) ! moist static energy (h bar)
282 : real(r8) shbs(pcols,pver) ! sat. specific humidity (sh bar star)
283 : real(r8) hbs(pcols,pver) ! sat. moist static energy (h bar star)
284 : real(r8) shbh(pcols,pverp) ! specific humidity on interfaces
285 : real(r8) sbh(pcols,pverp) ! s bar on interfaces
286 : real(r8) hbh(pcols,pverp) ! h bar on interfaces
287 : real(r8) cmrh(pcols,pverp) ! interface constituent mixing ratio
288 : real(r8) prec(pcols) ! instantaneous total precipitation
289 : real(r8) dzcld(pcols) ! depth of convective layer (m)
290 : real(r8) beta(pcols) ! overshoot parameter (fraction)
291 : real(r8) betamx(pcols) ! local maximum on overshoot
292 : real(r8) eta(pcols) ! convective mass flux (kg/m^2 s)
293 : real(r8) etagdt(pcols) ! eta*grav*dt
294 : real(r8) cldwtr(pcols) ! cloud water (mass)
295 : real(r8) rnwtr(pcols) ! rain water (mass)
296 : ! JJH extension to facilitate export of cloud liquid water
297 : real(r8) totcond(pcols) ! total condensate; mix of precip and cloud water (mass)
298 : real(r8) sc (pcols) ! dry static energy ("in-cloud")
299 : real(r8) shc (pcols) ! specific humidity ("in-cloud")
300 : real(r8) hc (pcols) ! moist static energy ("in-cloud")
301 : real(r8) cmrc(pcols) ! constituent mix rat ("in-cloud")
302 : real(r8) dq1(pcols) ! shb convective change (lower lvl)
303 : real(r8) dq2(pcols) ! shb convective change (mid level)
304 : real(r8) dq3(pcols) ! shb convective change (upper lvl)
305 : real(r8) ds1(pcols) ! sb convective change (lower lvl)
306 : real(r8) ds2(pcols) ! sb convective change (mid level)
307 : real(r8) ds3(pcols) ! sb convective change (upper lvl)
308 : real(r8) dcmr1(pcols) ! q convective change (lower lvl)
309 : real(r8) dcmr2(pcols) ! q convective change (mid level)
310 : real(r8) dcmr3(pcols) ! q convective change (upper lvl)
311 : real(r8) estemp(pcols,pver) ! saturation vapor pressure (scratch)
312 : real(r8) vtemp1(2*pcols) ! intermediate scratch vector
313 : real(r8) vtemp2(2*pcols) ! intermediate scratch vector
314 : real(r8) vtemp3(2*pcols) ! intermediate scratch vector
315 : real(r8) vtemp4(2*pcols) ! intermediate scratch vector
316 : real(r8) vtemp5(2*pcols) ! intermediate scratch vector
317 : integer indx1(pcols) ! longitude indices for condition true
318 : logical etagt0 ! true if eta > 0.0
319 : real(r8) sh1 ! dummy arg in qhalf statement func.
320 : real(r8) sh2 ! dummy arg in qhalf statement func.
321 : real(r8) shbs1 ! dummy arg in qhalf statement func.
322 : real(r8) shbs2 ! dummy arg in qhalf statement func.
323 : real(r8) cats ! modified characteristic adj. time
324 : real(r8) rtdt ! 1./ztodt
325 : real(r8) qprime ! modified specific humidity pert.
326 : real(r8) tprime ! modified thermal perturbation
327 : real(r8) pblhgt ! bounded pbl height (max[pblh,1m])
328 : real(r8) fac1 ! intermediate scratch variable
329 : real(r8) shprme ! intermediate specific humidity pert.
330 : real(r8) qsattp ! sat mix rat for thermally pert PBL parcels
331 : real(r8) dz ! local layer depth
332 : real(r8) temp1 ! intermediate scratch variable
333 : real(r8) b1 ! bouyancy measure in detrainment lvl
334 : real(r8) b2 ! bouyancy measure in condensation lvl
335 : real(r8) temp2 ! intermediate scratch variable
336 : real(r8) temp3 ! intermediate scratch variable
337 : real(r8) g ! bounded vertical gradient of hb
338 : real(r8) tmass ! total mass available for convective exch
339 : real(r8) denom ! intermediate scratch variable
340 : real(r8) qtest1 ! used in negative q test (middle lvl)
341 : real(r8) qtest2 ! used in negative q test (lower lvl)
342 : real(r8) fslkp ! flux lw static energy (bot interface)
343 : real(r8) fslkm ! flux lw static energy (top interface)
344 : real(r8) fqlkp ! flux total water (bottom interface)
345 : real(r8) fqlkm ! flux total water (top interface)
346 : real(r8) botflx ! bottom constituent mixing ratio flux
347 : real(r8) topflx ! top constituent mixing ratio flux
348 : real(r8) efac1 ! ratio q to convectively induced chg (btm lvl)
349 : real(r8) efac2 ! ratio q to convectively induced chg (mid lvl)
350 : real(r8) efac3 ! ratio q to convectively induced chg (top lvl)
351 : real(r8) tb(pcols,pver) ! working storage for temp (t bar)
352 : real(r8) shb(pcols,pver) ! working storage for spec hum (sh bar)
353 : real(r8) adjfac ! adjustment factor (relaxation related)
354 : real(r8) rktp
355 : real(r8) rk
356 : #if ( defined DIAGNS )
357 : !
358 : ! Following 7 real variables are used in diagnostics calculations
359 : !
360 : real(r8) rh ! relative humidity
361 : real(r8) es ! sat vapor pressure
362 : real(r8) hsum1 ! moist static energy integral
363 : real(r8) qsum1 ! total water integral
364 : real(r8) hsum2 ! final moist static energy integral
365 : real(r8) qsum2 ! final total water integral
366 : real(r8) fac ! intermediate scratch variable
367 : #endif
368 : integer i,k ! longitude, level indices
369 : integer ii ! index on "gathered" vectors
370 : integer len1 ! vector length of "gathered" vectors
371 : integer m ! constituent index
372 : integer ktp ! tmp indx used to track top of convective layer
373 : #if ( defined DIAGNS )
374 : integer n ! vertical index (diagnostics)
375 : integer kp ! vertical index (diagnostics)
376 : integer kpp ! index offset, kp+1 (diagnostics)
377 : integer kpm1 ! index offset, kp-1 (diagnostics)
378 : integer lat(pcols) ! latitude indices
379 : integer lon(pcols) ! longitude indices
380 : #endif
381 : !
382 : !---------------------------Statement functions-------------------------
383 : !
384 : real(r8) qhalf
385 : qhalf(sh1,sh2,shbs1,shbs2) = min(max(sh1,sh2),(shbs2*sh1 + shbs1*sh2)/(shbs1+shbs2))
386 : !
387 : !-----------------------------------------------------------------------
388 :
389 : !** BAB initialize output tendencies here
390 : ! copy q to dq; use dq below for passive tracer transport
391 650693736 : cmfdt(:ncol,:) = 0._r8
392 650693736 : cmfdq(:ncol,:) = 0._r8
393 1302882840 : dq(:ncol,:,2:) = q(:ncol,:,2:)
394 675662904 : cmfmc(:ncol,:) = 0._r8
395 650693736 : cmfdqr(:ncol,:) = 0._r8
396 650693736 : cmfsl(:ncol,:) = 0._r8
397 650693736 : cmflq(:ncol,:) = 0._r8
398 650693736 : qc(:ncol,:) = 0._r8
399 24969168 : rliq(:ncol) = 0._r8
400 : !
401 : #if ( defined DIAGNS )
402 : ! Determine chunk latitudes and longitudes
403 : call get_lat_all_p(lchnk, ncol, lat)
404 : call get_lon_all_p(lchnk, ncol, lon)
405 : #endif
406 : !
407 : ! Ensure that characteristic adjustment time scale (cmftau) assumed
408 : ! in estimate of eta isn't smaller than model time scale (ztodt)
409 : ! The time over which the convection is assumed to act (the adjustment
410 : ! time scale) can be applied with each application of the three-level
411 : ! cloud model, or applied to the column tendencies after a "hard"
412 : ! adjustment (i.e., on a 2-delta t time scale) is evaluated
413 : !
414 1495368 : if (rlxclm) then
415 1495368 : cats = ztodt ! relaxation applied to column
416 1495368 : adjfac = ztodt/(max(ztodt,cmftau))
417 : else
418 0 : cats = max(ztodt,cmftau) ! relaxation applied to triplet
419 0 : adjfac = 1.0_r8
420 : endif
421 1495368 : rtdt = 1.0_r8/ztodt
422 : !
423 : ! Move temperature and moisture into working storage
424 : !
425 34393464 : do k=limcnv,pver
426 550817064 : do i=1,ncol
427 516423600 : tb (i,k) = t(i,k)
428 549321696 : shb(i,k) = q(i,k,1)
429 : end do
430 : end do
431 40374936 : do k=1,pver
432 650693736 : do i=1,ncol
433 649198368 : icwmr(i,k) = 0._r8
434 : end do
435 : end do
436 : !
437 : ! Compute sb,hb,shbs,hbs
438 : !
439 34393464 : do k = limcnv,pver
440 : call qsat(tb(1:ncol,k), pmid(1:ncol,k), &
441 : estemp(1:ncol,k), shbs(1:ncol,k), ncol, &
442 34393464 : gam=gam(1:ncol,k))
443 : end do
444 : !
445 34393464 : do k=limcnv,pver
446 550817064 : do i=1,ncol
447 516423600 : sb (i,k) = cp*tb(i,k) + zm(i,k)*grav + phis(i)
448 516423600 : hb (i,k) = sb(i,k) + hlat*shb(i,k)
449 549321696 : hbs(i,k) = sb(i,k) + hlat*shbs(i,k)
450 : end do
451 : end do
452 : !
453 : ! Compute sbh, shbh
454 : !
455 32898096 : do k=limcnv+1,pver
456 525847896 : do i=1,ncol
457 492949800 : sbh (i,k) = 0.5_r8*(sb(i,k-1) + sb(i,k))
458 492949800 : shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
459 524352528 : hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
460 : end do
461 : end do
462 : !
463 : ! Specify properties at top of model (not used, but filling anyway)
464 : !
465 24969168 : do i=1,ncol
466 23473800 : sbh (i,limcnv) = sb(i,limcnv)
467 23473800 : shbh(i,limcnv) = shb(i,limcnv)
468 24969168 : hbh (i,limcnv) = hb(i,limcnv)
469 : end do
470 : !
471 : ! Zero vertically independent control, tendency & diagnostic arrays
472 : !
473 24969168 : do i=1,ncol
474 23473800 : prec(i) = 0.0_r8
475 23473800 : dzcld(i) = 0.0_r8
476 23473800 : cnb(i) = 0.0_r8
477 24969168 : cnt(i) = real(pver+1,r8)
478 : end do
479 : #if ( defined DIAGNS )
480 : !#######################################################################
481 : !# #
482 : !# output initial thermodynamic profile if debug diagnostics #
483 : !# #
484 : do i=1,ncol
485 : if ((lat(i).eq.jloc) .and. (lon(i).eq.iloc) &
486 : .and. (nstep.ge.nsloc)) then
487 : !# #
488 : !# approximate vertical integral of moist static energy #
489 : !# and total preciptable water #
490 : !# #
491 : hsum1 = 0.0_r8
492 : qsum1 = 0.0_r8
493 : do k=limcnv,pver
494 : hsum1 = hsum1 + pdel(i,k)*rgrav*hb(i,k)
495 : qsum1 = qsum1 + pdel(i,k)*rgrav*shb(i,k)
496 : end do
497 : !# #
498 : write(iulog,8010)
499 : fac = grav*864._r8
500 : do k=limcnv,pver
501 : rh = shb(i,k)/shbs(i,k)
502 : write(iulog,8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc(i,k),cmfsl(i,k), cmflq(i,k)
503 : write(iulog,8040) tb(i,k),shb(i,k),rh,sb(i,k),hb(i,k),hbs(i,k),ztodt*cmfdt(i,k), &
504 : ztodt*cmfdq(i,k),ztodt*cmfdqr(i,k)
505 : end do
506 : write(iulog, 8000) prec(i)
507 : end if
508 : enddo
509 : #endif
510 : !# #
511 : !# #
512 : !#######################################################################
513 : !
514 : ! Begin moist convective mass flux adjustment procedure.
515 : ! Formalism ensures that negative cloud liquid water can never occur
516 : !
517 31402728 : do 70 k=pver-1,limcnv+1,-1
518 499383360 : do 10 i=1,ncol
519 469476000 : etagdt(i) = 0.0_r8
520 469476000 : eta (i) = 0.0_r8
521 469476000 : beta (i) = 0.0_r8
522 469476000 : ds1 (i) = 0.0_r8
523 469476000 : ds2 (i) = 0.0_r8
524 469476000 : ds3 (i) = 0.0_r8
525 469476000 : dq1 (i) = 0.0_r8
526 469476000 : dq2 (i) = 0.0_r8
527 469476000 : dq3 (i) = 0.0_r8
528 : !
529 : ! Specification of "cloud base" conditions
530 : !
531 469476000 : qprime = 0.0_r8
532 469476000 : tprime = 0.0_r8
533 : !
534 : ! Assign tprime within the PBL to be proportional to the quantity
535 : ! tpert (which will be bounded by tpmax), passed to this routine by
536 : ! the PBL routine. Don't allow perturbation to produce a dry
537 : ! adiabatically unstable parcel. Assign qprime within the PBL to be
538 : ! an appropriately modified value of the quantity qpert (which will be
539 : ! bounded by shpmax) passed to this routine by the PBL routine. The
540 : ! quantity qprime should be less than the local saturation value
541 : ! (qsattp=qsat[t+tprime,p]). In both cases, tpert and qpert are
542 : ! linearly reduced toward zero as the PBL top is approached.
543 : !
544 469476000 : pblhgt = max(pblh(i),1.0_r8)
545 469476000 : if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0_r8 ) then
546 54440113 : fac1 = max(0.0_r8,1.0_r8-zm(i,k+1)/pblhgt)
547 54440113 : tprime = min(tpert(i),tpmax)*fac1
548 54440113 : qsattp = shbs(i,k+1) + cp*rhlat*gam(i,k+1)*tprime
549 54440113 : shprme = min(min(qpert(i,1),shpmax)*fac1,max(qsattp-shb(i,k+1),0.0_r8))
550 54440113 : qprime = max(qprime,shprme)
551 : else
552 : tprime = 0.0_r8
553 : qprime = 0.0_r8
554 : end if
555 : !
556 : ! Specify "updraft" (in-cloud) thermodynamic properties
557 : !
558 469476000 : sc (i) = sb (i,k+1) + cp*tprime
559 469476000 : shc(i) = shb(i,k+1) + qprime
560 469476000 : hc (i) = sc (i ) + hlat*shc(i)
561 469476000 : vtemp4(i) = hc(i) - hbs(i,k)
562 469476000 : dz = pdel(i,k)*rgas*tb(i,k)*rgrav/pmid(i,k)
563 469476000 : if (vtemp4(i) > 0.0_r8) then
564 32205264 : dzcld(i) = dzcld(i) + dz
565 : else
566 437270736 : dzcld(i) = 0.0_r8
567 : end if
568 29907360 : 10 continue
569 : #if ( defined DIAGNS )
570 : !#######################################################################
571 : !# #
572 : !# output thermodynamic perturbation information #
573 : !# #
574 : do i=1,ncol
575 : if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
576 : write(iulog,8090) k+1,sc(iloc),shc(iloc),hc(iloc)
577 : end if
578 : enddo
579 : !# #
580 : !#######################################################################
581 : #endif
582 : !
583 : ! Check on moist convective instability
584 : ! Build index vector of points where instability exists
585 : !
586 : len1 = 0
587 499383360 : do i=1,ncol
588 499383360 : if (vtemp4(i) > 0.0_r8) then
589 32205264 : len1 = len1 + 1
590 32205264 : indx1(len1) = i
591 : end if
592 : end do
593 29907360 : if (len1 <= 0) go to 70
594 : !
595 : ! Current level just below top level => no overshoot
596 : !
597 5008903 : if (k <= limcnv+1) then
598 0 : do ii=1,len1
599 0 : i = indx1(ii)
600 0 : temp1 = vtemp4(i)/(1.0_r8 + gam(i,k))
601 0 : cldwtr(i) = max(0.0_r8,(sb(i,k) - sc(i) + temp1))
602 0 : beta(i) = 0.0_r8
603 0 : vtemp3(i) = (1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k))
604 : end do
605 : else
606 : !
607 : ! First guess at overshoot parameter using crude buoyancy closure
608 : ! 10% overshoot assumed as a minimum and 1-c0*dz maximum to start
609 : ! If pre-existing supersaturation in detrainment layer, beta=0
610 : ! cldwtr is temporarily equal to hlat*l (l=> liquid water)
611 : !
612 37214167 : do ii=1,len1
613 32205264 : i = indx1(ii)
614 32205264 : temp1 = vtemp4(i)/(1.0_r8 + gam(i,k))
615 32205264 : cldwtr(i) = max(0.0_r8,(sb(i,k)-sc(i)+temp1))
616 32205264 : betamx(i) = 1.0_r8 - c0*max(0.0_r8,(dzcld(i)-dzmin))
617 32205264 : b1 = (hc(i) - hbs(i,k-1))*pdel(i,k-1)
618 32205264 : b2 = (hc(i) - hbs(i,k ))*pdel(i,k )
619 32205264 : beta(i) = max(betamn,min(betamx(i), 1.0_r8 + b1/b2))
620 32205264 : if (hbs(i,k-1) <= hb(i,k-1)) beta(i) = 0.0_r8
621 : !
622 : ! Bound maximum beta to ensure physically realistic solutions
623 : !
624 : ! First check constrains beta so that eta remains positive
625 : ! (assuming that eta is already positive for beta equal zero)
626 : !
627 32205264 : vtemp1(i) = -(hbh(i,k+1) - hc(i))*pdel(i,k)*rpdel(i,k+1)+ &
628 32205264 : (1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i))
629 32205264 : vtemp2(i) = (1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k))
630 32205264 : vtemp3(i) = vtemp2(i)
631 37214167 : if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
632 189 : betamx(i) = 0.99_r8*(vtemp1(i)/vtemp2(i))
633 189 : beta(i) = max(0.0_r8,min(betamx(i),beta(i)))
634 : end if
635 : end do
636 : !
637 : ! Second check involves supersaturation of "detrainment layer"
638 : ! small amount of supersaturation acceptable (by ssfac factor)
639 : !
640 37214167 : do ii=1,len1
641 32205264 : i = indx1(ii)
642 37214167 : if (hb(i,k-1) < hbs(i,k-1)) then
643 28918269 : vtemp1(i) = vtemp1(i)*rpdel(i,k)
644 28918269 : temp2 = gam(i,k-1)*(sbh(i,k) - sc(i) + cldwtr(i)) - &
645 28918269 : hbh(i,k) + hc(i) - sc(i) + sbh(i,k)
646 28918269 : temp3 = vtemp3(i)*rpdel(i,k)
647 : vtemp2(i) = (ztodt/cats)*(hc(i) - hbs(i,k))*temp2/ &
648 28918269 : (pdel(i,k-1)*(hbs(i,k-1) - hb(i,k-1))) + temp3
649 28918269 : if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
650 3169407 : betamx(i) = ssfac*(vtemp1(i)/vtemp2(i))
651 3169407 : beta(i) = max(0.0_r8,min(betamx(i),beta(i)))
652 : end if
653 : else
654 3286995 : beta(i) = 0.0_r8
655 : end if
656 : end do
657 : !
658 : ! Third check to avoid introducing 2 delta x thermodynamic
659 : ! noise in the vertical ... constrain adjusted h (or theta e)
660 : ! so that the adjustment doesn't contribute to "kinks" in h
661 : !
662 37214167 : do ii=1,len1
663 32205264 : i = indx1(ii)
664 32205264 : g = min(0.0_r8,hb(i,k) - hb(i,k-1))
665 32205264 : temp1 = (hb(i,k) - hb(i,k-1) - g)*(cats/ztodt)/(hc(i) - hbs(i,k))
666 32205264 : vtemp1(i) = temp1*vtemp1(i) + (hc(i) - hbh(i,k+1))*rpdel(i,k)
667 32205264 : vtemp2(i) = temp1*vtemp3(i)*rpdel(i,k) + (hc(i) - hbh(i,k) - cldwtr(i))* &
668 32205264 : (rpdel(i,k) + rpdel(i,k+1))
669 37214167 : if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._r8) then
670 2014555 : if (vtemp2(i) /= 0.0_r8) then
671 2014555 : betamx(i) = vtemp1(i)/vtemp2(i)
672 : else
673 0 : betamx(i) = 0.0_r8
674 : end if
675 2014555 : beta(i) = max(0.0_r8,min(betamx(i),beta(i)))
676 : end if
677 : end do
678 : end if
679 : !
680 : ! Calculate mass flux required for stabilization.
681 : !
682 : ! Ensure that the convective mass flux, eta, is positive by
683 : ! setting negative values of eta to zero..
684 : ! Ensure that estimated mass flux cannot move more than the
685 : ! minimum of total mass contained in either layer k or layer k+1.
686 : ! Also test for other pathological cases that result in non-
687 : ! physical states and adjust eta accordingly.
688 : !
689 37214167 : do ii=1,len1
690 32205264 : i = indx1(ii)
691 32205264 : beta(i) = max(0.0_r8,beta(i))
692 32205264 : temp1 = hc(i) - hbs(i,k)
693 32205264 : temp2 = ((1.0_r8 + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i)) - &
694 64410528 : beta(i)*vtemp3(i))*rpdel(i,k) - (hbh(i,k+1) - hc(i))*rpdel(i,k+1)
695 32205264 : eta(i) = temp1/(temp2*grav*cats)
696 32205264 : tmass = min(pdel(i,k),pdel(i,k+1))*rgrav
697 32205264 : if (eta(i) > tmass*rtdt .or. eta(i) <= 0.0_r8) eta(i) = 0.0_r8
698 : !
699 : ! Check on negative q in top layer (bound beta)
700 : !
701 32205264 : if (shc(i)-shbh(i,k) < 0.0_r8 .and. beta(i)*eta(i) /= 0.0_r8) then
702 19843 : denom = eta(i)*grav*ztodt*(shc(i) - shbh(i,k))*rpdel(i,k-1)
703 19843 : beta(i) = max(0.0_r8,min(-0.999_r8*shb(i,k-1)/denom,beta(i)))
704 : end if
705 : !
706 : ! Check on negative q in middle layer (zero eta)
707 : !
708 : qtest1 = shb(i,k) + eta(i)*grav*ztodt*((shc(i) - shbh(i,k+1)) - &
709 : (1.0_r8 - beta(i))*cldwtr(i)*rhlat - beta(i)*(shc(i) - shbh(i,k)))* &
710 32205264 : rpdel(i,k)
711 32205264 : if (qtest1 <= 0.0_r8) eta(i) = 0.0_r8
712 : !
713 : ! Check on negative q in lower layer (bound eta)
714 : !
715 32205264 : fac1 = -(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
716 32205264 : qtest2 = shb(i,k+1) - eta(i)*grav*ztodt*fac1
717 32205264 : if (qtest2 < 0.0_r8) then
718 0 : eta(i) = 0.99_r8*shb(i,k+1)/(grav*ztodt*fac1)
719 : end if
720 37214167 : etagdt(i) = eta(i)*grav*ztodt
721 : end do
722 : !
723 : #if ( defined DIAGNS )
724 : !#######################################################################
725 : !# #
726 : do i=1,ncol
727 : if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep >= nsloc)) then
728 : write(iulog,8080) beta(iloc), eta(iloc)
729 : end if
730 : enddo
731 : !# #
732 : !#######################################################################
733 : #endif
734 : !
735 : ! Calculate cloud water, rain water, and thermodynamic changes
736 : !
737 37214167 : do 30 ii=1,len1
738 32205264 : i = indx1(ii)
739 32205264 : icwmr(i,k) = cldwtr(i)*rhlat
740 32205264 : cldwtr(i) = etagdt(i)*cldwtr(i)*rhlat*rgrav
741 : ! JJH changes to facilitate export of cloud liquid water --------------------------------
742 32205264 : totcond(i) = (1.0_r8 - beta(i))*cldwtr(i)
743 32205264 : rnwtr(i) = min(totcond(i),c0*(dzcld(i)-dzmin)*cldwtr(i))
744 32205264 : ds1(i) = etagdt(i)*(sbh(i,k+1) - sc(i))*rpdel(i,k+1)
745 32205264 : dq1(i) = etagdt(i)*(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
746 : ds2(i) = (etagdt(i)*(sc(i) - sbh(i,k+1)) + &
747 32205264 : hlat*grav*cldwtr(i) - beta(i)*etagdt(i)*(sc(i) - sbh(i,k)))*rpdel(i,k)
748 : ! JJH change for export of cloud liquid water; must use total condensate
749 : ! since rainwater no longer represents total condensate
750 : dq2(i) = (etagdt(i)*(shc(i) - shbh(i,k+1)) - grav*totcond(i) - beta(i)* &
751 32205264 : etagdt(i)*(shc(i) - shbh(i,k)))*rpdel(i,k)
752 : ds3(i) = beta(i)*(etagdt(i)*(sc(i) - sbh(i,k)) - hlat*grav*cldwtr(i))* &
753 32205264 : rpdel(i,k-1)
754 32205264 : dq3(i) = beta(i)*etagdt(i)*(shc(i) - shbh(i,k))*rpdel(i,k-1)
755 : !
756 : ! Isolate convective fluxes for later diagnostics
757 : !
758 32205264 : fslkp = eta(i)*(sc(i) - sbh(i,k+1))
759 32205264 : fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) - hlat*cldwtr(i)*rtdt)
760 32205264 : fqlkp = eta(i)*(shc(i) - shbh(i,k+1))
761 32205264 : fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k))
762 : !
763 : ! Update thermodynamic profile (update sb, hb, & hbs later)
764 : !
765 32205264 : tb (i,k+1) = tb(i,k+1) + ds1(i)*rcp
766 32205264 : tb (i,k ) = tb(i,k ) + ds2(i)*rcp
767 32205264 : tb (i,k-1) = tb(i,k-1) + ds3(i)*rcp
768 32205264 : shb(i,k+1) = shb(i,k+1) + dq1(i)
769 32205264 : shb(i,k ) = shb(i,k ) + dq2(i)
770 32205264 : shb(i,k-1) = shb(i,k-1) + dq3(i)
771 : !
772 : ! ** Update diagnostic information for final budget **
773 : ! Tracking precipitation, temperature & specific humidity tendencies,
774 : ! rainout term, convective mass flux, convective liquid
775 : ! water static energy flux, and convective total water flux
776 : ! The variable afac makes the necessary adjustment to the
777 : ! diagnostic fluxes to account for adjustment time scale based on
778 : ! how relaxation time scale is to be applied (column vs. triplet)
779 : !
780 32205264 : prec(i) = prec(i) + (rnwtr(i)/rhoh2o)*adjfac
781 : !
782 : ! The following variables have units of "units"/second
783 : !
784 32205264 : cmfdt (i,k+1) = cmfdt (i,k+1) + ds1(i)*rtdt*adjfac
785 32205264 : cmfdt (i,k ) = cmfdt (i,k ) + ds2(i)*rtdt*adjfac
786 32205264 : cmfdt (i,k-1) = cmfdt (i,k-1) + ds3(i)*rtdt*adjfac
787 32205264 : cmfdq (i,k+1) = cmfdq (i,k+1) + dq1(i)*rtdt*adjfac
788 32205264 : cmfdq (i,k ) = cmfdq (i,k ) + dq2(i)*rtdt*adjfac
789 32205264 : cmfdq (i,k-1) = cmfdq (i,k-1) + dq3(i)*rtdt*adjfac
790 : ! JJH changes to export cloud liquid water --------------------------------
791 32205264 : qc (i,k ) = (grav*(totcond(i)-rnwtr(i))*rpdel(i,k))*rtdt*adjfac
792 32205264 : cmfdqr(i,k ) = cmfdqr(i,k ) + (grav*rnwtr(i)*rpdel(i,k))*rtdt*adjfac
793 32205264 : cmfmc (i,k+1) = cmfmc (i,k+1) + eta(i)*adjfac
794 32205264 : cmfmc (i,k ) = cmfmc (i,k ) + beta(i)*eta(i)*adjfac
795 : !
796 : ! The following variables have units of w/m**2
797 : !
798 32205264 : cmfsl (i,k+1) = cmfsl (i,k+1) + fslkp*adjfac
799 32205264 : cmfsl (i,k ) = cmfsl (i,k ) + fslkm*adjfac
800 32205264 : cmflq (i,k+1) = cmflq (i,k+1) + hlat*fqlkp*adjfac
801 32205264 : cmflq (i,k ) = cmflq (i,k ) + hlat*fqlkm*adjfac
802 5008903 : 30 continue
803 : !
804 : ! Next, convectively modify passive constituents
805 : ! For now, when applying relaxation time scale to thermal fields after
806 : ! entire column has undergone convective overturning, constituents will
807 : ! be mixed using a "relaxed" value of the mass flux determined above
808 : ! Although this will be inconsistant with the treatment of the thermal
809 : ! fields, it's computationally much cheaper, no more-or-less justifiable,
810 : ! and consistent with how the history tape mass fluxes would be used in
811 : ! an off-line mode (i.e., using an off-line transport model)
812 : !
813 15026709 : do 50 m=2,pcnst ! note: indexing assumes water is first field
814 10017806 : if (cnst_get_type_byind(m).eq.'dry') then
815 0 : pd(:ncol,:) = pdeldry(:ncol,:)
816 0 : rpd(:ncol,:) = rpdeldry(:ncol,:)
817 0 : pm(:ncol,:) = pmiddry(:ncol,:)
818 : else
819 4364572634 : pd(:ncol,:) = pdel(:ncol,:)
820 4364572634 : rpd(:ncol,:) = rpdel(:ncol,:)
821 4374590440 : pm(:ncol,:) = pmid(:ncol,:)
822 : endif
823 74428334 : do 40 ii=1,len1
824 64410528 : i = indx1(ii)
825 : !
826 : ! If any of the reported values of the constituent is negative in
827 : ! the three adjacent levels, nothing will be done to the profile
828 : !
829 64410528 : if ((dq(i,k+1,m) < 0.0_r8) .or. (dq(i,k,m) < 0.0_r8) .or. (dq(i,k-1,m) < 0.0_r8)) go to 40
830 : !
831 : ! Specify constituent interface values (linear interpolation)
832 : !
833 64410528 : cmrh(i,k ) = 0.5_r8*(dq(i,k-1,m) + dq(i,k ,m))
834 64410528 : cmrh(i,k+1) = 0.5_r8*(dq(i,k ,m) + dq(i,k+1,m))
835 : !
836 : ! Specify perturbation properties of constituents in PBL
837 : !
838 64410528 : pblhgt = max(pblh(i),1.0_r8)
839 64410528 : if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0_r8 ) then
840 0 : fac1 = max(0.0_r8,1.0_r8-zm(i,k+1)/pblhgt)
841 0 : cmrc(i) = dq(i,k+1,m) + qpert(i,m)*fac1
842 : else
843 64410528 : cmrc(i) = dq(i,k+1,m)
844 : end if
845 : !
846 : ! Determine fluxes, flux divergence => changes due to convection
847 : ! Logic must be included to avoid producing negative values. A bit
848 : ! messy since there are no a priori assumptions about profiles.
849 : ! Tendency is modified (reduced) when pending disaster detected.
850 : !
851 64410528 : botflx = etagdt(i)*(cmrc(i) - cmrh(i,k+1))*adjfac
852 64410528 : topflx = beta(i)*etagdt(i)*(cmrc(i)-cmrh(i,k))*adjfac
853 64410528 : dcmr1(i) = -botflx*rpd(i,k+1)
854 64410528 : efac1 = 1.0_r8
855 64410528 : efac2 = 1.0_r8
856 64410528 : efac3 = 1.0_r8
857 : !
858 64410528 : if (dq(i,k+1,m)+dcmr1(i) < 0.0_r8) then
859 0 : if ( abs(dcmr1(i)) > 1.e-300_r8 ) then
860 0 : efac1 = max(tiny,abs(dq(i,k+1,m)/dcmr1(i)) - eps)
861 : else
862 0 : efac1 = tiny
863 : endif
864 : end if
865 : !
866 64410528 : if (efac1 == tiny .or. efac1 > 1.0_r8) efac1 = 0.0_r8
867 64410528 : dcmr1(i) = -efac1*botflx*rpd(i,k+1)
868 64410528 : dcmr2(i) = (efac1*botflx - topflx)*rpd(i,k)
869 : !
870 64410528 : if (dq(i,k,m)+dcmr2(i) < 0.0_r8) then
871 1276174 : if ( abs(dcmr2(i)) > 1.e-300_r8 ) then
872 1276174 : efac2 = max(tiny,abs(dq(i,k ,m)/dcmr2(i)) - eps)
873 : else
874 : efac2 = tiny
875 : endif
876 : end if
877 : !
878 64410528 : if (efac2 == tiny .or. efac2 > 1.0_r8) efac2 = 0.0_r8
879 64410528 : dcmr2(i) = (efac1*botflx - efac2*topflx)*rpd(i,k)
880 64410528 : dcmr3(i) = efac2*topflx*rpd(i,k-1)
881 : !
882 64410528 : if ( (dq(i,k-1,m)+dcmr3(i) < 0.0_r8 ) ) then
883 7319589 : if ( abs(dcmr3(i)) > 1.e-300_r8 ) then
884 7319589 : efac3 = max(tiny,abs(dq(i,k-1,m)/dcmr3(i)) - eps)
885 : else
886 : efac3 = tiny
887 : endif
888 : end if
889 : !
890 64410528 : if (efac3 == tiny .or. efac3 > 1.0_r8) efac3 = 0.0_r8
891 64410528 : efac3 = min(efac2,efac3)
892 64410528 : dcmr2(i) = (efac1*botflx - efac3*topflx)*rpd(i,k)
893 64410528 : dcmr3(i) = efac3*topflx*rpd(i,k-1)
894 : !
895 64410528 : dq(i,k+1,m) = dq(i,k+1,m) + dcmr1(i)
896 64410528 : dq(i,k ,m) = dq(i,k ,m) + dcmr2(i)
897 64410528 : dq(i,k-1,m) = dq(i,k-1,m) + dcmr3(i)
898 10017806 : 40 continue
899 5008903 : 50 continue ! end of m=2,pcnst loop
900 : !
901 : ! Constituent modifications complete
902 : !
903 5008903 : if (k == limcnv+1) go to 60
904 : !
905 : ! Complete update of thermodynamic structure at integer levels
906 : ! gather/scatter points that need new values of shbs and gamma
907 : !
908 37214167 : do ii=1,len1
909 32205264 : i = indx1(ii)
910 32205264 : vtemp1(ii ) = tb(i,k)
911 32205264 : vtemp1(ii+len1) = tb(i,k-1)
912 32205264 : vtemp2(ii ) = pmid(i,k)
913 37214167 : vtemp2(ii+len1) = pmid(i,k-1)
914 : end do
915 0 : call qsat(vtemp1(1:2*len1), vtemp2(1:2*len1), &
916 5008903 : vtemp5(1:2*len1), vtemp3(1:2*len1), 2*len1, gam=vtemp4(1:2*len1))
917 37214167 : do ii=1,len1
918 32205264 : i = indx1(ii)
919 32205264 : shbs(i,k ) = vtemp3(ii )
920 32205264 : shbs(i,k-1) = vtemp3(ii+len1)
921 32205264 : gam(i,k ) = vtemp4(ii )
922 32205264 : gam(i,k-1) = vtemp4(ii+len1)
923 32205264 : sb (i,k ) = sb(i,k ) + ds2(i)
924 32205264 : sb (i,k-1) = sb(i,k-1) + ds3(i)
925 32205264 : hb (i,k ) = sb(i,k ) + hlat*shb(i,k )
926 32205264 : hb (i,k-1) = sb(i,k-1) + hlat*shb(i,k-1)
927 32205264 : hbs(i,k ) = sb(i,k ) + hlat*shbs(i,k )
928 37214167 : hbs(i,k-1) = sb(i,k-1) + hlat*shbs(i,k-1)
929 : end do
930 : !
931 : ! Update thermodynamic information at half (i.e., interface) levels
932 : !
933 37214167 : do ii=1,len1
934 32205264 : i = indx1(ii)
935 32205264 : sbh (i,k) = 0.5_r8*(sb(i,k) + sb(i,k-1))
936 32205264 : shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
937 32205264 : hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
938 32205264 : sbh (i,k-1) = 0.5_r8*(sb(i,k-1) + sb(i,k-2))
939 32205264 : shbh(i,k-1) = qhalf(shb(i,k-2),shb(i,k-1),shbs(i,k-2),shbs(i,k-1))
940 37214167 : hbh (i,k-1) = sbh(i,k-1) + hlat*shbh(i,k-1)
941 : end do
942 : !
943 : #if ( defined DIAGNS )
944 : !#######################################################################
945 : !# #
946 : !# this update necessary, only if debugging diagnostics requested #
947 : !# #
948 : do i=1,ncol
949 : if (lat(i) == jloc .and. nstep >= nsloc) then
950 : call qsat(tb(i,k+1), pmid(i,k+1), &
951 : es, shbs(i,k+1), gam=gam(i,k+1))
952 : sb (i,k+1) = sb(i,k+1) + ds1(i)
953 : hb (i,k+1) = sb(i,k+1) + hlat*shb(i,k+1)
954 : hbs(i,k+1) = sb(i,k+1) + hlat*shbs(i,k+1)
955 : kpp = k + 2
956 : if (k+1 == pver) kpp = k + 1
957 : do kp=k+1,kpp
958 : kpm1 = kp-1
959 : sbh(i,kp) = 0.5_r8*(sb(i,kpm1) + sb(i,kp))
960 : shbh(i,kp) = qhalf(shb(i,kpm1),shb(i,kp),shbs(i,kpm1),shbs(i,kp))
961 : hbh(i,kp) = sbh(i,kp) + hlat*shbh(i,kp)
962 : end do
963 : end if
964 : end do
965 : !# #
966 : !# diagnostic output #
967 : !# #
968 : do i=1,ncol
969 : if ((lat(i)== jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
970 : write(iulog, 8060) k
971 : fac = grav*864._r8
972 : do n=limcnv,pver
973 : rh = shb(i,n)/shbs(i,n)
974 : write(iulog,8020)shbh(i,n),sbh(i,n),hbh(i,n),fac*cmfmc(i,n), &
975 : cmfsl(i,n), cmflq(i,n)
976 : !--------------write(iulog, 8050)
977 : !--------------write(iulog, 8030) fac*cmfmc(i,n),cmfsl(i,n), cmflq(i,n)
978 : write(iulog, 8040) tb(i,n),shb(i,n),rh,sb(i,n),hb(i,n), &
979 : hbs(i,n), ztodt*cmfdt(i,n),ztodt*cmfdq(i,n), &
980 : ztodt*cmfdqr(i,n)
981 : end do
982 : write(iulog, 8000) prec(i)
983 : end if
984 : end do
985 : !# #
986 : !# #
987 : !#######################################################################
988 : #endif
989 : !
990 : ! Ensure that dzcld is reset if convective mass flux zero
991 : ! specify the current vertical extent of the convective activity
992 : ! top of convective layer determined by size of overshoot param.
993 : !
994 83741439 : 60 do i=1,ncol
995 78732536 : etagt0 = eta(i).gt.0.0_r8
996 78732536 : if ( .not. etagt0) dzcld(i) = 0.0_r8
997 78732536 : if (etagt0 .and. beta(i) > betamn) then
998 15727755 : ktp = k-1
999 : else
1000 : ktp = k
1001 : end if
1002 108639896 : if (etagt0) then
1003 32182355 : rk=k
1004 32182355 : rktp=ktp
1005 32182355 : cnt(i) = min(cnt(i),rktp)
1006 32182355 : cnb(i) = max(cnb(i),rk)
1007 : end if
1008 : end do
1009 1495368 : 70 continue ! end of k loop
1010 : !
1011 : ! ** apply final thermodynamic tendencies **
1012 : !
1013 : !**BAB don't update input profiles
1014 : !!$ do k=limcnv,pver
1015 : !!$ do i=1,ncol
1016 : !!$ t (i,k) = t (i,k) + cmfdt(i,k)*ztodt
1017 : !!$ q(i,k,1) = q(i,k,1) + cmfdq(i,k)*ztodt
1018 : !!$ end do
1019 : !!$ end do
1020 : ! Set output q tendencies
1021 650693736 : dq(:ncol,:,1 ) = cmfdq(:ncol,:)
1022 1302882840 : dq(:ncol,:,2:) = (dq(:ncol,:,2:) - q(:ncol,:,2:))/ztodt
1023 : !
1024 : ! Kludge to prevent cnb-cnt from being zero (in the event
1025 : ! someone decides that they want to divide by this quantity)
1026 : !
1027 24969168 : do i=1,ncol
1028 24969168 : if (cnb(i) /= 0.0_r8 .and. cnb(i) == cnt(i)) then
1029 5544475 : cnt(i) = cnt(i) - 1.0_r8
1030 : end if
1031 : end do
1032 : !
1033 24969168 : do i=1,ncol
1034 24969168 : precc(i) = prec(i)*rtdt
1035 : end do
1036 : !
1037 : ! Compute reserved liquid (not yet in cldliq) for energy integrals.
1038 : ! Treat rliq as flux out bottom, to be added back later.
1039 40374936 : do k = 1, pver
1040 650693736 : do i = 1, ncol
1041 649198368 : rliq(i) = rliq(i) + qc(i,k)*pdel(i,k)/grav
1042 : end do
1043 : end do
1044 24969168 : rliq(:ncol) = rliq(:ncol) /1000._r8
1045 :
1046 : #if ( defined DIAGNS )
1047 : !#######################################################################
1048 : !# #
1049 : !# we're done ... show final result if debug diagnostics requested #
1050 : !# #
1051 : do i=1,ncol
1052 : if ((lat(i)==jloc).and.(lon(i)==iloc).and.(nstep>=nsloc)) then
1053 : fac = grav*864._r8
1054 : write(iulog, 8010)
1055 : do k=limcnv,pver
1056 : rh = shb(i,k)/shbs(i,k)
1057 : write(iulog, 8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc(i,k), &
1058 : cmfsl(i,k), cmflq(i,k)
1059 : write(iulog, 8040) tb(i,k),shb(i,k),rh ,sb(i,k),hb(i,k), &
1060 : hbs(i,k), ztodt*cmfdt(i,k),ztodt*cmfdq(i,k), &
1061 : ztodt*cmfdqr(i,k)
1062 : end do
1063 : write(iulog, 8000) prec(i)
1064 : !# #
1065 : !# approximate vertical integral of moist static energy and #
1066 : !# total preciptable water after adjustment and output changes #
1067 : !# #
1068 : hsum2 = 0.0_r8
1069 : qsum2 = 0.0_r8
1070 : do k=limcnv,pver
1071 : hsum2 = hsum2 + pdel(i,k)*rgrav*hb(i,k)
1072 : qsum2 = qsum2 + pdel(i,k)*rgrav*shb(i,k)
1073 : end do
1074 : !# #
1075 : write(iulog,8070) hsum1, hsum2, abs(hsum2-hsum1)/hsum2, &
1076 : qsum1, qsum2, abs(qsum2-qsum1)/qsum2
1077 : end if
1078 : enddo
1079 : !# #
1080 : !# #
1081 : !#######################################################################
1082 : #endif
1083 1495368 : return ! we're all done ... return to calling procedure
1084 : #if ( defined DIAGNS )
1085 : !
1086 : ! Formats
1087 : !
1088 : 8000 format(///,10x,'PREC = ',3pf12.6,/)
1089 : 8010 format('1** TB SHB RH SB', &
1090 : ' HB HBS CAH CAM PRECC ', &
1091 : ' ETA FSL FLQ **', /)
1092 : 8020 format(' ----- ', 9x,3p,f7.3,2x,2p, 9x,-3p,f7.3,2x, &
1093 : f7.3, 37x, 0p,2x,f8.2, 0p,2x,f8.2,2x,f8.2, ' ----- ')
1094 : 8030 format(' ----- ', 0p,82x,f8.2, 0p,2x,f8.2,2x,f8.2, &
1095 : ' ----- ')
1096 : 8040 format(' - - - ',f7.3,2x,3p,f7.3,2x,2p,f7.3,2x,-3p,f7.3,2x, &
1097 : f7.3, 2x,f8.3,2x,0p,f7.3,3p,2x,f7.3,2x,f7.3,30x, &
1098 : ' - - - ')
1099 : 8050 format(' ----- ',110x,' ----- ')
1100 : 8060 format('1 K =>', i4,/, &
1101 : ' TB SHB RH SB', &
1102 : ' HB HBS CAH CAM PREC ', &
1103 : ' ETA FSL FLQ', /)
1104 : 8070 format(' VERTICALLY INTEGRATED MOIST STATIC ENERGY BEFORE, AFTER', &
1105 : ' AND PERCENTAGE DIFFERENCE => ',1p,2e15.7,2x,2p,f7.3,/, &
1106 : ' VERTICALLY INTEGRATED MOISTURE BEFORE, AFTER' &,
1107 : ' AND PERCENTAGE DIFFERENCE => ',1p,2e15_r8.7_r8,2x,2p,f7.3,/)
1108 : 8080 format(' BETA, ETA => ', 1p,2e12.3)
1109 : 8090 format (' k+1, sc, shc, hc => ', 1x, i2, 1p, 3e12.4)
1110 : #endif
1111 : !
1112 : end subroutine cmfmca
1113 : end module hk_conv
|
