Line data Source code
1 : !-----------------------------------------------------------------------
2 : ! $Id$
3 : !===============================================================================
4 : module diagnose_correlations_module
5 :
6 : use clubb_precision, only: &
7 : core_rknd
8 :
9 : implicit none
10 :
11 : public :: calc_mean, calc_varnce, calc_w_corr, &
12 : calc_cholesky_corr_mtx_approx, &
13 : cholesky_to_corr_mtx_approx, setup_corr_cholesky_mtx, &
14 : diagnose_correlations
15 :
16 :
17 : private :: diagnose_corr, rearrange_corr_array, &
18 : corr_array_assertion_checks
19 :
20 : private ! Default scope
21 : contains
22 :
23 : !-----------------------------------------------------------------------
24 0 : subroutine diagnose_correlations( pdf_dim, iiPDF_w, & ! Intent(in)
25 0 : corr_array_pre, & ! Intent(in)
26 : l_calc_w_corr, & ! Intent(in)
27 0 : corr_array ) ! Intent(out)
28 : ! Description:
29 : ! This subroutine diagnoses the correlation matrix in order to feed it
30 : ! into SILHS microphysics.
31 :
32 : ! References:
33 : ! Larson et al. (2011), J. of Geophysical Research, Vol. 116, D00T02
34 : ! (see CLUBB Trac ticket#514)
35 : !-----------------------------------------------------------------------
36 :
37 : use clubb_precision, only: &
38 : core_rknd ! Variable(s)
39 :
40 : use constants_clubb, only: &
41 : zero
42 :
43 : implicit none
44 :
45 : intrinsic :: max, sqrt, transpose
46 :
47 : ! Input Variables
48 : integer, intent(in) :: &
49 : pdf_dim, & ! number of diagnosed correlations
50 : iiPDF_w
51 :
52 : real( kind = core_rknd ), dimension(pdf_dim, pdf_dim), intent(in) :: &
53 : corr_array_pre ! Prescribed correlations
54 :
55 : logical, intent(in) :: &
56 : l_calc_w_corr ! Calculate the correlations between w and the hydrometeors
57 :
58 : ! Output variables
59 : real( kind = core_rknd ), dimension(pdf_dim, pdf_dim), intent(out) :: &
60 : corr_array
61 :
62 : ! Local Variables
63 : real( kind = core_rknd ), dimension(pdf_dim, pdf_dim) :: &
64 0 : corr_array_pre_swapped, &
65 0 : corr_array_swapped
66 :
67 : ! We actually don't need this right now
68 : real( kind = core_rknd ), dimension(pdf_dim) :: &
69 0 : sigma2_on_mu2_ip_array ! Ratios: sigma_x^2/mu_x^2 (ith PDF comp.) ip [-]
70 :
71 : integer :: i ! Loop iterator
72 :
73 : !-------------------- Begin code --------------------
74 :
75 : ! Initialize sigma2_on_mu2_ip_array
76 0 : do i = 1, pdf_dim
77 0 : sigma2_on_mu2_ip_array(i) = zero
78 : end do
79 :
80 : ! Swap the w-correlations to the first row for the prescribed correlations
81 : call rearrange_corr_array( pdf_dim, iiPDF_w, & ! Intent(in)
82 : corr_array_pre, & ! Intent(in)
83 0 : corr_array_pre_swapped ) ! Intent(inout)
84 :
85 : ! diagnose correlations
86 :
87 0 : if ( .not. l_calc_w_corr ) then
88 0 : corr_array_swapped = corr_array_pre_swapped
89 : endif
90 :
91 : call diagnose_corr( pdf_dim, sqrt(sigma2_on_mu2_ip_array), & ! intent(in)
92 : corr_array_pre_swapped, & ! intent(in)
93 0 : corr_array_swapped ) ! intent(inout)
94 :
95 : ! Swap rows back
96 : call rearrange_corr_array( pdf_dim, iiPDF_w, & ! Intent(in)
97 : corr_array_swapped, & ! Intent(in)
98 0 : corr_array) ! Intent(out)
99 :
100 0 : end subroutine diagnose_correlations
101 :
102 :
103 : !-----------------------------------------------------------------------
104 0 : subroutine diagnose_corr( n_variables, sqrt_sigma2_on_mu2_ip, & ! intent(in)
105 0 : corr_matrix_prescribed, & !intent(in)
106 0 : corr_matrix_approx ) ! intent(inout)
107 :
108 : ! Description:
109 : ! This subroutine diagnoses the correlation matrix for each timestep.
110 :
111 : ! References:
112 : ! Larson et al. (2011), J. of Geophysical Research, Vol. 116, D00T02
113 : ! (see CLUBB Trac ticket#514)
114 : !-----------------------------------------------------------------------
115 :
116 : use clubb_precision, only: &
117 : core_rknd ! Variable(s)
118 :
119 : use constants_clubb, only: &
120 : max_mag_correlation
121 :
122 : implicit none
123 :
124 : intrinsic :: &
125 : sqrt, abs, sign
126 :
127 : ! Input Variables
128 : integer, intent(in) :: &
129 : n_variables ! number of variables in the correlation matrix [-]
130 :
131 : real( kind = core_rknd ), dimension(n_variables), intent(in) :: &
132 : sqrt_sigma2_on_mu2_ip ! sqrt of sigma_x^2/mu_x^2 (ith PDF comp.) ip [-]
133 :
134 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(in) :: &
135 : corr_matrix_prescribed ! correlation matrix [-]
136 :
137 : ! Input/Output Variables
138 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(inout) :: &
139 : corr_matrix_approx ! correlation matrix [-]
140 :
141 :
142 : ! Local Variables
143 : integer :: i, j ! Loop iterator
144 :
145 : real( kind = core_rknd ) :: &
146 : f_ij
147 : ! f_ij_o
148 :
149 : real( kind = core_rknd ), dimension(n_variables) :: &
150 0 : s_1j ! s_1j = sqrt(1-c_1j^2)
151 :
152 :
153 : !-------------------- Begin code --------------------
154 :
155 : ! Remove compiler warnings about unused variables.
156 : if ( .false. ) then
157 : print *, "sqrt_sigma2_on_mu2_ip = ", sqrt_sigma2_on_mu2_ip
158 : endif
159 :
160 : ! calculate all square roots
161 0 : do i = 1, n_variables
162 :
163 0 : s_1j(i) = sqrt(1._core_rknd-corr_matrix_approx(i,1)**2)
164 :
165 : end do
166 :
167 :
168 : ! Diagnose the missing correlations (upper triangle)
169 0 : do j = 2, (n_variables-1)
170 0 : do i = (j+1), n_variables
171 :
172 : ! formula (16) in the ref. paper (Larson et al. (2011))
173 : !f_ij = alpha_corr * sqrt_sigma2_on_mu2_ip(i) * sqrt_sigma2_on_mu2_ip(j) &
174 : ! * sign(1.0_core_rknd,corr_matrix_approx(1,i)*corr_matrix_approx(1,j))
175 :
176 : ! If the predicting c1i's are small then cij will be closer to the prescribed value. If
177 : ! the c1i's are bigger, then cij will be closer to formular (15) from the ref. paper. See
178 : ! clubb:ticket:514:comment:61 for details.
179 : !f_ij = (1-abs(corr_matrix_approx(1,i)*corr_matrix_approx(1,j)))*corr_matrix_prescribed(i,j) &
180 : ! + abs(corr_matrix_approx(1,i)*corr_matrix_approx(1,j))*f_ij_o
181 :
182 0 : f_ij = corr_matrix_prescribed(i,j)
183 :
184 : ! make sure -1 < f_ij < 1
185 0 : if ( f_ij < -max_mag_correlation ) then
186 :
187 : f_ij = -max_mag_correlation
188 :
189 0 : else if ( f_ij > max_mag_correlation ) then
190 :
191 0 : f_ij = max_mag_correlation
192 :
193 : end if
194 :
195 :
196 : ! formula (15) in the ref. paper (Larson et al. (2011))
197 0 : corr_matrix_approx(i,j) = corr_matrix_approx(i,1) * corr_matrix_approx(j,1) &
198 0 : + f_ij * s_1j(i) * s_1j(j)
199 :
200 : end do ! do j
201 : end do ! do i
202 :
203 0 : end subroutine diagnose_corr
204 :
205 :
206 : !-----------------------------------------------------------------------
207 : ! subroutine approx_w_corr( nz, pdf_dim, pdf_params, & ! Intent(in)
208 : ! rrm, Nrm, Ncnm, &
209 : ! stdev_w, sigma_rr_1, &
210 : ! sigma_Nr_1, sigma_Ncn_1, &
211 : ! corr_array) ! Intent(out)
212 : ! ! Description:
213 : ! ! Approximate the correlations of w with the hydrometeors.
214 : !
215 : ! ! References:
216 : ! ! clubb:ticket:514
217 : ! !-----------------------------------------------------------------------
218 : !
219 : ! use clubb_precision, only: &
220 : ! core_rknd ! Variable(s)
221 : !
222 : ! use pdf_parameter_module, only: &
223 : ! pdf_parameter ! Type
224 : !
225 : ! use constants_clubb, only: &
226 : ! one, & ! Constant(s)
227 : ! rr_tol, &
228 : ! Nr_tol, &
229 : ! Ncn_tol, &
230 : ! w_tol, & ! [m/s]
231 : ! chi_tol ! [kg/kg]
232 : !
233 : ! implicit none
234 : !
235 : ! ! Input Variables
236 : ! integer, intent(in) :: &
237 : ! pdf_dim, & ! Number of diagnosed correlations
238 : ! nz ! Number of model vertical grid levels
239 : !
240 : ! type(pdf_parameter), dimension(nz), intent(in) :: &
241 : ! pdf_params ! PDF parameters [units vary]
242 : !
243 : ! real( kind = core_rknd ), dimension(nz), intent(in) :: &
244 : ! rrm, & ! Mean rain water mixing ratio, < r_r > [kg/kg]
245 : ! Nrm, & ! Mean rain drop concentration, < N_r > [num/kg]
246 : ! Ncnm, & ! Mean cloud nuclei conc., < N_cn > [num/kg]
247 : ! stdev_w ! Standard deviation of w [m/s]
248 : !
249 : ! real( kind = core_rknd ), intent(in) :: &
250 : ! sigma_Ncn_1, & ! Standard deviation of Ncn (1st PDF component) [num/kg]
251 : ! sigma_Nr_1, & ! Standard deviation of Nr (2nd PDF component) [num/kg]
252 : ! sigma_rr_1 ! Standard dev. of ln rr (1st PDF comp.) ip [ln(kg/kg)]
253 : !
254 : ! ! Output Variables
255 : ! real( kind = core_rknd ), dimension(pdf_dim, pdf_dim, nz), intent(out) :: &
256 : ! corr_array
257 : !
258 : ! ! Local Variables
259 : ! real( kind = core_rknd ), dimension(nz) :: &
260 : ! corr_chi_w, & ! Correlation between w and chi(s_mellor) (both components) [-]
261 : ! corr_wrr, & ! Correlation between w and rr (both components) [-]
262 : ! corr_wNr, & ! Correlation between w and Nr (both components) [-]
263 : ! corr_wNcn ! Correlation between w and Ncn (both components) [-]
264 : !
265 : ! real( kind = core_rknd ), dimension(nz) :: &
266 : ! wpchip_zt, & ! Covariance of chi and w on the zt-grid [(m/s)(kg/kg)]
267 : ! wprrp_zt, & ! Covariance of r_r and w on the zt-grid [(m/s)(kg/kg)]
268 : ! wpNrp_zt, & ! Covariance of N_r and w on the zt-grid [(m/s)(#/kg)]
269 : ! wpNcnp_zt ! Covariance of N_cn and w on the zt-grid [(m/s)(#/kg)]
270 : !
271 : ! real( kind = core_rknd ) :: &
272 : ! chi_m, & ! Mean of chi (s_mellor) [kg/kg]
273 : ! stdev_chi ! Standard deviation of chi (s_mellor) [kg/kg]
274 : !
275 : ! integer :: k ! vertical loop iterator
276 : !
277 : ! ! ----- Begin Code -----
278 : !
279 : ! call approx_w_covar( nz, pdf_params, rrm, Nrm, Ncnm, & ! Intent(in)
280 : ! wpchip_zt, wprrp_zt, wpNrp_zt, wpNcnp_zt ) ! Intent(out)
281 : !
282 : ! do k = 1, nz
283 : !
284 : ! chi_m &
285 : ! = calc_mean( pdf_params(k)%mixt_frac, pdf_params(k)%chi_1, &
286 : ! pdf_params(k)%chi_2 )
287 : !
288 : ! stdev_chi &
289 : ! = sqrt( pdf_params(k)%mixt_frac &
290 : ! * ( ( pdf_params(k)%chi_1 - chi_m )**2 &
291 : ! + pdf_params(k)%stdev_chi_1**2 ) &
292 : ! + ( one - pdf_params(k)%mixt_frac ) &
293 : ! * ( ( pdf_params(k)%chi_2 - chi_m )**2 &
294 : ! + pdf_params(k)%stdev_chi_2**2 ) &
295 : ! )
296 : !
297 : ! corr_chi_w(k) &
298 : ! = calc_w_corr( wpchip_zt(k), stdev_w(k), stdev_chi, &
299 : ! w_tol, chi_tol )
300 : !
301 : ! corr_wrr(k) &
302 : ! = calc_w_corr( wprrp_zt(k), stdev_w(k), sigma_rr_1, w_tol, rr_tol )
303 : !
304 : ! corr_wNr(k) &
305 : ! = calc_w_corr( wpNrp_zt(k), stdev_w(k), sigma_Nr_1, w_tol, Nr_tol )
306 : !
307 : ! corr_wNcn(k) &
308 : ! = calc_w_corr( wpNcnp_zt(k), stdev_w(k), sigma_Ncn_1, w_tol, Ncn_tol )
309 : !
310 : ! enddo
311 : !
312 : ! call set_w_corr( nz, pdf_dim, & ! Intent(in)
313 : ! corr_chi_w, corr_wrr, corr_wNr, corr_wNcn, &
314 : ! corr_array ) ! Intent(inout)
315 : !
316 : ! end subroutine approx_w_corr
317 :
318 :
319 : !-----------------------------------------------------------------------
320 : ! subroutine approx_w_covar( nz, pdf_params, rrm, Nrm, Ncnm, Kh_zm, & ! Intent(in)
321 : ! wpchip_zt, wprrp_zt, wpNrp_zt, wpNcnp_zt ) ! Intent(out)
322 : ! ! Description:
323 : ! ! Approximate the covariances of w with the hydrometeors using Eddy
324 : ! ! diffusivity.
325 : !
326 : ! ! References:
327 : ! ! clubb:ticket:514
328 : ! !-----------------------------------------------------------------------
329 : !
330 : ! use clubb_precision, only: &
331 : ! core_rknd ! Variable(s)
332 : !
333 : ! use grid_class, only: &
334 : ! gr, & ! Variable(s)
335 : ! zm2zt, & ! Procedure(s)
336 : ! zt2zm
337 : !
338 : ! use pdf_parameter_module, only: &
339 : ! pdf_parameter ! Type
340 : !
341 : ! use constants_clubb, only: &
342 : ! one ! Constant(s)
343 : !
344 : ! use advance_windm_edsclrm_module, only: &
345 : ! xpwp_fnc ! Procedure(s)
346 : !
347 : ! implicit none
348 : !
349 : ! ! Input Variables
350 : ! integer, intent(in) :: &
351 : ! nz ! Number of model vertical grid levels
352 : !
353 : ! type(pdf_parameter), dimension(nz), intent(in) :: &
354 : ! pdf_params ! PDF parameters [units vary]
355 : !
356 : ! real( kind = core_rknd ), dimension(nz), intent(in) :: &
357 : ! rrm, & ! Mean rain water mixing ratio, < r_r > [kg/kg]
358 : ! Nrm, & ! Mean rain drop concentration, < N_r > [num/kg]
359 : ! Ncnm, & ! Mean cloud nuclei concentration, < N_cn > [num/kg]
360 : ! Kh_zm ! Eddy diffusivity coef. on momentum levels [m^2/s]
361 : !
362 : ! ! Output Variables
363 : ! real( kind = core_rknd ), dimension(nz), intent(out) :: &
364 : ! wpchip_zt, & ! Covariance of chi(s) and w on the zt-grid [(m/s)(kg/kg)]
365 : ! wprrp_zt, & ! Covariance of r_r and w on the zt-grid [(m/s)(kg/kg)]
366 : ! wpNrp_zt, & ! Covariance of N_r and w on the zt-grid [(m/s)(#/kg)]
367 : ! wpNcnp_zt ! Covariance of N_cn and w on the zt-grid [(m/s)(#/kg)]
368 : !
369 : ! ! Local Variables
370 : ! real( kind = core_rknd ), dimension(nz) :: &
371 : ! wpchip_zm, & ! Covariance of chi(s) and w on the zm-grid [(m/s)(kg/kg)]
372 : ! wprrp_zm, & ! Covariance of r_r and w on the zm-grid [(m/s)(kg/kg)]
373 : ! wpNrp_zm, & ! Covariance of N_r and w on the zm-grid [(m/s)(#/kg)]
374 : ! wpNcnp_zm ! Covariance of N_cn and w on the zm-grid [(m/s)(#/kg)]
375 : !
376 : ! integer :: k ! vertical loop iterator
377 : !
378 : ! ! ----- Begin Code -----
379 : !
380 : ! ! calculate the covariances of w with the hydrometeors
381 : ! do k = 1, nz
382 : ! wpchip_zm(k) = pdf_params(k)%mixt_frac &
383 : ! * ( one - pdf_params(k)%mixt_frac ) &
384 : ! * ( pdf_params(k)%chi_1 - pdf_params(k)%chi_2 ) &
385 : ! * ( pdf_params(k)%w_1 - pdf_params(k)%w_2 )
386 : ! enddo
387 : !
388 : !! same for wpNrp
389 : !! wprrp_zm(1:nz-1) &
390 : !! = xpwp_fnc( -c_K_hm * Kh_zm(1:nz-1), &
391 : !! rrm(1:nz-1) / max( precip_frac(1:nz-1), eps ), &
392 : !! rrm(2:nz) / max( precip_frac(2:nz), eps ), &
393 : !! gr%invrs_dzm(1:nz-1) )
394 : !
395 : ! wprrp_zm(1:nz-1) &
396 : ! = xpwp_fnc( -c_K_hm * Kh_zm(1:nz-1), &
397 : ! rrm(1:nz-1), rrm(2:nz), &
398 : ! gr%invrs_dzm(1:nz-1) )
399 : !
400 : ! wpNrp_zm(1:nz-1) &
401 : ! = xpwp_fnc( -c_K_hm * Kh_zm(1:nz-1), &
402 : ! Nrm(1:nz-1), Nrm(2:nz), &
403 : ! gr%invrs_dzm(1:nz-1) )
404 : !
405 : ! wpNcnp_zm(1:nz-1) = xpwp_fnc( -c_K_hm * Kh_zm(1:nz-1), Ncnm(1:nz-1), &
406 : ! Ncnm(2:nz), gr%invrs_dzm(1:nz-1) )
407 : !
408 : ! ! Boundary conditions; We are assuming constant flux at the top.
409 : ! wprrp_zm(nz) = wprrp_zm(nz-1)
410 : ! wpNrp_zm(nz) = wpNrp_zm(nz-1)
411 : ! wpNcnp_zm(nz) = wpNcnp_zm(nz-1)
412 : !
413 : ! ! interpolate back to zt-grid
414 : ! wpchip_zt = zm2zt(wpchip_zm)
415 : ! wprrp_zt = zm2zt(wprrp_zm)
416 : ! wpNrp_zt = zm2zt(wpNrp_zm)
417 : ! wpNcnp_zt = zm2zt(wpNcnp_zm)
418 : !
419 : ! end subroutine approx_w_covar
420 :
421 : !-----------------------------------------------------------------------
422 0 : function calc_w_corr( wpxp, stdev_w, stdev_x, w_tol, x_tol )
423 : ! Description:
424 : ! Compute the correlations of w with the hydrometeors.
425 :
426 : ! References:
427 : ! clubb:ticket:514
428 : !-----------------------------------------------------------------------
429 :
430 : use clubb_precision, only: &
431 : core_rknd ! Variable(s)
432 :
433 : use constants_clubb, only: &
434 : max_mag_correlation
435 :
436 : implicit none
437 :
438 : intrinsic :: max
439 :
440 : ! Input Variables
441 : real( kind = core_rknd ), intent(in) :: &
442 : stdev_w, & ! standard deviation of w [m/s]
443 : stdev_x, & ! standard deviation of x [units vary]
444 : wpxp, & ! Covariances of w with the hydrometeors [units vary]
445 : w_tol, & ! tolerance for w [m/s]
446 : x_tol ! tolerance for x [units vary]
447 :
448 : real( kind = core_rknd ) :: &
449 : calc_w_corr
450 :
451 : ! --- Begin Code ---
452 :
453 0 : calc_w_corr = wpxp / ( max(stdev_x, x_tol) * max(stdev_w, w_tol) )
454 :
455 : ! Make sure the correlation is in [-1,1]
456 0 : if ( calc_w_corr < -max_mag_correlation ) then
457 :
458 : calc_w_corr = -max_mag_correlation
459 :
460 0 : else if ( calc_w_corr > max_mag_correlation ) then
461 :
462 0 : calc_w_corr = max_mag_correlation
463 :
464 : end if
465 :
466 0 : end function calc_w_corr
467 :
468 :
469 : !-----------------------------------------------------------------------
470 0 : function calc_varnce( mixt_frac, x1, x2, xm, x1p2, x2p2 )
471 :
472 : ! Description:
473 : ! Calculate the variance xp2 from the components x1, x2.
474 :
475 : ! References:
476 : ! Larson et al. (2011), J. of Geophysical Research, Vol. 116, D00T02,
477 : ! page 3535
478 : !-----------------------------------------------------------------------
479 :
480 : use clubb_precision, only: &
481 : core_rknd ! Variable(s)
482 :
483 : implicit none
484 :
485 : ! Input Variables
486 : real( kind = core_rknd ), intent(in) :: &
487 : mixt_frac, & ! mixing ratio [-]
488 : x1, & ! first component of the double gaussian [units vary]
489 : x2, & ! second component of the double gaussian [units vary]
490 : xm, & ! mean of x [units vary]
491 : x1p2, & ! variance of the first component [units vary]
492 : x2p2 ! variance of the second component [units vary]
493 :
494 : ! Return Variable
495 : real( kind = core_rknd ) :: &
496 : calc_varnce ! variance of x (both components) [units vary]
497 :
498 : ! --- Begin Code ---
499 :
500 : calc_varnce &
501 : = mixt_frac * ( ( x1 - xm )**2 + x1p2 ) &
502 0 : + ( 1.0_core_rknd - mixt_frac ) * ( ( x2 - xm )**2 + x2p2 )
503 :
504 : return
505 : end function calc_varnce
506 :
507 : !-----------------------------------------------------------------------
508 0 : function calc_mean( mixt_frac, x1, x2 )
509 :
510 : ! Description:
511 : ! Calculate the mean xm from the components x1, x2.
512 :
513 : ! References:
514 : ! Larson et al. (2011), J. of Geophysical Research, Vol. 116, D00T02,
515 : ! page 3535
516 : !-----------------------------------------------------------------------
517 :
518 : use clubb_precision, only: &
519 : core_rknd ! Variable(s)
520 :
521 : implicit none
522 :
523 : ! Input Variables
524 : real( kind = core_rknd ), intent(in) :: &
525 : mixt_frac, & ! mixing ratio [-]
526 : x1, & ! first component of the double gaussian [units vary]
527 : x2 ! second component of the double gaussian [units vary]
528 :
529 : ! Return Variable
530 : real( kind = core_rknd ) :: &
531 : calc_mean ! mean of x (both components) [units vary]
532 :
533 : ! --- Begin Code ---
534 :
535 0 : calc_mean = mixt_frac * x1 + (1.0_core_rknd - mixt_frac) * x2
536 :
537 : return
538 : end function calc_mean
539 :
540 :
541 : !-----------------------------------------------------------------------
542 0 : subroutine calc_cholesky_corr_mtx_approx( n_variables, iiPDF_w, & ! Intent(in)
543 0 : corr_matrix, & ! intent(in)
544 0 : corr_cholesky_mtx, corr_mtx_approx ) ! intent(out)
545 :
546 : ! Description:
547 : ! This subroutine calculates the transposed correlation cholesky matrix
548 : ! from the correlation matrix
549 : !
550 : ! References:
551 : ! 1 Larson et al. (2011), J. of Geophysical Research, Vol. 116, D00T02
552 : ! 2 CLUBB Trac ticket#514
553 : !-----------------------------------------------------------------------
554 :
555 : use clubb_precision, only: &
556 : core_rknd ! Variable(s)
557 :
558 : use constants_clubb, only: &
559 : zero ! Variable(s)
560 :
561 : implicit none
562 :
563 : ! Input Variables
564 : integer, intent(in) :: &
565 : n_variables, & ! number of variables in the correlation matrix [-]
566 : iiPDF_w
567 :
568 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(in) :: &
569 : corr_matrix ! correlation matrix [-]
570 :
571 : ! Output Variables
572 :
573 : ! correlation cholesky matrix transposed L', C = LL'; see reference 1 formula 10
574 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(out) :: &
575 : corr_cholesky_mtx, & ! Transposed correlation cholesky matrix [-]
576 : corr_mtx_approx ! Approximated correlation matrix (C = LL') [-]
577 :
578 : ! Local Variables
579 : integer :: i, j ! Loop iterators
580 :
581 : ! Swapped means that the w-correlations are swapped to the first row
582 : real( kind = core_rknd ), dimension(n_variables,n_variables) :: &
583 0 : corr_cholesky_mtx_swap, & ! Swapped correlation cholesky matrix [-]
584 0 : corr_mtx_approx_swap, & ! Swapped correlation matrix (approx.) [-]
585 0 : corr_mtx_swap ! Swapped correlation matrix [-]
586 :
587 : !-------------------- Begin code --------------------
588 :
589 : call rearrange_corr_array( n_variables, iiPDF_w, & ! Intent(in)
590 : corr_matrix, & ! Intent(in)
591 0 : corr_mtx_swap ) ! Intent(inout)
592 :
593 : call setup_corr_cholesky_mtx( n_variables, corr_mtx_swap, & ! intent(in)
594 0 : corr_cholesky_mtx_swap ) ! intent(out)
595 :
596 : call rearrange_corr_array( n_variables, iiPDF_w, & ! Intent(in)
597 : corr_cholesky_mtx_swap, & ! Intent(in)
598 0 : corr_cholesky_mtx ) ! Intent(inout)
599 :
600 : call cholesky_to_corr_mtx_approx( n_variables, corr_cholesky_mtx_swap, & ! intent(in)
601 0 : corr_mtx_approx_swap ) ! intent(out)
602 :
603 : call rearrange_corr_array( n_variables, iiPDF_w, & ! Intent(in)
604 : corr_mtx_approx_swap, & ! Intent(in)
605 0 : corr_mtx_approx ) ! Intent(inout)
606 :
607 0 : call corr_array_assertion_checks( n_variables, corr_mtx_approx ) ! intent(in)
608 :
609 : ! Set lower triangle to zero for conformity
610 0 : do i = 2, n_variables
611 0 : do j = 1, i-1
612 0 : corr_mtx_approx(j,i) = zero
613 : end do
614 : end do
615 :
616 0 : return
617 :
618 : end subroutine calc_cholesky_corr_mtx_approx
619 : !-----------------------------------------------------------------------
620 :
621 : !-----------------------------------------------------------------------
622 0 : subroutine setup_corr_cholesky_mtx( n_variables, corr_matrix, & ! intent(in)
623 0 : corr_cholesky_mtx_t ) ! intent(out)
624 :
625 : ! Description:
626 : ! This subroutine calculates the transposed correlation cholesky matrix
627 : ! from the correlation matrix
628 : !
629 : ! References:
630 : ! 1 Larson et al. (2011), J. of Geophysical Research, Vol. 116, D00T02
631 : ! 2 CLUBB Trac ticket#514
632 : !-----------------------------------------------------------------------
633 :
634 : use clubb_precision, only: &
635 : core_rknd ! Variable(s)
636 :
637 : use constants_clubb, only: &
638 : zero, & ! Variable(s)
639 : one
640 :
641 : implicit none
642 :
643 : intrinsic :: sqrt
644 :
645 : ! Input Variables
646 : integer, intent(in) :: &
647 : n_variables ! number of variables in the correlation matrix [-]
648 :
649 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(in) :: &
650 : corr_matrix ! correlation matrix [-]
651 :
652 : ! Output Variables
653 :
654 : ! correlation cholesky matrix transposed L', C = LL'; see reference 1 formula 10
655 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(out) :: &
656 : corr_cholesky_mtx_t ! transposed correlation cholesky matrix [-]
657 :
658 : ! Local Variables
659 : integer :: i, j, k ! Loop iterators
660 :
661 : real( kind = core_rknd ), dimension(n_variables, n_variables) :: &
662 0 : s ! s(i,j) = sqrt(1-c(i,j)^2); see ref 1
663 :
664 : !-------------------- Begin code --------------------
665 :
666 : ! calculate all necessary square roots
667 0 : do i = 1, n_variables-1
668 0 : do j = i+1, n_variables
669 :
670 0 : s(j,i) = sqrt(1._core_rknd - corr_matrix(j,i)**2)
671 :
672 : end do
673 : end do
674 :
675 : !!! calculate transposed correlation cholesky matrix; ref 1 formula 10
676 :
677 : ! initialize matrix to zero
678 0 : do i = 1, n_variables
679 0 : do j = 1, n_variables
680 :
681 0 : corr_cholesky_mtx_t(j,i) = zero
682 :
683 : end do
684 : end do
685 :
686 : ! initialize upper triangle and diagonal to one
687 0 : do i = 1, n_variables
688 0 : do j = i, n_variables
689 :
690 0 : corr_cholesky_mtx_t(j,i) = one
691 :
692 : end do
693 : end do
694 :
695 : ! set diagonal elements
696 0 : do j = 2, n_variables
697 0 : do i = 1, j-1
698 :
699 0 : corr_cholesky_mtx_t(j,j) = corr_cholesky_mtx_t(j,j)*s(j,i)
700 : ! print *, "s(", j, ",", i, ") = ", s(j,i)
701 :
702 : end do
703 : end do
704 :
705 : ! set first row
706 0 : do j = 2, n_variables
707 :
708 0 : corr_cholesky_mtx_t(j,1) = corr_matrix(j,1)
709 :
710 : end do
711 :
712 : ! set upper triangle
713 0 : do i = 2, n_variables-1
714 0 : do j = i+1, n_variables
715 0 : do k = 1, i-1
716 :
717 0 : corr_cholesky_mtx_t(j,i) = corr_cholesky_mtx_t(j,i)*s(j,k)
718 :
719 : end do
720 :
721 0 : corr_cholesky_mtx_t(j,i) = corr_cholesky_mtx_t(j,i)*corr_matrix(j,i)
722 :
723 : end do
724 : end do
725 :
726 0 : return
727 :
728 : end subroutine setup_corr_cholesky_mtx
729 : !-----------------------------------------------------------------------
730 :
731 :
732 : !-----------------------------------------------------------------------
733 0 : subroutine cholesky_to_corr_mtx_approx( n_variables, corr_cholesky_mtx_t, & ! intent(in)
734 0 : corr_matrix_approx ) ! intent(out)
735 :
736 : ! Description:
737 : ! This subroutine approximates the correlation matrix from the correlation
738 : ! cholesky matrix
739 : !
740 : ! References:
741 : ! 1 Larson et al. (2011), J. of Geophysical Research, Vol. 116, D00T02
742 : ! 2 CLUBB Trac ticket#514
743 : !-----------------------------------------------------------------------
744 :
745 : use clubb_precision, only: &
746 : core_rknd ! Variable(s)
747 :
748 : implicit none
749 :
750 : intrinsic :: matmul, transpose
751 :
752 : ! Input Variables
753 : integer, intent(in) :: &
754 : n_variables ! number of variables in the correlation matrix [-]
755 :
756 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(in) :: &
757 : corr_cholesky_mtx_t ! transposed correlation cholesky matrix [-]
758 :
759 : ! Output Variables
760 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(out) :: &
761 : corr_matrix_approx ! correlation matrix [-]
762 :
763 : !-------------------- Begin code --------------------
764 :
765 : ! approximate the correlation matrix; see ref 1 formula (8)
766 0 : corr_matrix_approx = matmul(corr_cholesky_mtx_t, transpose(corr_cholesky_mtx_t))
767 :
768 0 : return
769 :
770 : end subroutine cholesky_to_corr_mtx_approx
771 : !-----------------------------------------------------------------------
772 :
773 :
774 : !-----------------------------------------------------------------------
775 0 : subroutine corr_array_assertion_checks( n_variables, corr_array )
776 :
777 : ! Description:
778 : ! This subroutine does the assertion checks for the corr_array.
779 :
780 : ! References:
781 : !
782 : !
783 : !-----------------------------------------------------------------------
784 :
785 : use clubb_precision, only: &
786 : core_rknd ! Variable(s)
787 :
788 : use constants_clubb, only: &
789 : max_mag_correlation ! Variable(s)
790 :
791 : use constants_clubb, only: &
792 : one ! Variable(s)
793 :
794 : use error_code, only: &
795 : clubb_at_least_debug_level ! Procedure
796 :
797 : implicit none
798 :
799 : ! Input Variables
800 : integer, intent(in) :: &
801 : n_variables ! number of variables in the correlation matrix [-]
802 :
803 : real( kind = core_rknd ), dimension(n_variables,n_variables), intent(in) :: &
804 : corr_array ! correlation matrix [-]
805 :
806 : ! Local Variables
807 : integer :: i, j ! Loop iterator
808 :
809 : real( kind = core_rknd ), parameter :: &
810 : tol = 1.e-6_core_rknd ! Maximum acceptable tolerance for the difference of the diagonal
811 : ! elements of corr_array to one
812 :
813 : !-------------------- Begin code --------------------
814 :
815 0 : if ( clubb_at_least_debug_level( 1 ) ) then
816 :
817 0 : do i = 1, n_variables - 1
818 0 : do j = i+1, n_variables
819 :
820 : ! Check if upper and lower triangle values are within the correlation boundaries
821 0 : if ( ( corr_array(i,j) < -max_mag_correlation ) &
822 : .or. ( corr_array(i,j) > max_mag_correlation ) &
823 : .or. ( corr_array(j,i) < -max_mag_correlation ) &
824 0 : .or. ( corr_array(j,i) > max_mag_correlation ) ) &
825 0 : then
826 :
827 0 : error stop "Error: A value in the correlation matrix is out of range."
828 :
829 : endif
830 :
831 : enddo
832 : enddo
833 :
834 : endif
835 :
836 0 : if ( clubb_at_least_debug_level( 2 ) ) then
837 :
838 0 : do i = 1, n_variables
839 : ! Check if the diagonal elements are one (up to a tolerance)
840 0 : if ( ( corr_array(i,i) > one + tol ) .or. (corr_array(i,i) < one - tol ) ) then
841 :
842 0 : error stop "Error: Diagonal element(s) of the correlation matrix are unequal to one."
843 :
844 : endif
845 : enddo
846 :
847 : endif
848 :
849 0 : return
850 :
851 : end subroutine corr_array_assertion_checks
852 :
853 :
854 : !-----------------------------------------------------------------------
855 0 : subroutine rearrange_corr_array( pdf_dim, iiPDF_w, & ! Intent(in)
856 0 : corr_array, & ! Intent(in)
857 0 : corr_array_swapped) ! Intent(out)
858 : ! Description:
859 : ! This subroutine swaps the w-correlations to the first row if the input
860 : ! matrix is in the same order as the *_corr_array_cloud.in files. It swaps
861 : ! the rows back to the order of the *_corr_array_cloud.in files if the
862 : ! input matrix is already swapped (first row w-correlations).
863 : !
864 : ! References:
865 : !
866 : !-----------------------------------------------------------------------
867 :
868 : use clubb_precision, only: &
869 : core_rknd ! Variable(s)
870 :
871 : implicit none
872 :
873 : intrinsic :: max, sqrt, transpose
874 :
875 : ! Input Variables
876 : integer, intent(in) :: &
877 : pdf_dim, & ! number of diagnosed correlations
878 : iiPDF_w
879 :
880 : real( kind = core_rknd ), dimension(pdf_dim, pdf_dim), intent(in) :: &
881 : corr_array ! Correlation matrix
882 :
883 : ! Output variables
884 : real( kind = core_rknd ), dimension(pdf_dim, pdf_dim), intent(out) :: &
885 : corr_array_swapped ! Swapped correlation matrix
886 :
887 : ! Local Variables
888 : real( kind = core_rknd ), dimension(pdf_dim) :: &
889 0 : swap_array
890 :
891 : !-------------------- Begin code --------------------
892 :
893 : ! Swap the w-correlations to the first row for the prescribed correlations
894 0 : corr_array_swapped = corr_array
895 0 : swap_array = corr_array_swapped (:,1)
896 0 : corr_array_swapped(1:iiPDF_w, 1) = corr_array_swapped(iiPDF_w, iiPDF_w:1:-1)
897 0 : corr_array_swapped((iiPDF_w+1):pdf_dim, 1) = corr_array_swapped( &
898 0 : (iiPDF_w+1):pdf_dim, iiPDF_w)
899 0 : corr_array_swapped(iiPDF_w, 1:iiPDF_w) = swap_array(iiPDF_w:1:-1)
900 0 : corr_array_swapped((iiPDF_w+1):pdf_dim, iiPDF_w) = swap_array((iiPDF_w+1):pdf_dim)
901 :
902 0 : return
903 :
904 : end subroutine rearrange_corr_array
905 : !-----------------------------------------------------------------------
906 :
907 :
908 : !-----------------------------------------------------------------------
909 : ! subroutine set_w_corr( nz, pdf_dim, & ! Intent(in)
910 : ! corr_chi_w, corr_wrr, corr_wNr, corr_wNcn, &
911 : ! corr_array ) ! Intent(inout)
912 : !
913 : ! ! Description:
914 : ! ! Set the first row of corr_array to the according w-correlations.
915 : !
916 : ! ! References:
917 : ! ! clubb:ticket:514
918 : ! !-----------------------------------------------------------------------
919 : !
920 : ! use clubb_precision, only: &
921 : ! core_rknd ! Variable(s)
922 : !
923 : ! use array_index, only: &
924 : ! iiPDF_w, & ! Variable(s)
925 : ! iiPDF_chi, &
926 : ! iiPDF_rr, &
927 : ! iiPDF_Nr, &
928 : ! iiPDF_Ncn
929 : !
930 : ! implicit none
931 : !
932 : ! ! Input Variables
933 : ! integer, intent(in) :: &
934 : ! nz, & ! Number of model vertical grid levels
935 : ! pdf_dim ! Number of Variables to be diagnosed
936 : !
937 : ! real( kind = core_rknd ), dimension(nz), intent(in) :: &
938 : ! corr_chi_w, & ! Correlation between chi (s) & w (both components) [-]
939 : ! corr_wrr, & ! Correlation between rr & w (both components) [-]
940 : ! corr_wNr, & ! Correlation between Nr & w (both components) [-]
941 : ! corr_wNcn ! Correlation between Ncn & w (both components) [-]
942 : !
943 : ! ! Input/Output Variables
944 : ! real( kind = core_rknd ), dimension(pdf_dim, pdf_dim, nz), &
945 : ! intent(inout) :: &
946 : ! corr_array
947 : !
948 : ! ! ----- Begin Code -----
949 : !
950 : ! corr_array(iiPDF_w, iiPDF_chi, :) = corr_chi_w
951 : ! corr_array(iiPDF_w, iiPDF_rr, :) = corr_wrr
952 : ! corr_array(iiPDF_w, iiPDF_Nr, :) = corr_wNr
953 : ! corr_array(iiPDF_w, iiPDF_Ncn, :) = corr_wNcn
954 : !
955 : ! end subroutine set_w_corr
956 :
957 : !=============================================================================
958 : ! subroutine unpack_correlations( pdf_dim, corr_array, & ! Intent(in)
959 : ! corr_w_chi, corr_wrr, corr_wNr, corr_wNcn, &
960 : ! corr_chi_eta, corr_chi_rr, corr_chi_Nr, corr_chi_Ncn, &
961 : ! corr_eta_rr, corr_eta_Nr, corr_eta_Ncn, corr_rrNr )
962 : !
963 : ! ! Description:
964 : !
965 : ! ! References:
966 : ! !-----------------------------------------------------------------------
967 :
968 : ! use clubb_precision, only: &
969 : ! core_rknd ! Variable(s)
970 :
971 : ! use array_index, only: &
972 : ! iiPDF_w, & ! Variable(s)
973 : ! iiPDF_chi, &
974 : ! iiPDF_eta, &
975 : ! iiPDF_rr, &
976 : ! iiPDF_Nr, &
977 : ! iiPDF_Ncn
978 :
979 : ! implicit none
980 :
981 : ! intrinsic :: max, sqrt, transpose
982 :
983 : ! ! Input Variables
984 : ! integer, intent(in) :: &
985 : ! pdf_dim ! number of diagnosed correlations
986 :
987 : ! real( kind = core_rknd ), dimension(pdf_dim, pdf_dim), intent(in) :: &
988 : ! corr_array ! Prescribed correlations
989 :
990 : ! ! Output variables
991 : ! real( kind = core_rknd ), intent(out) :: &
992 : ! corr_w_chi, & ! Correlation between w and chi(s) (1st PDF component) [-]
993 : ! corr_wrr, & ! Correlation between w and rr (1st PDF component) ip [-]
994 : ! corr_wNr, & ! Correlation between w and Nr (1st PDF component) ip [-]
995 : ! corr_wNcn, & ! Correlation between w and Ncn (1st PDF component) [-]
996 : ! corr_chi_eta, & ! Correlation between chi(s) and eta(t) (1st PDF component) [-]
997 : ! corr_chi_rr, & ! Correlation between chi(s) and rr (1st PDF component) ip [-]
998 : ! corr_chi_Nr, & ! Correlation between chi(s) and Nr (1st PDF component) ip [-]
999 : ! corr_chi_Ncn, & ! Correlation between chi(s) and Ncn (1st PDF component) [-]
1000 : ! corr_eta_rr, & ! Correlation between eta(t) and rr (1st PDF component) ip [-]
1001 : ! corr_eta_Nr, & ! Correlation between eta(t) and Nr (1st PDF component) ip [-]
1002 : ! corr_eta_Ncn, & ! Correlation between (t) and Ncn (1st PDF component) [-]
1003 : ! corr_rrNr ! Correlation between rr & Nr (1st PDF component) ip [-]
1004 :
1005 : ! ---- Begin Code ----
1006 :
1007 : ! corr_w_chi = corr_array(iiPDF_w, iiPDF_chi)
1008 : ! corr_wrr = corr_array(iiPDF_w, iiPDF_rr)
1009 : ! corr_wNr = corr_array(iiPDF_w, iiPDF_Nr)
1010 : ! corr_wNcn = corr_array(iiPDF_w, iiPDF_Ncn)
1011 : ! corr_chi_eta = corr_array(iiPDF_chi, iiPDF_eta)
1012 : ! corr_chi_rr = corr_array(iiPDF_chi, iiPDF_rr)
1013 : ! corr_chi_Nr = corr_array(iiPDF_chi, iiPDF_Nr)
1014 : ! corr_chi_Ncn = corr_array(iiPDF_chi, iiPDF_Ncn)
1015 : ! corr_eta_rr = corr_array(iiPDF_eta, iiPDF_rr)
1016 : ! corr_eta_Nr = corr_array(iiPDF_eta, iiPDF_Nr)
1017 : ! corr_eta_Ncn = corr_array(iiPDF_eta, iiPDF_Ncn)
1018 : ! corr_rrNr = corr_array(iiPDF_rr, iiPDF_Nr)
1019 :
1020 : ! corr_w_chi = corr_array(iiPDF_chi, iiPDF_w)
1021 : ! corr_wrr = corr_array(iiPDF_rr, iiPDF_w)
1022 : ! corr_wNr = corr_array(iiPDF_Nr, iiPDF_w)
1023 : ! corr_wNcn = corr_array(iiPDF_Ncn, iiPDF_w)
1024 : ! corr_chi_eta = corr_array(iiPDF_eta, iiPDF_chi)
1025 : ! corr_chi_rr = corr_array(iiPDF_rr, iiPDF_chi)
1026 : ! corr_chi_Nr = corr_array(iiPDF_Nr, iiPDF_chi)
1027 : ! corr_chi_Ncn = corr_array(iiPDF_Ncn, iiPDF_chi)
1028 : ! corr_eta_rr = corr_array(iiPDF_rr, iiPDF_eta)
1029 : ! corr_eta_Nr = corr_array(iiPDF_Nr, iiPDF_eta)
1030 : ! corr_eta_Ncn = corr_array(iiPDF_Ncn, iiPDF_eta)
1031 : ! corr_rrNr = corr_array(iiPDF_rr, iiPDF_Nr)
1032 :
1033 : ! end subroutine unpack_correlations
1034 :
1035 : !===============================================================================
1036 :
1037 : end module diagnose_correlations_module
|