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