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