Line data Source code
1 : ! Copyright (C) 2024-2025 National Science Foundation-National Center for Atmospheric Research
2 : ! SPDX-License-Identifier: Apache-2.0
3 : !
4 : ! Hack shallow convective scheme.
5 : ! The main subroutine was formerly named "cmfmca", and its initialization "mfinti".
6 : !
7 : ! Original Author: J. Hack
8 : ! CCPPized: Haipeng Lin, October 2024
9 : module hack_convect_shallow
10 : use ccpp_kinds, only: kind_phys
11 : implicit none
12 : private
13 : save
14 :
15 : ! public CCPP-compliant subroutines
16 : public :: hack_convect_shallow_init
17 : public :: hack_convect_shallow_run
18 :
19 : ! namelist variables for tuning of hack shallow convective scheme
20 : real(kind_phys) :: cmftau ! characteristic adjustment time scale [s]
21 : real(kind_phys) :: c0 ! rain water autoconversion coefficient [m-1]
22 :
23 : ! host-model physical constants and shorthands
24 : real(kind_phys) :: cp ! specific heat of dry air [J K-1 kg-1]
25 : real(kind_phys) :: rgas ! gas constant for dry air [J K-1 kg-1]
26 : real(kind_phys) :: grav ! gravitational constant [m s-2]
27 : real(kind_phys) :: hlat ! latent heat of vaporization [J kg-1]
28 : real(kind_phys) :: rhoh2o ! density of liquid water at STP [kg m-3]
29 :
30 : real(kind_phys) :: rcp ! reciprocal of cp
31 : real(kind_phys) :: rgrav ! reciprocal of grav
32 : real(kind_phys) :: rhlat ! reciprocal of hlat
33 :
34 : integer :: limcnv ! top vertical interface level limit for convection [index]
35 : ! derived from reference pressures to below 40 mb
36 :
37 : ! internal parameters
38 : real(kind_phys) :: betamn = 0.10_kind_phys ! minimum overshoot parameter [???]
39 : real(kind_phys) :: dzmin = 0.0_kind_phys ! minimum convective depth for precipitation [m]
40 : logical :: rlxclm = .true. ! control for relaxing column versus cloud triplet (default: true)
41 : ! true: relaxation timescale should be applied to column as opposed to triplets individually
42 : real(kind_phys) :: ssfac = 1.001_kind_phys ! detrained air supersaturation bound [???]
43 :
44 : ! internal parameters for tolerance
45 : real(kind_phys) :: tiny = 1.0e-36_kind_phys ! arbitrary small num in scalar transport estimates
46 : real(kind_phys) :: eps = 1.0e-13_kind_phys ! machine dependent convergence criteria
47 : real(kind_phys) :: tpmax = 1.50_kind_phys ! maximum acceptable T perturbation [K]
48 : real(kind_phys) :: shpmax = 1.50e-3_kind_phys ! maximum acceptable Q perturbation [g g-1]
49 :
50 : ! diagnostic only
51 : logical :: debug_verbose = .false. ! control for debug messages
52 :
53 :
54 : contains
55 : ! Initialization of moist convective mass procedure including namelist read.
56 : !> \section arg_table_hack_convect_shallow_init Argument Table
57 : !! \htmlinclude hack_convect_shallow_init.html
58 0 : subroutine hack_convect_shallow_init( &
59 : pver, &
60 : amIRoot, iulog, &
61 : cmftau_in, c0_in, &
62 : rair, cpair, gravit, latvap, rhoh2o_in, &
63 0 : pref_edge, &
64 0 : use_shfrc, shfrc, &
65 : top_lev, &
66 0 : errmsg, errflg)
67 :
68 : integer, intent(in) :: pver ! number of vertical levels
69 : logical, intent(in) :: amIRoot
70 : integer, intent(in) :: iulog ! log output unit
71 : real(kind_phys), intent(in) :: cmftau_in ! characteristic adjustment time scale [s]
72 : real(kind_phys), intent(in) :: c0_in ! rain water autoconversion coefficient [m-1]
73 : real(kind_phys), intent(in) :: rair ! gas constant for dry air [J K-1 kg-1]
74 : real(kind_phys), intent(in) :: cpair ! specific heat of dry air [J K-1 kg-1]
75 : real(kind_phys), intent(in) :: gravit ! gravitational constant [m s-2]
76 : real(kind_phys), intent(in) :: latvap ! latent heat of vaporization [J kg-1]
77 : real(kind_phys), intent(in) :: rhoh2o_in ! density of liquid water [kg m-3]
78 : real(kind_phys), intent(in) :: pref_edge(:) ! reference pressures at interface [Pa]
79 :
80 : logical, intent(out) :: use_shfrc ! this shallow scheme provides convective cloud fractions? [flag]
81 : real(kind_phys), intent(out) :: shfrc(:,:) ! (dummy) shallow convective cloud fractions calculated in-scheme [fraction]
82 :
83 : integer, intent(out) :: top_lev ! top level for cloud fraction [index]
84 :
85 : character(len=512), intent(out) :: errmsg
86 : integer, intent(out) :: errflg
87 :
88 : ! local variables
89 : integer :: k
90 :
91 0 : errmsg = ''
92 0 : errflg = 0
93 :
94 : ! namelist variables
95 0 : cmftau = cmftau_in
96 0 : c0 = c0_in
97 :
98 0 : if(amIRoot) then
99 0 : write(iulog,*) 'tuning parameters hack_convect_shallow: cmftau',cmftau
100 0 : write(iulog,*) 'tuning parameters hack_convect_shallow: c0',c0
101 : endif
102 :
103 : ! host model physical constants
104 0 : cp = cpair
105 0 : rcp = 1.0_kind_phys/cp
106 0 : hlat = latvap
107 0 : rhlat = 1.0_kind_phys/hlat
108 0 : grav = gravit
109 0 : rgrav = 1.0_kind_phys/gravit
110 0 : rgas = rair
111 0 : rhoh2o = rhoh2o_in
112 :
113 : ! determine limit of shallow convection: regions below 40 mb
114 : ! logic ported from convect_shallow_init with note that this calculation is repeated in the deep
115 : ! convection interface.
116 0 : if(pref_edge(1) >= 4.e3_kind_phys) then
117 0 : limcnv = 1
118 : else
119 0 : limcnv = pver + 1
120 0 : do k = 1, pver
121 0 : if(pref_edge(k) < 4.e3_kind_phys .and. pref_edge(k+1) >= 4.e3_kind_phys) then
122 0 : limcnv = k
123 : endif
124 : enddo
125 : endif
126 :
127 0 : if(amIRoot) then
128 0 : write(iulog,*) "hack_convect_shallow_init: convection will be capped at interface ", limcnv, &
129 0 : "which is ", pref_edge(limcnv), " pascals"
130 : endif
131 :
132 : ! flags for whether this shallow convection scheme
133 : ! calculates and provides convective cloud fractions
134 : ! to convective cloud cover scheme.
135 : !
136 : ! the Hack scheme does not provide this.
137 : ! a dummy shfrc is provided and is never used.
138 0 : use_shfrc = .false.
139 0 : shfrc(:,:) = 0._kind_phys
140 :
141 : ! for Hack shallow convection (CAM4 physics), do not limit cloud fraction
142 : ! (extend all the way to model top)
143 0 : top_lev = 1
144 0 : end subroutine hack_convect_shallow_init
145 :
146 : ! Moist convective mass flux procedure.
147 : !
148 : ! If stratification is unstable to nonentraining parcel ascent,
149 : ! complete an adjustment making successive use of a simple cloud model
150 : ! consisting of three layers (sometimes referred to as a triplet)
151 : !
152 : ! Code generalized to allow specification of parcel ("updraft")
153 : ! properties, as well as convective transport of an arbitrary
154 : ! number of passive constituents (see q array). The code
155 : ! is written so the water vapor field is passed independently
156 : ! in the calling list from the block of other transported
157 : ! constituents, even though as currently designed, it is the
158 : ! first component in the constituents field.
159 : !
160 : ! Reports tendencies in cmfdt and dq instead of updating profiles.
161 : !
162 : ! Original author: J. Hack, BAB
163 : !> \section arg_table_hack_convect_shallow_run Argument Table
164 : !! \htmlinclude hack_convect_shallow_run.html
165 0 : subroutine hack_convect_shallow_run( &
166 : ncol, pver, pcnst, &
167 : iulog, &
168 0 : const_props, &
169 : ztodt, &
170 0 : pmid, pmiddry, &
171 0 : pdel, pdeldry, rpdel, rpdeldry, &
172 0 : zm, &
173 0 : qpert_in, &
174 0 : phis, &
175 0 : pblh, &
176 0 : t, &
177 0 : q, & ! ... below are output arguments:
178 0 : dq, &
179 0 : qc_sh, &
180 0 : cmfdt, &
181 0 : cmfmc_sh, &
182 0 : cmfdqr, &
183 0 : cmfsl, &
184 0 : cmflq, &
185 0 : precc, &
186 0 : cnt_sh, &
187 0 : cnb_sh, &
188 0 : icwmr, &
189 0 : rliq_sh, &
190 0 : scheme_name, &
191 0 : flx_cnd, &
192 0 : errmsg, errflg &
193 : )
194 : ! framework dependency for const_props
195 : use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
196 :
197 : ! dependency to get constituent index
198 : use ccpp_const_utils, only: ccpp_const_get_idx
199 :
200 : ! to_be_ccppized
201 : use wv_saturation, only: qsat
202 :
203 : ! Input arguments
204 : integer, intent(in) :: ncol ! number of atmospheric columns
205 : integer, intent(in) :: pver ! number of vertical levels
206 : integer, intent(in) :: pcnst ! number of ccpp constituents
207 : integer, intent(in) :: iulog ! log output unit
208 : type(ccpp_constituent_prop_ptr_t), &
209 : intent(in) :: const_props(:) ! ccpp constituent properties pointer
210 : real(kind_phys), intent(in) :: ztodt ! physics timestep [s]
211 :
212 : real(kind_phys), intent(in) :: pmid(:,:) ! midpoint pressures [Pa]
213 : real(kind_phys), intent(in) :: pmiddry(:,:) ! dry pressure at midpoints [Pa]
214 : real(kind_phys), intent(in) :: pdel(:,:) ! layer thickness (delta-p) [Pa]
215 : real(kind_phys), intent(in) :: pdeldry(:,:) ! dry layer thickness [Pa]
216 : real(kind_phys), intent(in) :: rpdel(:,:) ! 1.0 / pdel
217 : real(kind_phys), intent(in) :: rpdeldry(:,:) ! 1.0 / pdeldry
218 :
219 : real(kind_phys), intent(in) :: zm(:,:) ! geopotential height at midpoints [m]
220 : real(kind_phys), intent(in) :: qpert_in(:) ! PBL perturbation specific humidity (convective humidity excess) [kg kg-1]
221 : real(kind_phys), intent(in) :: phis(:) ! surface geopotential [m2 s-2]
222 : real(kind_phys), intent(in) :: pblh(:) ! PBL height [m]
223 : real(kind_phys), intent(in) :: t(:,:) ! temperature [K]
224 : real(kind_phys), intent(in) :: q(:,:,:) ! constituents [kg kg-1]
225 :
226 : ! Output arguments
227 : real(kind_phys), intent(out) :: dq(:,:,:) ! constituent tendencies [kg kg-1 s-1]
228 : real(kind_phys), intent(out) :: qc_sh(:,:) ! dq/dt due to export of cloud water / shallow reserved cloud condensate [kg kg-1 s-1]
229 : real(kind_phys), intent(out) :: cmfdt(:,:) ! heating rate (to ptend%s) [J kg-1 s-1]
230 : real(kind_phys), intent(out) :: cmfmc_sh(:,:) ! convective updraft mass flux, shallow [kg s-1 m-2]
231 : real(kind_phys), intent(out) :: cmfdqr(:,:) ! q tendency due to shallow convective rainout [kg kg-1 s-1]
232 : real(kind_phys), intent(out) :: cmfsl(:,:) ! moist shallow convection liquid water static energy flux [W m-2]
233 : real(kind_phys), intent(out) :: cmflq(:,:) ! moist shallow convection total water flux [W m-2]
234 : real(kind_phys), intent(out) :: precc(:) ! shallow convective precipitation rate [m s-1]
235 : integer, intent(out) :: cnt_sh(:) ! top level of shallow convective activity [index]
236 : integer, intent(out) :: cnb_sh(:) ! bottom level of shallow convective activity [index]
237 : real(kind_phys), intent(out) :: icwmr(:,:) ! shallow convection in-cloud water mixing ratio [kg kg-1]
238 : real(kind_phys), intent(out) :: rliq_sh(:) ! vertically-integrated shallow reserved cloud condensate [m s-1]
239 :
240 : character(len=64), intent(out) :: scheme_name ! scheme name
241 : real(kind_phys), intent(out) :: flx_cnd(:) ! net_liquid_and_lwe_ice_fluxes_through_top_and_bottom_of_atmosphere_column [m s-1] for check_energy_chng
242 :
243 : character(len=512), intent(out) :: errmsg
244 : integer, intent(out) :: errflg
245 :
246 : ! Local variables
247 0 : real(kind_phys) :: tpert(ncol) ! PBL perturbation temperature (convective temperature excess) [K]
248 :
249 : character(len=256) :: const_standard_name ! temp: constituent standard name
250 : logical :: const_is_dry ! temp: constituent is dry flag
251 : integer :: const_wv_idx ! temp: index of water vapor
252 :
253 0 : real(kind_phys) :: pm(ncol,pver) ! pressure [Pa]
254 0 : real(kind_phys) :: pd(ncol,pver) ! delta-p [Pa]
255 0 : real(kind_phys) :: rpd(ncol,pver) ! 1./pdel [Pa-1]
256 0 : real(kind_phys) :: cmfdq(ncol,pver) ! dq(wv)/dt due to moist convection (later copied to dq(:,:,const_wv_idx)) [kg kg-1 s-1]
257 0 : real(kind_phys) :: gam(ncol,pver) ! 1/cp (d(qsat)/dT) change in saturation mixing ratio with temp
258 0 : real(kind_phys) :: sb(ncol,pver) ! dry static energy (s bar) [J kg-1]
259 0 : real(kind_phys) :: hb(ncol,pver) ! moist static energy (h bar) [J kg-1]
260 0 : real(kind_phys) :: shbs(ncol,pver) ! sat. specific humidity (sh bar star)
261 0 : real(kind_phys) :: hbs(ncol,pver) ! sat. moist static energy (h bar star)
262 0 : real(kind_phys) :: shbh(ncol,pver+1) ! specific humidity on interfaces
263 0 : real(kind_phys) :: sbh(ncol,pver+1) ! s bar on interfaces
264 0 : real(kind_phys) :: hbh(ncol,pver+1) ! h bar on interfaces
265 0 : real(kind_phys) :: cmrh(ncol,pver+1) ! interface constituent mixing ratio
266 0 : real(kind_phys) :: prec(ncol) ! instantaneous total precipitation
267 0 : real(kind_phys) :: dzcld(ncol) ! depth of convective layer (m)
268 0 : real(kind_phys) :: beta(ncol) ! overshoot parameter (fraction)
269 0 : real(kind_phys) :: betamx(ncol) ! local maximum on overshoot
270 0 : real(kind_phys) :: eta(ncol) ! convective mass flux (kg/m^2 s)
271 0 : real(kind_phys) :: etagdt(ncol) ! eta*grav*dt
272 0 : real(kind_phys) :: cldwtr(ncol) ! cloud water (mass)
273 0 : real(kind_phys) :: rnwtr(ncol) ! rain water (mass)
274 0 : real(kind_phys) :: totcond(ncol) ! total condensate; mix of precip and cloud water (mass)
275 0 : real(kind_phys) :: sc (ncol) ! dry static energy ("in-cloud")
276 0 : real(kind_phys) :: shc (ncol) ! specific humidity ("in-cloud")
277 0 : real(kind_phys) :: hc (ncol) ! moist static energy ("in-cloud")
278 0 : real(kind_phys) :: cmrc(ncol) ! constituent mix rat ("in-cloud")
279 0 : real(kind_phys) :: dq1(ncol) ! shb convective change (lower lvl)
280 0 : real(kind_phys) :: dq2(ncol) ! shb convective change (mid level)
281 0 : real(kind_phys) :: dq3(ncol) ! shb convective change (upper lvl)
282 0 : real(kind_phys) :: ds1(ncol) ! sb convective change (lower lvl)
283 0 : real(kind_phys) :: ds2(ncol) ! sb convective change (mid level)
284 0 : real(kind_phys) :: ds3(ncol) ! sb convective change (upper lvl)
285 0 : real(kind_phys) :: dcmr1(ncol) ! q convective change (lower lvl)
286 0 : real(kind_phys) :: dcmr2(ncol) ! q convective change (mid level)
287 0 : real(kind_phys) :: dcmr3(ncol) ! q convective change (upper lvl)
288 0 : real(kind_phys) :: estemp(ncol,pver) ! saturation vapor pressure (scratch)
289 0 : real(kind_phys) :: vtemp1(2*ncol) ! intermediate scratch vector
290 0 : real(kind_phys) :: vtemp2(2*ncol) ! intermediate scratch vector
291 0 : real(kind_phys) :: vtemp3(2*ncol) ! intermediate scratch vector
292 0 : real(kind_phys) :: vtemp4(2*ncol) ! intermediate scratch vector
293 0 : real(kind_phys) :: vtemp5(2*ncol) ! intermediate scratch vector
294 0 : integer :: indx1(ncol) ! longitude indices for condition true
295 : logical :: etagt0 ! true if eta > 0.0
296 : real(kind_phys) :: cats ! modified characteristic adj. time
297 : real(kind_phys) :: rtdt ! 1./ztodt
298 : real(kind_phys) :: qprime ! modified specific humidity pert.
299 : real(kind_phys) :: tprime ! modified thermal perturbation
300 : real(kind_phys) :: pblhgt ! bounded pbl height (max[pblh,1m])
301 : real(kind_phys) :: fac1 ! intermediate scratch variable
302 : real(kind_phys) :: shprme ! intermediate specific humidity pert.
303 : real(kind_phys) :: qsattp ! sat mix rat for thermally pert PBL parcels
304 : real(kind_phys) :: dz ! local layer depth
305 : real(kind_phys) :: temp1 ! intermediate scratch variable
306 : real(kind_phys) :: b1 ! bouyancy measure in detrainment lvl
307 : real(kind_phys) :: b2 ! bouyancy measure in condensation lvl
308 : real(kind_phys) :: temp2 ! intermediate scratch variable
309 : real(kind_phys) :: temp3 ! intermediate scratch variable
310 : real(kind_phys) :: g ! bounded vertical gradient of hb
311 : real(kind_phys) :: tmass ! total mass available for convective exch
312 : real(kind_phys) :: denom ! intermediate scratch variable
313 : real(kind_phys) :: qtest1 ! used in negative q test (middle lvl)
314 : real(kind_phys) :: qtest2 ! used in negative q test (lower lvl)
315 : real(kind_phys) :: fslkp ! flux lw static energy (bot interface)
316 : real(kind_phys) :: fslkm ! flux lw static energy (top interface)
317 : real(kind_phys) :: fqlkp ! flux total water (bottom interface)
318 : real(kind_phys) :: fqlkm ! flux total water (top interface)
319 : real(kind_phys) :: botflx ! bottom constituent mixing ratio flux
320 : real(kind_phys) :: topflx ! top constituent mixing ratio flux
321 : real(kind_phys) :: efac1 ! ratio q to convectively induced chg (btm lvl)
322 : real(kind_phys) :: efac2 ! ratio q to convectively induced chg (mid lvl)
323 : real(kind_phys) :: efac3 ! ratio q to convectively induced chg (top lvl)
324 0 : real(kind_phys) :: tb(ncol,pver) ! working storage for temp (t bar)
325 0 : real(kind_phys) :: shb(ncol,pver) ! working storage for spec hum (sh bar)
326 : real(kind_phys) :: adjfac ! adjustment factor (relaxation related)
327 : integer :: i,k ! longitude, level indices
328 : integer :: ii ! index on "gathered" vectors
329 : integer :: len1 ! vector length of "gathered" vectors
330 : integer :: m ! constituent index
331 : integer :: ktp ! tmp indx used to track top of convective layer
332 :
333 : ! debug use quantities
334 : real(kind_phys) :: rh ! relative humidity
335 : real(kind_phys) :: es ! sat vapor pressure
336 : real(kind_phys) :: hsum1 ! moist static energy integral
337 : real(kind_phys) :: qsum1 ! total water integral
338 : real(kind_phys) :: hsum2 ! final moist static energy integral
339 : real(kind_phys) :: qsum2 ! final total water integral
340 : real(kind_phys) :: fac ! intermediate scratch variable
341 : integer :: n ! vertical index (diagnostics)
342 : integer :: kp ! vertical index (diagnostics)
343 : integer :: kpp ! index offset, kp+1 (diagnostics)
344 : integer :: kpm1 ! index offset, kp-1 (diagnostics)
345 :
346 0 : errmsg = ''
347 0 : errflg = 0
348 :
349 0 : scheme_name = 'hack_convect_shallow'
350 :
351 : !---------------------------------------------------
352 : ! Initialize output tendencies
353 : !---------------------------------------------------
354 0 : cmfdt (:ncol,:) = 0._kind_phys
355 0 : cmfdq (:ncol,:) = 0._kind_phys
356 0 : cmfmc_sh(:ncol,:) = 0._kind_phys
357 0 : cmfdqr (:ncol,:) = 0._kind_phys
358 0 : cmfsl (:ncol,:) = 0._kind_phys
359 0 : cmflq (:ncol,:) = 0._kind_phys
360 0 : qc_sh (:ncol,:) = 0._kind_phys
361 0 : rliq_sh (:ncol) = 0._kind_phys
362 :
363 : ! Check constituents list and locate water vapor index
364 : ! (not assumed to be 1)
365 : call ccpp_const_get_idx(const_props, &
366 : 'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water', &
367 0 : const_wv_idx, errmsg, errflg)
368 :
369 : !---------------------------------------------------
370 : ! copy q to dq for passive tracer transport.
371 : ! this is NOT an initialization. the dq at this point
372 : ! is not physical (used as temporary here) only at the end
373 : ! dq is updated to be an actual tendency.
374 : !---------------------------------------------------
375 0 : if(pcnst > 1) then
376 : ! set dq for passive tracer transport from q as temporary...
377 0 : dq(:ncol,:,:) = q(:ncol,:,:)
378 :
379 : ! except for water vapor
380 0 : dq(:ncol,:,const_wv_idx) = 0._kind_phys
381 : endif
382 :
383 : !---------------------------------------------------
384 : ! Quantity preparations from convect_shallow.F90.
385 : !---------------------------------------------------
386 :
387 : ! convect_shallow.F90 is not linked to pbuf tpert and always sets to zero.
388 : ! "This field probably should reference the pbuf tpert field but it doesnt"
389 0 : tpert(:ncol) = 0.0_kind_phys
390 :
391 : !---------------------------------------------------
392 : ! Preparation of working arrays
393 : !---------------------------------------------------
394 : ! Ensure that characteristic adjustment time scale (cmftau) assumed
395 : ! in estimate of eta isn't smaller than model time scale (ztodt)
396 : ! The time over which the convection is assumed to act (the adjustment
397 : ! time scale) can be applied with each application of the three-level
398 : ! cloud model, or applied to the column tendencies after a "hard"
399 : ! adjustment (i.e., on a 2-delta t time scale) is evaluated
400 0 : if (rlxclm) then
401 0 : cats = ztodt ! relaxation applied to column
402 0 : adjfac = ztodt/(max(ztodt,cmftau))
403 : else
404 0 : cats = max(ztodt,cmftau) ! relaxation applied to triplet
405 0 : adjfac = 1.0_kind_phys
406 : endif
407 0 : rtdt = 1.0_kind_phys/ztodt
408 :
409 : ! Move temperature and moisture into working storage
410 0 : do k=limcnv,pver
411 0 : do i=1,ncol
412 0 : tb (i,k) = t(i,k)
413 0 : shb(i,k) = q(i,k,const_wv_idx)
414 : end do
415 : end do
416 0 : do k=1,pver
417 0 : do i=1,ncol
418 0 : icwmr(i,k) = 0._kind_phys
419 : end do
420 : end do
421 :
422 : ! Compute sb,hb,shbs,hbs
423 0 : do k = limcnv,pver
424 0 : call qsat(tb(1:ncol,k), pmid(1:ncol,k), &
425 : estemp(1:ncol,k), shbs(1:ncol,k), ncol, &
426 0 : gam=gam(1:ncol,k))
427 : end do
428 :
429 0 : do k=limcnv,pver
430 0 : do i=1,ncol
431 0 : sb (i,k) = cp*tb(i,k) + zm(i,k)*grav + phis(i)
432 0 : hb (i,k) = sb(i,k) + hlat*shb(i,k)
433 0 : hbs(i,k) = sb(i,k) + hlat*shbs(i,k)
434 : end do
435 : end do
436 :
437 : ! Compute sbh, shbh
438 0 : do k=limcnv+1,pver
439 0 : do i=1,ncol
440 0 : sbh (i,k) = 0.5_kind_phys*(sb(i,k-1) + sb(i,k))
441 0 : shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
442 0 : hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
443 : end do
444 : end do
445 :
446 : ! Specify properties at top of model (not used, but filling anyway)
447 0 : do i=1,ncol
448 0 : sbh (i,limcnv) = sb(i,limcnv)
449 0 : shbh(i,limcnv) = shb(i,limcnv)
450 0 : hbh (i,limcnv) = hb(i,limcnv)
451 : end do
452 :
453 : ! Zero vertically independent control, tendency & diagnostic arrays
454 0 : do i=1,ncol
455 0 : prec(i) = 0.0_kind_phys
456 0 : dzcld(i) = 0.0_kind_phys
457 0 : cnb_sh(i)= 0
458 0 : cnt_sh(i)= pver+1
459 : end do
460 :
461 0 : if(debug_verbose) then
462 : ! DEBUG DIAGNOSTICS - Output initial thermodynamic profile
463 0 : do i=1,ncol
464 0 : if(i == 1) then
465 : ! Approximate vertical integral of moist static energy
466 : ! and total precipitable water
467 0 : hsum1 = 0.0_kind_phys
468 0 : qsum1 = 0.0_kind_phys
469 0 : do k=limcnv,pver
470 0 : hsum1 = hsum1 + pdel(i,k)*rgrav*hb(i,k)
471 0 : qsum1 = qsum1 + pdel(i,k)*rgrav*shb(i,k)
472 : end do
473 :
474 0 : write(iulog,8010)
475 0 : fac = grav*864._kind_phys
476 0 : do k=limcnv,pver
477 0 : rh = shb(i,k)/shbs(i,k)
478 0 : write(iulog,8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc_sh(i,k),cmfsl(i,k), cmflq(i,k)
479 0 : write(iulog,8040) tb(i,k),shb(i,k),rh,sb(i,k),hb(i,k),hbs(i,k),ztodt*cmfdt(i,k), &
480 0 : ztodt*cmfdq(i,k),ztodt*cmfdqr(i,k)
481 : end do
482 0 : write(iulog, 8000) prec(i)
483 : end if
484 : end do
485 : endif
486 :
487 : !---------------------------------------------------
488 : ! Begin moist convective mass flux adjustment procedure.
489 : ! Formalism ensures that negative cloud liquid water can never occur
490 : !---------------------------------------------------
491 0 : kloop: do k = pver-1,limcnv+1,-1
492 0 : do i = 1, ncol
493 0 : etagdt(i) = 0.0_kind_phys
494 0 : eta (i) = 0.0_kind_phys
495 0 : beta (i) = 0.0_kind_phys
496 0 : ds1 (i) = 0.0_kind_phys
497 0 : ds2 (i) = 0.0_kind_phys
498 0 : ds3 (i) = 0.0_kind_phys
499 0 : dq1 (i) = 0.0_kind_phys
500 0 : dq2 (i) = 0.0_kind_phys
501 0 : dq3 (i) = 0.0_kind_phys
502 : ! Specification of "cloud base" conditions
503 0 : qprime = 0.0_kind_phys
504 0 : tprime = 0.0_kind_phys
505 :
506 : ! Assign tprime within the PBL to be proportional to the quantity
507 : ! tpert (which will be bounded by tpmax), passed to this routine by
508 : ! the PBL routine. Don't allow perturbation to produce a dry
509 : ! adiabatically unstable parcel. Assign qprime within the PBL to be
510 : ! an appropriately modified value of the quantity qpert (which will be
511 : ! bounded by shpmax) passed to this routine by the PBL routine. The
512 : ! quantity qprime should be less than the local saturation value
513 : ! (qsattp=qsat[t+tprime,p]). In both cases, tpert and qpert are
514 : ! linearly reduced toward zero as the PBL top is approached.
515 0 : pblhgt = max(pblh(i),1.0_kind_phys)
516 0 : if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0_kind_phys ) then
517 0 : fac1 = max(0.0_kind_phys,1.0_kind_phys-zm(i,k+1)/pblhgt)
518 0 : tprime = min(tpert(i),tpmax)*fac1
519 0 : qsattp = shbs(i,k+1) + cp*rhlat*gam(i,k+1)*tprime
520 0 : shprme = min(min(qpert_in(i),shpmax)*fac1,max(qsattp-shb(i,k+1),0.0_kind_phys))
521 0 : qprime = max(qprime,shprme)
522 : else
523 : tprime = 0.0_kind_phys
524 : qprime = 0.0_kind_phys
525 : end if
526 :
527 : ! Specify "updraft" (in-cloud) thermodynamic properties
528 0 : sc (i) = sb (i,k+1) + cp*tprime
529 0 : shc(i) = shb(i,k+1) + qprime
530 0 : hc (i) = sc (i ) + hlat*shc(i)
531 0 : vtemp4(i) = hc(i) - hbs(i,k)
532 0 : dz = pdel(i,k)*rgas*tb(i,k)*rgrav/pmid(i,k)
533 0 : if (vtemp4(i) > 0.0_kind_phys) then
534 0 : dzcld(i) = dzcld(i) + dz
535 : else
536 0 : dzcld(i) = 0.0_kind_phys
537 : end if
538 : enddo
539 :
540 0 : if(debug_verbose) then
541 : ! DEBUG DIAGNOSTICS - output thermodynamic perturbation information
542 0 : do i=1,ncol
543 0 : if(i == 1) then
544 0 : write(iulog,8090) k+1,sc(i),shc(i),hc(i)
545 : end if
546 : enddo
547 : endif
548 :
549 :
550 : ! Check on moist convective instability
551 : ! Build index vector of points where instability exists
552 : len1 = 0
553 0 : do i=1,ncol
554 0 : if (vtemp4(i) > 0.0_kind_phys) then
555 0 : len1 = len1 + 1
556 0 : indx1(len1) = i
557 : end if
558 : end do
559 :
560 0 : if (len1 <= 0) cycle kloop
561 :
562 : ! Current level just below top level => no overshoot
563 0 : if (k <= limcnv+1) then
564 0 : do ii=1,len1
565 0 : i = indx1(ii)
566 0 : temp1 = vtemp4(i)/(1.0_kind_phys + gam(i,k))
567 0 : cldwtr(i) = max(0.0_kind_phys,(sb(i,k) - sc(i) + temp1))
568 0 : beta(i) = 0.0_kind_phys
569 0 : vtemp3(i) = (1.0_kind_phys + gam(i,k))*(sc(i) - sbh(i,k))
570 : end do
571 : else
572 : ! First guess at overshoot parameter using crude buoyancy closure
573 : ! 10% overshoot assumed as a minimum and 1-c0*dz maximum to start
574 : ! If pre-existing supersaturation in detrainment layer, beta=0
575 : ! cldwtr is temporarily equal to hlat*l (l=> liquid water)
576 0 : do ii=1,len1
577 0 : i = indx1(ii)
578 0 : temp1 = vtemp4(i)/(1.0_kind_phys + gam(i,k))
579 0 : cldwtr(i) = max(0.0_kind_phys,(sb(i,k)-sc(i)+temp1))
580 0 : betamx(i) = 1.0_kind_phys - c0*max(0.0_kind_phys,(dzcld(i)-dzmin))
581 0 : b1 = (hc(i) - hbs(i,k-1))*pdel(i,k-1)
582 0 : b2 = (hc(i) - hbs(i,k ))*pdel(i,k )
583 0 : beta(i) = max(betamn,min(betamx(i), 1.0_kind_phys + b1/b2))
584 0 : if (hbs(i,k-1) <= hb(i,k-1)) beta(i) = 0.0_kind_phys
585 :
586 : ! Bound maximum beta to ensure physically realistic solutions
587 : !
588 : ! First check constrains beta so that eta remains positive
589 : ! (assuming that eta is already positive for beta equal zero)
590 0 : vtemp1(i) = -(hbh(i,k+1) - hc(i))*pdel(i,k)*rpdel(i,k+1)+ &
591 0 : (1.0_kind_phys + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i))
592 0 : vtemp2(i) = (1.0_kind_phys + gam(i,k))*(sc(i) - sbh(i,k))
593 0 : vtemp3(i) = vtemp2(i)
594 0 : if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._kind_phys) then
595 0 : betamx(i) = 0.99_kind_phys*(vtemp1(i)/vtemp2(i))
596 0 : beta(i) = max(0.0_kind_phys,min(betamx(i),beta(i)))
597 : end if
598 : end do
599 :
600 : ! Second check involves supersaturation of "detrainment layer"
601 : ! small amount of supersaturation acceptable (by ssfac factor)
602 0 : do ii=1,len1
603 0 : i = indx1(ii)
604 0 : if (hb(i,k-1) < hbs(i,k-1)) then
605 0 : vtemp1(i) = vtemp1(i)*rpdel(i,k)
606 0 : temp2 = gam(i,k-1)*(sbh(i,k) - sc(i) + cldwtr(i)) - &
607 0 : hbh(i,k) + hc(i) - sc(i) + sbh(i,k)
608 0 : temp3 = vtemp3(i)*rpdel(i,k)
609 0 : vtemp2(i) = (ztodt/cats)*(hc(i) - hbs(i,k))*temp2/ &
610 0 : (pdel(i,k-1)*(hbs(i,k-1) - hb(i,k-1))) + temp3
611 0 : if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._kind_phys) then
612 0 : betamx(i) = ssfac*(vtemp1(i)/vtemp2(i))
613 0 : beta(i) = max(0.0_kind_phys,min(betamx(i),beta(i)))
614 : end if
615 : else
616 0 : beta(i) = 0.0_kind_phys
617 : end if
618 : end do
619 :
620 : ! Third check to avoid introducing 2 delta x thermodynamic
621 : ! noise in the vertical ... constrain adjusted h (or theta e)
622 : ! so that the adjustment doesn't contribute to "kinks" in h
623 0 : do ii=1,len1
624 0 : i = indx1(ii)
625 0 : g = min(0.0_kind_phys,hb(i,k) - hb(i,k-1))
626 0 : temp1 = (hb(i,k) - hb(i,k-1) - g)*(cats/ztodt)/(hc(i) - hbs(i,k))
627 0 : vtemp1(i) = temp1*vtemp1(i) + (hc(i) - hbh(i,k+1))*rpdel(i,k)
628 0 : vtemp2(i) = temp1*vtemp3(i)*rpdel(i,k) + (hc(i) - hbh(i,k) - cldwtr(i))* &
629 0 : (rpdel(i,k) + rpdel(i,k+1))
630 0 : if ((beta(i)*vtemp2(i) - vtemp1(i)) > 0._kind_phys) then
631 0 : if (vtemp2(i) /= 0.0_kind_phys) then
632 0 : betamx(i) = vtemp1(i)/vtemp2(i)
633 : else
634 0 : betamx(i) = 0.0_kind_phys
635 : end if
636 0 : beta(i) = max(0.0_kind_phys,min(betamx(i),beta(i)))
637 : end if
638 : end do
639 : end if ! (k <= limcnv+1) Current level just below top level => no overshoot
640 :
641 :
642 : ! Calculate mass flux required for stabilization.
643 : !
644 : ! Ensure that the convective mass flux, eta, is positive by
645 : ! setting negative values of eta to zero..
646 : ! Ensure that estimated mass flux cannot move more than the
647 : ! minimum of total mass contained in either layer k or layer k+1.
648 : ! Also test for other pathological cases that result in non-
649 : ! physical states and adjust eta accordingly.
650 0 : do ii=1,len1
651 0 : i = indx1(ii)
652 0 : beta(i) = max(0.0_kind_phys,beta(i))
653 0 : temp1 = hc(i) - hbs(i,k)
654 0 : temp2 = ((1.0_kind_phys + gam(i,k))*(sc(i) - sbh(i,k+1) + cldwtr(i)) - &
655 0 : beta(i)*vtemp3(i))*rpdel(i,k) - (hbh(i,k+1) - hc(i))*rpdel(i,k+1)
656 0 : eta(i) = temp1/(temp2*grav*cats)
657 0 : tmass = min(pdel(i,k),pdel(i,k+1))*rgrav
658 0 : if (eta(i) > tmass*rtdt .or. eta(i) <= 0.0_kind_phys) eta(i) = 0.0_kind_phys
659 :
660 : ! Check on negative q in top layer (bound beta)
661 0 : if (shc(i)-shbh(i,k) < 0.0_kind_phys .and. beta(i)*eta(i) /= 0.0_kind_phys) then
662 0 : denom = eta(i)*grav*ztodt*(shc(i) - shbh(i,k))*rpdel(i,k-1)
663 0 : beta(i) = max(0.0_kind_phys,min(-0.999_kind_phys*shb(i,k-1)/denom,beta(i)))
664 : end if
665 :
666 : ! Check on negative q in middle layer (zero eta)
667 : qtest1 = shb(i,k) + eta(i)*grav*ztodt*((shc(i) - shbh(i,k+1)) - &
668 : (1.0_kind_phys - beta(i))*cldwtr(i)*rhlat - beta(i)*(shc(i) - shbh(i,k)))* &
669 0 : rpdel(i,k)
670 0 : if (qtest1 <= 0.0_kind_phys) eta(i) = 0.0_kind_phys
671 :
672 : ! Check on negative q in lower layer (bound eta)
673 0 : fac1 = -(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
674 0 : qtest2 = shb(i,k+1) - eta(i)*grav*ztodt*fac1
675 0 : if (qtest2 < 0.0_kind_phys) then
676 0 : eta(i) = 0.99_kind_phys*shb(i,k+1)/(grav*ztodt*fac1)
677 : end if
678 0 : etagdt(i) = eta(i)*grav*ztodt
679 : end do
680 :
681 0 : if(debug_verbose) then
682 0 : do i=1,ncol
683 0 : if (i == 1) then
684 0 : write(iulog,8080) beta(i), eta(i)
685 : endif
686 : enddo
687 : endif
688 :
689 : ! Calculate cloud water, rain water, and thermodynamic changes
690 0 : do ii=1,len1
691 0 : i = indx1(ii)
692 0 : icwmr(i,k) = cldwtr(i)*rhlat
693 0 : cldwtr(i) = etagdt(i)*cldwtr(i)*rhlat*rgrav
694 :
695 : ! JJH changes to facilitate export of cloud liquid water --------------------------------
696 0 : totcond(i) = (1.0_kind_phys - beta(i))*cldwtr(i)
697 0 : rnwtr(i) = min(totcond(i),c0*(dzcld(i)-dzmin)*cldwtr(i))
698 0 : ds1(i) = etagdt(i)*(sbh(i,k+1) - sc(i))*rpdel(i,k+1)
699 0 : dq1(i) = etagdt(i)*(shbh(i,k+1) - shc(i))*rpdel(i,k+1)
700 : ds2(i) = (etagdt(i)*(sc(i) - sbh(i,k+1)) + &
701 0 : hlat*grav*cldwtr(i) - beta(i)*etagdt(i)*(sc(i) - sbh(i,k)))*rpdel(i,k)
702 :
703 : ! JJH change for export of cloud liquid water; must use total condensate
704 : ! since rainwater no longer represents total condensate
705 : dq2(i) = (etagdt(i)*(shc(i) - shbh(i,k+1)) - grav*totcond(i) - beta(i)* &
706 0 : etagdt(i)*(shc(i) - shbh(i,k)))*rpdel(i,k)
707 : ds3(i) = beta(i)*(etagdt(i)*(sc(i) - sbh(i,k)) - hlat*grav*cldwtr(i))* &
708 0 : rpdel(i,k-1)
709 0 : dq3(i) = beta(i)*etagdt(i)*(shc(i) - shbh(i,k))*rpdel(i,k-1)
710 :
711 : ! Isolate convective fluxes for later diagnostics
712 0 : fslkp = eta(i)*(sc(i) - sbh(i,k+1))
713 0 : fslkm = beta(i)*(eta(i)*(sc(i) - sbh(i,k)) - hlat*cldwtr(i)*rtdt)
714 0 : fqlkp = eta(i)*(shc(i) - shbh(i,k+1))
715 0 : fqlkm = beta(i)*eta(i)*(shc(i) - shbh(i,k))
716 :
717 : ! Update thermodynamic profile (update sb, hb, & hbs later)
718 0 : tb (i,k+1) = tb(i,k+1) + ds1(i)*rcp
719 0 : tb (i,k ) = tb(i,k ) + ds2(i)*rcp
720 0 : tb (i,k-1) = tb(i,k-1) + ds3(i)*rcp
721 0 : shb(i,k+1) = shb(i,k+1) + dq1(i)
722 0 : shb(i,k ) = shb(i,k ) + dq2(i)
723 0 : shb(i,k-1) = shb(i,k-1) + dq3(i)
724 :
725 : ! ** Update diagnostic information for final budget **
726 : ! Tracking precipitation, temperature & specific humidity tendencies,
727 : ! rainout term, convective mass flux, convective liquid
728 : ! water static energy flux, and convective total water flux
729 : ! The variable afac makes the necessary adjustment to the
730 : ! diagnostic fluxes to account for adjustment time scale based on
731 : ! how relaxation time scale is to be applied (column vs. triplet)
732 0 : prec(i) = prec(i) + (rnwtr(i)/rhoh2o)*adjfac
733 :
734 : ! The following variables have units of "units"/second
735 0 : cmfdt (i,k+1) = cmfdt (i,k+1) + ds1(i)*rtdt*adjfac
736 0 : cmfdt (i,k ) = cmfdt (i,k ) + ds2(i)*rtdt*adjfac
737 0 : cmfdt (i,k-1) = cmfdt (i,k-1) + ds3(i)*rtdt*adjfac
738 0 : cmfdq (i,k+1) = cmfdq (i,k+1) + dq1(i)*rtdt*adjfac
739 0 : cmfdq (i,k ) = cmfdq (i,k ) + dq2(i)*rtdt*adjfac
740 0 : cmfdq (i,k-1) = cmfdq (i,k-1) + dq3(i)*rtdt*adjfac
741 :
742 : ! JJH changes to export cloud liquid water --------------------------------
743 0 : qc_sh (i,k ) = (grav*(totcond(i)-rnwtr(i))*rpdel(i,k))*rtdt*adjfac
744 0 : cmfdqr (i,k ) = cmfdqr(i,k ) + (grav*rnwtr(i)*rpdel(i,k))*rtdt*adjfac
745 0 : cmfmc_sh(i,k+1) = cmfmc_sh(i,k+1) + eta(i)*adjfac
746 0 : cmfmc_sh(i,k ) = cmfmc_sh(i,k ) + beta(i)*eta(i)*adjfac
747 :
748 : ! The following variables have units of w/m**2
749 0 : cmfsl (i,k+1) = cmfsl (i,k+1) + fslkp*adjfac
750 0 : cmfsl (i,k ) = cmfsl (i,k ) + fslkm*adjfac
751 0 : cmflq (i,k+1) = cmflq (i,k+1) + hlat*fqlkp*adjfac
752 0 : cmflq (i,k ) = cmflq (i,k ) + hlat*fqlkm*adjfac
753 : enddo
754 :
755 : ! Next, convectively modify passive constituents
756 : ! For now, when applying relaxation time scale to thermal fields after
757 : ! entire column has undergone convective overturning, constituents will
758 : ! be mixed using a "relaxed" value of the mass flux determined above
759 : ! Although this will be inconsistant with the treatment of the thermal
760 : ! fields, it's computationally much cheaper, no more-or-less justifiable,
761 : ! and consistent with how the history tape mass fluxes would be used in
762 : ! an off-line mode (i.e., using an off-line transport model)
763 0 : const_modify_loop: do m = 1, pcnst
764 : ! Water vapor needs to be skipped in the loop.
765 0 : if (m == const_wv_idx) then
766 : cycle const_modify_loop
767 : endif
768 :
769 : ! assign pd, rpd, pm temporary properties based on constituent dry/moist mixing ratio
770 0 : call const_props(m)%is_dry(const_is_dry, errflg, errmsg)
771 0 : if(const_is_dry) then
772 0 : pd (:ncol,:) = pdeldry (:ncol,:)
773 0 : rpd(:ncol,:) = rpdeldry(:ncol,:)
774 0 : pm (:ncol,:) = pmiddry (:ncol,:)
775 : else
776 0 : pd (:ncol,:) = pdel (:ncol,:)
777 0 : rpd(:ncol,:) = rpdel (:ncol,:)
778 0 : pm (:ncol,:) = pmid (:ncol,:)
779 : endif
780 :
781 0 : pcl1loop: do ii=1,len1
782 0 : i = indx1(ii)
783 :
784 : ! If any of the reported values of the constituent is negative in
785 : ! the three adjacent levels, nothing will be done to the profile
786 0 : if ((dq(i,k+1,m) < 0.0_kind_phys) .or. (dq(i,k,m) < 0.0_kind_phys) .or. (dq(i,k-1,m) < 0.0_kind_phys)) cycle pcl1loop
787 :
788 : ! Specify constituent interface values (linear interpolation)
789 0 : cmrh(i,k ) = 0.5_kind_phys*(dq(i,k-1,m) + dq(i,k ,m))
790 0 : cmrh(i,k+1) = 0.5_kind_phys*(dq(i,k ,m) + dq(i,k+1,m))
791 :
792 : ! Specify perturbation properties of constituents in PBL
793 0 : pblhgt = max(pblh(i),1.0_kind_phys)
794 0 : if ( (zm(i,k+1) <= pblhgt) .and. dzcld(i) == 0.0_kind_phys ) then
795 0 : fac1 = max(0.0_kind_phys,1.0_kind_phys-zm(i,k+1)/pblhgt)
796 : ! cmrc(i) = dq(i,k+1,m) + qpert(i,m)*fac1
797 : ! hplin - qpert for m>1 is always zero
798 0 : cmrc(i) = dq(i,k+1,m)
799 : else
800 0 : cmrc(i) = dq(i,k+1,m)
801 : end if
802 :
803 : ! Determine fluxes, flux divergence => changes due to convection
804 : ! Logic must be included to avoid producing negative values. A bit
805 : ! messy since there are no a priori assumptions about profiles.
806 : ! Tendency is modified (reduced) when pending disaster detected.
807 :
808 0 : botflx = etagdt(i)*(cmrc(i) - cmrh(i,k+1))*adjfac
809 0 : topflx = beta(i)*etagdt(i)*(cmrc(i)-cmrh(i,k))*adjfac
810 0 : dcmr1(i) = -botflx*rpd(i,k+1)
811 0 : efac1 = 1.0_kind_phys
812 0 : efac2 = 1.0_kind_phys
813 0 : efac3 = 1.0_kind_phys
814 :
815 0 : if (dq(i,k+1,m)+dcmr1(i) < 0.0_kind_phys) then
816 0 : if ( abs(dcmr1(i)) > 1.e-300_kind_phys ) then
817 0 : efac1 = max(tiny,abs(dq(i,k+1,m)/dcmr1(i)) - eps)
818 : else
819 0 : efac1 = tiny
820 : endif
821 : end if
822 :
823 0 : if (efac1 == tiny .or. efac1 > 1.0_kind_phys) efac1 = 0.0_kind_phys
824 0 : dcmr1(i) = -efac1*botflx*rpd(i,k+1)
825 0 : dcmr2(i) = (efac1*botflx - topflx)*rpd(i,k)
826 :
827 0 : if (dq(i,k,m)+dcmr2(i) < 0.0_kind_phys) then
828 0 : if ( abs(dcmr2(i)) > 1.e-300_kind_phys ) then
829 0 : efac2 = max(tiny,abs(dq(i,k ,m)/dcmr2(i)) - eps)
830 : else
831 0 : efac2 = tiny
832 : endif
833 : end if
834 :
835 0 : if (efac2 == tiny .or. efac2 > 1.0_kind_phys) efac2 = 0.0_kind_phys
836 0 : dcmr2(i) = (efac1*botflx - efac2*topflx)*rpd(i,k)
837 0 : dcmr3(i) = efac2*topflx*rpd(i,k-1)
838 :
839 0 : if ( (dq(i,k-1,m)+dcmr3(i) < 0.0_kind_phys ) ) then
840 0 : if ( abs(dcmr3(i)) > 1.e-300_kind_phys ) then
841 0 : efac3 = max(tiny,abs(dq(i,k-1,m)/dcmr3(i)) - eps)
842 : else
843 0 : efac3 = tiny
844 : endif
845 : end if
846 :
847 0 : if (efac3 == tiny .or. efac3 > 1.0_kind_phys) efac3 = 0.0_kind_phys
848 0 : efac3 = min(efac2,efac3)
849 0 : dcmr2(i) = (efac1*botflx - efac3*topflx)*rpd(i,k)
850 0 : dcmr3(i) = efac3*topflx*rpd(i,k-1)
851 :
852 0 : dq(i,k+1,m) = dq(i,k+1,m) + dcmr1(i)
853 0 : dq(i,k ,m) = dq(i,k ,m) + dcmr2(i)
854 0 : dq(i,k-1,m) = dq(i,k-1,m) + dcmr3(i)
855 : end do pcl1loop
856 : end do const_modify_loop
857 : ! Constituent modifications complete
858 :
859 : ! This if restructured from a goto
860 0 : if (k /= limcnv+1) then
861 : ! Complete update of thermodynamic structure at integer levels
862 : ! gather/scatter points that need new values of shbs and gamma
863 0 : do ii=1,len1
864 0 : i = indx1(ii)
865 0 : vtemp1(ii ) = tb(i,k)
866 0 : vtemp1(ii+len1) = tb(i,k-1)
867 0 : vtemp2(ii ) = pmid(i,k)
868 0 : vtemp2(ii+len1) = pmid(i,k-1)
869 : end do
870 0 : call qsat(vtemp1(1:2*len1), vtemp2(1:2*len1), &
871 0 : vtemp5(1:2*len1), vtemp3(1:2*len1), 2*len1, gam=vtemp4(1:2*len1))
872 0 : do ii=1,len1
873 0 : i = indx1(ii)
874 0 : shbs(i,k ) = vtemp3(ii )
875 0 : shbs(i,k-1) = vtemp3(ii+len1)
876 0 : gam(i,k ) = vtemp4(ii )
877 0 : gam(i,k-1) = vtemp4(ii+len1)
878 0 : sb (i,k ) = sb(i,k ) + ds2(i)
879 0 : sb (i,k-1) = sb(i,k-1) + ds3(i)
880 0 : hb (i,k ) = sb(i,k ) + hlat*shb(i,k )
881 0 : hb (i,k-1) = sb(i,k-1) + hlat*shb(i,k-1)
882 0 : hbs(i,k ) = sb(i,k ) + hlat*shbs(i,k )
883 0 : hbs(i,k-1) = sb(i,k-1) + hlat*shbs(i,k-1)
884 : end do
885 :
886 : ! Update thermodynamic information at half (i.e., interface) levels
887 0 : do ii=1,len1
888 0 : i = indx1(ii)
889 0 : sbh (i,k) = 0.5_kind_phys*(sb(i,k) + sb(i,k-1))
890 0 : shbh(i,k) = qhalf(shb(i,k-1),shb(i,k),shbs(i,k-1),shbs(i,k))
891 0 : hbh (i,k) = sbh(i,k) + hlat*shbh(i,k)
892 0 : sbh (i,k-1) = 0.5_kind_phys*(sb(i,k-1) + sb(i,k-2))
893 0 : shbh(i,k-1) = qhalf(shb(i,k-2),shb(i,k-1),shbs(i,k-2),shbs(i,k-1))
894 0 : hbh (i,k-1) = sbh(i,k-1) + hlat*shbh(i,k-1)
895 : end do
896 : end if ! k /= limcnv+1
897 :
898 : ! Ensure that dzcld is reset if convective mass flux zero
899 : ! specify the current vertical extent of the convective activity
900 : ! top of convective layer determined by size of overshoot param.
901 0 : do i=1,ncol
902 0 : etagt0 = eta(i).gt.0.0_kind_phys
903 0 : if ( .not. etagt0) dzcld(i) = 0.0_kind_phys
904 0 : if (etagt0 .and. beta(i) > betamn) then
905 0 : ktp = k-1
906 : else
907 : ktp = k
908 : end if
909 0 : if (etagt0) then
910 0 : cnt_sh(i) = min(cnt_sh(i),ktp)
911 0 : cnb_sh(i) = max(cnb_sh(i),k)
912 : end if
913 : end do
914 : end do kloop
915 :
916 : !---------------------------------------------------
917 : ! apply final thermodynamic tendencies
918 : !---------------------------------------------------
919 : ! Set output q tendencies...
920 : ! ...for water vapor
921 0 : dq(:ncol,:,const_wv_idx) = cmfdq(:ncol,:)
922 :
923 : ! ...for other tracers from passive tracer transport
924 0 : do m = 1, pcnst
925 0 : if (m .ne. const_wv_idx) then
926 0 : dq(:ncol,:,m) = (dq(:ncol,:,m) - q(:ncol,:,m))/ztodt
927 : endif
928 : enddo
929 :
930 : ! Kludge to prevent cnb_sh-cnt_sh from being zero (in the event
931 : ! someone decides that they want to divide by this quantity)
932 0 : do i=1,ncol
933 0 : if (cnb_sh(i) /= 0 .and. cnb_sh(i) == cnt_sh(i)) then
934 0 : cnt_sh(i) = cnt_sh(i) - 1
935 : end if
936 : end do
937 :
938 0 : do i=1,ncol
939 0 : precc(i) = prec(i)*rtdt
940 : end do
941 :
942 : ! Compute reserved liquid (not yet in cldliq) for energy integrals.
943 : ! Treat rliq_sh as flux out bottom, to be added back later.
944 0 : do k = 1, pver
945 0 : do i = 1, ncol
946 0 : rliq_sh(i) = rliq_sh(i) + qc_sh(i,k)*pdel(i,k)/grav
947 : end do
948 : end do
949 :
950 : ! rliq_sh is converted to precipitation units [m s-1]
951 0 : rliq_sh(:ncol) = rliq_sh(:ncol) / 1000._kind_phys
952 :
953 : ! Prepare boundary fluxes for check_energy [m s-1]
954 0 : flx_cnd(:ncol) = precc(:ncol) + rliq_sh(:ncol)
955 :
956 0 : if(debug_verbose) then
957 : ! DEBUG DIAGNOSTICS - show final result
958 0 : do i=1,ncol
959 0 : if (i == 1) then
960 0 : fac = grav*864._kind_phys
961 0 : write(iulog, 8010)
962 0 : do k=limcnv,pver
963 0 : rh = shb(i,k)/shbs(i,k)
964 0 : write(iulog, 8020) shbh(i,k),sbh(i,k),hbh(i,k),fac*cmfmc_sh(i,k), &
965 0 : cmfsl(i,k), cmflq(i,k)
966 0 : write(iulog, 8040) tb(i,k),shb(i,k),rh ,sb(i,k),hb(i,k), &
967 0 : hbs(i,k), ztodt*cmfdt(i,k),ztodt*cmfdq(i,k), &
968 0 : ztodt*cmfdqr(i,k)
969 : end do
970 0 : write(iulog, 8000) prec(i)
971 :
972 : ! approximate vertical integral of moist static energy and
973 : ! total preciptable water after adjustment and output changes
974 0 : hsum2 = 0.0_kind_phys
975 0 : qsum2 = 0.0_kind_phys
976 0 : do k=limcnv,pver
977 0 : hsum2 = hsum2 + pdel(i,k)*rgrav*hb(i,k)
978 0 : qsum2 = qsum2 + pdel(i,k)*rgrav*shb(i,k)
979 : end do
980 0 : write(iulog,8070) hsum1, hsum2, abs(hsum2-hsum1)/hsum2, &
981 0 : qsum1, qsum2, abs(qsum2-qsum1)/qsum2
982 : end if
983 : enddo
984 : endif
985 :
986 : ! Diagnostic use format strings
987 : 8000 format(///,10x,'PREC = ',3pf12.6,/)
988 : 8010 format('1** TB SHB RH SB', &
989 : ' HB HBS CAH CAM PRECC ', &
990 : ' ETA FSL FLQ **', /)
991 : 8020 format(' ----- ', 9x,3p,f7.3,2x,2p, 9x,-3p,f7.3,2x, &
992 : f7.3, 37x, 0p,2x,f8.2, 0p,2x,f8.2,2x,f8.2, ' ----- ')
993 : 8030 format(' ----- ', 0p,82x,f8.2, 0p,2x,f8.2,2x,f8.2, &
994 : ' ----- ')
995 : 8040 format(' - - - ',f7.3,2x,3p,f7.3,2x,2p,f7.3,2x,-3p,f7.3,2x, &
996 : f7.3, 2x,f8.3,2x,0p,f7.3,3p,2x,f7.3,2x,f7.3,30x, &
997 : ' - - - ')
998 : 8050 format(' ----- ',110x,' ----- ')
999 : 8060 format('1 K =>', i4,/, &
1000 : ' TB SHB RH SB', &
1001 : ' HB HBS CAH CAM PREC ', &
1002 : ' ETA FSL FLQ', /)
1003 : 8070 format(' VERTICALLY INTEGRATED MOIST STATIC ENERGY BEFORE, AFTER', &
1004 : ' AND PERCENTAGE DIFFERENCE => ',1p,2e15.7,2x,2p,f7.3,/, &
1005 : ' VERTICALLY INTEGRATED MOISTURE BEFORE, AFTER', &
1006 : ' AND PERCENTAGE DIFFERENCE => ',1p,2e15.7,2x,2p,f7.3,/)
1007 : 8080 format(' BETA, ETA => ', 1p,2e12.3)
1008 : 8090 format (' k+1, sc, shc, hc => ', 1x, i2, 1p, 3e12.4)
1009 0 : end subroutine hack_convect_shallow_run
1010 :
1011 : ! qhalf computes the specific humidity at interface levels between two model layers (interpolate moisture)
1012 0 : pure function qhalf(sh1,sh2,shbs1,shbs2) result(qh)
1013 : real(kind_phys), intent(in) :: sh1 ! humidity of layer 1 [kg kg-1]
1014 : real(kind_phys), intent(in) :: sh2 ! humidity of layer 2 [kg kg-1]
1015 : real(kind_phys), intent(in) :: shbs1 ! saturation specific humidity of layer 1 [kg kg-1]
1016 : real(kind_phys), intent(in) :: shbs2 ! saturation specific humidity of layer 2 [kg kg-1]
1017 : real(kind_phys) :: qh
1018 0 : qh = min(max(sh1,sh2),(shbs2*sh1 + shbs1*sh2)/(shbs1+shbs2))
1019 0 : end function qhalf
1020 : end module hack_convect_shallow
|
