Line data Source code
1 : module zm_convr
2 :
3 : use ccpp_kinds, only: kind_phys
4 :
5 : implicit none
6 :
7 : save
8 : private ! Make default type private to the module
9 : !
10 : ! PUBLIC: interfaces
11 : !
12 : public zm_convr_init ! ZM schemea
13 : public zm_convr_run ! ZM schemea
14 :
15 : real(kind_phys) rl ! wg latent heat of vaporization.
16 : real(kind_phys) cpres ! specific heat at constant pressure in j/kg-degk.
17 : real(kind_phys) :: capelmt ! namelist configurable:
18 : ! threshold value for cape for deep convection.
19 : real(kind_phys) :: ke ! Tunable evaporation efficiency set from namelist input zmconv_ke
20 : real(kind_phys) :: ke_lnd
21 : real(kind_phys) :: c0_lnd ! set from namelist input zmconv_c0_lnd
22 : real(kind_phys) :: c0_ocn ! set from namelist input zmconv_c0_ocn
23 : integer :: num_cin ! set from namelist input zmconv_num_cin
24 : ! The number of negative buoyancy regions that are allowed
25 : ! before the convection top and CAPE calculations are completed.
26 : real(kind_phys) tau ! convective time scale
27 : real(kind_phys) :: tfreez
28 : real(kind_phys) :: eps1
29 : real(kind_phys) :: momcu
30 : real(kind_phys) :: momcd
31 :
32 : logical :: no_deep_pbl ! default = .false.
33 : ! no_deep_pbl = .true. eliminates deep convection entirely within PBL
34 :
35 :
36 : real(kind_phys) :: rgrav ! reciprocal of grav
37 : real(kind_phys) :: rgas ! gas constant for dry air
38 : real(kind_phys) :: grav ! = gravit
39 : real(kind_phys) :: cp ! = cpres = cpair
40 :
41 : integer limcnv ! top interface level limit for convection
42 :
43 : logical :: lparcel_pbl ! Switch to turn on mixing of parcel MSE air, and picking launch level to be the top of the PBL.
44 : real(kind_phys) :: parcel_hscale
45 :
46 : real(kind_phys) :: tiedke_add ! namelist configurable
47 : real(kind_phys) :: dmpdz_param ! namelist configurable
48 :
49 : contains
50 :
51 :
52 :
53 : !===============================================================================
54 : !> \section arg_table_zm_convr_init Argument Table
55 : !! \htmlinclude zm_convr_init.html
56 : !!
57 1536 : subroutine zm_convr_init(plev, plevp, cpair, epsilo, gravit, latvap, tmelt, rair, &
58 1536 : pref_edge, zmconv_c0_lnd, zmconv_c0_ocn, zmconv_ke, zmconv_ke_lnd, &
59 : zmconv_momcu, zmconv_momcd, zmconv_num_cin, &
60 : no_deep_pbl_in, zmconv_tiedke_add, &
61 : zmconv_capelmt, zmconv_dmpdz, zmconv_parcel_pbl, zmconv_parcel_hscale, zmconv_tau, &
62 1536 : masterproc, iulog, errmsg, errflg)
63 :
64 : integer, intent(in) :: plev
65 : integer, intent(in) :: plevp
66 :
67 : real(kind_phys), intent(in) :: cpair ! specific heat of dry air (J K-1 kg-1)
68 : real(kind_phys), intent(in) :: epsilo ! ratio of h2o to dry air molecular weights
69 : real(kind_phys), intent(in) :: gravit ! gravitational acceleration (m s-2)
70 : real(kind_phys), intent(in) :: latvap ! Latent heat of vaporization (J kg-1)
71 : real(kind_phys), intent(in) :: tmelt ! Freezing point of water (K)
72 : real(kind_phys), intent(in) :: rair ! Dry air gas constant (J K-1 kg-1)
73 : real(kind_phys), intent(in) :: pref_edge(:) ! reference pressures at interfaces
74 : integer, intent(in) :: zmconv_num_cin ! Number negative buoyancy regions that are allowed
75 : ! before the convection top and CAPE calculations are completed.
76 : real(kind_phys),intent(in) :: zmconv_c0_lnd
77 : real(kind_phys),intent(in) :: zmconv_c0_ocn
78 : real(kind_phys),intent(in) :: zmconv_ke
79 : real(kind_phys),intent(in) :: zmconv_ke_lnd
80 : real(kind_phys),intent(in) :: zmconv_momcu
81 : real(kind_phys),intent(in) :: zmconv_momcd
82 : logical, intent(in) :: no_deep_pbl_in ! no_deep_pbl = .true. eliminates ZM convection entirely within PBL
83 : real(kind_phys),intent(in) :: zmconv_tiedke_add
84 : real(kind_phys),intent(in) :: zmconv_capelmt
85 : real(kind_phys),intent(in) :: zmconv_dmpdz
86 : logical, intent(in) :: zmconv_parcel_pbl ! Should the parcel properties include PBL mixing?
87 : real(kind_phys),intent(in) :: zmconv_parcel_hscale ! Fraction of PBL over which to mix ZM parcel.
88 : real(kind_phys),intent(in) :: zmconv_tau
89 : logical, intent(in) :: masterproc
90 : integer, intent(in) :: iulog
91 : character(len=512), intent(out) :: errmsg
92 : integer, intent(out) :: errflg
93 :
94 : integer :: k
95 :
96 1536 : errmsg =''
97 1536 : errflg = 0
98 :
99 : ! Initialization of ZM constants
100 1536 : tfreez = tmelt
101 1536 : eps1 = epsilo
102 1536 : rl = latvap
103 1536 : cpres = cpair
104 1536 : rgrav = 1.0_kind_phys/gravit
105 1536 : rgas = rair
106 1536 : grav = gravit
107 1536 : cp = cpres
108 :
109 1536 : c0_lnd = zmconv_c0_lnd
110 1536 : c0_ocn = zmconv_c0_ocn
111 1536 : num_cin = zmconv_num_cin
112 1536 : ke = zmconv_ke
113 1536 : ke_lnd = zmconv_ke_lnd
114 1536 : momcu = zmconv_momcu
115 1536 : momcd = zmconv_momcd
116 :
117 1536 : tiedke_add = zmconv_tiedke_add
118 1536 : capelmt = zmconv_capelmt
119 1536 : dmpdz_param = zmconv_dmpdz
120 1536 : no_deep_pbl = no_deep_pbl_in
121 1536 : lparcel_pbl = zmconv_parcel_pbl
122 1536 : parcel_hscale = zmconv_parcel_hscale
123 :
124 1536 : tau = zmconv_tau
125 :
126 : !
127 : ! Limit deep convection to regions below 40 mb
128 : ! Note this calculation is repeated in the shallow convection interface
129 : !
130 1536 : limcnv = 0 ! null value to check against below
131 1536 : if (pref_edge(1) >= 4.e3_kind_phys) then
132 0 : limcnv = 1
133 : else
134 67584 : do k=1,plev
135 67584 : if (pref_edge(k) < 4.e3_kind_phys .and. pref_edge(k+1) >= 4.e3_kind_phys) then
136 1536 : limcnv = k
137 1536 : exit
138 : end if
139 : end do
140 1536 : if ( limcnv == 0 ) limcnv = plevp
141 : end if
142 :
143 1536 : if ( masterproc ) then
144 2 : write(iulog,*)'ZM_CONV_INIT: Deep convection will be capped at intfc ',limcnv, &
145 4 : ' which is ',pref_edge(limcnv),' pascals'
146 : endif
147 :
148 1536 : if (masterproc) write(iulog,*)'**** ZM: DILUTE Buoyancy Calculation ****'
149 :
150 1536 : end subroutine zm_convr_init
151 :
152 :
153 : !===============================================================================
154 : !> \section arg_table_zm_convr_run Argument Table
155 : !! \htmlinclude zm_convr_run.html
156 : !!
157 80640 : subroutine zm_convr_run( ncol ,pver , &
158 : pverp, gravit ,latice ,cpwv ,cpliq , rh2o, &
159 80640 : lat, long, &
160 80640 : t ,qh ,prec , &
161 161280 : pblh ,zm ,geos ,zi ,qtnd , &
162 161280 : heat ,pap ,paph ,dpp , &
163 161280 : delt ,mcon ,cme ,cape , &
164 241920 : tpert ,dlf ,dif ,zdu ,rprd , &
165 80640 : mu ,md ,du ,eu ,ed , &
166 322560 : dp ,dsubcld ,jt ,maxg ,ideep , &
167 241920 : ql ,rliq ,landfrac, &
168 161280 : rice ,lengath ,scheme_name, errmsg ,errflg)
169 : !-----------------------------------------------------------------------
170 : !
171 : ! Purpose:
172 : ! Main driver for zhang-mcfarlane convection scheme
173 : !
174 : ! Method:
175 : ! performs deep convective adjustment based on mass-flux closure
176 : ! algorithm.
177 : !
178 : ! Author:guang jun zhang, m.lazare, n.mcfarlane. CAM Contact: P. Rasch
179 : !
180 : ! This is contributed code not fully standardized by the CAM core group.
181 : ! All variables have been typed, where most are identified in comments
182 : ! The current procedure will be reimplemented in a subsequent version
183 : ! of the CAM where it will include a more straightforward formulation
184 : ! and will make use of the standard CAM nomenclature
185 : !
186 : !-----------------------------------------------------------------------
187 : !
188 : ! ************************ index of variables **********************
189 : !
190 : ! wg * alpha array of vertical differencing used (=1. for upstream).
191 : ! w * cape convective available potential energy.
192 : ! wg * capeg gathered convective available potential energy.
193 : ! c * capelmt threshold value for cape for deep convection.
194 : ! ic * cpres specific heat at constant pressure in j/kg-degk.
195 : ! i * dpp
196 : ! ic * delt length of model time-step in seconds.
197 : ! wg * dp layer thickness in mbs (between upper/lower interface).
198 : ! wg * dqdt mixing ratio tendency at gathered points.
199 : ! wg * dsdt dry static energy ("temp") tendency at gathered points.
200 : ! wg * dudt u-wind tendency at gathered points.
201 : ! wg * dvdt v-wind tendency at gathered points.
202 : ! wg * dsubcld layer thickness in mbs between lcl and maxi.
203 : ! ic * grav acceleration due to gravity in m/sec2.
204 : ! wg * du detrainment in updraft. specified in mid-layer
205 : ! wg * ed entrainment in downdraft.
206 : ! wg * eu entrainment in updraft.
207 : ! wg * hmn moist static energy.
208 : ! wg * hsat saturated moist static energy.
209 : ! w * ideep holds position of gathered points vs longitude index.
210 : ! ic * pver number of model levels.
211 : ! wg * j0 detrainment initiation level index.
212 : ! wg * jd downdraft initiation level index.
213 : ! ic * jlatpr gaussian latitude index for printing grids (if needed).
214 : ! wg * jt top level index of deep cumulus convection.
215 : ! w * lcl base level index of deep cumulus convection.
216 : ! wg * lclg gathered values of lcl.
217 : ! w * lel index of highest theoretical convective plume.
218 : ! wg * lelg gathered values of lel.
219 : ! w * lon index of onset level for deep convection.
220 : ! w * maxi index of level with largest moist static energy.
221 : ! wg * maxg gathered values of maxi.
222 : ! wg * mb cloud base mass flux.
223 : ! wg * mc net upward (scaled by mb) cloud mass flux.
224 : ! wg * md downward cloud mass flux (positive up).
225 : ! wg * mu upward cloud mass flux (positive up). specified
226 : ! at interface
227 : ! ic * msg number of missing moisture levels at the top of model.
228 : ! w * p grid slice of ambient mid-layer pressure in mbs.
229 : ! i * pblt row of pbl top indices.
230 : ! w * pcpdh scaled surface pressure.
231 : ! w * pf grid slice of ambient interface pressure in mbs.
232 : ! wg * pg grid slice of gathered values of p.
233 : ! w * q grid slice of mixing ratio.
234 : ! wg * qd grid slice of mixing ratio in downdraft.
235 : ! wg * qg grid slice of gathered values of q.
236 : ! i/o * qh grid slice of specific humidity.
237 : ! w * qh0 grid slice of initial specific humidity.
238 : ! wg * qhat grid slice of upper interface mixing ratio.
239 : ! wg * ql grid slice of cloud liquid water.
240 : ! wg * qs grid slice of saturation mixing ratio.
241 : ! w * qstp grid slice of parcel temp. saturation mixing ratio.
242 : ! wg * qstpg grid slice of gathered values of qstp.
243 : ! wg * qu grid slice of mixing ratio in updraft.
244 : ! ic * rgas dry air gas constant.
245 : ! wg * rl latent heat of vaporization.
246 : ! w * s grid slice of scaled dry static energy (t+gz/cp).
247 : ! wg * sd grid slice of dry static energy in downdraft.
248 : ! wg * sg grid slice of gathered values of s.
249 : ! wg * shat grid slice of upper interface dry static energy.
250 : ! wg * su grid slice of dry static energy in updraft.
251 : ! i/o * t
252 : ! wg * tg grid slice of gathered values of t.
253 : ! w * tl row of parcel temperature at lcl.
254 : ! wg * tlg grid slice of gathered values of tl.
255 : ! w * tp grid slice of parcel temperatures.
256 : ! wg * tpg grid slice of gathered values of tp.
257 : ! i/o * u grid slice of u-wind (real).
258 : ! wg * ug grid slice of gathered values of u.
259 : ! i/o * utg grid slice of u-wind tendency (real).
260 : ! i/o * v grid slice of v-wind (real).
261 : ! w * va work array re-used by called subroutines.
262 : ! wg * vg grid slice of gathered values of v.
263 : ! i/o * vtg grid slice of v-wind tendency (real).
264 : ! i * w grid slice of diagnosed large-scale vertical velocity.
265 : ! w * z grid slice of ambient mid-layer height in metres.
266 : ! w * zf grid slice of ambient interface height in metres.
267 : ! wg * zfg grid slice of gathered values of zf.
268 : ! wg * zg grid slice of gathered values of z.
269 : !
270 : !-----------------------------------------------------------------------
271 : !
272 : ! multi-level i/o fields:
273 : ! i => input arrays.
274 : ! i/o => input/output arrays.
275 : ! w => work arrays.
276 : ! wg => work arrays operating only on gathered points.
277 : ! ic => input data constants.
278 : ! c => data constants pertaining to subroutine itself.
279 : !
280 : ! input arguments
281 : !
282 : integer, intent(in) :: ncol ! number of atmospheric columns
283 : integer, intent(in) :: pver, pverp
284 :
285 : real(kind_phys), intent(in) :: gravit ! gravitational acceleration (m s-2)
286 : real(kind_phys), intent(in) :: latice ! Latent heat of fusion (J kg-1)
287 : real(kind_phys), intent(in) :: cpwv ! specific heat of water vapor (J K-1 kg-1)
288 : real(kind_phys), intent(in) :: cpliq ! specific heat of fresh h2o (J K-1 kg-1)
289 : real(kind_phys), intent(in) :: rh2o ! Water vapor gas constant (J K-1 kg-1)
290 :
291 : real(kind_phys), intent(in) :: lat(:)
292 : real(kind_phys), intent(in) :: long(:)
293 :
294 : real(kind_phys), intent(in) :: t(:,:) ! grid slice of temperature at mid-layer. (ncol,pver)
295 : real(kind_phys), intent(in) :: qh(:,:) ! grid slice of specific humidity. (ncol,pver)
296 : real(kind_phys), intent(in) :: pap(:,:) ! (ncol,pver)
297 : real(kind_phys), intent(in) :: paph(:,:) ! (ncol,pver+1)
298 : real(kind_phys), intent(in) :: dpp(:,:) ! local sigma half-level thickness (i.e. dshj). (ncol,pver)
299 : real(kind_phys), intent(in) :: zm(:,:) ! (ncol,pver)
300 : real(kind_phys), intent(in) :: geos(:) ! (ncol)
301 : real(kind_phys), intent(in) :: zi(:,:) ! (ncol,pver+1)
302 : real(kind_phys), intent(in) :: pblh(:) ! (ncol)
303 : real(kind_phys), intent(in) :: tpert(:) ! (ncol)
304 : real(kind_phys), intent(in) :: landfrac(:) ! RBN Landfrac (ncol)
305 :
306 : ! output arguments
307 : !
308 : real(kind_phys), intent(out) :: qtnd(:,:) ! specific humidity tendency (kg/kg/s) (ncol,pver)
309 : real(kind_phys), intent(out) :: heat(:,:) ! heating rate (dry static energy tendency, W/kg) (ncol,pver)
310 : real(kind_phys), intent(out) :: mcon(:,:) ! (ncol,pverp)
311 : real(kind_phys), intent(out) :: dif(:,:)
312 : real(kind_phys), intent(out) :: dlf(:,:) ! scattrd version of the detraining cld h2o tend (ncol,pver)
313 : real(kind_phys), intent(out) :: cme(:,:) ! (ncol,pver)
314 : real(kind_phys), intent(out) :: cape(:) ! w convective available potential energy. (ncol)
315 : real(kind_phys), intent(out) :: zdu(:,:) ! (ncol,pver)
316 : real(kind_phys), intent(out) :: rprd(:,:) ! rain production rate (ncol,pver)
317 :
318 : ! move these vars from local storage to output so that convective
319 : ! transports can be done in outside of conv_cam.
320 : real(kind_phys), intent(out) :: mu(:,:) ! (ncol,pver)
321 : real(kind_phys), intent(out) :: eu(:,:) ! (ncol,pver)
322 : real(kind_phys), intent(out) :: du(:,:) ! (ncol,pver)
323 : real(kind_phys), intent(out) :: md(:,:) ! (ncol,pver)
324 : real(kind_phys), intent(out) :: ed(:,:) ! (ncol,pver)
325 : real(kind_phys), intent(out) :: dp(:,:) ! wg layer thickness in mbs (between upper/lower interface). (ncol,pver)
326 : real(kind_phys), intent(out) :: dsubcld(:) ! wg layer thickness in mbs between lcl and maxi. (ncol)
327 : real(kind_phys), intent(out) :: prec(:) ! (ncol)
328 : real(kind_phys), intent(out) :: rliq(:) ! reserved liquid (not yet in cldliq) for energy integrals (ncol)
329 : real(kind_phys), intent(out) :: rice(:) ! reserved ice (not yet in cldce) for energy integrals (ncol)
330 :
331 : integer, intent(out) :: ideep(:) ! column indices of gathered points (ncol)
332 :
333 : integer, intent(out) :: jt(:) ! wg top level index of deep cumulus convection.
334 : integer, intent(out) :: maxg(:)! wg gathered values of maxi.
335 :
336 : integer, intent(out) :: lengath
337 :
338 : real(kind_phys),intent(out):: ql(:,:) ! wg grid slice of cloud liquid water.
339 :
340 : character(len=40), intent(out) :: scheme_name
341 : character(len=512), intent(out) :: errmsg
342 : integer, intent(out) :: errflg
343 :
344 :
345 : ! Local variables
346 :
347 :
348 161280 : real(kind_phys) zs(ncol)
349 161280 : real(kind_phys) dlg(ncol,pver) ! gathrd version of the detraining cld h2o tend
350 161280 : real(kind_phys) cug(ncol,pver) ! gathered condensation rate
351 :
352 161280 : real(kind_phys) evpg(ncol,pver) ! gathered evap rate of rain in downdraft
353 : real(kind_phys) dptot(ncol)
354 :
355 161280 : real(kind_phys) mumax(ncol)
356 :
357 : !
358 161280 : real(kind_phys) pblt(ncol) ! i row of pbl top indices.
359 :
360 :
361 :
362 :
363 : !
364 : !-----------------------------------------------------------------------
365 : !
366 : ! general work fields (local variables):
367 : !
368 161280 : real(kind_phys) q(ncol,pver) ! w grid slice of mixing ratio.
369 161280 : real(kind_phys) p(ncol,pver) ! w grid slice of ambient mid-layer pressure in mbs.
370 161280 : real(kind_phys) z(ncol,pver) ! w grid slice of ambient mid-layer height in metres.
371 161280 : real(kind_phys) s(ncol,pver) ! w grid slice of scaled dry static energy (t+gz/cp).
372 161280 : real(kind_phys) tp(ncol,pver) ! w grid slice of parcel temperatures.
373 161280 : real(kind_phys) zf(ncol,pver+1) ! w grid slice of ambient interface height in metres.
374 161280 : real(kind_phys) pf(ncol,pver+1) ! w grid slice of ambient interface pressure in mbs.
375 161280 : real(kind_phys) qstp(ncol,pver) ! w grid slice of parcel temp. saturation mixing ratio.
376 :
377 161280 : real(kind_phys) tl(ncol) ! w row of parcel temperature at lcl.
378 :
379 161280 : integer lcl(ncol) ! w base level index of deep cumulus convection.
380 161280 : integer lel(ncol) ! w index of highest theoretical convective plume.
381 161280 : integer lon(ncol) ! w index of onset level for deep convection.
382 161280 : integer maxi(ncol) ! w index of level with largest moist static energy.
383 :
384 : real(kind_phys) precip
385 : !
386 : ! gathered work fields:
387 : !
388 161280 : real(kind_phys) qg(ncol,pver) ! wg grid slice of gathered values of q.
389 161280 : real(kind_phys) tg(ncol,pver) ! w grid slice of temperature at interface.
390 161280 : real(kind_phys) pg(ncol,pver) ! wg grid slice of gathered values of p.
391 161280 : real(kind_phys) zg(ncol,pver) ! wg grid slice of gathered values of z.
392 161280 : real(kind_phys) sg(ncol,pver) ! wg grid slice of gathered values of s.
393 161280 : real(kind_phys) tpg(ncol,pver) ! wg grid slice of gathered values of tp.
394 161280 : real(kind_phys) zfg(ncol,pver+1) ! wg grid slice of gathered values of zf.
395 161280 : real(kind_phys) qstpg(ncol,pver) ! wg grid slice of gathered values of qstp.
396 161280 : real(kind_phys) ug(ncol,pver) ! wg grid slice of gathered values of u.
397 161280 : real(kind_phys) vg(ncol,pver) ! wg grid slice of gathered values of v.
398 161280 : real(kind_phys) cmeg(ncol,pver)
399 :
400 161280 : real(kind_phys) rprdg(ncol,pver) ! wg gathered rain production rate
401 161280 : real(kind_phys) capeg(ncol) ! wg gathered convective available potential energy.
402 161280 : real(kind_phys) tlg(ncol) ! wg grid slice of gathered values of tl.
403 161280 : real(kind_phys) landfracg(ncol) ! wg grid slice of landfrac
404 :
405 161280 : integer lclg(ncol) ! wg gathered values of lcl.
406 161280 : integer lelg(ncol)
407 : !
408 : ! work fields arising from gathered calculations.
409 : !
410 161280 : real(kind_phys) dqdt(ncol,pver) ! wg mixing ratio tendency at gathered points.
411 161280 : real(kind_phys) dsdt(ncol,pver) ! wg dry static energy ("temp") tendency at gathered points.
412 161280 : real(kind_phys) sd(ncol,pver) ! wg grid slice of dry static energy in downdraft.
413 161280 : real(kind_phys) qd(ncol,pver) ! wg grid slice of mixing ratio in downdraft.
414 161280 : real(kind_phys) mc(ncol,pver) ! wg net upward (scaled by mb) cloud mass flux.
415 161280 : real(kind_phys) qhat(ncol,pver) ! wg grid slice of upper interface mixing ratio.
416 161280 : real(kind_phys) qu(ncol,pver) ! wg grid slice of mixing ratio in updraft.
417 161280 : real(kind_phys) su(ncol,pver) ! wg grid slice of dry static energy in updraft.
418 161280 : real(kind_phys) qs(ncol,pver) ! wg grid slice of saturation mixing ratio.
419 161280 : real(kind_phys) shat(ncol,pver) ! wg grid slice of upper interface dry static energy.
420 161280 : real(kind_phys) hmn(ncol,pver) ! wg moist static energy.
421 161280 : real(kind_phys) hsat(ncol,pver) ! wg saturated moist static energy.
422 161280 : real(kind_phys) qlg(ncol,pver)
423 161280 : real(kind_phys) dudt(ncol,pver) ! wg u-wind tendency at gathered points.
424 161280 : real(kind_phys) dvdt(ncol,pver) ! wg v-wind tendency at gathered points.
425 :
426 161280 : real(kind_phys) qldeg(ncol,pver) ! cloud liquid water mixing ratio for detrainment (kg/kg)
427 161280 : real(kind_phys) mb(ncol) ! wg cloud base mass flux.
428 :
429 161280 : integer jlcl(ncol)
430 161280 : integer j0(ncol) ! wg detrainment initiation level index.
431 80640 : integer jd(ncol) ! wg downdraft initiation level index.
432 :
433 : real(kind_phys),intent(in):: delt ! length of model time-step in seconds.
434 :
435 : integer i
436 : integer ii
437 : integer k, kk, l, m
438 :
439 : integer msg ! ic number of missing moisture levels at the top of model.
440 : real(kind_phys) qdifr
441 : real(kind_phys) sdifr
442 :
443 : real(kind_phys), parameter :: dcon = 25.e-6_kind_phys
444 : real(kind_phys), parameter :: mucon = 5.3_kind_phys
445 : real(kind_phys) negadq
446 : logical doliq
447 :
448 :
449 : !
450 : !--------------------------Data statements------------------------------
451 :
452 80640 : scheme_name = "zm_convr_run"
453 80640 : errmsg = ''
454 80640 : errflg = 0
455 : !
456 : ! Set internal variable "msg" (convection limit) to "limcnv-1"
457 : !
458 80640 : msg = limcnv - 1
459 : !
460 : ! initialize necessary arrays.
461 : ! zero out variables not used in cam
462 : !
463 :
464 :
465 87010560 : qtnd(:,:) = 0._kind_phys
466 87010560 : heat(:,:) = 0._kind_phys
467 88252416 : mcon(:,:) = 0._kind_phys
468 1241856 : rliq(:ncol) = 0._kind_phys
469 1241856 : rice(:ncol) = 0._kind_phys
470 :
471 : !
472 : ! initialize convective tendencies
473 : !
474 1241856 : prec(:ncol) = 0._kind_phys
475 5725440 : do k = 1,pver
476 87010560 : do i = 1,ncol
477 81285120 : dqdt(i,k) = 0._kind_phys
478 81285120 : dsdt(i,k) = 0._kind_phys
479 81285120 : dudt(i,k) = 0._kind_phys
480 81285120 : dvdt(i,k) = 0._kind_phys
481 81285120 : cme(i,k) = 0._kind_phys
482 81285120 : rprd(i,k) = 0._kind_phys
483 81285120 : zdu(i,k) = 0._kind_phys
484 81285120 : ql(i,k) = 0._kind_phys
485 81285120 : qlg(i,k) = 0._kind_phys
486 81285120 : dlf(i,k) = 0._kind_phys
487 81285120 : dlg(i,k) = 0._kind_phys
488 81285120 : qldeg(i,k) = 0._kind_phys
489 :
490 86929920 : dif(i,k) = 0._kind_phys
491 :
492 : end do
493 : end do
494 :
495 1241856 : do i = 1,ncol
496 1161216 : pblt(i) = pver
497 1241856 : dsubcld(i) = 0._kind_phys
498 :
499 : end do
500 :
501 :
502 : !
503 : ! calculate local pressure (mbs) and height (m) for both interface
504 : ! and mid-layer locations.
505 : !
506 1241856 : do i = 1,ncol
507 1161216 : zs(i) = geos(i)*rgrav
508 1161216 : pf(i,pver+1) = paph(i,pver+1)*0.01_kind_phys
509 1241856 : zf(i,pver+1) = zi(i,pver+1) + zs(i)
510 : end do
511 5725440 : do k = 1,pver
512 87010560 : do i = 1,ncol
513 81285120 : p(i,k) = pap(i,k)*0.01_kind_phys
514 81285120 : pf(i,k) = paph(i,k)*0.01_kind_phys
515 81285120 : z(i,k) = zm(i,k) + zs(i)
516 86929920 : zf(i,k) = zi(i,k) + zs(i)
517 : end do
518 : end do
519 : !
520 2177280 : do k = pver - 1,msg + 1,-1
521 32368896 : do i = 1,ncol
522 32288256 : if (abs(z(i,k)-zs(i)-pblh(i)) < (zf(i,k)-zf(i,k+1))*0.5_kind_phys) pblt(i) = k
523 : end do
524 : end do
525 : !
526 : ! store incoming specific humidity field for subsequent calculation
527 : ! of precipitation (through change in storage).
528 : ! define dry static energy (normalized by cp).
529 : !
530 5725440 : do k = 1,pver
531 87010560 : do i = 1,ncol
532 81285120 : q(i,k) = qh(i,k)
533 81285120 : s(i,k) = t(i,k) + (grav/cpres)*z(i,k)
534 81285120 : tp(i,k)=0.0_kind_phys
535 81285120 : shat(i,k) = s(i,k)
536 86929920 : qhat(i,k) = q(i,k)
537 : end do
538 : end do
539 :
540 1241856 : do i = 1,ncol
541 1161216 : capeg(i) = 0._kind_phys
542 1161216 : lclg(i) = 1
543 1161216 : lelg(i) = pver
544 1161216 : maxg(i) = 1
545 1161216 : tlg(i) = 400._kind_phys
546 1241856 : dsubcld(i) = 0._kind_phys
547 : end do
548 :
549 :
550 : ! Evaluate Tparcel, qs(Tparcel), buoyancy and CAPE,
551 : ! lcl, lel, parcel launch level at index maxi()=hmax
552 :
553 : call buoyan_dilute( ncol ,pver , &
554 : cpliq ,latice ,cpwv ,rh2o ,&
555 : q ,t ,p ,z ,pf , &
556 : tp ,qstp ,tl ,rl ,cape , &
557 : pblt ,lcl ,lel ,lon ,maxi , &
558 : rgas ,grav ,cpres ,msg , &
559 : zi ,zs ,tpert , landfrac,&
560 175182336 : lat ,long ,errmsg ,errflg)
561 :
562 : !
563 : ! determine whether grid points will undergo some deep convection
564 : ! (ideep=1) or not (ideep=0), based on values of cape,lcl,lel
565 : ! (require cape.gt. 0 and lel<lcl as minimum conditions).
566 : !
567 80640 : lengath = 0
568 1241856 : ideep = 0
569 1241856 : do i=1,ncol
570 1241856 : if (cape(i) > capelmt) then
571 277765 : lengath = lengath + 1
572 277765 : ideep(lengath) = i
573 : end if
574 : end do
575 :
576 80640 : if (lengath.eq.0) return
577 : !
578 : ! obtain gathered arrays necessary for ensuing calculations.
579 : !
580 5681207 : do k = 1,pver
581 25124757 : do i = 1,lengath
582 19443550 : dp(i,k) = 0.01_kind_phys*dpp(ideep(i),k)
583 19443550 : qg(i,k) = q(ideep(i),k)
584 19443550 : tg(i,k) = t(ideep(i),k)
585 19443550 : pg(i,k) = p(ideep(i),k)
586 19443550 : zg(i,k) = z(ideep(i),k)
587 19443550 : sg(i,k) = s(ideep(i),k)
588 19443550 : tpg(i,k) = tp(ideep(i),k)
589 19443550 : zfg(i,k) = zf(ideep(i),k)
590 19443550 : qstpg(i,k) = qstp(ideep(i),k)
591 19443550 : ug(i,k) = 0._kind_phys
592 25044740 : vg(i,k) = 0._kind_phys
593 : end do
594 : end do
595 :
596 : !
597 357782 : do i = 1,lengath
598 357782 : zfg(i,pver+1) = zf(ideep(i),pver+1)
599 : end do
600 357782 : do i = 1,lengath
601 277765 : capeg(i) = cape(ideep(i))
602 277765 : lclg(i) = lcl(ideep(i))
603 277765 : lelg(i) = lel(ideep(i))
604 277765 : maxg(i) = maxi(ideep(i))
605 277765 : tlg(i) = tl(ideep(i))
606 357782 : landfracg(i) = landfrac(ideep(i))
607 : end do
608 : !
609 : ! calculate sub-cloud layer pressure "thickness" for use in
610 : ! closure and tendency routines.
611 : !
612 2240476 : do k = msg + 1,pver
613 9740131 : do i = 1,lengath
614 9660114 : if (k >= maxg(i)) then
615 298552 : dsubcld(i) = dsubcld(i) + dp(i,k)
616 : end if
617 : end do
618 : end do
619 : !
620 : ! define array of factors (alpha) which defines interfacial
621 : ! values, as well as interfacial values for (q,s) used in
622 : ! subsequent routines.
623 : !
624 2160459 : do k = msg + 2,pver
625 9382349 : do i = 1,lengath
626 7221890 : sdifr = 0._kind_phys
627 7221890 : qdifr = 0._kind_phys
628 7221890 : if (sg(i,k) > 0._kind_phys .or. sg(i,k-1) > 0._kind_phys) &
629 7221890 : sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k)))
630 7221890 : if (qg(i,k) > 0._kind_phys .or. qg(i,k-1) > 0._kind_phys) &
631 7221890 : qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k)))
632 7221890 : if (sdifr > 1.E-6_kind_phys) then
633 7219572 : shat(i,k) = log(sg(i,k-1)/sg(i,k))*sg(i,k-1)*sg(i,k)/(sg(i,k-1)-sg(i,k))
634 : else
635 2318 : shat(i,k) = 0.5_kind_phys* (sg(i,k)+sg(i,k-1))
636 : end if
637 9302332 : if (qdifr > 1.E-6_kind_phys) then
638 7221865 : qhat(i,k) = log(qg(i,k-1)/qg(i,k))*qg(i,k-1)*qg(i,k)/(qg(i,k-1)-qg(i,k))
639 : else
640 25 : qhat(i,k) = 0.5_kind_phys* (qg(i,k)+qg(i,k-1))
641 : end if
642 : end do
643 : end do
644 : !
645 : ! obtain cloud properties.
646 : !
647 :
648 : call cldprp(ncol ,pver ,pverp ,cpliq , &
649 : latice ,cpwv ,rh2o ,&
650 : qg ,tg ,ug ,vg ,pg , &
651 : zg ,sg ,mu ,eu ,du , &
652 : md ,ed ,sd ,qd ,mc , &
653 : qu ,su ,zfg ,qs ,hmn , &
654 : hsat ,shat ,qlg , &
655 : cmeg ,maxg ,lelg ,jt ,jlcl , &
656 : maxg ,j0 ,jd ,rl ,lengath , &
657 : rgas ,grav ,cpres ,msg , &
658 : evpg ,cug ,rprdg ,limcnv ,landfracg , &
659 431397617 : qldeg ,qhat )
660 :
661 :
662 : !
663 : ! convert detrainment from units of "1/m" to "1/mb".
664 : !
665 :
666 2240476 : do k = msg + 1,pver
667 9740131 : do i = 1,lengath
668 7499655 : du (i,k) = du (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
669 7499655 : eu (i,k) = eu (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
670 7499655 : ed (i,k) = ed (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
671 7499655 : cug (i,k) = cug (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
672 7499655 : cmeg (i,k) = cmeg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
673 7499655 : rprdg(i,k) = rprdg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
674 9660114 : evpg (i,k) = evpg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k)
675 : end do
676 : end do
677 :
678 : call closure(ncol ,pver , &
679 : qg ,tg ,pg ,zg ,sg , &
680 : tpg ,qs ,qu ,su ,mc , &
681 : du ,mu ,md ,qd ,sd , &
682 : qhat ,shat ,dp ,qstpg ,zfg , &
683 : qlg ,dsubcld ,mb ,capeg ,tlg , &
684 : lclg ,lelg ,jt ,maxg ,1 , &
685 : lengath ,rgas ,grav ,cpres ,rl , &
686 345134097 : msg ,capelmt )
687 : !
688 : ! limit cloud base mass flux to theoretical upper bound.
689 : !
690 357782 : do i=1,lengath
691 357782 : mumax(i) = 0
692 : end do
693 2160459 : do k=msg + 2,pver
694 9382349 : do i=1,lengath
695 9302332 : mumax(i) = max(mumax(i), mu(i,k)/dp(i,k))
696 : end do
697 : end do
698 :
699 357782 : do i=1,lengath
700 357782 : if (mumax(i) > 0._kind_phys) then
701 276413 : mb(i) = min(mb(i),1._kind_phys/(delt*mumax(i)))
702 : else
703 1352 : mb(i) = 0._kind_phys
704 : endif
705 : end do
706 : ! If no_deep_pbl = .true., don't allow convection entirely
707 : ! within PBL (suggestion of Bjorn Stevens, 8-2000)
708 :
709 80017 : if (no_deep_pbl) then
710 0 : do i=1,lengath
711 0 : if (zm(ideep(i),jt(i)) < pblh(ideep(i))) mb(i) = 0
712 : end do
713 : end if
714 :
715 2240476 : do k=msg+1,pver
716 9740131 : do i=1,lengath
717 7499655 : mu (i,k) = mu (i,k)*mb(i)
718 7499655 : md (i,k) = md (i,k)*mb(i)
719 7499655 : mc (i,k) = mc (i,k)*mb(i)
720 7499655 : du (i,k) = du (i,k)*mb(i)
721 7499655 : eu (i,k) = eu (i,k)*mb(i)
722 7499655 : ed (i,k) = ed (i,k)*mb(i)
723 7499655 : cmeg (i,k) = cmeg (i,k)*mb(i)
724 7499655 : rprdg(i,k) = rprdg(i,k)*mb(i)
725 7499655 : cug (i,k) = cug (i,k)*mb(i)
726 9660114 : evpg (i,k) = evpg (i,k)*mb(i)
727 :
728 : end do
729 : end do
730 : !
731 : ! compute temperature and moisture changes due to convection.
732 : !
733 : call q1q2_pjr(ncol ,pver ,latice , &
734 : dqdt ,dsdt ,qg ,qs ,qu , &
735 : su ,du ,qhat ,shat ,dp , &
736 : mu ,md ,sd ,qd ,qldeg , &
737 : dsubcld ,jt ,maxg ,1 ,lengath , &
738 : cpres ,rl ,msg , &
739 345134097 : dlg ,evpg ,cug)
740 :
741 : !
742 : ! gather back temperature and mixing ratio.
743 : !
744 :
745 2240476 : do k = msg + 1,pver
746 9740131 : do i = 1,lengath
747 : !
748 : ! q is updated to compute net precip.
749 : !
750 7499655 : q(ideep(i),k) = qh(ideep(i),k) + delt*dqdt(i,k)
751 7499655 : qtnd(ideep(i),k) = dqdt (i,k)
752 7499655 : cme (ideep(i),k) = cmeg (i,k)
753 7499655 : rprd(ideep(i),k) = rprdg(i,k)
754 7499655 : zdu (ideep(i),k) = du (i,k)
755 7499655 : mcon(ideep(i),k) = mc (i,k)
756 7499655 : heat(ideep(i),k) = dsdt (i,k)*cpres
757 7499655 : dlf (ideep(i),k) = dlg (i,k)
758 9660114 : ql (ideep(i),k) = qlg (i,k)
759 : end do
760 : end do
761 :
762 : ! Compute precip by integrating change in water vapor minus detrained cloud water
763 2240476 : do k = pver,msg + 1,-1
764 33353089 : do i = 1,ncol
765 33273072 : prec(i) = prec(i) - dpp(i,k)* (q(i,k)-qh(i,k)) - dpp(i,k)*(dlf(i,k)+dif(i,k))*delt
766 : end do
767 : end do
768 :
769 : ! obtain final precipitation rate in m/s.
770 1232336 : do i = 1,ncol
771 1232336 : prec(i) = rgrav*max(prec(i),0._kind_phys)/ delt/1000._kind_phys
772 : end do
773 :
774 : ! Compute reserved liquid (not yet in cldliq) for energy integrals.
775 : ! Treat rliq as flux out bottom, to be added back later.
776 5681207 : do k = 1, pver
777 86343537 : do i = 1, ncol
778 80662330 : rliq(i) = rliq(i) + (dlf(i,k)+dif(i,k))*dpp(i,k)/gravit
779 86263520 : rice(i) = rice(i) + dif(i,k)*dpp(i,k)/gravit
780 : end do
781 : end do
782 1232336 : rliq(:ncol) = rliq(:ncol) /1000._kind_phys ! Converted to precip units m s-1
783 1232336 : rice(:ncol) = rice(:ncol) /1000._kind_phys ! Converted to precip units m s-1
784 :
785 : ! Convert mass flux from reported mb s-1 to kg m-2 s-1
786 87575873 : mcon(:ncol,:pverp) = mcon(:ncol,:pverp) * 100._kind_phys / gravit
787 :
788 : return
789 : end subroutine zm_convr_run
790 :
791 : subroutine zm_convr_finalize
792 : end subroutine zm_convr_finalize
793 :
794 : !=========================================================================================
795 :
796 80640 : subroutine buoyan_dilute( ncol ,pver , &
797 : cpliq ,latice ,cpwv ,rh2o ,&
798 80640 : q ,t ,p ,z ,pf , &
799 80640 : tp ,qstp ,tl ,rl ,cape , &
800 80640 : pblt ,lcl ,lel ,lon ,mx , &
801 : rd ,grav ,cp ,msg , &
802 80640 : zi ,zs ,tpert , landfrac,&
803 80640 : lat ,long ,errmsg ,errflg)
804 : !-----------------------------------------------------------------------
805 : !
806 : ! Purpose:
807 : ! Calculates CAPE the lifting condensation level and the convective top
808 : ! where buoyancy is first -ve.
809 : !
810 : ! Method: Calculates the parcel temperature based on a simple constant
811 : ! entraining plume model. CAPE is integrated from buoyancy.
812 : ! 09/09/04 - Simplest approach using an assumed entrainment rate for
813 : ! testing (dmpdp).
814 : ! 08/04/05 - Swap to convert dmpdz to dmpdp
815 : !
816 : ! SCAM Logical Switches - DILUTE:RBN - Now Disabled
817 : ! ---------------------
818 : ! switch(1) = .T. - Uses the dilute parcel calculation to obtain tendencies.
819 : ! switch(2) = .T. - Includes entropy/q changes due to condensate loss and freezing.
820 : ! switch(3) = .T. - Adds the PBL Tpert for the parcel temperature at all levels.
821 : !
822 : ! References:
823 : ! Raymond and Blythe (1992) JAS
824 : !
825 : ! Author:
826 : ! Richard Neale - September 2004
827 : !
828 : !-----------------------------------------------------------------------
829 : implicit none
830 : !-----------------------------------------------------------------------
831 : !
832 : ! input arguments
833 : !
834 : integer, intent(in) :: ncol ! number of atmospheric columns
835 : integer, intent(in) :: pver
836 : real(kind_phys), intent(in) :: cpliq
837 : real(kind_phys), intent(in) :: latice
838 : real(kind_phys), intent(in) :: cpwv
839 : real(kind_phys), intent(in) :: rh2o
840 :
841 : real(kind_phys), intent(in) :: q(ncol,pver) ! spec. humidity
842 : real(kind_phys), intent(in) :: t(ncol,pver) ! temperature
843 : real(kind_phys), intent(in) :: p(ncol,pver) ! pressure
844 : real(kind_phys), intent(in) :: z(ncol,pver) ! height
845 : real(kind_phys), intent(in) :: pf(ncol,pver+1) ! pressure at interfaces
846 : real(kind_phys), intent(in) :: pblt(ncol) ! index of pbl depth
847 : real(kind_phys), intent(in) :: tpert(ncol) ! perturbation temperature by pbl processes
848 :
849 : ! Use z interface/surface relative values for PBL parcel calculations.
850 : real(kind_phys), intent(in) :: zi(ncol,pver+1)
851 : real(kind_phys), intent(in) :: zs(ncol)
852 :
853 : real(kind_phys), intent(in) :: lat(:)
854 : real(kind_phys), intent(in) :: long(:)
855 :
856 : !
857 : ! output arguments
858 : !
859 :
860 : real(kind_phys), intent(out) :: tp(ncol,pver) ! parcel temperature
861 : real(kind_phys), intent(out) :: qstp(ncol,pver) ! saturation mixing ratio of parcel (only above lcl, just q below).
862 : real(kind_phys), intent(out) :: tl(ncol) ! parcel temperature at lcl
863 : real(kind_phys), intent(out) :: cape(ncol) ! convective aval. pot. energy.
864 : integer lcl(ncol) !
865 : integer lel(ncol) !
866 : integer lon(ncol) ! level of onset of deep convection
867 : integer mx(ncol) ! level of max moist static energy
868 :
869 : real(kind_phys), intent(in) :: landfrac(ncol)
870 : character(len=512), intent(out) :: errmsg
871 : integer, intent(out) :: errflg
872 :
873 : !
874 : !--------------------------Local Variables------------------------------
875 : !
876 161280 : real(kind_phys) capeten(ncol,5) ! provisional value of cape
877 161280 : real(kind_phys) tv(ncol,pver) !
878 161280 : real(kind_phys) tpv(ncol,pver) !
879 161280 : real(kind_phys) buoy(ncol,pver)
880 :
881 : real(kind_phys) a1(ncol)
882 : real(kind_phys) a2(ncol)
883 : real(kind_phys) estp(ncol)
884 161280 : real(kind_phys) pl(ncol)
885 : real(kind_phys) plexp(ncol)
886 161280 : real(kind_phys) hmax(ncol)
887 161280 : real(kind_phys) hmn(ncol)
888 : real(kind_phys) y(ncol)
889 :
890 161280 : logical plge600(ncol)
891 161280 : integer knt(ncol)
892 161280 : integer lelten(ncol,5)
893 :
894 :
895 :
896 :
897 : ! Parcel property variables
898 :
899 161280 : real(kind_phys) :: hmn_lev(ncol,pver) ! Vertical profile of moist static energy for each column
900 161280 : real(kind_phys) :: dp_lev(ncol,pver) ! Level dpressure between interfaces
901 161280 : real(kind_phys) :: hmn_zdp(ncol,pver) ! Integrals of hmn_lev*dp_lev at each level
902 161280 : real(kind_phys) :: q_zdp(ncol,pver) ! Integrals of q*dp_lev at each level
903 : real(kind_phys) :: dp_zfrac ! Fraction of vertical grid box below mixing top (usually pblt)
904 161280 : real(kind_phys) :: parcel_dz(ncol) ! Depth of parcel mixing (usually parcel_hscale*parcel_dz)
905 161280 : real(kind_phys) :: parcel_ztop(ncol) ! Height of parcel mixing (usually parcel_ztop+zm(nlev))
906 161280 : real(kind_phys) :: parcel_dp(ncol) ! Pressure integral over parcel mixing depth (usually pblt)
907 161280 : real(kind_phys) :: parcel_hdp(ncol) ! Pressure*MSE integral over parcel mixing depth (usually pblt)
908 161280 : real(kind_phys) :: parcel_qdp(ncol) ! Pressure*q integral over parcel mixing depth (usually pblt)
909 161280 : real(kind_phys) :: pbl_dz(ncol) ! Previously diagnosed PBL height
910 161280 : real(kind_phys) :: hpar(ncol) ! Initial MSE of the parcel
911 161280 : real(kind_phys) :: qpar(ncol) ! Initial humidity of the parcel
912 80640 : real(kind_phys) :: ql(ncol) ! Initial parcel humidity (for ientropy routine)
913 : integer :: ipar ! Index for top of parcel mixing/launch level.
914 :
915 :
916 :
917 :
918 : real(kind_phys) cp
919 : real(kind_phys) e
920 : real(kind_phys) grav
921 :
922 : integer i
923 : integer k
924 : integer msg
925 : integer n
926 :
927 : real(kind_phys) rd
928 : real(kind_phys) rl
929 :
930 : !
931 : !-----------------------------------------------------------------------
932 : !
933 483840 : do n = 1,5
934 6289920 : do i = 1,ncol
935 5806080 : lelten(i,n) = pver
936 6209280 : capeten(i,n) = 0._kind_phys
937 : end do
938 : end do
939 : !
940 1241856 : do i = 1,ncol
941 1161216 : lon(i) = pver
942 1161216 : knt(i) = 0
943 1161216 : lel(i) = pver
944 1161216 : mx(i) = lon(i)
945 1161216 : cape(i) = 0._kind_phys
946 1161216 : hmax(i) = 0._kind_phys
947 1161216 : pbl_dz(i) = z(i,nint(pblt(i)))-zs(i) ! mid-point z (zm) reference to PBL depth
948 1161216 : parcel_dz(i) = max(zi(i,pver),parcel_hscale*pbl_dz(i)) ! PBL mixing depth [parcel_hscale*Boundary, but no thinner than zi(i,pver)]
949 1161216 : parcel_ztop(i) = parcel_dz(i)+zs(i) ! PBL mixing height ztop this is wrt zs=0
950 1161216 : parcel_hdp(i) = 0._kind_phys
951 1161216 : parcel_dp(i) = 0._kind_phys
952 1161216 : parcel_qdp(i) = 0._kind_phys
953 1161216 : hpar(i) = 0._kind_phys
954 1241856 : qpar(i) = 0._kind_phys
955 : end do
956 :
957 87010560 : tp(:ncol,:) = t(:ncol,:)
958 87010560 : qstp(:ncol,:) = q(:ncol,:)
959 87010560 : hmn_lev(:ncol,:) = 0._kind_phys
960 :
961 :
962 :
963 : !!! Initialize tv and buoy for output.
964 : !!! tv=tv : tpv=tpv : qstp=q : buoy=0.
965 87010560 : tv(:ncol,:) = t(:ncol,:) *(1._kind_phys+1.608_kind_phys*q(:ncol,:))/ (1._kind_phys+q(:ncol,:))
966 87010560 : tpv(:ncol,:) = tv(:ncol,:)
967 87010560 : buoy(:ncol,:) = 0._kind_phys
968 :
969 :
970 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
971 : ! Mix the parcel over a certain dp or dz and take the launch level as the top level
972 : ! of this mixing region and the parcel properties as this mixed value
973 : ! Should be well mixed by other processes in the very near PBL.
974 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
975 :
976 :
977 :
978 80640 : if (lparcel_pbl) then
979 :
980 : ! Vertical profile of MSE and pressure weighted of the same.
981 0 : hmn_lev(:ncol,1:pver) = cp*t(:ncol,1:pver) + grav*z(:ncol,1:pver) + rl*q(:ncol,1:pver)
982 0 : dp_lev(:ncol,1:pver) = pf(:ncol,2:pver+1)-pf(:ncol,1:pver)
983 0 : hmn_zdp(:ncol,1:pver) = hmn_lev(:ncol,1:pver)*dp_lev(:ncol,1:pver)
984 0 : q_zdp(:ncol,1:pver) = q(:ncol,1:pver)*dp_lev(:ncol,1:pver)
985 :
986 :
987 : ! Mix profile over vertical length scale of 0.5*PBLH.
988 :
989 0 : do i = 1,ncol ! Loop columns
990 0 : do k = pver,msg + 1,-1
991 :
992 0 : if (zi(i,k+1)<= parcel_dz(i)) then ! Has to be relative to near-surface layer center elevation
993 0 : ipar = k
994 :
995 0 : if (k == pver) then ! Always at least the full depth of lowest model layer.
996 : dp_zfrac = 1._kind_phys
997 : else
998 : ! Fraction of grid cell depth (mostly 1, except when parcel_ztop is in between levels.
999 0 : dp_zfrac = min(1._kind_phys,(parcel_dz(i)-zi(i,k+1))/(zi(i,k)-zi(i,k+1)))
1000 : end if
1001 :
1002 0 : parcel_hdp(i) = parcel_hdp(i)+hmn_zdp(i,k)*dp_zfrac ! Sum parcel profile up to a certain level.
1003 0 : parcel_qdp(i) = parcel_qdp(i)+q_zdp(i,k)*dp_zfrac ! Sum parcel profile up to a certain level.
1004 0 : parcel_dp(i) = parcel_dp(i)+dp_lev(i,k)*dp_zfrac ! SUM dp's for weighting of parcel_hdp
1005 :
1006 : end if
1007 : end do
1008 0 : hpar(i) = parcel_hdp(i)/parcel_dp(i)
1009 0 : qpar(i) = parcel_qdp(i)/parcel_dp(i)
1010 0 : mx(i) = ipar
1011 : end do
1012 :
1013 : else ! Default method finding level of MSE maximum (nlev sensitive though)
1014 : !
1015 : ! set "launching" level(mx) to be at maximum moist static energy.
1016 : ! search for this level stops at planetary boundary layer top.
1017 : !
1018 2257920 : do k = pver,msg + 1,-1
1019 33610752 : do i = 1,ncol
1020 31352832 : hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
1021 33530112 : if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then
1022 1757831 : hmax(i) = hmn(i)
1023 1757831 : mx(i) = k
1024 : end if
1025 : end do
1026 : end do
1027 :
1028 : end if ! Default method of determining parcel launch properties.
1029 :
1030 :
1031 :
1032 :
1033 :
1034 : ! LCL dilute calculation - initialize to mx(i)
1035 : ! Determine lcl in parcel_dilute and get pl,tl after parcel_dilute
1036 : ! Original code actually sets LCL as level above wher condensate forms.
1037 : ! Therefore in parcel_dilute lcl(i) will be at first level where qsmix < qtmix.
1038 :
1039 80640 : if (lparcel_pbl) then
1040 :
1041 : ! For parcel dilute need to invert hpar and qpar.
1042 : ! Now need to supply ql(i) as it is mixed parcel version, just q(i,max(i)) in default
1043 :
1044 0 : do i = 1,ncol ! Initialise LCL variables.
1045 0 : lcl(i) = mx(i)
1046 0 : tl(i) = (hpar(i)-rl*qpar(i)-grav*parcel_ztop(i))/cp
1047 0 : ql(i) = qpar(i)
1048 0 : pl(i) = p(i,mx(i))
1049 : end do
1050 :
1051 : else
1052 :
1053 1241856 : do i = 1,ncol
1054 1161216 : lcl(i) = mx(i)
1055 1161216 : tl(i) = t(i,mx(i))
1056 1161216 : ql(i) = q(i,mx(i))
1057 1241856 : pl(i) = p(i,mx(i))
1058 : end do
1059 :
1060 : end if ! Mixed parcel properties
1061 :
1062 :
1063 :
1064 : !
1065 : ! main buoyancy calculation.
1066 : !
1067 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1068 : !!! DILUTE PLUME CALCULATION USING ENTRAINING PLUME !!!
1069 : !!! RBN 9/9/04 !!!
1070 :
1071 : call parcel_dilute(ncol, pver, cpliq, cpwv, rh2o, latice, msg, mx, p, t, q, &
1072 : tpert, tp, tpv, qstp, pl, tl, ql, lcl, &
1073 80640 : landfrac, lat, long, errmsg, errflg)
1074 :
1075 :
1076 : ! If lcl is above the nominal level of non-divergence (600 mbs),
1077 : ! no deep convection is permitted (ensuing calculations
1078 : ! skipped and cape retains initialized value of zero).
1079 : !
1080 1241856 : do i = 1,ncol
1081 1241856 : plge600(i) = pl(i).ge.600._kind_phys ! Just change to always allow buoy calculation.
1082 : end do
1083 :
1084 : !
1085 : ! Main buoyancy calculation.
1086 : !
1087 2257920 : do k = pver,msg + 1,-1
1088 33610752 : do i=1,ncol
1089 33530112 : if (k <= mx(i) .and. plge600(i)) then ! Define buoy from launch level to cloud top.
1090 28583546 : tv(i,k) = t(i,k)* (1._kind_phys+1.608_kind_phys*q(i,k))/ (1._kind_phys+q(i,k))
1091 28583546 : buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add ! +0.5K or not?
1092 : else
1093 2769286 : qstp(i,k) = q(i,k)
1094 2769286 : tp(i,k) = t(i,k)
1095 2769286 : tpv(i,k) = tv(i,k)
1096 : endif
1097 : end do
1098 : end do
1099 :
1100 :
1101 :
1102 : !-------------------------------------------------------------------------------
1103 :
1104 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1105 :
1106 :
1107 :
1108 : !
1109 2177280 : do k = msg + 2,pver
1110 32368896 : do i = 1,ncol
1111 32288256 : if (k < lcl(i) .and. plge600(i)) then
1112 22770604 : if (buoy(i,k+1) > 0._kind_phys .and. buoy(i,k) <= 0._kind_phys) then
1113 920606 : knt(i) = min(num_cin,knt(i) + 1)
1114 920606 : lelten(i,knt(i)) = k
1115 : end if
1116 : end if
1117 : end do
1118 : end do
1119 : !
1120 : ! calculate convective available potential energy (cape).
1121 : !
1122 161280 : do n = 1,num_cin
1123 2338560 : do k = msg + 1,pver
1124 33610752 : do i = 1,ncol
1125 33530112 : if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then
1126 5641946 : capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k))
1127 : end if
1128 : end do
1129 : end do
1130 : end do
1131 : !
1132 : ! find maximum cape from all possible tentative capes from
1133 : ! one sounding,
1134 : ! and use it as the final cape, april 26, 1995
1135 : !
1136 161280 : do n = 1,num_cin
1137 1322496 : do i = 1,ncol
1138 1241856 : if (capeten(i,n) > cape(i)) then
1139 740506 : cape(i) = capeten(i,n)
1140 740506 : lel(i) = lelten(i,n)
1141 : end if
1142 : end do
1143 : end do
1144 : !
1145 : ! put lower bound on cape for diagnostic purposes.
1146 : !
1147 1241856 : do i = 1,ncol
1148 1241856 : cape(i) = max(cape(i), 0._kind_phys)
1149 : end do
1150 : !
1151 80640 : return
1152 : end subroutine buoyan_dilute
1153 :
1154 80640 : subroutine parcel_dilute (ncol, pver, cpliq, cpwv, rh2o, latice, msg, klaunch, p, t, q, &
1155 161280 : tpert, tp, tpv, qstp, pl, tl, ql, lcl, &
1156 80640 : landfrac,lat,long,errmsg,errflg)
1157 :
1158 : ! Routine to determine
1159 : ! 1. Tp - Parcel temperature
1160 : ! 2. qstp - Saturated mixing ratio at the parcel temperature.
1161 :
1162 : !--------------------
1163 : implicit none
1164 : !--------------------
1165 :
1166 : integer, intent(in) :: ncol
1167 : integer, intent(in) :: pver
1168 : real(kind_phys), intent(in) :: cpliq
1169 : real(kind_phys), intent(in) :: cpwv
1170 : real(kind_phys), intent(in) :: rh2o
1171 : real(kind_phys), intent(in) :: latice
1172 : integer, intent(in) :: msg
1173 :
1174 : integer, intent(in), dimension(ncol) :: klaunch(ncol)
1175 :
1176 : real(kind_phys), intent(in), dimension(ncol,pver) :: p
1177 : real(kind_phys), intent(in), dimension(ncol,pver) :: t
1178 : real(kind_phys), intent(in), dimension(ncol,pver) :: q
1179 : real(kind_phys), intent(in), dimension(ncol) :: tpert ! PBL temperature perturbation.
1180 :
1181 : real(kind_phys), intent(in) :: lat(:)
1182 : real(kind_phys), intent(in) :: long(:)
1183 :
1184 : real(kind_phys), intent(inout), dimension(ncol,pver) :: tp ! Parcel temp.
1185 : real(kind_phys), intent(inout), dimension(ncol,pver) :: qstp ! Parcel water vapour (sat value above lcl).
1186 : real(kind_phys), intent(inout), dimension(ncol) :: tl ! Actual temp of LCL.
1187 : real(kind_phys), intent(inout), dimension(ncol) :: ql ! Actual humidity of LCL
1188 : real(kind_phys), intent(inout), dimension(ncol) :: pl ! Actual pressure of LCL.
1189 :
1190 : integer, intent(inout), dimension(ncol) :: lcl ! Lifting condesation level (first model level with saturation).
1191 :
1192 : real(kind_phys), intent(out), dimension(ncol,pver) :: tpv ! Define tpv within this routine.
1193 :
1194 : character(len=512), intent(out) :: errmsg
1195 : integer, intent(out) :: errflg
1196 :
1197 :
1198 :
1199 : real(kind_phys), intent(in), dimension(ncol) :: landfrac
1200 : !--------------------
1201 :
1202 : ! Have to be careful as s is also dry static energy.
1203 :
1204 :
1205 : ! If we are to retain the fact that CAM loops over grid-points in the internal
1206 : ! loop then we need to dimension sp,atp,mp,xsh2o with ncol.
1207 :
1208 :
1209 161280 : real(kind_phys) tmix(ncol,pver) ! Tempertaure of the entraining parcel.
1210 161280 : real(kind_phys) qtmix(ncol,pver) ! Total water of the entraining parcel.
1211 161280 : real(kind_phys) qsmix(ncol,pver) ! Saturated mixing ratio at the tmix.
1212 161280 : real(kind_phys) smix(ncol,pver) ! Entropy of the entraining parcel.
1213 161280 : real(kind_phys) xsh2o(ncol,pver) ! Precipitate lost from parcel.
1214 161280 : real(kind_phys) ds_xsh2o(ncol,pver) ! Entropy change due to loss of condensate.
1215 161280 : real(kind_phys) ds_freeze(ncol,pver) ! Entropy change sue to freezing of precip.
1216 : real(kind_phys) dmpdz2d(ncol,pver) ! variable detrainment rate
1217 :
1218 161280 : real(kind_phys) mp(ncol) ! Parcel mass flux.
1219 161280 : real(kind_phys) qtp(ncol) ! Parcel total water.
1220 161280 : real(kind_phys) sp(ncol) ! Parcel entropy.
1221 :
1222 161280 : real(kind_phys) sp0(ncol) ! Parcel launch entropy.
1223 161280 : real(kind_phys) qtp0(ncol) ! Parcel launch total water.
1224 80640 : real(kind_phys) mp0(ncol) ! Parcel launch relative mass flux.
1225 :
1226 : real(kind_phys) lwmax ! Maximum condesate that can be held in cloud before rainout.
1227 : real(kind_phys) dmpdp ! Parcel fractional mass entrainment rate (/mb).
1228 : real(kind_phys) dmpdz ! Parcel fractional mass entrainment rate (/m)
1229 : real(kind_phys) dpdz,dzdp ! Hydrstatic relation and inverse of.
1230 : real(kind_phys) senv ! Environmental entropy at each grid point.
1231 : real(kind_phys) qtenv ! Environmental total water " " ".
1232 : real(kind_phys) penv ! Environmental total pressure " " ".
1233 : real(kind_phys) tenv ! Environmental total temperature " " ".
1234 : real(kind_phys) new_s ! Hold value for entropy after condensation/freezing adjustments.
1235 : real(kind_phys) new_q ! Hold value for total water after condensation/freezing adjustments.
1236 : real(kind_phys) dp ! Layer thickness (center to center)
1237 : real(kind_phys) tfguess ! First guess for entropy inversion - crucial for efficiency!
1238 : real(kind_phys) tscool ! Super cooled temperature offset (in degC) (eg -35).
1239 :
1240 : real(kind_phys) qxsk, qxskp1 ! LCL excess water (k, k+1)
1241 : real(kind_phys) dsdp, dqtdp, dqxsdp ! LCL s, qt, p gradients (k, k+1)
1242 : real(kind_phys) slcl,qtlcl,qslcl ! LCL s, qt, qs values.
1243 : real(kind_phys) dmpdz_lnd, dmpdz_mask
1244 :
1245 : integer rcall ! Number of ientropy call for errors recording
1246 : integer nit_lheat ! Number of iterations for condensation/freezing loop.
1247 : integer i,k,ii ! Loop counters.
1248 :
1249 : !======================================================================
1250 : ! SUMMARY
1251 : !
1252 : ! 9/9/04 - Assumes parcel is initiated from level of maxh (klaunch)
1253 : ! and entrains at each level with a specified entrainment rate.
1254 : !
1255 : ! 15/9/04 - Calculates lcl(i) based on k where qsmix is first < qtmix.
1256 : !
1257 : !======================================================================
1258 : !
1259 : ! Set some values that may be changed frequently.
1260 : !
1261 :
1262 80640 : nit_lheat = 2 ! iterations for ds,dq changes from condensation freezing.
1263 80640 : dmpdz=dmpdz_param ! Entrainment rate. (-ve for /m)
1264 80640 : dmpdz_lnd=-1.e-3_kind_phys
1265 80640 : lwmax = 1.e-3_kind_phys ! Need to put formula in for this.
1266 80640 : tscool = 0.0_kind_phys ! Temp at which water loading freezes in the cloud.
1267 :
1268 87010560 : qtmix=0._kind_phys
1269 87010560 : smix=0._kind_phys
1270 :
1271 80640 : qtenv = 0._kind_phys
1272 80640 : senv = 0._kind_phys
1273 80640 : tenv = 0._kind_phys
1274 80640 : penv = 0._kind_phys
1275 :
1276 1241856 : qtp0 = 0._kind_phys
1277 1241856 : sp0 = 0._kind_phys
1278 1241856 : mp0 = 0._kind_phys
1279 :
1280 1241856 : qtp = 0._kind_phys
1281 1241856 : sp = 0._kind_phys
1282 1241856 : mp = 0._kind_phys
1283 :
1284 80640 : new_q = 0._kind_phys
1285 80640 : new_s = 0._kind_phys
1286 :
1287 : ! **** Begin loops ****
1288 :
1289 2257920 : do k = pver, msg+1, -1
1290 33610752 : do i=1,ncol
1291 :
1292 : ! Initialize parcel values at launch level.
1293 :
1294 31352832 : if (k == klaunch(i)) then
1295 :
1296 1161216 : if (lparcel_pbl) then ! Modifcations to parcel properties if lparcel_pbl set.
1297 :
1298 0 : qtp0(i) = ql(i) ! Parcel launch q (PBL mixed value).
1299 0 : sp0(i) = entropy(tl(i),pl(i),qtp0(i),cpliq,cpwv,rh2o) ! Parcel launch entropy could be a mixed parcel.
1300 :
1301 : else
1302 :
1303 1161216 : qtp0(i) = q(i,k) ! Parcel launch total water (assuming subsaturated)
1304 1161216 : sp0(i) = entropy(t(i,k),p(i,k),qtp0(i),cpliq,cpwv,rh2o) ! Parcel launch entropy.
1305 :
1306 : end if
1307 :
1308 1161216 : mp0(i) = 1._kind_phys ! Parcel launch relative mass (i.e. 1 parcel stays 1 parcel for dmpdp=0, undilute).
1309 1161216 : smix(i,k) = sp0(i)
1310 1161216 : qtmix(i,k) = qtp0(i)
1311 1161216 : tfguess = t(i,k)
1312 1161216 : rcall = 1
1313 : call ientropy (rcall,i,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess,cpliq,cpwv,rh2o,&
1314 1161216 : lat(i), long(i), errmsg,errflg)
1315 : end if
1316 :
1317 : ! Entraining levels
1318 :
1319 33530112 : if (k < klaunch(i)) then
1320 :
1321 : ! Set environmental values for this level.
1322 :
1323 29482978 : dp = (p(i,k)-p(i,k+1)) ! In -ve mb as p decreasing with height - difference between center of layers.
1324 29482978 : qtenv = 0.5_kind_phys*(q(i,k)+q(i,k+1)) ! Total water of environment.
1325 29482978 : tenv = 0.5_kind_phys*(t(i,k)+t(i,k+1))
1326 29482978 : penv = 0.5_kind_phys*(p(i,k)+p(i,k+1))
1327 :
1328 29482978 : senv = entropy(tenv,penv,qtenv,cpliq,cpwv,rh2o) ! Entropy of environment.
1329 :
1330 : ! Determine fractional entrainment rate /pa given value /m.
1331 :
1332 29482978 : dpdz = -(penv*grav)/(rgas*tenv) ! in mb/m since p in mb.
1333 29482978 : dzdp = 1._kind_phys/dpdz ! in m/mb
1334 29482978 : dmpdp = dmpdz*dzdp
1335 :
1336 : ! Sum entrainment to current level
1337 : ! entrains q,s out of intervening dp layers, in which linear variation is assumed
1338 : ! so really it entrains the mean of the 2 stored values.
1339 :
1340 29482978 : sp(i) = sp(i) - dmpdp*dp*senv
1341 29482978 : qtp(i) = qtp(i) - dmpdp*dp*qtenv
1342 29482978 : mp(i) = mp(i) - dmpdp*dp
1343 :
1344 : ! Entrain s and qt to next level.
1345 :
1346 29482978 : smix(i,k) = (sp0(i) + sp(i)) / (mp0(i) + mp(i))
1347 29482978 : qtmix(i,k) = (qtp0(i) + qtp(i)) / (mp0(i) + mp(i))
1348 :
1349 : ! Invert entropy from s and q to determine T and saturation-capped q of mixture.
1350 : ! t(i,k) used as a first guess so that it converges faster.
1351 :
1352 29482978 : tfguess = tmix(i,k+1)
1353 29482978 : rcall = 2
1354 29482978 : call ientropy(rcall,i,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess,cpliq,cpwv,rh2o,lat(i),&
1355 29482978 : long(i),errmsg,errflg)
1356 :
1357 : !
1358 : ! Determine if this is lcl of this column if qsmix <= qtmix.
1359 : ! FIRST LEVEL where this happens on ascending.
1360 :
1361 29482978 : if (qsmix(i,k) <= qtmix(i,k) .and. qsmix(i,k+1) > qtmix(i,k+1)) then
1362 1093419 : lcl(i) = k
1363 1093419 : qxsk = qtmix(i,k) - qsmix(i,k)
1364 1093419 : qxskp1 = qtmix(i,k+1) - qsmix(i,k+1)
1365 1093419 : dqxsdp = (qxsk - qxskp1)/dp
1366 1093419 : pl(i) = p(i,k+1) - qxskp1/dqxsdp ! pressure level of actual lcl.
1367 1093419 : dsdp = (smix(i,k) - smix(i,k+1))/dp
1368 1093419 : dqtdp = (qtmix(i,k) - qtmix(i,k+1))/dp
1369 1093419 : slcl = smix(i,k+1) + dsdp* (pl(i)-p(i,k+1))
1370 1093419 : qtlcl = qtmix(i,k+1) + dqtdp*(pl(i)-p(i,k+1))
1371 :
1372 1093419 : tfguess = tmix(i,k)
1373 1093419 : rcall = 3
1374 1093419 : call ientropy (rcall,i,slcl,pl(i),qtlcl,tl(i),qslcl,tfguess,cpliq,cpwv,rh2o,lat(i), long(i), errmsg,errflg)
1375 :
1376 : endif
1377 : !
1378 : end if ! k < klaunch
1379 :
1380 :
1381 : end do ! Levels loop
1382 : end do ! Columns loop
1383 :
1384 : !!!!!!!!!!!!!!!!!!!!!!!!!!END ENTRAINMENT LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1385 :
1386 : !! Could stop now and test with this as it will provide some estimate of buoyancy
1387 : !! without the effects of freezing/condensation taken into account for tmix.
1388 :
1389 : !! So we now have a profile of entropy and total water of the entraining parcel
1390 : !! Varying with height from the launch level klaunch parcel=environment. To the
1391 : !! top allowed level for the existence of convection.
1392 :
1393 : !! Now we have to adjust these values such that the water held in vaopor is < or
1394 : !! = to qsmix. Therefore, we assume that the cloud holds a certain amount of
1395 : !! condensate (lwmax) and the rest is rained out (xsh2o). This, obviously
1396 : !! provides latent heating to the mixed parcel and so this has to be added back
1397 : !! to it. But does this also increase qsmix as well? Also freezing processes
1398 :
1399 :
1400 87010560 : xsh2o = 0._kind_phys
1401 87010560 : ds_xsh2o = 0._kind_phys
1402 87010560 : ds_freeze = 0._kind_phys
1403 :
1404 : !!!!!!!!!!!!!!!!!!!!!!!!!PRECIPITATION/FREEZING LOOP!!!!!!!!!!!!!!!!!!!!!!!!!!
1405 : !! Iterate solution twice for accuracy
1406 :
1407 :
1408 :
1409 2257920 : do k = pver, msg+1, -1
1410 33610752 : do i=1,ncol
1411 :
1412 : ! Initialize variables at k=klaunch
1413 :
1414 31352832 : if (k == klaunch(i)) then
1415 :
1416 : ! Set parcel values at launch level assume no liquid water.
1417 :
1418 1161216 : tp(i,k) = tmix(i,k)
1419 1161216 : qstp(i,k) = q(i,k)
1420 1161216 : tpv(i,k) = (tp(i,k) + tpert(i)) * (1._kind_phys+1.608_kind_phys*qstp(i,k)) / (1._kind_phys+qstp(i,k))
1421 :
1422 : end if
1423 :
1424 33530112 : if (k < klaunch(i)) then
1425 :
1426 : ! Initiaite loop if switch(2) = .T. - RBN:DILUTE - TAKEN OUT BUT COULD BE RETURNED LATER.
1427 :
1428 : ! Iterate nit_lheat times for s,qt changes.
1429 :
1430 88448934 : do ii=0,nit_lheat-1
1431 :
1432 : ! Rain (xsh2o) is excess condensate, bar LWMAX (Accumulated loss from qtmix).
1433 :
1434 58965956 : xsh2o(i,k) = max (0._kind_phys, qtmix(i,k) - qsmix(i,k) - lwmax)
1435 :
1436 : ! Contribution to ds from precip loss of condensate (Accumulated change from smix).(-ve)
1437 :
1438 58965956 : ds_xsh2o(i,k) = ds_xsh2o(i,k+1) - cpliq * log (tmix(i,k)/tfreez) * max(0._kind_phys,(xsh2o(i,k)-xsh2o(i,k+1)))
1439 : !
1440 : ! Entropy of freezing: latice times amount of water involved divided by T.
1441 : !
1442 :
1443 58965956 : if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) == 0._kind_phys) then ! One off freezing of condensate.
1444 5165052 : ds_freeze(i,k) = (latice/tmix(i,k)) * max(0._kind_phys,qtmix(i,k)-qsmix(i,k)-xsh2o(i,k)) ! Gain of LH
1445 : end if
1446 :
1447 58965956 : if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) /= 0._kind_phys) then ! Continual freezing of additional condensate.
1448 44017916 : ds_freeze(i,k) = ds_freeze(i,k+1)+(latice/tmix(i,k)) * max(0._kind_phys,(qsmix(i,k+1)-qsmix(i,k)))
1449 : end if
1450 :
1451 : ! Adjust entropy and accordingly to sum of ds (be careful of signs).
1452 :
1453 58965956 : new_s = smix(i,k) + ds_xsh2o(i,k) + ds_freeze(i,k)
1454 :
1455 : ! Adjust liquid water and accordingly to xsh2o.
1456 :
1457 58965956 : new_q = qtmix(i,k) - xsh2o(i,k)
1458 :
1459 : ! Invert entropy to get updated Tmix and qsmix of parcel.
1460 :
1461 58965956 : tfguess = tmix(i,k)
1462 58965956 : rcall =4
1463 : call ientropy (rcall,i,new_s, p(i,k), new_q, tmix(i,k), qsmix(i,k), tfguess,cpliq,cpwv,rh2o,&
1464 88448934 : lat(i), long(i), errmsg,errflg)
1465 :
1466 : end do ! Iteration loop for freezing processes.
1467 :
1468 : ! tp - Parcel temp is temp of mixture.
1469 : ! tpv - Parcel v. temp should be density temp with new_q total water.
1470 :
1471 29482978 : tp(i,k) = tmix(i,k)
1472 :
1473 : ! tpv = tprho in the presence of condensate (i.e. when new_q > qsmix)
1474 :
1475 29482978 : if (new_q > qsmix(i,k)) then ! Super-saturated so condensate present - reduces buoyancy.
1476 26356582 : qstp(i,k) = qsmix(i,k)
1477 : else ! Just saturated/sub-saturated - no condensate virtual effects.
1478 3126396 : qstp(i,k) = new_q
1479 : end if
1480 :
1481 29482978 : tpv(i,k) = (tp(i,k)+tpert(i))* (1._kind_phys+1.608_kind_phys*qstp(i,k)) / (1._kind_phys+ new_q)
1482 :
1483 : end if ! k < klaunch
1484 :
1485 : end do ! Loop for columns
1486 :
1487 : end do ! Loop for vertical levels.
1488 :
1489 :
1490 80640 : return
1491 : end subroutine parcel_dilute
1492 :
1493 : !-----------------------------------------------------------------------------------------
1494 566885567 : real(kind_phys) function entropy(TK,p,qtot,cpliq,cpwv,rh2o)
1495 : !-----------------------------------------------------------------------------------------
1496 : !
1497 : ! TK(K),p(mb),qtot(kg/kg)
1498 : ! from Raymond and Blyth 1992
1499 : !
1500 : real(kind_phys), intent(in) :: p,qtot,TK
1501 : real(kind_phys), intent(in) :: cpliq
1502 : real(kind_phys), intent(in) :: cpwv
1503 : real(kind_phys), intent(in) :: rh2o
1504 :
1505 : real(kind_phys) :: qv,qst,e,est,L
1506 : real(kind_phys), parameter :: pref = 1000._kind_phys
1507 :
1508 566885567 : L = rl - (cpliq - cpwv)*(TK-tfreez) ! T IN CENTIGRADE
1509 :
1510 566885567 : call qsat_hPa(TK, p, est, qst)
1511 :
1512 566885567 : qv = min(qtot,qst) ! Partition qtot into vapor part only.
1513 566885567 : e = qv*p / (eps1 +qv)
1514 :
1515 : entropy = (cpres + qtot*cpliq)*log( TK/tfreez) - rgas*log( (p-e)/pref ) + &
1516 566885567 : L*qv/TK - qv*rh2o*log(qv/qst)
1517 :
1518 566885567 : end FUNCTION entropy
1519 :
1520 : !
1521 : !-----------------------------------------------------------------------------------------
1522 90703569 : SUBROUTINE ientropy (rcall,icol,s,p,qt,T,qst,Tfg,cpliq,cpwv,rh2o,this_lat,this_lon,errmsg,errflg)
1523 : !-----------------------------------------------------------------------------------------
1524 : !
1525 : ! p(mb), Tfg/T(K), qt/qv(kg/kg), s(J/kg).
1526 : ! Inverts entropy, pressure and total water qt
1527 : ! for T and saturated vapor mixing ratio
1528 : !
1529 :
1530 : integer, intent(in) :: icol, rcall
1531 : real(kind_phys), intent(in) :: s, p, Tfg, qt
1532 : real(kind_phys), intent(in) :: cpliq
1533 : real(kind_phys), intent(in) :: cpwv
1534 : real(kind_phys), intent(in) :: rh2o
1535 :
1536 : real(kind_phys), intent(in) :: this_lat
1537 : real(kind_phys), intent(in) :: this_lon
1538 :
1539 : real(kind_phys), intent(out) :: qst, T
1540 : character(len=512), intent(out) :: errmsg
1541 : integer, intent(out) :: errflg
1542 :
1543 : real(kind_phys) :: est
1544 : real(kind_phys) :: a,b,c,d,ebr,fa,fb,fc,pbr,qbr,rbr,sbr,tol1,xm,tol
1545 : integer :: i
1546 :
1547 : logical :: converged
1548 :
1549 : ! Max number of iteration loops.
1550 : integer, parameter :: LOOPMAX = 100
1551 : real(kind_phys), parameter :: EPS = 3.e-8_kind_phys
1552 :
1553 90703569 : converged = .false.
1554 :
1555 : ! Invert the entropy equation -- use Brent's method
1556 : ! Brent, R. P. Ch. 3-4 in Algorithms for Minimization Without Derivatives. Englewood Cliffs, NJ: Prentice-Hall, 1973.
1557 :
1558 90703569 : T = Tfg ! Better first guess based on Tprofile from conv.
1559 :
1560 90703569 : a = Tfg-10 !low bracket
1561 90703569 : b = Tfg+10 !high bracket
1562 :
1563 90703569 : fa = entropy(a, p, qt,cpliq,cpwv,rh2o) - s
1564 90703569 : fb = entropy(b, p, qt,cpliq,cpwv,rh2o) - s
1565 :
1566 90703569 : c=b
1567 90703569 : fc=fb
1568 90703569 : tol=0.001_kind_phys
1569 :
1570 445537804 : converge: do i=0, LOOPMAX
1571 445537804 : if ((fb > 0.0_kind_phys .and. fc > 0.0_kind_phys) .or. &
1572 : (fb < 0.0_kind_phys .and. fc < 0.0_kind_phys)) then
1573 288738626 : c=a
1574 288738626 : fc=fa
1575 288738626 : d=b-a
1576 288738626 : ebr=d
1577 : end if
1578 445537804 : if (abs(fc) < abs(fb)) then
1579 105119017 : a=b
1580 105119017 : b=c
1581 105119017 : c=a
1582 105119017 : fa=fb
1583 105119017 : fb=fc
1584 105119017 : fc=fa
1585 : end if
1586 :
1587 445537804 : tol1=2.0_kind_phys*EPS*abs(b)+0.5_kind_phys*tol
1588 445537804 : xm=0.5_kind_phys*(c-b)
1589 445537804 : converged = (abs(xm) <= tol1 .or. fb == 0.0_kind_phys)
1590 445537804 : if (converged) exit converge
1591 :
1592 354834235 : if (abs(ebr) >= tol1 .and. abs(fa) > abs(fb)) then
1593 354834235 : sbr=fb/fa
1594 354834235 : if (a == c) then
1595 198037167 : pbr=2.0_kind_phys*xm*sbr
1596 198037167 : qbr=1.0_kind_phys-sbr
1597 : else
1598 156797068 : qbr=fa/fc
1599 156797068 : rbr=fb/fc
1600 156797068 : pbr=sbr*(2.0_kind_phys*xm*qbr*(qbr-rbr)-(b-a)*(rbr-1.0_kind_phys))
1601 156797068 : qbr=(qbr-1.0_kind_phys)*(rbr-1.0_kind_phys)*(sbr-1.0_kind_phys)
1602 : end if
1603 354834235 : if (pbr > 0.0_kind_phys) qbr=-qbr
1604 354834235 : pbr=abs(pbr)
1605 354834235 : if (2.0_kind_phys*pbr < min(3.0_kind_phys*xm*qbr-abs(tol1*qbr),abs(ebr*qbr))) then
1606 354717560 : ebr=d
1607 354717560 : d=pbr/qbr
1608 : else
1609 : d=xm
1610 : ebr=d
1611 : end if
1612 : else
1613 : d=xm
1614 : ebr=d
1615 : end if
1616 354834235 : a=b
1617 354834235 : fa=fb
1618 354834235 : b=b+merge(d,sign(tol1,xm), abs(d) > tol1 )
1619 :
1620 445537804 : fb = entropy(b, p, qt,cpliq,cpwv,rh2o) - s
1621 :
1622 : end do converge
1623 :
1624 90703569 : T = b
1625 90703569 : call qsat_hPa(T, p, est, qst)
1626 :
1627 90703569 : if (.not. converged) then
1628 0 : write(errmsg,100) ' ZM_CONV: IENTROPY. Details: call#,icol= ',rcall,icol, &
1629 0 : ' lat: ',this_lat,' lon: ',this_lon, &
1630 0 : ' P(mb)= ', p, ' Tfg(K)= ', Tfg, ' qt(g/kg) = ', 1000._kind_phys*qt, &
1631 0 : ' qst(g/kg) = ', 1000._kind_phys*qst,', s(J/kg) = ',s
1632 0 : errflg=1
1633 : end if
1634 :
1635 : 100 format (A,I4,I4,7(A,F6.2))
1636 :
1637 90703569 : end SUBROUTINE ientropy
1638 :
1639 80017 : subroutine cldprp(ncol ,pver ,pverp ,cpliq , &
1640 : latice ,cpwv ,rh2o ,&
1641 80017 : q ,t ,u ,v ,p , &
1642 80017 : z ,s ,mu ,eu ,du , &
1643 80017 : md ,ed ,sd ,qd ,mc , &
1644 80017 : qu ,su ,zf ,qst ,hmn , &
1645 80017 : hsat ,shat ,ql , &
1646 80017 : cmeg ,jb ,lel ,jt ,jlcl , &
1647 80017 : mx ,j0 ,jd ,rl ,il2g , &
1648 : rd ,grav ,cp ,msg , &
1649 80017 : evp ,cu ,rprd ,limcnv ,landfrac, &
1650 80017 : qcde ,qhat )
1651 : !----------------------------------------------
1652 : ! Purpose: Provide cloud properties
1653 : !----------------------------------------------
1654 :
1655 : implicit none
1656 :
1657 : !------------------------------------------------------------------------------
1658 : !
1659 : ! Input arguments
1660 : !
1661 : integer, intent(in) :: ncol
1662 : integer, intent(in) :: pver
1663 : integer, intent(in) :: pverp
1664 :
1665 : real(kind_phys), intent(in) :: cpliq
1666 : real(kind_phys), intent(in) :: latice
1667 : real(kind_phys), intent(in) :: cpwv
1668 : real(kind_phys), intent(in) :: rh2o
1669 :
1670 : real(kind_phys), intent(in) :: q(ncol,pver) ! spec. humidity of env
1671 : real(kind_phys), intent(in) :: t(ncol,pver) ! temp of env
1672 : real(kind_phys), intent(in) :: p(ncol,pver) ! pressure of env
1673 : real(kind_phys), intent(in) :: z(ncol,pver) ! height of env
1674 : real(kind_phys), intent(in) :: s(ncol,pver) ! normalized dry static energy of env
1675 : real(kind_phys), intent(in) :: zf(ncol,pverp) ! height of interfaces
1676 : real(kind_phys), intent(in) :: u(ncol,pver) ! zonal velocity of env
1677 : real(kind_phys), intent(in) :: v(ncol,pver) ! merid. velocity of env
1678 :
1679 : real(kind_phys), intent(in) :: landfrac(ncol) ! RBN Landfrac
1680 :
1681 : integer, intent(in) :: jb(ncol) ! updraft base level
1682 : integer, intent(in) :: lel(ncol) ! updraft launch level
1683 : integer, intent(in) :: mx(ncol) ! updraft base level (same is jb)
1684 : integer, intent(out) :: jt(ncol) ! updraft plume top
1685 : integer, intent(out) :: jlcl(ncol) ! updraft lifting cond level
1686 : integer, intent(out) :: j0(ncol) ! level where updraft begins detraining
1687 : integer, intent(out) :: jd(ncol) ! level of downdraft
1688 : integer, intent(in) :: limcnv ! convection limiting level
1689 : integer, intent(in) :: il2g !CORE GROUP REMOVE
1690 : integer, intent(in) :: msg ! missing moisture vals (always 0)
1691 : real(kind_phys), intent(in) :: rl ! latent heat of vap
1692 : real(kind_phys), intent(in) :: shat(ncol,pver) ! interface values of dry stat energy
1693 : real(kind_phys), intent(in) :: qhat(ncol,pver) ! wg grid slice of upper interface mixing ratio.
1694 :
1695 : !
1696 : ! output
1697 : !
1698 :
1699 : real(kind_phys), intent(out) :: rprd(ncol,pver) ! rate of production of precip at that layer
1700 : real(kind_phys), intent(out) :: du(ncol,pver) ! detrainement rate of updraft
1701 : real(kind_phys), intent(out) :: ed(ncol,pver) ! entrainment rate of downdraft
1702 : real(kind_phys), intent(out) :: eu(ncol,pver) ! entrainment rate of updraft
1703 : real(kind_phys), intent(out) :: hmn(ncol,pver) ! moist stat energy of env
1704 : real(kind_phys), intent(out) :: hsat(ncol,pver) ! sat moist stat energy of env
1705 : real(kind_phys), intent(out) :: mc(ncol,pver) ! net mass flux
1706 : real(kind_phys), intent(out) :: md(ncol,pver) ! downdraft mass flux
1707 : real(kind_phys), intent(out) :: mu(ncol,pver) ! updraft mass flux
1708 : real(kind_phys), intent(out) :: qd(ncol,pver) ! spec humidity of downdraft
1709 : real(kind_phys), intent(out) :: ql(ncol,pver) ! liq water of updraft
1710 : real(kind_phys), intent(out) :: qst(ncol,pver) ! saturation mixing ratio of env.
1711 : real(kind_phys), intent(out) :: qu(ncol,pver) ! spec hum of updraft
1712 : real(kind_phys), intent(out) :: sd(ncol,pver) ! normalized dry stat energy of downdraft
1713 : real(kind_phys), intent(out) :: su(ncol,pver) ! normalized dry stat energy of updraft
1714 : real(kind_phys), intent(out) :: qcde(ncol,pver) ! cloud water mixing ratio for detrainment (kg/kg)
1715 :
1716 : real(kind_phys) rd ! gas constant for dry air
1717 : real(kind_phys) grav ! gravity
1718 : real(kind_phys) cp ! heat capacity of dry air
1719 :
1720 : !
1721 : ! Local workspace
1722 : !
1723 160034 : real(kind_phys) gamma(ncol,pver)
1724 160034 : real(kind_phys) dz(ncol,pver)
1725 160034 : real(kind_phys) iprm(ncol,pver)
1726 160034 : real(kind_phys) hu(ncol,pver)
1727 160034 : real(kind_phys) hd(ncol,pver)
1728 160034 : real(kind_phys) eps(ncol,pver)
1729 160034 : real(kind_phys) f(ncol,pver)
1730 160034 : real(kind_phys) k1(ncol,pver)
1731 160034 : real(kind_phys) i2(ncol,pver)
1732 160034 : real(kind_phys) ihat(ncol,pver)
1733 160034 : real(kind_phys) i3(ncol,pver)
1734 160034 : real(kind_phys) idag(ncol,pver)
1735 160034 : real(kind_phys) i4(ncol,pver)
1736 160034 : real(kind_phys) qsthat(ncol,pver)
1737 160034 : real(kind_phys) hsthat(ncol,pver)
1738 160034 : real(kind_phys) gamhat(ncol,pver)
1739 : real(kind_phys) cu(ncol,pver)
1740 : real(kind_phys) evp(ncol,pver)
1741 : real(kind_phys) cmeg(ncol,pver)
1742 160034 : real(kind_phys) qds(ncol,pver)
1743 160034 : real(kind_phys) c0mask(ncol)
1744 :
1745 160034 : real(kind_phys) hmin(ncol)
1746 160034 : real(kind_phys) expdif(ncol)
1747 160034 : real(kind_phys) expnum(ncol)
1748 160034 : real(kind_phys) ftemp(ncol)
1749 160034 : real(kind_phys) eps0(ncol)
1750 160034 : real(kind_phys) rmue(ncol)
1751 160034 : real(kind_phys) zuef(ncol)
1752 160034 : real(kind_phys) zdef(ncol)
1753 160034 : real(kind_phys) epsm(ncol)
1754 160034 : real(kind_phys) ratmjb(ncol)
1755 160034 : real(kind_phys) est(ncol)
1756 160034 : real(kind_phys) totpcp(ncol)
1757 160034 : real(kind_phys) totevp(ncol)
1758 160034 : real(kind_phys) alfa(ncol)
1759 : real(kind_phys) ql1
1760 : real(kind_phys) tu
1761 : real(kind_phys) estu
1762 : real(kind_phys) qstu
1763 :
1764 : real(kind_phys) small
1765 : real(kind_phys) mdt
1766 :
1767 160034 : real(kind_phys) tug(ncol,pver)
1768 :
1769 160034 : real(kind_phys) tvuo(ncol,pver) ! updraft virtual T w/o freezing heating
1770 160034 : real(kind_phys) tvu(ncol,pver) ! updraft virtual T with freezing heating
1771 160034 : real(kind_phys) totfrz(ncol)
1772 160034 : real(kind_phys) frz (ncol,pver) ! rate of freezing
1773 160034 : integer jto(ncol) ! updraft plume old top
1774 160034 : integer tmplel(ncol)
1775 :
1776 : integer iter, itnum
1777 : integer m
1778 :
1779 : integer khighest
1780 : integer klowest
1781 : integer kount
1782 : integer i,k
1783 :
1784 160034 : logical doit(ncol)
1785 160034 : logical done(ncol)
1786 : !
1787 : !------------------------------------------------------------------------------
1788 : !
1789 :
1790 357782 : do i = 1,il2g
1791 277765 : ftemp(i) = 0._kind_phys
1792 277765 : expnum(i) = 0._kind_phys
1793 277765 : expdif(i) = 0._kind_phys
1794 357782 : c0mask(i) = c0_ocn * (1._kind_phys-landfrac(i)) + c0_lnd * landfrac(i)
1795 : end do
1796 : !
1797 : !jr Change from msg+1 to 1 to prevent blowup
1798 : !
1799 5681207 : do k = 1,pver
1800 25124757 : do i = 1,il2g
1801 25044740 : dz(i,k) = zf(i,k) - zf(i,k+1)
1802 : end do
1803 : end do
1804 :
1805 : !
1806 : ! initialize many output and work variables to zero
1807 : !
1808 5681207 : do k = 1,pver
1809 25124757 : do i = 1,il2g
1810 19443550 : k1(i,k) = 0._kind_phys
1811 19443550 : i2(i,k) = 0._kind_phys
1812 19443550 : i3(i,k) = 0._kind_phys
1813 19443550 : i4(i,k) = 0._kind_phys
1814 19443550 : mu(i,k) = 0._kind_phys
1815 19443550 : f(i,k) = 0._kind_phys
1816 19443550 : eps(i,k) = 0._kind_phys
1817 19443550 : eu(i,k) = 0._kind_phys
1818 19443550 : du(i,k) = 0._kind_phys
1819 19443550 : ql(i,k) = 0._kind_phys
1820 19443550 : cu(i,k) = 0._kind_phys
1821 19443550 : evp(i,k) = 0._kind_phys
1822 19443550 : cmeg(i,k) = 0._kind_phys
1823 19443550 : qds(i,k) = q(i,k)
1824 19443550 : md(i,k) = 0._kind_phys
1825 19443550 : ed(i,k) = 0._kind_phys
1826 19443550 : sd(i,k) = s(i,k)
1827 19443550 : qd(i,k) = q(i,k)
1828 19443550 : mc(i,k) = 0._kind_phys
1829 19443550 : qu(i,k) = q(i,k)
1830 19443550 : su(i,k) = s(i,k)
1831 19443550 : call qsat_hPa(t(i,k), p(i,k), est(i), qst(i,k))
1832 :
1833 19443550 : if ( p(i,k)-est(i) <= 0._kind_phys ) then
1834 6034892 : qst(i,k) = 1.0_kind_phys
1835 : end if
1836 :
1837 19443550 : gamma(i,k) = qst(i,k)*(1._kind_phys + qst(i,k)/eps1)*eps1*rl/(rd*t(i,k)**2)*rl/cp
1838 19443550 : hmn(i,k) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k)
1839 19443550 : hsat(i,k) = cp*t(i,k) + grav*z(i,k) + rl*qst(i,k)
1840 19443550 : hu(i,k) = hmn(i,k)
1841 19443550 : hd(i,k) = hmn(i,k)
1842 19443550 : rprd(i,k) = 0._kind_phys
1843 :
1844 19443550 : tug(i,k) = 0._kind_phys
1845 19443550 : qcde(i,k) = 0._kind_phys
1846 19443550 : tvuo(i,k) = (shat(i,k) - grav/cp*zf(i,k))*(1._kind_phys + 0.608_kind_phys*qhat(i,k))
1847 19443550 : tvu(i,k) = tvuo(i,k)
1848 44488290 : frz(i,k) = 0._kind_phys
1849 :
1850 : end do
1851 : end do
1852 :
1853 : !
1854 : !jr Set to zero things which make this routine blow up
1855 : !
1856 3520748 : do k=1,msg
1857 15464643 : do i=1,il2g
1858 15384626 : rprd(i,k) = 0._kind_phys
1859 : end do
1860 : end do
1861 : !
1862 : ! interpolate the layer values of qst, hsat and gamma to
1863 : ! layer interfaces
1864 : !
1865 3600765 : do k = 1, msg+1
1866 15822425 : do i = 1,il2g
1867 12221660 : hsthat(i,k) = hsat(i,k)
1868 12221660 : qsthat(i,k) = qst(i,k)
1869 15742408 : gamhat(i,k) = gamma(i,k)
1870 : end do
1871 : end do
1872 357782 : do i = 1,il2g
1873 277765 : totpcp(i) = 0._kind_phys
1874 357782 : totevp(i) = 0._kind_phys
1875 : end do
1876 2160459 : do k = msg + 2,pver
1877 9382349 : do i = 1,il2g
1878 7221890 : if (abs(qst(i,k-1)-qst(i,k)) > 1.E-6_kind_phys) then
1879 7065614 : qsthat(i,k) = log(qst(i,k-1)/qst(i,k))*qst(i,k-1)*qst(i,k)/ (qst(i,k-1)-qst(i,k))
1880 : else
1881 156276 : qsthat(i,k) = qst(i,k)
1882 : end if
1883 7221890 : hsthat(i,k) = cp*shat(i,k) + rl*qsthat(i,k)
1884 9302332 : if (abs(gamma(i,k-1)-gamma(i,k)) > 1.E-6_kind_phys) then
1885 : gamhat(i,k) = log(gamma(i,k-1)/gamma(i,k))*gamma(i,k-1)*gamma(i,k)/ &
1886 7221398 : (gamma(i,k-1)-gamma(i,k))
1887 : else
1888 492 : gamhat(i,k) = gamma(i,k)
1889 : end if
1890 : end do
1891 : end do
1892 : !
1893 : ! initialize cloud top to highest plume top.
1894 : !jr changed hard-wired 4 to limcnv+1 (not to exceed pver)
1895 : !
1896 1232336 : jt(:) = pver
1897 357782 : do i = 1,il2g
1898 277765 : jt(i) = max(lel(i),limcnv+1)
1899 277765 : jt(i) = min(jt(i),pver)
1900 277765 : jd(i) = pver
1901 277765 : jlcl(i) = lel(i)
1902 357782 : hmin(i) = 1.E6_kind_phys
1903 : end do
1904 : !
1905 : ! find the level of minimum hsat, where detrainment starts
1906 : !
1907 :
1908 2240476 : do k = msg + 1,pver
1909 9740131 : do i = 1,il2g
1910 9660114 : if (hsat(i,k) <= hmin(i) .and. k >= jt(i) .and. k <= jb(i)) then
1911 770301 : hmin(i) = hsat(i,k)
1912 770301 : j0(i) = k
1913 : end if
1914 : end do
1915 : end do
1916 357782 : do i = 1,il2g
1917 277765 : j0(i) = min(j0(i),jb(i)-2)
1918 277765 : j0(i) = max(j0(i),jt(i)+2)
1919 : !
1920 : ! Fix from Guang Zhang to address out of bounds array reference
1921 : !
1922 357782 : j0(i) = min(j0(i),pver)
1923 : end do
1924 : !
1925 : ! Initialize certain arrays inside cloud
1926 : !
1927 2240476 : do k = msg + 1,pver
1928 9740131 : do i = 1,il2g
1929 9660114 : if (k >= jt(i) .and. k <= jb(i)) then
1930 3004126 : hu(i,k) = hmn(i,mx(i)) + cp*tiedke_add
1931 3004126 : su(i,k) = s(i,mx(i)) + tiedke_add
1932 : end if
1933 : end do
1934 : end do
1935 : !
1936 : ! *********************************************************
1937 : ! compute taylor series for approximate eps(z) below
1938 : ! *********************************************************
1939 : !
1940 2160459 : do k = pver - 1,msg + 1,-1
1941 9382349 : do i = 1,il2g
1942 9302332 : if (k < jb(i) .and. k >= jt(i)) then
1943 2726361 : k1(i,k) = k1(i,k+1) + (hmn(i,mx(i))-hmn(i,k))*dz(i,k)
1944 2726361 : ihat(i,k) = 0.5_kind_phys* (k1(i,k+1)+k1(i,k))
1945 2726361 : i2(i,k) = i2(i,k+1) + ihat(i,k)*dz(i,k)
1946 2726361 : idag(i,k) = 0.5_kind_phys* (i2(i,k+1)+i2(i,k))
1947 2726361 : i3(i,k) = i3(i,k+1) + idag(i,k)*dz(i,k)
1948 2726361 : iprm(i,k) = 0.5_kind_phys* (i3(i,k+1)+i3(i,k))
1949 2726361 : i4(i,k) = i4(i,k+1) + iprm(i,k)*dz(i,k)
1950 : end if
1951 : end do
1952 : end do
1953 : !
1954 : ! re-initialize hmin array for ensuing calculation.
1955 : !
1956 357782 : do i = 1,il2g
1957 357782 : hmin(i) = 1.E6_kind_phys
1958 : end do
1959 2240476 : do k = msg + 1,pver
1960 9740131 : do i = 1,il2g
1961 9660114 : if (k >= j0(i) .and. k <= jb(i) .and. hmn(i,k) <= hmin(i)) then
1962 363396 : hmin(i) = hmn(i,k)
1963 363396 : expdif(i) = hmn(i,mx(i)) - hmin(i)
1964 : end if
1965 : end do
1966 : end do
1967 : !
1968 : ! *********************************************************
1969 : ! compute approximate eps(z) using above taylor series
1970 : ! *********************************************************
1971 : !
1972 2160459 : do k = msg + 2,pver
1973 9382349 : do i = 1,il2g
1974 7221890 : expnum(i) = 0._kind_phys
1975 7221890 : ftemp(i) = 0._kind_phys
1976 7221890 : if (k < jt(i) .or. k >= jb(i)) then
1977 4495529 : k1(i,k) = 0._kind_phys
1978 4495529 : expnum(i) = 0._kind_phys
1979 : else
1980 8179083 : expnum(i) = hmn(i,mx(i)) - (hsat(i,k-1)*(zf(i,k)-z(i,k)) + &
1981 10905444 : hsat(i,k)* (z(i,k-1)-zf(i,k)))/(z(i,k-1)-z(i,k))
1982 : end if
1983 7221890 : if ((expdif(i) > 100._kind_phys .and. expnum(i) > 0._kind_phys) .and. &
1984 2080442 : k1(i,k) > expnum(i)*dz(i,k)) then
1985 1888279 : ftemp(i) = expnum(i)/k1(i,k)
1986 : f(i,k) = ftemp(i) + i2(i,k)/k1(i,k)*ftemp(i)**2 + &
1987 : (2._kind_phys*i2(i,k)**2-k1(i,k)*i3(i,k))/k1(i,k)**2* &
1988 : ftemp(i)**3 + (-5._kind_phys*k1(i,k)*i2(i,k)*i3(i,k)+ &
1989 : 5._kind_phys*i2(i,k)**3+k1(i,k)**2*i4(i,k))/ &
1990 1888279 : k1(i,k)**3*ftemp(i)**4
1991 1888279 : f(i,k) = max(f(i,k),0._kind_phys)
1992 1888279 : f(i,k) = min(f(i,k),0.0002_kind_phys)
1993 : end if
1994 : end do
1995 : end do
1996 357782 : do i = 1,il2g
1997 357782 : if (j0(i) < jb(i)) then
1998 277765 : if (f(i,j0(i)) < 1.E-6_kind_phys .and. f(i,j0(i)+1) > f(i,j0(i))) j0(i) = j0(i) + 1
1999 : end if
2000 : end do
2001 2160459 : do k = msg + 2,pver
2002 9382349 : do i = 1,il2g
2003 9302332 : if (k >= jt(i) .and. k <= j0(i)) then
2004 958153 : f(i,k) = max(f(i,k),f(i,k-1))
2005 : end if
2006 : end do
2007 : end do
2008 357782 : do i = 1,il2g
2009 277765 : eps0(i) = f(i,j0(i))
2010 357782 : eps(i,jb(i)) = eps0(i)
2011 : end do
2012 : !
2013 : ! This is set to match the Rasch and Kristjansson paper
2014 : !
2015 2240476 : do k = pver,msg + 1,-1
2016 9740131 : do i = 1,il2g
2017 9660114 : if (k >= j0(i) .and. k <= jb(i)) then
2018 2323738 : eps(i,k) = f(i,j0(i))
2019 : end if
2020 : end do
2021 : end do
2022 2240476 : do k = pver,msg + 1,-1
2023 9740131 : do i = 1,il2g
2024 9660114 : if (k < j0(i) .and. k >= jt(i)) eps(i,k) = f(i,k)
2025 : end do
2026 : end do
2027 :
2028 : itnum = 1
2029 160034 : do iter=1, itnum
2030 :
2031 : !
2032 : ! specify the updraft mass flux mu, entrainment eu, detrainment du
2033 : ! and moist static energy hu.
2034 : ! here and below mu, eu,du, md and ed are all normalized by mb
2035 : !
2036 357782 : do i = 1,il2g
2037 277765 : if (eps0(i) > 0._kind_phys) then
2038 276413 : mu(i,jb(i)) = 1._kind_phys
2039 276413 : eu(i,jb(i)) = mu(i,jb(i))/dz(i,jb(i))
2040 : end if
2041 :
2042 357782 : tmplel(i) = jt(i)
2043 : end do
2044 2240476 : do k = pver,msg + 1,-1
2045 9740131 : do i = 1,il2g
2046 9660114 : if (eps0(i) > 0._kind_phys .and. (k >= tmplel(i) .and. k < jb(i))) then
2047 2710998 : zuef(i) = zf(i,k) - zf(i,jb(i))
2048 2710998 : rmue(i) = (1._kind_phys/eps0(i))* (exp(eps(i,k+1)*zuef(i))-1._kind_phys)/zuef(i)
2049 2710998 : mu(i,k) = (1._kind_phys/eps0(i))* (exp(eps(i,k )*zuef(i))-1._kind_phys)/zuef(i)
2050 2710998 : eu(i,k) = (rmue(i)-mu(i,k+1))/dz(i,k)
2051 2710998 : du(i,k) = (rmue(i)-mu(i,k))/dz(i,k)
2052 : end if
2053 : end do
2054 : end do
2055 :
2056 : khighest = pverp
2057 : klowest = 1
2058 357782 : do i=1,il2g
2059 277765 : khighest = min(khighest,lel(i))
2060 357782 : klowest = max(klowest,jb(i))
2061 : end do
2062 985131 : do k = klowest-1,khighest,-1
2063 4207720 : do i = 1,il2g
2064 4127703 : if (k <= jb(i)-1 .and. k >= lel(i) .and. eps0(i) > 0._kind_phys) then
2065 2710998 : if (mu(i,k) < 0.02_kind_phys) then
2066 182063 : hu(i,k) = hmn(i,k)
2067 182063 : mu(i,k) = 0._kind_phys
2068 182063 : eu(i,k) = 0._kind_phys
2069 182063 : du(i,k) = mu(i,k+1)/dz(i,k)
2070 : else
2071 2528935 : hu(i,k) = mu(i,k+1)/mu(i,k)*hu(i,k+1) + &
2072 5057870 : dz(i,k)/mu(i,k)* (eu(i,k)*hmn(i,k)- du(i,k)*hsat(i,k))
2073 : end if
2074 : end if
2075 : end do
2076 : end do
2077 : !
2078 : ! reset cloud top index beginning from two layers above the
2079 : ! cloud base (i.e. if cloud is only one layer thick, top is not reset
2080 : !
2081 357782 : do i=1,il2g
2082 277765 : doit(i) = .true.
2083 277765 : totfrz(i)= 0._kind_phys
2084 7857437 : do k = pver,msg + 1,-1
2085 7777420 : totfrz(i)= totfrz(i)+ frz(i,k)*dz(i,k)
2086 : end do
2087 : end do
2088 985131 : do k=klowest-2,khighest-1,-1
2089 4207720 : do i=1,il2g
2090 4127703 : if (doit(i) .and. k <= jb(i)-2 .and. k >= lel(i)-1) then
2091 4899660 : if (hu(i,k) <= hsthat(i,k) .and. hu(i,k+1) > hsthat(i,k+1) &
2092 7349490 : .and. mu(i,k) >= 0.02_kind_phys) then
2093 14113 : if (hu(i,k)-hsthat(i,k) < -2000._kind_phys) then
2094 8537 : jt(i) = k + 1
2095 8537 : doit(i) = .false.
2096 : else
2097 5576 : jt(i) = k
2098 5576 : doit(i) = .false.
2099 : end if
2100 2435717 : else if ( (hu(i,k) > hu(i,jb(i)) .and. totfrz(i)<=0._kind_phys) .or. mu(i,k) < 0.02_kind_phys) then
2101 263652 : jt(i) = k + 1
2102 263652 : doit(i) = .false.
2103 : end if
2104 : end if
2105 : end do
2106 : end do
2107 :
2108 1232336 : if (iter == 1) jto(:) = jt(:)
2109 :
2110 2240476 : do k = pver,msg + 1,-1
2111 9740131 : do i = 1,il2g
2112 7499655 : if (k >= lel(i) .and. k <= jt(i) .and. eps0(i) > 0._kind_phys) then
2113 533357 : mu(i,k) = 0._kind_phys
2114 533357 : eu(i,k) = 0._kind_phys
2115 533357 : du(i,k) = 0._kind_phys
2116 533357 : hu(i,k) = hmn(i,k)
2117 : end if
2118 9660114 : if (k == jt(i) .and. eps0(i) > 0._kind_phys) then
2119 276413 : du(i,k) = mu(i,k+1)/dz(i,k)
2120 276413 : eu(i,k) = 0._kind_phys
2121 276413 : mu(i,k) = 0._kind_phys
2122 : end if
2123 : end do
2124 : end do
2125 :
2126 357782 : do i = 1,il2g
2127 357782 : done(i) = .false.
2128 : end do
2129 : kount = 0
2130 495659 : do k = pver,msg + 2,-1
2131 2338972 : do i = 1,il2g
2132 1847768 : if (k == jb(i) .and. eps0(i) > 0._kind_phys) then
2133 276413 : qu(i,k) = q(i,mx(i))
2134 276413 : su(i,k) = (hu(i,k)-rl*qu(i,k))/cp
2135 : end if
2136 2338972 : if (( .not. done(i) .and. k > jt(i) .and. k < jb(i)) .and. eps0(i) > 0._kind_phys) then
2137 1563768 : su(i,k) = mu(i,k+1)/mu(i,k)*su(i,k+1) + &
2138 1563768 : dz(i,k)/mu(i,k)* (eu(i,k)-du(i,k))*s(i,k)
2139 : qu(i,k) = mu(i,k+1)/mu(i,k)*qu(i,k+1) + dz(i,k)/mu(i,k)* (eu(i,k)*q(i,k)- &
2140 781884 : du(i,k)*qst(i,k))
2141 781884 : tu = su(i,k) - grav/cp*zf(i,k)
2142 781884 : call qsat_hPa(tu, (p(i,k)+p(i,k-1))/2._kind_phys, estu, qstu)
2143 781884 : if (qu(i,k) >= qstu) then
2144 273180 : jlcl(i) = k
2145 273180 : kount = kount + 1
2146 273180 : done(i) = .true.
2147 : end if
2148 : end if
2149 : end do
2150 495659 : if (kount >= il2g) goto 690
2151 : end do
2152 : 690 continue
2153 2160459 : do k = msg + 2,pver
2154 9382349 : do i = 1,il2g
2155 9302332 : if ((k > jt(i) .and. k <= jlcl(i)) .and. eps0(i) > 0._kind_phys) then
2156 1668937 : su(i,k) = shat(i,k) + (hu(i,k)-hsthat(i,k))/(cp* (1._kind_phys+gamhat(i,k)))
2157 : qu(i,k) = qsthat(i,k) + gamhat(i,k)*(hu(i,k)-hsthat(i,k))/ &
2158 1668937 : (rl* (1._kind_phys+gamhat(i,k)))
2159 : end if
2160 : end do
2161 : end do
2162 :
2163 : ! compute condensation in updraft
2164 357782 : tmplel(:il2g) = jb(:il2g)
2165 :
2166 2160459 : do k = pver,msg + 2,-1
2167 9382349 : do i = 1,il2g
2168 9302332 : if (k >= jt(i) .and. k < tmplel(i) .and. eps0(i) > 0._kind_phys) then
2169 :
2170 4908108 : cu(i,k) = ((mu(i,k)*su(i,k)-mu(i,k+1)*su(i,k+1))/ &
2171 4908108 : dz(i,k)- (eu(i,k)-du(i,k))*s(i,k))/(rl/cp)
2172 2454054 : if (k == jt(i)) cu(i,k) = 0._kind_phys
2173 2454054 : cu(i,k) = max(0._kind_phys,cu(i,k))
2174 : end if
2175 : end do
2176 : end do
2177 :
2178 :
2179 : ! compute condensed liquid, rain production rate
2180 : ! accumulate total precipitation (condensation - detrainment of liquid)
2181 : ! Note ql1 = ql(k) + rprd(k)*dz(k)/mu(k)
2182 : ! The differencing is somewhat strange (e.g. du(i,k)*ql(i,k+1)) but is
2183 : ! consistently applied.
2184 : ! mu, ql are interface quantities
2185 : ! cu, du, eu, rprd are midpoint quantites
2186 :
2187 2240476 : do k = pver,msg + 2,-1
2188 9382349 : do i = 1,il2g
2189 7221890 : rprd(i,k) = 0._kind_phys
2190 9302332 : if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._kind_phys .and. mu(i,k) >= 0.0_kind_phys) then
2191 2454054 : if (mu(i,k) > 0._kind_phys) then
2192 2177641 : ql1 = 1._kind_phys/mu(i,k)* (mu(i,k+1)*ql(i,k+1)- &
2193 2177641 : dz(i,k)*du(i,k)*ql(i,k+1)+dz(i,k)*cu(i,k))
2194 2177641 : ql(i,k) = ql1/ (1._kind_phys+dz(i,k)*c0mask(i))
2195 : else
2196 276413 : ql(i,k) = 0._kind_phys
2197 : end if
2198 2454054 : totpcp(i) = totpcp(i) + dz(i,k)*(cu(i,k)-du(i,k)*ql(i,k+1))
2199 2454054 : rprd(i,k) = c0mask(i)*mu(i,k)*ql(i,k)
2200 2454054 : qcde(i,k) = ql(i,k)
2201 :
2202 :
2203 : end if
2204 : end do
2205 : end do
2206 : !
2207 : end do !iter
2208 : !
2209 : ! specify downdraft properties (no downdrafts if jd.ge.jb).
2210 : ! scale down downward mass flux profile so that net flux
2211 : ! (up-down) at cloud base in not negative.
2212 : !
2213 357782 : do i = 1,il2g
2214 : !
2215 : ! in normal downdraft strength run alfa=0.2. In test4 alfa=0.1
2216 : !
2217 277765 : alfa(i) = 0.1_kind_phys
2218 277765 : jt(i) = min(jt(i),jb(i)-1)
2219 277765 : jd(i) = max(j0(i),jt(i)+1)
2220 277765 : jd(i) = min(jd(i),jb(i))
2221 277765 : hd(i,jd(i)) = hmn(i,jd(i)-1)
2222 357782 : if (jd(i) < jb(i) .and. eps0(i) > 0._kind_phys) then
2223 274270 : epsm(i) = eps0(i)
2224 274270 : md(i,jd(i)) = -alfa(i)*epsm(i)/eps0(i)
2225 : end if
2226 : end do
2227 2240476 : do k = msg + 1,pver
2228 9740131 : do i = 1,il2g
2229 9660114 : if ((k > jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._kind_phys) then
2230 1949865 : zdef(i) = zf(i,jd(i)) - zf(i,k)
2231 1949865 : md(i,k) = -alfa(i)/ (2._kind_phys*eps0(i))*(exp(2._kind_phys*epsm(i)*zdef(i))-1._kind_phys)/zdef(i)
2232 : end if
2233 : end do
2234 : end do
2235 :
2236 2240476 : do k = msg + 1,pver
2237 9740131 : do i = 1,il2g
2238 9660114 : if ((k >= jt(i) .and. k <= jb(i)) .and. eps0(i) > 0._kind_phys .and. jd(i) < jb(i)) then
2239 2726181 : ratmjb(i) = min(abs(mu(i,jb(i))/md(i,jb(i))),1._kind_phys)
2240 2726181 : md(i,k) = md(i,k)*ratmjb(i)
2241 : end if
2242 : end do
2243 : end do
2244 :
2245 : small = 1.e-20_kind_phys
2246 2240476 : do k = msg + 1,pver
2247 9740131 : do i = 1,il2g
2248 9660114 : if ((k >= jt(i) .and. k <= pver) .and. eps0(i) > 0._kind_phys) then
2249 2749185 : ed(i,k-1) = (md(i,k-1)-md(i,k))/dz(i,k-1)
2250 2749185 : mdt = min(md(i,k),-small)
2251 2749185 : hd(i,k) = (md(i,k-1)*hd(i,k-1) - dz(i,k-1)*ed(i,k-1)*hmn(i,k-1))/mdt
2252 : end if
2253 : end do
2254 : end do
2255 : !
2256 : ! calculate updraft and downdraft properties.
2257 : !
2258 2160459 : do k = msg + 2,pver
2259 9382349 : do i = 1,il2g
2260 9302332 : if ((k >= jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._kind_phys .and. jd(i) < jb(i)) then
2261 2224135 : qds(i,k) = qsthat(i,k) + gamhat(i,k)*(hd(i,k)-hsthat(i,k))/ &
2262 4448270 : (rl*(1._kind_phys + gamhat(i,k)))
2263 : end if
2264 : end do
2265 : end do
2266 :
2267 357782 : do i = 1,il2g
2268 277765 : qd(i,jd(i)) = qds(i,jd(i))
2269 357782 : sd(i,jd(i)) = (hd(i,jd(i)) - rl*qd(i,jd(i)))/cp
2270 : end do
2271 : !
2272 2160459 : do k = msg + 2,pver
2273 9382349 : do i = 1,il2g
2274 9302332 : if (k >= jd(i) .and. k < jb(i) .and. eps0(i) > 0._kind_phys) then
2275 1949865 : qd(i,k+1) = qds(i,k+1)
2276 1949865 : evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k)-md(i,k+1)*qd(i,k+1))/dz(i,k)
2277 1949865 : evp(i,k) = max(evp(i,k),0._kind_phys)
2278 1949865 : mdt = min(md(i,k+1),-small)
2279 1949865 : sd(i,k+1) = ((rl/cp*evp(i,k)-ed(i,k)*s(i,k))*dz(i,k) + md(i,k)*sd(i,k))/mdt
2280 1949865 : totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k)
2281 : end if
2282 : end do
2283 : end do
2284 357782 : do i = 1,il2g
2285 357782 : totevp(i) = totevp(i) + md(i,jd(i))*qd(i,jd(i)) - md(i,jb(i))*qd(i,jb(i))
2286 : end do
2287 : !!$ if (.true.) then
2288 : if (.false.) then
2289 : do i = 1,il2g
2290 : k = jb(i)
2291 : if (eps0(i) > 0._kind_phys) then
2292 : evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k))/dz(i,k)
2293 : evp(i,k) = max(evp(i,k),0._kind_phys)
2294 : totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k)
2295 : end if
2296 : end do
2297 : endif
2298 :
2299 357782 : do i = 1,il2g
2300 277765 : totpcp(i) = max(totpcp(i),0._kind_phys)
2301 357782 : totevp(i) = max(totevp(i),0._kind_phys)
2302 : end do
2303 : !
2304 2160459 : do k = msg + 2,pver
2305 9382349 : do i = 1,il2g
2306 7221890 : if (totevp(i) > 0._kind_phys .and. totpcp(i) > 0._kind_phys) then
2307 7129564 : md(i,k) = md (i,k)*min(1._kind_phys, totpcp(i)/(totevp(i)+totpcp(i)))
2308 7129564 : ed(i,k) = ed (i,k)*min(1._kind_phys, totpcp(i)/(totevp(i)+totpcp(i)))
2309 7129564 : evp(i,k) = evp(i,k)*min(1._kind_phys, totpcp(i)/(totevp(i)+totpcp(i)))
2310 : else
2311 92326 : md(i,k) = 0._kind_phys
2312 92326 : ed(i,k) = 0._kind_phys
2313 92326 : evp(i,k) = 0._kind_phys
2314 : end if
2315 : ! cmeg is the cloud water condensed - rain water evaporated
2316 : ! rprd is the cloud water converted to rain - (rain evaporated)
2317 7221890 : cmeg(i,k) = cu(i,k) - evp(i,k)
2318 9302332 : rprd(i,k) = rprd(i,k)-evp(i,k)
2319 : end do
2320 : end do
2321 :
2322 : !
2323 2240476 : do k = msg + 1,pver
2324 9740131 : do i = 1,il2g
2325 9660114 : mc(i,k) = mu(i,k) + md(i,k)
2326 : end do
2327 : end do
2328 : !
2329 80017 : return
2330 : end subroutine cldprp
2331 :
2332 80017 : subroutine closure(ncol ,pver, &
2333 80017 : q ,t ,p ,z ,s , &
2334 80017 : tp ,qs ,qu ,su ,mc , &
2335 80017 : du ,mu ,md ,qd ,sd , &
2336 80017 : qhat ,shat ,dp ,qstp ,zf , &
2337 80017 : ql ,dsubcld ,mb ,cape ,tl , &
2338 80017 : lcl ,lel ,jt ,mx ,il1g , &
2339 : il2g ,rd ,grav ,cp ,rl , &
2340 : msg ,capelmt )
2341 : !
2342 : !-----------------------------Arguments---------------------------------
2343 : !
2344 : integer, intent(in) :: ncol
2345 : integer, intent(in) :: pver
2346 :
2347 : real(kind_phys), intent(inout) :: q(ncol,pver) ! spec humidity
2348 : real(kind_phys), intent(inout) :: t(ncol,pver) ! temperature
2349 : real(kind_phys), intent(inout) :: p(ncol,pver) ! pressure (mb)
2350 : real(kind_phys), intent(inout) :: mb(ncol) ! cloud base mass flux
2351 : real(kind_phys), intent(in) :: z(ncol,pver) ! height (m)
2352 : real(kind_phys), intent(in) :: s(ncol,pver) ! normalized dry static energy
2353 : real(kind_phys), intent(in) :: tp(ncol,pver) ! parcel temp
2354 : real(kind_phys), intent(in) :: qs(ncol,pver) ! sat spec humidity
2355 : real(kind_phys), intent(in) :: qu(ncol,pver) ! updraft spec. humidity
2356 : real(kind_phys), intent(in) :: su(ncol,pver) ! normalized dry stat energy of updraft
2357 : real(kind_phys), intent(in) :: mc(ncol,pver) ! net convective mass flux
2358 : real(kind_phys), intent(in) :: du(ncol,pver) ! detrainment from updraft
2359 : real(kind_phys), intent(in) :: mu(ncol,pver) ! mass flux of updraft
2360 : real(kind_phys), intent(in) :: md(ncol,pver) ! mass flux of downdraft
2361 : real(kind_phys), intent(in) :: qd(ncol,pver) ! spec. humidity of downdraft
2362 : real(kind_phys), intent(in) :: sd(ncol,pver) ! dry static energy of downdraft
2363 : real(kind_phys), intent(in) :: qhat(ncol,pver) ! environment spec humidity at interfaces
2364 : real(kind_phys), intent(in) :: shat(ncol,pver) ! env. normalized dry static energy at intrfcs
2365 : real(kind_phys), intent(in) :: dp(ncol,pver) ! pressure thickness of layers
2366 : real(kind_phys), intent(in) :: qstp(ncol,pver) ! spec humidity of parcel
2367 : real(kind_phys), intent(in) :: zf(ncol,pver+1) ! height of interface levels
2368 : real(kind_phys), intent(in) :: ql(ncol,pver) ! liquid water mixing ratio
2369 :
2370 : real(kind_phys), intent(in) :: cape(ncol) ! available pot. energy of column
2371 : real(kind_phys), intent(in) :: tl(ncol)
2372 : real(kind_phys), intent(in) :: dsubcld(ncol) ! thickness of subcloud layer
2373 :
2374 : integer, intent(in) :: lcl(ncol) ! index of lcl
2375 : integer, intent(in) :: lel(ncol) ! index of launch leve
2376 : integer, intent(in) :: jt(ncol) ! top of updraft
2377 : integer, intent(in) :: mx(ncol) ! base of updraft
2378 : !
2379 : !--------------------------Local variables------------------------------
2380 : !
2381 160034 : real(kind_phys) dtpdt(ncol,pver)
2382 160034 : real(kind_phys) dqsdtp(ncol,pver)
2383 160034 : real(kind_phys) dtmdt(ncol,pver)
2384 160034 : real(kind_phys) dqmdt(ncol,pver)
2385 160034 : real(kind_phys) dboydt(ncol,pver)
2386 160034 : real(kind_phys) thetavp(ncol,pver)
2387 160034 : real(kind_phys) thetavm(ncol,pver)
2388 :
2389 160034 : real(kind_phys) dtbdt(ncol),dqbdt(ncol),dtldt(ncol)
2390 : real(kind_phys) beta
2391 : real(kind_phys) capelmt
2392 : real(kind_phys) cp
2393 160034 : real(kind_phys) dadt(ncol)
2394 : real(kind_phys) debdt
2395 : real(kind_phys) dltaa
2396 : real(kind_phys) eb
2397 : real(kind_phys) grav
2398 :
2399 : integer i
2400 : integer il1g
2401 : integer il2g
2402 : integer k, kmin, kmax
2403 : integer msg
2404 :
2405 : real(kind_phys) rd
2406 : real(kind_phys) rl
2407 : ! change of subcloud layer properties due to convection is
2408 : ! related to cumulus updrafts and downdrafts.
2409 : ! mc(z)=f(z)*mb, mub=betau*mb, mdb=betad*mb are used
2410 : ! to define betau, betad and f(z).
2411 : ! note that this implies all time derivatives are in effect
2412 : ! time derivatives per unit cloud-base mass flux, i.e. they
2413 : ! have units of 1/mb instead of 1/sec.
2414 : !
2415 357782 : do i = il1g,il2g
2416 277765 : mb(i) = 0._kind_phys
2417 277765 : eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
2418 : dtbdt(i) = (1._kind_phys/dsubcld(i))* (mu(i,mx(i))*(shat(i,mx(i))-su(i,mx(i)))+ &
2419 277765 : md(i,mx(i))* (shat(i,mx(i))-sd(i,mx(i))))
2420 : dqbdt(i) = (1._kind_phys/dsubcld(i))* (mu(i,mx(i))*(qhat(i,mx(i))-qu(i,mx(i)))+ &
2421 277765 : md(i,mx(i))* (qhat(i,mx(i))-qd(i,mx(i))))
2422 277765 : debdt = eps1*p(i,mx(i))/ (eps1+q(i,mx(i)))**2*dqbdt(i)
2423 : dtldt(i) = -2840._kind_phys* (3.5_kind_phys/t(i,mx(i))*dtbdt(i)-debdt/eb)/ &
2424 357782 : (3.5_kind_phys*log(t(i,mx(i)))-log(eb)-4.805_kind_phys)**2
2425 : end do
2426 : !
2427 : ! dtmdt and dqmdt are cumulus heating and drying.
2428 : !
2429 2240476 : do k = msg + 1,pver
2430 9740131 : do i = il1g,il2g
2431 7499655 : dtmdt(i,k) = 0._kind_phys
2432 9660114 : dqmdt(i,k) = 0._kind_phys
2433 : end do
2434 : end do
2435 : !
2436 2160459 : do k = msg + 1,pver - 1
2437 9382349 : do i = il1g,il2g
2438 9302332 : if (k == jt(i)) then
2439 555530 : dtmdt(i,k) = (1._kind_phys/dp(i,k))*(mu(i,k+1)* (su(i,k+1)-shat(i,k+1)- &
2440 555530 : rl/cp*ql(i,k+1))+md(i,k+1)* (sd(i,k+1)-shat(i,k+1)))
2441 : dqmdt(i,k) = (1._kind_phys/dp(i,k))*(mu(i,k+1)* (qu(i,k+1)- &
2442 277765 : qhat(i,k+1)+ql(i,k+1))+md(i,k+1)*(qd(i,k+1)-qhat(i,k+1)))
2443 : end if
2444 : end do
2445 : end do
2446 : !
2447 80017 : beta = 0._kind_phys
2448 2160459 : do k = msg + 1,pver - 1
2449 9382349 : do i = il1g,il2g
2450 9302332 : if (k > jt(i) .and. k < mx(i)) then
2451 4355282 : dtmdt(i,k) = (mc(i,k)* (shat(i,k)-s(i,k))+mc(i,k+1)* (s(i,k)-shat(i,k+1)))/ &
2452 4355282 : dp(i,k) - rl/cp*du(i,k)*(beta*ql(i,k)+ (1-beta)*ql(i,k+1))
2453 :
2454 : dqmdt(i,k) = (mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)+cp/rl* (su(i,k+1)-s(i,k)))- &
2455 : mu(i,k)* (qu(i,k)-qhat(i,k)+cp/rl*(su(i,k)-s(i,k)))+md(i,k+1)* &
2456 : (qd(i,k+1)-qhat(i,k+1)+cp/rl*(sd(i,k+1)-s(i,k)))-md(i,k)* &
2457 : (qd(i,k)-qhat(i,k)+cp/rl*(sd(i,k)-s(i,k))))/dp(i,k) + &
2458 2177641 : du(i,k)* (beta*ql(i,k)+(1-beta)*ql(i,k+1))
2459 : end if
2460 : end do
2461 : end do
2462 : !
2463 2240476 : do k = msg + 1,pver
2464 9740131 : do i = il1g,il2g
2465 9660114 : if (k >= lel(i) .and. k <= lcl(i)) then
2466 2125841 : thetavp(i,k) = tp(i,k)* (1000._kind_phys/p(i,k))** (rd/cp)*(1._kind_phys+1.608_kind_phys*qstp(i,k)-q(i,mx(i)))
2467 2125841 : thetavm(i,k) = t(i,k)* (1000._kind_phys/p(i,k))** (rd/cp)*(1._kind_phys+0.608_kind_phys*q(i,k))
2468 2125841 : dqsdtp(i,k) = qstp(i,k)* (1._kind_phys+qstp(i,k)/eps1)*eps1*rl/(rd*tp(i,k)**2)
2469 : !
2470 : ! dtpdt is the parcel temperature change due to change of
2471 : ! subcloud layer properties during convection.
2472 : !
2473 : dtpdt(i,k) = tp(i,k)/ (1._kind_phys+rl/cp* (dqsdtp(i,k)-qstp(i,k)/tp(i,k)))* &
2474 : (dtbdt(i)/t(i,mx(i))+rl/cp* (dqbdt(i)/tl(i)-q(i,mx(i))/ &
2475 2125841 : tl(i)**2*dtldt(i)))
2476 : !
2477 : ! dboydt is the integrand of cape change.
2478 : !
2479 : dboydt(i,k) = ((dtpdt(i,k)/tp(i,k)+1._kind_phys/(1._kind_phys+1.608_kind_phys*qstp(i,k)-q(i,mx(i)))* &
2480 : (1.608_kind_phys * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608_kind_phys/ &
2481 2125841 : (1._kind_phys+0.608_kind_phys*q(i,k))*dqmdt(i,k)))*grav*thetavp(i,k)/thetavm(i,k)
2482 : end if
2483 : end do
2484 : end do
2485 : !
2486 2240476 : do k = msg + 1,pver
2487 9740131 : do i = il1g,il2g
2488 9660114 : if (k > lcl(i) .and. k < mx(i)) then
2489 602241 : thetavp(i,k) = tp(i,k)* (1000._kind_phys/p(i,k))** (rd/cp)*(1._kind_phys+0.608_kind_phys*q(i,mx(i)))
2490 602241 : thetavm(i,k) = t(i,k)* (1000._kind_phys/p(i,k))** (rd/cp)*(1._kind_phys+0.608_kind_phys*q(i,k))
2491 : !
2492 : ! dboydt is the integrand of cape change.
2493 : !
2494 : dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608_kind_phys/ (1._kind_phys+0.608_kind_phys*q(i,mx(i)))*dqbdt(i)- &
2495 : dtmdt(i,k)/t(i,k)-0.608_kind_phys/ (1._kind_phys+0.608_kind_phys*q(i,k))*dqmdt(i,k))* &
2496 602241 : grav*thetavp(i,k)/thetavm(i,k)
2497 : end if
2498 : end do
2499 : end do
2500 :
2501 : !
2502 : ! buoyant energy change is set to 2/3*excess cape per 3 hours
2503 : !
2504 357782 : dadt(il1g:il2g) = 0._kind_phys
2505 357782 : kmin = minval(lel(il1g:il2g))
2506 357782 : kmax = maxval(mx(il1g:il2g)) - 1
2507 985131 : do k = kmin, kmax
2508 4207720 : do i = il1g,il2g
2509 4127703 : if ( k >= lel(i) .and. k <= mx(i) - 1) then
2510 2726361 : dadt(i) = dadt(i) + dboydt(i,k)* (zf(i,k)-zf(i,k+1))
2511 : endif
2512 : end do
2513 : end do
2514 357782 : do i = il1g,il2g
2515 277765 : dltaa = -1._kind_phys* (cape(i)-capelmt)
2516 357782 : if (dadt(i) /= 0._kind_phys) mb(i) = max(dltaa/tau/dadt(i),0._kind_phys)
2517 : end do
2518 : !
2519 80017 : return
2520 : end subroutine closure
2521 :
2522 80017 : subroutine q1q2_pjr(ncol ,pver ,latice ,&
2523 80017 : dqdt ,dsdt ,q ,qs ,qu , &
2524 80017 : su ,du ,qhat ,shat ,dp , &
2525 80017 : mu ,md ,sd ,qd ,ql , &
2526 80017 : dsubcld ,jt ,mx ,il1g ,il2g , &
2527 : cp ,rl ,msg , &
2528 80017 : dl ,evp ,cu)
2529 :
2530 : implicit none
2531 :
2532 : !-----------------------------------------------------------------------
2533 : ! Purpose:
2534 : ! compute temperature and moisture changes due to convection.
2535 : !-----------------------------------------------------------------------
2536 :
2537 :
2538 : real(kind_phys), intent(in) :: cp
2539 :
2540 : integer, intent(in) :: ncol
2541 : integer, intent(in) :: pver
2542 : real(kind_phys), intent(in) :: latice
2543 : integer, intent(in) :: il1g
2544 : integer, intent(in) :: il2g
2545 : integer, intent(in) :: msg
2546 :
2547 : real(kind_phys), intent(in) :: q(ncol,pver)
2548 : real(kind_phys), intent(in) :: qs(ncol,pver)
2549 : real(kind_phys), intent(in) :: qu(ncol,pver)
2550 : real(kind_phys), intent(in) :: su(ncol,pver)
2551 : real(kind_phys), intent(in) :: du(ncol,pver)
2552 : real(kind_phys), intent(in) :: qhat(ncol,pver)
2553 : real(kind_phys), intent(in) :: shat(ncol,pver)
2554 : real(kind_phys), intent(in) :: dp(ncol,pver)
2555 : real(kind_phys), intent(in) :: mu(ncol,pver)
2556 : real(kind_phys), intent(in) :: md(ncol,pver)
2557 : real(kind_phys), intent(in) :: sd(ncol,pver)
2558 : real(kind_phys), intent(in) :: qd(ncol,pver)
2559 : real(kind_phys), intent(in) :: ql(ncol,pver)
2560 : real(kind_phys), intent(in) :: evp(ncol,pver)
2561 : real(kind_phys), intent(in) :: cu(ncol,pver)
2562 : real(kind_phys), intent(in) :: dsubcld(ncol)
2563 :
2564 : real(kind_phys),intent(out) :: dqdt(ncol,pver),dsdt(ncol,pver)
2565 : real(kind_phys),intent(out) :: dl(ncol,pver)
2566 :
2567 : integer kbm
2568 : integer ktm
2569 : integer jt(ncol)
2570 : integer mx(ncol)
2571 : !
2572 : ! work fields:
2573 : !
2574 : integer i
2575 : integer k
2576 :
2577 : real(kind_phys) emc
2578 : real(kind_phys) rl
2579 : !-------------------------------------------------------------------
2580 2240476 : do k = msg + 1,pver
2581 9740131 : do i = il1g,il2g
2582 7499655 : dsdt(i,k) = 0._kind_phys
2583 7499655 : dqdt(i,k) = 0._kind_phys
2584 9660114 : dl(i,k) = 0._kind_phys
2585 : end do
2586 : end do
2587 :
2588 : !
2589 : ! find the highest level top and bottom levels of convection
2590 : !
2591 80017 : ktm = pver
2592 80017 : kbm = pver
2593 357782 : do i = il1g, il2g
2594 277765 : ktm = min(ktm,jt(i))
2595 357782 : kbm = min(kbm,mx(i))
2596 : end do
2597 :
2598 946498 : do k = ktm,pver-1
2599 4037079 : do i = il1g,il2g
2600 6181162 : emc = -cu (i,k) & ! condensation in updraft
2601 6181162 : +evp(i,k) ! evaporating rain in downdraft
2602 :
2603 : dsdt(i,k) = -rl/cp*emc &
2604 3090581 : + (+mu(i,k+1)* (su(i,k+1)-shat(i,k+1)) &
2605 : -mu(i,k)* (su(i,k)-shat(i,k)) &
2606 : +md(i,k+1)* (sd(i,k+1)-shat(i,k+1)) &
2607 : -md(i,k)* (sd(i,k)-shat(i,k)) &
2608 3090581 : )/dp(i,k)
2609 :
2610 : dqdt(i,k) = emc + &
2611 : (+mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)) &
2612 : -mu(i,k)* (qu(i,k)-qhat(i,k)) &
2613 : +md(i,k+1)* (qd(i,k+1)-qhat(i,k+1)) &
2614 : -md(i,k)* (qd(i,k)-qhat(i,k)) &
2615 3090581 : )/dp(i,k)
2616 :
2617 3957062 : dl(i,k) = du(i,k)*ql(i,k+1)
2618 :
2619 : end do
2620 : end do
2621 :
2622 : !
2623 179335 : do k = kbm,pver
2624 541292 : do i = il1g,il2g
2625 461275 : if (k == mx(i)) then
2626 277765 : dsdt(i,k) = (1._kind_phys/dsubcld(i))* &
2627 : (-mu(i,k)* (su(i,k)-shat(i,k)) &
2628 : -md(i,k)* (sd(i,k)-shat(i,k)) &
2629 277765 : )
2630 : dqdt(i,k) = (1._kind_phys/dsubcld(i))* &
2631 : (-mu(i,k)*(qu(i,k)-qhat(i,k)) &
2632 : -md(i,k)*(qd(i,k)-qhat(i,k)) &
2633 277765 : )
2634 84192 : else if (k > mx(i)) then
2635 20787 : dsdt(i,k) = dsdt(i,k-1)
2636 20787 : dqdt(i,k) = dqdt(i,k-1)
2637 : end if
2638 : end do
2639 : end do
2640 : !
2641 80017 : return
2642 : end subroutine q1q2_pjr
2643 :
2644 :
2645 : ! Wrapper for qsat_water that does translation between Pa and hPa
2646 : ! qsat_water uses Pa internally, so get it right, need to pass in Pa.
2647 : ! Afterward, set es back to hPa.
2648 677814570 : subroutine qsat_hPa(t, p, es, qm)
2649 : use wv_saturation, only: qsat_water
2650 :
2651 : ! Inputs
2652 : real(kind_phys), intent(in) :: t ! Temperature (K)
2653 : real(kind_phys), intent(in) :: p ! Pressure (hPa)
2654 : ! Outputs
2655 : real(kind_phys), intent(out) :: es ! Saturation vapor pressure (hPa)
2656 : real(kind_phys), intent(out) :: qm ! Saturation mass mixing ratio
2657 : ! (vapor mass over dry mass, kg/kg)
2658 :
2659 677814570 : call qsat_water(t, p*100._kind_phys, es, qm)
2660 :
2661 677814570 : es = es*0.01_kind_phys
2662 :
2663 677814570 : end subroutine qsat_hPa
2664 :
2665 : end module zm_convr
|
