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