Line data Source code
1 : !-----------------------------------------------------------------------
2 : !$Id$
3 : !-------------------------------------------------------------------------------
4 : module corr_varnce_module
5 :
6 : use clubb_precision, only: &
7 : core_rknd
8 :
9 : use error_code, only: &
10 : clubb_at_least_debug_level, & ! Procedure
11 : clubb_fatal_error, & ! Constant
12 : err_code ! Error indicator
13 :
14 : implicit none
15 :
16 : type hmp2_ip_on_hmm2_ip_ratios_type
17 :
18 : ! Prescribed parameters for hydrometeor values of <hm|_ip'^2> / <hm|_ip>^2,
19 : ! where <hm|_ip'^2> is the in-precip. variance of the hydrometeor and
20 : ! <hm|_ip> is the in-precip. mean of the hydrometeor
21 : !
22 : ! These values are dependent on the horizontal grid spacing of the run, and are calculated
23 : ! using a slope and intercept corresponding to each hydrometer.
24 : real( kind = core_rknd ) :: &
25 : rr = 1.0_core_rknd, & ! For rain water mixing ratio [-]
26 : Nr = 1.0_core_rknd, & ! For rain drop concentration [-]
27 : ri = 1.0_core_rknd, & ! For ice mixing ratio [-]
28 : Ni = 1.0_core_rknd, & ! For ice concentration [-]
29 : rs = 1.0_core_rknd, & ! For snow mixing ratio [-]
30 : Ns = 1.0_core_rknd, & ! For snow concentration [-]
31 : rg = 1.0_core_rknd, & ! For graupel mixing ratio [-]
32 : Ng = 1.0_core_rknd ! For graupel concentration [-]
33 :
34 : end type hmp2_ip_on_hmm2_ip_ratios_type
35 :
36 : ! These slopes and intercepts below are used to calculate the hmp2_ip_on_hmm2_ip_ratios_type
37 : ! values that are defined above. This functionality is described by equations 8, 10, & 11
38 : ! in ``Parameterization of the Spatial Variability of Rain for Large-Scale Models and
39 : ! Remote Sensing`` Lebo, et al, October 2015
40 : ! https://journals.ametsoc.org/doi/pdf/10.1175/JAMC-D-15-0066.1
41 : ! see clubb:ticket:830 for more detail
42 : !
43 : ! hmp2_ip_on_hmm2_ip(iirr) = hmp2_ip_on_hmm2_ip_intrcpt%rr + &
44 : ! hmp2_ip_on_hmm2_ip_slope%rr * max( host_dx, host_dy )
45 : !
46 : ! In Lebo et al. the suggested values were
47 : ! slope = 2.12e-5 [1/m]
48 : ! intercept = 0.54 [-]
49 : !
50 : ! In CLUBB standalone, these parameters can be set based on the value for a
51 : ! given case in the CASE_model.in file.
52 : type hmp2_ip_on_hmm2_ip_slope_type
53 :
54 : real( kind = core_rknd ) :: &
55 : rr = 2.12e-5_core_rknd, & ! For rain water mixing ratio [1/m]
56 : Nr = 2.12e-5_core_rknd, & ! For rain drop concentration [1/m]
57 : ri = 2.12e-5_core_rknd, & ! For ice mixing ratio [1/m]
58 : Ni = 2.12e-5_core_rknd, & ! For ice concentration [1/m]
59 : rs = 2.12e-5_core_rknd, & ! For snow mixing ratio [1/m]
60 : Ns = 2.12e-5_core_rknd, & ! For snow concentration [1/m]
61 : rg = 2.12e-5_core_rknd, & ! For graupel mixing ratio [1/m]
62 : Ng = 2.12e-5_core_rknd ! For graupel concentration [1/m]
63 :
64 : end type hmp2_ip_on_hmm2_ip_slope_type
65 :
66 : type hmp2_ip_on_hmm2_ip_intrcpt_type
67 :
68 : real( kind = core_rknd ) :: &
69 : rr = 0.54_core_rknd, & ! For rain water mixing ratio [-]
70 : Nr = 0.54_core_rknd, & ! For rain drop concentration [-]
71 : ri = 0.54_core_rknd, & ! For ice mixing ratio [-]
72 : Ni = 0.54_core_rknd, & ! For ice concentration [-]
73 : rs = 0.54_core_rknd, & ! For snow mixing ratio [-]
74 : Ns = 0.54_core_rknd, & ! For snow concentration [-]
75 : rg = 0.54_core_rknd, & ! For graupel mixing ratio [-]
76 : Ng = 0.54_core_rknd ! For graupel concentration [-]
77 :
78 : end type hmp2_ip_on_hmm2_ip_intrcpt_type
79 :
80 :
81 :
82 : type hm_metadata_type
83 :
84 : ! Variables
85 : ! Microphysics mixing ratios
86 : integer :: &
87 : iirr, & ! Hydrometeor array index for rain water mixing ratio, rr
88 : iirs, & ! Hydrometeor array index for snow mixing ratio, rs
89 : iiri, & ! Hydrometeor array index for ice mixing ratio, ri
90 : iirg ! Hydrometeor array index for graupel mixing ratio, rg
91 :
92 : ! Microphysics concentrations
93 : integer :: &
94 : iiNr, & ! Hydrometeor array index for rain drop concentration, Nr
95 : iiNs, & ! Hydrometeor array index for snow concentration, Ns
96 : iiNi, & ! Hydrometeor array index for ice concentration, Ni
97 : iiNg ! Hydrometeor array index for graupel concentration, Ng
98 :
99 : ! Logical fields
100 : logical, dimension(:), allocatable :: &
101 : l_frozen_hm, & ! if true, then the hydrometeor is frozen; otherwise liquid
102 : l_mix_rat_hm ! if true, then the quantity is a hydrometeor mixing ratio
103 :
104 : character(len=10), dimension(:), allocatable :: &
105 : hydromet_list
106 :
107 : real( kind = core_rknd ), dimension(:), allocatable :: &
108 : hydromet_tol ! Tolerance values for all hydrometeors [units vary]
109 :
110 : ! Latin hypercube indices / Correlation array indices
111 : integer :: &
112 : iiPDF_chi = -1, &
113 : iiPDF_eta = -1, &
114 : iiPDF_w = -1
115 :
116 : integer :: &
117 : iiPDF_rr = -1, &
118 : iiPDF_rs = -1, &
119 : iiPDF_ri = -1, &
120 : iiPDF_rg = -1
121 :
122 : integer :: &
123 : iiPDF_Nr = -1, &
124 : iiPDF_Ns = -1, &
125 : iiPDF_Ni = -1, &
126 : iiPDF_Ng = -1, &
127 : iiPDF_Ncn = -1
128 :
129 : real( kind = core_rknd ), dimension(:), allocatable :: &
130 : hmp2_ip_on_hmm2_ip
131 :
132 : ! Prescribed parameter for <N_cn'^2> / <N_cn>^2.
133 : ! NOTE: In the case that l_const_Nc_in_cloud is true, Ncn is constant
134 : ! throughout the entire grid box, so the parameter below should be
135 : ! ignored.
136 : real( kind = core_rknd ) :: &
137 : Ncnp2_on_Ncnm2 = 1.0_core_rknd ! Prescribed ratio <N_cn'^2> / <N_cn>^2 [-]
138 :
139 : end type hm_metadata_type
140 :
141 : ! This "renaming" is used to shorten the matrix declarations below.
142 : integer, parameter :: c = core_rknd
143 :
144 : integer, parameter, public :: &
145 : d_var_total = 12 ! Size of the default correlation arrays
146 :
147 : real( kind = core_rknd ), public, dimension(d_var_total,d_var_total) :: &
148 : corr_array_n_cloud_def = reshape( &
149 : ! chi eta w Ncn rr Nr ri Ni rs Ns rg Ng
150 : (/1._c,-.6_c, .09_c , .09_c , .788_c, .675_c, .240_c, .222_c, .240_c, .222_c, .240_c, .222_c, & ! chi
151 : 0._c, 1._c, .027_c, .027_c, .114_c, .115_c,-.029_c, .093_c, .022_c, .013_c, 0._c , 0._c , & ! eta
152 : 0._c, 0._c, 1._c , .34_c , .315_c, .270_c, .120_c, .167_c, 0._c , 0._c , 0._c , 0._c , & ! w
153 : 0._c, 0._c, 0._c , 1._c , 0._c , 0._c , .464_c, .320_c, .168_c, .232_c, 0._c , 0._c , & ! Ncn
154 : 0._c, 0._c, 0._c , 0._c , 1._c , .821_c, 0._c , 0._c , .173_c, .164_c, .319_c, .308_c, & ! rr
155 : 0._c, 0._c, 0._c , 0._c , 0._c , 1._c , .152_c, .143_c, 0._c , 0._c , .285_c, .273_c, & ! Nr
156 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 1._c , .758_c, .585_c, .571_c, .379_c, .363_c, & ! ri
157 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 1._c , .571_c, .550_c, .363_c, .345_c, & ! Ni
158 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 1._c , .758_c, .485_c, .470_c, & ! rs
159 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 1._c , .470_c, .450_c, & ! Ns
160 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 1._c , .758_c, & ! rg
161 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 1._c/), & ! Ng
162 : (/d_var_total,d_var_total/) ), &
163 :
164 : corr_array_n_below_def = reshape( &
165 : ! chi eta w Ncn rr Nr ri Ni rs Ns rg Ng
166 : (/1._c, .3_c, .09_c , .09_c , .788_c, .675_c, .240_c, .222_c, .240_c, .222_c, .240_c, .222_c, &! chi
167 : 0._c, 1._c, .027_c, .027_c, .114_c, .115_c,-.029_c, .093_c, .022_c, .013_c, 0._c , 0._c , &! eta
168 : 0._c, 0._c, 1._c , .34_c , .315_c, .270_c, .120_c, .167_c, 0._c , 0._c , 0._c , 0._c , &! w
169 : 0._c, 0._c, 0._c , 1._c , 0._c , 0._c , .464_c, .320_c, .168_c, .232_c, 0._c , 0._c , &! Ncn
170 : 0._c, 0._c, 0._c , 0._c , 1._c , .821_c, 0._c , 0._c , .173_c, .164_c, .319_c, .308_c, &! rr
171 : 0._c, 0._c, 0._c , 0._c , 0._c , 1._c , .152_c, .143_c, 0._c , 0._c , .285_c, .273_c, &! Nr
172 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 1._c , .758_c, .585_c, .571_c, .379_c, .363_c, &! ri
173 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 1._c , .571_c, .550_c, .363_c, .345_c, &! Ni
174 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 1._c , .758_c, .485_c, .470_c, &! rs
175 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 1._c , .470_c, .450_c, &! Ns
176 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 1._c , .758_c, &! rg
177 : 0._c, 0._c, 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 0._c , 1._c/), &! Ng
178 : (/d_var_total,d_var_total/) )
179 :
180 : private
181 :
182 : public :: hmp2_ip_on_hmm2_ip_ratios_type, &
183 : hmp2_ip_on_hmm2_ip_slope_type, &
184 : hmp2_ip_on_hmm2_ip_intrcpt_type, &
185 : read_correlation_matrix, &
186 : setup_corr_varnce_array, &
187 : assert_corr_symmetric, print_corr_matrix, &
188 : hm_metadata_type, &
189 : init_pdf_hydromet_arrays
190 :
191 :
192 : private :: get_corr_var_index, def_corr_idx
193 :
194 :
195 : contains
196 :
197 : !-----------------------------------------------------------------------------
198 0 : function def_corr_idx( iiPDF_x, hm_metadata ) result(ii_def_corr)
199 :
200 : ! Description:
201 : ! Map from a iiPDF index to the corresponding index in the default
202 : ! correlation arrays.
203 : !-----------------------------------------------------------------------------
204 :
205 : implicit none
206 :
207 : ! Constant Parameters
208 :
209 : ! Indices that represent the order in the default corr arrays
210 : ! (chi (old s), eta (old t), w, Ncn, rr, Nr, ri, Ni, rs, Ns, rg, Ng)
211 : integer, parameter :: &
212 : ii_chi = 1, &
213 : ii_eta = 2, &
214 : ii_w = 3, &
215 : ii_Ncn = 4, &
216 : ii_rr = 5, &
217 : ii_Nr = 6, &
218 : ii_ri = 7, &
219 : ii_Ni = 8, &
220 : ii_rs = 9, &
221 : ii_Ns = 10, &
222 : ii_rg = 11, &
223 : ii_Ng = 12
224 :
225 : ! Input Variables
226 :
227 : integer, intent(in) :: iiPDF_x
228 :
229 : ! Return Variable
230 :
231 : integer :: ii_def_corr
232 :
233 : ! ---- Local Variables ----
234 : type (hm_metadata_type) :: &
235 : hm_metadata
236 :
237 : ! ---- Begin Code ----
238 :
239 0 : ii_def_corr = -1
240 :
241 0 : if (iiPDF_x == hm_metadata%iiPDF_chi) then
242 : ii_def_corr = ii_chi
243 :
244 0 : elseif (iiPDF_x == hm_metadata%iiPDF_eta) then
245 : ii_def_corr = ii_eta
246 :
247 0 : elseif (iiPDF_x == hm_metadata%iiPDF_w) then
248 : ii_def_corr = ii_w
249 :
250 0 : elseif (iiPDF_x == hm_metadata%iiPDF_Ncn) then
251 : ii_def_corr = ii_Ncn
252 :
253 0 : elseif (iiPDF_x == hm_metadata%iiPDF_rr) then
254 : ii_def_corr = ii_rr
255 :
256 0 : elseif (iiPDF_x == hm_metadata%iiPDF_Nr) then
257 : ii_def_corr = ii_Nr
258 :
259 0 : elseif (iiPDF_x == hm_metadata%iiPDF_ri) then
260 : ii_def_corr = ii_ri
261 :
262 0 : elseif (iiPDF_x == hm_metadata%iiPDF_Ni) then
263 : ii_def_corr = ii_Ni
264 :
265 0 : elseif (iiPDF_x == hm_metadata%iiPDF_rs) then
266 : ii_def_corr = ii_rs
267 :
268 0 : elseif (iiPDF_x == hm_metadata%iiPDF_Ns) then
269 : ii_def_corr = ii_Ns
270 :
271 0 : elseif (iiPDF_x == hm_metadata%iiPDF_rg) then
272 : ii_def_corr = ii_rg
273 :
274 0 : elseif (iiPDF_x == hm_metadata%iiPDF_Ng) then
275 0 : ii_def_corr = ii_Ng
276 :
277 : endif
278 :
279 0 : end function def_corr_idx
280 :
281 : !-----------------------------------------------------------------------------
282 0 : subroutine set_corr_arrays_to_default( pdf_dim, hm_metadata, &
283 0 : corr_array_n_cloud, corr_array_n_below )
284 :
285 : ! Description:
286 : ! If there are no corr_array.in files for the current case, default
287 : ! correlations are used.
288 : !-----------------------------------------------------------------------------
289 :
290 : use constants_clubb, only: &
291 : zero, &
292 : one
293 :
294 : implicit none
295 :
296 : !------------------------ Input Variables ------------------------
297 : integer, intent(in) :: &
298 : pdf_dim
299 :
300 : type (hm_metadata_type), intent(in) :: &
301 : hm_metadata
302 :
303 : !------------------------ Output Variables ------------------------
304 : real( kind = core_rknd ), intent(out), dimension(pdf_dim,pdf_dim) :: &
305 : corr_array_n_cloud, &
306 : corr_array_n_below
307 :
308 : !------------------------ Local Variables ------------------------
309 : integer :: i, j, idx_i, idx_j ! Loop iterators
310 :
311 : !------------------------ Begin Code ------------------------
312 :
313 0 : corr_array_n_cloud = zero
314 0 : corr_array_n_below = zero
315 :
316 0 : do i = 1, pdf_dim
317 0 : corr_array_n_cloud(i,i) = one
318 0 : corr_array_n_below(i,i) = one
319 : enddo
320 :
321 0 : do i = 1, pdf_dim-1
322 0 : do j = i+1, pdf_dim
323 :
324 0 : idx_i = def_corr_idx(i,hm_metadata)
325 0 : idx_j = def_corr_idx(j,hm_metadata)
326 :
327 0 : if ( idx_i > idx_j ) then
328 0 : corr_array_n_cloud(j, i) = corr_array_n_cloud_def( idx_i, idx_j )
329 0 : corr_array_n_below(j, i) = corr_array_n_below_def( idx_i, idx_j )
330 : else
331 0 : corr_array_n_cloud(j, i) = corr_array_n_cloud_def( idx_j, idx_i )
332 0 : corr_array_n_below(j, i) = corr_array_n_below_def( idx_j, idx_i )
333 : end if
334 : end do
335 : end do
336 :
337 0 : return
338 :
339 : end subroutine set_corr_arrays_to_default
340 :
341 :
342 : !-----------------------------------------------------------------------------
343 0 : subroutine read_correlation_matrix( iunit, input_file, &
344 : pdf_dim, hm_metadata, &
345 0 : corr_array_n )
346 :
347 : ! Description:
348 : ! Reads a correlation variance array from a file and stores it in an array.
349 : !-----------------------------------------------------------------------------
350 :
351 : use input_reader, only: &
352 : one_dim_read_var, & ! Variable(s)
353 : read_one_dim_file, deallocate_one_dim_vars, count_columns ! Procedure(s)
354 :
355 : use matrix_operations, only: set_lower_triangular_matrix ! Procedure(s)
356 :
357 : use constants_clubb, only: fstderr ! Variable(s)
358 :
359 : use clubb_precision, only: &
360 : core_rknd ! Variable(s)
361 :
362 : implicit none
363 :
364 : ! Input Variable(s)
365 : integer, intent(in) :: &
366 : iunit, & ! File I/O unit
367 : pdf_dim! number of variables in the array
368 :
369 : type (hm_metadata_type), intent(in) :: &
370 : hm_metadata
371 :
372 : character(len=*), intent(in) :: input_file ! Path to the file
373 :
374 : ! Input/Output Variable(s)
375 : real( kind = core_rknd ), dimension(pdf_dim,pdf_dim), intent(inout) :: &
376 : corr_array_n ! Normal space correlation array
377 :
378 : ! Local Variable(s)
379 :
380 : type(one_dim_read_var), allocatable, dimension(:) :: &
381 0 : retVars ! stores the variables read in from the corr_varnce.in file
382 :
383 : integer :: &
384 : var_index1, & ! variable index
385 : var_index2, & ! variable index
386 : nCols, & ! the number of columns in the file
387 : i, j ! Loop index
388 :
389 :
390 : !--------------------------- BEGIN CODE -------------------------
391 :
392 0 : nCols = count_columns( iunit, input_file )
393 :
394 : ! Allocate all arrays based on pdf_dim
395 0 : allocate( retVars(1:nCols) )
396 :
397 : ! Initializing to zero means that correlations we don't have are assumed to be 0.
398 0 : corr_array_n(:,:) = 0.0_core_rknd
399 :
400 : ! Set main diagonal to 1
401 0 : do i=1, pdf_dim
402 0 : corr_array_n(i,i) = 1.0_core_rknd
403 : end do
404 :
405 : ! Read the values from the specified file
406 : call read_one_dim_file( iunit, nCols, input_file, & ! intent(in)
407 0 : retVars ) ! intent(out)
408 :
409 0 : if( size( retVars(1)%values ) /= nCols ) then
410 0 : write(fstderr, *) "Correlation matrix must have an equal number of rows and cols in file ", &
411 0 : input_file
412 0 : error stop "Bad data in correlation file."
413 : end if
414 :
415 : ! Start at 2 because the first index is always just 1.0 in the first row
416 : ! and the rest of the rows are ignored
417 0 : do i=2, nCols
418 0 : var_index1 = get_corr_var_index( retVars(i)%name, hm_metadata )
419 0 : if( var_index1 > -1 ) then
420 0 : do j=1, (i-1)
421 0 : var_index2 = get_corr_var_index( retVars(j)%name, hm_metadata )
422 0 : if( var_index2 > -1 ) then
423 : call set_lower_triangular_matrix &
424 0 : ( pdf_dim, var_index1, var_index2, retVars(i)%values(j), & ! intent(in)
425 0 : corr_array_n ) ! intent(inout)
426 : end if
427 : end do
428 : end if
429 : end do
430 :
431 : call deallocate_one_dim_vars( nCols, & ! intent(in)
432 0 : retVars ) ! intent(inout)
433 :
434 0 : return
435 0 : end subroutine read_correlation_matrix
436 :
437 : !--------------------------------------------------------------------------
438 0 : function get_corr_var_index( var_name, hm_metadata ) result( i )
439 :
440 : ! Definition:
441 : ! Returns the index for a variable based on its name.
442 : !--------------------------------------------------------------------------
443 :
444 : implicit none
445 :
446 : character(len=*), intent(in) :: var_name ! The name of the variable
447 :
448 : type (hm_metadata_type), intent(in) :: &
449 : hm_metadata
450 :
451 : ! Output variable
452 : integer :: i
453 :
454 : !------------------ BEGIN CODE -----------------------------
455 0 : i = -1
456 :
457 0 : select case( trim(var_name) )
458 :
459 : case( "chi" )
460 0 : i = hm_metadata%iiPDF_chi
461 :
462 : case( "eta" )
463 0 : i = hm_metadata%iiPDF_eta
464 :
465 : case( "w" )
466 0 : i = hm_metadata%iiPDF_w
467 :
468 : case( "Ncn" )
469 0 : i = hm_metadata%iiPDF_Ncn
470 :
471 : case( "rr" )
472 0 : i = hm_metadata%iiPDF_rr
473 :
474 : case( "Nr" )
475 0 : i = hm_metadata%iiPDF_Nr
476 :
477 : case( "ri" )
478 0 : i = hm_metadata%iiPDF_ri
479 :
480 : case( "Ni" )
481 0 : i = hm_metadata%iiPDF_Ni
482 :
483 : case( "rs" )
484 0 : i = hm_metadata%iiPDF_rs
485 :
486 : case( "Ns" )
487 0 : i = hm_metadata%iiPDF_Ns
488 :
489 : case( "rg" )
490 0 : i = hm_metadata%iiPDF_rg
491 :
492 : case( "Ng" )
493 0 : i = hm_metadata%iiPDF_Ng
494 :
495 : end select
496 :
497 : return
498 :
499 : end function get_corr_var_index
500 :
501 : !================================================================================================
502 : ! subroutine init_pdf_hydromet_arrays
503 : !
504 : ! DESCRIPTION:
505 : ! This subroutine intializes the hydromet arrays(iirr, iiNr, etc.) to the values specified by
506 : ! the input arguements, this determines which hyrometeors are to be used by the microphysics
507 : ! scheme. It also sets up the corresponding pdf and hydromet arrays, and calculates the
508 : ! subgrid variance ratio for each hydrometeor.
509 : !
510 : ! OPTIONAL FUNCTIONALITY:
511 : ! The subgrid variance ratio for each hydrometeor is calculated based on the grid spacing
512 : ! defined by the host model. The calculation is a linear equation defined by a slope and
513 : ! intercept, each of which may or may not be passed in to this subroutine. If the slope
514 : ! and/or intercept are not passed in through the arguement list the default values, which
515 : ! are set in the corresponding type definitions, will be used. Otherwise the values
516 : ! specified by the aruements will be used.
517 : !
518 : ! NOTES:
519 : ! 'hmp2_ip_on_hmm2_ip_slope_in' is of type 'hmp2_ip_on_hmm2_ip_slope_type' and
520 : ! 'hmp2_ip_on_hmm2_ip_intrcpt_in' is of type 'hmp2_ip_on_hmm2_ip_intrcpt_in', both of which
521 : ! are deinfed in corr_vrnce_module.F90, and made public through this API.
522 : !
523 : ! If full control over the hydrometeor variance ratios is desired, pass in slopes that are
524 : ! initialized to 0.0, this causes the ratios to no longer depend on the grid spacing. Then
525 : ! pass in the intercepts set to the values of the desired ratios.
526 : !
527 : ! ARGUEMENTS:
528 : ! host_dx (real) - Horizontal grid spacings
529 : ! host_dy (real)
530 : !
531 : ! hydromet_dim (integer) - Number of enabled hydrometeors
532 : !
533 : ! Each of these is an index value corresponding to a hydrometeor,
534 : ! used to index the hydrometeor arrays. Each index has to be unqiue
535 : ! for each different hyrometeor that is enabled. Setting one of these
536 : ! indices to -1 disables that hydrometeor
537 : ! iirr (integer) - Index of rain water mixing ratio
538 : ! iiri (integer) - Index of rain drop concentration
539 : ! iirs (integer) - Index of ice mixing ratio
540 : ! iirg (integer) - Index of ice crystal concentration
541 : ! iiNr (integer) - Index of snow mixing ratio
542 : ! iiNi (integer) - Index of snow flake concentration
543 : ! iiNs (integer) - Index of graupel mixing ratio
544 : ! iiNg (integer) - Index of graupel concentration
545 : !
546 : ! hmp2_ip_on_hmm2_ip_slope_in (hmp2_ip_on_hmm2_ip_slope_type) - Custom slope values
547 : ! hmp2_ip_on_hmm2_ip_intrcpt_in (hmp2_ip_on_hmm2_ip_intrcpt_type) - Custom intercept values
548 : !
549 : !================================================================================================
550 0 : subroutine init_pdf_hydromet_arrays( host_dx, host_dy, hydromet_dim, & ! intent(in)
551 : iirr, iiNr, iiri, iiNi, & ! intent(in)
552 : iirs, iiNs, iirg, iiNg, & ! intent(in)
553 : Ncnp2_on_Ncnm2, & ! intent(in)
554 : hm_metadata, pdf_dim, & ! intent(out)
555 : hmp2_ip_on_hmm2_ip_slope_in, & ! optional(in)
556 : hmp2_ip_on_hmm2_ip_intrcpt_in ) ! optional(in)
557 : ! Description:
558 : !
559 : ! Initialization for the hydromet arrays. How the arrays are initialized is
560 : ! determined by how the hydromet indicies (iirr, iiNr, etc.) have been set,
561 : ! which are determined by the microphysics scheme.
562 : !-------------------------------------------------------------------------------
563 :
564 : use constants_clubb, only: &
565 : rr_tol, & ! Constants, tolerances for each hydrometeor
566 : ri_tol, &
567 : rs_tol, &
568 : rg_tol, &
569 : Nr_tol, &
570 : Ni_tol, &
571 : Ns_tol, &
572 : Ng_tol
573 :
574 : implicit none
575 :
576 : !------------------------ Input Variables ------------------------
577 : real( kind = core_rknd ), intent(in) :: &
578 : host_dx, host_dy ! Horizontal grid spacing, defined by host model [m]
579 :
580 : integer, intent(in) :: &
581 : hydromet_dim, & ! Total number of hydrometeor species.
582 : iirr, & ! Index of rain water mixing ratio
583 : iiNr, & ! Index of rain drop concentration
584 : iiri, & ! Index of ice mixing ratio
585 : iiNi, & ! Index of ice crystal concentration
586 : iirs, & ! Index of snow mixing ratio
587 : iiNs, & ! Index of snow flake concentration
588 : iirg, & ! Index of graupel mixing ratio
589 : iiNg ! Index of graupel concentration
590 :
591 : real( kind = core_rknd ), intent(in) :: &
592 : Ncnp2_on_Ncnm2
593 :
594 : !------------------------ InOut Variables ------------------------
595 : type (hm_metadata_type), intent(out) :: &
596 : hm_metadata
597 :
598 : !------------------------ Output Variables ------------------------
599 : integer, intent(out) :: &
600 : pdf_dim
601 :
602 : !------------------------ Optional Input Variables ------------------------
603 :
604 : ! Used to overwrite default values of slope and intercept
605 : type( hmp2_ip_on_hmm2_ip_slope_type ), optional, intent(in) :: &
606 : hmp2_ip_on_hmm2_ip_slope_in ! Custom slopes to overwrite defaults [1/m]
607 :
608 : type( hmp2_ip_on_hmm2_ip_intrcpt_type ), optional, intent(in) :: &
609 : hmp2_ip_on_hmm2_ip_intrcpt_in ! Custom intercepts to overwrite defaults [-]
610 :
611 : !------------------------ Local Variables ------------------------
612 : integer :: i, pdf_count
613 :
614 : type( hmp2_ip_on_hmm2_ip_slope_type ) :: &
615 : hmp2_ip_on_hmm2_ip_slope
616 :
617 : type( hmp2_ip_on_hmm2_ip_intrcpt_type ) :: &
618 : hmp2_ip_on_hmm2_ip_intrcpt
619 :
620 : !------------------------ Begin Code ------------------------
621 :
622 : ! Save indices
623 0 : hm_metadata%iirr = iirr
624 0 : hm_metadata%iiri = iiri
625 0 : hm_metadata%iirs = iirs
626 0 : hm_metadata%iirg = iirg
627 0 : hm_metadata%iiNr = iiNr
628 0 : hm_metadata%iiNi = iiNi
629 0 : hm_metadata%iiNs = iiNs
630 0 : hm_metadata%iiNg = iiNg
631 :
632 : ! Overwrite default Ncnp2_on_Ncnm2
633 0 : hm_metadata%Ncnp2_on_Ncnm2 = Ncnp2_on_Ncnm2
634 :
635 : !-----------------------------------------------------------
636 : ! Calculate the subgrid variances of the hydrometeors
637 : !-----------------------------------------------------------
638 :
639 : ! If slope and intercept are present in call, then overwrite default values
640 0 : if ( present( hmp2_ip_on_hmm2_ip_slope_in ) ) then
641 0 : hmp2_ip_on_hmm2_ip_slope = hmp2_ip_on_hmm2_ip_slope_in
642 : end if
643 :
644 0 : if ( present( hmp2_ip_on_hmm2_ip_intrcpt_in ) ) then
645 0 : hmp2_ip_on_hmm2_ip_intrcpt = hmp2_ip_on_hmm2_ip_intrcpt_in
646 : end if
647 :
648 0 : allocate( hm_metadata%hmp2_ip_on_hmm2_ip(hydromet_dim) )
649 0 : hm_metadata%hmp2_ip_on_hmm2_ip = 0.0_core_rknd
650 :
651 0 : if ( iirr > 0 ) then
652 0 : hm_metadata%hmp2_ip_on_hmm2_ip(iirr) &
653 : = hmp2_ip_on_hmm2_ip_intrcpt%rr &
654 0 : + hmp2_ip_on_hmm2_ip_slope%rr * max( host_dx, host_dy )
655 : endif
656 :
657 0 : if ( iirs > 0 ) then
658 0 : hm_metadata%hmp2_ip_on_hmm2_ip(iirs) &
659 : = hmp2_ip_on_hmm2_ip_intrcpt%rs &
660 0 : + hmp2_ip_on_hmm2_ip_slope%rs * max( host_dx, host_dy )
661 : endif
662 :
663 0 : if ( iiri > 0 ) then
664 0 : hm_metadata%hmp2_ip_on_hmm2_ip(iiri) &
665 : = hmp2_ip_on_hmm2_ip_intrcpt%ri &
666 0 : + hmp2_ip_on_hmm2_ip_slope%ri * max( host_dx, host_dy )
667 : endif
668 :
669 0 : if ( iirg > 0 ) then
670 0 : hm_metadata%hmp2_ip_on_hmm2_ip(iirg) &
671 : = hmp2_ip_on_hmm2_ip_intrcpt%rg &
672 0 : + hmp2_ip_on_hmm2_ip_slope%rg * max( host_dx, host_dy )
673 : endif
674 :
675 0 : if ( iiNr > 0 ) then
676 0 : hm_metadata%hmp2_ip_on_hmm2_ip(iiNr) &
677 : = hmp2_ip_on_hmm2_ip_intrcpt%Nr &
678 0 : + hmp2_ip_on_hmm2_ip_slope%Nr * max( host_dx, host_dy )
679 : endif
680 :
681 0 : if ( iiNs > 0 ) then
682 0 : hm_metadata%hmp2_ip_on_hmm2_ip(iiNs) &
683 : = hmp2_ip_on_hmm2_ip_intrcpt%Ns &
684 0 : + hmp2_ip_on_hmm2_ip_slope%Ns * max( host_dx, host_dy )
685 : endif
686 :
687 0 : if ( iiNi > 0 ) then
688 0 : hm_metadata%hmp2_ip_on_hmm2_ip(iiNi) &
689 : = hmp2_ip_on_hmm2_ip_intrcpt%Ni &
690 0 : + hmp2_ip_on_hmm2_ip_slope%Ni * max( host_dx, host_dy )
691 : endif
692 :
693 0 : if ( iiNg > 0 ) then
694 0 : hm_metadata%hmp2_ip_on_hmm2_ip(iiNg) &
695 : = hmp2_ip_on_hmm2_ip_intrcpt%Ng &
696 0 : + hmp2_ip_on_hmm2_ip_slope%Ng * max( host_dx, host_dy )
697 : endif
698 :
699 :
700 : !-----------------------------------------------------------
701 : ! Set up predictive precipitating hydrometeor arrays.
702 : !-----------------------------------------------------------
703 0 : allocate( hm_metadata%hydromet_list(hydromet_dim) )
704 0 : allocate( hm_metadata%hydromet_tol(hydromet_dim) )
705 0 : allocate( hm_metadata%l_mix_rat_hm(hydromet_dim) )
706 0 : allocate( hm_metadata%l_frozen_hm(hydromet_dim) )
707 :
708 0 : if ( iirr > 0 ) then
709 : ! The microphysics scheme predicts rain water mixing ratio, rr.
710 0 : hm_metadata%hydromet_list(iirr) = "rrm"
711 0 : hm_metadata%l_mix_rat_hm(iirr) = .true.
712 0 : hm_metadata%l_frozen_hm(iirr) = .false.
713 0 : hm_metadata%hydromet_tol(iirr) = rr_tol
714 : endif
715 :
716 0 : if ( iiri > 0 ) then
717 : ! The microphysics scheme predicts ice mixing ratio, ri.
718 0 : hm_metadata%hydromet_list(iiri) = "rim"
719 0 : hm_metadata%l_mix_rat_hm(iiri) = .true.
720 0 : hm_metadata%l_frozen_hm(iiri) = .true.
721 0 : hm_metadata%hydromet_tol(iiri) = ri_tol
722 : endif
723 :
724 0 : if ( iirs > 0 ) then
725 : ! The microphysics scheme predicts snow mixing ratio, rs.
726 0 : hm_metadata%hydromet_list(iirs) = "rsm"
727 0 : hm_metadata%l_mix_rat_hm(iirs) = .true.
728 0 : hm_metadata%l_frozen_hm(iirs) = .true.
729 0 : hm_metadata%hydromet_tol(iirs) = rs_tol
730 : endif
731 :
732 0 : if ( iirg > 0 ) then
733 : ! The microphysics scheme predicts graupel mixing ratio, rg.
734 0 : hm_metadata%hydromet_list(iirg) = "rgm"
735 0 : hm_metadata%l_mix_rat_hm(iirg) = .true.
736 0 : hm_metadata%l_frozen_hm(iirg) = .true.
737 0 : hm_metadata%hydromet_tol(iirg) = rg_tol
738 : endif
739 :
740 0 : if ( iiNr > 0 ) then
741 : ! The microphysics scheme predicts rain drop concentration, Nr.
742 0 : hm_metadata%hydromet_list(iiNr) = "Nrm"
743 0 : hm_metadata%l_frozen_hm(iiNr) = .false.
744 0 : hm_metadata%l_mix_rat_hm(iiNr) = .false.
745 0 : hm_metadata%hydromet_tol(iiNr) = Nr_tol
746 : endif
747 :
748 0 : if ( iiNi > 0 ) then
749 : ! The microphysics scheme predicts ice concentration, Ni.
750 0 : hm_metadata%hydromet_list(iiNi) = "Nim"
751 0 : hm_metadata%l_mix_rat_hm(iiNi) = .false.
752 0 : hm_metadata%l_frozen_hm(iiNi) = .true.
753 0 : hm_metadata%hydromet_tol(iiNi) = Ni_tol
754 : endif
755 :
756 0 : if ( iiNs > 0 ) then
757 : ! The microphysics scheme predicts snowflake concentration, Ns.
758 0 : hm_metadata%hydromet_list(iiNs) = "Nsm"
759 0 : hm_metadata%l_mix_rat_hm(iiNs) = .false.
760 0 : hm_metadata%l_frozen_hm(iiNs) = .true.
761 0 : hm_metadata%hydromet_tol(iiNs) = Ns_tol
762 : endif
763 :
764 0 : if ( iiNg > 0 ) then
765 : ! The microphysics scheme predicts graupel concentration, Ng.
766 0 : hm_metadata%hydromet_list(iiNg) = "Ngm"
767 0 : hm_metadata%l_mix_rat_hm(iiNg) = .false.
768 0 : hm_metadata%l_frozen_hm(iiNg) = .true.
769 0 : hm_metadata%hydromet_tol(iiNg) = Ng_tol
770 : endif
771 :
772 :
773 : !-----------------------------------------------------------
774 : ! Set up the PDF indices
775 : !-----------------------------------------------------------
776 :
777 0 : hm_metadata%iiPDF_chi = 1 ! Extended liquid water mixing ratio, chi
778 0 : hm_metadata%iiPDF_eta = 2 ! 'eta' orthogonal to 'chi'
779 0 : hm_metadata%iiPDF_w = 3 ! vertical velocity
780 0 : hm_metadata%iiPDF_Ncn = 4 ! Simplified cloud nuclei concentration or extended Nc.
781 :
782 0 : pdf_count = hm_metadata%iiPDF_Ncn
783 :
784 : ! Loop over hydrometeors.
785 : ! Hydrometeor indices in the PDF arrays should be in the same order as
786 : ! found in the hydrometeor arrays.
787 0 : if ( hydromet_dim > 0 ) then
788 :
789 0 : do i = 1, hydromet_dim, 1
790 :
791 0 : if ( i == iirr ) then
792 0 : pdf_count = pdf_count + 1
793 0 : hm_metadata%iiPDF_rr = pdf_count
794 : endif
795 :
796 0 : if ( i == iiNr ) then
797 0 : pdf_count = pdf_count + 1
798 0 : hm_metadata%iiPDF_Nr = pdf_count
799 : endif
800 :
801 0 : if ( i == iiri ) then
802 0 : pdf_count = pdf_count + 1
803 0 : hm_metadata%iiPDF_ri = pdf_count
804 : endif
805 :
806 0 : if ( i == iiNi ) then
807 0 : pdf_count = pdf_count + 1
808 0 : hm_metadata%iiPDF_Ni = pdf_count
809 : endif
810 :
811 0 : if ( i == iirs ) then
812 0 : pdf_count = pdf_count + 1
813 0 : hm_metadata%iiPDF_rs = pdf_count
814 : endif
815 :
816 0 : if ( i == iiNs ) then
817 0 : pdf_count = pdf_count + 1
818 0 : hm_metadata%iiPDF_Ns = pdf_count
819 : endif
820 :
821 0 : if ( i == iirg ) then
822 0 : pdf_count = pdf_count + 1
823 0 : hm_metadata%iiPDF_rg = pdf_count
824 : endif
825 :
826 0 : if ( i == iiNg ) then
827 0 : pdf_count = pdf_count + 1
828 0 : hm_metadata%iiPDF_Ng = pdf_count
829 : endif
830 :
831 : enddo ! i = 1, hydromet_dim, 1
832 :
833 : endif ! hydromet_dim > 0
834 :
835 0 : pdf_dim = pdf_count
836 :
837 0 : return
838 :
839 0 : end subroutine init_pdf_hydromet_arrays
840 :
841 : !===============================================================================
842 0 : subroutine setup_corr_varnce_array( input_file_cloud, input_file_below, &
843 : pdf_dim, hm_metadata, iunit, &
844 : l_fix_w_chi_eta_correlations, &
845 0 : corr_array_n_cloud, corr_array_n_below )
846 :
847 : ! Description:
848 : ! Setup an array with the x'^2/xm^2 variables on the diagonal and the other
849 : ! elements to be correlations between various variables.
850 :
851 : ! References:
852 : ! None.
853 : !-------------------------------------------------------------------------------
854 :
855 : use matrix_operations, only: mirror_lower_triangular_matrix ! Procedure
856 :
857 : use constants_clubb, only: &
858 : fstderr, & ! Constant(s)
859 : eps
860 :
861 : implicit none
862 :
863 : ! External
864 : intrinsic :: max, epsilon, trim
865 :
866 :
867 : !------------------------ Input Variables ------------------------
868 : character(len=*), intent(in) :: &
869 : input_file_cloud, & ! Path to the in cloud correlation file
870 : input_file_below ! Path to the out of cloud correlation file
871 :
872 : integer, intent(in) :: &
873 : pdf_dim
874 :
875 : type (hm_metadata_type), intent(in) :: &
876 : hm_metadata
877 :
878 : integer, intent(in) :: &
879 : iunit ! The file unit
880 :
881 : logical, intent(in) :: &
882 : l_fix_w_chi_eta_correlations ! Use a fixed correlation for s and t Mellor(chi/eta)
883 :
884 : !------------------------ Output Variables ------------------------
885 : real( kind = core_rknd ), intent(out), dimension(pdf_dim,pdf_dim) :: &
886 : corr_array_n_cloud, &
887 : corr_array_n_below
888 :
889 : !------------------------ Local variables ------------------------
890 : logical :: l_warning, l_corr_file_1_exist, l_corr_file_2_exist
891 : integer :: i
892 :
893 : !------------------------ Begin Code ------------------------
894 :
895 0 : inquire( file = input_file_cloud, exist = l_corr_file_1_exist )
896 0 : inquire( file = input_file_below, exist = l_corr_file_2_exist )
897 :
898 0 : if ( l_corr_file_1_exist .and. l_corr_file_2_exist ) then
899 :
900 : call read_correlation_matrix( iunit, trim( input_file_cloud ), &
901 : pdf_dim, hm_metadata, & ! In
902 0 : corr_array_n_cloud ) ! Out
903 :
904 : call read_correlation_matrix( iunit, trim( input_file_below ), &
905 : pdf_dim, hm_metadata, & ! In
906 0 : corr_array_n_below ) ! Out
907 :
908 : else ! Read in default correlation matrices
909 :
910 0 : if ( clubb_at_least_debug_level( 1 ) ) then
911 : write(fstderr,*) "Warning: "//trim( input_file_cloud )//" was not found! " // &
912 0 : "The default correlation arrays will be used."
913 : end if
914 :
915 : call set_corr_arrays_to_default( pdf_dim, hm_metadata, &
916 0 : corr_array_n_cloud, corr_array_n_below )
917 :
918 : endif
919 :
920 : ! Mirror the correlation matrices
921 : call mirror_lower_triangular_matrix( pdf_dim, & ! intent(in)
922 0 : corr_array_n_cloud ) ! intent(inout)
923 :
924 : call mirror_lower_triangular_matrix( pdf_dim, & ! intent(in)
925 0 : corr_array_n_below ) ! intent(inout)
926 :
927 :
928 : ! Sanity check to avoid confusing non-convergence results.
929 0 : if ( clubb_at_least_debug_level( 2 ) ) then
930 :
931 0 : if ( .not. l_fix_w_chi_eta_correlations .and. hm_metadata%iiPDF_Ncn > 0 ) then
932 : l_warning = .false.
933 0 : do i = 1, pdf_dim
934 0 : if ( ( (abs(corr_array_n_cloud(i,hm_metadata%iiPDF_Ncn)) > eps) .or. &
935 0 : (abs(corr_array_n_below(i,hm_metadata%iiPDF_Ncn))) > eps) .and. &
936 0 : i /= hm_metadata%iiPDF_Ncn ) then
937 0 : l_warning = .true.
938 : end if
939 : end do ! 1..pdf_dim
940 0 : if ( l_warning ) then
941 : write(fstderr,*) "Warning: the specified correlations for chi" &
942 0 : // " (old s) and Ncn are non-zero."
943 : write(fstderr,*) "The latin hypercube code will not converge to" &
944 0 : // " the analytic solution using these settings."
945 : end if
946 : end if ! l_fix_w_chi_eta_correlations .and. iiPDF_Ncn > 0
947 :
948 : end if ! clubb_at_least_debug_level( 2 )
949 :
950 :
951 0 : return
952 :
953 : end subroutine setup_corr_varnce_array
954 :
955 : !-----------------------------------------------------------------------------
956 0 : subroutine assert_corr_symmetric( corr_array_n, & ! intent(in)
957 : pdf_dim) ! intent(in)
958 :
959 : ! Description:
960 : ! Asserts that corr_matrix(i,j) == corr_matrix(j,i) for all indeces
961 : ! in the correlation array. If this is not the case, stops the program.
962 : ! References:
963 : ! None
964 : !---------------------------------------------------------------------------
965 :
966 : use constants_clubb, only: fstderr, eps, one ! Constant(s)
967 :
968 : implicit none
969 :
970 : ! Input Variables
971 : integer, intent(in) :: &
972 : pdf_dim ! Number of variables in the correlation array
973 :
974 : real( kind = core_rknd ), dimension(pdf_dim, pdf_dim), &
975 : intent(in) :: corr_array_n ! Normal space correlation array to be checked
976 :
977 : ! Local Variables
978 :
979 : ! tolerance used for real precision testing
980 : real( kind = core_rknd ), parameter :: tol = 1.0e-6_core_rknd
981 :
982 : integer :: n_row, n_col ! Indices
983 :
984 : !----- Begin Code -----
985 :
986 : !Do the check
987 0 : do n_col = 1, pdf_dim
988 0 : do n_row = 1, pdf_dim
989 :
990 0 : if ( abs(corr_array_n(n_col, n_row) - corr_array_n(n_row, n_col)) > tol ) then
991 0 : err_code = clubb_fatal_error
992 0 : write(fstderr,*) "in subroutine assert_corr_symmetric: ", &
993 0 : "Correlation array is non symmetric."
994 0 : write(fstderr,*) corr_array_n
995 0 : return
996 : end if
997 :
998 0 : if ( n_col == n_row .and. abs(corr_array_n(n_col, n_row)-one) > eps ) then
999 0 : err_code = clubb_fatal_error
1000 0 : write(fstderr,*) "in subroutine assert_corr_symmetric: ", &
1001 0 : "Correlation array is formatted incorrectly."
1002 0 : write(fstderr,*) corr_array_n
1003 0 : return
1004 : end if
1005 : end do
1006 : end do
1007 :
1008 : end subroutine assert_corr_symmetric
1009 :
1010 : !-----------------------------------------------------------------------------
1011 0 : subroutine print_corr_matrix( pdf_dim, & ! intent(in)
1012 0 : corr_array_n ) ! intent(in)
1013 :
1014 : ! Description:
1015 : ! Prints the correlation matrix to the console.
1016 : ! References:
1017 : ! None
1018 : !---------------------------------------------------------------------------
1019 :
1020 : use clubb_precision, only: core_rknd
1021 :
1022 : implicit none
1023 :
1024 : ! Input Variables
1025 : integer, intent(in) :: &
1026 : pdf_dim ! Number of variables in the correlation array
1027 :
1028 : real( kind = core_rknd ), dimension(pdf_dim, pdf_dim), &
1029 : intent(in) :: corr_array_n ! Normal space correlation array to be printed
1030 :
1031 : ! Local Variables
1032 : integer :: n, & ! Loop indeces
1033 : m, &
1034 : current_character_index ! keeps track of the position in the string
1035 :
1036 : character(LEN=72) :: current_line ! The current line to be printed
1037 : character(LEN=10) :: str_array_value
1038 :
1039 : !----- Begin Code -----
1040 :
1041 0 : current_character_index = 0
1042 :
1043 0 : do n = 1, pdf_dim
1044 0 : do m = 1, pdf_dim
1045 0 : write(str_array_value,'(F5.2)') corr_array_n(m,n)
1046 0 : current_line = current_line(1:current_character_index)//str_array_value
1047 0 : current_character_index = current_character_index + 6
1048 : end do
1049 0 : write(*, *) current_line
1050 0 : current_line = ""
1051 0 : current_character_index = 0
1052 : end do
1053 :
1054 0 : end subroutine print_corr_matrix
1055 : !-----------------------------------------------------------------------------
1056 :
1057 0 : end module corr_varnce_module
|