Line data Source code
1 : module cloud_fraction
2 :
3 : ! Cloud fraction parameterization.
4 :
5 :
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use ppgrid, only: pcols, pver, pverp
8 : use ref_pres, only: pref_mid
9 : use spmd_utils, only: masterproc
10 : use cam_logfile, only: iulog
11 : use cam_abortutils, only: endrun
12 : use ref_pres, only: trop_cloud_top_lev
13 :
14 : implicit none
15 : private
16 : save
17 :
18 : ! Public interfaces
19 : public &
20 : cldfrc_readnl, &! read cldfrc_nl namelist
21 : cldfrc_register, &! add fields to pbuf
22 : cldfrc_init, &! Inititialization of cloud_fraction run-time parameters
23 : cldfrc_getparams, &! public access of tuning parameters
24 : cldfrc, &! Computation of cloud fraction
25 : cldfrc_fice, &! Calculate fraction of condensate in ice phase (radiation partitioning)
26 : dp1, &! parameter for deep convection cloud fraction needed in clubb_intr
27 : dp2 ! parameter for deep convection cloud fraction needed in clubb_intr
28 :
29 : ! Private data
30 : real(r8), parameter :: unset_r8 = huge(1.0_r8)
31 :
32 : ! Top level
33 : integer :: top_lev = 1
34 :
35 : ! Physics buffer indices
36 : integer :: sh_frac_idx = 0
37 : integer :: dp_frac_idx = 0
38 :
39 : ! Namelist variables
40 : logical :: cldfrc_freeze_dry ! switch for Vavrus correction
41 : logical :: cldfrc_ice ! switch to compute ice cloud fraction
42 : real(r8) :: cldfrc_rhminl = unset_r8 ! minimum rh for low stable clouds
43 : real(r8) :: cldfrc_rhminl_adj_land = unset_r8 ! rhminl adjustment for snowfree land
44 : real(r8) :: cldfrc_rhminh = unset_r8 ! minimum rh for high stable clouds
45 : real(r8) :: cldfrc_sh1 = unset_r8 ! parameter for shallow convection cloud fraction
46 : real(r8) :: cldfrc_sh2 = unset_r8 ! parameter for shallow convection cloud fraction
47 : real(r8) :: cldfrc_dp1 = unset_r8 ! parameter for deep convection cloud fraction
48 : real(r8) :: cldfrc_dp2 = unset_r8 ! parameter for deep convection cloud fraction
49 : real(r8) :: cldfrc_premit = unset_r8 ! top pressure bound for mid level cloud
50 : real(r8) :: cldfrc_premib = unset_r8 ! bottom pressure bound for mid level cloud
51 : integer :: cldfrc_iceopt ! option for ice cloud closure
52 : ! 1=wang & sassen 2=schiller (iciwc)
53 : ! 3=wood & field, 4=Wilson (based on smith)
54 : real(r8) :: cldfrc_icecrit = unset_r8 ! Critical RH for ice clouds in Wilson & Ballard closure (smaller = more ice clouds)
55 :
56 : real(r8) :: rhminl ! set from namelist input cldfrc_rhminl
57 : real(r8) :: rhminl_adj_land ! set from namelist input cldfrc_rhminl_adj_land
58 : real(r8) :: rhminh ! set from namelist input cldfrc_rhminh
59 : real(r8) :: sh1, sh2 ! set from namelist input cldfrc_sh1, cldfrc_sh2
60 : real(r8) :: dp1,dp2 ! set from namelist input cldfrc_dp1, cldfrc_dp2
61 : real(r8) :: premit ! set from namelist input cldfrc_premit
62 : real(r8) :: premib ! set from namelist input cldfrc_premib
63 : integer :: iceopt ! set from namelist input cldfrc_iceopt
64 : real(r8) :: icecrit ! set from namelist input cldfrc_icecrit
65 :
66 : ! constants
67 : real(r8), parameter :: pnot = 1.e5_r8 ! reference pressure
68 : real(r8), parameter :: lapse = 6.5e-3_r8 ! U.S. Standard Atmosphere lapse rate
69 : real(r8), parameter :: pretop = 1.0e2_r8 ! pressure bounding high cloud
70 :
71 : integer count
72 :
73 : logical :: inversion_cld_off ! Turns off stratification-based cld frc
74 :
75 : integer :: k700 ! model level nearest 700 mb
76 :
77 : !================================================================================================
78 : contains
79 : !================================================================================================
80 :
81 1536 : subroutine cldfrc_readnl(nlfile)
82 :
83 : use namelist_utils, only: find_group_name
84 : use units, only: getunit, freeunit
85 : use mpishorthand
86 :
87 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
88 :
89 : ! Local variables
90 : integer :: unitn, ierr
91 : character(len=*), parameter :: subname = 'cldfrc_readnl'
92 :
93 : namelist /cldfrc_nl/ cldfrc_freeze_dry, cldfrc_ice, cldfrc_rhminl, &
94 : cldfrc_rhminl_adj_land, cldfrc_rhminh, cldfrc_sh1, &
95 : cldfrc_sh2, cldfrc_dp1, cldfrc_dp2, &
96 : cldfrc_premit, cldfrc_premib, cldfrc_iceopt, &
97 : cldfrc_icecrit
98 : !-----------------------------------------------------------------------------
99 :
100 1536 : if (masterproc) then
101 2 : unitn = getunit()
102 2 : open( unitn, file=trim(nlfile), status='old' )
103 2 : call find_group_name(unitn, 'cldfrc_nl', status=ierr)
104 2 : if (ierr == 0) then
105 2 : read(unitn, cldfrc_nl, iostat=ierr)
106 2 : if (ierr /= 0) then
107 0 : call endrun(subname // ':: ERROR reading namelist')
108 : end if
109 : end if
110 2 : close(unitn)
111 2 : call freeunit(unitn)
112 :
113 : ! set local variables
114 2 : rhminl = cldfrc_rhminl
115 2 : rhminl_adj_land = cldfrc_rhminl_adj_land
116 2 : rhminh = cldfrc_rhminh
117 2 : sh1 = cldfrc_sh1
118 2 : sh2 = cldfrc_sh2
119 2 : dp1 = cldfrc_dp1
120 2 : dp2 = cldfrc_dp2
121 2 : premit = cldfrc_premit
122 2 : premib = cldfrc_premib
123 2 : iceopt = cldfrc_iceopt
124 2 : icecrit = cldfrc_icecrit
125 :
126 : end if
127 :
128 : #ifdef SPMD
129 : ! Broadcast namelist variables
130 1536 : call mpibcast(cldfrc_freeze_dry, 1, mpilog, 0, mpicom)
131 1536 : call mpibcast(cldfrc_ice, 1, mpilog, 0, mpicom)
132 1536 : call mpibcast(rhminl, 1, mpir8, 0, mpicom)
133 1536 : call mpibcast(rhminl_adj_land, 1, mpir8, 0, mpicom)
134 1536 : call mpibcast(rhminh, 1, mpir8, 0, mpicom)
135 1536 : call mpibcast(sh1 , 1, mpir8, 0, mpicom)
136 1536 : call mpibcast(sh2 , 1, mpir8, 0, mpicom)
137 1536 : call mpibcast(dp1 , 1, mpir8, 0, mpicom)
138 1536 : call mpibcast(dp2 , 1, mpir8, 0, mpicom)
139 1536 : call mpibcast(premit, 1, mpir8, 0, mpicom)
140 1536 : call mpibcast(premib, 1, mpir8, 0, mpicom)
141 1536 : call mpibcast(iceopt, 1, mpiint, 0, mpicom)
142 1536 : call mpibcast(icecrit, 1, mpir8, 0, mpicom)
143 : #endif
144 :
145 1536 : end subroutine cldfrc_readnl
146 :
147 : !================================================================================================
148 :
149 1536 : subroutine cldfrc_register
150 :
151 : ! Register fields in the physics buffer.
152 :
153 : use physics_buffer, only : pbuf_add_field, dtype_r8
154 :
155 : !-----------------------------------------------------------------------
156 :
157 1536 : call pbuf_add_field('SH_FRAC', 'physpkg', dtype_r8, (/pcols,pver/), sh_frac_idx)
158 1536 : call pbuf_add_field('DP_FRAC', 'physpkg', dtype_r8, (/pcols,pver/), dp_frac_idx)
159 :
160 1536 : end subroutine cldfrc_register
161 :
162 : !================================================================================================
163 :
164 1536 : subroutine cldfrc_getparams(rhminl_out, rhminl_adj_land_out, rhminh_out, premit_out, &
165 : premib_out, iceopt_out, icecrit_out)
166 : !-----------------------------------------------------------------------
167 : ! Purpose: Return cldfrc tuning parameters
168 : !-----------------------------------------------------------------------
169 :
170 : real(r8), intent(out), optional :: rhminl_out
171 : real(r8), intent(out), optional :: rhminl_adj_land_out
172 : real(r8), intent(out), optional :: rhminh_out
173 : real(r8), intent(out), optional :: premit_out
174 : real(r8), intent(out), optional :: premib_out
175 : integer, intent(out), optional :: iceopt_out
176 : real(r8), intent(out), optional :: icecrit_out
177 :
178 1536 : if ( present(rhminl_out) ) rhminl_out = rhminl
179 1536 : if ( present(rhminl_adj_land_out) ) rhminl_adj_land_out = rhminl_adj_land
180 1536 : if ( present(rhminh_out) ) rhminh_out = rhminh
181 1536 : if ( present(premit_out) ) premit_out = premit
182 1536 : if ( present(premib_out) ) premib_out = premib
183 1536 : if ( present(iceopt_out) ) iceopt_out = iceopt
184 1536 : if ( present(icecrit_out) ) icecrit_out = icecrit
185 :
186 1536 : end subroutine cldfrc_getparams
187 :
188 : !===============================================================================
189 :
190 1536 : subroutine cldfrc_init
191 :
192 : ! Initialize cloud fraction run-time parameters
193 :
194 : use cam_history, only: addfld
195 : use phys_control, only: phys_getopts
196 :
197 : ! query interfaces for scheme settings
198 : character(len=16) :: shallow_scheme, eddy_scheme, macrop_scheme
199 :
200 : integer :: k
201 : !-----------------------------------------------------------------------------
202 :
203 : call phys_getopts(shallow_scheme_out = shallow_scheme ,&
204 : eddy_scheme_out = eddy_scheme ,&
205 1536 : macrop_scheme_out = macrop_scheme )
206 :
207 : ! Limit CAM5 cloud physics to below top cloud level.
208 1536 : if ( .not. (macrop_scheme == "rk" .or. macrop_scheme == "SPCAM_sam1mom")) top_lev = trop_cloud_top_lev
209 :
210 : ! Turn off inversion_cld if any UW PBL scheme is being used
211 1536 : if ( (eddy_scheme .eq. 'diag_TKE' ) .or. (shallow_scheme .eq. 'UW' ) .or.&
212 : (shallow_scheme .eq. 'SPCAM_m2005') ) then
213 0 : inversion_cld_off = .true.
214 : else
215 1536 : inversion_cld_off = .false.
216 : endif
217 :
218 1536 : if ( masterproc ) then
219 2 : write(iulog,*)'tuning parameters cldfrc_init: inversion_cld_off',inversion_cld_off
220 2 : write(iulog,*)'tuning parameters cldfrc_init: dp1',dp1,'dp2',dp2,'sh1',sh1,'sh2',sh2
221 2 : if (shallow_scheme .ne. 'UW' .or. shallow_scheme .eq. 'SPCAM_m2005' ) then
222 2 : write(iulog,*)'tuning parameters cldfrc_init: rhminl',rhminl,'rhminl_adj_land',rhminl_adj_land, &
223 4 : 'rhminh',rhminh,'premit',premit,'premib',premib
224 2 : write(iulog,*)'tuning parameters cldfrc_init: iceopt',iceopt,'icecrit',icecrit
225 : endif
226 : endif
227 :
228 1536 : if (pref_mid(top_lev) > 7.e4_r8) &
229 0 : call endrun ('cldfrc_init: model levels bracketing 700 mb not found')
230 :
231 : ! Find vertical level nearest 700 mb.
232 130560 : k700 = minloc(abs(pref_mid(top_lev:pver) - 7.e4_r8), 1)
233 :
234 1536 : if (masterproc) then
235 2 : write(iulog,*)'cldfrc_init: model level nearest 700 mb is',k700,'which is',pref_mid(k700),'pascals'
236 : end if
237 :
238 3072 : call addfld ('SH_CLD', (/ 'lev' /), 'A', 'fraction', 'Shallow convective cloud cover' )
239 3072 : call addfld ('DP_CLD', (/ 'lev' /), 'A', 'fraction', 'Deep convective cloud cover' )
240 :
241 1536 : end subroutine cldfrc_init
242 :
243 : !===============================================================================
244 :
245 0 : subroutine cldfrc(lchnk ,ncol , pbuf, &
246 : pmid ,temp ,q ,omga , phis, &
247 : shfrc ,use_shfrc, &
248 : cloud ,rhcloud, clc ,pdel , &
249 : cmfmc ,cmfmc2 ,landfrac,snowh ,concld ,cldst , &
250 : ts ,sst ,ps ,zdu ,ocnfrac ,&
251 : rhu00 ,cldice ,icecldf ,liqcldf ,relhum ,dindex )
252 : !-----------------------------------------------------------------------
253 : !
254 : ! Purpose:
255 : ! Compute cloud fraction
256 : !
257 : !
258 : ! Method:
259 : ! This calculate cloud fraction using a relative humidity threshold
260 : ! The threshold depends upon pressure, and upon the presence or absence
261 : ! of convection as defined by a reasonably large vertical mass flux
262 : ! entering that layer from below.
263 : !
264 : ! Author: Many. Last modified by Jim McCaa
265 : !
266 : !-----------------------------------------------------------------------
267 1536 : use cam_history, only: outfld
268 : use physconst, only: cappa, gravit, rair, tmelt
269 : use wv_saturation, only: qsat, qsat_water, svp_ice_vect
270 : use phys_grid, only: get_rlat_all_p, get_rlon_all_p
271 :
272 :
273 : !RBN - Need this to write shallow,deep fraction to phys buffer.
274 : !PJR - we should probably make seperate modules for determining convective
275 : ! clouds and make this one just responsible for relative humidity clouds
276 :
277 : use physics_buffer, only: physics_buffer_desc, pbuf_get_field
278 :
279 : ! Arguments
280 : integer, intent(in) :: lchnk ! chunk identifier
281 : integer, intent(in) :: ncol ! number of atmospheric columns
282 : integer, intent(in) :: dindex ! 0 or 1 to perturb rh
283 :
284 : type(physics_buffer_desc), pointer :: pbuf(:)
285 : real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures
286 : real(r8), intent(in) :: temp(pcols,pver) ! temperature
287 : real(r8), intent(in) :: q(pcols,pver) ! specific humidity
288 : real(r8), intent(in) :: omga(pcols,pver) ! vertical pressure velocity
289 : real(r8), intent(in) :: cmfmc(pcols,pverp) ! convective mass flux--m sub c
290 : real(r8), intent(in) :: cmfmc2(pcols,pverp) ! shallow convective mass flux--m sub c
291 : real(r8), intent(in) :: snowh(pcols) ! snow depth (liquid water equivalent)
292 : real(r8), intent(in) :: pdel(pcols,pver) ! pressure depth of layer
293 : real(r8), intent(in) :: landfrac(pcols) ! Land fraction
294 : real(r8), intent(in) :: ocnfrac(pcols) ! Ocean fraction
295 : real(r8), intent(in) :: ts(pcols) ! surface temperature
296 : real(r8), intent(in) :: sst(pcols) ! sea surface temperature
297 : real(r8), intent(in) :: ps(pcols) ! surface pressure
298 : real(r8), intent(in) :: zdu(pcols,pver) ! detrainment rate from deep convection
299 : real(r8), intent(in) :: phis(pcols) ! surface geopotential
300 : real(r8), intent(in) :: shfrc(pcols,pver) ! cloud fraction from convect_shallow
301 : real(r8), intent(in) :: cldice(pcols,pver) ! cloud ice mixing ratio
302 : logical, intent(in) :: use_shfrc
303 :
304 : ! Output arguments
305 : real(r8), intent(out) :: cloud(pcols,pver) ! cloud fraction
306 : real(r8), intent(out) :: rhcloud(pcols,pver) ! cloud fraction
307 : real(r8), intent(out) :: clc(pcols) ! column convective cloud amount
308 : real(r8), intent(out) :: cldst(pcols,pver) ! cloud fraction
309 : real(r8), intent(out) :: rhu00(pcols,pver) ! RH threshold for cloud
310 : real(r8), intent(out) :: relhum(pcols,pver) ! RH
311 : real(r8), intent(out) :: icecldf(pcols,pver) ! ice cloud fraction
312 : real(r8), intent(out) :: liqcldf(pcols,pver) ! liquid cloud fraction (combined into cloud)
313 :
314 : !---------------------------Local workspace-----------------------------
315 : !
316 : real(r8) concld(pcols,pver) ! convective cloud cover
317 : real(r8) cld ! intermediate scratch variable (low cld)
318 : real(r8) dthdpmn(pcols) ! most stable lapse rate below 750 mb
319 : real(r8) dthdp ! lapse rate (intermediate variable)
320 : real(r8) es(pcols,pver) ! saturation vapor pressure
321 : real(r8) qs(pcols,pver) ! saturation specific humidity
322 : real(r8) rhwght ! weighting function for rhlim transition
323 : real(r8) rh(pcols,pver) ! relative humidity
324 : real(r8) rhdif ! intermediate scratch variable
325 : real(r8) strat ! intermediate scratch variable
326 : real(r8) theta(pcols,pver) ! potential temperature
327 : real(r8) rhlim ! local rel. humidity threshold estimate
328 : real(r8) coef1 ! coefficient to convert mass flux to mb/d
329 : real(r8) clrsky(pcols) ! temporary used in random overlap calc
330 : real(r8) rpdeli(pcols,pver-1) ! 1./(pmid(k+1)-pmid(k))
331 : real(r8) rhpert !the specified perturbation to rh
332 :
333 0 : real(r8), pointer, dimension(:,:) :: deepcu ! deep convection cloud fraction
334 0 : real(r8), pointer, dimension(:,:) :: shallowcu ! shallow convection cloud fraction
335 :
336 : logical cldbnd(pcols) ! region below high cloud boundary
337 :
338 : integer i, ierror, k ! column, level indices
339 : integer kp1, ifld
340 : integer kdthdp(pcols)
341 : integer numkcld ! number of levels in which to allow clouds
342 :
343 : ! In Cloud Ice Content variables
344 : real(r8) :: a,b,c,as,bs,cs !fit parameters
345 : real(r8) :: Kc !constant for ice cloud calc (wood & field)
346 : real(r8) :: ttmp !limited temperature
347 : real(r8) :: icicval !empirical iwc value
348 : real(r8) :: rho !local air density
349 : real(r8) :: esl(pcols,pver) !liq sat vapor pressure
350 : real(r8) :: esi(pcols,pver) !ice sat vapor pressure
351 : real(r8) :: ncf,phi !Wilson and Ballard parameters
352 :
353 : real(r8) thetas(pcols) ! ocean surface potential temperature
354 : real(r8) :: clat(pcols) ! current latitudes(radians)
355 : real(r8) :: clon(pcols) ! current longitudes(radians)
356 :
357 : ! Statement functions
358 : logical land
359 : land(i) = nint(landfrac(i)) == 1
360 :
361 0 : call get_rlat_all_p(lchnk, ncol, clat)
362 0 : call get_rlon_all_p(lchnk, ncol, clon)
363 :
364 0 : call pbuf_get_field(pbuf, sh_frac_idx, shallowcu )
365 0 : call pbuf_get_field(pbuf, dp_frac_idx, deepcu )
366 :
367 : ! Initialise cloud fraction
368 0 : shallowcu = 0._r8
369 0 : deepcu = 0._r8
370 :
371 : !==================================================================================
372 : ! PHILOSOPHY OF PRESENT IMPLEMENTATION
373 : ! Modification to philosophy for ice supersaturation
374 : ! philosophy below is based on RH water only. This is 'liquid condensation'
375 : ! or liquid cloud (even though it will freeze immediately to ice)
376 : ! The idea is that the RH limits for condensation are strict only for
377 : ! water saturation
378 : !
379 : ! Ice clouds are formed by explicit parameterization of ice nucleation.
380 : ! Closure for ice cloud fraction is done on available cloud ice, such that
381 : ! the in-cloud ice content matches an empirical fit
382 : ! thus, icecldf = min(cldice/icicval,1) where icicval = f(temp,cldice,numice)
383 : ! for a first cut, icicval=f(temp) only.
384 : ! Combined cloud fraction is maximum overlap cloud=max(1,max(icecldf,liqcldf))
385 : ! No dA/dt term for ice?
386 : !
387 : ! There are three co-existing cloud types: convective, inversion related low-level
388 : ! stratocumulus, and layered cloud (based on relative humidity). Layered and
389 : ! stratocumulus clouds do not compete with convective cloud for which one creates
390 : ! the most cloud. They contribute collectively to the total grid-box average cloud
391 : ! amount. This is reflected in the way in which the total cloud amount is evaluated
392 : ! (a sum as opposed to a logical "or" operation)
393 : !
394 : !==================================================================================
395 : ! set defaults for rhu00
396 0 : rhu00(:,:) = 2.0_r8
397 : ! define rh perturbation in order to estimate rhdfda
398 0 : rhpert = 0.01_r8
399 :
400 : !set Wang and Sassen IWC paramters
401 0 : a=26.87_r8
402 0 : b=0.569_r8
403 0 : c=0.002892_r8
404 : !set schiller parameters
405 0 : as=-68.4202_r8
406 0 : bs=0.983917_r8
407 0 : cs=2.81795_r8
408 : !set wood and field paramters...
409 0 : Kc=75._r8
410 :
411 : ! Evaluate potential temperature and relative humidity
412 : ! If not computing ice cloud fraction then hybrid RH, if MG then water RH
413 0 : if ( cldfrc_ice ) then
414 0 : do k = top_lev,pver
415 0 : call qsat_water(temp(1:ncol,k), pmid(1:ncol,k), esl(1:ncol,k), qs(1:ncol,k), ncol)
416 0 : call svp_ice_vect(temp(1:ncol,k), esi(1:ncol,k), ncol)
417 : end do
418 : else
419 0 : do k = top_lev,pver
420 0 : call qsat(temp(1:ncol,k), pmid(1:ncol,k), es(1:ncol,k), qs(1:ncol,k), ncol)
421 : end do
422 : end if
423 :
424 0 : cloud = 0._r8
425 0 : icecldf = 0._r8
426 0 : liqcldf = 0._r8
427 0 : rhcloud = 0._r8
428 0 : cldst = 0._r8
429 0 : concld = 0._r8
430 :
431 0 : do k=top_lev,pver
432 0 : theta(:ncol,k) = temp(:ncol,k)*(pnot/pmid(:ncol,k))**cappa
433 :
434 0 : do i=1,ncol
435 0 : rh(i,k) = q(i,k)/qs(i,k)*(1.0_r8+real(dindex,r8)*rhpert)
436 : ! record relhum, rh itself will later be modified related with concld
437 0 : relhum(i,k) = rh(i,k)
438 : end do
439 : end do
440 :
441 : ! Initialize other temporary variables
442 0 : ierror = 0
443 0 : do i=1,ncol
444 : ! Adjust thetas(i) in the presence of non-zero ocean heights.
445 : ! This reduces the temperature for positive heights according to a standard lapse rate.
446 0 : if(ocnfrac(i).gt.0.01_r8) thetas(i) = &
447 0 : ( sst(i) - lapse * phis(i) / gravit) * (pnot/ps(i))**cappa
448 0 : if(ocnfrac(i).gt.0.01_r8.and.sst(i).lt.260._r8) ierror = i
449 0 : clc(i) = 0.0_r8
450 : end do
451 0 : coef1 = gravit*864.0_r8 ! conversion to millibars/day
452 :
453 0 : if (ierror > 0) then
454 0 : write(iulog,*) 'COLDSST: encountered in cldfrc:', lchnk,ierror,ocnfrac(ierror),sst(ierror)
455 : endif
456 :
457 0 : do k=top_lev,pver-1
458 0 : rpdeli(:ncol,k) = 1._r8/(pmid(:ncol,k+1) - pmid(:ncol,k))
459 : end do
460 :
461 : !
462 : ! Estimate of local convective cloud cover based on convective mass flux
463 : ! Modify local large-scale relative humidity to account for presence of
464 : ! convective cloud when evaluating relative humidity based layered cloud amount
465 : !
466 0 : concld(:ncol,top_lev:pver) = 0.0_r8
467 : !
468 : ! cloud mass flux in SI units of kg/m2/s; should produce typical numbers of 20%
469 : ! shallow and deep convective cloudiness are evaluated separately (since processes
470 : ! are evaluated separately) and summed
471 : !
472 : #ifndef PERGRO
473 0 : do k=top_lev,pver
474 0 : do i=1,ncol
475 0 : if ( .not. use_shfrc ) then
476 0 : shallowcu(i,k) = max(0.0_r8,min(sh1*log(1.0_r8+sh2*cmfmc2(i,k+1)),0.30_r8))
477 : else
478 0 : shallowcu(i,k) = shfrc(i,k)
479 : endif
480 0 : deepcu(i,k) = max(0.0_r8,min(dp1*log(1.0_r8+dp2*(cmfmc(i,k+1)-cmfmc2(i,k+1))),0.60_r8))
481 0 : concld(i,k) = min(shallowcu(i,k) + deepcu(i,k),0.80_r8)
482 0 : rh(i,k) = (rh(i,k) - concld(i,k))/(1.0_r8 - concld(i,k))
483 : end do
484 : end do
485 : #endif
486 : !==================================================================================
487 : !
488 : ! ****** Compute layer cloudiness ******
489 : !
490 : !====================================================================
491 : ! Begin the evaluation of layered cloud amount based on (modified) RH
492 : !====================================================================
493 : !
494 0 : numkcld = pver
495 0 : do k=top_lev+1,numkcld
496 0 : kp1 = min(k + 1,pver)
497 0 : do i=1,ncol
498 :
499 : ! This is now designed to apply FOR LIQUID CLOUDS (condensation > RH water)
500 :
501 0 : cldbnd(i) = pmid(i,k).ge.pretop
502 :
503 0 : if ( pmid(i,k).ge.premib ) then
504 : !==============================================================
505 : ! This is the low cloud (below premib) block
506 : !==============================================================
507 : ! enhance low cloud activation over land with no snow cover
508 0 : if (land(i) .and. (snowh(i) <= 0.000001_r8)) then
509 0 : rhlim = rhminl - rhminl_adj_land
510 : else
511 0 : rhlim = rhminl
512 : endif
513 :
514 0 : rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
515 0 : rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
516 :
517 : ! SJV: decrease cloud amount if very low water vapor content
518 : ! (thus very cold): "freeze dry"
519 0 : if (cldfrc_freeze_dry) then
520 0 : rhcloud(i,k) = rhcloud(i,k)*max(0.15_r8,min(1.0_r8,q(i,k)/0.0030_r8))
521 : endif
522 :
523 0 : else if ( pmid(i,k).lt.premit ) then
524 : !==============================================================
525 : ! This is the high cloud (above premit) block
526 : !==============================================================
527 : !
528 0 : rhlim = rhminh
529 : !
530 0 : rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
531 0 : rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
532 : else
533 : !==============================================================
534 : ! This is the middle cloud block
535 : !==============================================================
536 : !
537 : ! linear rh threshold transition between thresholds for low & high cloud
538 : !
539 0 : rhwght = (premib-(max(pmid(i,k),premit)))/(premib-premit)
540 :
541 0 : if (land(i) .and. (snowh(i) <= 0.000001_r8)) then
542 0 : rhlim = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
543 : else
544 0 : rhlim = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
545 : endif
546 0 : rhdif = (rh(i,k) - rhlim)/(1.0_r8-rhlim)
547 0 : rhcloud(i,k) = min(0.999_r8,(max(rhdif,0.0_r8))**2)
548 : end if
549 : !==================================================================================
550 : ! WE NEED TO DOCUMENT THE PURPOSE OF THIS TYPE OF CODE (ASSOCIATED WITH 2ND CALL)
551 : !==================================================================================
552 : ! !
553 : ! ! save rhlim to rhu00, it handles well by itself for low/high cloud
554 : ! !
555 0 : rhu00(i,k)=rhlim
556 : !==================================================================================
557 :
558 0 : if (cldfrc_ice) then
559 :
560 : ! Evaluate ice cloud fraction based on in-cloud ice content
561 :
562 : !--------ICE CLOUD OPTION 1--------Wang & Sassen 2002
563 : ! Evaluate desired in-cloud water content
564 : ! icicval = f(temp,cldice,numice)
565 : ! Start with a function of temperature.
566 : ! Wang & Sassen 2002 (JAS), based on ARM site MMCR (midlat cirrus)
567 : ! parameterization valid for 203-253K
568 : ! icival > 0 for t>195K
569 0 : if (iceopt.lt.3) then
570 0 : if (iceopt.eq.1) then
571 0 : ttmp=max(195._r8,min(temp(i,k),253._r8)) - 273.16_r8
572 0 : icicval=a + b * ttmp + c * ttmp**2._r8
573 : !convert units
574 0 : rho=pmid(i,k)/(rair*temp(i,k))
575 0 : icicval= icicval * 1.e-6_r8 / rho
576 : else
577 : !--------ICE CLOUD OPTION 2--------Schiller 2008 (JGR)
578 : ! Use a curve based on FISH measurements in
579 : ! tropics, mid-lats and arctic. Curve is for 180-250K (raise to 273K?)
580 : ! use median all flights
581 :
582 0 : ttmp=max(190._r8,min(temp(i,k),273.16_r8))
583 0 : icicval = 10._r8 **(as * bs**ttmp + cs)
584 : !convert units from ppmv to kg/kg
585 0 : icicval= icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
586 : endif
587 : !set icecldfraction for OPTION 1 or OPTION2
588 0 : icecldf(i,k) = max(0._r8,min(cldice(i,k)/icicval,1._r8))
589 :
590 0 : else if (iceopt.eq.3) then
591 :
592 : !--------ICE CLOUD OPTION 3--------Wood & Field 2000 (JAS)
593 : ! eq 6: cloud fraction = 1 - exp (-K * qc/qsati)
594 :
595 0 : icecldf(i,k)=1._r8 - exp(-Kc*cldice(i,k)/(qs(i,k)*(esi(i,k)/esl(i,k))))
596 0 : icecldf(i,k)=max(0._r8,min(icecldf(i,k),1._r8))
597 : else
598 : !--------ICE CLOUD OPTION 4--------Wilson and ballard 1999
599 : ! inversion of smith....
600 : ! ncf = cldice / ((1-RHcrit)*qs)
601 : ! then a function of ncf....
602 0 : ncf =cldice(i,k)/((1._r8 - icecrit)*qs(i,k))
603 0 : if (ncf.le.0._r8) then
604 0 : icecldf(i,k)=0._r8
605 0 : else if (ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8) then
606 0 : icecldf(i,k)=0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
607 0 : else if (ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8) then
608 0 : phi=(acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
609 0 : icecldf(i,k)=(1._r8 - 4._r8 * cos(phi) * cos(phi))
610 : else
611 0 : icecldf(i,k)=1._r8
612 : endif
613 0 : icecldf(i,k)=max(0._r8,min(icecldf(i,k),1._r8))
614 : endif
615 : !TEST: if ice present, icecldf=1.
616 : ! if (cldice(i,k).ge.1.e-8_r8) then
617 : ! icecldf(i,k) = 0.99_r8
618 : ! endif
619 :
620 : !! if ((cldice(i,k) .gt. icicval) .or. ((cldice(i,k) .gt. 0._r8) .and. (icecldf(i,k) .eq. 0._r8))) then
621 : ! if (cldice(i,k) .gt. 1.e-8_r8) then
622 : ! write(iulog,*) 'i,k,pmid,rho,t,cldice,icicval,icecldf,rhcloud: ', &
623 : ! i,k,pmid(i,k),rho,temp(i,k),cldice(i,k),icicval,icecldf(i,k),rhcloud(i,k)
624 : ! endif
625 :
626 : ! Combine ice and liquid cloud fraction assuming maximum overlap.
627 : ! Combined cloud fraction is maximum overlap
628 : ! cloud(i,k)=min(1._r8,max(icecldf(i,k),rhcloud(i,k)))
629 :
630 0 : liqcldf(i,k)=(1._r8 - icecldf(i,k))* rhcloud(i,k)
631 0 : cloud(i,k)=liqcldf(i,k) + icecldf(i,k)
632 : else
633 : ! For RK microphysics
634 0 : cloud(i,k) = rhcloud(i,k)
635 : end if
636 : end do
637 : end do
638 : !
639 : ! Add in the marine strat
640 : ! MARINE STRATUS SHOULD BE A SPECIAL CASE OF LAYERED CLOUD
641 : ! CLOUD CURRENTLY CONTAINS LAYERED CLOUD DETERMINED BY RH CRITERIA
642 : ! TAKE THE MAXIMUM OF THE DIAGNOSED LAYERED CLOUD OR STRATOCUMULUS
643 : !
644 : !===================================================================================
645 : !
646 : ! SOME OBSERVATIONS ABOUT THE FOLLOWING SECTION OF CODE (missed in earlier look)
647 : ! K700 IS SET AS A CONSTANT BASED ON HYBRID COORDINATE: IT DOES NOT DEPEND ON
648 : ! LOCAL PRESSURE; THERE IS NO PRESSURE RAMP => LOOKS LEVEL DEPENDENT AND
649 : ! DISCONTINUOUS IN SPACE (I.E., STRATUS WILL END SUDDENLY WITH NO TRANSITION)
650 : !
651 : ! IT APPEARS THAT STRAT IS EVALUATED ACCORDING TO KLEIN AND HARTMANN; HOWEVER,
652 : ! THE ACTUAL STRATUS AMOUNT (CLDST) APPEARS TO DEPEND DIRECTLY ON THE RH BELOW
653 : ! THE STRONGEST PART OF THE LOW LEVEL INVERSION.
654 : !PJR answers: 1) the rh limitation is a physical/mathematical limitation
655 : ! cant have more cloud than there is RH
656 : ! allowed the cloud to exist two layers below the inversion
657 : ! because the numerics frequently make 50% relative humidity
658 : ! in level below the inversion which would allow no cloud
659 : ! 2) since the cloud is only allowed over ocean, it should
660 : ! be very insensitive to surface pressure (except due to
661 : ! spectral ringing, which also causes so many other problems
662 : ! I didnt worry about it.
663 : !
664 : !==================================================================================
665 0 : if (.not.inversion_cld_off) then
666 : !
667 : ! Find most stable level below 750 mb for evaluating stratus regimes
668 : !
669 0 : do i=1,ncol
670 : ! Nothing triggers unless a stability greater than this minimum threshold is found
671 0 : dthdpmn(i) = -0.125_r8
672 0 : kdthdp(i) = 0
673 : end do
674 : !
675 0 : do k=top_lev+1,pver
676 0 : do i=1,ncol
677 0 : if (pmid(i,k) >= premib .and. ocnfrac(i).gt. 0.01_r8) then
678 : ! I think this is done so that dtheta/dp is in units of dg/mb (JJH)
679 0 : dthdp = 100.0_r8*(theta(i,k) - theta(i,k-1))*rpdeli(i,k-1)
680 0 : if (dthdp < dthdpmn(i)) then
681 0 : dthdpmn(i) = dthdp
682 0 : kdthdp(i) = k ! index of interface of max inversion
683 : end if
684 : end if
685 : end do
686 : end do
687 :
688 : ! Also check between the bottom layer and the surface
689 : ! Only perform this check if the criteria were not met above
690 :
691 0 : do i = 1,ncol
692 0 : if ( kdthdp(i) .eq. 0 .and. ocnfrac(i).gt.0.01_r8) then
693 0 : dthdp = 100.0_r8 * (thetas(i) - theta(i,pver)) / (ps(i)-pmid(i,pver))
694 0 : if (dthdp < dthdpmn(i)) then
695 0 : dthdpmn(i) = dthdp
696 0 : kdthdp(i) = pver ! index of interface of max inversion
697 : endif
698 : endif
699 : enddo
700 :
701 0 : do i=1,ncol
702 0 : if (kdthdp(i) /= 0) then
703 0 : k = kdthdp(i)
704 0 : kp1 = min(k+1,pver)
705 : ! Note: strat will be zero unless ocnfrac > 0.01
706 0 : strat = min(1._r8,max(0._r8, ocnfrac(i) * ((theta(i,k700)-thetas(i))*.057_r8-.5573_r8) ) )
707 : !
708 : ! assign the stratus to the layer just below max inversion
709 : ! the relative humidity changes so rapidly across the inversion
710 : ! that it is not safe to just look immediately below the inversion
711 : ! so limit the stratus cloud by rh in both layers below the inversion
712 : !
713 0 : cldst(i,k) = min(strat,max(rh(i,k),rh(i,kp1)))
714 : end if
715 : end do
716 : end if ! .not.inversion_cld_off
717 :
718 0 : do k=top_lev,pver
719 0 : do i=1,ncol
720 : !
721 : ! which is greater; standard layered cloud amount or stratocumulus diagnosis
722 : !
723 0 : cloud(i,k) = max(rhcloud(i,k),cldst(i,k))
724 : !
725 : ! add in the contributions of convective cloud (determined separately and accounted
726 : ! for by modifications to the large-scale relative humidity.
727 : !
728 0 : cloud(i,k) = min(cloud(i,k)+concld(i,k), 1.0_r8)
729 : end do
730 : end do
731 :
732 0 : call outfld( 'SH_CLD ', shallowcu , pcols, lchnk )
733 0 : call outfld( 'DP_CLD ', deepcu , pcols, lchnk )
734 :
735 : !
736 0 : return
737 0 : end subroutine cldfrc
738 :
739 : !================================================================================================
740 :
741 1495368 : subroutine cldfrc_fice(ncol, t, fice, fsnow)
742 : !
743 : ! Compute the fraction of the total cloud water which is in ice phase.
744 : ! The fraction depends on temperature only.
745 : ! This is the form that was used for radiation, the code came from cldefr originally
746 : !
747 : ! Author: B. A. Boville Sept 10, 2002
748 : ! modified: PJR 3/13/03 (added fsnow to ascribe snow production for convection )
749 : !-----------------------------------------------------------------------
750 0 : use physconst, only: tmelt
751 :
752 : ! Arguments
753 : integer, intent(in) :: ncol ! number of active columns
754 : real(r8), intent(in) :: t(:,:) ! temperature
755 :
756 : real(r8), intent(out) :: fice(:,:) ! Fractional ice content within cloud
757 : real(r8), intent(out) :: fsnow(:,:) ! Fractional snow content for convection
758 :
759 : ! Local variables
760 : real(r8) :: tmax_fice ! max temperature for cloud ice formation
761 : real(r8) :: tmin_fice ! min temperature for cloud ice formation
762 : real(r8) :: tmax_fsnow ! max temperature for transition to convective snow
763 : real(r8) :: tmin_fsnow ! min temperature for transition to convective snow
764 :
765 : integer :: i,k ! loop indexes
766 :
767 : !-----------------------------------------------------------------------
768 :
769 1495368 : tmax_fice = tmelt - 10._r8 ! max temperature for cloud ice formation
770 1495368 : tmin_fice = tmax_fice - 30._r8 ! min temperature for cloud ice formation
771 1495368 : tmax_fsnow = tmelt ! max temperature for transition to convective snow
772 1495368 : tmin_fsnow = tmelt - 5._r8 ! min temperature for transition to convective snow
773 :
774 227713248 : fice(:,:top_lev-1) = 0._r8
775 227713248 : fsnow(:,:top_lev-1) = 0._r8
776 :
777 : ! Define fractional amount of cloud that is ice
778 127106280 : do k=top_lev,pver
779 2098905480 : do i=1,ncol
780 :
781 : ! If warmer than tmax then water phase
782 1971799200 : if (t(i,k) > tmax_fice) then
783 502589268 : fice(i,k) = 0.0_r8
784 :
785 : ! If colder than tmin then ice phase
786 1469209932 : else if (t(i,k) < tmin_fice) then
787 1082589970 : fice(i,k) = 1.0_r8
788 :
789 : ! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax
790 : else
791 386619962 : fice(i,k) =(tmax_fice - t(i,k)) / (tmax_fice - tmin_fice)
792 : end if
793 :
794 : ! snow fraction partitioning
795 :
796 : ! If warmer than tmax then water phase
797 2097410112 : if (t(i,k) > tmax_fsnow) then
798 375834569 : fsnow(i,k) = 0.0_r8
799 :
800 : ! If colder than tmin then ice phase
801 1595964631 : else if (t(i,k) < tmin_fsnow) then
802 1531620515 : fsnow(i,k) = 1.0_r8
803 :
804 : ! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax
805 : else
806 64344116 : fsnow(i,k) =(tmax_fsnow - t(i,k)) / (tmax_fsnow - tmin_fsnow)
807 : end if
808 :
809 : end do
810 : end do
811 :
812 1495368 : end subroutine cldfrc_fice
813 :
814 : end module cloud_fraction
|
