Line data Source code
1 : module subcol_SILHS
2 : !---------------------------------------------------------------------------
3 : ! Purpose:
4 : !
5 : ! Implement a subcolumn scheme based on the Subgrid Importance Latin Hypercube
6 : ! Sampling (SILHS) functionality of the CLUBB moist turbulence parameterization.
7 : !
8 : !---------------------------------------------------------------------------
9 :
10 : use shr_kind_mod, only: r8=>shr_kind_r8, r4=>shr_kind_r4, i4=>shr_kind_i4
11 : use physics_types, only: physics_state, physics_tend, physics_ptend
12 : use ppgrid, only: pcols, psubcols, pver, pverp, begchunk, endchunk
13 : use constituents, only: pcnst, cnst_get_ind
14 : use cam_abortutils, only: endrun
15 : use cam_logfile, only: iulog
16 : use cam_history, only: addfld, add_default, outfld, horiz_only
17 : use ref_pres, only: top_lev => trop_cloud_top_lev
18 : #ifdef CLUBB_SGS
19 : #ifdef SILHS
20 : use clubb_intr, only: &
21 : clubb_config_flags, &
22 : clubb_params_single_col, &
23 : stats_metadata, &
24 : stats_zt, stats_zm, stats_sfc, &
25 : pdf_params_chnk, &
26 : hm_metadata, &
27 : hydromet_dim, &
28 : pdf_dim
29 :
30 : use clubb_api_module, only: &
31 : hmp2_ip_on_hmm2_ip_slope_type, &
32 : hmp2_ip_on_hmm2_ip_intrcpt_type, &
33 : precipitation_fractions, &
34 : stats, &
35 : core_rknd
36 :
37 : use silhs_api_module, only: &
38 : silhs_config_flags_type
39 : #endif
40 : #endif
41 : use physconst, only: cpair, gravit, latvap, latice, rair, rga, cappa
42 :
43 : implicit none
44 : private
45 : save
46 :
47 : public :: subcol_register_SILHS !
48 : public :: subcol_init_SILHS ! Initialize
49 : public :: subcol_gen_SILHS ! Generate subcolumn fields by calling SILHS
50 : public :: subcol_readnl_SILHS ! SILHS namelist reader
51 : public :: subcol_ptend_avg_SILHS
52 : public :: subcol_SILHS_var_covar_driver
53 : public :: subcol_SILHS_fill_holes_conserv
54 : public :: subcol_SILHS_hydromet_conc_tend_lim
55 : public :: init_state_subcol
56 : private :: fill_holes_sedimentation
57 : private :: fill_holes_same_phase_vert
58 : #ifdef SILHS
59 : ! Calc subcol mean ! Calc subcol variance
60 : private :: meansc
61 : private :: stdsc
62 :
63 : type (stats), target :: stats_lh_zt, &
64 : stats_lh_sfc
65 : !$omp threadprivate(stats_lh_zt, stats_lh_sfc)
66 :
67 : real( kind = core_rknd ), dimension(:,:), allocatable :: &
68 : corr_array_n_cloud, &
69 : corr_array_n_below
70 :
71 : #endif
72 :
73 : !-----
74 : ! Private module vars
75 : !-----
76 :
77 : ! constituent indicies
78 : integer :: &
79 : ixq = 0, &
80 : ixcldliq = 0, &
81 : ixnumliq = 0, &
82 : ixcldice = 0, &
83 : ixnumice = 0, &
84 : ixrain = 0, &
85 : ixnumrain= 0, &
86 : ixsnow = 0, &
87 : ixnumsnow= 0
88 :
89 : ! Pbuf indicies
90 : integer :: thlm_idx, rtm_idx, ice_supersat_idx, &
91 : alst_idx, cld_idx, qrain_idx, qsnow_idx, &
92 : nrain_idx, nsnow_idx, ztodt_idx, tke_idx, kvh_idx, &
93 : prec_pcw_idx, snow_pcw_idx, prec_str_idx, snow_str_idx, &
94 : qcsedten_idx, qrsedten_idx, qisedten_idx, qssedten_idx, &
95 : vtrmc_idx, umr_idx, vtrmi_idx, ums_idx, qcsevap_idx, qisevap_idx
96 :
97 : logical :: subcol_SILHS_weight ! if set, sets up weights for averaging subcolumns for SILHS
98 : integer :: subcol_SILHS_numsubcol ! number of subcolumns for this run
99 : logical :: docldfracscaling = .false. ! Weight tendencies by cloud fraction
100 :
101 : character(len=256) :: subcol_SILHS_corr_file_path
102 : character(len=16) :: subcol_SILHS_corr_file_name
103 :
104 : logical :: subcol_SILHS_q_to_micro, &
105 : subcol_SILHS_n_to_micro, &
106 : subcol_SILHS_use_clear_col, &
107 : subcol_SILHS_meanice, &
108 : subcol_SILHS_constrainmn
109 :
110 : logical :: subcol_SILHS_var_covar_src
111 :
112 : real(r8) :: subcol_SILHS_ncnp2_on_ncnm2
113 :
114 : ! There may or may not be a better place to put this.
115 : real(r8), parameter :: p0_clubb = 100000._r8
116 :
117 :
118 : ! real(r8) :: subcol_SILHS_c6rt, subcol_SILHS_c7, subcol_SILHS_c8, subcol_SILHS_c11, &
119 : ! subcol_SILHS_c11b, subcol_SILHS_gamma_coef, &
120 : ! subcol_SILHS_mult_coef, subcol_SILHS_mu
121 :
122 : real(r8) :: ztodt ! model timestep
123 : #ifdef CLUBB_SGS
124 : #ifdef SILHS
125 : type(hmp2_ip_on_hmm2_ip_slope_type) :: subcol_SILHS_hmp2_ip_on_hmm2_ip_slope
126 : type(hmp2_ip_on_hmm2_ip_intrcpt_type) :: subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt
127 :
128 : type(silhs_config_flags_type) :: silhs_config_flags
129 : #endif
130 : #endif
131 :
132 : contains
133 :
134 0 : subroutine subcol_register_SILHS()
135 :
136 : !--------------------------------
137 : ! Register fields needed by SILHS in the physics buffer
138 : ! Currently, most fields needed by SILHS but calculated by CLUBB are registered
139 : ! by clubb in clubb_intr.F90.
140 : !
141 : !--------------------------------
142 : #ifdef CLUBB_SGS
143 : #ifdef SILHS
144 : #endif
145 : #endif
146 0 : end subroutine subcol_register_SILHS
147 :
148 :
149 0 : subroutine subcol_readnl_SILHS(nlfile)
150 : #ifdef CLUBB_SGS
151 : #ifdef SILHS
152 : use namelist_utils, only: find_group_name
153 : use spmd_utils, only: masterproc, masterprocid, mpicom
154 : use spmd_utils, only: mpi_integer, mpi_logical, mpi_character, mpir8, iam
155 : use clubb_api_module, only: core_rknd
156 : use silhs_api_module, only: set_default_silhs_config_flags_api, &
157 : initialize_silhs_config_flags_type_api, &
158 : print_silhs_config_flags_api
159 : #endif
160 : #endif
161 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
162 :
163 : ! Local variables
164 : integer :: unitn, ierr
165 : #ifdef CLUBB_SGS
166 : #ifdef SILHS
167 :
168 : integer :: &
169 : cluster_allocation_strategy
170 :
171 : logical :: &
172 : subcol_silhs_l_lh_importance_sampling, &
173 : subcol_silhs_l_Lscale_vert_avg, &
174 : subcol_silhs_l_lh_straight_mc, &
175 : subcol_silhs_l_lh_clustered_sampling, &
176 : subcol_silhs_l_rcm_in_cloud_k_lh_start, &
177 : subcol_silhs_l_random_k_lh_start, &
178 : subcol_silhs_l_max_overlap_in_cloud, &
179 : subcol_silhs_l_lh_instant_var_covar_src, &
180 : subcol_silhs_l_lh_limit_weights, &
181 : subcol_silhs_l_lh_var_frac, &
182 : subcol_silhs_l_lh_normalize_weights
183 :
184 : namelist /subcol_SILHS_nl/ subcol_SILHS_weight, &
185 : subcol_SILHS_numsubcol, &
186 : subcol_SILHS_corr_file_path, &
187 : subcol_SILHS_corr_file_name, &
188 : subcol_SILHS_q_to_micro, &
189 : subcol_SILHS_n_to_micro, &
190 : subcol_SILHS_ncnp2_on_ncnm2, &
191 : subcol_SILHS_hmp2_ip_on_hmm2_ip_slope, &
192 : subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt, &
193 : subcol_SILHS_meanice, &
194 : subcol_SILHS_use_clear_col, &
195 : subcol_SILHS_constrainmn, &
196 : subcol_SILHS_var_covar_src
197 : ! subcol_SILHS_c6rt, subcol_SILHS_c7, &
198 : ! subcol_SILHS_c8, subcol_SILHS_c11, subcol_SILHS_c11b, &
199 : ! subcol_SILHS_gamma_coef, subcol_SILHS_mult_coef, subcol_SILHS_mu
200 :
201 : namelist /silhs_config_flags_nl/ subcol_silhs_l_lh_importance_sampling, &
202 : subcol_silhs_l_Lscale_vert_avg, &
203 : subcol_silhs_l_lh_straight_mc, &
204 : subcol_silhs_l_lh_clustered_sampling, &
205 : subcol_silhs_l_rcm_in_cloud_k_lh_start, &
206 : subcol_silhs_l_random_k_lh_start, &
207 : subcol_silhs_l_max_overlap_in_cloud, &
208 : subcol_silhs_l_lh_instant_var_covar_src, &
209 : subcol_silhs_l_lh_limit_weights, &
210 : subcol_silhs_l_lh_var_frac, &
211 : subcol_silhs_l_lh_normalize_weights
212 :
213 : !-----------------------------------------------------------------------------
214 : ! Set defaults
215 :
216 : ! Eric Raut changed a default.
217 : subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%Ni = 0.0_core_rknd
218 : subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%Ni = 0.5_core_rknd
219 :
220 : if (masterproc) then
221 : open( newunit=unitn, file=trim(nlfile), status='old' )
222 : call find_group_name(unitn, 'subcol_SILHS_nl', status=ierr)
223 : if (ierr == 0) then
224 : read(unitn, subcol_SILHS_nl, iostat=ierr)
225 : if (ierr /= 0) then
226 : call endrun('subcol_readnl_SILHS: ERROR reading namelist')
227 : end if
228 : end if
229 : close(unitn)
230 : end if
231 :
232 : ! Set default silhs_config_flags entires
233 : call set_default_silhs_config_flags_api( cluster_allocation_strategy, &
234 : subcol_silhs_l_lh_importance_sampling, &
235 : subcol_silhs_l_Lscale_vert_avg, &
236 : subcol_silhs_l_lh_straight_mc, &
237 : subcol_silhs_l_lh_clustered_sampling, &
238 : subcol_silhs_l_rcm_in_cloud_k_lh_start, &
239 : subcol_silhs_l_random_k_lh_start, &
240 : subcol_silhs_l_max_overlap_in_cloud, &
241 : subcol_silhs_l_lh_instant_var_covar_src, &
242 : subcol_silhs_l_lh_limit_weights, &
243 : subcol_silhs_l_lh_var_frac, &
244 : subcol_silhs_l_lh_normalize_weights )
245 :
246 : ! Get silhs_config_flags entries from namelist
247 : if (masterproc) then
248 : open( newunit=unitn, file=trim(nlfile), status='old' )
249 : call find_group_name(unitn, 'silhs_config_flags_nl', status=ierr)
250 : if (ierr == 0) then
251 : read(unitn, silhs_config_flags_nl, iostat=ierr)
252 : if (ierr /= 0) then
253 : call endrun('silhs_config_flags_nl: ERROR reading namelist')
254 : end if
255 : end if
256 : close(unitn)
257 : end if
258 :
259 : ! Save silhs_config_flags entries into module variable silhs_config_flags
260 : call initialize_silhs_config_flags_type_api( cluster_allocation_strategy, &
261 : subcol_silhs_l_lh_importance_sampling, &
262 : subcol_silhs_l_Lscale_vert_avg, &
263 : subcol_silhs_l_lh_straight_mc, &
264 : subcol_silhs_l_lh_clustered_sampling, &
265 : subcol_silhs_l_rcm_in_cloud_k_lh_start, &
266 : subcol_silhs_l_random_k_lh_start, &
267 : subcol_silhs_l_max_overlap_in_cloud, &
268 : subcol_silhs_l_lh_instant_var_covar_src, &
269 : subcol_silhs_l_lh_limit_weights, &
270 : subcol_silhs_l_lh_var_frac, &
271 : subcol_silhs_l_lh_normalize_weights, &
272 : silhs_config_flags )
273 :
274 : ! Print the SILHS configurable flags
275 : call print_silhs_config_flags_api( iulog, silhs_config_flags ) ! Intent(in)
276 :
277 : #ifdef SPMD
278 : ! Broadcast namelist variables
279 : call mpi_bcast(subcol_SILHS_weight, 1, mpi_logical, masterprocid, mpicom, ierr)
280 : call mpi_bcast(subcol_SILHS_numsubcol, 1, mpi_integer, masterprocid, mpicom, ierr)
281 : call mpi_bcast(subcol_SILHS_corr_file_path, len(subcol_SILHS_corr_file_path), &
282 : mpi_character, masterprocid, mpicom, ierr)
283 : call mpi_bcast(subcol_SILHS_corr_file_name, len(subcol_SILHS_corr_file_name), &
284 : mpi_character, masterprocid, mpicom, ierr)
285 : call mpi_bcast(subcol_SILHS_use_clear_col, 1, mpi_logical, masterprocid, mpicom, ierr)
286 : call mpi_bcast(subcol_SILHS_constrainmn, 1, mpi_logical, masterprocid, mpicom, ierr)
287 : call mpi_bcast(subcol_SILHS_meanice, 1, mpi_logical, masterprocid, mpicom, ierr)
288 : call mpi_bcast(subcol_SILHS_q_to_micro, 1, mpi_logical, masterprocid, mpicom, ierr)
289 : call mpi_bcast(subcol_SILHS_n_to_micro, 1, mpi_logical, masterprocid, mpicom, ierr)
290 : call mpi_bcast(subcol_SILHS_var_covar_src,1,mpi_logical,masterprocid, mpicom, ierr)
291 : call mpi_bcast(subcol_SILHS_ncnp2_on_ncnm2, 1, mpir8, masterprocid, mpicom, ierr)
292 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%rr, 1, mpir8, masterprocid, mpicom, ierr)
293 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%Nr, 1, mpir8, masterprocid, mpicom, ierr)
294 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%ri, 1, mpir8, masterprocid, mpicom, ierr)
295 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%Ni, 1, mpir8, masterprocid, mpicom, ierr)
296 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%rs, 1, mpir8, masterprocid, mpicom, ierr)
297 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_slope%Ns, 1, mpir8, masterprocid, mpicom, ierr)
298 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%rr, 1, mpir8, masterprocid, mpicom, ierr)
299 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%Nr, 1, mpir8, masterprocid, mpicom, ierr)
300 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%ri, 1, mpir8, masterprocid, mpicom, ierr)
301 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%Ni, 1, mpir8, masterprocid, mpicom, ierr)
302 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%rs, 1, mpir8, masterprocid, mpicom, ierr)
303 : call mpi_bcast(subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt%Ns, 1, mpir8, masterprocid, mpicom, ierr)
304 : ! call mpi_bcast(subcol_SILHS_c6rt, 1, mpir8, masterprocid, mpicom, ierr)
305 : ! call mpi_bcast(subcol_SILHS_c7, 1, mpir8, masterprocid, mpicom, ierr)
306 : ! call mpi_bcast(subcol_SILHS_c8, 1, mpir8, masterprocid, mpicom, ierr)
307 : ! call mpi_bcast(subcol_SILHS_c11, 1, mpir8, masterprocid, mpicom, ierr)
308 : ! call mpi_bcast(subcol_SILHS_c11b, 1, mpir8, masterprocid, mpicom, ierr)
309 : ! call mpi_bcast(subcol_SILHS_gamma_coef, 1, mpir8, masterprocid, mpicom, ierr)
310 : ! call mpi_bcast(subcol_SILHS_mult_coef, 1, mpir8, masterprocid, mpicom, ierr)
311 : ! call mpi_bcast(subcol_SILHS_mu, 1, mpir8, masterprocid, mpicom, ierr)
312 : call mpi_bcast(silhs_config_flags%l_lh_importance_sampling, 1, mpi_logical, masterprocid, mpicom, ierr)
313 : call mpi_bcast(silhs_config_flags%l_Lscale_vert_avg, 1, mpi_logical, masterprocid, mpicom, ierr)
314 : call mpi_bcast(silhs_config_flags%l_lh_straight_mc, 1, mpi_logical, masterprocid, mpicom, ierr)
315 : call mpi_bcast(silhs_config_flags%l_lh_clustered_sampling, 1, mpi_logical, masterprocid, mpicom, ierr)
316 : call mpi_bcast(silhs_config_flags%l_rcm_in_cloud_k_lh_start, 1, mpi_logical, masterprocid, mpicom, ierr)
317 : call mpi_bcast(silhs_config_flags%l_random_k_lh_start, 1, mpi_logical, masterprocid, mpicom, ierr)
318 : call mpi_bcast(silhs_config_flags%l_max_overlap_in_cloud, 1, mpi_logical, masterprocid, mpicom, ierr)
319 : call mpi_bcast(silhs_config_flags%l_lh_instant_var_covar_src, 1, mpi_logical, masterprocid, mpicom, ierr)
320 : call mpi_bcast(silhs_config_flags%l_lh_limit_weights, 1, mpi_logical, masterprocid, mpicom, ierr)
321 : call mpi_bcast(silhs_config_flags%l_lh_var_frac, 1, mpi_logical, masterprocid, mpicom, ierr)
322 : call mpi_bcast(silhs_config_flags%l_lh_normalize_weights, 1, mpi_logical, masterprocid, mpicom, ierr)
323 :
324 : ! SPMD
325 : #endif
326 : ! SILHS
327 : #endif
328 : ! CLUBB_SGS
329 : #endif
330 0 : end subroutine subcol_readnl_SILHS
331 :
332 :
333 0 : subroutine subcol_init_SILHS(pbuf2d)
334 :
335 : !--------------------------------
336 : ! Read in parameters and initialize SILHS PDF fields.
337 : ! Set up indexes into Pbuf fields.
338 : ! Register history outputs.
339 : !--------------------------------
340 :
341 : use physics_buffer, only: physics_buffer_desc, pbuf_get_field, &
342 : dtype_r8, pbuf_get_index
343 : #ifdef CLUBB_SGS
344 : #ifdef SILHS
345 : use clubb_api_module, only: core_rknd, &
346 : setup_corr_varnce_array_api, &
347 : init_pdf_hydromet_arrays_api, &
348 : set_clubb_debug_level_api
349 :
350 : #endif
351 : #endif
352 :
353 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
354 :
355 : #ifdef CLUBB_SGS
356 : #ifdef SILHS
357 :
358 : integer :: iunit = 501 ! Default value, will get iunit from CAM
359 : !character(len=*), parameter :: default_corr_case = "arm_97"
360 : character(len=*), parameter :: &
361 : cloud_file_ext = "_corr_array_cloud.in", & ! File extensions for corr files
362 : below_file_ext = "_corr_array_below.in"
363 : character(len=256) :: corr_file_path_cloud, corr_file_path_below
364 :
365 : ! To set up CLUBB hydromet indices
366 : integer :: &
367 : iirr, & ! Hydrometeor array index for rain water mixing ratio, rr
368 : iirs, & ! Hydrometeor array index for snow mixing ratio, rs
369 : iiri, & ! Hydrometeor array index for ice mixing ratio, ri
370 : iirg, & ! Hydrometeor array index for graupel mixing ratio, rg
371 : iiNr, & ! Hydrometeor array index for rain drop concentration, Nr
372 : iiNs, & ! Hydrometeor array index for snow concentration, Ns
373 : iiNi, & ! Hydrometeor array index for ice concentration, Ni
374 : iiNg ! Hydrometeor array index for graupel concentration, Ng
375 :
376 : integer :: l, ierr=0 ! Loop variable, error check
377 :
378 : ! Set CLUBB's debug level
379 : ! This is called in module clubb_intr; no need to do it here.
380 : ! call set_clubb_debug_level_api( 0 )
381 :
382 : !-------------------------------
383 : ! CLUBB-SILHS Parameters (global module variables)
384 : !-------------------------------
385 : clubb_config_flags%l_fix_w_chi_eta_correlations = .true.
386 : clubb_config_flags%l_diagnose_correlations = .false.
387 : clubb_config_flags%l_calc_w_corr = .false.
388 : ! l_prescribed_avg_deltaz = .false.
389 : clubb_config_flags%l_use_cloud_cover = .false.
390 : clubb_config_flags%l_const_Nc_in_cloud = .true.
391 :
392 : ! Values from the namelist
393 : docldfracscaling = subcol_SILHS_use_clear_col
394 :
395 : ! Namelist "tuning" or set correlations
396 : ! KTC Todo: Move these to a tuning "in" file or into the namelist
397 : ! JHTODO: we might want these on CLUBB's API and ultimatively on a namelist for tuning
398 : ! C6rt = subcol_SILHS_c6rt
399 : ! C7 = subcol_SILHS_c7 ! to all ice clouds
400 : ! C8 = subcol_SILHS_c8
401 : ! C11 = subcol_SILHS_c11
402 : ! C11b = subcol_SILHS_c11b
403 : ! gamma_coef = subcol_SILHS_gamma_coef
404 : ! mult_coef = subcol_SILHS_mult_coef
405 : ! mu = subcol_SILHS_mu
406 :
407 : !call set_clubb_debug_level( 0 ) !#KTCtodo: Add a namelist variable to set debug level
408 :
409 : ! Get constituent indices
410 : call cnst_get_ind('Q', ixq)
411 : call cnst_get_ind('CLDLIQ', ixcldliq)
412 : call cnst_get_ind('NUMLIQ', ixnumliq)
413 : call cnst_get_ind('CLDICE', ixcldice)
414 : call cnst_get_ind('NUMICE', ixnumice)
415 : call cnst_get_ind('RAINQM', ixrain, abort=.false.)
416 : call cnst_get_ind('NUMRAI', ixnumrain, abort=.false.)
417 : call cnst_get_ind('SNOWQM', ixsnow, abort=.false.)
418 : call cnst_get_ind('NUMSNO', ixnumsnow, abort=.false.)
419 :
420 : ! Get physics buffer indexes
421 : thlm_idx = pbuf_get_index('THLM')
422 : rtm_idx = pbuf_get_index('RTM')
423 : cld_idx = pbuf_get_index('CLD')
424 : alst_idx = pbuf_get_index('ALST') ! SILHS expects clubb's cloud_frac liq stratus fraction
425 : ztodt_idx = pbuf_get_index('ZTODT')
426 : ice_supersat_idx = pbuf_get_index('ISS_FRAC')
427 : tke_idx = pbuf_get_index('tke')
428 : kvh_idx = pbuf_get_index('kvh')
429 : prec_pcw_idx = pbuf_get_index('PREC_PCW')
430 : snow_pcw_idx = pbuf_get_index('SNOW_PCW')
431 : prec_str_idx = pbuf_get_index('PREC_STR')
432 : snow_str_idx = pbuf_get_index('SNOW_STR')
433 : qcsedten_idx = pbuf_get_index('QCSEDTEN')
434 : qrsedten_idx = pbuf_get_index('QRSEDTEN')
435 : qisedten_idx = pbuf_get_index('QISEDTEN')
436 : qssedten_idx = pbuf_get_index('QSSEDTEN')
437 : vtrmc_idx = pbuf_get_index('VTRMC')
438 : umr_idx = pbuf_get_index('UMR')
439 : vtrmi_idx = pbuf_get_index('VTRMI')
440 : ums_idx = pbuf_get_index('UMS')
441 : qcsevap_idx = pbuf_get_index('QCSEVAP')
442 : qisevap_idx = pbuf_get_index('QISEVAP')
443 : qrain_idx = pbuf_get_index('QRAIN')
444 : qsnow_idx = pbuf_get_index('QSNOW')
445 : nrain_idx = pbuf_get_index('NRAIN')
446 : nsnow_idx = pbuf_get_index('NSNOW')
447 :
448 : !-------------------------------
449 : ! Set up SILHS hydrometeors #KTCtodo: move microphys specification to config time,
450 : ! Steve wants to set up a microphysics query so I can ask the microphysics
451 : ! scheme which hydrometeors to use. For the future.
452 : !-------------------------------
453 : iirr = 1
454 : iirs = 3
455 : iiri = 5
456 : iirg = -1
457 :
458 : iiNr = 2
459 : iiNs = 4
460 : iiNi = 6
461 : iiNg = -1
462 :
463 : hydromet_dim = 6
464 :
465 : ! Set up pdf indices, hydromet indicies, hydromet arrays, and hydromet variance ratios
466 : call init_pdf_hydromet_arrays_api( 1.0_core_rknd, 1.0_core_rknd, hydromet_dim, & ! intent(in)
467 : iirr, iiNr, iiri, iiNi, & ! intent(in)
468 : iirs, iiNs, iirg, iiNg, & ! intent(in)
469 : subcol_SILHS_ncnp2_on_ncnm2, & ! intent(in)
470 : hm_metadata, pdf_dim, & ! intent(out)
471 : subcol_SILHS_hmp2_ip_on_hmm2_ip_slope, & ! optional(in)
472 : subcol_SILHS_hmp2_ip_on_hmm2_ip_intrcpt ) ! optional(in)
473 :
474 : !-------------------------------
475 : ! Set up hydrometeors and correlation arrays for SILHS
476 : !-------------------------------
477 : allocate( corr_array_n_cloud(pdf_dim,pdf_dim), corr_array_n_below(pdf_dim,pdf_dim), stat=ierr)
478 : if( ierr /= 0 ) call endrun(' subcol_init_SILHS: failed to allocate corr_array fields ')
479 :
480 : corr_file_path_cloud = trim( subcol_SILHS_corr_file_path )//trim( subcol_SILHS_corr_file_name )//cloud_file_ext
481 : corr_file_path_below = trim( subcol_SILHS_corr_file_path )//trim( subcol_SILHS_corr_file_name )//below_file_ext
482 :
483 : call setup_corr_varnce_array_api( corr_file_path_cloud, corr_file_path_below, &
484 : pdf_dim, hm_metadata, newunit(iunit), &
485 : clubb_config_flags%l_fix_w_chi_eta_correlations, & ! In
486 : corr_array_n_cloud, corr_array_n_below )
487 :
488 : !-------------------------------
489 : ! Register output fields from SILHS
490 : !-------------------------------
491 : call addfld('SILHS_NCLD_SCOL', (/'psubcols', 'ilev '/), 'I', 'm^-3', &
492 : 'Subcolumn Cloud Number Concentration', flag_xyfill=.true., fill_value=1.e30_r8)
493 : call addfld('SILHS_NRAIN_SCOL', (/'psubcols', 'ilev '/), 'I', 'm^-3', &
494 : 'Subcolumn Number Concentration of Rain from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
495 : call addfld('SILHS_OMEGA_SCOL', (/'psubcols', 'ilev '/), 'I', 'Pa/s', &
496 : 'Subcolumn vertical pressure velocity', flag_xyfill=.true., fill_value=1.e30_r8)
497 : call addfld('SILHS_RCM_SCOL', (/'psubcols', 'ilev '/), 'I', 'kg/kg', &
498 : 'Subcolumn Cloud Liquid Water from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
499 : call addfld('SILHS_RICLD_SCOL', (/'psubcols', 'ilev '/), 'I', 'kg/kg', &
500 : 'Subcolumn Cloud Ice Water from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
501 : call addfld('SILHS_NICLD_SCOL', (/'psubcols', 'ilev '/), 'I', 'kg/kg', &
502 : 'Subcolumn Cloud Ice Number Conc from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
503 : call addfld('SILHS_RRAIN_SCOL', (/'psubcols', 'ilev '/), 'I', 'kg/kg', &
504 : 'Subcolumn Precipitating Liquid Water from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
505 : call addfld('SILHS_RT_SCOL', (/'psubcols', 'ilev '/), 'I', 'kg/kg ', &
506 : 'Subcolumn Total Water from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
507 : call addfld('SILHS_THLM_SCOL', (/'psubcols', 'ilev '/), 'I', 'K', &
508 : 'Subcolumn liquid water pot temperature', flag_xyfill=.true., fill_value=1.e30_r8)
509 : call addfld('SILHS_WEIGHT_SCOL', (/'psubcols'/), 'I', 'frac', &
510 : 'Weights for each subcolumn', flag_xyfill=.true., fill_value=1.e30_r8)
511 : call addfld('SILHS_WM_SCOL', (/'psubcols', 'ilev '/), 'I', 'm/s', &
512 : 'Subcolumn vertical velocity from SILHS', flag_xyfill=.true., fill_value=1.e30_r8)
513 :
514 : call addfld('NR_IN_LH', (/ 'lev' /), 'I', 'm^-3', &
515 : 'Num Rain Conc as input to SILHS')
516 : call addfld('SILHS_RTM', (/ 'ilev' /), 'I', 'kg/kg', &
517 : 'Input total water mixing ratio')
518 : call addfld('SILHS_THLM', (/ 'ilev' /), 'I', 'K', &
519 : 'Input liquid water potential temperature')
520 : call addfld('SILHS_QC_IN', (/ 'lev' /), 'I', 'kg/kg', &
521 : 'Input cloud water mixing ratio')
522 : call addfld('SILHS_QI_IN', (/ 'lev' /), 'I', 'kg/kg', &
523 : 'Input cloud ice mixing ratio')
524 : call addfld('SILHS_NC_IN', (/ 'lev' /), 'I', '#/kg', &
525 : 'Input cloud water number concentration')
526 : call addfld('SILHS_NI_IN', (/ 'lev' /), 'I', '#/kg', &
527 : 'Input cloud ice number concentration')
528 : call addfld('AKM_CLUBB', (/ 'ilev' /), 'I', '(kg/kg)/s', &
529 : 'Exact Kessler autoconversion')
530 : call addfld('AKM_LH_CLUBB', (/ 'ilev' /), 'I', '(kg/kg)/s', &
531 : 'Monte Carlo estimate of Kessler autoconversion')
532 : call addfld('INVS_EXNER', (/ 'lev' /), 'I', 'none', &
533 : 'inverse EXNER function from state in subcol_SILHS')
534 : call addfld('SILHS_ZTODT', horiz_only, 'I', 's', &
535 : 'Length of Physics timestep (for debugging)')
536 : if ( subcol_SILHS_constrainmn ) then
537 : call addfld('SILHS_MSC_CLDICE', (/ 'lev' /), 'A', 'kg/kg', &
538 : 'Mean Cloud Ice across subcolumns')
539 : call addfld('SILHS_STDSC_CLDICE', (/ 'lev' /), 'A', 'kg/kg', &
540 : 'Standard deviation of Ice across subcolumns')
541 : if ( ixsnow > 0 ) then
542 : call addfld('SILHS_MSC_CLDLIQ', (/ 'lev' /), 'A', 'kg/kg', &
543 : 'Mean Cloud Liquid across subcolumns')
544 : call addfld('SILHS_STDSC_CLDLIQ', (/ 'lev' /), 'A', 'kg/kg', &
545 : 'Standard deviation of Liquid across subcolumns')
546 : call addfld('SILHS_MSC_Q', (/ 'lev' /), 'A', 'kg/kg', &
547 : 'Mean water vapor across subcolumns')
548 : call addfld('SILHS_STDSC_Q', (/ 'lev' /), 'A', 'kg/kg', &
549 : 'Standard deviation of water vapor across subcolumns')
550 : endif ! ixsnow > 0
551 : endif ! subcol_SILHS_constrainmn
552 : call addfld('SILHS_EFF_CLDFRAC', (/ 'lev' /), 'A', 'frac', &
553 : 'Calculated cloud fraction from subcolumn liq or ice')
554 :
555 : call addfld('SILHS_CLUBB_PRECIP_FRAC', (/ 'lev' /), 'A', 'frac', &
556 : 'Precipitation fraction from CLUBB (set_up_pdf_params_incl_hydromet)')
557 : call addfld('SILHS_CLUBB_ICE_SS_FRAC', (/ 'lev' /), 'A', 'frac', &
558 : 'Ice supersaturation fraction from CLUBB')
559 :
560 : call addfld ('QVHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Water vapor mixing ratio tendency from hole filling')
561 : call addfld ('QCHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Cloud water mixing ratio tendency from hole filling')
562 : call addfld ('QRHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Rain water mixing ratio tendency from hole filling')
563 : call addfld ('QIHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Cloud ice mixing ratio tendency from hole filling')
564 : call addfld ('QSHFTEN', (/ 'lev' /), 'A', 'kg/kg/s', 'Snow mixing ratio tendency from hole filling')
565 : call addfld ('THFTEN', (/ 'lev' /), 'A', 'K/s', 'Temperature tendency from hole filling')
566 :
567 : #endif
568 : #endif
569 0 : end subroutine subcol_init_SILHS
570 : !==============================================================!
571 0 : subroutine init_state_subcol(state, tend, state_sc, tend_sc)
572 :
573 0 : use ppgrid, only : pver, pverp, pcols
574 :
575 : use subcol_utils, only : subcol_set_subcols
576 :
577 : implicit none
578 :
579 : type(physics_state), intent(inout) :: state
580 : type(physics_tend), intent(inout) :: tend
581 : type(physics_state), intent(inout) :: state_sc ! sub-column state
582 : type(physics_tend), intent(inout) :: tend_sc ! sub-column tend
583 :
584 : integer, dimension(pcols) :: numsubcol_arr ! To set up the state struct
585 :
586 0 : numsubcol_arr(:) = 0 ! Start over each chunk
587 0 : numsubcol_arr(:state%ngrdcol) = subcol_SILHS_numsubcol ! Only set for valid grid columns
588 0 : call subcol_set_subcols(state, tend, numsubcol_arr, state_sc, tend_sc)
589 :
590 0 : end subroutine init_state_subcol
591 : !==================================================================!
592 0 : subroutine subcol_gen_SILHS(state, tend, state_sc, tend_sc, pbuf)
593 : !-------------------------------
594 : ! This is where the subcolumns are created, and the call to
595 : ! generate_silhs_sample_mod_api
596 : ! goes out. Variables needed to make this call are pulled from the
597 : ! pbuf, from module data, and calculated based on the CAM state.
598 : !-------------------------------
599 :
600 0 : use physics_buffer, only : physics_buffer_desc, pbuf_get_index, &
601 : pbuf_get_field
602 : use time_manager, only : get_nstep
603 : use subcol_utils, only : subcol_set_subcols, subcol_set_weight
604 : use subcol_pack_mod, only : subcol_pack
605 : use phys_control, only : phys_getopts
606 : use spmd_utils, only : masterproc
607 : use shr_const_mod, only : SHR_CONST_PI, SHR_CONST_RHOFW
608 :
609 : #ifdef CLUBB_SGS
610 : #ifdef SILHS
611 : use clubb_api_module, only : setup_pdf_parameters_api, &
612 :
613 : zm2zt_api, setup_grid_heights_api, &
614 :
615 : core_rknd, &
616 :
617 : w_tol_sqd, zero_threshold, &
618 : em_min, cloud_frac_min, & ! rc_tol, &
619 :
620 : genrand_intg, genrand_init_api, &
621 :
622 : nparams, ic_K, &
623 : read_parameters_api, &
624 : Cp, Lv, &
625 : grid, setup_grid_api, &
626 : init_precip_fracs_api
627 :
628 : use silhs_api_module, only : generate_silhs_sample_api, & ! Ncn_to_Nc, &
629 : clip_transform_silhs_output_api, &
630 : est_kessler_microphys_api, &
631 : vert_decorr_coef
632 :
633 : #endif
634 : #endif
635 :
636 : ! CAM data structures
637 : type(physics_state), intent(inout) :: state
638 : type(physics_tend), intent(inout) :: tend
639 : type(physics_state), intent(inout) :: state_sc ! sub-column state
640 : type(physics_tend), intent(inout) :: tend_sc ! sub-column tend
641 : type(physics_buffer_desc), pointer :: pbuf(:)
642 :
643 : #ifdef CLUBB_SGS
644 : #ifdef SILHS
645 : !----------------
646 : ! Local variables
647 : !----------------
648 : logical, parameter :: &
649 : l_implemented = .true. ! Implemented in a host model
650 : logical, parameter :: rx_Nc = .false. ! Use NC calculated based on grid mean effective radius
651 : integer, parameter :: &
652 : grid_type = 3 ! The 3 grid centered on momentum levels
653 : real(r8), parameter :: cldmin = 0.001_r8 ! To use when cld frac = 0.0 to be consistant with micro_mg
654 : real(r8), parameter :: min_num_conc = 1.0e-12_r8
655 : real(r8), parameter :: qsmall = 1.0e-18_r8 ! Microphysics cut-off for cloud
656 :
657 : integer :: i, j, k, ngrdcol, ncol, lchnk, stncol
658 : real(r8) :: sfc_elevation(state%ngrdcol) ! Surface elevation
659 :
660 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: zt_g, zi_g ! Thermo & Momentum grids for clubb
661 :
662 : real(r8), dimension(pverp) :: scfrac ! cloud fraction based on sc distributions
663 : real(r8) :: msc, std, maxcldfrac, maxsccldfrac
664 : real(r8) :: scale = 1.0_r8
665 :
666 : real(r8) :: c_K ! CLUBB parameter c_K (for eddy diffusivity)
667 :
668 : integer( kind = genrand_intg ) :: &
669 : lh_seed ! Seed used in random number generator that will be different
670 : ! for each column, yet reproducible for a restart run
671 :
672 : !----------------
673 : ! Required for set_up_pdf_params_incl_hydromet
674 : !----------------
675 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: cld_frac_in ! Cloud fraction
676 :
677 : real(r8), dimension(state%ngrdcol, pverp-top_lev+1, pdf_dim, pdf_dim) :: &
678 : corr_array_1, corr_array_2 ! Correlation matrix for pdf components
679 :
680 : real(r8), dimension(state%ngrdcol, pverp-top_lev+1, pdf_dim) :: &
681 : mu_x_1, mu_x_2, & ! Mean array for PDF components
682 : sigma_x_1, sigma_x_2 ! Std dev arr for PDF components
683 :
684 : real(r8), dimension(state%ngrdcol, pverp-top_lev+1, pdf_dim, pdf_dim) :: &
685 : corr_cholesky_mtx_1, corr_cholesky_mtx_2 ! Transposed corr cholesky mtx
686 :
687 : real(r8), dimension(state%ngrdcol, pverp-top_lev+1) :: Nc_in_cloud
688 : real(r8), dimension(state%ngrdcol, pverp-top_lev+1) :: ice_supersat_frac_in
689 : real(r8), dimension(state%ngrdcol, pverp-top_lev+1, hydromet_dim) :: hydrometp2
690 :
691 :
692 : !----------------
693 : ! Input to generate_silhs_sample
694 : !----------------
695 : integer :: iter ! CLUBB iteration
696 : integer :: num_subcols ! Number of subcolumns
697 : integer, parameter :: sequence_length = 1 ! Number of timesteps btn subcol calls
698 :
699 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: rho_ds_zt ! Dry static density (kg/m^3) on thermo levs
700 : real(r8), dimension(state%ngrdcol,pver) :: dz_g ! thickness of layer
701 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: delta_zm ! Difference in u wind altitudes
702 :
703 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1,hydromet_dim) :: hydromet ! Hydrometeor species
704 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1,hydromet_dim) :: wphydrometp ! Hydrometeor flux
705 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: Ncm ! Mean cloud droplet concentration, <N_c>
706 :
707 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: tke ! TKE
708 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: khzm ! Eddy diffusivity coef
709 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: Lscale_zm ! CLUBB's length scale on momentum (zm) levels
710 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: Lscale ! CLUBB's length scale
711 :
712 : logical, parameter :: &
713 : l_calc_weights_all_levs = .false. ! .false. if all time steps use the same
714 : ! weights at all vertical grid levels
715 : logical :: &
716 : l_calc_weights_all_levs_itime, & ! .true. if we calculate sample weights separately at all
717 : ! grid levels at the current time step
718 : l_rad_itime ! .true. if we calculate radiation at the current time step
719 :
720 : !---------------
721 : !Output from generate_silhs_sample
722 : !--------------
723 : real(r8), dimension(state%ngrdcol,subcol_SILHS_numsubcol,pverp-top_lev+1,pdf_dim) :: X_nl_all_levs ! Sample transformed to normal-lognormal
724 : real(r8), dimension(state%ngrdcol,subcol_SILHS_numsubcol,pverp-top_lev+1) :: lh_sample_point_weights ! Subcolumn weights
725 : integer, dimension(state%ngrdcol,subcol_SILHS_numsubcol,pverp-top_lev+1) :: X_mixt_comp_all_levs ! Which Mixture Component
726 :
727 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1, subcol_SILHS_numsubcol) :: &
728 : rc_all_points, & ! Calculate RCM from LH output
729 : rain_all_pts, & ! Calculate Rain from LH output
730 : nrain_all_pts, & ! Calculate Rain Conc from LH
731 : snow_all_pts, & ! Calculate Snow from LH output
732 : nsnow_all_pts, & ! Calculate Snow Conc from LH
733 : w_all_points, & ! Calculate W from LH output
734 : ice_all_pts, & ! Calculate Cld Ice from LH output
735 : nice_all_pts, & ! Calculate Num cld ice from LH
736 : nclw_all_pts ! Calculate Num cld wat from LH
737 :
738 : !----------------
739 : ! Output from clip_transform_silhs_output_api
740 : !----------------
741 : real( kind = core_rknd ), dimension(state%ngrdcol,subcol_SILHS_numsubcol,pverp-top_lev+1) :: &
742 : lh_rt_clipped, & ! rt generated from silhs sample points
743 : lh_thl_clipped, & ! thl generated from silhs sample points
744 : lh_rc_clipped, & ! rc generated from silhs sample points
745 : lh_rv_clipped, & ! rv generated from silhs sample points
746 : lh_Nc_clipped ! Nc generated from silhs sample points
747 :
748 : logical, parameter :: &
749 : l_use_Ncn_to_Nc = .true. ! Whether to call Ncn_to_Nc (.true.) or not (.false.);
750 : ! Ncn_to_Nc might cause problems with the MG microphysics
751 : ! since the changes made here (Nc-tendency) are not fed into
752 : ! the microphysics
753 :
754 :
755 : !----------------
756 : ! Output to history
757 : !----------------
758 : ! V. Larson note: These variables are on the zt (full) levels: why do they
759 : ! have dimension pverp? The pverp level corresponds to the CLUBB
760 : ! below-ground level.
761 : ! The variables in this paragraph are oriented like CAM variables (k=1 is
762 : ! the model top).
763 : ! They are flipped versions of CLUBB variables, for the entire chunk.
764 : real(r8), dimension(pcols*psubcols, pverp) :: RT_lh_out
765 : real(r8), dimension(pcols*psubcols, pverp) :: THL_lh_out
766 : real(r8), dimension(pcols*psubcols, pverp) :: OMEGA_lh_out
767 : real(r8), dimension(pcols*psubcols, pverp) :: WM_lh_out
768 : real(r8), dimension(pcols*psubcols, pverp) :: RVM_lh_out
769 : real(r8), dimension(pcols*psubcols, pverp) :: RCM_lh_out
770 : real(r8), dimension(pcols*psubcols, pverp) :: NCLW_lh_out
771 : real(r8), dimension(pcols*psubcols, pverp) :: ICE_lh_out
772 : real(r8), dimension(pcols*psubcols, pverp) :: NICE_lh_out
773 : real(r8), dimension(pcols*psubcols, pverp) :: RAIN_lh_out
774 : real(r8), dimension(pcols*psubcols, pverp) :: NRAIN_lh_out
775 : real(r8), dimension(pcols*psubcols, pverp) :: SNOW_lh_out
776 : real(r8), dimension(pcols*psubcols, pverp) :: NSNOW_lh_out
777 :
778 : real(r8), dimension(state_sc%psetcols) :: weights ! Subcol weights
779 :
780 : real(r8), dimension(pcols, pver) :: meansc_ice
781 : real(r8), dimension(pcols, pver) :: stdsc_ice
782 :
783 : real(r8), dimension(pcols, pver) :: meansc_liq
784 : real(r8), dimension(pcols, pver) :: stdsc_liq
785 :
786 : real(r8), dimension(pcols, pver) :: meansc_vap
787 : real(r8), dimension(pcols, pver) :: stdsc_vap
788 : real(r8), dimension(pcols, pver) :: grmn_eff_rad
789 : real(r8), dimension(pcols, pver) :: eff_cldfrac
790 : real(r8), dimension(pcols, pver) :: precip_frac_out
791 :
792 : real(r8) :: tmp_mean, diff_mean, rcubed
793 :
794 : !----------------
795 : ! Output from Est_Kessler_microphys
796 : !----------------
797 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: lh_Akm ! Monte Carlo estimate of Kessler Autoconversion
798 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKm ! Exact Kessler autoconversion
799 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKstd ! Exact Stdev of gba Kessler
800 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKstd_cld ! Exact w/in cloud stdev of gba Kessler
801 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKm_rcm ! Exact local gba Kessler auto based on rcm
802 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: AKm_rcc ! Exact local gba Kessler based on w/in cloud rc
803 : real(r8), dimension(state%ngrdcol,pverp-top_lev+1) :: lh_rcm_avg ! LH estimate of grid box avg liquid water
804 : real(r8), dimension(pcols,pverp) :: lh_AKm_out, AKm_out
805 :
806 : !----------------
807 : ! Needed to update State
808 : !----------------
809 : real(r8), dimension(pver) :: Temp_prof ! Subcolumn LWPT converted to Abs Temp
810 : real(r8), dimension(pver) :: SE_prof ! Static Energy calculated from Abs Temp
811 : real(r8), dimension(pver) :: No_cloud = 0.0_r8 ! Clear air condensate profile
812 : real(r8), dimension(pcols, pver) :: invs_exner ! inverse exner sent to conversion codw
813 : ! pcols for output to history
814 : real(r8) :: eff_rad_coef = 1.0_r8/(4.0_r8/3.0_r8*SHR_CONST_RHOFW*SHR_CONST_PI)
815 :
816 : !----------------
817 : ! Pointers
818 : !----------------
819 : real(r8), pointer, dimension(:) :: ztodt_ptr
820 : real(r8), pointer, dimension(:,:) :: thlm ! Mean temperature
821 : real(r8), pointer, dimension(:,:) :: ice_supersat_frac ! ice cloud fraction
822 : real(r8), pointer, dimension(:,:) :: rtm ! mean moisture mixing ratio
823 : real(r8), pointer, dimension(:,:) :: cld ! CAM cloud fraction
824 : real(r8), pointer, dimension(:,:) :: alst ! CLUBB liq cloud fraction
825 : real(r8), pointer, dimension(:,:) :: qrain ! micro_mg rain from previous step
826 : real(r8), pointer, dimension(:,:) :: qsnow
827 : real(r8), pointer, dimension(:,:) :: nrain ! micro_mg rain num conc
828 : real(r8), pointer, dimension(:,:) :: nsnow
829 :
830 : real(r8), pointer, dimension(:,:) :: tke_in ! TKE
831 : real(r8), pointer, dimension(:,:) :: khzm_in ! Eddy diffusivity coef
832 :
833 : logical, parameter :: l_est_kessler_microphys = .false.
834 : logical, parameter :: l_outfld_subcol = .false.
835 :
836 : type(grid) :: gr
837 :
838 : type(precipitation_fractions) :: precip_fracs
839 :
840 : ! Used as shortcuts to avoid typing hm_metadata%iiPDF_xx
841 : integer :: &
842 : iiPDF_chi, iiPDF_rr, iiPDF_w, iiPDF_Nr, &
843 : iiPDF_ri, iiPDF_Ni, iiPDF_Ncn, iiPDF_rs, iiPDF_Ns, &
844 : iirr, iiNr, iirs, iiri, &
845 : iirg, iiNs, iiNi, iiNg
846 :
847 : !------------------------------------------------
848 : ! Begin Code
849 : !------------------------------------------------
850 :
851 : #ifdef SILHS_OPENACC
852 : if ( l_est_kessler_microphys ) then
853 : call endrun('subcol_gen error: compilation with OpenACC requires l_est_kessler_microphys = .false.')
854 : end if
855 :
856 : if ( subcol_SILHS_constrainmn ) then
857 : call endrun('subcol_gen error: compilation with OpenACC requires subcol_SILHS_constrainmn = .false.')
858 : end if
859 :
860 : if ( subcol_SILHS_weight ) then
861 : call endrun('subcol_gen error: Importance sampling is not enabled for SILHS when using OpenACC. Set subcol_SILHS_weight to false.')
862 : end if
863 : #endif
864 :
865 : if (.not. allocated(state_sc%lat)) then
866 : call endrun('subcol_gen error: state_sc must be allocated before calling subcol_gen')
867 : end if
868 :
869 : if( rx_Nc ) then
870 : call endrun('subcol_gen_SILHS: rx_Nc not enabled')
871 : endif
872 :
873 : if (subcol_SILHS_meanice) then
874 : call endrun('subcol_gen_SILHS: subcol_SILHS_meanice = T not currently available')
875 : end if
876 :
877 : ! Determine num of columns and which chunk we're working on and what timestep
878 : ngrdcol = state%ngrdcol
879 : ncol = state%ncol
880 : lchnk = state%lchnk
881 : iter = get_nstep() ! #KTCtodo: The model iteration is passed into SILHS without taking
882 : ! substepping into account. I may need to change this in
883 : ! the future. Also, why does SILHS need an iter, but CLUBB
884 : ! does not?
885 : ! #ERDBG: The model iteration number is not used in SILHS unless
886 : ! sequence_length > 1, but nobody runs with that option.
887 :
888 : ! Copy hm_metadata indices to shortcuts
889 : iiPDF_chi = hm_metadata%iiPDF_chi
890 : iiPDF_Ncn = hm_metadata%iiPDF_Ncn
891 : iiPDF_rr = hm_metadata%iiPDF_rr
892 : iiPDF_w = hm_metadata%iiPDF_w
893 : iiPDF_Nr = hm_metadata%iiPDF_Nr
894 : iiPDF_ri = hm_metadata%iiPDF_ri
895 : iiPDF_Ni = hm_metadata%iiPDF_Ni
896 : iiPDF_rs = hm_metadata%iiPDF_rs
897 : iiPDF_Ns = hm_metadata%iiPDF_Ns
898 : iirr = hm_metadata%iirr
899 : iiNr = hm_metadata%iiNr
900 : iirs = hm_metadata%iirs
901 : iiri = hm_metadata%iiri
902 : iirg = hm_metadata%iirg
903 : iiNs = hm_metadata%iiNs
904 : iiNi = hm_metadata%iiNi
905 : iiNg = hm_metadata%iiNg
906 :
907 : !----------------
908 : ! Establish associations between pointers and physics buffer fields
909 : !----------------
910 : call pbuf_get_field(pbuf, thlm_idx, thlm)
911 : call pbuf_get_field(pbuf, ztodt_idx, ztodt_ptr)
912 : call pbuf_get_field(pbuf, ice_supersat_idx, ice_supersat_frac)
913 : call pbuf_get_field(pbuf, rtm_idx, rtm)
914 : call pbuf_get_field(pbuf, alst_idx, alst)
915 : call pbuf_get_field(pbuf, cld_idx, cld)
916 : call pbuf_get_field(pbuf, qrain_idx, qrain)
917 : call pbuf_get_field(pbuf, qsnow_idx, qsnow)
918 : call pbuf_get_field(pbuf, nrain_idx, nrain)
919 : call pbuf_get_field(pbuf, nsnow_idx, nsnow)
920 : call pbuf_get_field(pbuf, tke_idx, tke_in)
921 : call pbuf_get_field(pbuf, kvh_idx, khzm_in)
922 :
923 : ! Pull c_K from clubb parameters.
924 : c_K = clubb_params_single_col(ic_K)
925 :
926 : !----------------
927 : ! Copy state and populate numbers and values of sub-columns
928 : !----------------
929 : ztodt = ztodt_ptr(1)
930 : num_subcols = subcol_SILHS_numsubcol
931 :
932 : ! The number of vertical grid levels used in CLUBB is pverp, which is originally
933 : ! set in the call to setup_clubb_core_api from subroutine clubb_ini_cam. This
934 : ! is stored in CLUBB in the object gr%nz. This isn't changed in CLUBB.
935 : ! However, when SILHS is used, SILHS only uses pverp - top_lev + 1 vertical grid
936 : ! levels and also uses the gr%nz object. The value of gr%nz needs to be reset
937 : ! for SILHS here and then set again for CLUBB in subroutine clubb_tend_cam.
938 : gr%nz = pverp - top_lev + 1
939 :
940 : ! Calculate sample weights separately at all grid levels when
941 : ! radiation is not called
942 : l_calc_weights_all_levs_itime = .false. ! subcol_utils cannot compute weighted avgs
943 : ! when the weights vary with height.
944 : ! Don't set to true until this is fixed!!
945 :
946 :
947 : ! Setup the CLUBB vertical grid object. This must be done for each
948 : ! column as the z-distance between hybrid pressure levels can
949 : ! change easily.
950 : ! Define the CLUBB momentum grid (in height, units of m)
951 : do k = 1, pverp-top_lev+1
952 : do i = 1, ngrdcol
953 : zi_g(i,k) = state%zi(i,pverp-k+1)-state%zi(i,pverp)
954 : end do
955 : end do
956 :
957 :
958 : ! Define the CLUBB thermodynamic grid (in units of m)
959 : do k = 1, pver-top_lev+1
960 : do i = 1, ngrdcol
961 : zt_g(i,k+1) = state%zm(i,pver-k+1)-state%zi(i,pverp)
962 :
963 : ! Thermodynamic ghost point is below surface
964 : zt_g(i,1) = -1._r8*zt_g(i,2)
965 : end do
966 : end do
967 :
968 : do i=1, ncol
969 : ! Set the elevation of the surface
970 : sfc_elevation(i) = state%zi(i,pver+1)
971 : end do
972 :
973 : ! Heights need to be set at each timestep.
974 : call setup_grid_api( pverp+1-top_lev, ncol, sfc_elevation(1:ncol), l_implemented, & ! intent(in)
975 : grid_type, zi_g(1:ncol,2), zi_g(1:ncol,1), zi_g(1:ncol,pverp+1-top_lev), & ! intent(in)
976 : zi_g(1:ncol,:), zt_g(1:ncol,:), & ! intent(in)
977 : gr )
978 :
979 : ! Calculate the distance between grid levels on the host model grid,
980 : ! using host model grid indices.
981 : do k = top_lev, pver
982 : do i = 1, ngrdcol
983 : dz_g(i,k) = state%zi(i,k)-state%zi(i,k+1)
984 : end do
985 : end do
986 :
987 : ! Inverse delta_zm is required for the 3-level L-scale averaging
988 : do k = 1, pver-top_lev+1
989 : do i = 1, ngrdcol
990 : delta_zm(i,k+1) = state%zi(i,pverp-k)-state%zi(i,pverp-k+1)
991 :
992 : ! Handle CLUBB sub-sfc ghost point as done in clubb grid_class.F90
993 : delta_zm(i,1) = delta_zm(i,2)
994 : end do
995 : end do
996 :
997 : ! Compute dry static density on CLUBB vertical grid
998 : do k = 1, pver-top_lev+1
999 : do i = 1, ngrdcol
1000 : rho_ds_zt(i,k+1) = (rga)*state%pdel(i,pverp-k)/dz_g(i,pverp-k)
1001 :
1002 : ! CLUBB ghost point under the surface
1003 : rho_ds_zt(i,1) = rho_ds_zt(i,2)
1004 : end do
1005 : end do
1006 :
1007 : ! Set up hydromet array, flipped from CAM vert grid to CLUBB
1008 : if ( iirr > 0 ) then
1009 : ! If ixrain and family are greater than zero, then MG2 is
1010 : ! being used, and rain and snow are part of state. Otherwise,
1011 : ! diagnostic rain and snow from MG1 are used in hydromet.
1012 : if (ixrain > 0) then
1013 : do k = 1, pver-top_lev+1
1014 : do i = 1, ngrdcol
1015 : hydromet(i,k+1,iirr) = state%q(i,pver-k+1,ixrain)
1016 : end do
1017 : end do
1018 : else
1019 : do k = 1, pver-top_lev+1
1020 : do i = 1, ngrdcol
1021 : hydromet(i,k+1,iirr) = qrain(i,pver-k+1)
1022 : end do
1023 : end do
1024 : endif
1025 : endif
1026 :
1027 : if ( iiNr > 0 ) then
1028 : if (ixnumrain > 0) then
1029 : do k = 1, pver-top_lev+1
1030 : do i = 1, ngrdcol
1031 : hydromet(i,k+1,iiNr) = state%q(i,pver-k+1,ixnumrain)
1032 : end do
1033 : end do
1034 : else
1035 : do k = 1, pver-top_lev+1
1036 : do i = 1, ngrdcol
1037 : hydromet(i,k+1,iiNr) = nrain(i,pver-k+1)
1038 : end do
1039 : end do
1040 : endif
1041 : endif
1042 :
1043 : if ( iirs > 0 ) then
1044 : if (ixsnow > 0) then
1045 : do k = 1, pver-top_lev+1
1046 : do i = 1, ngrdcol
1047 : hydromet(i,k+1,iirs) = state%q(i,pver-k+1,ixsnow)
1048 : end do
1049 : end do
1050 : else
1051 : do k = 1, pver-top_lev+1
1052 : do i = 1, ngrdcol
1053 : hydromet(i,k+1,iirs) = qsnow(i,pver-k+1)
1054 : end do
1055 : end do
1056 : endif
1057 : endif
1058 :
1059 : if ( iiNs > 0 ) then
1060 : if (ixnumsnow > 0) then
1061 : do k = 1, pver-top_lev+1
1062 : do i = 1, ngrdcol
1063 : hydromet(i,k+1,iiNs) = state%q(i,pver-k+1,ixnumsnow)
1064 : end do
1065 : end do
1066 : else
1067 : do k = 1, pver-top_lev+1
1068 : do i = 1, ngrdcol
1069 : hydromet(i,k+1,iiNs) = nsnow(i,pver-k+1)
1070 : end do
1071 : end do
1072 : endif
1073 : endif
1074 :
1075 : if ( iiri > 0 ) then
1076 : do k = 1, pver-top_lev+1
1077 : do i = 1, ngrdcol
1078 : hydromet(i,k+1,iiri) = state%q(i,pver-k+1,ixcldice)
1079 : end do
1080 : end do
1081 : endif
1082 :
1083 : if ( iiNi > 0 ) then
1084 : do k = 1, pver-top_lev+1
1085 : do i = 1, ngrdcol
1086 : hydromet(i,k+1,iiNi) = state%q(i,pver-k+1,ixnumice)
1087 : end do
1088 : end do
1089 : endif
1090 :
1091 : do k = 1, hydromet_dim ! ghost point below the surface
1092 : do i = 1, ngrdcol
1093 : hydromet(i,1,k) = hydromet(i,2,k)
1094 : end do
1095 : end do
1096 :
1097 : do k = 1, pver-top_lev+1
1098 : do i = 1, ngrdcol
1099 : Ncm(i,k+1) = state%q(i,pver-k+1,ixnumliq)
1100 : Ncm(i,1) = Ncm(i,2)
1101 : end do
1102 : end do
1103 :
1104 : ! Convert from CAM vertical grid to CLUBB
1105 : do k = 1, pverp-top_lev+1
1106 : do i = 1, ngrdcol
1107 : ice_supersat_frac_in(i,k) = ice_supersat_frac(i,pverp-k+1)
1108 : end do
1109 : end do
1110 :
1111 :
1112 : do k = 1, pver-top_lev+1
1113 : do i = 1, ngrdcol
1114 : cld_frac_in(i,k+1) = alst(i,pver-k+1)
1115 : cld_frac_in(i,1) = cld_frac_in(i,2) ! Ghost pt below surface
1116 : end do
1117 : end do
1118 :
1119 : ! Calculate a clubb-specific exner function
1120 : ! (This is grid mean, as pressure levels do not change in
1121 : ! the subcolumn state)
1122 : do k = 1, pver-top_lev+1
1123 : do i = 1, ngrdcol
1124 : invs_exner(i,k) = ((state%pmid(i,k)/p0_clubb)**(cappa))
1125 : end do
1126 : end do
1127 :
1128 : ! Call setup_pdf_parameters to get the CLUBB PDF ready for SILHS
1129 : ! Compute Num concentration of cloud nuclei
1130 : do k = 1, pverp-top_lev+1
1131 : do i = 1, ngrdcol
1132 : Nc_in_cloud(i,k) = Ncm(i,k) / max( cld_frac_in(i,k), cloud_frac_min )
1133 : end do
1134 : end do
1135 :
1136 : ! The variable wphydrometp is only used when l_calc_w_corr is enabled.
1137 : ! The l_calc_w_corr flag is turned off by default, so wphydrometp will
1138 : ! simply be set to 0 to simplify matters.
1139 : wphydrometp = 0.0_r8
1140 :
1141 : do k = 1, pverp-top_lev+1
1142 : do i = 1, ngrdcol
1143 : khzm(i,k) = khzm_in(i,pverp-k+1)
1144 : end do
1145 : end do
1146 :
1147 : ! Allocate 2D arrays in precip_fracs for all grid columns and vertical levels
1148 : call init_precip_fracs_api( pverp-top_lev+1, ngrdcol, &
1149 : precip_fracs )
1150 :
1151 : call setup_pdf_parameters_api( gr, pverp-top_lev+1, ngrdcol, pdf_dim, hydromet_dim, ztodt, & ! In
1152 : Nc_in_cloud, cld_frac_in, khzm, & ! In
1153 : ice_supersat_frac_in, hydromet, wphydrometp, & ! In
1154 : corr_array_n_cloud, corr_array_n_below, & ! In
1155 : hm_metadata, & ! In
1156 : pdf_params_chnk(lchnk), & ! In
1157 : clubb_params_single_col, & ! In
1158 : clubb_config_flags%iiPDF_type, & ! In
1159 : clubb_config_flags%l_use_precip_frac, & ! In
1160 : clubb_config_flags%l_predict_upwp_vpwp, & ! In
1161 : clubb_config_flags%l_diagnose_correlations, & ! In
1162 : clubb_config_flags%l_calc_w_corr, & ! In
1163 : clubb_config_flags%l_const_Nc_in_cloud, & ! In
1164 : clubb_config_flags%l_fix_w_chi_eta_correlations, & ! In
1165 : stats_metadata, & ! In
1166 : stats_zt, stats_zm, stats_sfc, & ! In
1167 : hydrometp2, & ! Inout
1168 : mu_x_1, mu_x_2, & ! Out
1169 : sigma_x_1, sigma_x_2, & ! Out
1170 : corr_array_1, corr_array_2, & ! Out
1171 : corr_cholesky_mtx_1, corr_cholesky_mtx_2, & ! Out
1172 : precip_fracs ) ! Inout
1173 :
1174 : ! In order for Lscale to be used properly, it needs to be passed out of
1175 : ! advance_clubb_core, saved to the pbuf, and then pulled out of the
1176 : ! pbuf for use here. The profile of Lscale is passed into subroutine
1177 : ! generate_silhs_sample_api for use in calculating the vertical
1178 : ! correlation coefficient. Rather than output Lscale directly, its
1179 : ! value can be calculated from other fields that are already output to
1180 : ! pbuf. The equation relating Lscale to eddy diffusivity is:
1181 : !
1182 : ! Kh = c_K * Lscale * sqrt( TKE ).
1183 : !
1184 : ! Both Kh and TKE are written to the pbuf, and c_K is easily extracted
1185 : ! from CLUBB's tunable parameters. The equation for Lscale is:
1186 : !
1187 : ! Lscale = Kh / ( c_K * sqrt( TKE ) ).
1188 : !
1189 : ! Since Kh and TKE are output on momentum (interface) grid levels, the
1190 : ! resulting calculation of Lscale is also found on momentum levels. It
1191 : ! needs to be interpolated back to thermodynamic (midpoint) grid levels
1192 : ! for further use.
1193 : do k = 1, pverp-top_lev+1
1194 : do i = 1, ngrdcol
1195 : tke(i,k) = tke_in(i,pverp-k+1)
1196 : end do
1197 : end do
1198 :
1199 : do k = 1, pverp-top_lev+1
1200 : do i = 1, ngrdcol
1201 : Lscale_zm(i,k) = khzm(i,k) / ( c_K * sqrt( max( tke(i,k), em_min ) ) )
1202 : end do
1203 : end do
1204 :
1205 : do i = 1, ngrdcol
1206 : Lscale(i,1) = Lscale_zm(i,1) + ( Lscale_zm(i,2) - Lscale_zm(i,1) ) &
1207 : * ( zt_g(i,1) - zi_g(i,1) ) / ( zi_g(i,2) - zi_g(i,1) )
1208 : end do
1209 :
1210 : do k = 2, pverp-top_lev+1
1211 : do i = 1, ngrdcol
1212 : Lscale(i,k) = Lscale_zm(i,k-1) + ( Lscale_zm(i,k) - Lscale_zm(i,k-1) ) &
1213 : * ( zt_g(i,k) - zi_g(i,k-1) ) / ( zi_g(i,k) - zi_g(i,k-1) )
1214 : end do
1215 : end do
1216 :
1217 : do k = 2, pverp-top_lev+1
1218 : do i = 1, ngrdcol
1219 : Lscale(i,:) = max( Lscale(i,:), 0.01_r8 )
1220 : end do
1221 : end do
1222 :
1223 : !$acc data create( X_mixt_comp_all_levs, X_nl_all_levs, lh_rc_clipped, lh_Nc_clipped, &
1224 : !$acc& lh_sample_point_weights, lh_rt_clipped, lh_rt_clipped, &
1225 : !$acc& lh_rv_clipped, lh_thl_clipped, THL_lh_out, &
1226 : !$acc& RT_lh_out, RCM_lh_out, NCLW_lh_out, ICE_lh_out, &
1227 : !$acc& NICE_lh_out, RVM_lh_out, THL_lh_out, RAIN_lh_out, &
1228 : !$acc& NRAIN_lh_out, SNOW_lh_out, NSNOW_lh_out, WM_lh_out, &
1229 : !$acc& OMEGA_lh_out ) &
1230 : !$acc& copyin( state, state%zm, state%phis, rho_ds_zt, invs_exner ) &
1231 : !$acc& copyout( state%t, state%s, state%omega, state_sc%q )
1232 : !$acc& async(1)
1233 :
1234 : ! Set the seed to the random number generator based on a quantity that
1235 : ! will be reproducible for restarts.
1236 : lh_seed = int( 1.0e4_r8 * rtm(1,pver), kind = genrand_intg )
1237 :
1238 : ! Let's generate some subcolumns!!!!!
1239 : call generate_silhs_sample_api( &
1240 : iter, pdf_dim, num_subcols, sequence_length, pverp-top_lev+1, ngrdcol, & ! In
1241 : l_calc_weights_all_levs_itime, & ! In
1242 : pdf_params_chnk(lchnk), delta_zm, Lscale, & ! In
1243 : lh_seed, hm_metadata, & ! In
1244 : rho_ds_zt, & ! In
1245 : mu_x_1, mu_x_2, sigma_x_1, sigma_x_2, & ! In
1246 : corr_cholesky_mtx_1, corr_cholesky_mtx_2, & ! In
1247 : precip_fracs, silhs_config_flags, & ! In
1248 : vert_decorr_coef, & ! In
1249 : stats_metadata, & ! In
1250 : stats_lh_zt, stats_lh_sfc, & ! InOut
1251 : X_nl_all_levs, X_mixt_comp_all_levs, & ! Out
1252 : lh_sample_point_weights) ! Out
1253 :
1254 : ! Extract clipped variables from subcolumns
1255 : call clip_transform_silhs_output_api( gr, pverp-top_lev+1, ngrdcol, num_subcols, & ! In
1256 : pdf_dim, hydromet_dim, hm_metadata, & ! In
1257 : X_mixt_comp_all_levs, & ! In
1258 : X_nl_all_levs, & ! In
1259 : pdf_params_chnk(lchnk), & ! In
1260 : l_use_Ncn_to_Nc, & ! In
1261 : lh_rt_clipped, lh_thl_clipped, & ! Out
1262 : lh_rc_clipped, lh_rv_clipped, & ! Out
1263 : lh_Nc_clipped ) ! Out
1264 : !$acc wait
1265 :
1266 : if ( l_est_kessler_microphys ) then
1267 : call endrun('subcol_SILHS: l_est_kessler_microphys = T is not currently supported')
1268 : end if
1269 :
1270 : !-------------------------------------------------------------------------
1271 : ! Convert from CLUBB vertical grid to CAM grid
1272 : !------------------------------------------------------------------------
1273 : ! This kernel is executed in stream 1:
1274 : !$acc parallel loop collapse(3) default(present) async(1)
1275 : do k = top_lev, pverp
1276 : do j = 1, num_subcols
1277 : do i = 1, ngrdcol
1278 : RT_lh_out( num_subcols*(i-1)+j,k ) = lh_rt_clipped(i,j,pverp-k+1)
1279 : RCM_lh_out( num_subcols*(i-1)+j,k ) = lh_rc_clipped(i,j,pverp-k+1)
1280 : NCLW_lh_out( num_subcols*(i-1)+j,k ) = lh_Nc_clipped(i,j,pverp-k+1)
1281 : RVM_lh_out( num_subcols*(i-1)+j,k ) = lh_rv_clipped(i,j,pverp-k+1)
1282 : THL_lh_out( num_subcols*(i-1)+j,k ) = lh_thl_clipped(i,j,pverp-k+1)
1283 : end do
1284 : end do
1285 : end do
1286 :
1287 : ! This kernel is executed in stream 2:
1288 : !$acc parallel loop collapse(3) default(present) async(2)
1289 : do k = top_lev, pverp
1290 : do j = 1, num_subcols
1291 : do i = 1, ngrdcol
1292 : ICE_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_ri)
1293 : NICE_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_Ni)
1294 : RAIN_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_rr)
1295 : NRAIN_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_Nr)
1296 : SNOW_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_rs)
1297 : NSNOW_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_Ns)
1298 : WM_lh_out( num_subcols*(i-1)+j,k ) = X_nl_all_levs(i,j,pverp-k+1,iiPDF_w)
1299 : end do
1300 : end do
1301 : end do
1302 :
1303 : ! This kernel is executed in stream 2 because WM_lh_out comes from stream 2:
1304 : !$acc parallel loop collapse(3) default(present) async(2)
1305 : do k = top_lev, pverp
1306 : do j = 1, num_subcols
1307 : do i = 1, ngrdcol
1308 : OMEGA_lh_out( num_subcols*(i-1)+j,k ) = -1._r8 * WM_lh_out(num_subcols*(i-1)+j,k) &
1309 : * rho_ds_zt(i,pverp-k+1) * gravit
1310 : end do
1311 : end do
1312 : end do
1313 :
1314 : if ( l_est_kessler_microphys ) then
1315 : do k = top_lev, pverp
1316 : do j = 1, num_subcols
1317 : do i = 1, ngrdcol
1318 : AKm_out(i,k) = AKm(i,pverp-k+1)
1319 : lh_AKm_out(i,k) = lh_AKm(i,pverp-k+1)
1320 : end do
1321 : end do
1322 : end do
1323 : end if
1324 :
1325 : ! Pack up weights
1326 : ! Using grid level 2 always won't work if weights vary with height.
1327 : call subcol_pack(lchnk, lh_sample_point_weights(:,:,2), weights )
1328 : call subcol_set_weight(lchnk, weights)
1329 :
1330 : ! Constrain the sample distribution of cloud water and ice to the same mean
1331 : ! as the grid to prevent negative condensate errors
1332 : if(subcol_SILHS_constrainmn) then
1333 :
1334 : do i = 1, ngrdcol
1335 :
1336 : stncol = num_subcols*(i-1)
1337 :
1338 : call subcol_constrainmn( num_subcols, ICE_lh_out(stncol+1:stncol+num_subcols,:), &
1339 : weights(stncol+1:stncol+num_subcols), &
1340 : state%q(i,:,ixcldice), meansc_ice(i,:), stdsc_ice(i,:) )
1341 : if ( ixrain > 0 ) &
1342 : call subcol_constrainmn( num_subcols, RAIN_lh_out(stncol+1:stncol+num_subcols,:), &
1343 : weights(stncol+1:stncol+num_subcols), &
1344 : state%q(i,:,ixrain) )
1345 : if ( ixsnow > 0 ) &
1346 : call subcol_constrainmn( num_subcols, SNOW_lh_out(stncol+1:stncol+num_subcols,:), &
1347 : weights(stncol+1:stncol+num_subcols), &
1348 : state%q(i,:,ixsnow) )
1349 : call subcol_constrainmn( num_subcols, RCM_lh_out(stncol+1:stncol+num_subcols,:), &
1350 : weights(stncol+1:stncol+num_subcols), &
1351 : state%q(i,:,ixcldliq), meansc_liq(i,:), stdsc_liq(i,:) )
1352 : call subcol_constrainmn( num_subcols, RVM_lh_out(stncol+1:stncol+num_subcols,:), &
1353 : weights(stncol+1:stncol+num_subcols), &
1354 : state%q(i,:,ixq), meansc_vap(i,:), stdsc_vap(i,:) )
1355 : call subcol_constrainmn( num_subcols, NICE_lh_out(stncol+1:stncol+num_subcols,:), &
1356 : weights(stncol+1:stncol+num_subcols), &
1357 : state%q(i,:,ixnumice) )
1358 : if ( ixnumrain > 0 ) &
1359 : call subcol_constrainmn( num_subcols, NRAIN_lh_out(stncol+1:stncol+num_subcols,:), &
1360 : weights(stncol+1:stncol+num_subcols), &
1361 : state%q(i,:,ixnumrain) )
1362 : if ( ixnumsnow > 0 ) &
1363 : call subcol_constrainmn( num_subcols, NSNOW_lh_out(stncol+1:stncol+num_subcols,:), &
1364 : weights(stncol+1:stncol+num_subcols), &
1365 : state%q(i,:,ixnumsnow) )
1366 : call subcol_constrainmn( num_subcols, NCLW_lh_out(stncol+1:stncol+num_subcols,:), &
1367 : weights(stncol+1:stncol+num_subcols), &
1368 : state%q(i,:,ixnumliq) )
1369 : do k = top_lev, pver
1370 : ! Look for exceptionally large values of condensate
1371 : if(ANY(ICE_lh_out(stncol+1:stncol+num_subcols,k) .gt. 0.01_r8)) then
1372 : ! Clip the large values
1373 : where(ICE_lh_out(stncol+1:stncol+num_subcols,k) .gt. 0.01_r8)
1374 : ICE_lh_out(stncol+1:stncol+num_subcols,k) = 0.01_r8
1375 : NICE_lh_out(stncol+1:stncol+num_subcols,k) = 1.5e+7_r8
1376 : end where
1377 : ! Recalculate the weighted subcolumn mean
1378 : tmp_mean = meansc( ICE_lh_out( stncol+1:stncol+num_subcols, k ), &
1379 : weights(stncol+1:stncol+num_subcols), &
1380 : real(num_subcols,r8) )
1381 : ! Calculate the difference between the weighted mean and grid mean
1382 : diff_mean = state%q(i,k,ixcldice)-tmp_mean
1383 : ! Add the difference to each subcolumn
1384 : ICE_lh_out(stncol+1:stncol+num_subcols,k) = &
1385 : ICE_lh_out(stncol+1:stncol+num_subcols,k)+diff_mean
1386 : ! Recalculate the weight subcolumn mean for ice num conc
1387 : tmp_mean = meansc( NICE_lh_out( stncol+1:stncol+num_subcols, k ), &
1388 : weights(stncol+1:stncol+num_subcols), &
1389 : real(num_subcols,r8) )
1390 : ! Calculate the difference between the weighted mean and grid mean
1391 : diff_mean = state%q(i,k,ixnumice)-tmp_mean
1392 : ! Add the difference to each subcolumn
1393 : if(diff_mean.gt.0.0_r8) then
1394 : NICE_lh_out(stncol+1:stncol+num_subcols,k) = &
1395 : NICE_lh_out(stncol+1:stncol+num_subcols,k)+diff_mean
1396 : else ! just use the grid mean in each subcolumn
1397 : NICE_lh_out(stncol+1:stncol+num_subcols,k) = &
1398 : state%q(i,k,ixnumice)
1399 : end if
1400 : ! Test adjusted means for debugging
1401 : tmp_mean = meansc( ICE_lh_out( stncol+1:stncol+num_subcols, k ), &
1402 : weights(stncol+1:stncol+num_subcols), &
1403 : real(num_subcols,r8) )
1404 : diff_mean = state%q(i,k,ixcldice)-tmp_mean
1405 : tmp_mean = meansc( NICE_lh_out( stncol+1:stncol+num_subcols, k ), &
1406 : weights(stncol+1:stncol+num_subcols), &
1407 : real(num_subcols,r8) )
1408 : diff_mean = state%q(i,k,ixnumice)-tmp_mean
1409 : endif
1410 : end do ! k = top_lev, pver
1411 : end do
1412 : endif ! subcol_silhs_constrainm
1413 :
1414 :
1415 : !---------------------------------------------------
1416 : ! Updating state variables
1417 : !---------------------------------------------------
1418 : ! Code to update the state variables for interactive runs
1419 : ! This kernel is executed in stream 3, but waits for stream 1
1420 : ! because THL_lh_out and RCM_lh_out come from stream 1:
1421 : !$acc parallel loop collapse(3) default(present) wait(1) async(3)
1422 : do k = 1, pver-top_lev+1
1423 : do j = 1, num_subcols
1424 : do i = 1, ngrdcol
1425 :
1426 : state_sc%t(num_subcols*(i-1)+j,k) = THL_lh_out(num_subcols*(i-1)+j,k) * invs_exner(i,k) &
1427 : + Lv * RCM_lh_out(num_subcols*(i-1)+j,k) / Cp
1428 :
1429 : state_sc%s(num_subcols*(i-1)+j,k) = cpair * state_sc%t(num_subcols*(i-1)+j,k) &
1430 : + gravit * state%zm(i,k) + state%phis(i)
1431 : end do
1432 : end do
1433 : end do
1434 :
1435 : ! This kernel is executed in stream 4, but waits for stream 1 and 2
1436 : ! because RVM_lh_out is from stream 1 and OMEGA_lh_out is from stream 2:
1437 : !$acc parallel loop collapse(3) default(present) wait(1,2) async(4)
1438 : do k = 1, pver-top_lev+1
1439 : do j = 1, num_subcols
1440 : do i = 1, ngrdcol
1441 : ! Vertical Velocity is not part of the energy conservation checks, but
1442 : ! we need to be careful here, because the SILHS output VV is noisy.
1443 : state_sc%omega(num_subcols*(i-1)+j,k) = OMEGA_lh_out(num_subcols*(i-1)+j,k)
1444 : state_sc%q(num_subcols*(i-1)+j,k,ixq) = RVM_lh_out(num_subcols*(i-1)+j,k)
1445 : end do
1446 : end do
1447 : end do
1448 :
1449 :
1450 : if (subcol_SILHS_q_to_micro) then ! Send SILHS predicted constituents to microp
1451 :
1452 : ! This kernel is executed in stream 5, but waits for stream 1 and 2
1453 : ! because RCM_lh_out is from stream 1 and ICE_lh_out is from stream 2:
1454 : !$acc parallel loop collapse(3) default(present) wait(1,2) async(5)
1455 : do k = 1, pver-top_lev+1
1456 : do j = 1, num_subcols
1457 : do i = 1, ngrdcol
1458 : state_sc%q(num_subcols*(i-1)+j,k,ixcldliq) = RCM_lh_out(num_subcols*(i-1)+j,k)
1459 : state_sc%q(num_subcols*(i-1)+j,k,ixcldice) = ICE_lh_out(num_subcols*(i-1)+j,k)
1460 : end do
1461 : end do
1462 : end do
1463 :
1464 : if (ixrain > 0) then
1465 : ! This kernel is executed in stream 6, but waits for stream 2
1466 : ! because RAIN_lh_out is from stream 2:
1467 : !$acc parallel loop collapse(3) default(present) wait(2) async(6)
1468 : do k = 1, pver-top_lev+1
1469 : do j = 1, num_subcols
1470 : do i = 1, ngrdcol
1471 : state_sc%q(num_subcols*(i-1)+j,k,ixrain) = RAIN_lh_out(num_subcols*(i-1)+j,k)
1472 : end do
1473 : end do
1474 : end do
1475 : end if
1476 :
1477 : if (ixsnow > 0) then
1478 : ! This kernel is executed in stream 7, but waits for stream 2
1479 : ! because SNOW_lh_out is from stream 2:
1480 : !$acc parallel loop collapse(3) default(present) wait(2) async(7)
1481 : do k = 1, pver-top_lev+1
1482 : do j = 1, num_subcols
1483 : do i = 1, ngrdcol
1484 : state_sc%q(num_subcols*(i-1)+j,k,ixsnow) = SNOW_lh_out(num_subcols*(i-1)+j,k)
1485 : end do
1486 : end do
1487 : end do
1488 : end if
1489 :
1490 : else
1491 :
1492 : do k = 1, pver-top_lev+1
1493 : do j = 1, num_subcols
1494 : do i = 1, ngrdcol
1495 : state_sc%q(num_subcols*(i-1)+j,k,ixcldliq) = state%q(i,k,ixcldliq)
1496 : state_sc%q(num_subcols*(i-1)+j,k,ixcldice) = state%q(i,k,ixcldice)
1497 : if (ixrain > 0) then
1498 : state_sc%q(num_subcols*(i-1)+j,k,ixrain) = state%q(i,k,ixrain)
1499 : end if
1500 : if (ixsnow > 0) then
1501 : state_sc%q(num_subcols*(i-1)+j,k,ixsnow) = state%q(i,k,ixsnow)
1502 : end if
1503 : end do
1504 : end do
1505 : end do
1506 :
1507 : endif
1508 :
1509 : if (subcol_SILHS_n_to_micro) then ! Send SILHS predicted number conc to microp
1510 :
1511 : ! This kernel is executed in stream 8, but waits for stream 1 and 2
1512 : ! because NCLW_lh_out is from stream 1 and NICE_lh_out is from stream 2:
1513 : !$acc parallel loop collapse(3) default(present) wait(1,2) async(8)
1514 : do k = 1, pver-top_lev+1
1515 : do j = 1, num_subcols
1516 : do i = 1, ngrdcol
1517 : state_sc%q(num_subcols*(i-1)+j,k,ixnumice) = NICE_lh_out(num_subcols*(i-1)+j,k)
1518 : state_sc%q(num_subcols*(i-1)+j,k,ixnumliq) = NCLW_lh_out(num_subcols*(i-1)+j,k)
1519 : end do
1520 : end do
1521 : end do
1522 :
1523 : if (ixnumrain > 0) then
1524 : ! This kernel is executed in stream 9, but waits for stream 2
1525 : ! because NRAIN_lh_out is from stream 2:
1526 : !$acc parallel loop collapse(3) default(present) wait(2) async(9)
1527 : do k = 1, pver-top_lev+1
1528 : do j = 1, num_subcols
1529 : do i = 1, ngrdcol
1530 : state_sc%q(num_subcols*(i-1)+j,k,ixnumrain) = NRAIN_lh_out(num_subcols*(i-1)+j,k)
1531 : end do
1532 : end do
1533 : end do
1534 : end if
1535 :
1536 : if (ixnumsnow > 0) then
1537 : ! This kernel is executed in stream 10, but waits for stream 2
1538 : ! because NSNOW_lh_out is from stream 2:
1539 : !$acc parallel loop collapse(3) default(present) wait(2) async(10)
1540 : do k = 1, pver-top_lev+1
1541 : do j = 1, num_subcols
1542 : do i = 1, ngrdcol
1543 : state_sc%q(num_subcols*(i-1)+j,k,ixnumsnow) = NSNOW_lh_out(num_subcols*(i-1)+j,k)
1544 : end do
1545 : end do
1546 : end do
1547 : end if
1548 :
1549 : else
1550 :
1551 : do k = 1, pver-top_lev+1
1552 : do j = 1, num_subcols
1553 : do i = 1, ngrdcol
1554 : state_sc%q(num_subcols*(i-1)+j,k,ixnumliq) = state%q(i,k,ixnumliq)
1555 : state_sc%q(num_subcols*(i-1)+j,k,ixnumice) = state%q(i,k,ixnumice)
1556 : if (ixnumrain > 0) then
1557 : state_sc%q(num_subcols*(i-1)+j,k,ixnumrain) = state%q(i,k,ixnumrain)
1558 : end if
1559 : if (ixnumsnow > 0) then
1560 : state_sc%q(num_subcols*(i-1)+j,k,ixnumsnow) = state%q(i,k,ixnumsnow)
1561 : end if
1562 : end do
1563 : end do
1564 : end do
1565 :
1566 : endif
1567 :
1568 : ! This kernel is executed in stream 8, because state_sc%q(:,:,ixnumliq) and
1569 : ! state_sc%q(:,:,ixnumice) are from stream 8
1570 : !$acc parallel loop collapse(3) default(present) async(8)
1571 : do k = 1, pver-top_lev+1
1572 : do j = 1, num_subcols
1573 : do i = 1, ngrdcol
1574 : ! Change liq and ice (and rain and snow) num conc zeros to min values (1e-12)
1575 : if (state_sc%q(num_subcols*(i-1)+j,k,ixnumliq) .lt. min_num_conc) then
1576 : state_sc%q(num_subcols*(i-1)+j,k,ixnumliq) = min_num_conc
1577 : end if
1578 :
1579 : if (state_sc%q(num_subcols*(i-1)+j,k,ixnumice) .lt. min_num_conc) then
1580 : state_sc%q(num_subcols*(i-1)+j,k,ixnumice) = min_num_conc
1581 : end if
1582 : end do
1583 : end do
1584 : end do
1585 :
1586 : if (ixnumrain > 0) then
1587 : ! This kernel is executed in stream 9, because state_sc%q(:,:,ixnumrain) is
1588 : ! from stream 9
1589 : !$acc parallel loop collapse(3) default(present) async(9)
1590 : do k = 1, pver-top_lev+1
1591 : do j = 1, num_subcols
1592 : do i = 1, ngrdcol
1593 : if(state_sc%q(num_subcols*(i-1)+j,k,ixnumrain) .lt. min_num_conc) then
1594 : state_sc%q(num_subcols*(i-1)+j,k,ixnumrain) = min_num_conc
1595 : end if
1596 : end do
1597 : end do
1598 : end do
1599 : endif
1600 :
1601 : if (ixnumsnow > 0) then
1602 : ! This kernel is executed in stream 10, because state_sc%q(:,:,ixnumsnow) is
1603 : ! from stream 10
1604 : !$acc parallel loop collapse(3) default(present) async(10)
1605 : do k = 1, pver-top_lev+1
1606 : do j = 1, num_subcols
1607 : do i = 1, ngrdcol
1608 : if(state_sc%q(num_subcols*(i-1)+j,k,ixnumsnow) .lt. min_num_conc) then
1609 : state_sc%q(num_subcols*(i-1)+j,k,ixnumsnow) = min_num_conc
1610 : end if
1611 : end do
1612 : end do
1613 : end do
1614 : endif
1615 :
1616 : if ( l_outfld_subcol ) then
1617 :
1618 : do k = 1, pver-top_lev+1
1619 : do i = 1, ngrdcol
1620 : do j = 1, num_subcols
1621 :
1622 : ! Calc effective cloud fraction for testing
1623 : if ( ( lh_rc_clipped(i,j,pverp-k+1) .gt. qsmall ) &
1624 : .or. ( X_nl_all_levs(i,j,pverp-k+1,iiPDF_ri) .gt. qsmall ) ) then
1625 : eff_cldfrac(i,k) = eff_cldfrac(i,k) + lh_sample_point_weights(i,j,pverp-k+1)
1626 : else
1627 : eff_cldfrac(i,k) = 0.0_r8
1628 : endif
1629 :
1630 : end do
1631 :
1632 : eff_cldfrac(i,k) = eff_cldfrac(i,k)/real(num_subcols, kind=r8)
1633 :
1634 : end do
1635 : end do
1636 :
1637 : ! Pack precip_frac for output
1638 : do k = 2, pverp-top_lev+1
1639 : do i = 1, ngrdcol
1640 : precip_frac_out(i,pver-k+2) = precip_fracs%precip_frac(i,k)
1641 : end do
1642 : end do
1643 :
1644 : call outfld( 'SILHS_THLM_SCOL', THL_lh_out, pcols*psubcols, lchnk )
1645 : call outfld( 'SILHS_RT_SCOL', RT_lh_out, pcols*psubcols, lchnk )
1646 : call outfld( 'SILHS_OMEGA_SCOL', OMEGA_lh_out, pcols*psubcols, lchnk )
1647 : call outfld( 'SILHS_WM_SCOL', WM_lh_out, pcols*psubcols, lchnk )
1648 : call outfld( 'SILHS_RCM_SCOL', RCM_lh_out, pcols*psubcols, lchnk )
1649 : call outfld( 'SILHS_RICLD_SCOL', ICE_lh_out, pcols*psubcols, lchnk )
1650 : call outfld( 'SILHS_NICLD_SCOL', NICE_lh_out, pcols*psubcols, lchnk )
1651 : call outfld( 'SILHS_NCLD_SCOL', NCLW_lh_out, pcols*psubcols, lchnk )
1652 : call outfld( 'SILHS_RRAIN_SCOL', RAIN_lh_out, pcols*psubcols, lchnk )
1653 : call outfld( 'SILHS_NRAIN_SCOL', NRAIN_lh_out, pcols*psubcols, lchnk )
1654 : call outfld( 'SILHS_WEIGHT_SCOL', weights, pcols*psubcols, lchnk )
1655 : call outfld( 'NR_IN_LH', nrain, pcols, lchnk )
1656 : call outfld( 'SILHS_RTM', rtm, pcols, lchnk )
1657 : call outfld( 'SILHS_THLM', thlm, pcols, lchnk )
1658 : call outfld( 'SILHS_QC_IN', state%q(:,:,ixcldliq), pcols, lchnk )
1659 : call outfld( 'SILHS_QI_IN', state%q(:,:,ixcldice), pcols, lchnk )
1660 : call outfld( 'SILHS_NC_IN', state%q(:,:,ixnumliq), pcols, lchnk )
1661 : call outfld( 'SILHS_NI_IN', state%q(:,:,ixnumice), pcols, lchnk )
1662 : if ( l_est_kessler_microphys ) then
1663 : call outfld( 'AKM_CLUBB', AKm_out, pcols, lchnk )
1664 : call outfld( 'AKM_LH_CLUBB', lh_AKm_out, pcols, lchnk )
1665 : end if
1666 : call outfld( 'INVS_EXNER', invs_exner, pcols, lchnk )
1667 : call outfld( 'SILHS_ZTODT', ztodt_ptr, pcols, lchnk )
1668 : if ( subcol_SILHS_constrainmn ) then
1669 : call outfld( 'SILHS_MSC_CLDICE', meansc_ice, pcols, lchnk )
1670 : call outfld( 'SILHS_STDSC_CLDICE', stdsc_ice, pcols, lchnk )
1671 : if ( ixsnow > 0 ) then
1672 : call outfld( 'SILHS_MSC_CLDLIQ', meansc_liq, pcols, lchnk )
1673 : call outfld( 'SILHS_STDSC_CLDLIQ', stdsc_liq, pcols, lchnk )
1674 : call outfld( 'SILHS_MSC_Q', meansc_vap, pcols, lchnk )
1675 : call outfld( 'SILHS_STDSC_Q', stdsc_vap, pcols, lchnk )
1676 : endif ! ixsnow > 0
1677 : endif ! subcol_SILHS_constrainmn
1678 : call outfld( 'SILHS_EFF_CLDFRAC', eff_cldfrac, pcols, lchnk )
1679 : call outfld( 'SILHS_CLUBB_PRECIP_FRAC', precip_frac_out, pcols, lchnk )
1680 : call outfld( 'SILHS_CLUBB_ICE_SS_FRAC', ice_supersat_frac, pcols, lchnk )
1681 : end if
1682 :
1683 : !$acc end data
1684 : !$acc wait
1685 :
1686 : #endif
1687 : #endif
1688 0 : end subroutine subcol_gen_SILHS
1689 :
1690 0 : subroutine subcol_ptend_avg_SILHS(ptend_sc, ngrdcol, lchnk, ptend)
1691 0 : use physics_buffer, only: physics_buffer_desc
1692 : use subcol_utils, only: subcol_ptend_get_firstsubcol, subcol_ptend_avg_shr, &
1693 : subcol_get_weight, subcol_get_filter, &
1694 : is_filter_set, is_weight_set
1695 :
1696 : !-----------------------------------
1697 : ! Average the subcolumns dimension (pcols*psubcols) to the grid dimension (pcols)
1698 : !-----------------------------------
1699 :
1700 : type(physics_ptend), intent(in) :: ptend_sc ! intent in
1701 : integer, intent(in) :: ngrdcol ! # grid cols
1702 : integer, intent(in) :: lchnk ! chunk index
1703 : type(physics_ptend), intent(inout) :: ptend
1704 : ! Because we can't get a state passed in here, we might have to use values from the
1705 : ! subcolumn generation. This would make any conservation checks invalid if this
1706 : ! function is called after another parameterization... hmm.
1707 :
1708 0 : call subcol_ptend_avg_shr(ptend_sc, ngrdcol, lchnk, ptend, is_filter_set(), is_weight_set())
1709 :
1710 0 : end subroutine subcol_ptend_avg_SILHS
1711 :
1712 0 : subroutine subcol_SILHS_var_covar_driver &
1713 : ( ztodt, state_sc, ptend_sc, &
1714 : pbuf )
1715 :
1716 : ! This subroutine calculates microphysical effects on five variances and
1717 : ! covariances: rtp2, thlp2, wprtp, wpthlp, and rtpthlp.
1718 : !
1719 : ! This code is experimental!!
1720 :
1721 0 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field
1722 : #ifdef CLUBB_SGS
1723 : #ifdef SILHS
1724 : use subcol_utils, only: subcol_get_weight
1725 : use subcol_pack_mod, only: subcol_unpack, subcol_get_nsubcol
1726 : use clubb_api_module, only: T_in_K2thlm_api, &
1727 : init_pdf_params_api, &
1728 : copy_multi_pdf_params_to_single,&
1729 : pdf_parameter
1730 : use silhs_api_module, only: lh_microphys_var_covar_driver_api
1731 : #endif
1732 : #endif
1733 :
1734 : implicit none
1735 :
1736 : ! Parameters
1737 : ! This fill value is set to catch errors; it should not be read.
1738 : real(r8), parameter :: fillvalue = -999._r8
1739 :
1740 : ! Input Variables
1741 : real(r8), intent(in) :: ztodt ! model time increment
1742 : type(physics_state), intent(in) :: state_sc ! state for sub-columns
1743 : type(physics_ptend), intent(in) :: ptend_sc ! ptend for sub-columns
1744 :
1745 : ! Pointers
1746 : type(physics_buffer_desc), pointer :: pbuf(:)
1747 :
1748 : #ifdef CLUBB_SGS
1749 : #ifdef SILHS
1750 : ! Local Variables
1751 : integer :: lchnk, ngrdcol, igrdcol, isubcol, ns, k
1752 : integer, dimension(pcols) :: nsubcol
1753 : real(r8), dimension(pcols*psubcols) :: weights_packed
1754 : real(r8), dimension(pcols,psubcols) :: weights
1755 : real(r8), dimension(pcols,psubcols,pverp) :: rc_all, rv_all, rt_all, w_all, thl_all
1756 : real(r8), dimension(pcols,psubcols,pver ) :: s_all, t_all, zm_all, omega_all, pmid_all
1757 : real(r8), dimension(pcols,psubcols) :: phis_all
1758 : real(r8), dimension(pcols,psubcols,pver ) :: stend, ttend
1759 : real(r8), dimension(pcols,psubcols,pverp) :: thltend, qctend, qvtend
1760 :
1761 : real(r8), dimension(pcols,psubcols,pver) :: dz_g, pdel_all, rho
1762 : real(r8), dimension(pcols,psubcols,pverp) :: zi_all
1763 :
1764 : real(r8), dimension(pcols,psubcols,pver ) :: exner
1765 :
1766 : ! Inputs to lh_microphys_var_covar_driver
1767 : real(r8), dimension(pcols,psubcols,pverp) :: rt_all_clubb, thl_all_clubb, w_all_clubb, &
1768 : qctend_clubb, qvtend_clubb, thltend_clubb
1769 : real(r8), dimension(pcols,psubcols,pverp-top_lev+1) :: height_depndt_weights
1770 :
1771 : ! Outputs from lh_microphys_var_covar_driver
1772 : real(r8), dimension(:,:), pointer :: rtp2_mc_zt, thlp2_mc_zt, wprtp_mc_zt, &
1773 : wpthlp_mc_zt, rtpthlp_mc_zt
1774 :
1775 : ! pbuf indices
1776 : integer :: &
1777 : rtp2_mc_zt_idx, &
1778 : thlp2_mc_zt_idx, &
1779 : wprtp_mc_zt_idx, &
1780 : wpthlp_mc_zt_idx, &
1781 : rtpthlp_mc_zt_idx
1782 :
1783 : type(pdf_parameter) :: pdf_params_single_col
1784 :
1785 : !----- Begin Code -----
1786 :
1787 : call init_pdf_params_api( pverp+1-top_lev, 1, pdf_params_single_col )
1788 :
1789 : ! Don't do anything if this option isn't enabled.
1790 : if ( .not. subcol_SILHS_var_covar_src ) return
1791 :
1792 : lchnk = state_sc%lchnk
1793 : ngrdcol = state_sc%ngrdcol
1794 :
1795 : ! Obtain indices
1796 : rtp2_mc_zt_idx = pbuf_get_index('rtp2_mc_zt')
1797 : thlp2_mc_zt_idx = pbuf_get_index('thlp2_mc_zt')
1798 : wprtp_mc_zt_idx = pbuf_get_index('wprtp_mc_zt')
1799 : wpthlp_mc_zt_idx = pbuf_get_index('wpthlp_mc_zt')
1800 : rtpthlp_mc_zt_idx = pbuf_get_index('rtpthlp_mc_zt')
1801 :
1802 : ! Obtain pbuf fields for output
1803 : call pbuf_get_field(pbuf, rtp2_mc_zt_idx, rtp2_mc_zt)
1804 : call pbuf_get_field(pbuf, thlp2_mc_zt_idx, thlp2_mc_zt)
1805 : call pbuf_get_field(pbuf, wprtp_mc_zt_idx, wprtp_mc_zt)
1806 : call pbuf_get_field(pbuf, wpthlp_mc_zt_idx, wpthlp_mc_zt)
1807 : call pbuf_get_field(pbuf, rtpthlp_mc_zt_idx, rtpthlp_mc_zt)
1808 :
1809 : ! Unpack needed tendencies from subcolumn ptends
1810 : call subcol_unpack(lchnk, ptend_sc%s(:,:), stend, fillvalue)
1811 : call subcol_unpack(lchnk, ptend_sc%q(:,:,ixcldliq), qctend(:,:,1:pver), fillvalue)
1812 : call subcol_unpack(lchnk, ptend_sc%q(:,:,ixq), qvtend(:,:,1:pver), fillvalue)
1813 :
1814 : ! Unpack sample point values from subcolumn states
1815 : call subcol_unpack(lchnk, state_sc%q(:,:,ixcldliq), rc_all(:,:,1:pver), fillvalue)
1816 : call subcol_unpack(lchnk, state_sc%q(:,:,ixq), rv_all(:,:,1:pver), fillvalue)
1817 : call subcol_unpack(lchnk, state_sc%omega(:,:), omega_all (:,:,:), fillvalue)
1818 : call subcol_unpack(lchnk, state_sc%s(:,:), s_all, fillvalue)
1819 : call subcol_unpack(lchnk, state_sc%zm, zm_all, fillvalue)
1820 : call subcol_unpack(lchnk, state_sc%phis, phis_all, fillvalue)
1821 : call subcol_unpack(lchnk, state_sc%zi, zi_all, fillvalue)
1822 : call subcol_unpack(lchnk, state_sc%pdel, pdel_all, fillvalue)
1823 : call subcol_unpack(lchnk, state_sc%pmid, pmid_all, fillvalue)
1824 :
1825 : ! Initialize fields to fillvalue.
1826 : rt_all = fillvalue
1827 : thl_all = fillvalue
1828 : w_all = fillvalue
1829 : qctend = fillvalue
1830 : qvtend = fillvalue
1831 : thltend = fillvalue
1832 :
1833 : ! How many subcolumns in each column?
1834 : call subcol_get_nsubcol(lchnk, nsubcol)
1835 :
1836 : do igrdcol = 1, ngrdcol
1837 : do isubcol = 1, nsubcol(igrdcol)
1838 :
1839 : rt_all(igrdcol,isubcol,top_lev:pver) = rc_all(igrdcol,isubcol,top_lev:pver) &
1840 : + rv_all(igrdcol,isubcol,top_lev:pver)
1841 :
1842 : ! Compute dry static density on CLUBB vertical grid
1843 : do k = top_lev, pver
1844 : dz_g(igrdcol,isubcol,k) = zi_all(igrdcol,isubcol,k) - zi_all(igrdcol,isubcol,k+1) ! thickness
1845 : rho(igrdcol,isubcol,k) = (rga)*pdel_all(igrdcol,isubcol,k)/dz_g(igrdcol,isubcol,k)
1846 : end do
1847 :
1848 : ! Compute w from omega
1849 : w_all(igrdcol,isubcol,top_lev:pver) = -omega_all(igrdcol,isubcol,top_lev:pver) &
1850 : / ( rho(igrdcol,isubcol,top_lev:pver) * gravit )
1851 :
1852 : ! Convert stend and s_all to ttend and t_all
1853 : ! Note 1: With subcolumns, cpair is truly a constant (I think).
1854 : ! Note 2: For tendencies, the extra terns zm and phis should
1855 : ! not be included in the calculation.
1856 : ttend(igrdcol,isubcol,top_lev:pver) = stend(igrdcol,isubcol,top_lev:pver) / cpair
1857 :
1858 : do k = top_lev, pver
1859 : t_all(igrdcol,isubcol,k) = ( s_all(igrdcol,isubcol,k) &
1860 : - gravit * zm_all(igrdcol,isubcol,k) &
1861 : - phis_all(igrdcol,isubcol) ) / cpair
1862 : end do ! k = 1, pver
1863 :
1864 : ! This formula is taken from earlier in this file.
1865 : exner(igrdcol,isubcol,top_lev:pver) &
1866 : = ( pmid_all(igrdcol,isubcol,top_lev:pver) / p0_clubb )**(rair/cpair)
1867 :
1868 : ! Note: all tendencies or all means should be used in the call to
1869 : ! T_in_K2thlm_api (with the exception of exner)
1870 : do k = top_lev, pver
1871 : thltend(igrdcol,isubcol,k) &
1872 : = T_in_K2thlm_api( ttend(igrdcol,isubcol,k), exner(igrdcol,isubcol,k), &
1873 : qctend(igrdcol,isubcol,k) )
1874 : thl_all(igrdcol,isubcol,k) &
1875 : = T_in_K2thlm_api( t_all(igrdcol,isubcol,k), exner(igrdcol,isubcol,k), &
1876 : rc_all(igrdcol,isubcol,k) )
1877 : end do ! k = 1, pver
1878 :
1879 : ! Add ghost points
1880 : rt_all (igrdcol,isubcol,pverp) = rt_all (igrdcol,isubcol,pver)
1881 : thl_all(igrdcol,isubcol,pverp) = thl_all(igrdcol,isubcol,pver)
1882 : w_all (igrdcol,isubcol,pverp) = w_all (igrdcol,isubcol,pver)
1883 : qctend (igrdcol,isubcol,pverp) = qctend (igrdcol,isubcol,pver)
1884 : qvtend (igrdcol,isubcol,pverp) = qvtend (igrdcol,isubcol,pver)
1885 : thltend(igrdcol,isubcol,pverp) = thltend(igrdcol,isubcol,pver)
1886 :
1887 : ! Flip inputs to CLUBB's grid. Note the dimension ordering change.
1888 : rt_all_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( rt_all(igrdcol,isubcol,1:pverp) )
1889 : thl_all_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( thl_all(igrdcol,isubcol,1:pverp) )
1890 : w_all_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( w_all(igrdcol,isubcol,1:pverp) )
1891 : qctend_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( qctend(igrdcol,isubcol,1:pverp) )
1892 : qvtend_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( qvtend(igrdcol,isubcol,1:pverp) )
1893 : thltend_clubb(igrdcol,isubcol,1:pverp) = clubb_flip_grid( thltend(igrdcol,isubcol,1:pverp) )
1894 :
1895 : end do ! isubcol = 1, nsubcol(igrdcol)
1896 : end do ! igrdcol = 1, ngrdcol
1897 :
1898 : ! Obtain weights
1899 : call subcol_get_weight(lchnk, weights_packed)
1900 : call subcol_unpack(lchnk, weights_packed, weights, fillvalue)
1901 :
1902 : ! Call lh_microphys_var_covar_driver for each column
1903 : do igrdcol=1, ngrdcol
1904 : ns = nsubcol(igrdcol)
1905 :
1906 : ! This code assumes that the weights are height independent.
1907 : ! It will have to change once the weights vary with altitude!
1908 : ! I'm not sure whether the grid will need to be flipped.
1909 : do k = 1, pverp-top_lev+1
1910 : height_depndt_weights(igrdcol,1:ns,k) = weights(igrdcol,1:ns)
1911 : end do
1912 :
1913 : ! Copy the igrdcol column from the multicolumn pdf_params_chnk to the single column
1914 : ! version of pdf_params_single_col since lh_microphys_var_covar_driver_api only
1915 : ! works over 1 column currently
1916 : call copy_multi_pdf_params_to_single( pdf_params_chnk(lchnk), igrdcol, &
1917 : pdf_params_single_col )
1918 :
1919 : ! Make the call!!!!!
1920 : call lh_microphys_var_covar_driver_api &
1921 : ( pverp-top_lev+1, ns, ztodt, height_depndt_weights(igrdcol,1:ns,1:pverp-top_lev+1), &
1922 : pdf_params_single_col, &
1923 : rt_all_clubb(igrdcol,1:ns,1:pverp-top_lev+1), thl_all_clubb(igrdcol,1:ns,1:pverp-top_lev+1), &
1924 : w_all_clubb(igrdcol,1:ns,1:pverp-top_lev+1), qctend_clubb(igrdcol,1:ns,1:pverp-top_lev+1), &
1925 : qvtend_clubb(igrdcol,1:ns,1:pverp-top_lev+1), thltend_clubb(igrdcol,1:ns,1:pverp-top_lev+1), &
1926 : silhs_config_flags%l_lh_instant_var_covar_src, &
1927 : rtp2_mc_zt(igrdcol,1:pverp-top_lev+1), thlp2_mc_zt(igrdcol,1:pverp-top_lev+1), &
1928 : wprtp_mc_zt(igrdcol,1:pverp-top_lev+1), wpthlp_mc_zt(igrdcol,1:pverp-top_lev+1), &
1929 : rtpthlp_mc_zt(igrdcol,1:pverp-top_lev+1) )
1930 :
1931 : ! The *_mc_zt microphysics tendencies are passed out of SILHS and back
1932 : ! to CLUBB without being used at all in the rest of the host model code.
1933 : ! The arrays aren't flipped for the *_mc_zt microphysics tendencies, and
1934 : ! they don't need to be.
1935 :
1936 : ! CLUBB used pverp vertical levels, but SILHS only uses
1937 : ! pverp - top_lev + 1 vertical levels.
1938 : ! Fill the upper levels with 0s when necessary.
1939 : if ( pverp > pverp-top_lev+1 ) then
1940 : rtp2_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
1941 : thlp2_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
1942 : wprtp_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
1943 : wpthlp_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
1944 : rtpthlp_mc_zt(igrdcol,pverp-top_lev+2:pverp) = 0.0_r8
1945 : endif ! pverp > pverp-top_lev+1
1946 :
1947 : end do ! igrdcol = 1, ngrdcol
1948 : #endif
1949 : #endif
1950 :
1951 0 : return
1952 0 : end subroutine subcol_SILHS_var_covar_driver
1953 : #ifdef SILHS
1954 : real(r8) function meansc(arr_in, w_in, ns) result(val)
1955 : real(r8), intent(in) :: ns ! Length of Array
1956 : real(r8), dimension(int(ns)), intent(in) :: arr_in ! Input array
1957 : real(r8), dimension(int(ns)), intent(in) :: w_in ! Weights
1958 : real(r8) :: acc ! accumulator
1959 : integer :: i
1960 : acc = 0
1961 : val = 0
1962 : do i=1,ns
1963 : acc = acc + arr_in(i)*w_in(i)
1964 : end do
1965 : val = acc/ns
1966 : end function
1967 :
1968 : real(r8) function stdsc(arr_in, w_in, mn_in, ns) result(val)
1969 : real(r8), intent(in) :: ns ! Number of elements (subcolumns)
1970 : real(r8), dimension(int(ns)), intent(in) :: arr_in, w_in !Input array and weights
1971 : real(r8), intent(in) :: mn_in ! The mean of arr_in
1972 : real(r8) :: accvar, var
1973 : integer :: i
1974 : accvar = 0
1975 : do i=1,ns
1976 : accvar = accvar + ((arr_in(i)-mn_in)**2)*w_in(i)
1977 : end do
1978 : var = accvar/ns
1979 : val = sqrt(var)
1980 : end function
1981 :
1982 : subroutine THL_profile(nz, ABST_prof, ex_prof, rcm_prof, THL_prof)
1983 :
1984 : use clubb_api_module, only : T_in_K2thlm_api
1985 :
1986 : integer, intent(in) :: nz ! Num vert levels
1987 : real(r8), dimension(nz), intent(in) :: ABST_prof ! Abs Temp prof
1988 : real(r8), dimension(nz), intent(in) :: ex_prof ! Profile of Exner func
1989 : real(r8), dimension(nz), intent(in) :: rcm_prof ! Profile of Cld Wat MR
1990 : real(r8), dimension(nz), intent(out) :: THL_prof ! LWPT prof
1991 : integer :: i
1992 :
1993 : do i=1,nz
1994 : THL_prof(i) = T_in_K2thlm_api(ABST_prof(i), ex_prof(i), rcm_prof(i))
1995 : end do
1996 :
1997 : end subroutine
1998 :
1999 : subroutine subcol_constrainmn( num_subcols, samples, weights, grid_mean, mean_sc, std_sc )
2000 :
2001 : ! Input/Output Variables
2002 : integer, intent(in) :: num_subcols
2003 : real(r8), dimension(num_subcols, pverp), intent(inout) :: samples
2004 : real(r8), dimension(num_subcols), intent(in) :: weights
2005 : real(r8), dimension(pverp), intent(in) :: grid_mean
2006 : real(r8), dimension(pver), intent(out), optional :: mean_sc, std_sc
2007 :
2008 : ! Local Variables
2009 : real(r8) :: meansc_loc, adj_rat
2010 : integer :: k
2011 : !------------------------------------------------------------------
2012 : !----- Begin Code -----
2013 : do k=1, pver
2014 : meansc_loc = meansc( samples(:,k), weights(:), real(num_subcols, r8) )
2015 :
2016 : if (present(mean_sc)) &
2017 : mean_sc(k) = meansc_loc
2018 : if (present(std_sc)) &
2019 : std_sc(k) = stdsc( samples(:,k), weights(:), meansc_loc, &
2020 : real(num_subcols, r8) )
2021 :
2022 : if ( meansc_loc > 0.0_r8 ) then
2023 : adj_rat = grid_mean(k)/meansc_loc
2024 : else
2025 : ! If the mean is zero, then zero out all subcolumns to avoid
2026 : ! negative samples
2027 : adj_rat = 0.0_r8
2028 : end if
2029 : samples(:,k) = samples(:,k) * adj_rat
2030 : end do
2031 : end subroutine subcol_constrainmn
2032 :
2033 : ! =============================================================================== !
2034 : ! !
2035 : ! =============================================================================== !
2036 : function clubb_flip_grid ( profile ) result( profile_flipped )
2037 :
2038 : ! Description:
2039 : ! Swaps the elements in profile so they are in reverse order. CAM and
2040 : ! CLUBB's grids are flipped with respect to each other.
2041 : !
2042 : ! Usage:
2043 : ! clubb_var = clubb_flip_grid( cam_var )
2044 : ! cam_var = clubb_flip_grid( clubb_var )
2045 :
2046 : implicit none
2047 :
2048 : ! Input Variable
2049 : real(r8), dimension(pverp), intent(in) :: profile
2050 :
2051 : ! Output Variable
2052 : real(r8), dimension(pverp) :: profile_flipped
2053 :
2054 : ! Local Variable
2055 : integer :: k
2056 :
2057 : do k=1, pverp
2058 : profile_flipped(k) = profile(pverp-k+1)
2059 : end do ! k=1, pverp
2060 :
2061 : return
2062 : end function clubb_flip_grid
2063 : ! =============================================================================== !
2064 : ! !
2065 : ! =============================================================================== !
2066 : #endif
2067 : !============================================================================
2068 0 : subroutine subcol_SILHS_fill_holes_conserv( state, dt, ptend, pbuf )
2069 :
2070 : ! The William F. Buckley Jr. Conservative Hole Filler.
2071 :
2072 : ! Description:
2073 : ! Stops holes from forming in a hydrometeor mixing ratio by reducing the
2074 : ! microphysics tendency of that hydrometeor mixing ratio which would
2075 : ! otherwise cause that hydrometeor mixing ratio to have a negative value
2076 : ! once the microphysics tendency is applied. This code is used to prevent
2077 : ! holes in water mass, not number concentration.
2078 : !
2079 : ! This subroutine is called after microphysics has completed and after
2080 : ! microphysics fields from subcolumns have been averaged back to grid
2081 : ! columns, but before the grid-column microphysics tendencies have been
2082 : ! applied in physics_update. This code is meant for use with the SILHS
2083 : ! subcolumn approach. This code needs to be applied to grid columns, not
2084 : ! subcolumns.
2085 : !
2086 : ! This code adjusts the tendencies (ptend) before they are used to update
2087 : ! the grid mean fields (state variables).
2088 : !
2089 : ! The column-integrated total water needs to be conserved during
2090 : ! microphysics. The conserved amount includes the amount of water that
2091 : ! precipitated to the ground from sedimentation during microphysics.
2092 : ! The conservation equation for each grid column is:
2093 : !
2094 : ! SUM(k=top_lev:pver) ( rv_start(k) + rc_start(k) + rr_start(k)
2095 : ! + ri_start(k) + rs_start(k) ) * pdel(k) / g
2096 : ! = SUM(k=top_lev:pver) ( rv(k) + rc(k) + rr(k) + ri(k) + rs(k) )
2097 : ! * pdel(k) / g
2098 : ! + prect * dt * 1000;
2099 : !
2100 : ! where rv_start, rc_start, rr_start, ri_start, and rs_start are water
2101 : ! vapor, cloud water, rain water, cloud ice, and snow mixing ratios before
2102 : ! microphysics is called; rv, rc, rr, ri, and rs are water vapor, cloud
2103 : ! water, rain water, cloud ice, and snow mixing ratios after being updated
2104 : ! by microphysics; pdel is the pressure difference between vertical levels,
2105 : ! g is gravity, and prect * dt * 1000 is the total amount of water (from
2106 : ! all precipitating hydrometeors) that sedimented to the ground during
2107 : ! microphysics (dt is the timestep used for microphysics). The units of
2108 : ! column-integrated total water are kg (water) / m^2.
2109 : !
2110 : ! All the updated hydrometeor fields are related to the hydrometeor fields
2111 : ! at the start by:
2112 : !
2113 : ! rv(k) = rv_start(k) + rv_tend(k) * dt;
2114 : ! rc(k) = rc_start(k) + rc_tend(k) * dt;
2115 : ! rr(k) = rr_start(k) + rr_tend(k) * dt;
2116 : ! ri(k) = ri_start(k) + ri_tend(k) * dt; and
2117 : ! rs(k) = rs_start(k) + rs_tend(k) * dt;
2118 : !
2119 : ! where rv_tend, rc_tend, rr_tend, ri_tend, and rs_tend are water vapor,
2120 : ! cloud water, rain water, cloud ice, and snow mixing ratio tendencies
2121 : ! from microphysics, which includes the sum of microphysics process rates
2122 : ! and sedimentation. When these equations are applied to the equation
2123 : ! for column-integrated total water, that equation becomes:
2124 : !
2125 : ! SUM(k=top_lev:pver) ( rv_tend(k) + rc_tend(k) + rr_tend(k)
2126 : ! + ri_tend(k) + rs_tend(k) ) * dt * pdel(k) / g
2127 : ! + prect * dt * 1000 = 0.
2128 : !
2129 : ! As stated above, the hydrometeor tendencies are the sum of tendencies
2130 : ! from microphysics process rates and tendencies from sedimentation:
2131 : !
2132 : ! rv_tend(k) = rv_mc_tend(k);
2133 : ! rc_tend(k) = rc_mc_tend(k) + rc_sed_tend(k);
2134 : ! rr_tend(k) = rr_mc_tend(k) + rr_sed_tend(k);
2135 : ! ri_tend(k) = ri_mc_tend(k) + ri_sed_tend(k); and
2136 : ! rs_tend(k) = rs_mc_tend(k) + rs_sed_tend(k);
2137 : !
2138 : ! where rv_mc_tend, rc_mc_tend, rr_mc_tend, ri_mc_tend, and rs_mc_tend are
2139 : ! the tendencies of water vapor, cloud water, rain water, cloud ice, and
2140 : ! snow from microphysics process rates, and rc_sed_tend, rr_sed_tend,
2141 : ! ri_sed_tend, and rs_sed_tend are the tendencies of cloud water,
2142 : ! rain water, cloud ice, and snow from sedimentation. When these equations
2143 : ! are applied to the equation for column-integrated total water, that
2144 : ! equation becomes:
2145 : !
2146 : ! SUM(k=top_lev:pver) ( rv_mc_tend(k) + rc_mc_tend(k) + rr_mc_tend(k)
2147 : ! + ri_mc_tend(k) + rs_mc_tend(k) )
2148 : ! * dt * pdel(k) / g
2149 : ! + SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) + ri_sed_tend(k)
2150 : ! + rs_sed_tend(k) ) * dt * pdel(k) / g
2151 : ! + prect * dt * 1000 = 0.
2152 : !
2153 : ! At any vertical level, the tendencies from microphysics process rates
2154 : ! (mc_tend variables) must balance:
2155 : !
2156 : ! rv_mc_tend(k) + rc_mc_tend(k) + rr_mc_tend(k)
2157 : ! + ri_mc_tend(k) + rs_mc_tend(k) = 0; for all k from top_lev to pver.
2158 : !
2159 : ! The column-integrated total water equation can be applied to
2160 : ! sedimentation:
2161 : !
2162 : ! SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) + ri_sed_tend(k)
2163 : ! + rs_sed_tend(k) ) * dt * pdel(k) / g
2164 : ! + prect * dt * 1000 = 0.
2165 : !
2166 : ! The total precipitation rate, prect, can be split into liquid
2167 : ! precipitation rate, precl, and frozen precipitation rate, preci:
2168 : !
2169 : ! prect = precl + preci.
2170 : !
2171 : ! The microphysics code outputs prect and preci, so precl can be calculated
2172 : ! by precl = prect - preci. The column-integrated total water equation can
2173 : ! be split into:
2174 : !
2175 : ! SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) )
2176 : ! * dt * pdel(k) / g
2177 : ! + precl * dt * 1000 = 0; and
2178 : !
2179 : ! SUM(k=top_lev:pver) ( ri_sed_tend(k) + rs_sed_tend(k) )
2180 : ! * dt * pdel(k) / g
2181 : ! + preci * dt * 1000 = 0.
2182 : !
2183 : ! Overall, the conservation methods used in this subroutine are:
2184 : !
2185 : ! 1) When adjusting the tendencies from microphysics process rates,
2186 : ! conserve:
2187 : !
2188 : ! rv_mc_tend(k) + rc_mc_tend(k) + rr_mc_tend(k)
2189 : ! + ri_mc_tend(k) + rs_mc_tend(k) = 0; for all k from top_lev to pver.
2190 : !
2191 : ! 2) When adjusting the tendencies from microphysics process rates, adjust
2192 : ! dry static energy appropriately. The change in dry static energy
2193 : ! is necessary because of phase changes. This "puts back" the extra dry
2194 : ! static energy that was "taken out" when an excessive phase-changing
2195 : ! process rate was produced by microphysics.
2196 : !
2197 : ! 3) When adjusting the hydrometeor tendency from sedimentation of a
2198 : ! liquid hydrometeor (cloud water or rain water), conserve:
2199 : !
2200 : ! SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) )
2201 : ! * dt * pdel(k) / g
2202 : ! + precl * dt * 1000 = 0.
2203 : !
2204 : ! 4) When adjusting the hydrometeor tendency from sedimentation of a
2205 : ! frozen hydrometeor (cloud ice or snow), conserve:
2206 : !
2207 : ! SUM(k=top_lev:pver) ( ri_sed_tend(k) + rs_sed_tend(k) )
2208 : ! * dt * pdel(k) / g
2209 : ! + preci * dt * 1000 = 0.
2210 : !
2211 : ! The conservative hole filler works as follows. The total microphysics
2212 : ! tendency for each hydrometeor is provided in ptend. This is the sum of
2213 : ! the microphysics process rate tendency and sedimentation tendency for
2214 : ! each hydrometeor. The sedimentation tendency is provided in pbuf. The
2215 : ! sedimentation tendency is subtracted off the total microphysics tendency
2216 : ! to produce the microphysics process rate tendency for each hydrometeor.
2217 : ! The microphysics process rate tendency is adjusted when necessary so that
2218 : ! holes in the hydrometeor are not produced by microphysics process rates.
2219 : ! When a hydrometeor's negative microphysics process rate tendency needs to
2220 : ! be made smaller in magnitude to avoid a hole, all hydrometeor tendencies
2221 : ! that are positive at that grid level are also decreased proportionately
2222 : ! to maintain a balance. Dry static energy tendency is also adjusted
2223 : ! appropriately when necessary. After this, the vertical integral of each
2224 : ! hydrometeor species is greater than or equal to 0.
2225 : !
2226 : ! The sedimentation tendency is then added back onto the new microphysics
2227 : ! process rate tendency to produce a new total microphysics tendency for
2228 : ! each hydrometeor. Since the sedimentation tendency was based on the old
2229 : ! value of hydrometeor, before the hole-filling adjustment, it is possible
2230 : ! that the new total microphysics tendency may produce holes. When this
2231 : ! happens, sedimentation hole filling fills holes in the vertical profile
2232 : ! of each hydrometeor. Holes are filled using mass from other vertical
2233 : ! levels for the same hydrometeor (or from a same-phase hydrometeor when
2234 : ! necessary). Since the vertical integral of sedimentation tendency
2235 : ! (including surface precipitation rate) is 0, the vertical integral of the
2236 : ! hydrometeor must be greater than or equal to 0, which means that all
2237 : ! holes can be filled. The result is that all holes in any hydrometeor
2238 : ! mixing ratio are filled completely and conservatively. The value of
2239 : ! ptend is updated appropriately so that it can be applied later in
2240 : ! physics_update.
2241 :
2242 : !----------------------------------------------------------------------
2243 :
2244 0 : use physics_buffer, only: &
2245 : physics_buffer_desc, &
2246 : pbuf_get_field
2247 :
2248 : use ppgrid, only: &
2249 : pcols
2250 :
2251 : use constituents, only: &
2252 : qmin
2253 :
2254 : use ref_pres, only: &
2255 : top_lev => trop_cloud_top_lev
2256 :
2257 : implicit none
2258 :
2259 : ! Input Variables
2260 : type(physics_state), intent(in) :: state ! Physics state variables
2261 : real(r8), intent(in) :: dt ! Time step duration
2262 :
2263 : ! Input/Output Variables
2264 : type(physics_ptend), intent(inout) :: ptend ! Parameterization tendencies
2265 : type(physics_buffer_desc), pointer :: pbuf(:) ! Physics buffer
2266 :
2267 : ! Local Variables
2268 : real(r8), dimension(pcols,pver) :: &
2269 : rv_start, & ! Water vapor mixing ratio at start of microphysics [kg/kg]
2270 : rc_start, & ! Cloud water mixing ratio at start of microphysics [kg/kg]
2271 : rr_start, & ! Rain water mixing ratio at start of microphysics [kg/kg]
2272 : ri_start, & ! Cloud ice mixing ratio at start of microphysics [kg/kg]
2273 : rs_start ! Snow mixing ratio at start of microphysics [kg/kg]
2274 :
2275 : real(r8), dimension(pcols,pver) :: &
2276 : rv_tend, & ! Water vapor mixing ratio tendency [kg/kg/s]
2277 : rc_tend, & ! Cloud water mixing ratio tendency [kg/kg/s]
2278 : rr_tend, & ! Rain water mixing ratio tendency [kg/kg/s]
2279 : ri_tend, & ! Cloud ice mixing ratio tendency [kg/kg/s]
2280 : rs_tend, & ! Snow mixing ratio tendency [kg/kg/s]
2281 : stend ! Dry static energy tendency [J/kg/s]
2282 :
2283 : real(r8), dimension(:), pointer :: &
2284 0 : prect, & ! Total microphysics precipitation rate (surface) [m/s]
2285 0 : preci, & ! Ice-phase microphysics precipitation rate (surface) [m/s]
2286 0 : prec_str, & ! Total surface precipitation rate from stratoform [m/s]
2287 0 : snow_str ! Snow surface precipitation rate from stratoform [m/s]
2288 :
2289 : real(r8), dimension(:,:), pointer :: &
2290 0 : rc_sed_tend, & ! Mean cloud water sedimentation tendency [kg/kg/s]
2291 0 : rr_sed_tend, & ! Mean rain water sedimentation tendency [kg/kg/s]
2292 0 : ri_sed_tend, & ! Mean cloud ice sedimentation tendency [kg/kg/s]
2293 0 : rs_sed_tend, & ! Mean snow sedimentation tendency [kg/kg/s]
2294 0 : vtrmc, & ! Mean cloud water sedimentation velocity [m/s]
2295 0 : umr, & ! Mean rain water sedimentation velocity [m/s]
2296 0 : vtrmi, & ! Mean cloud ice sedimentation velocity [m/s]
2297 0 : ums, & ! Mean snow sedimentation velocity [m/s]
2298 0 : rc_sed_evap, & ! Mean evap of cloud water during sedimentation [kg/kg/s]
2299 0 : ri_sed_subl ! Mean subl of cloud ice during sedimentation [kg/kg/s]
2300 :
2301 : real(r8), dimension(pcols,pver) :: &
2302 : rv_mc_tend, & ! Water vapor mixing ratio microphysics tendency [kg/kg/s]
2303 : rc_mc_tend, & ! Cloud water mixing ratio microphysics tendency [kg/kg/s]
2304 : rr_mc_tend, & ! Rain water mixing ratio microphysics tendency [kg/kg/s]
2305 : ri_mc_tend, & ! Cloud ice mixing ratio microphysics tendency [kg/kg/s]
2306 : rs_mc_tend ! Snow mixing ratio microphysics tendency [kg/kg/s]
2307 :
2308 : real(r8) :: &
2309 : rv_curr, & ! Current water vapor mixing ratio [kg/kg]
2310 : rc_curr, & ! Current cloud water mixing ratio [kg/kg]
2311 : rr_curr, & ! Current rain water mixing ratio [kg/kg]
2312 : ri_curr, & ! Current cloud ice mixing ratio [kg/kg]
2313 : rs_curr ! Current snow mixing ratio [kg/kg]
2314 :
2315 : logical :: &
2316 : l_pos_rv_mc_tend, & ! Flag for positive water vapor mixing ratio mc tend.
2317 : l_pos_rc_mc_tend, & ! Flag for positive cloud water mixing ratio mc tend.
2318 : l_pos_rr_mc_tend, & ! Flag for positive rain water mixing ratio mc tend.
2319 : l_pos_ri_mc_tend, & ! Flag for positive cloud ice mixing ratio mc tend.
2320 : l_pos_rs_mc_tend ! Flag for positive snow mixing ratio mc tend.
2321 :
2322 : real(r8) :: &
2323 : mc_tend_max_mag, & ! Max. allowable mag. of (neg.) mc tend [kg/kg/s]
2324 : mc_tend_correction, & ! Amnt. correction necessary to mc tend [kg/kg/s]
2325 : total_mc_positive, & ! Total of all positive mc tendencies [kg/kg/s]
2326 : mc_correction_ratio ! Ratio: mc_tend_correction/total_mc_positive [-]
2327 :
2328 : real(r8), dimension(pcols) :: &
2329 : precl ! Liquid-phase precipitation rate (surface) [m/s]
2330 :
2331 : ! Budgeting terms for hole filling.
2332 : ! These variables are for use in stats output.
2333 : real(r8), dimension(pcols,pver) :: &
2334 : rv_hf_tend, & ! Water vapor mixing ratio hole-filling tendency [kg/kg/s]
2335 : rc_hf_tend, & ! Cloud water mixing ratio hole-filling tendency [kg/kg/s]
2336 : rr_hf_tend, & ! Rain water mixing ratio hole-filling tendency [kg/kg/s]
2337 : ri_hf_tend, & ! Cloud ice mixing ratio hole-filling tendency [kg/kg/s]
2338 : rs_hf_tend, & ! Snow mixing ratio hole-filling tendency [kg/kg/s]
2339 : s_hf_tend ! Dry static energy hole-filling tendency [J/kg/s]
2340 :
2341 : integer :: ncol ! Number of grid columns
2342 :
2343 : integer :: icol, k ! Loop indices
2344 :
2345 : ! Flag to perform hole filling after the original sedimentation tendency
2346 : ! is added back on to the new microphysics process tendency. This calls
2347 : ! the sedimentation hole filler.
2348 : logical, parameter :: &
2349 : l_sed_hole_fill = .true.
2350 :
2351 : logical, parameter :: &
2352 : l_check_conservation = .true. ! Flag to perform water conservation check
2353 :
2354 : ! Vertically-integrated grand total water (rv+rc+rr+ri+rs) [kg/m^2]
2355 : real(r8), dimension(pcols) :: &
2356 : grand_total_water_column_start, & ! Column integral at start
2357 : grand_total_water_column_finish ! Column integral at finish
2358 :
2359 : ! Vertically-integrated total water energy [J/m^2]
2360 : real(r8), dimension(pcols) :: &
2361 : total_energy_column_start, & ! Column integral at start
2362 : total_energy_column_finish ! Column integral at finish
2363 :
2364 : real(r8), dimension(pcols) :: &
2365 : tot_water_rel_err, & ! Relative error: vert-integrated grand total water
2366 : tot_energy_rel_err ! Relative error: vert-integrated total energy
2367 :
2368 : real(r8), parameter :: &
2369 : err_thresh = 1.0e-14_r8 ! Threshold of relative error
2370 :
2371 :
2372 : ! Get the number of grid columns.
2373 0 : ncol = state%ncol
2374 :
2375 : ! Get fields from the pbuf.
2376 0 : call pbuf_get_field(pbuf, prec_pcw_idx, prect)
2377 0 : call pbuf_get_field(pbuf, snow_pcw_idx, preci)
2378 0 : call pbuf_get_field(pbuf, prec_str_idx, prec_str)
2379 0 : call pbuf_get_field(pbuf, snow_str_idx, snow_str)
2380 0 : call pbuf_get_field(pbuf, qcsedten_idx, rc_sed_tend)
2381 0 : call pbuf_get_field(pbuf, qrsedten_idx, rr_sed_tend)
2382 0 : call pbuf_get_field(pbuf, qisedten_idx, ri_sed_tend)
2383 0 : call pbuf_get_field(pbuf, qssedten_idx, rs_sed_tend)
2384 0 : call pbuf_get_field(pbuf, vtrmc_idx, vtrmc)
2385 0 : call pbuf_get_field(pbuf, umr_idx, umr)
2386 0 : call pbuf_get_field(pbuf, vtrmi_idx, vtrmi)
2387 0 : call pbuf_get_field(pbuf, ums_idx, ums)
2388 0 : call pbuf_get_field(pbuf, qcsevap_idx, rc_sed_evap)
2389 0 : call pbuf_get_field(pbuf, qisevap_idx, ri_sed_subl)
2390 :
2391 : ! Calculate liquid precipitation rate (precl) from the total precipitation
2392 : ! rate (prect) and the frozen preciptation rate (preci). This should never
2393 : ! be negative, but just to be safe, threshold at 0.
2394 0 : precl(:ncol) = max( prect(:ncol) - preci(:ncol), 0.0_r8 )
2395 :
2396 : ! Perform total water and total energy conservation checks.
2397 : if ( l_check_conservation ) then
2398 :
2399 : ! Calculate total water in each column.
2400 : ! This calculation is the vertically-integrated grand total water (where
2401 : ! grand total water is the sum of water vapor, cloud water, rain water,
2402 : ! cloud ice, and snow, as well as the amount of water that precipitated
2403 : ! to the surface) in each grid column after microphysics, but at the
2404 : ! start of hole filling.
2405 0 : do icol = 1, ncol
2406 0 : grand_total_water_column_start(icol) = 0.0_r8
2407 0 : do k = top_lev, pver
2408 : grand_total_water_column_start(icol) &
2409 : = grand_total_water_column_start(icol) &
2410 0 : + ( state%q(icol,k,1) + ptend%q(icol,k,1) * dt &
2411 0 : + state%q(icol,k,ixcldliq) &
2412 0 : + ptend%q(icol,k,ixcldliq) * dt &
2413 0 : + state%q(icol,k,ixcldice) &
2414 0 : + ptend%q(icol,k,ixcldice) * dt ) &
2415 0 : * state%pdel(icol,k) * rga
2416 0 : if ( ixrain > 0 ) then
2417 : grand_total_water_column_start(icol) &
2418 : = grand_total_water_column_start(icol) &
2419 0 : + ( state%q(icol,k,ixrain) + ptend%q(icol,k,ixrain) * dt ) &
2420 0 : * state%pdel(icol,k) * rga
2421 : endif
2422 0 : if ( ixsnow > 0 ) then
2423 : grand_total_water_column_start(icol) &
2424 : = grand_total_water_column_start(icol) &
2425 0 : + ( state%q(icol,k,ixsnow) + ptend%q(icol,k,ixsnow) * dt ) &
2426 0 : * state%pdel(icol,k) * rga
2427 : endif
2428 : end do ! k = top_lev, pver
2429 : grand_total_water_column_start(icol) &
2430 : = grand_total_water_column_start(icol) &
2431 0 : + prect(icol) * dt * 1000.0_r8
2432 : end do ! icol = 1, ncol
2433 :
2434 : ! Calculate total energy in each column.
2435 : ! This calculation is the vertically-integrated total energy in each
2436 : ! grid column after microphysics, but at the start of hole filling.
2437 : ! Since the microphysics and hole filling code does not directly change
2438 : ! kinetic energy, 0.5 * ( u^2 + v^2 ), it can be skipped as part of the
2439 : ! energy conservation check.
2440 0 : do icol = 1, ncol
2441 0 : total_energy_column_start(icol) = 0.0_r8
2442 0 : do k = top_lev, pver
2443 : total_energy_column_start(icol) &
2444 : = total_energy_column_start(icol) &
2445 0 : + ( state%s(icol,k) + ptend%s(icol,k) * dt &
2446 : + ( latvap + latice ) &
2447 0 : * ( state%q(icol,k,1) + ptend%q(icol,k,1) * dt ) &
2448 0 : + latice * ( state%q(icol,k,ixcldliq) &
2449 0 : + ptend%q(icol,k,ixcldliq) * dt ) ) &
2450 0 : * state%pdel(icol,k) * rga
2451 0 : if ( ixrain > 0 ) then
2452 : total_energy_column_start(icol) &
2453 : = total_energy_column_start(icol) &
2454 0 : + latice * ( state%q(icol,k,ixrain) &
2455 0 : + ptend%q(icol,k,ixrain) * dt ) &
2456 0 : * state%pdel(icol,k) * rga
2457 : endif
2458 : end do ! k = top_lev, pver
2459 : total_energy_column_start(icol) &
2460 : = total_energy_column_start(icol) &
2461 0 : + latice * precl(icol) * dt * 1000.0_r8
2462 : end do ! icol = 1, ncol
2463 :
2464 : endif ! l_check_conservation
2465 :
2466 : ! The fields within state haven't been updated yet, since this is before
2467 : ! the call to physics_update.
2468 0 : rv_start = state%q(:,:,1)
2469 0 : rc_start = state%q(:,:,ixcldliq)
2470 0 : if ( ixrain > 0 ) then
2471 0 : rr_start = state%q(:,:,ixrain)
2472 : endif
2473 0 : ri_start = state%q(:,:,ixcldice)
2474 0 : if ( ixsnow > 0 ) then
2475 0 : rs_start = state%q(:,:,ixsnow)
2476 : endif
2477 :
2478 : ! Unpack the current total tendencies for hydrometeor mixing ratio fields.
2479 0 : rv_tend = ptend%q(:,:,1)
2480 0 : rc_tend = ptend%q(:,:,ixcldliq)
2481 0 : if ( ixrain > 0 ) then
2482 0 : rr_tend = ptend%q(:,:,ixrain)
2483 : endif
2484 0 : ri_tend = ptend%q(:,:,ixcldice)
2485 0 : if ( ixsnow > 0 ) then
2486 0 : rs_tend = ptend%q(:,:,ixsnow)
2487 : endif
2488 :
2489 : ! Unpack the current tendency for dry static energy.
2490 0 : stend = ptend%s
2491 :
2492 : ! The total hydrometeor tendencies are the sum of microphysics process
2493 : ! rates and sedimentation rates. Calculate the microphysics process
2494 : ! tendencies by subtracting the sedimentation tendencies from the overall
2495 : ! tendencies.
2496 : ! The sedimentation tendencies for cloud water (rc_sed_tend) and cloud ice
2497 : ! (ri_sed_tend) include the evaporation of cloud water during sedimentation
2498 : ! and the sublimation of cloud ice during sedimentation, respectively. The
2499 : ! true sedimentation of cloud water is the sum of rc_sed_tend and
2500 : ! rc_sed_evap, and the true sedimentation of cloud ice is the sum of
2501 : ! ri_sed_tend and ri_sed_subl. Subtract off only the true sedimentation
2502 : ! rates, as evaporation and sublimation need to be included in the
2503 : ! microphysics process rates.
2504 0 : rv_mc_tend(:ncol,:) = rv_tend(:ncol,:)
2505 0 : rc_mc_tend(:ncol,:) = rc_tend(:ncol,:) - ( rc_sed_tend(:ncol,:) + rc_sed_evap(:ncol,:) )
2506 0 : if ( ixrain > 0 ) then
2507 0 : rr_mc_tend(:ncol,:) = rr_tend(:ncol,:) - rr_sed_tend(:ncol,:)
2508 : endif
2509 0 : ri_mc_tend(:ncol,:) = ri_tend(:ncol,:) - ( ri_sed_tend(:ncol,:) + ri_sed_subl(:ncol,:) )
2510 0 : if ( ixsnow > 0 ) then
2511 0 : rs_mc_tend(:ncol,:) = rs_tend(:ncol,:) - rs_sed_tend(:ncol,:)
2512 : endif
2513 :
2514 : ! This section adjusts microphysics process rate tendencies so that the
2515 : ! resulting values of all hydrometeor mixing ratios are greater than or
2516 : ! equal to qmin after this section is complete. Once sedimentation is
2517 : ! added back on after this section, some of the hydrometeor mixing ratios
2518 : ! may become less than qmin again.
2519 : !
2520 : ! This section, which again is concerned only with adjusting microphysics
2521 : ! process rates, makes use of the following two principles:
2522 : !
2523 : ! 1) When adjusting the tendencies from microphysics process rates,
2524 : ! conserve:
2525 : !
2526 : ! rv_mc_tend(k) + rc_mc_tend(k) + rr_mc_tend(k)
2527 : ! + ri_mc_tend(k) + rs_mc_tend(k) = 0; for all k from top_lev to pver.
2528 : !
2529 : ! 2) When adjusting the tendencies from microphysics process rates, adjust
2530 : ! dry static energy appropriately. The change in dry static energy
2531 : ! is necessary because of phase changes. This "puts back" the extra dry
2532 : ! static energy that was "taken out" when an excessive phase-changing
2533 : ! process rate was produced by microphysics.
2534 :
2535 : ! Loop over all columns, performing any tendency adjustments one column
2536 : ! at a time.
2537 0 : do icol = 1, ncol
2538 :
2539 : ! Loop over all vertical levels, performing any microphysics process
2540 : ! tendency adjustments one level at a time.
2541 0 : do k = top_lev, pver
2542 :
2543 : ! Find which hydrometeors have positive microphysics process
2544 : ! tendencies at this level.
2545 0 : if ( rv_mc_tend(icol,k) >= 0.0_r8 ) then
2546 : l_pos_rv_mc_tend = .true.
2547 : else
2548 0 : l_pos_rv_mc_tend = .false.
2549 : endif
2550 0 : if ( rc_mc_tend(icol,k) >= 0.0_r8 ) then
2551 : l_pos_rc_mc_tend = .true.
2552 : else
2553 0 : l_pos_rc_mc_tend = .false.
2554 : endif
2555 0 : if ( ixrain > 0 ) then
2556 0 : if ( rr_mc_tend(icol,k) >= 0.0_r8 ) then
2557 : l_pos_rr_mc_tend = .true.
2558 : else
2559 0 : l_pos_rr_mc_tend = .false.
2560 : endif
2561 : endif
2562 0 : if ( ri_mc_tend(icol,k) >= 0.0_r8 ) then
2563 : l_pos_ri_mc_tend = .true.
2564 : else
2565 0 : l_pos_ri_mc_tend = .false.
2566 : endif
2567 0 : if ( ixsnow > 0 ) then
2568 0 : if ( rs_mc_tend(icol,k) >= 0.0_r8 ) then
2569 : l_pos_rs_mc_tend = .true.
2570 : else
2571 0 : l_pos_rs_mc_tend = .false.
2572 : endif
2573 : endif
2574 :
2575 : !!! Check for holes in water vapor mixing ratio
2576 0 : if ( .not. l_pos_rv_mc_tend ) then
2577 :
2578 : ! Calculate the water vapor mixing ratio as it would be with the
2579 : ! current microphysics process tendency.
2580 0 : rv_curr = rv_start(icol,k) + rv_mc_tend(icol,k) * dt
2581 :
2582 0 : if ( rv_curr < qmin(1) ) then
2583 :
2584 : ! Microphysics processes are causing a hole in water vapor
2585 : ! mixing ratio.
2586 :
2587 : ! Calculate the maximum allowable magnitude of (negative) water
2588 : ! vapor microphysics process tendency.
2589 0 : mc_tend_max_mag = ( qmin(1) - rv_start(icol,k) ) / dt
2590 :
2591 : ! Calculate the amount of the correction that needs to be made
2592 : ! to the water vapor mixing ratio microphysics process
2593 : ! tendency. This number is positive.
2594 0 : mc_tend_correction = mc_tend_max_mag - rv_mc_tend(icol,k)
2595 :
2596 : ! Calculate the total amount of positive microphysics process
2597 : ! tendencies for all hydrometeor mixing ratios.
2598 0 : total_mc_positive = 0.0_r8
2599 0 : if ( l_pos_rc_mc_tend ) then
2600 0 : total_mc_positive = total_mc_positive + rc_mc_tend(icol,k)
2601 : endif
2602 0 : if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
2603 0 : total_mc_positive = total_mc_positive + rr_mc_tend(icol,k)
2604 : endif
2605 0 : if ( l_pos_ri_mc_tend ) then
2606 0 : total_mc_positive = total_mc_positive + ri_mc_tend(icol,k)
2607 : endif
2608 0 : if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
2609 0 : total_mc_positive = total_mc_positive + rs_mc_tend(icol,k)
2610 : endif
2611 :
2612 : ! Calculate the correction ratio.
2613 : ! In principle, this should never be greater than 1 outside of
2614 : ! numerical round-off errors. This is limited at 1 to be safe.
2615 : mc_correction_ratio &
2616 : = min( mc_tend_correction &
2617 0 : / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
2618 :
2619 : ! Adjust (decrease) the tendencies of all positive hydrometeor
2620 : ! mixing ratio tendencies to balance the adjustment (increase)
2621 : ! to the excessively negative water vapor mixing ratio.
2622 : ! Transfer dry static energy appropriately (in response to the
2623 : ! excessive depletion of water vapor).
2624 0 : if ( l_pos_rc_mc_tend ) then
2625 : ! Changing cloud water to water vapor cools and reduces
2626 : ! dry static energy.
2627 : stend(icol,k) &
2628 : = stend(icol,k) &
2629 0 : - latvap * mc_correction_ratio * rc_mc_tend(icol,k)
2630 : ! Update cloud water mixing ratio microphysics tendency.
2631 : rc_mc_tend(icol,k) &
2632 0 : = rc_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2633 : endif
2634 0 : if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
2635 : ! Changing rain water to water vapor cools and reduces
2636 : ! dry static energy.
2637 : stend(icol,k) &
2638 : = stend(icol,k) &
2639 0 : - latvap * mc_correction_ratio * rr_mc_tend(icol,k)
2640 : ! Update rain water mixing ratio microphysics tendency.
2641 : rr_mc_tend(icol,k) &
2642 0 : = rr_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2643 : endif
2644 0 : if ( l_pos_ri_mc_tend ) then
2645 : ! Changing cloud ice to water vapor cools and reduces
2646 : ! dry static energy.
2647 : stend(icol,k) &
2648 : = stend(icol,k) &
2649 : - ( latvap + latice ) &
2650 0 : * mc_correction_ratio * ri_mc_tend(icol,k)
2651 : ! Update cloud ice mixing ratio microphysics tendency.
2652 : ri_mc_tend(icol,k) &
2653 0 : = ri_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2654 : endif
2655 0 : if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
2656 : ! Changing snow to water vapor cools and reduces dry
2657 : ! static energy.
2658 : stend(icol,k) &
2659 : = stend(icol,k) &
2660 : - ( latvap + latice ) &
2661 0 : * mc_correction_ratio * rs_mc_tend(icol,k)
2662 : ! Update snow mixing ratio microphysics tendency.
2663 : rs_mc_tend(icol,k) &
2664 0 : = rs_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2665 : endif
2666 :
2667 : ! Calculate the new water vapor mixing ratio microphysics
2668 : ! process tendency. This should be equal to the maximum
2669 : ! magnitude (negative) amount allowed, mc_tend_max_mag.
2670 : rv_mc_tend(icol,k) &
2671 : = rv_mc_tend(icol,k) &
2672 0 : + mc_correction_ratio * total_mc_positive
2673 :
2674 : endif ! rv_curr < qmin(1)
2675 :
2676 : endif ! .not. l_pos_rv_mc_tend
2677 :
2678 : !!! Check for holes in cloud water mixing ratio
2679 0 : if ( .not. l_pos_rc_mc_tend ) then
2680 :
2681 : ! Calculate the cloud water mixing ratio as it would be with the
2682 : ! current microphysics process tendency.
2683 0 : rc_curr = rc_start(icol,k) + rc_mc_tend(icol,k) * dt
2684 :
2685 0 : if ( rc_curr < qmin(ixcldliq) ) then
2686 :
2687 : ! Microphysics processes are causing a hole in cloud water
2688 : ! mixing ratio.
2689 :
2690 : ! Calculate the maximum allowable magnitude of (negative) cloud
2691 : ! water microphysics process tendency.
2692 0 : mc_tend_max_mag = ( qmin(ixcldliq) - rc_start(icol,k) ) / dt
2693 :
2694 : ! Calculate the amount of the correction that needs to be made
2695 : ! to the cloud water mixing ratio microphysics process
2696 : ! tendency. This number is positive.
2697 0 : mc_tend_correction = mc_tend_max_mag - rc_mc_tend(icol,k)
2698 :
2699 : ! Calculate the total amount of positive microphysics process
2700 : ! tendencies for all hydrometeor mixing ratios.
2701 0 : total_mc_positive = 0.0_r8
2702 0 : if ( l_pos_rv_mc_tend ) then
2703 0 : total_mc_positive = total_mc_positive + rv_mc_tend(icol,k)
2704 : endif
2705 0 : if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
2706 0 : total_mc_positive = total_mc_positive + rr_mc_tend(icol,k)
2707 : endif
2708 0 : if ( l_pos_ri_mc_tend ) then
2709 0 : total_mc_positive = total_mc_positive + ri_mc_tend(icol,k)
2710 : endif
2711 0 : if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
2712 0 : total_mc_positive = total_mc_positive + rs_mc_tend(icol,k)
2713 : endif
2714 :
2715 : ! Calculate the correction ratio.
2716 : ! In principle, this should never be greater than 1 outside of
2717 : ! numerical round-off errors. This is limited at 1 to be safe.
2718 : mc_correction_ratio &
2719 : = min( mc_tend_correction &
2720 0 : / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
2721 :
2722 : ! Adjust (decrease) the tendencies of all positive hydrometeor
2723 : ! mixing ratio tendencies to balance the adjustment (increase)
2724 : ! to the excessively negative cloud water mixing ratio.
2725 : ! Transfer dry static energy appropriately (in response to the
2726 : ! excessive depletion of cloud water).
2727 0 : if ( l_pos_rv_mc_tend ) then
2728 : ! Changing water vapor to cloud water heats and increases
2729 : ! dry static energy.
2730 : stend(icol,k) &
2731 : = stend(icol,k) &
2732 0 : + latvap * mc_correction_ratio * rv_mc_tend(icol,k)
2733 : ! Update water vapor mixing ratio microphysics tendency.
2734 : rv_mc_tend(icol,k) &
2735 0 : = rv_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2736 : endif
2737 0 : if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
2738 : ! Changing rain water to cloud water does not change
2739 : ! dry static energy.
2740 : ! Update rain water mixing ratio microphysics tendency.
2741 : rr_mc_tend(icol,k) &
2742 0 : = rr_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2743 : endif
2744 0 : if ( l_pos_ri_mc_tend ) then
2745 : ! Changing cloud ice to cloud water cools and reduces
2746 : ! dry static energy.
2747 : stend(icol,k) &
2748 : = stend(icol,k) &
2749 0 : - latice * mc_correction_ratio * ri_mc_tend(icol,k)
2750 : ! Update cloud ice mixing ratio microphysics tendency.
2751 : ri_mc_tend(icol,k) &
2752 0 : = ri_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2753 : endif
2754 0 : if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
2755 : ! Changing snow to cloud water cools and reduces dry
2756 : ! static energy.
2757 : stend(icol,k) &
2758 : = stend(icol,k) &
2759 0 : - latice * mc_correction_ratio * rs_mc_tend(icol,k)
2760 : ! Update snow mixing ratio microphysics tendency.
2761 : rs_mc_tend(icol,k) &
2762 0 : = rs_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2763 : endif
2764 :
2765 : ! Calculate the new cloud water mixing ratio microphysics
2766 : ! process tendency. This should be equal to the maximum
2767 : ! magnitude (negative) amount allowed, mc_tend_max_mag.
2768 : rc_mc_tend(icol,k) &
2769 : = rc_mc_tend(icol,k) &
2770 0 : + mc_correction_ratio * total_mc_positive
2771 :
2772 : endif ! rc_curr < qmin(ixcldliq)
2773 :
2774 : endif ! .not. l_pos_rc_mc_tend
2775 :
2776 : !!! Check for holes in rain water mixing ratio
2777 0 : if ( ixrain > 0 .and. ( .not. l_pos_rr_mc_tend ) ) then
2778 :
2779 : ! Calculate the rain water mixing ratio as it would be with the
2780 : ! current microphysics process tendency.
2781 0 : rr_curr = rr_start(icol,k) + rr_mc_tend(icol,k) * dt
2782 :
2783 0 : if ( rr_curr < qmin(ixrain) ) then
2784 :
2785 : ! Microphysics processes are causing a hole in rain water
2786 : ! mixing ratio.
2787 :
2788 : ! Calculate the maximum allowable magnitude of (negative) rain
2789 : ! water microphysics process tendency.
2790 0 : mc_tend_max_mag = ( qmin(ixrain) - rr_start(icol,k) ) / dt
2791 :
2792 : ! Calculate the amount of the correction that needs to be made
2793 : ! to the rain water mixing ratio microphysics process
2794 : ! tendency. This number is positive.
2795 0 : mc_tend_correction = mc_tend_max_mag - rr_mc_tend(icol,k)
2796 :
2797 : ! Calculate the total amount of positive microphysics process
2798 : ! tendencies for all hydrometeor mixing ratios.
2799 0 : total_mc_positive = 0.0_r8
2800 0 : if ( l_pos_rv_mc_tend ) then
2801 0 : total_mc_positive = total_mc_positive + rv_mc_tend(icol,k)
2802 : endif
2803 0 : if ( l_pos_rc_mc_tend ) then
2804 0 : total_mc_positive = total_mc_positive + rc_mc_tend(icol,k)
2805 : endif
2806 0 : if ( l_pos_ri_mc_tend ) then
2807 0 : total_mc_positive = total_mc_positive + ri_mc_tend(icol,k)
2808 : endif
2809 0 : if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
2810 0 : total_mc_positive = total_mc_positive + rs_mc_tend(icol,k)
2811 : endif
2812 :
2813 : ! Calculate the correction ratio.
2814 : ! In principle, this should never be greater than 1 outside of
2815 : ! numerical round-off errors. This is limited at 1 to be safe.
2816 : mc_correction_ratio &
2817 : = min( mc_tend_correction &
2818 0 : / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
2819 :
2820 : ! Adjust (decrease) the tendencies of all positive hydrometeor
2821 : ! mixing ratio tendencies to balance the adjustment (increase)
2822 : ! to the excessively negative rain water mixing ratio.
2823 : ! Transfer dry static energy appropriately (in response to the
2824 : ! excessive depletion of rain water).
2825 0 : if ( l_pos_rv_mc_tend ) then
2826 : ! Changing water vapor to rain water heats and increases
2827 : ! dry static energy.
2828 : stend(icol,k) &
2829 : = stend(icol,k) &
2830 0 : + latvap * mc_correction_ratio * rv_mc_tend(icol,k)
2831 : ! Update water vapor mixing ratio microphysics tendency.
2832 : rv_mc_tend(icol,k) &
2833 0 : = rv_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2834 : endif
2835 0 : if ( l_pos_rc_mc_tend ) then
2836 : ! Changing cloud water to rain water does not change
2837 : ! dry static energy.
2838 : ! Update cloud water mixing ratio microphysics tendency.
2839 : rc_mc_tend(icol,k) &
2840 0 : = rc_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2841 : endif
2842 0 : if ( l_pos_ri_mc_tend ) then
2843 : ! Changing cloud ice to rain water cools and reduces
2844 : ! dry static energy.
2845 : stend(icol,k) &
2846 : = stend(icol,k) &
2847 0 : - latice * mc_correction_ratio * ri_mc_tend(icol,k)
2848 : ! Update cloud ice mixing ratio microphysics tendency.
2849 : ri_mc_tend(icol,k) &
2850 0 : = ri_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2851 : endif
2852 0 : if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
2853 : ! Changing snow to rain water cools and reduces dry
2854 : ! static energy.
2855 : stend(icol,k) &
2856 : = stend(icol,k) &
2857 0 : - latice * mc_correction_ratio * rs_mc_tend(icol,k)
2858 : ! Update snow mixing ratio microphysics tendency.
2859 : rs_mc_tend(icol,k) &
2860 0 : = rs_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2861 : endif
2862 :
2863 : ! Calculate the new rain water mixing ratio microphysics
2864 : ! process tendency. This should be equal to the maximum
2865 : ! magnitude (negative) amount allowed, mc_tend_max_mag.
2866 : rr_mc_tend(icol,k) &
2867 : = rr_mc_tend(icol,k) &
2868 0 : + mc_correction_ratio * total_mc_positive
2869 :
2870 : endif ! rr_curr < qmin(ixrain)
2871 :
2872 : endif ! ixrain > 0 .and. ( .not. l_pos_rr_mc_tend )
2873 :
2874 : !!! Check for holes in cloud ice mixing ratio
2875 0 : if ( .not. l_pos_ri_mc_tend ) then
2876 :
2877 : ! Calculate the cloud ice mixing ratio as it would be with the
2878 : ! current microphysics process tendency.
2879 0 : ri_curr = ri_start(icol,k) + ri_mc_tend(icol,k) * dt
2880 :
2881 0 : if ( ri_curr < qmin(ixcldice) ) then
2882 :
2883 : ! Microphysics processes are causing a hole in cloud ice
2884 : ! mixing ratio.
2885 :
2886 : ! Calculate the maximum allowable magnitude of (negative) cloud
2887 : ! ice microphysics process tendency.
2888 0 : mc_tend_max_mag = ( qmin(ixcldice) - ri_start(icol,k) ) / dt
2889 :
2890 : ! Calculate the amount of the correction that needs to be made
2891 : ! to the cloud ice mixing ratio microphysics process
2892 : ! tendency. This number is positive.
2893 0 : mc_tend_correction = mc_tend_max_mag - ri_mc_tend(icol,k)
2894 :
2895 : ! Calculate the total amount of positive microphysics process
2896 : ! tendencies for all hydrometeor mixing ratios.
2897 0 : total_mc_positive = 0.0_r8
2898 0 : if ( l_pos_rv_mc_tend ) then
2899 0 : total_mc_positive = total_mc_positive + rv_mc_tend(icol,k)
2900 : endif
2901 0 : if ( l_pos_rc_mc_tend ) then
2902 0 : total_mc_positive = total_mc_positive + rc_mc_tend(icol,k)
2903 : endif
2904 0 : if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
2905 0 : total_mc_positive = total_mc_positive + rr_mc_tend(icol,k)
2906 : endif
2907 0 : if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
2908 0 : total_mc_positive = total_mc_positive + rs_mc_tend(icol,k)
2909 : endif
2910 :
2911 : ! Calculate the correction ratio.
2912 : ! In principle, this should never be greater than 1 outside of
2913 : ! numerical round-off errors. This is limited at 1 to be safe.
2914 : mc_correction_ratio &
2915 : = min( mc_tend_correction &
2916 0 : / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
2917 :
2918 : ! Adjust (decrease) the tendencies of all positive hydrometeor
2919 : ! mixing ratio tendencies to balance the adjustment (increase)
2920 : ! to the excessively negative cloud ice mixing ratio.
2921 : ! Transfer dry static energy appropriately (in response to the
2922 : ! excessive depletion of cloud ice).
2923 0 : if ( l_pos_rv_mc_tend ) then
2924 : ! Changing water vapor to cloud ice heats and increases
2925 : ! dry static energy.
2926 : stend(icol,k) &
2927 : = stend(icol,k) &
2928 : + ( latvap + latice ) &
2929 0 : * mc_correction_ratio * rv_mc_tend(icol,k)
2930 : ! Update water vapor mixing ratio microphysics tendency.
2931 : rv_mc_tend(icol,k) &
2932 0 : = rv_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2933 : endif
2934 0 : if ( l_pos_rc_mc_tend ) then
2935 : ! Changing cloud water to cloud ice heats and increases
2936 : ! dry static energy.
2937 : stend(icol,k) &
2938 : = stend(icol,k) &
2939 0 : + latice * mc_correction_ratio * rc_mc_tend(icol,k)
2940 : ! Update cloud water mixing ratio microphysics tendency.
2941 : rc_mc_tend(icol,k) &
2942 0 : = rc_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2943 : endif
2944 0 : if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
2945 : ! Changing rain water to cloud ice heats and increases
2946 : ! dry static energy.
2947 : stend(icol,k) &
2948 : = stend(icol,k) &
2949 0 : + latice * mc_correction_ratio * rr_mc_tend(icol,k)
2950 : ! Update rain water mixing ratio microphysics tendency.
2951 : rr_mc_tend(icol,k) &
2952 0 : = rr_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2953 : endif
2954 0 : if ( ixsnow > 0 .and. l_pos_rs_mc_tend ) then
2955 : ! Changing snow to cloud ice does not change dry static
2956 : ! energy.
2957 : ! Update snow mixing ratio microphysics tendency.
2958 : rs_mc_tend(icol,k) &
2959 0 : = rs_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
2960 : endif
2961 :
2962 : ! Calculate the new cloud ice mixing ratio microphysics
2963 : ! process tendency. This should be equal to the maximum
2964 : ! magnitude (negative) amount allowed, mc_tend_max_mag.
2965 : ri_mc_tend(icol,k) &
2966 : = ri_mc_tend(icol,k) &
2967 0 : + mc_correction_ratio * total_mc_positive
2968 :
2969 : endif ! ri_curr < qmin(ixcldice)
2970 :
2971 : endif ! .not. l_pos_ri_mc_tend
2972 :
2973 : !!! Check for holes in snow mixing ratio
2974 0 : if ( ixsnow > 0 .and. ( .not. l_pos_rs_mc_tend ) ) then
2975 :
2976 : ! Calculate the snow mixing ratio as it would be with the
2977 : ! current microphysics process tendency.
2978 0 : rs_curr = rs_start(icol,k) + rs_mc_tend(icol,k) * dt
2979 :
2980 0 : if ( rs_curr < qmin(ixsnow) ) then
2981 :
2982 : ! Microphysics processes are causing a hole in snow mixing
2983 : ! ratio.
2984 :
2985 : ! Calculate the maximum allowable magnitude of (negative) snow
2986 : ! microphysics process tendency.
2987 0 : mc_tend_max_mag = ( qmin(ixsnow) - rs_start(icol,k) ) / dt
2988 :
2989 : ! Calculate the amount of the correction that needs to be made
2990 : ! to the snow mixing ratio microphysics process tendency.
2991 : ! This number is positive.
2992 0 : mc_tend_correction = mc_tend_max_mag - rs_mc_tend(icol,k)
2993 :
2994 : ! Calculate the total amount of positive microphysics process
2995 : ! tendencies for all hydrometeor mixing ratios.
2996 0 : total_mc_positive = 0.0_r8
2997 0 : if ( l_pos_rv_mc_tend ) then
2998 0 : total_mc_positive = total_mc_positive + rv_mc_tend(icol,k)
2999 : endif
3000 0 : if ( l_pos_rc_mc_tend ) then
3001 0 : total_mc_positive = total_mc_positive + rc_mc_tend(icol,k)
3002 : endif
3003 0 : if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
3004 0 : total_mc_positive = total_mc_positive + rr_mc_tend(icol,k)
3005 : endif
3006 0 : if ( l_pos_ri_mc_tend ) then
3007 0 : total_mc_positive = total_mc_positive + ri_mc_tend(icol,k)
3008 : endif
3009 :
3010 : ! Calculate the correction ratio.
3011 : ! In principle, this should never be greater than 1 outside of
3012 : ! numerical round-off errors. This is limited at 1 to be safe.
3013 : mc_correction_ratio &
3014 : = min( mc_tend_correction &
3015 0 : / max( total_mc_positive, 1.0e-30_r8 ), 1.0_r8 )
3016 :
3017 : ! Adjust (decrease) the tendencies of all positive hydrometeor
3018 : ! mixing ratio tendencies to balance the adjustment (increase)
3019 : ! to the excessively negative snow mixing ratio.
3020 : ! Transfer dry static energy appropriately (in response to the
3021 : ! excessive depletion of snow).
3022 0 : if ( l_pos_rv_mc_tend ) then
3023 : ! Changing water vapor to snow heats and increases dry
3024 : ! static energy.
3025 : stend(icol,k) &
3026 : = stend(icol,k) &
3027 : + ( latvap + latice ) &
3028 0 : * mc_correction_ratio * rv_mc_tend(icol,k)
3029 : ! Update water vapor mixing ratio microphysics tendency.
3030 : rv_mc_tend(icol,k) &
3031 0 : = rv_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
3032 : endif
3033 0 : if ( l_pos_rc_mc_tend ) then
3034 : ! Changing cloud water to snow heats and increases dry
3035 : ! static energy.
3036 : stend(icol,k) &
3037 : = stend(icol,k) &
3038 0 : + latice * mc_correction_ratio * rc_mc_tend(icol,k)
3039 : ! Update cloud water mixing ratio microphysics tendency.
3040 : rc_mc_tend(icol,k) &
3041 0 : = rc_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
3042 : endif
3043 0 : if ( ixrain > 0 .and. l_pos_rr_mc_tend ) then
3044 : ! Changing rain water to snow heats and increases dry
3045 : ! static energy.
3046 : stend(icol,k) &
3047 : = stend(icol,k) &
3048 0 : + latice * mc_correction_ratio * rr_mc_tend(icol,k)
3049 : ! Update rain water mixing ratio microphysics tendency.
3050 : rr_mc_tend(icol,k) &
3051 0 : = rr_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
3052 : endif
3053 0 : if ( l_pos_ri_mc_tend ) then
3054 : ! Changing cloud ice to snow does not change dry static
3055 : ! energy.
3056 : ! Update cloud ice mixing ratio microphysics tendency.
3057 : ri_mc_tend(icol,k) &
3058 0 : = ri_mc_tend(icol,k) * ( 1.0_r8 - mc_correction_ratio )
3059 : endif
3060 :
3061 : ! Calculate the new snow mixing ratio microphysics process
3062 : ! tendency. This should be equal to the maximum magnitude
3063 : ! (negative) amount allowed, mc_tend_max_mag.
3064 : rs_mc_tend(icol,k) &
3065 : = rs_mc_tend(icol,k) &
3066 0 : + mc_correction_ratio * total_mc_positive
3067 :
3068 : endif ! rs_curr < qmin(ixsnow)
3069 :
3070 : endif ! ixsnow > 0 .and. ( .not. l_pos_rs_mc_tend )
3071 :
3072 : end do ! k = top_lev, pver
3073 :
3074 : end do ! icol = 1, ncol
3075 :
3076 : ! Calculate the new overall tendencies by adding the sedimentation
3077 : ! tendencies back onto the new microphysics process tendencies.
3078 : ! For cloud water and cloud ice, the sedimentation tendencies that are
3079 : ! added back on are the true sedimentation tendencies. For cloud water,
3080 : ! this is the sum of rc_sed_tend and rc_sed_evap, and for cloud ice, this
3081 : ! is the sum of ri_sed_tend and ri_sed_subl.
3082 0 : rv_tend(:ncol,:) = rv_mc_tend(:ncol,:)
3083 0 : rc_tend(:ncol,:) = rc_mc_tend(:ncol,:) + ( rc_sed_tend(:ncol,:) + rc_sed_evap(:ncol,:) )
3084 0 : if ( ixrain > 0 ) then
3085 0 : rr_tend(:ncol,:) = rr_mc_tend(:ncol,:) + rr_sed_tend(:ncol,:)
3086 : endif
3087 0 : ri_tend(:ncol,:) = ri_mc_tend(:ncol,:) + ( ri_sed_tend(:ncol,:) + ri_sed_subl(:ncol,:) )
3088 0 : if ( ixsnow > 0 ) then
3089 0 : rs_tend(:ncol,:) = rs_mc_tend(:ncol,:) + rs_sed_tend(:ncol,:)
3090 : endif
3091 :
3092 : ! Now that the original sedimentation tendency has been added to the
3093 : ! new microphysics process tendency, the new total microphysics tendency
3094 : ! can still cause holes to form. After the microphysics process rates were
3095 : ! adjusted, the values of the hydrometeor fields were greater than or equal
3096 : ! to 0 at all grid levels, which means their vertical integrals were also
3097 : ! greater than or equal to 0. Sedimentation by itself has a vertical
3098 : ! integral of 0 (including the amount that sedimented to the surface).
3099 : ! This means that after the microphysics process rates have been adjusted
3100 : ! and sedimentation has been added back on, the resulting hydrometeor
3101 : ! fields all still have vertical integrals that are greater than or equal
3102 : ! to 0. Holes that develop at any particular grid level can be filled.
3103 : ! These holes can be filled conservatively using the sedimentation hole
3104 : ! filler.
3105 : if ( l_sed_hole_fill ) then
3106 :
3107 : ! This section makes use of the following principle:
3108 : !
3109 : ! 3) When adjusting the hydrometeor tendency from sedimentation of a
3110 : ! liquid hydrometeor (cloud water or rain water), conserve:
3111 : !
3112 : ! SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) )
3113 : ! * dt * pdel(k) / g
3114 : ! + precl * dt * 1000 = 0.
3115 :
3116 : ! Call the sedimentation hole filler for rain water mixing ratio.
3117 : ! This can update rr_tend and precl.
3118 0 : if ( ixrain > 0 ) then
3119 : call fill_holes_sedimentation( dt, ncol, rr_start, state%pdel, &
3120 0 : umr, state%zi, qmin(ixrain), &
3121 0 : rr_tend, precl )
3122 : endif ! ixrain > 0
3123 :
3124 : ! Call the sedimentation hole filler for cloud water mixing ratio.
3125 : ! This can update rc_tend and precl.
3126 : call fill_holes_sedimentation( dt, ncol, rc_start, state%pdel, &
3127 0 : vtrmc, state%zi, qmin(ixcldliq), &
3128 0 : rc_tend, precl )
3129 :
3130 : ! Occasionally, a situation can occur where filling a hole in rain can
3131 : ! deplete all the surface liquid-phase precipitation (precl), resulting
3132 : ! in not enough water mass in the vertical profile of cloud water to
3133 : ! fill a hole in cloud water. When this happens, there must be liquid
3134 : ! water found in the vertical profile of rain, so pull the water from
3135 : ! rain to fill any remaining holes in cloud water.
3136 0 : if ( ixrain > 0 ) then
3137 : call fill_holes_same_phase_vert( dt, ncol, rc_start, rr_start, &
3138 0 : state%pdel, qmin(ixcldliq), &
3139 0 : qmin(ixrain), &
3140 0 : rc_tend, rr_tend )
3141 : endif ! ixrain > 0
3142 :
3143 : ! This section makes use of the following principle:
3144 : !
3145 : ! 4) When adjusting the hydrometeor tendency from sedimentation of a
3146 : ! frozen hydrometeor (cloud ice or snow), conserve:
3147 : !
3148 : ! SUM(k=top_lev:pver) ( ri_sed_tend(k) + rs_sed_tend(k) )
3149 : ! * dt * pdel(k) / g
3150 : ! + preci * dt * 1000 = 0.
3151 :
3152 : ! Call the sedimentation hole filler for snow mixing ratio.
3153 : ! This can update rs_tend and preci.
3154 0 : if ( ixsnow > 0 ) then
3155 : call fill_holes_sedimentation( dt, ncol, rs_start, state%pdel, &
3156 0 : ums, state%zi, qmin(ixsnow), &
3157 0 : rs_tend, preci )
3158 : endif ! ixsnow > 0
3159 :
3160 : ! Call the sedimentation hole filler for cloud ice mixing ratio.
3161 : ! This can update ri_tend and preci.
3162 : call fill_holes_sedimentation( dt, ncol, ri_start, state%pdel, &
3163 0 : vtrmi, state%zi, qmin(ixcldice), &
3164 0 : ri_tend, preci )
3165 :
3166 : ! Occasionally, a situation can occur where filling a hole in snow can
3167 : ! deplete all the surface ice-phase precipitation (preci), resulting
3168 : ! in not enough water mass in the vertical profile of cloud ice to
3169 : ! fill a hole in cloud ice. When this happens, there must be ice-phase
3170 : ! water found in the vertical profile of snow, so pull the water from
3171 : ! snow to fill any remaining holes in cloud ice.
3172 0 : if ( ixsnow > 0 ) then
3173 : call fill_holes_same_phase_vert( dt, ncol, ri_start, rs_start, &
3174 0 : state%pdel, qmin(ixcldice), &
3175 0 : qmin(ixsnow), &
3176 0 : ri_tend, rs_tend )
3177 : endif ! ixsnow > 0
3178 :
3179 : ! Update the total precipitation rate (prect) from the updated liquid
3180 : ! precipitation rate (precl) and the updated frozen preciptation rate
3181 : ! (preci).
3182 0 : prect(:ncol) = precl(:ncol) + preci(:ncol)
3183 :
3184 : ! The MG code sets prec_str equal to prect (prec_pcw) and snow_str equal
3185 : ! to preci (snow_pcw). The prec_str and snow_str variables are used
3186 : ! in the calculations for energy and water conservation. Since prect
3187 : ! and preci are adjusted here, when necessary, prec_str and snow_str
3188 : ! also need to be adjusted.
3189 0 : prec_str(:ncol) = prect(:ncol)
3190 0 : snow_str(:ncol) = preci(:ncol)
3191 :
3192 : endif ! l_sed_hole_fill
3193 :
3194 : ! The updated total microphysics tendencies after hole filling have not
3195 : ! been used to update ptend yet, so record the budget terms for hole
3196 : ! filling first.
3197 0 : rv_hf_tend = rv_tend - ptend%q(:,:,1)
3198 0 : rc_hf_tend = rc_tend - ptend%q(:,:,ixcldliq)
3199 0 : if ( ixrain > 0 ) then
3200 0 : rr_hf_tend = rr_tend - ptend%q(:,:,ixrain)
3201 : endif ! ixrain > 0
3202 0 : ri_hf_tend = ri_tend - ptend%q(:,:,ixcldice)
3203 0 : if ( ixsnow > 0 ) then
3204 0 : rs_hf_tend = rs_tend - ptend%q(:,:,ixsnow)
3205 : endif ! ixsnow > 0
3206 :
3207 : ! The updated dry static energy tendency after hole filling has not been
3208 : ! used to update ptend yet, so record the budget term for hole filling
3209 : ! first.
3210 0 : s_hf_tend = stend - ptend%s
3211 :
3212 : ! Pack the current total tendencies for hydrometeor mixing ratio fields.
3213 0 : ptend%q(:,:,1) = rv_tend
3214 0 : ptend%q(:,:,ixcldliq) = rc_tend
3215 0 : if ( ixrain > 0 ) then
3216 0 : ptend%q(:,:,ixrain) = rr_tend
3217 : endif
3218 0 : ptend%q(:,:,ixcldice) = ri_tend
3219 0 : if ( ixsnow > 0 ) then
3220 0 : ptend%q(:,:,ixsnow) = rs_tend
3221 : endif
3222 :
3223 : ! Pack the current tendency for dry static energy.
3224 0 : ptend%s = stend
3225 :
3226 : ! Output stats for hole filling tendencies.
3227 0 : call outfld( 'QVHFTEN', rv_hf_tend, pcols, state%lchnk )
3228 0 : call outfld( 'QCHFTEN', rc_hf_tend, pcols, state%lchnk )
3229 0 : call outfld( 'QRHFTEN', rr_hf_tend, pcols, state%lchnk )
3230 0 : call outfld( 'QIHFTEN', ri_hf_tend, pcols, state%lchnk )
3231 0 : call outfld( 'QSHFTEN', rs_hf_tend, pcols, state%lchnk )
3232 0 : call outfld( 'THFTEN', s_hf_tend / cpair, pcols, state%lchnk )
3233 :
3234 : ! Perform total water and total energy conservation checks.
3235 : if ( l_check_conservation ) then
3236 :
3237 : ! Calculate total water in each grid column.
3238 : ! This calculation is the vertically-integrated grand total water
3239 : ! in each grid column updated for all microphysics and hole filling.
3240 : ! This includes the amount that precipitated to the surface.
3241 0 : do icol = 1, ncol
3242 0 : grand_total_water_column_finish(icol) = 0.0_r8
3243 0 : do k = top_lev, pver
3244 : grand_total_water_column_finish(icol) &
3245 : = grand_total_water_column_finish(icol) &
3246 0 : + ( state%q(icol,k,1) &
3247 0 : + ptend%q(icol,k,1) * dt &
3248 0 : + state%q(icol,k,ixcldliq) &
3249 0 : + ptend%q(icol,k,ixcldliq) * dt &
3250 0 : + state%q(icol,k,ixcldice) &
3251 0 : + ptend%q(icol,k,ixcldice) * dt ) &
3252 0 : * state%pdel(icol,k) * rga
3253 0 : if ( ixrain > 0 ) then
3254 : grand_total_water_column_finish(icol) &
3255 : = grand_total_water_column_finish(icol) &
3256 0 : + ( state%q(icol,k,ixrain) + ptend%q(icol,k,ixrain) * dt ) &
3257 0 : * state%pdel(icol,k) * rga
3258 : endif
3259 0 : if ( ixsnow > 0 ) then
3260 : grand_total_water_column_finish(icol) &
3261 : = grand_total_water_column_finish(icol) &
3262 0 : + ( state%q(icol,k,ixsnow) + ptend%q(icol,k,ixsnow) * dt ) &
3263 0 : * state%pdel(icol,k) * rga
3264 : endif
3265 : end do ! k = top_lev, pver
3266 : grand_total_water_column_finish(icol) &
3267 : = grand_total_water_column_finish(icol) &
3268 0 : + prect(icol) * dt * 1000.0_r8
3269 : end do ! icol = 1, ncol
3270 :
3271 : ! Calculate total energy in each column.
3272 : ! This calculation is the vertically-integrated total energy in each
3273 : ! grid column updated for all microphysics and hole filling. This
3274 : ! includes the amount that precipitated to the surface. Since, the
3275 : ! microphysics code does not directly change kinetic energy,
3276 : ! 0.5 * ( u^2 + v^2 ), it can be skipped as part of the energy
3277 : ! conservation check.
3278 0 : do icol = 1, ncol
3279 0 : total_energy_column_finish(icol) = 0.0_r8
3280 0 : do k = top_lev, pver
3281 : total_energy_column_finish(icol) &
3282 : = total_energy_column_finish(icol) &
3283 0 : + ( state%s(icol,k) + ptend%s(icol,k) * dt &
3284 : + ( latvap + latice ) &
3285 0 : * ( state%q(icol,k,1) + ptend%q(icol,k,1) * dt ) &
3286 0 : + latice * ( state%q(icol,k,ixcldliq) &
3287 0 : + ptend%q(icol,k,ixcldliq) * dt ) ) &
3288 0 : * state%pdel(icol,k) * rga
3289 0 : if ( ixrain > 0 ) then
3290 : total_energy_column_finish(icol) &
3291 : = total_energy_column_finish(icol) &
3292 0 : + latice * ( state%q(icol,k,ixrain) &
3293 0 : + ptend%q(icol,k,ixrain) * dt ) &
3294 0 : * state%pdel(icol,k) * rga
3295 : endif
3296 : end do ! k = top_lev, pver
3297 : total_energy_column_finish(icol) &
3298 : = total_energy_column_finish(icol) &
3299 0 : + latice * precl(icol) * dt * 1000.0_r8
3300 : end do ! icol = 1, ncol
3301 :
3302 : ! Calculate the total relative error in each grid column.
3303 0 : do icol = 1, ncol
3304 :
3305 0 : tot_water_rel_err(icol) &
3306 : = abs( ( grand_total_water_column_finish(icol) &
3307 : - grand_total_water_column_start(icol) ) ) &
3308 : / min( grand_total_water_column_finish(icol), &
3309 0 : grand_total_water_column_start(icol) )
3310 :
3311 : tot_energy_rel_err(icol) &
3312 : = abs( ( total_energy_column_finish(icol) &
3313 : - total_energy_column_start(icol) ) ) &
3314 : / min( total_energy_column_finish(icol), &
3315 0 : total_energy_column_start(icol) )
3316 :
3317 : end do ! icol = 1, ncol
3318 :
3319 : ! Print an error message if any total water relative error is found to
3320 : ! be greater than the threshold.
3321 0 : if ( any( tot_water_rel_err(:ncol) >= err_thresh ) ) then
3322 0 : write(iulog,*) "Water conservation error reported in hole filling"
3323 0 : do icol = 1, ncol
3324 0 : if ( tot_water_rel_err(icol) >= err_thresh ) then
3325 0 : write(iulog,*) "Column = ", icol, &
3326 0 : "Relative error = ", tot_water_rel_err(icol), &
3327 0 : "Column-integrated grand total water at start = ", &
3328 0 : grand_total_water_column_start(icol), &
3329 0 : "Column-integrated grand total water at finish = ", &
3330 0 : grand_total_water_column_finish(icol)
3331 : endif ! tot_water_rel_err(icol) >= err_thresh
3332 : end do ! icol = 1, ncol
3333 : endif ! any( tot_water_rel_err >= err_thresh )
3334 :
3335 : ! Print an error message if any total energy relative error is found to
3336 : ! be greater than the threshold.
3337 0 : if ( any( tot_energy_rel_err(:ncol) >= err_thresh ) ) then
3338 0 : write(iulog,*) "Energy conservation error reported in hole filling"
3339 0 : do icol = 1, ncol
3340 0 : if ( tot_energy_rel_err(icol) >= err_thresh ) then
3341 0 : write(iulog,*) "Column = ", icol, &
3342 0 : "Relative error = ", tot_energy_rel_err(icol), &
3343 0 : "Column-integrated total energy at start = ", &
3344 0 : total_energy_column_start(icol), &
3345 0 : "Column-integrated total energy at finish = ", &
3346 0 : total_energy_column_finish(icol)
3347 : endif ! tot_energy_rel_err(icol) >= err_thresh
3348 : end do ! icol = 1, ncol
3349 : endif ! any( tot_energy_rel_err >= err_thresh )
3350 :
3351 : endif ! l_check_conservation
3352 :
3353 :
3354 0 : return
3355 :
3356 0 : end subroutine subcol_SILHS_fill_holes_conserv
3357 :
3358 : !============================================================================
3359 0 : subroutine fill_holes_sedimentation( dt, ncol, hm_start, pdel, &
3360 : fallspeed_m_per_s, zi, qmin_hm, &
3361 : hm_tend, prec )
3362 :
3363 : ! Description:
3364 : ! After hydrometeor tendencies from microphysics processes were adjusted
3365 : ! so that holes don't form in a hydrometeor field from microphysics
3366 : ! processes, the sedimentation tendency was added back on to produce an
3367 : ! updated total microphysics tendency. The first-order "up-gravity"
3368 : ! sedimentation method that was originally used is positive definite.
3369 : ! However, since the microphysics process tendencies were altered so that
3370 : ! holes don't form, it is possible that adding the old sedimentation
3371 : ! tendencies back onto the new microphysics process tendencies could
3372 : ! produce new total microphysics tendencies that cause holes to form.
3373 : !
3374 : ! In this subroutine, holes in a hydrometeor field are checked for after
3375 : ! the updated microphysics tendency is applied. If any are found, they are
3376 : ! filled from positive hydrometeor mass found at grid levels below where
3377 : ! the hole is found. The levels that are used to fill are within range
3378 : ! based on fallspeed of the hydrometeor. If the level that contains the
3379 : ! hole is within proximity to the surface, then the water that sedimented
3380 : ! to the surface can be used to fill the hole, as well.
3381 : !
3382 : ! If there isn't enough total hydrometeor mass within the fall range to
3383 : ! fill the hole, then positive hydrometeor mass from levels below the fall
3384 : ! range is to be added to the total available mass to fill the hole. Mass
3385 : ! is added one level at a time until enough mass is found to fill the hole
3386 : ! or until the surface is reached and the surface precipitation is added to
3387 : ! the total available fill mass.
3388 : !
3389 : ! After this, if there still isn't enough available mass to fill the hole,
3390 : ! then positive hydrometeor mass is added from all levels above the hole to
3391 : ! the total mass that is available to fill the hole.
3392 :
3393 : !----------------------------------------------------------------------
3394 :
3395 0 : use ppgrid, only: &
3396 : pcols
3397 :
3398 : use ref_pres, only: &
3399 : top_lev => trop_cloud_top_lev
3400 :
3401 : implicit none
3402 :
3403 : ! Input Variables
3404 : real(r8), intent(in) :: dt ! Time step duration
3405 :
3406 : integer, intent(in) :: ncol ! Number of grid columns
3407 :
3408 : real(r8), dimension(pcols,pver), intent(in) :: &
3409 : hm_start, & ! Hydrometeor mixing ratio at start of microphysics [kg/kg]
3410 : pdel ! Pressure difference between grid levels [Pa]
3411 :
3412 : real(r8), dimension(pcols,pver), intent(in) :: &
3413 : fallspeed_m_per_s ! Hydrometeor mixing ratio fall speed [m/s]
3414 :
3415 : real(r8), dimension(pcols,pverp), intent(in) :: &
3416 : zi ! Height of momentum (interface) grid levels [m]
3417 :
3418 : real(r8), intent(in) :: &
3419 : qmin_hm ! Minimum threshold value of hydrometeor mixing ratio [kg/kg]
3420 :
3421 : ! Input/Output Variables
3422 : real(r8), dimension(pcols,pver), intent(inout) :: &
3423 : hm_tend ! Hydrometeor mixing ratio tendency [kg/kg/s]
3424 :
3425 : real(r8), dimension(pcols), intent(inout) :: &
3426 : prec ! Precipitation rate (surface) [m/s]
3427 :
3428 : ! Local Variables
3429 : real(r8), dimension(pver) :: &
3430 : hm_update, & ! Hydrometeor mixing ratio; start of sed. hole fill [kg/kg]
3431 : hm_curr ! Current value of hydrometeor mixing ratio [kg/kg]
3432 :
3433 : real(r8) :: &
3434 : total_hole, & ! Total mass of hole in hydrometeor [kg/m^2]
3435 : total_fill_mass, & ! Total mass available to fill hole [kg/m^2]
3436 : hole_fillmass_ratio, & ! Ratio: total_hole / total_fill_mass [-]
3437 : fallspeed_Pa_per_s, & ! Hydrometeor mixing ratio fall speed [Pa/s]
3438 : total_fall_Pa, & ! Pressure "distance" hydrometeor fell [Pa]
3439 : sum_pdel ! Sum of pdel over levels [Pa]
3440 :
3441 : logical, dimension(pver) :: &
3442 : l_pos_hm ! Flag for a hydrometeor having a positive (>= qmin_hm) value
3443 :
3444 : ! Flag for whether surface precipitation mass needs to be included in
3445 : ! the total_fill_mass for hole filling.
3446 : logical :: l_reached_surface
3447 :
3448 : ! Flag for whether hydrometeor mass from levels above the hole needs to be
3449 : ! included in the total_fill_mass for hole filling.
3450 : logical :: l_fill_from_above
3451 :
3452 : integer :: icol ! Grid column index
3453 :
3454 : integer :: k, idx ! Vertical grid level indices
3455 :
3456 : ! Index of the lowest vertical grid level that needs to be included in the
3457 : ! total_fill_mass for hole filling.
3458 : integer :: lowest_level_idx
3459 :
3460 :
3461 : ! Loop over all columns, performing any adjustments one column at a time.
3462 0 : do icol = 1, ncol
3463 :
3464 : ! Calculate the updated value of the hydrometeor field based on the
3465 : ! updated microphysics tendency. Since the original sedimentation
3466 : ! tendency has been added to the updated microphysics process tendency
3467 : ! to produce the updated total microphysics tendency (hm_tend), the
3468 : ! updated value of the hydrometeor field (hm_update) could be negative.
3469 0 : hm_update = hm_start(icol,:) + hm_tend(icol,:) * dt
3470 0 : hm_curr = hm_update
3471 :
3472 : ! Check for any holes in the vertical profile
3473 0 : if ( any( hm_curr(top_lev:pver) < qmin_hm ) ) then
3474 :
3475 : ! At least one hole is found in this hydrometeor species in this
3476 : ! grid column. The holes must be filled conservatively.
3477 :
3478 : ! Check which levels have values of the hydrometeor that are at or
3479 : ! above the minimum threshold value.
3480 0 : do k = top_lev, pver
3481 0 : if ( hm_curr(k) >= qmin_hm ) then
3482 0 : l_pos_hm(k) = .true.
3483 : else ! hm_curr < qmin_hm
3484 0 : l_pos_hm(k) = .false.
3485 : endif ! hm_curr >= qmin_hm
3486 : end do ! k = top_lev, pver
3487 :
3488 0 : do k = pver, top_lev, -1
3489 :
3490 0 : if ( .not. l_pos_hm(k) ) then
3491 :
3492 : ! A hole is found in the hydrometeor at this grid level.
3493 :
3494 : ! Calculate the total hydrometeor mass of the hole that needs
3495 : ! to be filled.
3496 : ! The value of the hydrometeor mixing ratio is negative, but
3497 : ! the value of total_hole is positive.
3498 0 : total_hole = ( qmin_hm - hm_curr(k) ) * pdel(icol,k) * rga
3499 :
3500 : ! Calculate the total hydrometeor mass available from below
3501 : ! to fill the hole.
3502 0 : if ( k == pver ) then
3503 :
3504 : ! A hole is found at the lowermost level.
3505 : ! The only place the hydrometeor could have sedimented
3506 : ! to is the surface, so fill from only the surface.
3507 0 : l_reached_surface = .true.
3508 :
3509 : ! Calculate the available amount of hydrometeor mass to
3510 : ! fill the hole.
3511 0 : total_fill_mass = prec(icol) * dt * 1000.0_r8
3512 :
3513 : else ! top_lev <= k < pver
3514 :
3515 : ! Calculate the hydrometeor fallspeed in Pa/s.
3516 : ! In MG2, the equation for this is given by:
3517 : !
3518 : ! fallspeed([Pa/s]) = g * rho * fallspeed([m/s]).
3519 : !
3520 : ! The value of rho is typically calculated from the
3521 : ! hydrostatic approximation:
3522 : !
3523 : ! rho = - ( 1 / g ) * dp/dz.
3524 : !
3525 : ! The equation for fallspeed in Pa/s becomes:
3526 : !
3527 : ! fallspeed([Pa/s]) = - dp/dz * fallspeed([m/s]).
3528 : fallspeed_Pa_per_s &
3529 : = fallspeed_m_per_s(icol,k) &
3530 0 : * pdel(icol,k) / ( zi(icol,k) - zi(icol,k+1) )
3531 :
3532 : ! Calculate the fall "distance" in Pa.
3533 0 : total_fall_Pa = fallspeed_Pa_per_s * dt
3534 :
3535 : ! Find the index of the vertical level that the hydrometeor
3536 : ! sedimented to in one timestep. It must sediment at least
3537 : ! one level.
3538 0 : sum_pdel = 0.0_r8
3539 0 : idx = k + 1
3540 0 : do
3541 : ! Update the total pressure difference between the
3542 : ! level of origin and the current level.
3543 0 : sum_pdel = sum_pdel + pdel(icol,idx)
3544 0 : if ( sum_pdel >= total_fall_Pa ) then
3545 : ! The total pressure difference between the level of
3546 : ! origin and the current level exceeds the total
3547 : ! hydrometeor fall "distance" (in Pa).
3548 : lowest_level_idx = idx
3549 : l_reached_surface = .false.
3550 : exit
3551 : else ! sum_pdel < total_fall_Pa
3552 : ! The total hydrometeor fall "distance" (in Pa)
3553 : ! exceeds the total pressure difference between the
3554 : ! level of origin and the current level.
3555 0 : if ( idx == pver ) then
3556 : ! The lowest level of the model has been reached.
3557 : ! The hydrometeor sedimented to the surface.
3558 : lowest_level_idx = pver
3559 : l_reached_surface = .true.
3560 : exit
3561 : else ! idx < pver
3562 : ! Increment idx and keep going.
3563 0 : idx = idx + 1
3564 : endif ! idx == pver
3565 : endif ! sum_pdel >= total_fall_Pa
3566 : end do
3567 :
3568 : ! Calculate the available amount of hydrometeor mass to
3569 : ! fill the hole.
3570 : total_fill_mass = 0.0_r8
3571 : if ( l_reached_surface ) then
3572 : ! The hydrometeor sedimented to the surface, so
3573 : ! automatically loop down to pver and include the
3574 : ! surface mass.
3575 0 : do idx = k+1, pver, 1
3576 0 : if ( l_pos_hm(idx) ) then
3577 : total_fill_mass &
3578 : = total_fill_mass &
3579 : + ( hm_curr(idx) - qmin_hm ) &
3580 0 : * pdel(icol,idx) * rga
3581 : endif ! l_pos_hm(idx)
3582 : end do ! idx = k+1, pver, 1
3583 : ! Contribution to total fill mass from the surface.
3584 : total_fill_mass &
3585 0 : = total_fill_mass + prec(icol) * dt * 1000.0_r8
3586 : else ! .not. l_reached_surface
3587 : ! The hydrometeor sedimented to lowest_level_idx.
3588 : idx = k + 1
3589 : do
3590 0 : if ( l_pos_hm(idx) ) then
3591 : total_fill_mass &
3592 : = total_fill_mass &
3593 : + ( hm_curr(idx) - qmin_hm ) &
3594 0 : * pdel(icol,idx) * rga
3595 : endif ! l_pos_hm(idx)
3596 0 : if ( idx >= lowest_level_idx ) then
3597 : ! Check if enough mass has been gathered in
3598 : ! total_fill_mass to fill the hole.
3599 0 : if ( total_fill_mass >= total_hole ) then
3600 : ! There has been enough total_fill_mass
3601 : ! gathered to completely fill the hole.
3602 : lowest_level_idx = idx
3603 : exit
3604 : else ! total_fill_mass < total_hole
3605 : ! Even though lowest_level_idx has been reached,
3606 : ! more total_fill_mass needs to be added in
3607 : ! order to completely fill the hole, so keep
3608 : ! going.
3609 0 : if ( idx == pver ) then
3610 : ! The lowest vertical level has already been
3611 : ! reached, so go to the surface.
3612 0 : lowest_level_idx = pver
3613 0 : l_reached_surface = .true.
3614 : ! Contribution to total fill mass from the
3615 : ! surface.
3616 : total_fill_mass &
3617 : = total_fill_mass &
3618 0 : + prec(icol) * dt * 1000.0_r8
3619 0 : exit
3620 : else ! idx < pver
3621 : ! Haven't reached pver yet, so increment
3622 : ! and keep going.
3623 0 : idx = idx + 1
3624 : endif ! idx == pver
3625 : endif ! total_fill_mass >= total_hole
3626 : else ! idx < lowest_level_idx
3627 : ! Haven't reached lowest_level_idx yet, so
3628 : ! increment and keep going.
3629 0 : idx = idx + 1
3630 : endif ! idx >= lowest_level_idx
3631 : end do
3632 : endif ! l_reached_surface
3633 :
3634 : endif ! k == pver
3635 :
3636 : ! If mass has been added all the way down to the surface and
3637 : ! there's still not enough mass to fill the hole, then fill the
3638 : ! hole pulling mass from above.
3639 0 : if ( total_fill_mass >= total_hole ) then
3640 : l_fill_from_above = .false.
3641 : else ! total_fill_mass < total_hole
3642 0 : l_fill_from_above = .true.
3643 0 : do idx = top_lev, k-1, 1
3644 0 : if ( l_pos_hm(idx) ) then
3645 : total_fill_mass &
3646 : = total_fill_mass &
3647 : + ( hm_curr(idx) - qmin_hm ) &
3648 0 : * pdel(icol,idx) * rga
3649 : endif ! l_pos_hm(idx)
3650 : end do ! idx = top_lev, k-1, 1
3651 : endif ! total_fill_mass >= total_hole
3652 :
3653 : ! Calculate the ratio of total hole to total fill mass. This
3654 : ! should not exceed 1 except as a result of numerical round-off
3655 : ! errors. Use thresholding to be safe.
3656 : hole_fillmass_ratio &
3657 : = min( total_hole / max( total_fill_mass, 1.0e-30_r8 ), &
3658 0 : 1.0_r8 )
3659 :
3660 0 : if ( k < pver ) then
3661 : ! Modify (reduce) the amount of the hydrometeor at levels
3662 : ! that were used to fill the hole.
3663 0 : do idx = k+1, lowest_level_idx
3664 0 : if ( l_pos_hm(idx) ) then
3665 : ! Since pdel at a grid level does not change and
3666 : ! gravit is constant, the only variable that needs to
3667 : ! be modified proportionately is hm_curr.
3668 : hm_curr(idx) &
3669 : = qmin_hm &
3670 : + ( hm_curr(idx) - qmin_hm ) &
3671 0 : * ( 1.0_r8 - hole_fillmass_ratio )
3672 : endif ! l_pos_hm(idx)
3673 : end do ! idx = k+1, lowest_level_idx
3674 : endif ! k < pver
3675 :
3676 0 : if ( l_reached_surface ) then
3677 : ! Modify (reduce) the amount of surface precipitation in
3678 : ! order to fill the hole. Since dt and 1000 are constants,
3679 : ! the only variable that needs to be modified
3680 : ! proportionately is prec.
3681 0 : prec(icol) = prec(icol) * ( 1.0_r8 - hole_fillmass_ratio )
3682 : endif ! l_reached_surface
3683 :
3684 0 : if ( l_fill_from_above ) then
3685 : ! Modify (reduce) the amount of the hydrometeor at levels
3686 : ! that were used to fill the hole.
3687 0 : do idx = top_lev, k-1
3688 0 : if ( l_pos_hm(idx) ) then
3689 : ! Since pdel at a grid level does not change and
3690 : ! gravit is constant, the only variable that needs to
3691 : ! be modified proportionately is hm_curr.
3692 : hm_curr(idx) &
3693 : = qmin_hm &
3694 : + ( hm_curr(idx) - qmin_hm ) &
3695 0 : * ( 1.0_r8 - hole_fillmass_ratio )
3696 : endif ! l_pos_hm(idx)
3697 : end do ! idx = top_lev, k-1
3698 : endif ! l_fill_from_above
3699 :
3700 : ! Update the value of the hydrometeor at the level where the
3701 : ! hole was found. Mathematically, as long as the available
3702 : ! mass was able to fill the entire hole, the new value of the
3703 : ! hydrometeor mixing ratio (hm_curr) should be qmin_hm.
3704 : hm_curr(k) &
3705 : = hm_curr(k) &
3706 : + hole_fillmass_ratio * total_fill_mass &
3707 0 : * gravit / pdel(icol,k)
3708 :
3709 : endif ! .not. l_pos_hm(k)
3710 :
3711 : end do ! k = pver, top_lev, -1
3712 :
3713 : endif ! any( hm_curr(top_lev:pver) < qmin_hm )
3714 :
3715 : ! Update the value of total microphysics tendency after hole filling.
3716 0 : hm_tend(icol,:) = hm_tend(icol,:) + ( hm_curr - hm_update ) / dt
3717 :
3718 : end do ! icol = 1, ncol
3719 :
3720 :
3721 0 : return
3722 :
3723 0 : end subroutine fill_holes_sedimentation
3724 :
3725 : !============================================================================
3726 0 : subroutine fill_holes_same_phase_vert( dt, ncol, hm_start, hm_start_filler, &
3727 : pdel, qmin_hm, qmin_hm_filler, &
3728 : hm_tend, hm_tend_filler )
3729 :
3730 : ! Description:
3731 : ! Fills remaining holes in a hydrometeor with mass from the the vertical
3732 : ! profile of another hydrometeor of the same phase. Remaining holes in
3733 : ! cloud water are filled with rain water and remaining holes in snow are
3734 : ! filled with cloud ice.
3735 : !
3736 : ! This subroutine, combined with subroutine fill_holes_sedimentation, fill
3737 : ! holes making use of the following principles:
3738 : !
3739 : ! 3) When adjusting the hydrometeor tendency from sedimentation of a
3740 : ! liquid hydrometeor (cloud water or rain water), conserve:
3741 : !
3742 : ! SUM(k=top_lev:pver) ( rc_sed_tend(k) + rr_sed_tend(k) )
3743 : ! * dt * pdel(k) / g
3744 : ! + precl * dt * 1000 = 0.
3745 : !
3746 : ! 4) When adjusting the hydrometeor tendency from sedimentation of a
3747 : ! frozen hydrometeor (cloud ice or snow), conserve:
3748 : !
3749 : ! SUM(k=top_lev:pver) ( ri_sed_tend(k) + rs_sed_tend(k) )
3750 : ! * dt * pdel(k) / g
3751 : ! + preci * dt * 1000 = 0.
3752 : !
3753 : ! These two equations (one for liquid-phase hydrometeors and one for
3754 : ! ice-phase hydrometeors) could be further split into one equation for
3755 : ! each hydrometeor if there was prec output for each hydrometeor. However,
3756 : ! there's only prec output for ice-phase precipitation rate and total
3757 : ! precipitation rate (liquid preciptation rate is total rate minus
3758 : ! ice-phase rate).
3759 : !
3760 : ! Since only liquid-phase precipitation rate (precl) and ice-phase
3761 : ! precipitation rate (preci) are available, and there are two hydrometeors
3762 : ! in each category, one hydrometeor from each category must fill before
3763 : ! the other hydrometeor from its category and get priority access to precl
3764 : ! or preci. Since a vast majority of liquid precipitation comes from rain
3765 : ! rather than sedimenting cloud water, rain is filled before cloud water
3766 : ! and gets priority access to precl. Likewise, since a vast majority of
3767 : ! frozen precipitation comes from snow rather than sedimenting cloud ice,
3768 : ! snow is filled before cloud ice and gets priority access to preci.
3769 : !
3770 : ! The order of sedimentation hole filling is as follows. First, a level
3771 : ! with a hole in it is identified. The fall distance for the hydrometeor
3772 : ! that originated at a level is calculated. Total mass to fill the hole is
3773 : ! calculated from all levels within the fall range that have positive
3774 : ! values of the hydrometeor. The amount that precipitated to the surface
3775 : ! is also included if the hydrometeor fell that far. If that isn't enough
3776 : ! mass to fill the hole, then levels that are lower in the profile are
3777 : ! included (if the hydrometeor has a positive value) until enough mass is
3778 : ! found to fill the hole or until the surface is reached. If there isn't
3779 : ! enough mass found in all levels below the hole, including the amount that
3780 : ! precipitated to the ground, to fill the hole, then the hydrometeor mass
3781 : ! from all levels above the hole (again, where a positive value of the
3782 : ! hydrometeor is found) are included in the total available mass to fill
3783 : ! the hole.
3784 : !
3785 : ! Occasionally, a situation can occur where both hydrometeors in a category
3786 : ! contributed to surface precipitation rate, and filling a hole in rain
3787 : ! (or snow) can deplete all the surface precl (or preci), resulting in not
3788 : ! enough water mass in the vertical profile (including the surface) of
3789 : ! cloud water (or cloud ice) to fill a hole in cloud water (or cloud ice).
3790 : ! When this happens, there must still be liquid water (or frozen water)
3791 : ! found in the vertical profile of rain (or snow), so pull the water from
3792 : ! rain (or snow) to fill any remaining holes in cloud water (or cloud ice).
3793 :
3794 : !----------------------------------------------------------------------
3795 :
3796 0 : use ppgrid, only: &
3797 : pcols
3798 :
3799 : use ref_pres, only: &
3800 : top_lev => trop_cloud_top_lev
3801 :
3802 : implicit none
3803 :
3804 : ! Input Variables
3805 : real(r8), intent(in) :: dt ! Time step duration
3806 :
3807 : integer, intent(in) :: ncol ! Number of grid columns
3808 :
3809 : real(r8), dimension(pcols,pver), intent(in) :: &
3810 : hm_start, & ! Hydrometeor mixing ratio (microphys start) [kg/kg]
3811 : hm_start_filler, & ! Filler hydromet mix ratio (microphys start) [kg/kg]
3812 : pdel ! Pressure difference between grid levels [Pa]
3813 :
3814 : real(r8), intent(in) :: &
3815 : qmin_hm, & ! Minimum threshold hydrometeor mixing ratio [kg/kg]
3816 : qmin_hm_filler ! Min threshold filler hydromet mixing ratio [kg/kg]
3817 :
3818 : ! Input/Output Variables
3819 : real(r8), dimension(pcols,pver), intent(inout) :: &
3820 : hm_tend, & ! Hydrometeor mixing ratio tendency [kg/kg/s]
3821 : hm_tend_filler ! Filler hydrometeor mixing ratio tendency [kg/kg/s]
3822 :
3823 : ! Local Variables
3824 : real(r8), dimension(pver) :: &
3825 : hm_update, & ! Hydrometeor mixing ratio; start [kg/kg]
3826 : hm_update_filler, & ! Filler Hydrometeor mixing ratio; start [kg/kg]
3827 : hm_curr, & ! Current hydrometeor mixing ratio [kg/kg]
3828 : hm_curr_filler ! Current filler hydrometeor mixing ratio [kg/kg]
3829 :
3830 : real(r8) :: &
3831 : total_hole, & ! Total mass of hole in hydrometeor [kg/m^2]
3832 : total_fill_mass, & ! Total mass available to fill hole [kg/m^2]
3833 : hole_fillmass_ratio ! Ratio: total_hole / total_fill_mass [-]
3834 :
3835 : logical, dimension(pver) :: &
3836 : l_pos_hm, & ! Flag: hydrometeor has positive (>= qmin_hm) value
3837 : l_pos_hm_filler ! Flag: filler hydrometeor has positive value
3838 :
3839 : integer :: icol ! Grid column index
3840 :
3841 : integer :: k, idx ! Vertical grid level indices
3842 :
3843 :
3844 : ! Loop over all columns, performing any adjustments one column at a time.
3845 0 : do icol = 1, ncol
3846 :
3847 : ! Calculate the updated value of the hydrometeor field based on the
3848 : ! updated microphysics tendency.
3849 0 : hm_update = hm_start(icol,:) + hm_tend(icol,:) * dt
3850 0 : hm_curr = hm_update
3851 :
3852 : ! Calculate the updated value of the filler hydrometeor field based on
3853 : ! the updated microphysics tendency.
3854 0 : hm_update_filler = hm_start_filler(icol,:) + hm_tend_filler(icol,:) * dt
3855 0 : hm_curr_filler = hm_update_filler
3856 :
3857 : ! Check for any holes in the vertical profile
3858 0 : if ( any( hm_curr(top_lev:pver) < qmin_hm ) ) then
3859 :
3860 : ! At least one hole is found in this hydrometeor species in this
3861 : ! grid column. The holes must be filled conservatively.
3862 :
3863 : ! Check which levels have values of the hydrometeor that are at or
3864 : ! above the minimum threshold value.
3865 0 : do k = top_lev, pver
3866 : ! Check for the hydrometeor that might need to be filled.
3867 0 : if ( hm_curr(k) >= qmin_hm ) then
3868 0 : l_pos_hm(k) = .true.
3869 : else ! hm_curr < qmin_hm
3870 0 : l_pos_hm(k) = .false.
3871 : endif ! hm_curr >= qmin_hm
3872 : ! Check for the filler hydrometeor, as some levels might have
3873 : ! numerical round-off level, small negative values.
3874 0 : if ( hm_curr_filler(k) >= qmin_hm_filler ) then
3875 0 : l_pos_hm_filler(k) = .true.
3876 : else ! hm_curr_filler < qmin_hm_filler
3877 0 : l_pos_hm_filler(k) = .false.
3878 : endif ! hm_curr_filler >= qmin_hm_filler
3879 : end do ! k = top_lev, pver
3880 :
3881 0 : do k = top_lev, pver
3882 :
3883 0 : if ( .not. l_pos_hm(k) ) then
3884 :
3885 : ! A hole is found in the hydrometeor at this grid level.
3886 :
3887 : ! Calculate the total hydrometeor mass of the hole that needs
3888 : ! to be filled.
3889 : ! The value of the hydrometeor mixing ratio is negative, but
3890 : ! the value of total_hole is positive.
3891 0 : total_hole = ( qmin_hm - hm_curr(k) ) * pdel(icol,k) * rga
3892 :
3893 : ! Calculate the total hydrometeor mass available from the
3894 : ! filler hydrometeor to fill the hole.
3895 0 : total_fill_mass = 0.0_r8
3896 0 : do idx = top_lev, pver, 1
3897 0 : if ( l_pos_hm_filler(idx) ) then
3898 : total_fill_mass &
3899 : = total_fill_mass &
3900 : + ( hm_curr_filler(idx) - qmin_hm_filler ) &
3901 0 : * pdel(icol,idx) * rga
3902 : endif ! l_pos_hm_filler(idx)
3903 : end do ! idx = top_lev, pver, 1
3904 :
3905 : ! Calculate the ratio of total hole to total fill mass. This
3906 : ! should not exceed 1 except as a result of numerical round-off
3907 : ! errors. Use thresholding to be safe.
3908 : hole_fillmass_ratio &
3909 : = min( total_hole / max( total_fill_mass, 1.0e-30_r8 ), &
3910 0 : 1.0_r8 )
3911 :
3912 : ! Modify (reduce) the amount of the filler hydrometeor.
3913 0 : do idx = top_lev, pver
3914 0 : if ( l_pos_hm_filler(idx) ) then
3915 : ! Since pdel at a grid level does not change and gravit
3916 : ! is constant, the only variable that needs to be
3917 : ! modified proportionately is hm_curr_filler.
3918 : hm_curr_filler(idx) &
3919 : = qmin_hm_filler &
3920 : + ( hm_curr_filler(idx) - qmin_hm_filler ) &
3921 0 : * ( 1.0_r8 - hole_fillmass_ratio )
3922 : endif ! l_pos_hm_filler(idx)
3923 : end do ! idx = top_lev, pver
3924 :
3925 : ! Update the value of the hydrometeor at the level where the
3926 : ! hole was found. Mathematically, as long as the available
3927 : ! mass was able to fill the entire hole, the new value of the
3928 : ! hydrometeor mixing ratio (hm_curr) should be qmin_hm.
3929 : hm_curr(k) &
3930 : = hm_curr(k) &
3931 : + hole_fillmass_ratio * total_fill_mass &
3932 0 : * gravit / pdel(icol,k)
3933 :
3934 : endif ! .not. l_pos_hm(k)
3935 :
3936 : end do ! k = top_lev, pver
3937 :
3938 : endif ! any( hm_curr(top_lev:pver) < qmin_hm )
3939 :
3940 : ! Update the value of total microphysics tendency after hole filling.
3941 0 : hm_tend(icol,:) = hm_tend(icol,:) + ( hm_curr - hm_update ) / dt
3942 :
3943 : ! Update the value of total microphysics tendency after hole filling for
3944 : ! the filler hydrometeor.
3945 : hm_tend_filler(icol,:) &
3946 0 : = hm_tend_filler(icol,:) + ( hm_curr_filler - hm_update_filler ) / dt
3947 :
3948 : end do ! icol = 1, ncol
3949 :
3950 :
3951 0 : return
3952 :
3953 0 : end subroutine fill_holes_same_phase_vert
3954 :
3955 : !============================================================================
3956 0 : subroutine subcol_SILHS_hydromet_conc_tend_lim( state, dt, ptend )
3957 :
3958 : ! Description:
3959 : ! Limits the values of mean hydrometeor concentrations so that the mean
3960 : ! drop size for the hydrometeor type remains reasonable and does not become
3961 : ! too large.
3962 :
3963 : !----------------------------------------------------------------------
3964 :
3965 0 : use shr_const_mod, only: &
3966 : shr_const_pi, &
3967 : shr_const_rhofw
3968 :
3969 : use constituents, only: &
3970 : qmin
3971 :
3972 : use ref_pres, only: &
3973 : top_lev => trop_cloud_top_lev
3974 :
3975 : implicit none
3976 :
3977 : ! Input Variables
3978 : type(physics_state), intent(in) :: state ! Physics state variables
3979 : real(r8), intent(in) :: dt ! Time step duration
3980 :
3981 : ! Input/Output Variable
3982 : type(physics_ptend), intent(inout) :: ptend ! Parameterization tendencies
3983 :
3984 : ! Local Variables
3985 : real( r8 ) :: &
3986 : rcm_update, & ! New value of mean cloud water mixing ratio [kg/kg]
3987 : rrm_update, & ! New value of mean rain water mixing ratio [kg/kg]
3988 : rim_update, & ! New value of mean ice mixing ratio [kg/kg]
3989 : rsm_update ! New value of mean snow mixing ratio [kg/kg]
3990 :
3991 : real( r8 ) :: &
3992 : Nc_tend_min, & ! Minimum value of cloud droplet conc. tendency [num/kg/s]
3993 : Nr_tend_min, & ! Minimum value of rain drop conc. tendency [num/kg/s]
3994 : Ni_tend_min, & ! Minimum value of ice conc. tendency [num/kg/s]
3995 : Ns_tend_min ! Minimum value of snow conc. tendency [num/kg/s]
3996 :
3997 : real( r8 ), parameter :: &
3998 : four_thirds = 4.0_r8/3.0_r8, & ! 4/3
3999 : rho_ice = 917.0_r8, & ! Density of ice [kg/m^3]
4000 : rho_lw = shr_const_rhofw, & ! Density of liquid water [kg/m^3]
4001 : pi = shr_const_pi ! Pi
4002 :
4003 : real( r8 ), parameter :: &
4004 : mvr_cloud_max = 1.6E-5_r8, & ! Max. avg. mean vol. rad. cloud [m]
4005 : mvr_rain_max = 5.0E-3_r8, & ! Max. avg. mean vol. rad. rain [m]
4006 : mvr_ice_max = 1.3E-4_r8, & ! Max. avg. mean vol. rad. ice [m]
4007 : mvr_snow_max = 1.0E-2_r8 ! Max. avg. mean vol. rad. snow [m]
4008 :
4009 : ! Calculate the coefficient for the minimum mean cloud droplet
4010 : ! concentration, where <Nc>|_min = Ncm_min_coef * <rc> and has units of
4011 : ! 1/kg.
4012 : real( r8 ), parameter :: &
4013 : Ncm_min_coef = 1.0_r8 / ( four_thirds * pi * rho_lw * mvr_cloud_max**3 )
4014 :
4015 : ! Calculate the coefficient for the minimum mean rain drop concentration,
4016 : ! where <Nr>|_min = Nrm_min_coef * <rr> and has units of 1/kg.
4017 : real( r8 ), parameter :: &
4018 : Nrm_min_coef = 1.0_r8 / ( four_thirds * pi * rho_lw * mvr_rain_max**3 )
4019 :
4020 : ! Calculate the coefficient for the minimum mean ice crystal concentration,
4021 : ! where <Ni>|_min = Nim_min_coef * <ri> and has units of 1/kg.
4022 : real( r8 ), parameter :: &
4023 : Nim_min_coef = 1.0_r8 / ( four_thirds * pi * rho_ice * mvr_ice_max**3 )
4024 :
4025 : ! Calculate the coefficient for the minimum mean snow flake concentration,
4026 : ! where <Ns>|_min = Nsm_min_coef * <rs> and has units of 1/kg.
4027 : real( r8 ), parameter :: &
4028 : Nsm_min_coef = 1.0_r8 / ( four_thirds * pi * rho_ice * mvr_snow_max**3 )
4029 :
4030 : integer :: ncol ! Number of grid columns
4031 :
4032 : integer :: icol ! Column loop index
4033 :
4034 : integer :: k ! Vertical level loop index
4035 :
4036 :
4037 : ! Get the number of grid columns.
4038 0 : ncol = state%ncol
4039 :
4040 : ! Loop over all grid columns.
4041 0 : do icol = 1, ncol
4042 :
4043 : ! Loop over all vertical levels from top_lev to pver.
4044 0 : do k = top_lev, pver
4045 :
4046 : ! Cloud droplet concentration
4047 0 : if ( ixcldliq > 0 .and. ixnumliq > 0 ) then
4048 :
4049 : ! Calculate the value of cloud water mixing ratio after the
4050 : ! update.
4051 : rcm_update &
4052 0 : = max( state%q(icol,k,ixcldliq) + ptend%q(icol,k,ixcldliq) * dt, &
4053 0 : qmin(ixcldliq) )
4054 :
4055 : ! Calculate the limiting cloud droplet concentration tendency so
4056 : ! that cloud maintains a reasonable (not too big) mean volume
4057 : ! radius.
4058 : Nc_tend_min &
4059 0 : = ( Ncm_min_coef * rcm_update - state%q(icol,k,ixnumliq) ) / dt
4060 :
4061 : ! The cloud droplet concentration tendency needs to be the greater
4062 : ! of the current Nc_tend and Nc_tend_min.
4063 0 : ptend%q(icol,k,ixnumliq) &
4064 0 : = max( ptend%q(icol,k,ixnumliq), Nc_tend_min )
4065 :
4066 : endif ! ixcldliq > 0 .and. ixnumliq > 0
4067 :
4068 : ! Rain drop concentration
4069 0 : if ( ixrain > 0 .and. ixnumrain > 0 ) then
4070 :
4071 : ! Calculate the value of rain water mixing ratio after the update.
4072 : rrm_update &
4073 0 : = max( state%q(icol,k,ixrain) + ptend%q(icol,k,ixrain) * dt, &
4074 0 : qmin(ixrain) )
4075 :
4076 : ! Calculate the limiting rain drop concentration tendency so that
4077 : ! rain maintains a reasonable (not too big) mean volume radius.
4078 : Nr_tend_min &
4079 0 : = ( Nrm_min_coef * rrm_update - state%q(icol,k,ixnumrain) ) / dt
4080 :
4081 : ! The rain drop concentration tendency needs to be the greater of
4082 : ! the current Nr_tend and Nr_tend_min.
4083 0 : ptend%q(icol,k,ixnumrain) &
4084 0 : = max( ptend%q(icol,k,ixnumrain), Nr_tend_min )
4085 :
4086 : endif ! ixrain > 0 .and. ixnumrain > 0
4087 :
4088 : ! Ice crystal concentration
4089 0 : if ( ixcldice > 0 .and. ixnumice > 0 ) then
4090 :
4091 : ! Calculate the value of ice mixing ratio after the update.
4092 : rim_update &
4093 0 : = max( state%q(icol,k,ixcldice) + ptend%q(icol,k,ixcldice) * dt, &
4094 0 : qmin(ixcldice) )
4095 :
4096 : ! Calculate the limiting ice crystal concentration tendency so
4097 : ! that ice maintains a reasonable (not too big) mean volume
4098 : ! radius.
4099 : Ni_tend_min &
4100 0 : = ( Nim_min_coef * rim_update - state%q(icol,k,ixnumice) ) / dt
4101 :
4102 : ! The ice crystal concentration tendency needs to be the greater
4103 : ! of the current Ni_tend and Ni_tend_min.
4104 0 : ptend%q(icol,k,ixnumice) &
4105 0 : = max( ptend%q(icol,k,ixnumice), Ni_tend_min )
4106 :
4107 : endif ! ixcldice > 0 .and. ixnumice > 0
4108 :
4109 : ! Snow flake concentration
4110 0 : if ( ixsnow > 0 .and. ixnumsnow > 0 ) then
4111 :
4112 : ! Calculate the value of snow mixing ratio after the update.
4113 : rsm_update &
4114 0 : = max( state%q(icol,k,ixsnow) + ptend%q(icol,k,ixsnow) * dt, &
4115 0 : qmin(ixsnow) )
4116 :
4117 : ! Calculate the limiting snow flake concentration tendency so that
4118 : ! snow maintains a reasonable (not too big) mean volume radius.
4119 : Ns_tend_min &
4120 0 : = ( Nsm_min_coef * rsm_update - state%q(icol,k,ixnumsnow) ) / dt
4121 :
4122 : ! The snow flake concentration tendency needs to be the greater of
4123 : ! the current Ns_tend and Ns_tend_min.
4124 0 : ptend%q(icol,k,ixnumsnow) &
4125 0 : = max( ptend%q(icol,k,ixnumsnow), Ns_tend_min )
4126 :
4127 : endif ! ixsnow > 0 .and. ixnumsnow > 0
4128 :
4129 : end do ! k = top_lev, pver
4130 :
4131 : end do ! icol = 1, ncol
4132 :
4133 :
4134 0 : return
4135 :
4136 0 : end subroutine subcol_SILHS_hydromet_conc_tend_lim
4137 :
4138 : !============================================================================
4139 :
4140 : ! Getunit and Freeunit are depreciated in Fortran going forward, so this is a
4141 : ! small function to get an unused stream identifier to send to setup_corr_varnce_array_api
4142 : ! or any other silhs/clubb functions that require a unit number argument
4143 : ! This comes directly from the Fortran wiki
4144 : integer function newunit(unit)
4145 : integer, intent(out), optional :: unit
4146 :
4147 : integer, parameter :: LUN_MIN=10, LUN_MAX=1000
4148 : logical :: opened
4149 : integer :: lun
4150 :
4151 : newunit=-1
4152 : do lun=LUN_MIN,LUN_MAX
4153 : inquire(unit=lun,opened=opened)
4154 : if (.not. opened) then
4155 : newunit=lun
4156 : exit
4157 : end if
4158 : end do
4159 : if (present(unit)) unit=newunit
4160 0 : end function newunit
4161 :
4162 : end module subcol_SILHS
|
