Line data Source code
1 : module clubb_mf
2 :
3 : ! =============================================================================== !
4 : ! Mass-flux module for use with CLUBB !
5 : ! Together (CLUBB+MF) they comprise a eddy-diffusivity mass-flux approach (EDMF) !
6 : ! =============================================================================== !
7 :
8 : use shr_kind_mod, only: r8=>shr_kind_r8
9 : use spmd_utils, only: masterproc
10 : use cam_logfile, only: iulog
11 :
12 : implicit none
13 : private
14 : save
15 :
16 : public :: integrate_mf, &
17 : clubb_mf_readnl, &
18 : do_clubb_mf, &
19 : do_clubb_mf_diag
20 :
21 : real(r8) :: clubb_mf_L0 = 0._r8
22 : real(r8) :: clubb_mf_ent0 = 0._r8
23 : integer :: clubb_mf_nup = 0
24 : logical, protected :: do_clubb_mf = .false.
25 : logical, protected :: do_clubb_mf_diag = .false.
26 :
27 : contains
28 :
29 1536 : subroutine clubb_mf_readnl(nlfile)
30 :
31 : ! =============================================================================== !
32 : ! MF namelists !
33 : ! =============================================================================== !
34 :
35 : use namelist_utils, only: find_group_name
36 : use cam_abortutils, only: endrun
37 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_real8, mpi_integer, mpi_logical
38 :
39 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
40 :
41 : character(len=*), parameter :: sub = 'clubb_mf_readnl'
42 :
43 : integer :: iunit, read_status, ierr
44 :
45 :
46 : namelist /clubb_mf_nl/ clubb_mf_L0, clubb_mf_ent0, clubb_mf_nup, do_clubb_mf, do_clubb_mf_diag
47 :
48 1536 : if (masterproc) then
49 2 : open( newunit=iunit, file=trim(nlfile), status='old' )
50 2 : call find_group_name(iunit, 'clubb_mf_nl', status=read_status)
51 2 : if (read_status == 0) then
52 2 : read(iunit, clubb_mf_nl, iostat=ierr)
53 2 : if (ierr /= 0) then
54 0 : call endrun('clubb_mf_readnl: ERROR reading namelist')
55 : end if
56 : end if
57 2 : close(iunit)
58 : end if
59 :
60 1536 : call mpi_bcast(clubb_mf_L0, 1, mpi_real8, mstrid, mpicom, ierr)
61 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_mf_L0")
62 1536 : call mpi_bcast(clubb_mf_ent0, 1, mpi_real8, mstrid, mpicom, ierr)
63 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_mf_ent0")
64 1536 : call mpi_bcast(clubb_mf_nup, 1, mpi_integer, mstrid, mpicom, ierr)
65 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: clubb_mf_nup")
66 1536 : call mpi_bcast(do_clubb_mf, 1, mpi_logical, mstrid, mpicom, ierr)
67 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_clubb_mf")
68 1536 : call mpi_bcast(do_clubb_mf_diag, 1, mpi_logical, mstrid, mpicom, ierr)
69 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: do_clubb_mf_diag")
70 :
71 1536 : if ((.not. do_clubb_mf) .and. do_clubb_mf_diag ) then
72 0 : call endrun('clubb_mf_readnl: Error - cannot turn on do_clubb_mf_diag without also turning on do_clubb_mf')
73 : end if
74 :
75 :
76 1536 : end subroutine clubb_mf_readnl
77 :
78 0 : subroutine integrate_mf( nz, dzt, zm, p_zm, iexner_zm, & ! input
79 0 : p_zt, iexner_zt, & ! input
80 0 : u, v, thl, qt, thv, & ! input
81 0 : thl_zm, qt_zm, & ! input
82 : wthl, wqt, pblh, & ! input
83 0 : dry_a, moist_a, & ! output - plume diagnostics
84 0 : dry_w, moist_w, & ! output - plume diagnostics
85 0 : dry_qt, moist_qt, & ! output - plume diagnostics
86 0 : dry_thl, moist_thl, & ! output - plume diagnostics
87 0 : dry_u, moist_u, & ! output - plume diagnostics
88 0 : dry_v, moist_v, & ! output - plume diagnostics
89 0 : moist_qc, & ! output - plume diagnostics
90 0 : ae, aw, & ! output - diagnosed fluxes BEFORE mean field update
91 0 : awthl, awqt, & ! output - diagnosed fluxes BEFORE mean field update
92 0 : awql, awqi, & ! output - diagnosed fluxes BEFORE mean field update
93 0 : awu, awv, & ! output - diagnosed fluxes BEFORE mean field update
94 0 : thlflx, qtflx ) ! output - variables needed for solver
95 :
96 : ! ================================================================================= !
97 : ! Mass-flux algorithm !
98 : ! !
99 : ! Provides rtm and thl fluxes due to mass flux ensemble, !
100 : ! which are fed into the mixed explicit/implicit clubb solver as explicit terms !
101 : ! !
102 : ! Variables needed for solver !
103 : ! ae = sum_i (1-a_i) !
104 : ! aw3 = sum (a_i w_i) !
105 : ! aws3 = sum (a_i w_i*s_i); s=thl*cp !
106 : ! aws3,awqv3,awql3,awqi3,awu3,awv3 similar as above except for different variables !
107 : ! !
108 : ! Mass flux variables are computed on edges (i.e. momentum grid): !
109 : ! upa,upw,upqt,... !
110 : ! dry_a,moist_a,dry_w,moist_w, ... !
111 : ! !
112 : ! In CLUBB (unlike CAM) nlevs of momentum grid = nlevs of thermodynamic grid, !
113 : ! due to a subsurface thermodynamic layer. To avoid confusion, below the variables !
114 : ! are grouped by the grid they are on. !
115 : ! !
116 : ! *note that state on the lowest thermo level is equal to state on the lowest !
117 : ! momentum level due to state_zt(1) = state_zt(2), and lowest momentum level !
118 : ! is a weighted combination of the lowest two thermodynamic levels. !
119 : ! !
120 : ! ---------------------------------Authors---------------------------------------- !
121 : ! Marcin Kurowski, JPL !
122 : ! Modified heavily by Mikael Witte, UCLA/JPL for implementation in CESM2/E3SM !
123 : ! Additional modifications by Adam Herrington, NCAR !
124 : ! ================================================================================= !
125 :
126 : use physconst, only: rair, cpair, gravit, latvap, latice, zvir
127 :
128 : integer, intent(in) :: nz
129 : real(r8), dimension(nz), intent(in) :: u, v, & ! thermodynamic grid
130 : thl, thv, & ! thermodynamic grid
131 : qt, & ! thermodynamic grid
132 : dzt, & ! thermodynamic grid
133 : p_zt, iexner_zt, & ! thermodynamic grid
134 : thl_zm, & ! momentum grid
135 : qt_zm, & ! momentum grid
136 : zm, & ! momentum grid
137 : p_zm, iexner_zm ! momentum grid
138 :
139 : real(r8), intent(in) :: wthl,wqt
140 : real(r8), intent(inout) :: pblh
141 :
142 : real(r8),dimension(nz), intent(out) :: dry_a, moist_a, & ! momentum grid
143 : dry_w, moist_w, & ! momentum grid
144 : dry_qt, moist_qt, & ! momentum grid
145 : dry_thl, moist_thl, & ! momentum grid
146 : dry_u, moist_u, & ! momentum grid
147 : dry_v, moist_v, & ! momentum grid
148 : moist_qc, & ! momentum grid
149 : ae, aw, & ! momentum grid
150 : awthl, awqt, & ! momentum grid
151 : awql, awqi, & ! momentum grid
152 : awu, awv, & ! momentum grid
153 : thlflx, qtflx ! momentum grid
154 :
155 : ! =============================================================================== !
156 : ! INTERNAL VARIABLES
157 : !
158 : ! sums over all plumes
159 : real(r8), dimension(nz) :: moist_th, dry_th, & ! momentum grid
160 0 : awqv, awth ! momentum grid
161 : !
162 : ! updraft properties
163 0 : real(r8), dimension(nz,clubb_mf_nup) :: upw, upa, & ! momentum grid
164 0 : upqt, upqc, & ! momentum grid
165 0 : upqv, upqs, & ! momentum grid
166 0 : upql, upqi, & ! momentum grid
167 0 : upth, upthv, & ! momentum grid
168 0 : upthl, & ! momentum grid
169 0 : upu, upv ! momentum grid
170 : !
171 : ! entrainment profiles
172 0 : real(r8), dimension(nz,clubb_mf_nup) :: ent, entf ! thermodynamic grid
173 0 : integer, dimension(nz,clubb_mf_nup) :: enti ! thermodynamic grid
174 : !
175 : ! other variables
176 : integer :: k,i,ih
177 0 : real(r8), dimension(clubb_mf_nup) :: zcb
178 : real(r8) :: zcb_unset, &
179 : wthv, &
180 : wstar, qstar, thvstar, &
181 : sigmaw, sigmaqt, sigmathv,&
182 : wmin, wmax, &
183 : wlv, wtv, wp, &
184 : B, & ! thermodynamic grid
185 : entexp, entexpu, entw, & ! thermodynamic grid
186 : thln, thvn, thn, & ! momentum grid
187 : qtn, qsn, & ! momentum grid
188 : qcn, qln, qin, & ! momentum grid
189 : un, vn, wn2, & ! momentum grid
190 : lmixn, & ! momentum grid
191 : supqt, supthl ! thermodynamic grid
192 : !
193 : ! parameters defining initial conditions for updrafts
194 : real(r8),parameter :: pwmin = 1.5_r8, &
195 : pwmax = 3._r8
196 :
197 : !
198 : ! alpha, z-scores after Suselj etal 2019
199 : real(r8),parameter :: alphw = 0.572_r8, &
200 : alphqt = 2.890_r8, &
201 : alphthv = 2.890_r8
202 : !
203 : ! w' covariance after Suselj etal 2019
204 : real(r8),parameter :: cwqt = 0.32_r8, &
205 : cwthv = 0.58_r8
206 : !
207 : ! virtual mass coefficients for w-eqn after Suselj etal 2019
208 : real(r8),parameter :: wa = 1.0_r8, &
209 : wb = 1.5_r8
210 : !
211 : ! min values to avoid singularities
212 : real(r8),parameter :: wstarmin = 1.e-3_r8, &
213 : pblhmin = 100._r8
214 : !
215 : ! to condensate or not to condensate
216 : logical :: do_condensation = .true.
217 : !
218 : ! to precip or not to precip
219 : logical :: do_precip = .false.
220 : !
221 : ! to debug flag (overides stochastic entrainment)
222 : logical :: debug = .false.
223 : real(r8),parameter :: fixent = 0.004_r8
224 :
225 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
226 : !!!!!!!!!!!!!!!!!!!!!! BEGIN CODE !!!!!!!!!!!!!!!!!!!!!!!
227 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
228 :
229 : ! INITIALIZE OUTPUT VARIABLES
230 : ! set updraft properties to zero
231 0 : dry_a = 0._r8
232 0 : moist_a = 0._r8
233 0 : dry_w = 0._r8
234 0 : moist_w = 0._r8
235 0 : dry_qt = 0._r8
236 0 : moist_qt = 0._r8
237 0 : dry_thl = 0._r8
238 0 : moist_thl = 0._r8
239 0 : dry_u = 0._r8
240 0 : moist_u = 0._r8
241 0 : dry_v = 0._r8
242 0 : moist_v = 0._r8
243 0 : moist_qc = 0._r8
244 : ! outputs - variables needed for solver
245 0 : aw = 0._r8
246 0 : awthl = 0._r8
247 0 : awqt = 0._r8
248 0 : awqv = 0._r8
249 0 : awql = 0._r8
250 0 : awqi = 0._r8
251 0 : awu = 0._r8
252 0 : awv = 0._r8
253 0 : thlflx = 0._r8
254 0 : qtflx = 0._r8
255 :
256 0 : ent = 0._r8
257 0 : entf = 0._r8
258 0 : enti = 0
259 :
260 : ! this is the environmental area - by default 1.
261 0 : ae = 1._r8
262 :
263 : ! START MAIN COMPUTATION
264 0 : upw = 0._r8
265 0 : upthl = 0._r8
266 0 : upthv = 0._r8
267 0 : upqt = 0._r8
268 0 : upa = 0._r8
269 0 : upu = 0._r8
270 0 : upv = 0._r8
271 0 : upqc = 0._r8
272 0 : upth = 0._r8
273 0 : upql = 0._r8
274 0 : upqi = 0._r8
275 0 : upqv = 0._r8
276 0 : upqs = 0._r8
277 :
278 : ! unique identifier
279 : zcb_unset = 9999999._r8
280 0 : zcb = zcb_unset
281 :
282 0 : pblh = max(pblh,pblhmin)
283 0 : wthv = wthl+zvir*thv(1)*wqt
284 :
285 : ! if surface buoyancy is positive then do mass-flux
286 0 : if ( wthv > 0._r8 ) then
287 :
288 0 : if (debug) then
289 : ! overide stochastic entrainment with fixent
290 0 : ent(:,:) = fixent
291 : else
292 :
293 : ! get entrainment coefficient, dz/L0
294 0 : do i=1,clubb_mf_nup
295 0 : do k=1,nz
296 0 : entf(k,i) = dzt(k) / clubb_mf_L0
297 : enddo
298 : enddo
299 :
300 : ! get poisson, P(dz/L0)
301 0 : call poisson( nz, clubb_mf_nup, entf, enti, u(2:5))
302 :
303 : ! get entrainment, ent=ent0/dz*P(dz/L0)
304 0 : do i=1,clubb_mf_nup
305 0 : do k=1,nz
306 0 : ent(k,i) = real( enti(k,i))*clubb_mf_ent0/dzt(k)
307 : enddo
308 : enddo
309 :
310 : end if
311 :
312 : ! get surface conditions
313 0 : wstar = max( wstarmin, (gravit/thv(1)*wthv*pblh)**(1._r8/3._r8) )
314 0 : qstar = wqt / wstar
315 0 : thvstar = wthv / wstar
316 :
317 0 : sigmaw = alphw * wstar
318 0 : sigmaqt = alphqt * abs(qstar)
319 0 : sigmathv = alphthv * abs(thvstar)
320 :
321 0 : wmin = sigmaw * pwmin
322 0 : wmax = sigmaw * pwmax
323 :
324 0 : do i=1,clubb_mf_nup
325 :
326 0 : wlv = wmin + (wmax-wmin) / (real(clubb_mf_nup,r8)) * (real(i-1, r8))
327 0 : wtv = wmin + (wmax-wmin) / (real(clubb_mf_nup,r8)) * real(i,r8)
328 :
329 0 : upw(1,i) = 0.5_r8 * (wlv+wtv)
330 : upa(1,i) = 0.5_r8 * erf( wtv/(sqrt(2._r8)*sigmaw) ) &
331 0 : - 0.5_r8 * erf( wlv/(sqrt(2._r8)*sigmaw) )
332 :
333 0 : upu(1,i) = u(1)
334 0 : upv(1,i) = v(1)
335 :
336 0 : upqt(1,i) = qt(1) + cwqt * upw(1,i) * sigmaqt/sigmaw
337 0 : upthv(1,i) = thv(1) + cwthv * upw(1,i) * sigmathv/sigmaw
338 0 : upthl(1,i) = upthv(1,i) / (1._r8+zvir*upqt(1,i))
339 :
340 : ! get cloud, lowest momentum level
341 0 : if (do_condensation) then
342 : call condensation_mf(upqt(1,i), upthl(1,i), p_zm(1), iexner_zm(1), &
343 0 : thvn, qcn, thn, qln, qin, qsn, lmixn)
344 0 : upthv(1,i) = thvn
345 0 : upqc(1,i) = qcn
346 0 : upql(1,i) = qln
347 0 : upqi(1,i) = qin
348 0 : upqs(1,i) = qsn
349 0 : if (qcn > 0._r8) zcb(i) = zm(1)
350 : else
351 : ! assume no cldliq
352 0 : upqc(1,i) = 0._r8
353 : end if
354 :
355 : enddo
356 :
357 : ! get updraft properties
358 0 : do i=1,clubb_mf_nup
359 0 : do k=1,nz-1
360 :
361 : ! get microphysics, autoconversion
362 0 : if (do_precip .and. upqc(k,i) > 0._r8) then
363 0 : call precip_mf(upqs(k,i),upqt(k,i),upw(k,i),dzt(k+1),zm(k+1)-zcb(i),supqt)
364 :
365 0 : supthl = -1._r8*lmixn*supqt*iexner_zt(k+1)/cpair
366 : else
367 0 : supqt = 0._r8
368 0 : supthl = 0._r8
369 : end if
370 :
371 : ! integrate updraft
372 0 : entexp = exp(-ent(k+1,i)*dzt(k+1))
373 0 : entexpu = exp(-ent(k+1,i)*dzt(k+1)/3._r8)
374 :
375 0 : qtn = qt(k+1) *(1._r8-entexp ) + upqt (k,i)*entexp + supqt
376 0 : thln = thl(k+1)*(1._r8-entexp ) + upthl(k,i)*entexp + supthl
377 0 : un = u(k+1) *(1._r8-entexpu) + upu (k,i)*entexpu
378 0 : vn = v(k+1) *(1._r8-entexpu) + upv (k,i)*entexpu
379 :
380 : ! get cloud, momentum levels
381 0 : if (do_condensation) then
382 : call condensation_mf(qtn, thln, p_zm(k+1), iexner_zm(k+1), &
383 0 : thvn, qcn, thn, qln, qin, qsn, lmixn)
384 0 : if (zcb(i).eq.zcb_unset .and. qcn > 0._r8) zcb(i) = zm(k+1)
385 : else
386 0 : thvn = thln*(1._r8+zvir*qtn)
387 : end if
388 :
389 : ! get buoyancy
390 0 : B=gravit*(0.5_r8*(thvn + upthv(k,i))/thv(k+1)-1._r8)
391 0 : if (debug) then
392 0 : if ( masterproc ) then
393 0 : write(iulog,*) "B(k,i), k, i ", B, k, i
394 : end if
395 : end if
396 :
397 : ! get wn^2
398 0 : wp = wb*ent(k+1,i)
399 0 : if (wp==0._r8) then
400 0 : wn2 = upw(k,i)**2._r8+2._r8*wa*B*dzt(k+1)
401 : else
402 0 : entw = exp(-2._r8*wp*dzt(k+1))
403 0 : wn2 = entw*upw(k,i)**2._r8+wa*B/(wb*ent(k+1,i))*(1._r8-entw)
404 : end if
405 :
406 0 : if (wn2>0._r8) then
407 0 : upw(k+1,i) = sqrt(wn2)
408 0 : upthv(k+1,i) = thvn
409 0 : upthl(k+1,i) = thln
410 0 : upqt(k+1,i) = qtn
411 0 : upqc(k+1,i) = qcn
412 0 : upqs(k+1,i) = qsn
413 0 : upu(k+1,i) = un
414 0 : upv(k+1,i) = vn
415 0 : upa(k+1,i) = upa(k,i)
416 0 : upql(k+1,i) = qln
417 0 : upqi(k+1,i) = qin
418 0 : upqv(k+1,i) = qtn - qcn
419 : else
420 : exit
421 : end if
422 : enddo
423 : enddo
424 :
425 : ! writing updraft properties for output
426 0 : do k=1,nz
427 :
428 : ! first sum over all i-updrafts
429 0 : do i=1,clubb_mf_nup
430 0 : if (upqc(k,i)>0._r8) then
431 0 : moist_a(k) = moist_a(k) + upa(k,i)
432 0 : moist_w(k) = moist_w(k) + upa(k,i)*upw(k,i)
433 0 : moist_qt(k) = moist_qt(k) + upa(k,i)*upqt(k,i)
434 0 : moist_thl(k) = moist_thl(k) + upa(k,i)*upthl(k,i)
435 0 : moist_u(k) = moist_u(k) + upa(k,i)*upu(k,i)
436 0 : moist_v(k) = moist_v(k) + upa(k,i)*upv(k,i)
437 0 : moist_qc(k) = moist_qc(k) + upa(k,i)*upqc(k,i)
438 : else
439 0 : dry_a(k) = dry_a(k) + upa(k,i)
440 0 : dry_w(k) = dry_w(k) + upa(k,i)*upw(k,i)
441 0 : dry_qt(k) = dry_qt(k) + upa(k,i)*upqt(k,i)
442 0 : dry_thl(k) = dry_thl(k) + upa(k,i)*upthl(k,i)
443 0 : dry_u(k) = dry_u(k) + upa(k,i)*upu(k,i)
444 0 : dry_v(k) = dry_v(k) + upa(k,i)*upv(k,i)
445 : endif
446 : enddo
447 :
448 0 : if ( dry_a(k) > 0._r8 ) then
449 0 : dry_w(k) = dry_w(k) / dry_a(k)
450 0 : dry_qt(k) = dry_qt(k) / dry_a(k)
451 0 : dry_thl(k) = dry_thl(k) / dry_a(k)
452 0 : dry_u(k) = dry_u(k) / dry_a(k)
453 0 : dry_v(k) = dry_v(k) / dry_a(k)
454 : else
455 0 : dry_w(k) = 0._r8
456 0 : dry_qt(k) = 0._r8
457 0 : dry_thl(k) = 0._r8
458 0 : dry_u(k) = 0._r8
459 0 : dry_v(k) = 0._r8
460 : endif
461 :
462 0 : if ( moist_a(k) > 0._r8 ) then
463 0 : moist_w(k) = moist_w(k) / moist_a(k)
464 0 : moist_qt(k) = moist_qt(k) / moist_a(k)
465 0 : moist_thl(k) = moist_thl(k) / moist_a(k)
466 0 : moist_u(k) = moist_u(k) / moist_a(k)
467 0 : moist_v(k) = moist_v(k) / moist_a(k)
468 0 : moist_qc(k) = moist_qc(k) / moist_a(k)
469 : else
470 0 : moist_w(k) = 0._r8
471 0 : moist_qt(k) = 0._r8
472 0 : moist_thl(k) = 0._r8
473 0 : moist_u(k) = 0._r8
474 0 : moist_v(k) = 0._r8
475 0 : moist_qc(k) = 0._r8
476 : endif
477 :
478 : enddo
479 :
480 0 : do k=1,nz
481 0 : do i=1,clubb_mf_nup
482 0 : ae (k) = ae (k) - upa(k,i)
483 0 : aw (k) = aw (k) + upa(k,i)*upw(k,i)
484 0 : awu (k) = awu (k) + upa(k,i)*upw(k,i)*upu(k,i)
485 0 : awv (k) = awv (k) + upa(k,i)*upw(k,i)*upv(k,i)
486 0 : awthl(k)= awthl(k)+ upa(k,i)*upw(k,i)*upthl(k,i)
487 0 : awqt(k) = awqt(k) + upa(k,i)*upw(k,i)*upqt(k,i)
488 0 : awqv(k) = awqv(k) + upa(k,i)*upw(k,i)*upqv(k,i)
489 0 : awql(k) = awql(k) + upa(k,i)*upw(k,i)*upql(k,i)
490 0 : awqi(k) = awqi(k) + upa(k,i)*upw(k,i)*upqi(k,i)
491 : enddo
492 : enddo
493 :
494 0 : do k=1,nz
495 0 : thlflx(k)= awthl(k) - aw(k)*thl_zm(k)
496 0 : qtflx(k)= awqt(k) - aw(k)*qt_zm(k)
497 : enddo
498 :
499 : end if ! ( wthv > 0.0 )
500 :
501 0 : end subroutine integrate_mf
502 :
503 0 : subroutine condensation_mf( qt, thl, p, iex, thv, qc, th, ql, qi, qs, lmix )
504 : ! =============================================================================== !
505 : ! zero or one condensation for edmf: calculates thv and qc !
506 : ! =============================================================================== !
507 : use physconst, only: cpair, zvir, h2otrip
508 : use wv_saturation, only : qsat
509 :
510 : real(r8),intent(in) :: qt,thl,p,iex
511 : real(r8),intent(out):: thv,qc,th,ql,qi,qs,lmix
512 :
513 : !local variables
514 : integer :: niter,i
515 : real(r8) :: diff,t,qstmp,qcold,es,wf
516 :
517 : ! max number of iterations
518 0 : niter=50
519 : ! minimum difference
520 0 : diff=2.e-5_r8
521 :
522 0 : qc=0._r8
523 0 : t=thl/iex
524 :
525 : !by definition:
526 : ! T = Th*Exner, Exner=(p/p0)^(R/cp) (1)
527 : ! Thl = Th - L/cp*ql/Exner (2)
528 : !so:
529 : ! Th = Thl + L/cp*ql/Exner (3)
530 : ! T = Th*Exner=(Thl+L/cp*ql/Exner)*Exner (4)
531 : ! = Thl*Exner + L/cp*ql
532 0 : do i=1,niter
533 0 : wf = get_watf(t)
534 0 : t = thl/iex+get_alhl(wf)/cpair*qc !as in (4)
535 :
536 : ! qsat, p is in pascal (check!)
537 0 : call qsat(t,p,es,qstmp)
538 0 : qcold = qc
539 0 : qc = max(0.5_r8*qc+0.5_r8*(qt-qstmp),0._r8)
540 0 : if (abs(qc-qcold)<diff) exit
541 : enddo
542 :
543 0 : wf = get_watf(t)
544 0 : t = thl/iex+get_alhl(wf)/cpair*qc
545 :
546 0 : call qsat(t,p,es,qs)
547 0 : qc = max(qt-qs,0._r8)
548 0 : thv = (thl+get_alhl(wf)/cpair*iex*qc)*(1._r8+zvir*(qt-qc)-qc)
549 0 : lmix = get_alhl(wf)
550 0 : th = t*iex
551 0 : qi = qc*(1._r8-wf)
552 0 : ql = qc*wf
553 :
554 : contains
555 :
556 0 : function get_watf(t)
557 : real(r8) :: t,get_watf,tc
558 : real(r8), parameter :: &
559 : tmax=-10._r8, &
560 : tmin=-40._r8
561 :
562 0 : tc=t-h2otrip
563 :
564 0 : if (tc>tmax) then
565 : get_watf=1._r8
566 0 : else if (tc<tmin) then
567 : get_watf=0._r8
568 : else
569 0 : get_watf=(tc-tmin)/(tmax-tmin);
570 : end if
571 :
572 0 : end function get_watf
573 :
574 :
575 0 : function get_alhl(wf)
576 : !latent heat of the mixture based on water fraction
577 : use physconst, only : latvap , latice
578 : real(r8) :: get_alhl,wf
579 :
580 0 : get_alhl = wf*latvap+(1._r8-wf)*(latvap+latice)
581 :
582 0 : end function get_alhl
583 :
584 : end subroutine condensation_mf
585 :
586 0 : subroutine precip_mf(qs,qt,w,dz,dzcld,Supqt)
587 : !**********************************************************************
588 : ! Precipitation microphysics
589 : ! By Adam Herrington, after Kay Suselj
590 : !**********************************************************************
591 :
592 : real(r8),intent(in) :: qs,qt,w,dz,dzcld
593 : real(r8),intent(out) :: Supqt
594 : !
595 : ! local vars
596 : real(r8) :: tauwgt, tau, & ! time-scale vars
597 : qstar ! excess cloud liquid
598 :
599 : real(r8),parameter :: tau0 = 15._r8, & ! base time-scale
600 : zmin = 300._r8, & ! small cloud thick
601 : zmax = 3000._r8, & ! large cloud thick
602 : qcmin = 0.00125_r8 ! supersat threshold
603 :
604 0 : qstar = qs+qcmin
605 :
606 0 : if (qt > qstar) then
607 : ! get precip efficiency
608 0 : tauwgt = (dzcld-zmin)/(zmax-zmin)
609 0 : tauwgt = min(max(tauwgt,0._r8),1._r8)
610 0 : tau = tauwgt/tau0
611 :
612 : ! get source for updraft
613 0 : Supqt = (qstar-qt)*(1._r8 - exp(-1._r8*tau*dz/w))
614 : else
615 0 : Supqt = 0._r8
616 : end if
617 :
618 0 : end subroutine precip_mf
619 :
620 0 : subroutine poisson(nz,nup,lambda,poi,state)
621 : !**********************************************************************
622 : ! Set a unique (but reproduceble) seed for the kiss RNG
623 : ! Call Poisson deviate
624 : ! By Adam Herrington
625 : !**********************************************************************
626 : use shr_RandNum_mod, only: ShrKissRandGen
627 :
628 : integer, intent(in) :: nz,nup
629 : real(r8), dimension(4), intent(in) :: state
630 : real(r8), dimension(nz,nup), intent(in) :: lambda
631 : integer, dimension(nz,nup), intent(out) :: poi
632 : integer, dimension(1,4) :: tmpseed
633 : integer :: i,j
634 0 : type(ShrKissRandGen) :: kiss_gen
635 :
636 : ! Compute seed
637 0 : tmpseed(1,1) = int((state(1) - int(state(1))) * 1000000000._r8)
638 0 : tmpseed(1,2) = int((state(2) - int(state(2))) * 1000000000._r8)
639 0 : tmpseed(1,3) = int((state(3) - int(state(3))) * 1000000000._r8)
640 0 : tmpseed(1,4) = int((state(4) - int(state(4))) * 1000000000._r8)
641 :
642 : ! Set seed
643 0 : kiss_gen = ShrKissRandGen(tmpseed)
644 :
645 0 : do i=1,nz
646 0 : do j=1,nup
647 0 : call knuth(kiss_gen,lambda(i,j),poi(i,j))
648 : enddo
649 : enddo
650 :
651 0 : end subroutine poisson
652 :
653 0 : subroutine knuth(kiss_gen,lambda,kout)
654 : !**********************************************************************
655 : ! Discrete random poisson from Knuth
656 : ! The Art of Computer Programming, v2, 137-138
657 : ! By Adam Herrington
658 : !**********************************************************************
659 : use shr_RandNum_mod, only: ShrKissRandGen
660 :
661 : type(ShrKissRandGen), intent(inout) :: kiss_gen
662 : real(r8), intent(in) :: lambda
663 : integer, intent(out) :: kout
664 :
665 : ! Local variables
666 : real(r8), dimension(1,1) :: tmpuni
667 : real(r8) :: puni, explam
668 : integer :: k
669 :
670 0 : k = 0
671 0 : explam = exp(-1._r8*lambda)
672 0 : puni = 1._r8
673 0 : do while (puni > explam)
674 0 : k = k + 1
675 0 : call kiss_gen%random(tmpuni)
676 0 : puni = puni*tmpuni(1,1)
677 : end do
678 0 : kout = k - 1
679 :
680 0 : end subroutine knuth
681 :
682 : end module clubb_mf
|
