Line data Source code
1 : module cldfrc2m
2 :
3 : ! cloud fraction calculations
4 :
5 : use shr_kind_mod, only: r8=>shr_kind_r8
6 : use spmd_utils, only: masterproc
7 : use ppgrid, only: pcols
8 : use physconst, only: rair
9 : use wv_saturation, only: qsat_water, svp_water, svp_ice, &
10 : svp_water_vect, svp_ice_vect
11 : use cam_logfile, only: iulog
12 : use cam_abortutils, only: endrun
13 :
14 : implicit none
15 : private
16 : save
17 :
18 : public :: &
19 : cldfrc2m_readnl, &
20 : cldfrc2m_init, &
21 : astG_PDF_single, &
22 : astG_PDF, &
23 : astG_RHU_single, &
24 : astG_RHU, &
25 : aist_single, &
26 : aist_vector, &
27 : CAMstfrac, &
28 : rhmini_const, &
29 : rhmaxi_const, &
30 : rhminis_const, &
31 : rhmaxis_const
32 :
33 : ! Namelist variables
34 : real(r8) :: cldfrc2m_rhmini ! Minimum rh for ice cloud fraction > 0.
35 : real(r8) :: cldfrc2m_rhmaxi
36 : real(r8) :: cldfrc2m_rhminis ! Minimum rh for ice cloud fraction > 0 in the stratsophere.
37 : real(r8) :: cldfrc2m_rhmaxis
38 : real(r8) :: cldfrc2m_qist_min ! Minimum in-stratus IWC constraint [ kg/kg ]
39 : real(r8) :: cldfrc2m_qist_max ! Maximum in-stratus IWC constraint [ kg/kg ]
40 : logical :: cldfrc2m_do_subgrid_growth = .false.
41 : logical :: cldfrc2m_do_avg_aist_algs = .false.
42 : ! -------------------------- !
43 : ! Parameters for Ice Stratus !
44 : ! -------------------------- !
45 : real(r8), protected :: rhmini_const ! Minimum rh for ice cloud fraction > 0.
46 : real(r8), protected :: rhmaxi_const
47 : real(r8), protected :: rhminis_const ! Minimum rh for ice cloud fraction > 0.
48 : real(r8), protected :: rhmaxis_const
49 :
50 : ! ----------------------------- !
51 : ! Parameters for Liquid Stratus !
52 : ! ----------------------------- !
53 :
54 : logical, parameter :: CAMstfrac = .false. ! If .true. (.false.),
55 : ! use Slingo (triangular PDF-based) liquid stratus fraction
56 : logical, parameter :: freeze_dry = .false. ! If .true., use 'freeze dry' in liquid stratus fraction formula
57 : real(r8) :: rhminl_const ! Critical RH for low-level liquid stratus clouds
58 : real(r8) :: rhminl_adj_land_const ! rhminl adjustment for snowfree land
59 : real(r8) :: rhminh_const ! Critical RH for high-level liquid stratus clouds
60 : real(r8) :: premit ! Top height for mid-level liquid stratus fraction
61 : real(r8) :: premib ! Bottom height for mid-level liquid stratus fraction
62 : integer :: iceopt ! option for ice cloud closure
63 : ! 1=wang & sassen 2=schiller (iciwc)
64 : ! 3=wood & field, 4=Wilson (based on smith)
65 : ! 5=modified slingo (ssat & empyt cloud)
66 : real(r8) :: icecrit ! Critical RH for ice clouds in Wilson & Ballard closure
67 : ! ( smaller = more ice clouds )
68 :
69 : !================================================================================================
70 : contains
71 : !================================================================================================
72 :
73 1536 : subroutine cldfrc2m_readnl(nlfile)
74 :
75 : use namelist_utils, only: find_group_name
76 : use units, only: getunit, freeunit
77 : use spmd_utils, only: mpicom, masterprocid, mpi_logical, mpi_real8
78 :
79 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
80 :
81 : ! Local variables
82 : integer :: unitn, ierr
83 : character(len=*), parameter :: subname = 'cldfrc2m_readnl'
84 :
85 : namelist /cldfrc2m_nl/ cldfrc2m_rhmini, cldfrc2m_rhmaxi, cldfrc2m_rhminis, cldfrc2m_rhmaxis, cldfrc2m_do_subgrid_growth, &
86 : cldfrc2m_qist_min, cldfrc2m_qist_max, cldfrc2m_do_avg_aist_algs
87 : !-----------------------------------------------------------------------------
88 :
89 1536 : if (masterproc) then
90 2 : unitn = getunit()
91 2 : open( unitn, file=trim(nlfile), status='old' )
92 2 : call find_group_name(unitn, 'cldfrc2m_nl', status=ierr)
93 2 : if (ierr == 0) then
94 2 : read(unitn, cldfrc2m_nl, iostat=ierr)
95 2 : if (ierr /= 0) then
96 0 : call endrun(subname // ':: ERROR reading namelist')
97 : end if
98 : end if
99 2 : close(unitn)
100 2 : call freeunit(unitn)
101 :
102 : ! set local variables
103 2 : rhmini_const = cldfrc2m_rhmini
104 2 : rhmaxi_const = cldfrc2m_rhmaxi
105 2 : rhminis_const = cldfrc2m_rhminis
106 2 : rhmaxis_const = cldfrc2m_rhmaxis
107 : end if
108 :
109 : ! Broadcast namelist variables
110 1536 : call mpi_bcast(rhmini_const, 1, mpi_real8, masterprocid, mpicom, ierr)
111 1536 : call mpi_bcast(rhmaxi_const, 1, mpi_real8, masterprocid, mpicom, ierr)
112 1536 : call mpi_bcast(rhminis_const, 1, mpi_real8, masterprocid, mpicom, ierr)
113 1536 : call mpi_bcast(rhmaxis_const, 1, mpi_real8, masterprocid, mpicom, ierr)
114 1536 : call mpi_bcast(cldfrc2m_qist_min, 1, mpi_real8, masterprocid, mpicom, ierr)
115 1536 : call mpi_bcast(cldfrc2m_qist_max, 1, mpi_real8, masterprocid, mpicom, ierr)
116 1536 : call mpi_bcast(cldfrc2m_do_subgrid_growth, 1, mpi_logical,masterprocid, mpicom, ierr)
117 1536 : call mpi_bcast(cldfrc2m_do_avg_aist_algs, 1, mpi_logical,masterprocid, mpicom, ierr)
118 :
119 1536 : end subroutine cldfrc2m_readnl
120 :
121 : !================================================================================================
122 :
123 1536 : subroutine cldfrc2m_init()
124 :
125 : use cloud_fraction, only: cldfrc_getparams
126 :
127 : call cldfrc_getparams(rhminl_out=rhminl_const, rhminl_adj_land_out=rhminl_adj_land_const, &
128 : rhminh_out=rhminh_const, premit_out=premit, premib_out=premib, &
129 1536 : iceopt_out=iceopt, icecrit_out=icecrit)
130 :
131 1536 : if( masterproc ) then
132 2 : write(iulog,*) 'cldfrc2m parameters:'
133 2 : write(iulog,*) ' rhminl = ', rhminl_const
134 2 : write(iulog,*) ' rhminl_adj_land = ', rhminl_adj_land_const
135 2 : write(iulog,*) ' rhminh = ', rhminh_const
136 2 : write(iulog,*) ' premit = ', premit
137 2 : write(iulog,*) ' premib = ', premib
138 2 : write(iulog,*) ' iceopt = ', iceopt
139 2 : write(iulog,*) ' icecrit = ', icecrit
140 2 : write(iulog,*) ' rhmini = ', rhmini_const
141 2 : write(iulog,*) ' rhmaxi = ', rhmaxi_const
142 2 : write(iulog,*) ' rhminis = ', rhminis_const
143 2 : write(iulog,*) ' rhmaxis = ', rhmaxis_const
144 2 : write(iulog,*) ' do_subgrid_growth = ', cldfrc2m_do_subgrid_growth
145 2 : write(iulog,*) ' do_avg_aist_algs = ', cldfrc2m_do_avg_aist_algs
146 2 : write(iulog,*) ' cldfrc2m_qist_min = ', cldfrc2m_qist_min
147 2 : write(iulog,*) ' cldfrc2m_qist_max = ', cldfrc2m_qist_max
148 : end if
149 :
150 1536 : end subroutine cldfrc2m_init
151 :
152 : !================================================================================================
153 :
154 :
155 0 : subroutine astG_PDF_single(U, p, qv, landfrac, snowh, a, Ga, orhmin, &
156 : rhminl_in, rhminl_adj_land_in, rhminh_in )
157 :
158 : ! --------------------------------------------------------- !
159 : ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the !
160 : ! analytical formulation of triangular PDF. !
161 : ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', !
162 : ! so using constant 'dV' assume that width is proportional !
163 : ! to the saturation specific humidity. !
164 : ! dV ~ 0.1. !
165 : ! cldrh : RH of in-stratus( = 1 if no supersaturation) !
166 : ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
167 : ! G is discontinuous across U = 1. In fact, it does not !
168 : ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that !
169 : ! they will produce the same results. !
170 : ! --------------------------------------------------------- !
171 :
172 : real(r8), intent(in) :: U ! Relative humidity
173 : real(r8), intent(in) :: p ! Pressure [Pa]
174 : real(r8), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg]
175 : real(r8), intent(in) :: landfrac ! Land fraction
176 : real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
177 :
178 : real(r8), intent(out) :: a ! Stratus fraction
179 : real(r8), intent(out) :: Ga ! dU/da
180 : real(r8), optional, intent(out) :: orhmin ! Critical RH
181 :
182 : real(r8), optional, intent(in) :: rhminl_in ! Critical relative humidity for low-level liquid stratus
183 : real(r8), optional, intent(in) :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land
184 : real(r8), optional, intent(in) :: rhminh_in ! Critical relative humidity for high-level liquid stratus
185 :
186 : ! Local variables
187 : integer :: i ! Loop indexes
188 : real(r8) dV ! Width of triangular PDF
189 : real(r8) cldrh ! RH of stratus cloud
190 : real(r8) rhmin ! Critical RH
191 : real(r8) rhwght
192 :
193 : real(r8) :: rhminl
194 : real(r8) :: rhminl_adj_land
195 : real(r8) :: rhminh
196 :
197 : ! Statement functions
198 : logical land
199 0 : land = nint(landfrac) == 1
200 :
201 : ! ---------- !
202 : ! Parameters !
203 : ! ---------- !
204 :
205 0 : cldrh = 1.0_r8
206 :
207 0 : rhminl = rhminl_const
208 0 : if (present(rhminl_in)) rhminl = rhminl_in
209 0 : rhminl_adj_land = rhminl_adj_land_const
210 0 : if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in
211 0 : rhminh = rhminh_const
212 0 : if (present(rhminh_in)) rhminh = rhminh_in
213 :
214 : ! ---------------- !
215 : ! Main computation !
216 : ! ---------------- !
217 :
218 0 : if( p .ge. premib ) then
219 :
220 0 : if( land .and. (snowh.le.0.000001_r8) ) then
221 0 : rhmin = rhminl - rhminl_adj_land
222 : else
223 : rhmin = rhminl
224 : endif
225 :
226 0 : dV = cldrh - rhmin
227 :
228 0 : if( U .ge. 1._r8 ) then
229 0 : a = 1._r8
230 0 : Ga = 1.e10_r8
231 0 : elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
232 0 : a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
233 0 : Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
234 0 : elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
235 : a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
236 0 : (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
237 0 : Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
238 0 : elseif( U .le. (cldrh-dV) ) then
239 0 : a = 0._r8
240 0 : Ga = 1.e10_r8
241 : endif
242 :
243 : if( freeze_dry ) then
244 : a = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
245 : Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
246 : endif
247 :
248 0 : elseif( p .lt. premit ) then
249 :
250 0 : rhmin = rhminh
251 0 : dV = cldrh - rhmin
252 :
253 0 : if( U .ge. 1._r8 ) then
254 0 : a = 1._r8
255 0 : Ga = 1.e10_r8
256 0 : elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
257 0 : a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
258 0 : Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
259 0 : elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
260 : a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
261 0 : (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
262 0 : Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
263 0 : elseif( U .le. (cldrh-dV) ) then
264 0 : a = 0._r8
265 0 : Ga = 1.e10_r8
266 : endif
267 :
268 : else
269 :
270 0 : rhwght = (premib-(max(p,premit)))/(premib-premit)
271 :
272 : ! if( land .and. (snowh.le.0.000001_r8) ) then
273 : ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
274 : ! else
275 0 : rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
276 : ! endif
277 :
278 0 : dV = cldrh - rhmin
279 :
280 0 : if( U .ge. 1._r8 ) then
281 0 : a = 1._r8
282 0 : Ga = 1.e10_r8
283 0 : elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
284 0 : a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
285 0 : Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
286 0 : elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
287 : a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
288 0 : (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
289 0 : Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
290 0 : elseif( U .le. (cldrh-dV) ) then
291 0 : a = 0._r8
292 0 : Ga = 1.e10_r8
293 : endif
294 :
295 : endif
296 :
297 0 : if (present(orhmin)) orhmin = rhmin
298 :
299 1536 : end subroutine astG_PDF_single
300 :
301 : !================================================================================================
302 :
303 0 : subroutine astG_PDF(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, &
304 : rhminl_in, rhminl_adj_land_in, rhminh_in )
305 :
306 : ! --------------------------------------------------------- !
307 : ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the !
308 : ! analytical formulation of triangular PDF. !
309 : ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', !
310 : ! so using constant 'dV' assume that width is proportional !
311 : ! to the saturation specific humidity. !
312 : ! dV ~ 0.1. !
313 : ! cldrh : RH of in-stratus( = 1 if no supersaturation) !
314 : ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
315 : ! G is discontinuous across U = 1. In fact, it does not !
316 : ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that !
317 : ! they will produce the same results. !
318 : ! --------------------------------------------------------- !
319 :
320 : real(r8), intent(in) :: U_in(pcols) ! Relative humidity
321 : real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa]
322 : real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor specific humidity [kg/kg]
323 : real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction
324 : real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent)
325 :
326 : real(r8), intent(out) :: a_out(pcols) ! Stratus fraction
327 : real(r8), intent(out) :: Ga_out(pcols) ! dU/da
328 : integer, intent(in) :: ncol
329 :
330 : real(r8), optional, intent(in) :: rhminl_in(pcols) ! Critical relative humidity for low-level liquid stratus
331 : real(r8), optional, intent(in) :: rhminl_adj_land_in(pcols) ! Adjustment drop of rhminl over the land
332 : real(r8), optional, intent(in) :: rhminh_in(pcols) ! Critical relative humidity for high-level liquid stratus
333 :
334 : real(r8) :: rhminl ! Critical relative humidity for low-level liquid stratus
335 : real(r8) :: rhminl_adj_land ! Adjustment drop of rhminl over the land
336 : real(r8) :: rhminh ! Critical relative humidity for high-level liquid stratus
337 :
338 : real(r8) :: U ! Relative humidity
339 : real(r8) :: p ! Pressure [Pa]
340 : real(r8) :: qv ! Grid-mean water vapor specific humidity [kg/kg]
341 : real(r8) :: landfrac ! Land fraction
342 : real(r8) :: snowh ! Snow depth (liquid water equivalent)
343 :
344 : real(r8) :: a ! Stratus fraction
345 : real(r8) :: Ga ! dU/da
346 :
347 : ! Local variables
348 : integer :: i ! Loop indexes
349 : real(r8) dV ! Width of triangular PDF
350 : real(r8) cldrh ! RH of stratus cloud
351 : real(r8) rhmin ! Critical RH
352 : real(r8) rhwght
353 :
354 : ! Statement functions
355 : logical land
356 : land(i) = nint(landfrac_in(i)) == 1
357 :
358 : ! ---------- !
359 : ! Parameters !
360 : ! ---------- !
361 :
362 0 : cldrh = 1.0_r8
363 :
364 0 : rhminl = rhminl_const
365 0 : rhminl_adj_land = rhminl_adj_land_const
366 0 : rhminh = rhminh_const
367 :
368 : ! ---------------- !
369 : ! Main computation !
370 : ! ---------------- !
371 :
372 0 : a_out(:) = 0._r8
373 0 : Ga_out(:) = 0._r8
374 :
375 0 : do i = 1, ncol
376 :
377 0 : U = U_in(i)
378 0 : p = p_in(i)
379 0 : qv = qv_in(i)
380 0 : landfrac = landfrac_in(i)
381 0 : snowh = snowh_in(i)
382 :
383 0 : if (present(rhminl_in)) rhminl = rhminl_in(i)
384 0 : if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i)
385 0 : if (present(rhminh_in)) rhminh = rhminh_in(i)
386 :
387 0 : if( p .ge. premib ) then
388 :
389 0 : if( land(i) .and. (snowh.le.0.000001_r8) ) then
390 0 : rhmin = rhminl - rhminl_adj_land
391 : else
392 : rhmin = rhminl
393 : endif
394 :
395 0 : dV = cldrh - rhmin
396 :
397 0 : if( U .ge. 1._r8 ) then
398 : a = 1._r8
399 : Ga = 1.e10_r8
400 0 : elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
401 0 : a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
402 0 : Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
403 0 : elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
404 : a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
405 0 : (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
406 0 : Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
407 0 : elseif( U .le. (cldrh-dV) ) then
408 0 : a = 0._r8
409 0 : Ga = 1.e10_r8
410 : endif
411 :
412 : if( freeze_dry ) then
413 : a = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
414 : Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
415 : endif
416 :
417 0 : elseif( p .lt. premit ) then
418 :
419 0 : rhmin = rhminh
420 0 : dV = cldrh - rhmin
421 :
422 0 : if( U .ge. 1._r8 ) then
423 : a = 1._r8
424 : Ga = 1.e10_r8
425 0 : elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
426 0 : a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
427 0 : Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
428 0 : elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
429 : a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
430 0 : (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
431 0 : Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
432 0 : elseif( U .le. (cldrh-dV) ) then
433 0 : a = 0._r8
434 0 : Ga = 1.e10_r8
435 : endif
436 :
437 : else
438 :
439 0 : rhwght = (premib-(max(p,premit)))/(premib-premit)
440 :
441 : ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
442 : ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
443 : ! else
444 0 : rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
445 : ! endif
446 :
447 0 : dV = cldrh - rhmin
448 :
449 0 : if( U .ge. 1._r8 ) then
450 : a = 1._r8
451 : Ga = 1.e10_r8
452 0 : elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
453 0 : a = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
454 0 : Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
455 0 : elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
456 : a = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
457 0 : (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
458 0 : Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
459 0 : elseif( U .le. (cldrh-dV) ) then
460 0 : a = 0._r8
461 0 : Ga = 1.e10_r8
462 : endif
463 :
464 : endif
465 :
466 0 : a_out(i) = a
467 0 : Ga_out(i) = Ga
468 :
469 : enddo
470 :
471 0 : end subroutine astG_PDF
472 : !================================================================================================
473 :
474 0 : subroutine astG_RHU_single(U, p, qv, landfrac, snowh, a, Ga, orhmin, &
475 : rhminl_in, rhminl_adj_land_in, rhminh_in )
476 :
477 : ! --------------------------------------------------------- !
478 : ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the !
479 : ! CAM35 cloud fraction formula. !
480 : ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core !
481 : ! For the other cases, I should re-define 'rhminl,rhminh' & !
482 : ! 'premib,premit'. !
483 : ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
484 : ! G is discontinuous across U = 1. !
485 : ! --------------------------------------------------------- !
486 :
487 : real(r8), intent(in) :: U ! Relative humidity
488 : real(r8), intent(in) :: p ! Pressure [Pa]
489 : real(r8), intent(in) :: qv ! Grid-mean water vapor specific humidity [kg/kg]
490 : real(r8), intent(in) :: landfrac ! Land fraction
491 : real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
492 :
493 : real(r8), intent(out) :: a ! Stratus fraction
494 : real(r8), intent(out) :: Ga ! dU/da
495 : real(r8), optional, intent(out) :: orhmin ! Critical RH
496 :
497 : real(r8), optional, intent(in) :: rhminl_in ! Critical relative humidity for low-level liquid stratus
498 : real(r8), optional, intent(in) :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land
499 : real(r8), optional, intent(in) :: rhminh_in ! Critical relative humidity for high-level liquid stratus
500 :
501 : ! Local variables
502 : real(r8) rhmin ! Critical RH
503 : real(r8) rhdif ! Factor for stratus fraction
504 : real(r8) rhwght
505 :
506 : real(r8) :: rhminl
507 : real(r8) :: rhminl_adj_land
508 : real(r8) :: rhminh
509 :
510 : ! Statement functions
511 : logical land
512 0 : land = nint(landfrac) == 1
513 :
514 0 : rhminl = rhminl_const
515 0 : if (present(rhminl_in)) rhminl = rhminl_in
516 0 : rhminl_adj_land = rhminl_adj_land_const
517 0 : if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in
518 0 : rhminh = rhminh_const
519 0 : if (present(rhminh_in)) rhminh = rhminh_in
520 :
521 : ! ---------------- !
522 : ! Main computation !
523 : ! ---------------- !
524 :
525 0 : if( p .ge. premib ) then
526 :
527 0 : if( land .and. (snowh.le.0.000001_r8) ) then
528 0 : rhmin = rhminl - rhminl_adj_land
529 : else
530 : rhmin = rhminl
531 : endif
532 0 : rhdif = (U-rhmin)/(1.0_r8-rhmin)
533 0 : a = min(1._r8,(max(rhdif,0.0_r8))**2)
534 0 : if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
535 0 : Ga = 1.e20_r8
536 : else
537 0 : Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
538 : endif
539 : if( freeze_dry ) then
540 : a = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
541 : Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
542 : endif
543 :
544 0 : elseif( p .lt. premit ) then
545 :
546 0 : rhmin = rhminh
547 0 : rhdif = (U-rhmin)/(1.0_r8-rhmin)
548 0 : a = min(1._r8,(max(rhdif,0._r8))**2)
549 0 : if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
550 0 : Ga = 1.e20_r8
551 : else
552 0 : Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
553 : endif
554 :
555 : else
556 :
557 0 : rhwght = (premib-(max(p,premit)))/(premib-premit)
558 :
559 : ! if( land .and. (snowh.le.0.000001_r8) ) then
560 : ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
561 : ! else
562 0 : rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
563 : ! endif
564 :
565 0 : rhdif = (U-rhmin)/(1.0_r8-rhmin)
566 0 : a = min(1._r8,(max(rhdif,0._r8))**2)
567 0 : if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
568 0 : Ga = 1.e10_r8
569 : else
570 0 : Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
571 : endif
572 :
573 : endif
574 :
575 0 : if (present(orhmin)) orhmin = rhmin
576 :
577 0 : end subroutine astG_RHU_single
578 :
579 : !================================================================================================
580 :
581 0 : subroutine astG_RHU(U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol, &
582 : rhminl_in, rhminl_adj_land_in, rhminh_in )
583 :
584 : ! --------------------------------------------------------- !
585 : ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the !
586 : ! CAM35 cloud fraction formula. !
587 : ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core !
588 : ! For the other cases, I should re-define 'rhminl,rhminh' & !
589 : ! 'premib,premit'. !
590 : ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
591 : ! G is discontinuous across U = 1. !
592 : ! --------------------------------------------------------- !
593 :
594 : real(r8), intent(in) :: U_in(pcols) ! Relative humidity
595 : real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa]
596 : real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor specific humidity [kg/kg]
597 : real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction
598 : real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent)
599 :
600 : real(r8), intent(out) :: a_out(pcols) ! Stratus fraction
601 : real(r8), intent(out) :: Ga_out(pcols) ! dU/da
602 : integer, intent(in) :: ncol
603 :
604 : real(r8), optional, intent(in) :: rhminl_in(pcols) ! Critical relative humidity for low-level liquid stratus
605 : real(r8), optional, intent(in) :: rhminl_adj_land_in(pcols) ! Adjustment drop of rhminl over the land
606 : real(r8), optional, intent(in) :: rhminh_in(pcols) ! Critical relative humidity for high-level liquid stratus
607 :
608 : real(r8) :: U ! Relative humidity
609 : real(r8) :: p ! Pressure [Pa]
610 : real(r8) :: qv ! Grid-mean water vapor specific humidity [kg/kg]
611 : real(r8) :: landfrac ! Land fraction
612 : real(r8) :: snowh ! Snow depth (liquid water equivalent)
613 :
614 : real(r8) :: rhminl ! Critical relative humidity for low-level liquid stratus
615 : real(r8) :: rhminl_adj_land ! Adjustment drop of rhminl over the land
616 : real(r8) :: rhminh ! Critical relative humidity for high-level liquid stratus
617 :
618 : real(r8) :: a ! Stratus fraction
619 : real(r8) :: Ga ! dU/da
620 :
621 : ! Local variables
622 : integer i
623 : real(r8) rhmin ! Critical RH
624 : real(r8) rhdif ! Factor for stratus fraction
625 : real(r8) rhwght
626 :
627 : ! Statement functions
628 : logical land
629 : land(i) = nint(landfrac_in(i)) == 1
630 :
631 0 : rhminl = rhminl_const
632 0 : rhminl_adj_land = rhminl_adj_land_const
633 0 : rhminh = rhminh_const
634 :
635 : ! ---------------- !
636 : ! Main computation !
637 : ! ---------------- !
638 :
639 0 : a_out(:) = 0._r8
640 0 : Ga_out(:) = 0._r8
641 :
642 0 : do i = 1, ncol
643 :
644 0 : U = U_in(i)
645 0 : p = p_in(i)
646 0 : qv = qv_in(i)
647 0 : landfrac = landfrac_in(i)
648 0 : snowh = snowh_in(i)
649 :
650 0 : if (present(rhminl_in)) rhminl = rhminl_in(i)
651 0 : if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i)
652 0 : if (present(rhminh_in)) rhminh = rhminh_in(i)
653 :
654 0 : if( p .ge. premib ) then
655 :
656 0 : if( land(i) .and. (snowh.le.0.000001_r8) ) then
657 0 : rhmin = rhminl - rhminl_adj_land
658 : else
659 : rhmin = rhminl
660 : endif
661 0 : rhdif = (U-rhmin)/(1.0_r8-rhmin)
662 0 : a = min(1._r8,(max(rhdif,0.0_r8))**2)
663 0 : if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
664 : Ga = 1.e20_r8
665 : else
666 0 : Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
667 : endif
668 : if( freeze_dry ) then
669 : a = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
670 : Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
671 : endif
672 :
673 0 : elseif( p .lt. premit ) then
674 :
675 0 : rhmin = rhminh
676 0 : rhdif = (U-rhmin)/(1.0_r8-rhmin)
677 0 : a = min(1._r8,(max(rhdif,0._r8))**2)
678 0 : if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
679 : Ga = 1.e20_r8
680 : else
681 0 : Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
682 : endif
683 :
684 : else
685 :
686 0 : rhwght = (premib-(max(p,premit)))/(premib-premit)
687 :
688 : ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
689 : ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
690 : ! else
691 0 : rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
692 : ! endif
693 :
694 0 : rhdif = (U-rhmin)/(1.0_r8-rhmin)
695 0 : a = min(1._r8,(max(rhdif,0._r8))**2)
696 0 : if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
697 : Ga = 1.e10_r8
698 : else
699 0 : Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
700 : endif
701 :
702 : endif
703 :
704 0 : a_out(i) = a
705 0 : Ga_out(i) = Ga
706 :
707 : enddo
708 :
709 0 : end subroutine astG_RHU
710 :
711 : !================================================================================================
712 :
713 0 : subroutine aist_single(qv, T, p, qi, landfrac, snowh, aist, &
714 : rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, &
715 : qsatfac_out)
716 :
717 : ! --------------------------------------------------------- !
718 : ! Compute non-physical ice stratus fraction !
719 : ! --------------------------------------------------------- !
720 :
721 : real(r8), intent(in) :: qv ! Grid-mean water vapor[kg/kg]
722 : real(r8), intent(in) :: T ! Temperature
723 : real(r8), intent(in) :: p ! Pressure [Pa]
724 : real(r8), intent(in) :: qi ! Grid-mean ice water content [kg/kg]
725 : real(r8), intent(in) :: landfrac ! Land fraction
726 : real(r8), intent(in) :: snowh ! Snow depth (liquid water equivalent)
727 :
728 : real(r8), intent(out) :: aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
729 :
730 : real(r8), optional, intent(in) :: rhmaxi_in
731 : real(r8), optional, intent(in) :: rhmini_in ! Critical relative humidity for ice stratus
732 : real(r8), optional, intent(in) :: rhminl_in ! Critical relative humidity for low-level liquid stratus
733 : real(r8), optional, intent(in) :: rhminl_adj_land_in ! Adjustment drop of rhminl over the land
734 : real(r8), optional, intent(in) :: rhminh_in ! Critical relative humidity for high-level liquid stratus
735 : real(r8), optional, intent(out) :: qsatfac_out ! Subgrid scaling factor for qsat
736 :
737 : ! Local variables
738 : real(r8) rhmin ! Critical RH
739 : real(r8) rhwght
740 :
741 : real(r8) a,b,c,as,bs,cs ! Fit parameters
742 : real(r8) Kc ! Constant for ice cloud calc (wood & field)
743 : real(r8) ttmp ! Limited temperature
744 : real(r8) icicval ! Empirical IWC value [ kg/kg ]
745 : real(r8) rho ! Local air density
746 : real(r8) esl ! Liq sat vapor pressure
747 : real(r8) esi ! Ice sat vapor pressure
748 : real(r8) ncf,phi ! Wilson and Ballard parameters
749 : real(r8) es, qs
750 :
751 : real(r8) rhi ! grid box averaged relative humidity over ice
752 : real(r8) minice ! minimum grid box avg ice for having a 'cloud'
753 : real(r8) mincld ! minimum ice cloud fraction threshold
754 : real(r8) icimr ! in cloud ice mixing ratio
755 : real(r8) rhdif ! working variable for slingo scheme
756 :
757 : real(r8) :: rhmaxi
758 : real(r8) :: rhmini
759 : real(r8) :: rhminl
760 : real(r8) :: rhminl_adj_land
761 : real(r8) :: rhminh
762 :
763 : ! Statement functions
764 : logical land
765 0 : land = nint(landfrac) == 1
766 :
767 : ! --------- !
768 : ! Constants !
769 : ! --------- !
770 :
771 : ! Wang and Sassen IWC paramters ( Option.1 )
772 0 : a = 26.87_r8
773 0 : b = 0.569_r8
774 0 : c = 0.002892_r8
775 : ! Schiller parameters ( Option.2 )
776 0 : as = -68.4202_r8
777 0 : bs = 0.983917_r8
778 0 : cs = 2.81795_r8
779 : ! Wood and Field parameters ( Option.3 )
780 0 : Kc = 75._r8
781 : ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds)
782 : ! Slingo modified (option 5)
783 0 : minice = 1.e-12_r8
784 0 : mincld = 1.e-4_r8
785 :
786 0 : rhmaxi = rhmaxi_const
787 0 : if (present(rhmaxi_in)) rhmaxi = rhmaxi_in
788 0 : rhmini = rhmini_const
789 0 : if (present(rhmini_in)) rhmini = rhmini_in
790 0 : rhminl = rhminl_const
791 0 : if (present(rhminl_in)) rhminl = rhminl_in
792 0 : rhminl_adj_land = rhminl_adj_land_const
793 0 : if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in
794 0 : rhminh = rhminh_const
795 0 : if (present(rhminh_in)) rhminh = rhminh_in
796 0 : if (present(qsatfac_out)) qsatfac_out = 1.0_r8
797 :
798 :
799 : ! ---------------- !
800 : ! Main computation !
801 : ! ---------------- !
802 :
803 0 : call qsat_water(T, p, es, qs)
804 0 : esl = svp_water(T)
805 0 : esi = svp_ice(T)
806 :
807 0 : if( iceopt.lt.3 ) then
808 0 : if( iceopt.eq.1 ) then
809 0 : ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8
810 0 : icicval = a + b * ttmp + c * ttmp**2._r8
811 0 : rho = p/(rair*T)
812 0 : icicval = icicval * 1.e-6_r8 / rho
813 : else
814 0 : ttmp = max(190._r8,min(T,273.16_r8))
815 0 : icicval = 10._r8 **(as * bs**ttmp + cs)
816 0 : icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
817 : endif
818 0 : aist = max(0._r8,min(qi/icicval,1._r8))
819 0 : elseif( iceopt.eq.3 ) then
820 0 : aist = 1._r8 - exp(-Kc*qi/(qs*(esi/esl)))
821 0 : aist = max(0._r8,min(aist,1._r8))
822 0 : elseif( iceopt.eq.4) then
823 0 : if( p .ge. premib ) then
824 0 : if( land .and. (snowh.le.0.000001_r8) ) then
825 : rhmin = rhminl - rhminl_adj_land
826 : else
827 : rhmin = rhminl
828 : endif
829 0 : elseif( p .lt. premit ) then
830 : rhmin = rhminh
831 : else
832 0 : rhwght = (premib-(max(p,premit)))/(premib-premit)
833 : ! if( land .and. (snowh.le.0.000001_r8) ) then
834 : ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
835 : ! else
836 0 : rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
837 : ! endif
838 : endif
839 0 : ncf = qi/((1._r8 - icecrit)*qs)
840 0 : if( ncf.le.0._r8 ) then
841 0 : aist = 0._r8
842 0 : elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then
843 0 : aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
844 0 : elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then
845 0 : phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
846 0 : aist = (1._r8 - 4._r8 * cos(phi) * cos(phi))
847 : else
848 0 : aist = 1._r8
849 : endif
850 0 : aist = max(0._r8,min(aist,1._r8))
851 0 : elseif (iceopt.eq.5) then
852 : ! set rh ice cloud fraction
853 0 : rhi= (qv+qi)/qs * (esl/esi)
854 0 : if (rhmaxi .eq. rhmini) then
855 0 : if (rhi .gt. rhmini) then
856 : rhdif = 1._r8
857 : else
858 0 : rhdif = 0._r8
859 : end if
860 : else
861 0 : rhdif = (rhi-rhmini) / (rhmaxi - rhmini)
862 : end if
863 0 : aist = min(1.0_r8, max(rhdif,0._r8)**2)
864 :
865 : ! Similar to alpha in Wilson & Ballard (1999), determine a
866 : ! scaling factor for saturation vapor pressure that reflects
867 : ! the cloud fraction, rhmini, and rhmaxi.
868 : !
869 : ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less.
870 0 : if (present(qsatfac_out) .and. cldfrc2m_do_subgrid_growth) then
871 0 : qsatfac_out = max(min(qv / qs, 1._r8), (1._r8 - aist) * rhmini + aist * rhmaxi)
872 : end if
873 :
874 : ! limiter to remove empty cloud and ice with no cloud
875 : ! and set icecld fraction to mincld if ice exists
876 :
877 0 : if (qi.lt.minice) then
878 0 : aist=0._r8
879 : else
880 0 : aist=max(mincld,aist)
881 : endif
882 :
883 : ! enforce limits on icimr
884 0 : if (qi.ge.minice) then
885 0 : icimr=qi/aist
886 :
887 : !minimum
888 0 : if (icimr.lt.cldfrc2m_qist_min) then
889 0 : if (cldfrc2m_do_avg_aist_algs) then
890 : !
891 : ! Take the geometric mean of the iceopt=4 and iceopt=5 values.
892 : ! Mods developed by Thomas Toniazzo for NorESM.
893 0 : aist = max(0._r8,min(1._r8,sqrt(aist*qi/cldfrc2m_qist_min)))
894 : else
895 : !
896 : ! Default for iceopt=5
897 0 : aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_min))
898 : end if
899 : endif
900 : !maximum
901 0 : if (icimr.gt.cldfrc2m_qist_max) then
902 0 : aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_max))
903 : endif
904 :
905 : endif
906 : endif
907 :
908 : ! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate
909 : ! computed after updating 'qi_st'.
910 :
911 0 : aist = max(0._r8,min(aist,0.999_r8))
912 :
913 0 : end subroutine aist_single
914 :
915 : !================================================================================================
916 :
917 375272352 : subroutine aist_vector(qv_in, T_in, p_in, qi_in, ni_in, landfrac_in, snowh_in, aist_out, ncol, &
918 : rhmaxi_in, rhmini_in, rhminl_in, rhminl_adj_land_in, rhminh_in, &
919 : qsatfac_out )
920 :
921 : ! --------------------------------------------------------- !
922 : ! Compute non-physical ice stratus fraction !
923 : ! --------------------------------------------------------- !
924 :
925 : real(r8), intent(in) :: qv_in(pcols) ! Grid-mean water vapor[kg/kg]
926 : real(r8), intent(in) :: T_in(pcols) ! Temperature
927 : real(r8), intent(in) :: p_in(pcols) ! Pressure [Pa]
928 : real(r8), intent(in) :: qi_in(pcols) ! Grid-mean ice water content [kg/kg]
929 : real(r8), intent(in) :: ni_in(pcols) ! Grid-mean ice water number concentration [#/kg]
930 : real(r8), intent(in) :: landfrac_in(pcols) ! Land fraction
931 : real(r8), intent(in) :: snowh_in(pcols) ! Snow depth (liquid water equivalent)
932 :
933 : real(r8), intent(out) :: aist_out(pcols) ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
934 : integer, intent(in) :: ncol
935 :
936 : real(r8), optional, intent(in) :: rhmaxi_in(pcols)
937 : real(r8), optional, intent(in) :: rhmini_in(pcols) ! Critical relative humidity for ice stratus
938 : real(r8), optional, intent(in) :: rhminl_in(pcols) ! Critical relative humidity for low-level liquid stratus
939 : real(r8), optional, intent(in) :: rhminl_adj_land_in(pcols) ! Adjustment drop of rhminl over the land
940 : real(r8), optional, intent(in) :: rhminh_in(pcols) ! Critical relative humidity for high-level liquid stratus
941 : real(r8), optional, intent(out) :: qsatfac_out(pcols) ! Subgrid scaling factor for qsat
942 :
943 : ! Local variables
944 :
945 : real(r8) qv ! Grid-mean water vapor[kg/kg]
946 : real(r8) T ! Temperature
947 : real(r8) p ! Pressure [Pa]
948 : real(r8) qi ! Grid-mean ice water content [kg/kg]
949 : real(r8) ni
950 : real(r8) landfrac ! Land fraction
951 : real(r8) snowh ! Snow depth (liquid water equivalent)
952 :
953 : real(r8) rhmaxi ! Critical relative humidity for ice stratus
954 : real(r8) rhmini ! Critical relative humidity for ice stratus
955 : real(r8) rhminl ! Critical relative humidity for low-level liquid stratus
956 : real(r8) rhminl_adj_land ! Adjustment drop of rhminl over the land
957 : real(r8) rhminh ! Critical relative humidity for high-level liquid stratus
958 :
959 : real(r8) aist ! Non-physical ice stratus fraction ( 0<= aist <= 1 )
960 :
961 : real(r8) rhmin ! Critical RH
962 : real(r8) rhwght
963 :
964 : real(r8) a,b,c,as,bs,cs,ah,bh,ch ! Fit parameters
965 : real(r8) nil
966 : real(r8) Kc ! Constant for ice cloud calc (wood & field)
967 : real(r8) ttmp ! Limited temperature
968 : real(r8) icicval ! Empirical IWC value [ kg/kg ]
969 : real(r8) rho ! Local air density
970 : real(r8) esl(pcols) ! Liq sat vapor pressure
971 : real(r8) esi(pcols) ! Ice sat vapor pressure
972 : real(r8) ncf,phi ! Wilson and Ballard parameters
973 : real(r8) qs
974 : real(r8) esat_in(pcols)
975 : real(r8) qsat_in(pcols)
976 :
977 : real(r8) rhi ! grid box averaged relative humidity over ice
978 : real(r8) minice ! minimum grid box avg ice for having a 'cloud'
979 : real(r8) mincld ! minimum ice cloud fraction threshold
980 : real(r8) icimr ! in cloud ice mixing ratio
981 : real(r8) rhdif ! working variable for slingo scheme
982 :
983 : integer i
984 :
985 :
986 : ! Statement functions
987 : logical land
988 : land(i) = nint(landfrac_in(i)) == 1
989 :
990 : ! --------- !
991 : ! Constants !
992 : ! --------- !
993 :
994 : ! Wang and Sassen IWC paramters ( Option.1 )
995 375272352 : a = 26.87_r8
996 375272352 : b = 0.569_r8
997 375272352 : c = 0.002892_r8
998 : ! Schiller parameters ( Option.2 )
999 375272352 : as = -68.4202_r8
1000 375272352 : bs = 0.983917_r8
1001 375272352 : cs = 2.81795_r8
1002 : ! Wood and Field parameters ( Option.3 )
1003 375272352 : Kc = 75._r8
1004 : ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds)
1005 : ! Slingo modified (option 5)
1006 375272352 : minice = 1.e-12_r8
1007 375272352 : mincld = 1.e-4_r8
1008 :
1009 375272352 : rhmaxi = rhmaxi_const
1010 :
1011 375272352 : rhmini = rhmini_const
1012 375272352 : rhminl = rhminl_const
1013 375272352 : rhminl_adj_land = rhminl_adj_land_const
1014 375272352 : rhminh = rhminh_const
1015 :
1016 6379629984 : if (present(qsatfac_out)) qsatfac_out = 1.0_r8
1017 :
1018 : ! ---------------- !
1019 : ! Main computation !
1020 : ! ---------------- !
1021 :
1022 375272352 : aist_out(:) = 0._r8
1023 375272352 : esat_in(:) = 0._r8
1024 375272352 : qsat_in(:) = 0._r8
1025 :
1026 375272352 : call qsat_water(T_in(1:ncol), p_in(1:ncol), esat_in(1:ncol), qsat_in(1:ncol), ncol)
1027 375272352 : call svp_water_vect(T_in(1:ncol), esl(1:ncol), ncol)
1028 375272352 : call svp_ice_vect(T_in(1:ncol), esi(1:ncol), ncol)
1029 :
1030 6266175552 : do i = 1, ncol
1031 :
1032 5890903200 : landfrac = landfrac_in(i)
1033 5890903200 : snowh = snowh_in(i)
1034 5890903200 : T = T_in(i)
1035 5890903200 : qv = qv_in(i)
1036 5890903200 : p = p_in(i)
1037 5890903200 : qi = qi_in(i)
1038 5890903200 : ni = ni_in(i)
1039 5890903200 : qs = qsat_in(i)
1040 :
1041 5890903200 : if (present(rhmaxi_in)) rhmaxi = rhmaxi_in(i)
1042 5890903200 : if (present(rhmini_in)) rhmini = rhmini_in(i)
1043 5890903200 : if (present(rhminl_in)) rhminl = rhminl_in(i)
1044 5890903200 : if (present(rhminl_adj_land_in)) rhminl_adj_land = rhminl_adj_land_in(i)
1045 5890903200 : if (present(rhminh_in)) rhminh = rhminh_in(i)
1046 :
1047 5890903200 : if( iceopt.lt.3 ) then
1048 0 : if( iceopt.eq.1 ) then
1049 0 : ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8
1050 0 : icicval = a + b * ttmp + c * ttmp**2._r8
1051 0 : rho = p/(rair*T)
1052 0 : icicval = icicval * 1.e-6_r8 / rho
1053 : else
1054 0 : ttmp = max(190._r8,min(T,273.16_r8))
1055 0 : icicval = 10._r8 **(as * bs**ttmp + cs)
1056 0 : icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
1057 : endif
1058 0 : aist = max(0._r8,min(qi/icicval,1._r8))
1059 5890903200 : elseif( iceopt.eq.3 ) then
1060 0 : aist = 1._r8 - exp(-Kc*qi/(qs*(esi(i)/esl(i))))
1061 0 : aist = max(0._r8,min(aist,1._r8))
1062 5890903200 : elseif( iceopt.eq.4) then
1063 0 : if( p .ge. premib ) then
1064 0 : if( land(i) .and. (snowh.le.0.000001_r8) ) then
1065 : rhmin = rhminl - rhminl_adj_land
1066 : else
1067 : rhmin = rhminl
1068 : endif
1069 0 : elseif( p .lt. premit ) then
1070 : rhmin = rhminh
1071 : else
1072 0 : rhwght = (premib-(max(p,premit)))/(premib-premit)
1073 : ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
1074 : ! rhmin = rhminh*rhwght + (rhminl - rhminl_adj_land)*(1.0_r8-rhwght)
1075 : ! else
1076 0 : rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
1077 : ! endif
1078 : endif
1079 0 : ncf = qi/((1._r8 - icecrit)*qs)
1080 0 : if( ncf.le.0._r8 ) then
1081 : aist = 0._r8
1082 0 : elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then
1083 0 : aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
1084 0 : elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then
1085 0 : phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
1086 0 : aist = (1._r8 - 4._r8 * cos(phi) * cos(phi))
1087 : else
1088 : aist = 1._r8
1089 : endif
1090 0 : aist = max(0._r8,min(aist,1._r8))
1091 5890903200 : elseif (iceopt.eq.5) then
1092 : ! set rh ice cloud fraction
1093 5890903200 : rhi= (qv+qi)/qs * (esl(i)/esi(i))
1094 5890903200 : if (rhmaxi .eq. rhmini) then
1095 3033111842 : if (rhi .gt. rhmini) then
1096 : rhdif = 1._r8
1097 : else
1098 3014101397 : rhdif = 0._r8
1099 : end if
1100 : else
1101 2857791358 : rhdif = (rhi-rhmini) / (rhmaxi - rhmini)
1102 : end if
1103 5890903200 : aist = min(1.0_r8, max(rhdif,0._r8)**2)
1104 :
1105 0 : elseif (iceopt.eq.6) then
1106 : !----- ICE CLOUD OPTION 6: fit based on T and Number (Gettelman: based on Heymsfield obs)
1107 : ! Use observations from Heymsfield et al 2012 of IWC and Ni v. Temp
1108 : ! Multivariate fit follows form of Boudala 2002: ICIWC = a * exp(b*T) * N^c
1109 : ! a=6.73e-8, b=0.05, c=0.349
1110 : ! N is #/L, so need to convert Ni_L=N*rhoa/1000.
1111 0 : ah= 6.73834e-08_r8
1112 0 : bh= 0.0533110_r8
1113 0 : ch= 0.3493813_r8
1114 0 : rho=p/(rair*T)
1115 0 : nil=ni*rho/1000._r8
1116 0 : icicval = ah * exp(bh*T) * nil**ch
1117 : !result is in g m-3, convert to kg H2O / kg air (icimr...)
1118 0 : icicval = icicval / rho / 1000._r8
1119 0 : aist = max(0._r8,min(qi/icicval,1._r8))
1120 0 : aist = min(aist,1._r8)
1121 :
1122 : endif
1123 :
1124 5890903200 : if (iceopt.eq.5 .or. iceopt.eq.6) then
1125 :
1126 : ! Similar to alpha in Wilson & Ballard (1999), determine a
1127 : ! scaling factor for saturation vapor pressure that reflects
1128 : ! the cloud fraction, rhmini, and rhmaxi.
1129 : !
1130 : ! NOTE: Limit qsatfac so that adjusted RHliq would be 1. or less.
1131 5890903200 : if (present(qsatfac_out) .and. cldfrc2m_do_subgrid_growth) then
1132 5890903200 : qsatfac_out(i) = max(min(qv / qs, 1._r8), (1._r8 - aist) * rhmini + aist * rhmaxi)
1133 : end if
1134 :
1135 : ! limiter to remove empty cloud and ice with no cloud
1136 : ! and set icecld fraction to mincld if ice exists
1137 :
1138 5890903200 : if (qi.lt.minice) then
1139 : aist=0._r8
1140 : else
1141 1070716011 : aist=max(mincld,aist)
1142 : endif
1143 :
1144 : ! enforce limits on icimr
1145 5890903200 : if (qi.ge.minice) then
1146 1070716011 : icimr=qi/aist
1147 :
1148 : !minimum
1149 1070716011 : if (icimr.lt.cldfrc2m_qist_min) then
1150 171125593 : if (cldfrc2m_do_avg_aist_algs) then
1151 : !
1152 : ! Take the geometric mean of the iceopt=4 and iceopt=5 values.
1153 : ! Mods developed by Thomas Toniazzo for NorESM.
1154 0 : aist = max(0._r8,min(1._r8,sqrt(aist*qi/cldfrc2m_qist_min)))
1155 : else
1156 : !
1157 : ! Default for iceopt=5
1158 171125593 : aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_min))
1159 : end if
1160 : endif
1161 : !maximum
1162 1070716011 : if (icimr.gt.cldfrc2m_qist_max) then
1163 27542498 : aist = max(0._r8,min(1._r8,qi/cldfrc2m_qist_max))
1164 : endif
1165 :
1166 : endif
1167 : endif
1168 :
1169 : ! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate
1170 : ! computed after updating 'qi_st'.
1171 :
1172 5890903200 : aist = max(0._r8,min(aist,0.999_r8))
1173 :
1174 6266175552 : aist_out(i) = aist
1175 :
1176 : enddo
1177 :
1178 375272352 : end subroutine aist_vector
1179 :
1180 : !================================================================================================
1181 :
1182 : end module cldfrc2m
|
