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