CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
sgscloud_radpre.F90
1
8
10
11 contains
12
37 subroutine sgscloud_radpre_run( &
38 im,dt,fhswr,levs, &
39 flag_init,flag_restart, &
40 con_g, con_pi, eps, epsm1, &
41 r_v, cpv, rcp, &
42 xlv, xlf, cp, &
43 do_mynnedmf, &
44 qc, qi, qv, T3D, P3D, exner, &
45 qr, qs, qg, &
46 qci_conv,qlc,qli,ud_mf, &
47 imfdeepcnv, imfdeepcnv_gf, &
48 imfdeepcnv_c3, &
49 imfdeepcnv_sas, &
50 qc_save, qi_save, qs_save, &
51 qc_bl,qi_bl,cldfra_bl, &
52 delp,clouds1,clouds2,clouds3, &
53 clouds4,clouds5, &
54 clouds8,clouds9,slmsk, &
55 nlay, plyr, xlat, dz,de_lgth, &
56 cldsa,mtopa,mbota, &
57 imp_physics, imp_physics_gfdl,&
58 imp_physics_fa, &
59 iovr, &
60 errmsg, errflg )
61
62 use machine , only : kind_phys
64 use radcons, only: qmin ! Minimum values for various calculations
65 use funcphys, only: fpvs ! Function to compute sat. vapor pressure over liq.
66!-------------------------------------------------------------------
67 implicit none
68!-------------------------------------------------------------------
69 ! Interface variables
70 real(kind=kind_phys), intent(in) :: con_g, con_pi, eps, epsm1
71 real(kind=kind_phys), intent(in) :: r_v, cpv, rcp
72 real(kind=kind_phys), intent(in) :: xlv, xlf, cp
73 real(kind=kind_phys), intent(in) :: dt,fhswr
74 real :: xls, xlvcp, xlscp !derived below
75 real(kind=kind_phys) :: gfac
76 integer, intent(in) :: im, levs, imfdeepcnv, imfdeepcnv_gf, &
77 & nlay, imfdeepcnv_sas, imfdeepcnv_c3, imp_physics, &
78 & imp_physics_gfdl, imp_physics_fa
79 logical, intent(in) :: flag_init, flag_restart, do_mynnedmf
80
81 real(kind=kind_phys), dimension(:,:), intent(inout) :: qc, qi
82 real(kind=kind_phys), dimension(:,:), intent(inout) :: qr, qs, qg
83 ! note: qci_conv only allocated if GF is used
84 real(kind=kind_phys), dimension(:,:), intent(inout), optional :: qci_conv
85 real(kind=kind_phys), dimension(:,:), intent(inout) :: qlc, qli !for SAS
86 real(kind=kind_phys), dimension(:,:), intent(in), optional :: ud_mf
87 real(kind=kind_phys), dimension(:,:), intent(in) :: t3d,delp
88 real(kind=kind_phys), dimension(:,:), intent(in) :: qv,p3d,exner
89 real(kind=kind_phys), dimension(:,:), intent(inout) :: &
90 & clouds1,clouds2,clouds3,clouds4,clouds5, &
91 & clouds8,clouds9
92 real(kind=kind_phys), dimension(:,:), intent(inout) :: qc_save, qi_save, qs_save
93 real(kind=kind_phys), dimension(:,:), intent(in), optional :: qc_bl, qi_bl, cldfra_bl
94 real(kind=kind_phys), dimension(:), intent(in) :: slmsk, xlat, de_lgth
95 real(kind=kind_phys), dimension(:,:), intent(in) :: plyr, dz
96 real(kind=kind_phys), dimension(:,:), intent(inout) :: cldsa
97 integer, dimension(:,:), intent(inout) :: mbota, mtopa
98 integer, intent(in) :: iovr
99 character(len=*), intent(out) :: errmsg
100 integer, intent(out) :: errflg
101 ! Local variables
102 ! pressure limits of cloud domain interfaces (low,mid,high) in mb (0.1kPa)
103 real(kind=kind_phys) :: ptop1(im,3+1)
104 real(kind=kind_phys) :: ptopc(3+1,2 )
106 data ptopc / 1050., 650., 400., 0.0, 1050., 750., 500., 0.0 /
107 real(kind=kind_phys), dimension(im,nlay) :: cldcnv
108 real(kind=kind_phys), dimension(im) :: rxlat
109 real(kind=kind_phys) :: tc, tk, liqfrac, iwc, ice_frac, snow_frac
110 integer :: i, k, id
111 ! DH* 20200723 - see comment at the end of this routine around 'gethml'
112 real(kind=kind_phys), dimension(im,nlay) :: alpha_dummy
113 ! *DH
114
115 ! PARAMETERS FOR RANDALL AND XU (1996) CLOUD FRACTION
116 real, parameter :: coef_p = 0.25, coef_gamm = 0.49, coef_alph = 100.
117 real :: rhgrid,h2oliq,qsat,tem1,tem2,clwt,es,onemrh,value
118
119 !Chaboureau and Bechtold (2002 and 2005)
120 real :: a, f, sigq, qmq, qt, xl, th, thl, rsl, cpm, cb_cf
121 real(kind=kind_phys) :: tlk
122
123 !Option to convective cloud fraction
124 integer, parameter :: conv_cf_opt = 0 !0: C-B, 1: X-R
125
126 ! Initialize CCPP error handling variables
127 errmsg = ''
128 errflg = 0
129
130 ! some derived variables from incoming constants:
131 xls=xlv+xlf
132 xlvcp=xlv/cp
133 xlscp=(xlv+xlf)/cp
134
135 !write(0,*)"=============================================="
136 !write(0,*)"in SGSCLoud_RadPre"
137 gfac=1.0e5/con_g
138 if (flag_init .and. (.not. flag_restart)) then
139 !write (0,*) 'Skip this flag_init = ', flag_init
140 ! return
141 ! Need default cloud fraction when MYNN is not used: Resort to
142 ! Xu-Randall (1996).
143 ! cloud fraction =
144 ! {1-exp[-100.0*qc/((1-RH)*qsat)**0.49]}*RH**0.25
145 do k = 1, levs
146 do i = 1, im
147 if ( qi(i,k) > 1e-7 .OR. qc(i,k) > 1e-7 ) then
148 es = min( p3d(i,k), fpvs( t3d(i,k) ) ) ! fpvs and prsl in pa
149 qsat = max( qmin, eps * es / (p3d(i,k) + epsm1*es) )
150 rhgrid = max( 0., min( 1., qv(i,k)/qsat ) )
151 h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k) ! g/kg
152 clwt = 1.0e-6 * (p3d(i,k)*0.00001)
153
154 if (h2oliq > clwt) then
155 onemrh= max( 1.e-10, 1.0-rhgrid )
156 tem1 = min(max((onemrh*qsat)**0.49,0.0001),1.0) !jhan
157 tem1 = 100.0 / tem1
158 value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
159 tem2 = sqrt( sqrt(rhgrid) )
160
161 clouds1(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
162 endif
163
164 endif
165 enddo
166 enddo
167
168 else ! timestep > 1 or restart
169
170 ! Back-up microphysics cloud information:
171 do k = 1, levs
172 do i = 1, im
173 qc_save(i,k) = qc(i,k)
174 qi_save(i,k) = qi(i,k)
175 qs_save(i,k) = qs(i,k)
176 end do
177 end do
178
179 if ( do_mynnedmf ) then
180
181 ! add boundary layer clouds - Note: now the temperature-dependent sorting of
182 ! ice and water subgrid-scale clouds is done inside the MYNN-EDMF
183
184 do k = 1, levs
185 do i = 1, im
186
187 !if (imp_physics == imp_physics_gfdl) then
188 ! ! only complement the GFDL cloud fractions
189 ! if (clouds1(i,k) < 0.01 .and. cldfra_bl(i,k) > 0.01) then
190 ! clouds1(i,k) = cldfra_bl(i,k)
191 ! endif
192 !else
193 clouds1(i,k) = cldfra_bl(i,k)
194 !endif
195
196 if (qc(i,k) < 1.e-6 .and. cldfra_bl(i,k)>0.001) then
197 qc(i,k) = qc_bl(i,k)
198
199 !eff radius cloud water (microns) from Miles et al. (2007)
200 if (nint(slmsk(i)) == 1) then !land
201 if(qc(i,k)>1.e-8)clouds3(i,k)=5.4
202 else
203 if(qc(i,k)>1.e-8)clouds3(i,k)=9.6
204 endif
205
206 !calculate the liquid water path using additional BL clouds
207 clouds2(i,k) = max(0.0, qc(i,k) * gfac * delp(i,k))
208 endif
209
210 tc = t3d(i,k) - 273.15
211 !crudely split frozen species into 50% ice and 50% snow below
212 !~700 mb and decrease snow to zero by ~300 mb
213 snow_frac = min(0.5, max((p3d(i,k)-30000.0),0.0)/140000.0)
214 ice_frac = 1.0 - snow_frac
215 if (qi(i,k) < 1.e-9 .and. cldfra_bl(i,k)>0.001) then
216 qi(i,k) = ice_frac*qi_bl(i,k)
217
218 !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b)
219 clouds5(i,k)=max(173.45 + 2.14*tc, 20.)
220 !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 8b)
221 !iwc = qi(i,k)*1.0e6*rho(i,k)
222 !clouds5(i,k)=MAX(139.7 + 1.76*Tc + 13.49*LOG(iwc), 20.)
223
224 !calculate the ice water path using additional BL clouds
225 clouds4(i,k) = max(0.0, qi(i,k) * gfac * delp(i,k))
226 endif
227
228 if (qs(i,k) < 1.e-9 .and. cldfra_bl(i,k)>0.001) then
229 qs(i,k) = snow_frac*qi_bl(i,k)
230
231 !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b)
232 clouds9(i,k)=max(2.*(173.45 + 2.14*tc), 50.)
233
234 !calculate the snow water path using additional BL clouds
235 clouds8(i,k) = max(0.0, qs(i,k) * gfac * delp(i,k))
236 endif
237 enddo
238 enddo
239
240 elseif (imp_physics /= imp_physics_gfdl) then
241
242 ! Non-MYNN cloud fraction AND non-GFDL microphysics, since both
243 ! have their own cloud fractions. In this case, we resort to
244 ! Xu-Randall (1996).
245 ! cloud fraction =
246 ! {1-exp[-100.0*qc/((1-RH)*qsat)**0.49]}*RH**0.25
247 do k = 1, levs
248 do i = 1, im
249 if ( qi(i,k) > 1e-7 .OR. qc(i,k) > 1e-7 ) then
250
251 es = min( p3d(i,k), fpvs( t3d(i,k) ) ) ! fpvs and prsl in pa
252 qsat = max( qmin, eps * es / (p3d(i,k) + epsm1*es) )
253 rhgrid = max( 0., min( 1., qv(i,k)/qsat ) )
254 h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k) ! g/kg
255 clwt = 1.0e-6 * (p3d(i,k)*0.00001)
256
257 if (h2oliq > clwt) then
258 onemrh= max( 1.e-10, 1.0-rhgrid )
259 tem1 = min(max((onemrh*qsat)**0.49,0.0001),1.0) !jhan
260 tem1 = 100.0 / tem1
261 value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
262 tem2 = sqrt( sqrt(rhgrid) )
263
264 clouds1(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
265 endif
266
267 endif
268 enddo
269 enddo
270
271 endif ! end MYNN or OTHER choice for background clouds fractions
272
273 ! At this point, we have cloud properties for all non-deep convective clouds.
274 ! So now we add the convective clouds:
275
276 if (imfdeepcnv == imfdeepcnv_gf .or. imfdeepcnv == imfdeepcnv_c3) then
277 do k = 1, levs
278 do i = 1, im
279 if ( qci_conv(i,k) > 0. ) then
280 tk = t3d(i,k)
281 tc = tk - 273.15
282
283 !Partition the convective clouds into water & frozen species
284 liqfrac = min(1., max(0., (tk-244.)/29.))
285 qc(i,k) = qc(i,k)+qci_conv(i,k)*liqfrac
286 !split ice & snow 50-50%
287 qi(i,k) = qi(i,k)+0.5*qci_conv(i,k)*(1. - liqfrac)
288 qs(i,k) = qs(i,k)+0.5*qci_conv(i,k)*(1. - liqfrac)
289
290 !eff radius cloud water (microns)
291 if (nint(slmsk(i)) == 1) then !land
292 if(qc(i,k)>1.e-8)clouds3(i,k)=5.4
293 else
294 !from Miles et al.
295 if(qc(i,k)>1.e-8)clouds3(i,k)=9.6
296 endif
297 !from Mishra et al. (2014, JGR Atmos), assume R_sno = 2*R_ice
298 if(qi(i,k)>1.e-8)clouds5(i,k)=max( 173.45 + 2.14*tc , 20.)
299 if(qs(i,k)>1.e-8)clouds9(i,k)=max(2.0*(173.45 + 2.14*tc), 50.)
300
301 if ( conv_cf_opt .eq. 0 ) then
302 !print *,'Chab-Bechtold cloud fraction used'
303 ! clouds1(i,k) = cldfra_bl(i,k)
304
305 !Alternatively, use Chaboureau-Bechtold (CB) convective component
306 !Based on both CB2002 and CB2005.
307 xl = xlv*liqfrac + xls*(1.-liqfrac) ! blended heat capacity
308 tlk = t3d(i,k) - xlvcp/exner(i,k)*qc(i,k) &
309 & - xlscp/exner(i,k)*qi(i,k)! liquid temp
310 ! get saturation water vapor mixing ratio at tl and p
311 es = min( p3d(i,k), fpvs( tlk ) ) ! fpvs and prsl in pa
312 qsat= max( qmin, eps*es / (p3d(i,k) + epsm1*es) )
313 rsl = xl*qsat / (r_v*tlk**2) ! slope of C-C curve at t = tl
314 ! CB02, Eqn. 4
315 qt = qc(i,k) + qi(i,k) + qv(i,k) !total water
316 cpm = cp + qt*cpv ! CB02, sec. 2, para. 1
317 a = 1./(1. + xl*rsl/cpm) ! CB02 variable "a"
318 !Now calculate convective component of the cloud fraction:
319 if (a > 0.0) then
320 f = min(1.0/a, 4.0) ! f is the vertical profile
321 else ! scaling function (CB2005)
322 f = 1.0
323 endif
324 sigq = 1.5e-3 * ud_mf(i,k)/dt * f
325 !sigq = 3.E-3 * ud_mf(i,k)/dt * f
326 sigq = sqrt(sigq**2 + 1e-10) ! combined conv + background components
327 qmq = a * (qt - qsat) ! saturation deficit/excess;
328 ! the numerator of Q1
329 cb_cf= min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.0),0.99)
330 if (qci_conv(i,k) .lt. 1e-9) cb_cf = 0.0
331 if (do_mynnedmf .and. qmq .ge. 0.0) then
332 ! leverage C-B stratus clouds from MYNN in saturated conditions
333 if (cb_cf .gt. 0.0) then
334 clouds1(i,k) = 0.5*(clouds1(i,k) + cb_cf)
335 else
336 !default to MYNN clouds - already specified
337 endif
338 else ! unsaturated
339 clouds1(i,k) = cb_cf
340 endif
341 else
342 !print *,'GF with Xu-Randall cloud fraction'
343 ! Xu-Randall (1996) cloud fraction
344 es = min( p3d(i,k), fpvs( t3d(i,k) ) ) ! fpvs and prsl in pa
345 qsat = max( qmin, eps*es / (p3d(i,k) + epsm1*es) )
346 rhgrid = max( 0., min( 1.00, qv(i,k)/qsat ) )
347 h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k) ! g/kg
348 clwt = 1.0e-6 * (p3d(i,k)*0.00001)
349
350 if (h2oliq > clwt) then
351 onemrh= max( 1.e-10, 1.0-rhgrid )
352 tem1 = min(max((onemrh*qsat)**0.49,0.0001),1.0) !jhan
353 tem1 = 100.0 / tem1
354 value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
355 tem2 = sqrt( sqrt(rhgrid) )
356
357 clouds1(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
358 else
359 clouds1(i,k) = 0.0
360 endif
361 !print*,"XuRandla- cf:",clouds1(i,k)," rh:",rhgrid," qt:",h2oliq
362 !print*,"XuRandlb- clwt:",clwt," qsat:",qsat," p:",p3d(i,k)
363 endif ! end convective cf choice
364 endif ! qci_conv
365 enddo
366 enddo
367
368 elseif (imfdeepcnv == imfdeepcnv_sas) then
369
370 do k = 1, levs
371 do i = 1, im
372 h2oliq = qlc(i,k)+qli(i,k)
373 if ( h2oliq > 0. ) then
374 tk = t3d(i,k)
375 tc = tk - 273.15
376
377 !Partition the convective clouds into water & frozen species
378 liqfrac = min(1., max(0., (tk-244.)/29.))
379
380 qc(i,k) = qc(i,k)+qlc(i,k)
381 !split ice & snow 50-50%
382 qi(i,k) = qi(i,k)+0.5*qli(i,k)
383 qs(i,k) = qs(i,k)+0.5*qli(i,k)
384
385 !eff radius cloud water (microns)
386 if (nint(slmsk(i)) == 1) then !land
387 if(qc(i,k)>1.e-8)clouds3(i,k)=5.4
388 else
389 !from Miles et al.
390 if(qc(i,k)>1.e-8)clouds3(i,k)=9.6
391 endif
392 !from Mishra et al. (2014, JGR Atmos), assume R_sno = 2*R_ice
393 if(qi(i,k)>1.e-8)clouds5(i,k)=max( 173.45 + 2.14*tc , 20.)
394 if(qs(i,k)>1.e-8)clouds9(i,k)=max(2.0*(173.45 + 2.14*tc), 50.)
395
396 if ( conv_cf_opt .eq. 0 ) then
397 !print *,'Chab-Bechtold cloud fraction used'
398 !Alternatively, use Chaboureau-Bechtold (CB) convective component
399 !Based on both CB2002 and CB2005.
400 xl = xlv*liqfrac + xls*(1.-liqfrac) ! blended heat capacity
401 tlk = t3d(i,k) - xlvcp/exner(i,k)*qc(i,k) &
402 & - xlscp/exner(i,k)*qi(i,k)! liquid temp
403 ! get saturation water vapor mixing ratio at tl and p
404 es = min( p3d(i,k), fpvs( tlk ) ) ! fpvs and prsl in pa
405 qsat= max( qmin, eps*es / (p3d(i,k) + epsm1*es) )
406 rsl = xl*qsat / (r_v*tlk**2) ! slope of C-C curve at t = tl
407 ! CB02, Eqn. 4
408 qt = qc(i,k) + qi(i,k) + qv(i,k) !total water
409 cpm = cp + qt*cpv ! CB02, sec. 2, para. 1
410 a = 1./(1. + xl*rsl/cpm) ! CB02 variable "a"
411 !Now calculate convective component of the cloud fraction:
412 if (a > 0.0) then
413 f = min(1.0/a, 4.0) ! f is the vertical profile
414 else ! scaling function (CB2005)
415 f = 1.0
416 endif
417 sigq = 1.5e-3 * ud_mf(i,k)/dt * f
418 !sigq = 3.E-3 * ud_mf(i,k)/dt * f
419 sigq = sqrt(sigq**2 + 1e-10) ! combined conv + background components
420 qmq = a * (qt - qsat) ! saturation deficit/excess;
421 ! the numerator of Q1
422 cb_cf= min(max(0.5 + 0.36 * atan(1.55*(qmq/sigq)),0.0),0.99)
423 if (h2oliq .lt. 1e-9) cb_cf = 0.0
424 if (do_mynnedmf .and. qmq .ge. 0.0) then
425 ! leverage C-B stratus clouds from MYNN in saturated conditions
426 if (cb_cf .gt. 0.0) then
427 clouds1(i,k) = 0.5*(clouds1(i,k) + cb_cf)
428 else
429 !default to MYNN clouds - already specified
430 endif
431 else ! unsaturated
432 clouds1(i,k) = cb_cf
433 endif
434 else
435 !print *,'SAS with Xu-Randall cloud fraction'
436 ! Xu-Randall (1996) cloud fraction
437 es = min( p3d(i,k), fpvs( t3d(i,k) ) ) ! fpvs and prsl in pa
438 qsat = max( qmin, eps*es / (p3d(i,k) + epsm1*es) )
439 rhgrid = max( 0., min( 1.00, qv(i,k)/qsat ) )
440 h2oliq = qc(i,k) + qi(i,k) + qr(i,k) + qs(i,k) + qg(i,k) ! g/kg
441 clwt = 1.0e-6 * (p3d(i,k)*0.00001)
442
443 if (h2oliq > clwt) then
444 onemrh= max( 1.e-10, 1.0-rhgrid )
445 tem1 = min(max((onemrh*qsat)**0.49,0.0001),1.0) !jhan
446 tem1 = 100.0 / tem1
447 value = max( min( tem1*(h2oliq-clwt), 50.0 ), 0.0 )
448 tem2 = sqrt( sqrt(rhgrid) )
449
450 clouds1(i,k) = max( tem2*(1.0-exp(-value)), 0.0 )
451 else
452 clouds1(i,k) = 0.0
453 endif
454 !print*,"XuRandla- cf:",clouds1(i,k)," rh:",rhgrid," qt:",h2oliq
455 !print*,"XuRandlb- clwt:",clwt," qsat:",qsat," p:",p3d(i,k)
456 endif ! end convective cf choice
457 endif ! qlc/qli check
458 enddo
459 enddo
460
461 endif ! convection scheme check
462
463 endif ! timestep > 1
464
465 end subroutine sgscloud_radpre_run
466 end module sgscloud_radpre
real(kind=kind_phys), dimension(nk_clds+1, 2), save ptopc
pressure limits of cloud domain interfaces (low, mid, high) in mb (0.1kPa)
real(kind=kind_phys) gfac
subroutine, public gethml(plyr, ptop1, cldtot, cldcnv, dz, de_lgth, alpha, ix, nlay, iovr, iovr_rand, iovr_maxrand, iovr_max, iovr_dcorr, iovr_exp, iovr_exprand, top_at_1, si, clds, mtop, mbot)
This subroutine computes high, mid, low, total, and boundary cloud fractions and cloud top/bottom lay...
This module provides an Application Program Interface (API) for computing basic thermodynamic physics...
Definition funcphys.f90:26
This module computes cloud related quantities for radiation computations.
This module contains some of the most frequently used math and physics constants for RRTMG.
Definition radcons.f90:7