| blob | e16e0c79324507ff24d09950fcf7b9938c69d814 |
1 subroutine da_setup_be_regional(xb, be, grid)
3 !---------------------------------------------------------------------------
4 ! Purpose: Define and allocate components of background errors
5 !
6 ! Updates:
7 !
8 ! Implementation of multi-variate BE
9 ! Syed RH Rizvi, MMM/NESL/NCAR, Date: 02/01/2010
10 !---------------------------------------------------------------------------
12 #if (WRF_CHEM == 1)
13 USE module_state_description, ONLY : p_bc1, p_bc2, p_oc1, p_oc2 & !inline gocart
14 ,p_dust_1, p_dust_2, p_dust_3 ,p_dust_4 &
15 ,p_seas_1, p_seas_2, p_seas_3 ,p_seas_4 &
16 ,p_co, p_o3, p_no2, p_so2 &
17 ,p_p25, p_p10, p_sulf
19 USE module_state_description, ONLY : p_bc_a01, p_bc_a02, p_bc_a03, p_bc_a04 & !inline mosaic
20 ,p_oc_a01, p_oc_a02, p_oc_a03, p_oc_a04 &
21 ,p_so4_a01, p_so4_a02, p_so4_a03, p_so4_a04 &
22 ,p_no3_a01, p_no3_a02, p_no3_a03, p_no3_a04 &
23 ,p_nh4_a01, p_nh4_a02, p_nh4_a03, p_nh4_a04 &
24 ,p_cl_a01, p_cl_a02, p_cl_a03, p_cl_a04 &
25 ,p_na_a01, p_na_a02, p_na_a03, p_na_a04 &
26 ,p_oin_a01, p_oin_a02, p_oin_a03, p_oin_a04
28 USE module_state_description, ONLY : p_so4aj, p_so4ai, p_nh4aj, p_nh4ai, & ! racm_soa_vbs_da
29 p_no3aj, p_no3ai, p_naaj, p_naai, &
30 p_asoa1j, p_asoa1i, p_asoa2j, p_asoa2i, &
31 p_asoa3j, p_asoa3i, p_asoa4j, p_asoa4i, &
32 p_bsoa1j, p_bsoa1i, p_bsoa2j, p_bsoa2i, &
33 p_bsoa3j, p_bsoa3i, p_bsoa4j, p_bsoa4i, &
34 p_orgpaj, p_orgpai, p_ecj, p_eci, &
35 p_p25j, p_p25i, p_antha, p_seas, &
36 p_claj, p_clai, p_soila
37 ! p_so2, p_no2, p_o3, p_co ! gas species
39 #endif
41 implicit none
43 type (xb_type), intent(in) :: xb ! First guess structure.
44 type (be_type), intent(inout) :: be ! Back. errors structure.
45 type (domain), intent(in) :: grid
47 logical do_normalize1 ! Test do_normalize.
48 integer :: i, j, k, m ! Loop counters.
49 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point
50 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point
51 integer :: bin_type ! Type of bin to average over.
52 integer :: num_bins ! Number of bins (3D fields).
53 integer :: num_bins2d ! Number of bins (3D fields).
54 real*8 :: junk ! For reading truncated sd, wsd
55 real*8 :: lat_min, lat_max, binwidth_lat ! Used if bin_type = 2 (degrees)..
56 real*8 :: hgt_min, hgt_max, binwidth_hgt ! Used if bin_type = 2 (m). .
58 integer :: num_cv_3d_basic ! number of 3D mandatory control variables
59 integer :: num_cv_3d_extra ! number of 3D optional control variables
60 integer :: num_cv_2d ! number of 2D control variables
62 real*8, allocatable :: be1_eval_loc(:,:) ! Temp arrays.
63 real*8, allocatable :: be2_eval_loc(:,:) ! Temp arrays.
64 real*8, allocatable :: be3_eval_loc(:,:) ! Temp arrays.
65 real*8, allocatable :: be4_eval_loc(:,:) ! Temp arrays.
66 real*8, allocatable :: be5_eval_loc(:,:) ! Temp arrays.
68 real*8, allocatable :: be6_eval_loc(:,:) ! Temp arrays.
69 real*8, allocatable :: be7_eval_loc(:,:) ! Temp arrays.
70 real*8, allocatable :: be8_eval_loc(:,:) ! Temp arrays.
71 real*8, allocatable :: be9_eval_loc(:,:) ! Temp arrays.
72 real*8, allocatable :: be10_eval_loc(:,:) ! Temp arrays.
73 real*8, allocatable :: be11_eval_loc(:,:) ! Temp arrays.
75 real*8, allocatable :: be1_eval_glo(:) ! Global Eigenvalues.
76 real*8, allocatable :: be2_eval_glo(:) ! Global Eigenvalues.
77 real*8, allocatable :: be3_eval_glo(:) ! Global Eigenvalues.
78 real*8, allocatable :: be4_eval_glo(:) ! Global Eigenvalues.
79 real*8, allocatable :: be5_eval_glo(:) ! Global Eigenvalues.
81 real*8, allocatable :: be6_eval_glo(:) ! Global Eigenvalues.
82 real*8, allocatable :: be7_eval_glo(:) ! Global Eigenvalues.
83 real*8, allocatable :: be8_eval_glo(:) ! Global Eigenvalues.
84 real*8, allocatable :: be9_eval_glo(:) ! Global Eigenvalues.
85 real*8, allocatable :: be10_eval_glo(:) ! Global Eigenvalues.
86 real*8, allocatable :: be11_eval_glo(:) ! Global Eigenvalues.
88 real*8, allocatable :: alpha_val(:) ! Global Eigenvalues.
90 real*8, allocatable :: be1_evec_loc(:,:,:) ! Local Eigenvectors.
91 real*8, allocatable :: be2_evec_loc(:,:,:) ! Local Eigenvectors.
92 real*8, allocatable :: be3_evec_loc(:,:,:) ! Local Eigenvectors.
93 real*8, allocatable :: be4_evec_loc(:,:,:) ! Local Eigenvectors.
94 real*8, allocatable :: be5_evec_loc(:,:,:) ! Local Eigenvectors.
96 real*8, allocatable :: be6_evec_loc(:,:,:) ! Local Eigenvectors.
97 real*8, allocatable :: be7_evec_loc(:,:,:) ! Local Eigenvectors.
98 real*8, allocatable :: be8_evec_loc(:,:,:) ! Local Eigenvectors.
99 real*8, allocatable :: be9_evec_loc(:,:,:) ! Local Eigenvectors.
100 real*8, allocatable :: be10_evec_loc(:,:,:) ! Local Eigenvectors.
101 real*8, allocatable :: be11_evec_loc(:,:,:) ! Local Eigenvectors.
103 real*8, allocatable :: be1_evec_glo(:,:) ! Global Eigenvectors.
104 real*8, allocatable :: be2_evec_glo(:,:) ! Global Eigenvectors.
105 real*8, allocatable :: be3_evec_glo(:,:) ! Global Eigenvectors.
106 real*8, allocatable :: be4_evec_glo(:,:) ! Global Eigenvectors.
107 real*8, allocatable :: be5_evec_glo(:,:) ! Global Eigenvectors.
109 real*8, allocatable :: be6_evec_glo(:,:) ! Global Eigenvectors.
110 real*8, allocatable :: be7_evec_glo(:,:) ! Global Eigenvectors.
111 real*8, allocatable :: be8_evec_glo(:,:) ! Global Eigenvectors.
112 real*8, allocatable :: be9_evec_glo(:,:) ! Global Eigenvectors.
113 real*8, allocatable :: be10_evec_glo(:,:) ! Global Eigenvectors.
114 real*8, allocatable :: be11_evec_glo(:,:) ! Global Eigenvectors.
117 real*8, allocatable :: alpha_evec(:,:) ! Global Eigenvectors.
119 real*8, allocatable :: be1_rf_lengthscale(:) ! RF lengthscale.
120 real*8, allocatable :: be2_rf_lengthscale(:) ! RF lengthscale.
121 real*8, allocatable :: be3_rf_lengthscale(:) ! RF lengthscale.
122 real*8, allocatable :: be4_rf_lengthscale(:) ! RF lengthscale.
123 real*8, allocatable :: be5_rf_lengthscale(:)
125 real*8, allocatable :: be6_rf_lengthscale(:)
126 real*8, allocatable :: be7_rf_lengthscale(:)
127 real*8, allocatable :: be8_rf_lengthscale(:)
128 real*8, allocatable :: be9_rf_lengthscale(:)
129 real*8, allocatable :: be10_rf_lengthscale(:)
130 real*8, allocatable :: be11_rf_lengthscale(:)
132 real*8, allocatable :: alpha_rf_lengthscale(:)
133 real*8, allocatable :: alpha_rf_scale_factor(:)
135 real*8, allocatable :: evec_loc(:,:,:) ! Latitudinally varying eigenvectors.
136 real*8, allocatable :: eval_loc(:,:) ! Latitudinally varying eigenvalues.
138 #if (WRF_CHEM == 1)
139 integer :: ispcc ! chem species index.
140 integer :: num_cv_3d_chem ! number of 3D chem control variables
141 real*8, allocatable :: be12_eval_loc(:,:,:) ! Temp arrays for chem.
142 real*8, allocatable :: be12_eval_glo(:,:) ! Global Eigenvalues for chem.
143 real*8, allocatable :: be12_evec_loc(:,:,:,:) ! Local Eigenvectors for chem.
144 real*8, allocatable :: be12_evec_glo(:,:,:) ! Global Eigenvectors for chem.
145 real*8, allocatable :: be12_rf_lengthscale(:,:) ! for chem
146 #endif
148 character*5 :: namv(5)=(/"psi ","chi_u","t_u ","rh ","ps_u "/)
149 character*10 :: variable
150 character*80 :: var80c
151 character(len=15) :: write_fmt
152 integer :: ni, nj, nk, nk_2d, b
153 integer :: ix, jy, kz, kzs(5), mz, mzs(5)
154 real, allocatable :: regcoeff_psi_chi(:) ! psi/chi regression cooefficient.
155 real, allocatable :: regcoeff_psi_t(:,:,:) ! psi/t regression cooefficient.
156 real, allocatable :: regcoeff_psi_ps(:,:) ! psi/ps regression cooefficient.
157 real, allocatable :: regcoeff_psi_rh(:,:,:) ! psi/rh regression cooefficient.
158 real, allocatable :: regcoeff_chi_u_t(:,:,:) ! chi_u/t regression coefficient
159 real, allocatable :: regcoeff_chi_u_ps(:,:) ! chi_u/ps regression coefficient
160 real, allocatable :: regcoeff_chi_u_rh(:,:,:) ! chi_u/rh regression coefficient
161 real, allocatable :: regcoeff_t_u_rh(:,:,:) ! t_u/rh regression coefficient
162 real, allocatable :: regcoeff_ps_u_rh(:,:) ! ps_u/rh regression coefficient
163 !real :: qrain_th_low, qrain_th_high
165 integer :: be_unit, ier, be_rf_unit, be_print_unit, it, idummy
167 !-----------for interpolating CV5--------------------------------------------------------------
168 REAL, ALLOCATABLE :: reg_psi_ps0(:,:), reg_psi_chi0(:,:), reg_psi_t0(:,:,:), &
169 reg_psi_ps (:), reg_psi_chi (:), reg_psi_t (:,:), &
170 reg_psi_ps_be(:), reg_psi_chi_be(:), reg_psi_t_be(:,:)
171 REAL, DIMENSION(:,:), ALLOCATABLE:: covm1_be, covm2_be, covm3_be, covm4_be,&
172 covm1 , covm2 , covm3 , covm4
173 REAL, DIMENSION(:), ALLOCATABLE:: rfls1_be, rfls2_be, rfls3_be, rfls4_be
174 REAL, DIMENSION(:), ALLOCATABLE:: eta_be
175 ! For write out the Interpolated CV5 BE (bin_type=5) file:
176 integer :: bin_type_out, num_bins_out, num_bins2d_out
177 integer, allocatable :: bin_out(:,:,:), bin2d_out(:,:)
178 integer, parameter :: be_out_unit = 135
180 !===========for interpolating CV5==============================================================
182 if (trace_use) call da_trace_entry("da_setup_be_regional")
184 ix = xb % mix
185 jy = xb % mjy
186 kz = xb % mkz
188 if ( jb_factor > 0.0 ) then
189 if( (cv_options == 5) .or. (cv_options == 6) .or. (cv_options == 7) ) then
190 write (unit=message(1),fmt='(A,I3)') &
191 'Set up background errors for regional application for cv_options = ',cv_options
192 write (unit=message(2),fmt='(a)') ' '
193 call da_message(message(1:2))
194 if ( (cv_options == 7 .or. cv_options == 6) .and. interpolate_stats ) then
195 write (unit=message(1),fmt='(A,I3)') 'interpolate_stats is not implemented for cv_options=',cv_options
196 call da_error(__FILE__,__LINE__,message(1:1))
197 end if
198 if ( ( cloud_cv_options >=2 ) .and. interpolate_stats ) then
199 write (unit=message(1),fmt='(A,I3)') 'interpolate_stats is not implemented for cloud_cv_options=',cloud_cv_options
200 call da_error(__FILE__,__LINE__,message(1:1))
201 end if
202 #if (WRF_CHEM == 1)
203 if ( ( chem_cv_options >=10 ) .and. interpolate_stats ) then
204 write (unit=message(1),fmt='(A,I3)') 'interpolate_stats is not implemented for chem_cv_options=',chem_cv_options
205 call da_error(__FILE__,__LINE__,message(1:1))
206 end if
207 #endif
208 else
209 write (unit=message(1),fmt='(a)') 'cv_options should be 5, 6, or 7'
210 call da_error(__FILE__,__LINE__,message(1:1))
211 end if
213 !be_rf_unit = unit_end + 1
214 !be_print_unit = unit_end + 2
216 call da_get_unit(be_unit)
217 call da_get_unit(be_rf_unit)
218 call da_get_unit(be_print_unit)
219 open(unit=be_unit,file="be.dat", status="old",form="unformatted")
220 open(unit=be_rf_unit ,file="be_rf.dat" , status="unknown",form="unformatted")
221 open(unit=be_print_unit,file="be_print.dat", status="unknown")
223 rewind (be_unit)
224 read (be_unit, iostat=ier) ni, nj, nk
225 print *, 'ni, nj, nk = ', ni, nj, nk
226 if (ier /= 0) then
227 write (unit=message(1),fmt='(a,i3)') 'Error in reading be.dat, unit= ',be_unit
228 call da_error(__FILE__,__LINE__,message(1:1))
229 end if
230 read (be_unit) bin_type
231 print *, 'bin_type = ', bin_type
233 !-----------for interpolating CV5--------------------------------------------------------------
234 if ( .not. interpolate_stats ) then
235 if (ni /= ix .or. nj /= jy .or. nk /= kz) then
236 write(unit=message(1),fmt='(a)') &
237 'Dimensions of the assimilation domain do not match those in the BE file.'
238 write(unit=message(2),fmt='(3x,a,3i4)') "in be.dat: ni, nj, nk = ",ni, nj, nk
239 write(unit=message(3),fmt='(3x,a,3i4)') "in fg: ix, jy, kz = ",ix, jy, kz
240 write(message(4),'(a)') 'be.dat need to be generated using data consistent with the assimilation domain.'
241 call da_error(__FILE__,__LINE__,message(1:4))
242 end if
243 else ! when interpolate_stats = true
244 ! Must use the domain averaged Reg. Coeff.:
245 lat_stats_option = .FALSE.
246 ! Must use the global Eigenvector/eigenvalue
247 vert_evalue = 1
248 write (unit=message(1),fmt='(A,I3,A)') 'interpolate_stats is .true., BE for cv_options = ',cv_options,' will be interpolated.'
249 call da_message(message(1:1))
250 if (rootproc) write(be_out_unit) ix, jy, kz
251 write(6,'(a/("k=",i3,2x,f10.4))') 'xb%sigmah=',(k,xb%sigmah(k),k=1,kz)
252 allocate (eta_be(1:nk))
253 read(be_unit,iostat=ier) eta_be
254 ! If not available from BE file, get it from namelist
255 if (ier /= 0) then
256 write(message(1),'("ier=",i5,2x,a)') ier, " NO Eta values existed in BE file."
257 write(message(2),'(a)') "Will use BE Eta values from the namelist be_eta."
258 call da_message(message(1:2))
259 rewind(be_unit)
260 read(be_unit) idummy, idummy, idummy
261 read(be_unit) idummy
262 eta_be = be_eta(1:nk)
263 endif
264 write(6,'(a/("k=",i3,2x,f10.4))') 'eta_be=',(k,eta_be(k),k=1,nk)
265 ! If eta values not defined correctly, stop
266 if (eta_be(1) == eta_be(nk)) then
267 write(message(1),'(a)') 'Wrong eta values found.'
268 write(message(2),'(a)') 'The eta levels used in generating be.dat should be specified in namelist &wrfvar10 be_eta.'
269 call da_error(__FILE__,__LINE__,message(1:2))
270 endif
271 end if
273 !===========for interpolating CV5==============================================================
275 allocate (bin(1:ni,1:nj,1:nk))
276 allocate (bin2d(1:ni,1:nj))
278 !if(cloud_cv_options.eq.2)then
279 ! read (be_unit)num_bins, num_bins2d
280 ! read (be_unit)lat_min, lat_max, binwidth_lat
281 ! read (be_unit)hgt_min, hgt_max, binwidth_hgt
282 ! read (be_unit)qrain_th_low, qrain_th_high
283 ! read (be_unit)bin(1:ni,1:nj,1:nk)
284 ! read (be_unit)bin2d(1:ni,1:nj)
285 !else
286 read (be_unit)lat_min, lat_max, binwidth_lat
287 read (be_unit)hgt_min, hgt_max, binwidth_hgt
288 read (be_unit)num_bins, num_bins2d
289 read (be_unit)bin(1:ni,1:nj,1:nk)
290 read (be_unit)bin2d(1:ni,1:nj)
291 !end if
293 print *, lat_min, lat_max, binwidth_lat
294 print *, hgt_min, hgt_max, binwidth_hgt
295 print *, 'num_bins, num_bins2d = ', num_bins, num_bins2d
296 print *, 'bin = ', bin(1:1,1:1,1:1)
297 print *, 'bin2d = ', bin2d(1:1,1:1)
299 num_cv_3d_basic = 4 ! (Psi, Chi_u, T_u, Pseudo RH) or (U,V,T,Pseudo RH)
300 num_cv_3d_extra = 0 ! default: no cloud and w
301 num_cv_2d = 1 ! Ps_u
302 if ( cloud_cv_options >= 2 ) then
303 num_cv_3d_extra = 5 ! 5 hydrometeors, qc,qr,qi,qs,qg
304 end if
305 #if (WRF_CHEM == 1)
306 select case ( chem_cv_options )
307 case ( 10 )
308 num_cv_3d_chem = 19
309 case ( 20 )
310 num_cv_3d_chem = 36
311 case ( 108 ) ! racm_soa_vbs_da
312 num_cv_3d_chem = 39
313 case default ! mozcart_da
314 num_cv_3d_chem = 19
315 end select
316 if ( chem_cv_options >= 1 ) then
317 allocate (be % v12 (num_cv_3d_chem) )
318 num_cv_3d_extra = num_cv_3d_extra + num_cv_3d_chem
319 endif ! ( chem_cv_options >= 1 )
320 #endif
321 if ( use_cv_w ) then
322 num_cv_3d_extra = num_cv_3d_extra + 1 ! clouds+w
323 end if
325 !-----------for interpolating CV5--------------------------------------------------------------
326 ! Write out the interpolated CV5 BE:
327 if (interpolate_stats) then
329 ! No matter which bin_type BE inputed, after Interpolated, the BE always
330 ! became bin_type = 1 BE, i.e. the global BE.
331 bin_type_out = 5
332 num_bins_out = kz
333 num_bins2d_out = 1
334 if (rootproc) then
335 write(be_out_unit)bin_type_out
336 write(be_out_unit)lat_min, lat_max, binwidth_lat
337 write(be_out_unit)hgt_min, hgt_max, binwidth_hgt
338 write(be_out_unit)num_bins_out, num_bins2d_out
339 endif
340 allocate (bin_out(1:ix,1:jy,1:kz))
341 allocate (bin2d_out(1:ix,1:jy))
342 do k = 1, kz
343 bin_out (:,:,k) = k
344 enddo
345 bin2d_out = 1
346 if (rootproc) then
347 write(be_out_unit)bin_out
348 write(be_out_unit)bin2d_out
349 endif
350 deallocate (bin_out)
351 deallocate (bin2d_out)
352 endif
353 !===========for interpolating CV5==============================================================
354 !
355 if ( cv_options /= 7 ) then ! No regression coefficients for cv_options == 7
357 ! 1.1 Read in regression coefficients
358 allocate (regcoeff_psi_chi(1:num_bins))
359 allocate (regcoeff_psi_t(1:nk,1:nk,1:num_bins2d))
360 allocate (regcoeff_psi_ps(1:nk,1:num_bins2d))
361 if ( cv_options == 6 ) then
362 allocate (regcoeff_psi_rh(1:nk,1:nk,1:num_bins2d))
363 allocate (regcoeff_chi_u_t(1:nk,1:nk,1:num_bins2d))
364 allocate (regcoeff_chi_u_ps(1:nk,1:num_bins2d))
365 allocate (regcoeff_chi_u_rh(1:nk,1:nk,1:num_bins2d))
366 allocate (regcoeff_t_u_rh(1:nk,1:nk,1:num_bins2d))
367 allocate (regcoeff_ps_u_rh(1:nk,1:num_bins2d))
368 end if
370 if ( cv_options == 5 ) then
371 read (be_unit) regcoeff_psi_chi
372 read (be_unit) regcoeff_psi_ps
373 read (be_unit) regcoeff_psi_t
374 else
375 do i = 1 , 9
376 read (be_unit) var80c
377 select case( trim(adjustl(var80c)) )
378 case ('regcoeff_psi_chi')
379 read (be_unit) regcoeff_psi_chi
380 case ('regcoeff_psi_t')
381 read (be_unit) regcoeff_psi_t
382 case ('regcoeff_psi_ps')
383 read (be_unit) regcoeff_psi_ps
384 case ('regcoeff_psi_rh')
385 read (be_unit) regcoeff_psi_rh
386 case ('regcoeff_chi_u_t')
387 read (be_unit) regcoeff_chi_u_t
388 case ('regcoeff_chi_u_ps')
389 read (be_unit) regcoeff_chi_u_ps
390 case ('regcoeff_chi_u_rh')
391 read (be_unit) regcoeff_chi_u_rh
392 case ('regcoeff_t_u_rh')
393 read (be_unit) regcoeff_t_u_rh
394 case ('regcoeff_ps_u_rh')
395 read (be_unit) regcoeff_ps_u_rh
396 case default;
397 message(1)=' Read problem in regression coefficients in BE file '
398 write (unit=message(2),fmt='(A,A)') ' Trying to read regression coefficients for variable: ',trim(adjustl(var80c))
399 call da_error(__FILE__,__LINE__,message(1:2))
400 end select
401 end do
402 end if
403 ! 1.2 Fill regression coeff. array for model BE, vertical dimension is kz:
405 allocate (be%reg_psi_chi (1:jy,1:kz))
406 allocate (be%reg_psi_t (1:jy,1:kz,1:kz))
407 allocate (be%reg_psi_ps (1:jy,1:kz))
408 if ( cv_options == 6 ) then
409 allocate (be%reg_psi_rh (1:jy,1:kz,1:kz))
410 allocate (be%reg_chi_u_t (1:jy,1:kz,1:kz))
411 allocate (be%reg_chi_u_ps (1:jy,1:kz))
412 allocate (be%reg_chi_u_rh (1:jy,1:kz,1:kz))
413 allocate (be%reg_t_u_rh (1:jy,1:kz,1:kz))
414 allocate (be%reg_ps_u_rh (1:jy,1:kz))
415 end if
417 !-----------for interpolating CV5--------------------------------------------------------------
418 if ( interpolate_stats ) then
419 allocate (reg_psi_chi0(1:nj,1:nk))
420 allocate (reg_psi_ps0 (1:nj,1:nk))
421 allocate (reg_psi_t0 (1:nj,1:nk,1:nk))
422 end if
424 be%reg_psi_chi = 0.
425 be%reg_psi_t = 0.
426 be%reg_psi_ps = 0.
427 if ( cv_options == 6 ) then
428 be%reg_psi_rh = 0.
429 be%reg_chi_u_t = 0.
430 be%reg_chi_u_ps = 0.
431 be%reg_chi_u_rh = 0.
432 be%reg_t_u_rh = 0.
433 end if
435 if ( interpolate_stats ) then
436 do k=1,nk
437 do j=1,nj
438 reg_psi_chi0(j,k) = regcoeff_psi_chi(bin(1,j,k))
439 b = bin2d(1,j)
440 reg_psi_ps0(j,k) = regcoeff_psi_ps(k,b)
441 end do
442 end do
443 do k=1,nk
444 do i=1,nk
445 do j=1,nj
446 b = bin2d(1,j)
447 reg_psi_t0(j,i,k) = regcoeff_psi_t(i,k,b)
448 end do
449 end do
450 end do
451 end if ! if interpolate_stats
453 if ( .not. interpolate_stats ) then
454 do k=1,nk
455 do j=1,nj
456 be%reg_psi_chi(j,k) = psi_chi_factor * regcoeff_psi_chi(bin(1,j,k))
457 b = bin2d(1,j)
458 be%reg_psi_ps(j,k) = psi_ps_factor * regcoeff_psi_ps(k,b)
459 end do
460 end do
461 do k=1,nk
462 do i=1,nk
463 do j=1,nj
464 b = bin2d(1,j)
465 be%reg_psi_t(j,i,k) = psi_t_factor * regcoeff_psi_t(i,k,b)
466 end do
467 end do
468 end do
470 if ( cv_options == 6 ) then
471 do k=1,nk
472 do j=1,nj
473 b = bin2d(1,j)
474 be%reg_ps_u_rh(j,k) = ps_u_rh_factor * regcoeff_ps_u_rh(k,b)
475 be%reg_chi_u_ps(j,k) = chi_u_ps_factor * regcoeff_chi_u_ps(k,b)
476 end do
477 end do
478 do k=1,nk
479 do i=1,nk
480 do j=1,nj
481 b = bin2d(1,j)
482 be%reg_psi_rh(j,i,k) = psi_rh_factor * regcoeff_psi_rh(i,k,b)
483 be%reg_chi_u_t(j,i,k) = chi_u_t_factor * regcoeff_chi_u_t(i,k,b)
484 be%reg_chi_u_rh(j,i,k) = chi_u_rh_factor * regcoeff_chi_u_rh(i,k,b)
485 be%reg_t_u_rh(j,i,k) = t_u_rh_factor * regcoeff_t_u_rh(i,k,b)
486 end do
487 end do
488 end do
489 end if ! if cv_options 6
490 end if ! if not interpolate_stats
492 deallocate (regcoeff_psi_chi)
493 deallocate (regcoeff_psi_t)
494 deallocate (regcoeff_psi_ps)
495 if ( cv_options == 6 ) then
496 deallocate (regcoeff_psi_rh)
497 deallocate (regcoeff_chi_u_t)
498 deallocate (regcoeff_chi_u_ps)
499 deallocate (regcoeff_chi_u_rh)
500 deallocate (regcoeff_t_u_rh)
501 deallocate (regcoeff_ps_u_rh)
502 end if
503 ! 1.3 Domain_averaged regression coefficients
505 if (.not.lat_stats_option) then
506 write (unit=message(1), fmt='(a)') ' '
507 write (unit=message(2), fmt='(3x, a)') &
508 'Using the averaged regression coefficients for unbalanced part'
510 !-----------for interpolating CV5--------------------------------------------------------------
511 if ( interpolate_stats ) then
512 allocate (reg_psi_ps_be (1:nk))
513 allocate (reg_psi_chi_be(1:nk))
514 allocate (reg_psi_t_be (1:nk,1:nk))
515 do k=1,nk
516 reg_psi_ps_be(k) = sum(reg_psi_ps0 (:,k))/dble(nj)
517 reg_psi_chi_be(k) = sum(reg_psi_chi0(:,k))/dble(nj)
518 end do
519 do m=1,nk
520 do k=1,nk
521 reg_psi_t_be(k,m)=sum(reg_psi_t0(:,k,m))/dble(nj)
522 end do
523 end do
524 deallocate (reg_psi_chi0)
525 deallocate (reg_psi_ps0 )
526 deallocate (reg_psi_t0 )
527 end if ! if interpolate_stats
528 !===========for interpolating CV5==============================================================
530 if ( (.not. interpolate_stats) .and. (bin_type /= 5) ) then
531 do k=1,kz
532 be%reg_psi_ps (:,k) = sum(be%reg_psi_ps (:,k))/dble(jy)
533 be%reg_psi_chi (:,k) = sum(be%reg_psi_chi(:,k))/dble(jy)
534 end do
535 do m=1,kz
536 do k=1,kz
537 be%reg_psi_t (:,k,m) = sum(be%reg_psi_t(:,k,m))/dble(jy)
538 end do
539 end do
540 if ( cv_options == 6 ) then
541 do k=1,kz
542 be%reg_ps_u_rh (:,k) = sum(be%reg_ps_u_rh (:,k))/dble(jy)
543 be%reg_chi_u_ps(:,k) = sum(be%reg_chi_u_ps(:,k))/dble(jy)
544 end do
545 do m=1,kz
546 do k=1,kz
547 be%reg_psi_rh (:,k,m) = sum(be%reg_psi_rh (:,k,m))/dble(jy)
548 be%reg_chi_u_t (:,k,m) = sum(be%reg_chi_u_t (:,k,m))/dble(jy)
549 be%reg_chi_u_rh(:,k,m) = sum(be%reg_chi_u_rh(:,k,m))/dble(jy)
550 be%reg_t_u_rh (:,k,m) = sum(be%reg_t_u_rh (:,k,m))/dble(jy)
551 end do
552 end do
553 end if ! if cv_options 6
554 end if ! if not interpolate_stats
556 else
557 write (unit=message(1), fmt='(a)') ' '
558 write (unit=message(2), fmt='(3x, a)') &
559 'Using the geographically-dependent regression coefficients for unbalanced part'
560 end if
562 call da_message(message(1:2))
564 end if ! cv_options /= 7
566 ! 2.0 Load the eigenvector and eigenvalue
568 allocate (be1_eval_loc (1:nj,1:nk))
569 allocate (be2_eval_loc (1:nj,1:nk))
570 allocate (be3_eval_loc (1:nj,1:nk))
571 allocate (be4_eval_loc (1:nj,1:nk))
572 allocate (be5_eval_loc (1:nj,1:1))
574 if ( cloud_cv_options >= 2 ) then
575 allocate (be6_eval_loc (1:nj,1:nk))
576 allocate (be7_eval_loc (1:nj,1:nk))
577 allocate (be8_eval_loc (1:nj,1:nk))
578 allocate (be9_eval_loc (1:nj,1:nk))
579 allocate (be10_eval_loc (1:nj,1:nk))
580 end if
581 #if (WRF_CHEM == 1)
582 if ( chem_cv_options >= 10 ) then
583 allocate (be12_eval_loc (num_cv_3d_chem,1:nj,1:nk))
584 end if
585 #endif
586 if ( use_cv_w ) then
587 allocate (be11_eval_loc (1:nj,1:nk))
588 end if
590 if (vert_corr == vert_corr_2) then
592 allocate (be1_eval_glo(1:nk))
593 allocate (be2_eval_glo(1:nk))
594 allocate (be3_eval_glo(1:nk))
595 allocate (be4_eval_glo(1:nk))
596 allocate (be5_eval_glo(1:1))
598 if ( cloud_cv_options >= 2 ) then
599 allocate (be6_eval_glo(1:nk))
600 allocate (be7_eval_glo(1:nk))
601 allocate (be8_eval_glo(1:nk))
602 allocate (be9_eval_glo(1:nk))
603 allocate (be10_eval_glo(1:nk))
604 end if
606 #if (WRF_CHEM == 1)
607 if ( chem_cv_options >= 10 ) then
608 allocate (be12_eval_glo(num_cv_3d_chem,1:nk))
609 end if
610 #endif
612 if ( use_cv_w ) then
613 allocate (be11_eval_glo(1:nk))
614 end if
616 allocate (be1_evec_loc(1:nj,1:nk,1:nk))
617 allocate (be2_evec_loc(1:nj,1:nk,1:nk))
618 allocate (be3_evec_loc(1:nj,1:nk,1:nk))
619 allocate (be4_evec_loc(1:nj,1:nk,1:nk))
620 allocate (be5_evec_loc(1:nj,1: 1,1: 1))
622 if ( cloud_cv_options >= 2 ) then
623 allocate (be6_evec_loc (1:nj,1:nk,1:nk))
624 allocate (be7_evec_loc (1:nj,1:nk,1:nk))
625 allocate (be8_evec_loc (1:nj,1:nk,1:nk))
626 allocate (be9_evec_loc (1:nj,1:nk,1:nk))
627 allocate (be10_evec_loc(1:nj,1:nk,1:nk))
628 end if
629 #if (WRF_CHEM == 1)
630 if ( chem_cv_options >= 10 ) then
631 allocate (be12_evec_loc (num_cv_3d_chem,1:nj,1:nk,1:nk))
632 end if
633 #endif
634 if ( use_cv_w ) then
635 allocate (be11_evec_loc(1:nj,1:nk,1:nk))
636 end if
638 allocate (be1_evec_glo(1:nk,1:nk))
639 allocate (be2_evec_glo(1:nk,1:nk))
640 allocate (be3_evec_glo(1:nk,1:nk))
641 allocate (be4_evec_glo(1:nk,1:nk))
642 allocate (be5_evec_glo(1:1,1:1))
644 if ( cloud_cv_options >= 2 ) then
645 allocate (be6_evec_glo (1:nk,1:nk))
646 allocate (be7_evec_glo (1:nk,1:nk))
647 allocate (be8_evec_glo (1:nk,1:nk))
648 allocate (be9_evec_glo (1:nk,1:nk))
649 allocate (be10_evec_glo(1:nk,1:nk))
650 end if
652 #if (WRF_CHEM == 1)
653 if ( chem_cv_options >= 10 ) then
654 allocate (be12_evec_glo(num_cv_3d_chem,1:nk,1:nk))
655 end if
656 #endif
658 if ( use_cv_w ) then
659 allocate (be11_evec_glo(1:nk,1:nk))
660 end if
662 end if
664 if ( use_rf ) then
665 if ( cloud_cv_options >= 2 ) then
666 allocate ( be6_rf_lengthscale(1:kz))
667 allocate ( be7_rf_lengthscale(1:kz))
668 allocate ( be8_rf_lengthscale(1:kz))
669 allocate ( be9_rf_lengthscale(1:kz))
670 allocate (be10_rf_lengthscale(1:kz))
671 end if
672 #if (WRF_CHEM == 1)
673 if ( chem_cv_options >= 10 ) then
674 allocate (be12_rf_lengthscale(num_cv_3d_chem,1:kz))
675 end if
676 #endif
677 if ( use_cv_w ) then
678 allocate (be11_rf_lengthscale(1:kz))
679 end if
680 end if
682 if ( cloud_cv_options == 3 ) then
683 ! hard-coded the v6-v11 BE values here
684 be % v6 % name = "qcloud"
685 be % v7 % name = "qrain"
686 be % v8 % name = "qice"
687 be % v9 % name = "qsnow"
688 be % v10 % name = "qgraup"
689 be6_eval_glo = 1.0e-6
690 be7_eval_glo = 1.0e-6
691 be8_eval_glo = 1.0e-6
692 be9_eval_glo = 1.0e-6
693 be10_eval_glo = 1.0e-6
694 if ( use_cv_w ) then
695 be % v11 % name = "w"
696 be11_eval_glo = 1.0
697 end if
698 if ( use_rf ) then
699 be6_rf_lengthscale = 1.0
700 be7_rf_lengthscale = 1.0
701 be8_rf_lengthscale = 1.0
702 be9_rf_lengthscale = 1.0
703 be10_rf_lengthscale = 1.0
704 if ( use_cv_w ) then
705 be11_rf_lengthscale = 1.0
706 end if
707 end if
708 end if
710 ! 2.2 Read in the eigenvector and eigenvalue
712 print *, '-------- reading eigen vector/value -------'
713 do i = 1 , num_cv_3d_basic
714 read (be_unit) variable
715 read (be_unit) nk, num_bins2d
716 print *, trim(adjustl(variable)), nk, num_bins2d
717 if ( i == 1 ) then
718 allocate (evec_loc(1:nk,1:nk,1:num_bins2d))
719 allocate (eval_loc(1:nk, 1:num_bins2d))
720 end if
721 select case( trim(adjustl(variable)) )
722 case ('psi', 'u')
723 be % v1 % name = trim(adjustl(variable))
724 read (be_unit) be1_evec_glo
725 read (be_unit) be1_eval_glo
726 read (be_unit) evec_loc
727 read (be_unit) eval_loc
728 do j=1,nj
729 b = bin2d(1,j)
730 be1_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
731 be1_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
732 end do
734 case ('chi_u', 'v')
735 be % v2 % name = trim(adjustl(variable))
736 read (be_unit) be2_evec_glo
737 read (be_unit) be2_eval_glo
738 read (be_unit) evec_loc
739 read (be_unit) eval_loc
740 do j=1,nj
741 b = bin2d(1,j)
742 be2_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
743 be2_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
744 end do
746 case ('t_u', 't')
747 be % v3 % name = trim(adjustl(variable))
748 read (be_unit) be3_evec_glo
749 read (be_unit) be3_eval_glo
750 read (be_unit) evec_loc
751 read (be_unit) eval_loc
752 do j=1,nj
753 b = bin2d(1,j)
754 be3_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
755 be3_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
756 end do
758 case ('rh_u' , 'rh' )
759 be % v4 % name = trim(adjustl(variable))
760 read (be_unit) be4_evec_glo
761 read (be_unit) be4_eval_glo
762 read (be_unit) evec_loc
763 read (be_unit) eval_loc
764 do j=1,nj
765 b = bin2d(1,j)
766 be4_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
767 be4_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
768 end do
770 case default;
771 message(1)=' Read problem in eigen vectors/values in BE file '
772 write (unit=message(2),fmt='(A,A)') ' Trying to read Eigenvectors for variable: ',trim(adjustl(variable))
773 write (unit=message(3),fmt='(A)') ' Make sure you are using the correct be.dat file for your cv_options setting!'
774 call da_error(__FILE__,__LINE__,message(1:3))
775 end select
777 end do !1-num_cv_3d_basic
779 if ( cloud_cv_options == 2 ) then
780 print*,"cloud_cv_options == 2 variables from id ",num_cv_3d_basic+1 ,"to id", num_cv_3d_basic+num_cv_3d_extra
781 do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra
783 read (be_unit) variable
784 read (be_unit) nk, num_bins2d
785 print *, trim(adjustl(variable)), nk, num_bins2d
787 select case( trim(adjustl(variable)) )
789 case ('qcloud' )
790 be % v6 % name = trim(adjustl(variable))
791 read (be_unit) be6_evec_glo
792 read (be_unit) be6_eval_glo
793 read (be_unit) evec_loc
794 read (be_unit) eval_loc
795 do j=1,nj
796 b = bin2d(1,j)
797 be6_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
798 be6_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
799 end do
801 case ('qrain' )
802 be % v7 % name = trim(adjustl(variable))
803 read (be_unit) be7_evec_glo
804 read (be_unit) be7_eval_glo
805 read (be_unit) evec_loc
806 read (be_unit) eval_loc
807 do j=1,nj
808 b = bin2d(1,j)
809 be7_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
810 be7_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
811 end do
813 case ('qice' )
814 be % v8 % name = trim(adjustl(variable))
815 read (be_unit) be8_evec_glo
816 read (be_unit) be8_eval_glo
817 read (be_unit) evec_loc
818 read (be_unit) eval_loc
819 do j=1,nj
820 b = bin2d(1,j)
821 be8_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
822 be8_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
823 end do
825 case ('qsnow' )
826 be % v9 % name = trim(adjustl(variable))
827 read (be_unit) be9_evec_glo
828 read (be_unit) be9_eval_glo
829 read (be_unit) evec_loc
830 read (be_unit) eval_loc
831 do j=1,nj
832 b = bin2d(1,j)
833 be9_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
834 be9_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
835 end do
837 case ('qgraup' )
838 be % v10 % name = trim(adjustl(variable))
839 !read (be_unit) nk, num_bins2d
840 read (be_unit) be10_evec_glo
841 read (be_unit) be10_eval_glo
842 read (be_unit) evec_loc
843 read (be_unit) eval_loc
844 do j=1,nj
845 b = bin2d(1,j)
846 be10_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
847 be10_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
848 end do
850 case ('w' )
851 !if ( use_cv_w ) then
852 ! if ( trim(adjustl(variable)) == 'w' ) then
853 be % v11 % name = trim(adjustl(variable))
854 read (be_unit) be11_evec_glo
855 read (be_unit) be11_eval_glo
856 read (be_unit) evec_loc
857 read (be_unit) eval_loc
858 do j=1,nj
859 b = bin2d(1,j)
860 be11_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
861 be11_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
862 end do
863 ! end if
864 !end if
866 case default;
867 message(1)=' Read problem in eigen vectors/values in BE file '
868 write (unit=message(2),fmt='(A,A)') ' Trying to read Eigenvectors for variable: ',trim(adjustl(variable))
869 write (unit=message(3),fmt='(A)') ' Make sure you are using the correct be.dat file for your cv_options setting!'
870 call da_error(__FILE__,__LINE__,message(1:3))
871 end select
873 end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra-1
874 print*,"reading variables for cloud_cv_options=2 is done"
875 end if ! cloud_cv_options=2
877 #if (WRF_CHEM == 1)
878 if ( chem_cv_options == 10 ) then
879 print*,"chem_cv_options == 10 variables from id ",num_cv_3d_basic+1 ,"to id", num_cv_3d_basic+num_cv_3d_extra
880 do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra
882 read (be_unit) variable
883 read (be_unit) nk, num_bins2d
884 print *, trim(adjustl(variable)), nk, num_bins2d
886 select case( trim(adjustl(variable)) )
888 case ('dust_1' )
889 ispcc=p_dust_1-1
890 be % v12(ispcc) % name = trim(adjustl(variable))
891 read (be_unit) be12_evec_glo(ispcc,:,:)
892 read (be_unit) be12_eval_glo(ispcc,:)
893 read (be_unit) evec_loc
894 read (be_unit) eval_loc
895 do j=1,nj
896 b = bin2d(1,j)
897 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
898 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
899 end do
901 case ('dust_2' )
902 ispcc=p_dust_2-1
903 be % v12(ispcc) % name = trim(adjustl(variable))
904 read (be_unit) be12_evec_glo(ispcc,:,:)
905 read (be_unit) be12_eval_glo(ispcc,:)
906 read (be_unit) evec_loc
907 read (be_unit) eval_loc
908 do j=1,nj
909 b = bin2d(1,j)
910 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
911 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
912 end do
914 case ('dust_3' )
915 ispcc=p_dust_3-1
916 be % v12(ispcc) % name = trim(adjustl(variable))
917 read (be_unit) be12_evec_glo(ispcc,:,:)
918 read (be_unit) be12_eval_glo(ispcc,:)
919 read (be_unit) evec_loc
920 read (be_unit) eval_loc
921 do j=1,nj
922 b = bin2d(1,j)
923 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
924 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
925 end do
927 case ('dust_4' )
928 ispcc=p_dust_4-1
929 be % v12(ispcc) % name = trim(adjustl(variable))
930 read (be_unit) be12_evec_glo(ispcc,:,:)
931 read (be_unit) be12_eval_glo(ispcc,:)
932 read (be_unit) evec_loc
933 read (be_unit) eval_loc
934 do j=1,nj
935 b = bin2d(1,j)
936 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
937 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
938 end do
940 case ('seas_1' )
941 ispcc=p_seas_1-1
942 be % v12(ispcc) % name = trim(adjustl(variable))
943 read (be_unit) be12_evec_glo(ispcc,:,:)
944 read (be_unit) be12_eval_glo(ispcc,:)
945 read (be_unit) evec_loc
946 read (be_unit) eval_loc
947 do j=1,nj
948 b = bin2d(1,j)
949 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
950 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
951 end do
953 case ('seas_2' )
954 ispcc=p_seas_2-1
955 be % v12(ispcc) % name = trim(adjustl(variable))
956 read (be_unit) be12_evec_glo(ispcc,:,:)
957 read (be_unit) be12_eval_glo(ispcc,:)
958 read (be_unit) evec_loc
959 read (be_unit) eval_loc
960 do j=1,nj
961 b = bin2d(1,j)
962 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
963 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
964 end do
966 case ('seas_3' )
967 ispcc=p_seas_3-1
968 be % v12(ispcc) % name = trim(adjustl(variable))
969 read (be_unit) be12_evec_glo(ispcc,:,:)
970 read (be_unit) be12_eval_glo(ispcc,:)
971 read (be_unit) evec_loc
972 read (be_unit) eval_loc
973 do j=1,nj
974 b = bin2d(1,j)
975 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
976 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
977 end do
979 case ('seas_4' )
980 ispcc=p_seas_4-1
981 be % v12(ispcc) % name = trim(adjustl(variable))
982 read (be_unit) be12_evec_glo(ispcc,:,:)
983 read (be_unit) be12_eval_glo(ispcc,:)
984 read (be_unit) evec_loc
985 read (be_unit) eval_loc
986 do j=1,nj
987 b = bin2d(1,j)
988 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
989 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
990 end do
992 case ('oc1' )
993 ispcc=p_oc1-1
994 be % v12(ispcc) % name = trim(adjustl(variable))
995 read (be_unit) be12_evec_glo(ispcc,:,:)
996 read (be_unit) be12_eval_glo(ispcc,:)
997 read (be_unit) evec_loc
998 read (be_unit) eval_loc
999 do j=1,nj
1000 b = bin2d(1,j)
1001 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1002 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1003 end do
1005 case ('oc2' )
1006 ispcc=p_oc2-1
1007 be % v12(ispcc) % name = trim(adjustl(variable))
1008 read (be_unit) be12_evec_glo(ispcc,:,:)
1009 read (be_unit) be12_eval_glo(ispcc,:)
1010 read (be_unit) evec_loc
1011 read (be_unit) eval_loc
1012 do j=1,nj
1013 b = bin2d(1,j)
1014 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1015 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1016 end do
1018 case ('bc1' )
1019 ispcc=p_bc1-1
1020 be % v12(ispcc) % name = trim(adjustl(variable))
1021 read (be_unit) be12_evec_glo(ispcc,:,:)
1022 read (be_unit) be12_eval_glo(ispcc,:)
1023 read (be_unit) evec_loc
1024 read (be_unit) eval_loc
1025 do j=1,nj
1026 b = bin2d(1,j)
1027 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1028 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1029 end do
1031 case ('bc2' )
1032 ispcc=p_bc2-1
1033 be % v12(ispcc) % name = trim(adjustl(variable))
1034 read (be_unit) be12_evec_glo(ispcc,:,:)
1035 read (be_unit) be12_eval_glo(ispcc,:)
1036 read (be_unit) evec_loc
1037 read (be_unit) eval_loc
1038 do j=1,nj
1039 b = bin2d(1,j)
1040 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1041 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1042 end do
1044 case ('so2' )
1045 ispcc=p_so2-1
1046 be % v12(ispcc) % name = trim(adjustl(variable))
1047 read (be_unit) be12_evec_glo(ispcc,:,:)
1048 read (be_unit) be12_eval_glo(ispcc,:)
1049 read (be_unit) evec_loc
1050 read (be_unit) eval_loc
1051 do j=1,nj
1052 b = bin2d(1,j)
1053 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1054 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1055 end do
1057 case ('no2' )
1058 ispcc=p_no2-1
1059 be % v12(ispcc) % name = trim(adjustl(variable))
1060 read (be_unit) be12_evec_glo(ispcc,:,:)
1061 read (be_unit) be12_eval_glo(ispcc,:)
1062 read (be_unit) evec_loc
1063 read (be_unit) eval_loc
1064 do j=1,nj
1065 b = bin2d(1,j)
1066 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1067 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1068 end do
1070 case ('o3' )
1071 ispcc=p_o3-1
1072 be % v12(ispcc) % name = trim(adjustl(variable))
1073 read (be_unit) be12_evec_glo(ispcc,:,:)
1074 read (be_unit) be12_eval_glo(ispcc,:)
1075 read (be_unit) evec_loc
1076 read (be_unit) eval_loc
1077 do j=1,nj
1078 b = bin2d(1,j)
1079 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1080 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1081 end do
1083 case ('co' )
1084 ispcc=p_co-1
1085 be % v12(ispcc) % name = trim(adjustl(variable))
1086 read (be_unit) be12_evec_glo(ispcc,:,:)
1087 read (be_unit) be12_eval_glo(ispcc,:)
1088 read (be_unit) evec_loc
1089 read (be_unit) eval_loc
1090 do j=1,nj
1091 b = bin2d(1,j)
1092 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1093 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1094 end do
1096 case ('sulf' )
1097 ispcc=p_sulf-1
1098 be % v12(ispcc) % name = trim(adjustl(variable))
1099 read (be_unit) be12_evec_glo(ispcc,:,:)
1100 read (be_unit) be12_eval_glo(ispcc,:)
1101 read (be_unit) evec_loc
1102 read (be_unit) eval_loc
1103 do j=1,nj
1104 b = bin2d(1,j)
1105 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1106 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1107 end do
1109 case ('p25' )
1110 ispcc=p_p25-1
1111 be % v12(ispcc) % name = trim(adjustl(variable))
1112 read (be_unit) be12_evec_glo(ispcc,:,:)
1113 read (be_unit) be12_eval_glo(ispcc,:)
1114 read (be_unit) evec_loc
1115 read (be_unit) eval_loc
1116 do j=1,nj
1117 b = bin2d(1,j)
1118 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1119 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1120 end do
1122 case ('p10' )
1123 ispcc=p_p10-1
1124 be % v12(ispcc) % name = trim(adjustl(variable))
1125 read (be_unit) be12_evec_glo(ispcc,:,:)
1126 read (be_unit) be12_eval_glo(ispcc,:)
1127 read (be_unit) evec_loc
1128 read (be_unit) eval_loc
1129 do j=1,nj
1130 b = bin2d(1,j)
1131 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1132 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1133 end do
1135 case ('w' )
1136 !if ( use_cv_w ) then
1137 ! if ( trim(adjustl(variable)) == 'w' ) then
1138 be % v11 % name = trim(adjustl(variable))
1139 read (be_unit) be11_evec_glo
1140 read (be_unit) be11_eval_glo
1141 read (be_unit) evec_loc
1142 read (be_unit) eval_loc
1143 do j=1,nj
1144 b = bin2d(1,j)
1145 be11_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1146 be11_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
1147 end do
1148 ! end if
1149 !end if
1151 case default;
1152 message(1)=' Read problem in eigen vectors/values in BE file '
1153 write (unit=message(2),fmt='(A,A)') ' Trying to read Eigenvectors for variable: ',trim(adjustl(variable))
1154 write (unit=message(3),fmt='(A)') ' Make sure you are using the correct be.dat file for your cv_options setting!'
1155 call da_error(__FILE__,__LINE__,message(1:3))
1156 end select
1158 end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra-1
1159 print*,"reading variables for chem_cv_options=10 is done"
1160 end if ! chem_cv_options=10
1163 if ( chem_cv_options == 20 ) then
1164 print*,"chem_cv_options == 20 variables from id ",num_cv_3d_basic+1 ,"to id", num_cv_3d_basic+num_cv_3d_extra
1165 do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra
1166 read (be_unit) variable
1167 read (be_unit) nk, num_bins2d
1168 print *, trim(adjustl(variable)), nk, num_bins2d
1169 select case( trim(adjustl(variable)) )
1171 case ('bc_1' )
1172 ispcc=p_bc_a01-1
1173 be % v12(ispcc) % name = trim(adjustl(variable))
1174 read (be_unit) be12_evec_glo(ispcc,:,:)
1175 read (be_unit) be12_eval_glo(ispcc,:)
1176 read (be_unit) evec_loc
1177 read (be_unit) eval_loc
1178 do j=1,nj
1179 b = bin2d(1,j)
1180 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1181 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1182 end do
1184 case ('bc_2' )
1185 ispcc=p_bc_a02-1
1186 be % v12(ispcc) % name = trim(adjustl(variable))
1187 read (be_unit) be12_evec_glo(ispcc,:,:)
1188 read (be_unit) be12_eval_glo(ispcc,:)
1189 read (be_unit) evec_loc
1190 read (be_unit) eval_loc
1191 do j=1,nj
1192 b = bin2d(1,j)
1193 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1194 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1195 end do
1197 case ('bc_3' )
1198 ispcc=p_bc_a03-1
1199 be % v12(ispcc) % name = trim(adjustl(variable))
1200 read (be_unit) be12_evec_glo(ispcc,:,:)
1201 read (be_unit) be12_eval_glo(ispcc,:)
1202 read (be_unit) evec_loc
1203 read (be_unit) eval_loc
1204 do j=1,nj
1205 b = bin2d(1,j)
1206 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1207 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1208 end do
1210 case ('bc_4' )
1211 ispcc=p_bc_a04-1
1212 be % v12(ispcc) % name = trim(adjustl(variable))
1213 read (be_unit) be12_evec_glo(ispcc,:,:)
1214 read (be_unit) be12_eval_glo(ispcc,:)
1215 read (be_unit) evec_loc
1216 read (be_unit) eval_loc
1217 do j=1,nj
1218 b = bin2d(1,j)
1219 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1220 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1221 end do
1223 case ('oc_1' )
1224 ispcc=p_oc_a01-1
1225 be % v12(ispcc) % name = trim(adjustl(variable))
1226 read (be_unit) be12_evec_glo(ispcc,:,:)
1227 read (be_unit) be12_eval_glo(ispcc,:)
1228 read (be_unit) evec_loc
1229 read (be_unit) eval_loc
1230 do j=1,nj
1231 b = bin2d(1,j)
1232 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1233 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1234 end do
1236 case ('oc_2' )
1237 ispcc=p_oc_a02-1
1238 be % v12(ispcc) % name = trim(adjustl(variable))
1239 read (be_unit) be12_evec_glo(ispcc,:,:)
1240 read (be_unit) be12_eval_glo(ispcc,:)
1241 read (be_unit) evec_loc
1242 read (be_unit) eval_loc
1243 do j=1,nj
1244 b = bin2d(1,j)
1245 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1246 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1247 end do
1249 case ('oc_3' )
1250 ispcc=p_oc_a03-1
1251 be % v12(ispcc) % name = trim(adjustl(variable))
1252 read (be_unit) be12_evec_glo(ispcc,:,:)
1253 read (be_unit) be12_eval_glo(ispcc,:)
1254 read (be_unit) evec_loc
1255 read (be_unit) eval_loc
1256 do j=1,nj
1257 b = bin2d(1,j)
1258 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1259 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1260 end do
1262 case ('oc_4' )
1263 ispcc=p_oc_a04-1
1264 be % v12(ispcc) % name = trim(adjustl(variable))
1265 read (be_unit) be12_evec_glo(ispcc,:,:)
1266 read (be_unit) be12_eval_glo(ispcc,:)
1267 read (be_unit) evec_loc
1268 read (be_unit) eval_loc
1269 do j=1,nj
1270 b = bin2d(1,j)
1271 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1272 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1273 end do
1275 case ('so4_1' )
1276 ispcc=p_so4_a01-1
1277 be % v12(ispcc) % name = trim(adjustl(variable))
1278 read (be_unit) be12_evec_glo(ispcc,:,:)
1279 read (be_unit) be12_eval_glo(ispcc,:)
1280 read (be_unit) evec_loc
1281 read (be_unit) eval_loc
1282 do j=1,nj
1283 b = bin2d(1,j)
1284 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1285 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1286 end do
1288 case ('so4_2' )
1289 ispcc=p_so4_a02-1
1290 be % v12(ispcc) % name = trim(adjustl(variable))
1291 read (be_unit) be12_evec_glo(ispcc,:,:)
1292 read (be_unit) be12_eval_glo(ispcc,:)
1293 read (be_unit) evec_loc
1294 read (be_unit) eval_loc
1295 do j=1,nj
1296 b = bin2d(1,j)
1297 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1298 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1299 end do
1301 case ('so4_3' )
1302 ispcc=p_so4_a03-1
1303 be % v12(ispcc) % name = trim(adjustl(variable))
1304 read (be_unit) be12_evec_glo(ispcc,:,:)
1305 read (be_unit) be12_eval_glo(ispcc,:)
1306 read (be_unit) evec_loc
1307 read (be_unit) eval_loc
1308 do j=1,nj
1309 b = bin2d(1,j)
1310 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1311 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1312 end do
1314 case ('so4_4' )
1315 ispcc=p_so4_a04-1
1316 be % v12(ispcc) % name = trim(adjustl(variable))
1317 read (be_unit) be12_evec_glo(ispcc,:,:)
1318 read (be_unit) be12_eval_glo(ispcc,:)
1319 read (be_unit) evec_loc
1320 read (be_unit) eval_loc
1321 do j=1,nj
1322 b = bin2d(1,j)
1323 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1324 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1325 end do
1327 case ('no3_1' )
1328 ispcc=p_no3_a01-1
1329 be % v12(ispcc) % name = trim(adjustl(variable))
1330 read (be_unit) be12_evec_glo(ispcc,:,:)
1331 read (be_unit) be12_eval_glo(ispcc,:)
1332 read (be_unit) evec_loc
1333 read (be_unit) eval_loc
1334 do j=1,nj
1335 b = bin2d(1,j)
1336 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1337 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1338 end do
1340 case ('no3_2' )
1341 ispcc=p_no3_a02-1
1342 be % v12(ispcc) % name = trim(adjustl(variable))
1343 read (be_unit) be12_evec_glo(ispcc,:,:)
1344 read (be_unit) be12_eval_glo(ispcc,:)
1345 read (be_unit) evec_loc
1346 read (be_unit) eval_loc
1347 do j=1,nj
1348 b = bin2d(1,j)
1349 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1350 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1351 end do
1353 case ('no3_3' )
1354 ispcc=p_no3_a03-1
1355 be % v12(ispcc) % name = trim(adjustl(variable))
1356 read (be_unit) be12_evec_glo(ispcc,:,:)
1357 read (be_unit) be12_eval_glo(ispcc,:)
1358 read (be_unit) evec_loc
1359 read (be_unit) eval_loc
1360 do j=1,nj
1361 b = bin2d(1,j)
1362 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1363 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1364 end do
1366 case ('no3_4' )
1367 ispcc=p_no3_a04-1
1368 be % v12(ispcc) % name = trim(adjustl(variable))
1369 read (be_unit) be12_evec_glo(ispcc,:,:)
1370 read (be_unit) be12_eval_glo(ispcc,:)
1371 read (be_unit) evec_loc
1372 read (be_unit) eval_loc
1373 do j=1,nj
1374 b = bin2d(1,j)
1375 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1376 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1377 end do
1379 case ('nh4_1' )
1380 ispcc=p_nh4_a01-1
1381 be % v12(ispcc) % name = trim(adjustl(variable))
1382 read (be_unit) be12_evec_glo(ispcc,:,:)
1383 read (be_unit) be12_eval_glo(ispcc,:)
1384 read (be_unit) evec_loc
1385 read (be_unit) eval_loc
1386 do j=1,nj
1387 b = bin2d(1,j)
1388 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1389 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1390 end do
1392 case ('nh4_2' )
1393 ispcc=p_nh4_a02-1
1394 be % v12(ispcc) % name = trim(adjustl(variable))
1395 read (be_unit) be12_evec_glo(ispcc,:,:)
1396 read (be_unit) be12_eval_glo(ispcc,:)
1397 read (be_unit) evec_loc
1398 read (be_unit) eval_loc
1399 do j=1,nj
1400 b = bin2d(1,j)
1401 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1402 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1403 end do
1405 case ('nh4_3' )
1406 ispcc=p_nh4_a03-1
1407 be % v12(ispcc) % name = trim(adjustl(variable))
1408 read (be_unit) be12_evec_glo(ispcc,:,:)
1409 read (be_unit) be12_eval_glo(ispcc,:)
1410 read (be_unit) evec_loc
1411 read (be_unit) eval_loc
1412 do j=1,nj
1413 b = bin2d(1,j)
1414 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1415 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1416 end do
1418 case ('nh4_4' )
1419 ispcc=p_nh4_a04-1
1420 be % v12(ispcc) % name = trim(adjustl(variable))
1421 read (be_unit) be12_evec_glo(ispcc,:,:)
1422 read (be_unit) be12_eval_glo(ispcc,:)
1423 read (be_unit) evec_loc
1424 read (be_unit) eval_loc
1425 do j=1,nj
1426 b = bin2d(1,j)
1427 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1428 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1429 end do
1431 case ('cl_1' )
1432 ispcc=p_cl_a01-1
1433 be % v12(ispcc) % name = trim(adjustl(variable))
1434 read (be_unit) be12_evec_glo(ispcc,:,:)
1435 read (be_unit) be12_eval_glo(ispcc,:)
1436 read (be_unit) evec_loc
1437 read (be_unit) eval_loc
1438 do j=1,nj
1439 b = bin2d(1,j)
1440 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1441 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1442 end do
1444 case ('cl_2' )
1445 ispcc=p_cl_a02-1
1446 be % v12(ispcc) % name = trim(adjustl(variable))
1447 read (be_unit) be12_evec_glo(ispcc,:,:)
1448 read (be_unit) be12_eval_glo(ispcc,:)
1449 read (be_unit) evec_loc
1450 read (be_unit) eval_loc
1451 do j=1,nj
1452 b = bin2d(1,j)
1453 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1454 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1455 end do
1457 case ('cl_3' )
1458 ispcc=p_cl_a03-1
1459 be % v12(ispcc) % name = trim(adjustl(variable))
1460 read (be_unit) be12_evec_glo(ispcc,:,:)
1461 read (be_unit) be12_eval_glo(ispcc,:)
1462 read (be_unit) evec_loc
1463 read (be_unit) eval_loc
1464 do j=1,nj
1465 b = bin2d(1,j)
1466 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1467 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1468 end do
1470 case ('cl_4' )
1471 ispcc=p_cl_a04-1
1472 be % v12(ispcc) % name = trim(adjustl(variable))
1473 read (be_unit) be12_evec_glo(ispcc,:,:)
1474 read (be_unit) be12_eval_glo(ispcc,:)
1475 read (be_unit) evec_loc
1476 read (be_unit) eval_loc
1477 do j=1,nj
1478 b = bin2d(1,j)
1479 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1480 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1481 end do
1483 case ('na_1' )
1484 ispcc=p_na_a01-1
1485 be % v12(ispcc) % name = trim(adjustl(variable))
1486 read (be_unit) be12_evec_glo(ispcc,:,:)
1487 read (be_unit) be12_eval_glo(ispcc,:)
1488 read (be_unit) evec_loc
1489 read (be_unit) eval_loc
1490 do j=1,nj
1491 b = bin2d(1,j)
1492 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1493 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1494 end do
1496 case ('na_2' )
1497 ispcc=p_na_a02-1
1498 be % v12(ispcc) % name = trim(adjustl(variable))
1499 read (be_unit) be12_evec_glo(ispcc,:,:)
1500 read (be_unit) be12_eval_glo(ispcc,:)
1501 read (be_unit) evec_loc
1502 read (be_unit) eval_loc
1503 do j=1,nj
1504 b = bin2d(1,j)
1505 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1506 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1507 end do
1509 case ('na_3' )
1510 ispcc=p_na_a03-1
1511 be % v12(ispcc) % name = trim(adjustl(variable))
1512 read (be_unit) be12_evec_glo(ispcc,:,:)
1513 read (be_unit) be12_eval_glo(ispcc,:)
1514 read (be_unit) evec_loc
1515 read (be_unit) eval_loc
1516 do j=1,nj
1517 b = bin2d(1,j)
1518 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1519 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1520 end do
1522 case ('na_4' )
1523 ispcc=p_na_a04-1
1524 be % v12(ispcc) % name = trim(adjustl(variable))
1525 read (be_unit) be12_evec_glo(ispcc,:,:)
1526 read (be_unit) be12_eval_glo(ispcc,:)
1527 read (be_unit) evec_loc
1528 read (be_unit) eval_loc
1529 do j=1,nj
1530 b = bin2d(1,j)
1531 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1532 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1533 end do
1535 case ('oin_1' )
1536 ispcc=p_oin_a01-1
1537 be % v12(ispcc) % name = trim(adjustl(variable))
1538 read (be_unit) be12_evec_glo(ispcc,:,:)
1539 read (be_unit) be12_eval_glo(ispcc,:)
1540 read (be_unit) evec_loc
1541 read (be_unit) eval_loc
1542 do j=1,nj
1543 b = bin2d(1,j)
1544 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1545 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1546 end do
1548 case ('oin_2' )
1549 ispcc=p_oin_a02-1
1550 be % v12(ispcc) % name = trim(adjustl(variable))
1551 read (be_unit) be12_evec_glo(ispcc,:,:)
1552 read (be_unit) be12_eval_glo(ispcc,:)
1553 read (be_unit) evec_loc
1554 read (be_unit) eval_loc
1555 do j=1,nj
1556 b = bin2d(1,j)
1557 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1558 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1559 end do
1561 case ('oin_3' )
1562 ispcc=p_oin_a03-1
1563 be % v12(ispcc) % name = trim(adjustl(variable))
1564 read (be_unit) be12_evec_glo(ispcc,:,:)
1565 read (be_unit) be12_eval_glo(ispcc,:)
1566 read (be_unit) evec_loc
1567 read (be_unit) eval_loc
1568 do j=1,nj
1569 b = bin2d(1,j)
1570 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1571 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1572 end do
1574 case ('oin_4' )
1575 ispcc=p_oin_a04-1
1576 be % v12(ispcc) % name = trim(adjustl(variable))
1577 read (be_unit) be12_evec_glo(ispcc,:,:)
1578 read (be_unit) be12_eval_glo(ispcc,:)
1579 read (be_unit) evec_loc
1580 read (be_unit) eval_loc
1581 do j=1,nj
1582 b = bin2d(1,j)
1583 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1584 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1585 end do
1587 case ('so2' )
1588 ispcc=p_so2-1
1589 be % v12(ispcc) % name = trim(adjustl(variable))
1590 read (be_unit) be12_evec_glo(ispcc,:,:)
1591 read (be_unit) be12_eval_glo(ispcc,:)
1592 read (be_unit) evec_loc
1593 read (be_unit) eval_loc
1594 do j=1,nj
1595 b = bin2d(1,j)
1596 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1597 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1598 end do
1600 case ('no2' )
1601 ispcc=p_no2-1
1602 be % v12(ispcc) % name = trim(adjustl(variable))
1603 read (be_unit) be12_evec_glo(ispcc,:,:)
1604 read (be_unit) be12_eval_glo(ispcc,:)
1605 read (be_unit) evec_loc
1606 read (be_unit) eval_loc
1607 do j=1,nj
1608 b = bin2d(1,j)
1609 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1610 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1611 end do
1613 case ('o3' )
1614 ispcc=p_o3-1
1615 be % v12(ispcc) % name = trim(adjustl(variable))
1616 read (be_unit) be12_evec_glo(ispcc,:,:)
1617 read (be_unit) be12_eval_glo(ispcc,:)
1618 read (be_unit) evec_loc
1619 read (be_unit) eval_loc
1620 do j=1,nj
1621 b = bin2d(1,j)
1622 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1623 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1624 end do
1626 case ('co' )
1627 ispcc=p_co-1
1628 be % v12(ispcc) % name = trim(adjustl(variable))
1629 read (be_unit) be12_evec_glo(ispcc,:,:)
1630 read (be_unit) be12_eval_glo(ispcc,:)
1631 read (be_unit) evec_loc
1632 read (be_unit) eval_loc
1633 do j=1,nj
1634 b = bin2d(1,j)
1635 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1636 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1637 end do
1639 case ('w' )
1640 if ( use_cv_w ) then
1641 be % v11 % name = trim(adjustl(variable))
1642 read (be_unit) be11_evec_glo
1643 read (be_unit) be11_eval_glo
1644 read (be_unit) evec_loc
1645 read (be_unit) eval_loc
1646 do j=1,nj
1647 b = bin2d(1,j)
1648 be11_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1649 be11_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
1650 end do
1651 end if
1653 case default;
1654 message(1)=' Read problem in eigen vectors/values in BE file '
1655 write (unit=message(2),fmt='(A,A)') ' Trying to read Eigenvectors for variable: ',trim(adjustl(variable))
1656 write (unit=message(3),fmt='(A)') ' Make sure you are using the correct be.dat file for your cv_options setting!'
1657 call da_error(__FILE__,__LINE__,message(1:3))
1658 end select
1660 end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra-1
1661 print*,"reading variables for chem_cv_options=20 is done"
1662 end if ! chem_cv_options=20
1664 if ( chem_cv_options == 108 ) then ! racm_soa_vbs_da
1665 print*,"chem_cv_options == 108 variables from id ",num_cv_3d_basic+1 ,"to id", &
1666 num_cv_3d_basic+num_cv_3d_extra
1667 do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra
1668 read (be_unit) variable
1669 read (be_unit) nk, num_bins2d
1670 !print *, trim(adjustl(variable)), nk, num_bins2d
1672 select case( trim(adjustl(variable)) )
1674 case ('so4aj' )
1675 ispcc=p_so4aj-1
1676 be % v12(ispcc) % name = trim(adjustl(variable))
1677 read (be_unit) be12_evec_glo(ispcc,:,:)
1678 read (be_unit) be12_eval_glo(ispcc,:)
1679 read (be_unit) evec_loc
1680 read (be_unit) eval_loc
1681 do j=1,nj
1682 b = bin2d(1,j)
1683 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1684 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1685 end do
1687 case ('so4ai' )
1688 ispcc=p_so4ai-1
1689 be % v12(ispcc) % name = trim(adjustl(variable))
1690 read (be_unit) be12_evec_glo(ispcc,:,:)
1691 read (be_unit) be12_eval_glo(ispcc,:)
1692 read (be_unit) evec_loc
1693 read (be_unit) eval_loc
1694 do j=1,nj
1695 b = bin2d(1,j)
1696 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1697 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1698 end do
1700 case ('nh4aj' )
1701 ispcc=p_nh4aj-1
1702 be % v12(ispcc) % name = trim(adjustl(variable))
1703 read (be_unit) be12_evec_glo(ispcc,:,:)
1704 read (be_unit) be12_eval_glo(ispcc,:)
1705 read (be_unit) evec_loc
1706 read (be_unit) eval_loc
1707 do j=1,nj
1708 b = bin2d(1,j)
1709 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1710 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1711 end do
1713 case ('nh4ai' )
1714 ispcc=p_nh4ai-1
1715 be % v12(ispcc) % name = trim(adjustl(variable))
1716 read (be_unit) be12_evec_glo(ispcc,:,:)
1717 read (be_unit) be12_eval_glo(ispcc,:)
1718 read (be_unit) evec_loc
1719 read (be_unit) eval_loc
1720 do j=1,nj
1721 b = bin2d(1,j)
1722 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1723 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1724 end do
1726 case ('no3aj' )
1727 ispcc=p_no3aj-1
1728 be % v12(ispcc) % name = trim(adjustl(variable))
1729 read (be_unit) be12_evec_glo(ispcc,:,:)
1730 read (be_unit) be12_eval_glo(ispcc,:)
1731 read (be_unit) evec_loc
1732 read (be_unit) eval_loc
1733 do j=1,nj
1734 b = bin2d(1,j)
1735 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1736 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1737 end do
1739 case ('no3ai' )
1740 ispcc=p_no3ai-1
1741 be % v12(ispcc) % name = trim(adjustl(variable))
1742 read (be_unit) be12_evec_glo(ispcc,:,:)
1743 read (be_unit) be12_eval_glo(ispcc,:)
1744 read (be_unit) evec_loc
1745 read (be_unit) eval_loc
1746 do j=1,nj
1747 b = bin2d(1,j)
1748 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1749 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1750 end do
1752 case ('naaj' )
1753 ispcc=p_naaj-1
1754 be % v12(ispcc) % name = trim(adjustl(variable))
1755 read (be_unit) be12_evec_glo(ispcc,:,:)
1756 read (be_unit) be12_eval_glo(ispcc,:)
1757 read (be_unit) evec_loc
1758 read (be_unit) eval_loc
1759 do j=1,nj
1760 b = bin2d(1,j)
1761 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1762 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1763 end do
1765 case ('naai' )
1766 ispcc=p_naai-1
1767 be % v12(ispcc) % name = trim(adjustl(variable))
1768 read (be_unit) be12_evec_glo(ispcc,:,:)
1769 read (be_unit) be12_eval_glo(ispcc,:)
1770 read (be_unit) evec_loc
1771 read (be_unit) eval_loc
1772 do j=1,nj
1773 b = bin2d(1,j)
1774 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1775 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1776 end do
1778 case ('claj' )
1779 ispcc=p_claj-1
1780 be % v12(ispcc) % name = trim(adjustl(variable))
1781 read (be_unit) be12_evec_glo(ispcc,:,:)
1782 read (be_unit) be12_eval_glo(ispcc,:)
1783 read (be_unit) evec_loc
1784 read (be_unit) eval_loc
1785 do j=1,nj
1786 b = bin2d(1,j)
1787 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1788 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1789 end do
1791 case ('clai' )
1792 ispcc=p_clai-1
1793 be % v12(ispcc) % name = trim(adjustl(variable))
1794 read (be_unit) be12_evec_glo(ispcc,:,:)
1795 read (be_unit) be12_eval_glo(ispcc,:)
1796 read (be_unit) evec_loc
1797 read (be_unit) eval_loc
1798 do j=1,nj
1799 b = bin2d(1,j)
1800 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1801 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1802 end do
1804 case ('asoa1j' )
1805 ispcc=p_asoa1j-1
1806 be % v12(ispcc) % name = trim(adjustl(variable))
1807 read (be_unit) be12_evec_glo(ispcc,:,:)
1808 read (be_unit) be12_eval_glo(ispcc,:)
1809 read (be_unit) evec_loc
1810 read (be_unit) eval_loc
1811 do j=1,nj
1812 b = bin2d(1,j)
1813 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1814 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1815 end do
1817 case ('asoa1i' )
1818 ispcc=p_asoa1i-1
1819 be % v12(ispcc) % name = trim(adjustl(variable))
1820 read (be_unit) be12_evec_glo(ispcc,:,:)
1821 read (be_unit) be12_eval_glo(ispcc,:)
1822 read (be_unit) evec_loc
1823 read (be_unit) eval_loc
1824 do j=1,nj
1825 b = bin2d(1,j)
1826 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1827 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1828 end do
1830 case ('asoa2j' )
1831 ispcc=p_asoa2j-1
1832 be % v12(ispcc) % name = trim(adjustl(variable))
1833 read (be_unit) be12_evec_glo(ispcc,:,:)
1834 read (be_unit) be12_eval_glo(ispcc,:)
1835 read (be_unit) evec_loc
1836 read (be_unit) eval_loc
1837 do j=1,nj
1838 b = bin2d(1,j)
1839 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1840 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1841 end do
1843 case ('asoa2i' )
1844 ispcc=p_asoa2i-1
1845 be % v12(ispcc) % name = trim(adjustl(variable))
1846 read (be_unit) be12_evec_glo(ispcc,:,:)
1847 read (be_unit) be12_eval_glo(ispcc,:)
1848 read (be_unit) evec_loc
1849 read (be_unit) eval_loc
1850 do j=1,nj
1851 b = bin2d(1,j)
1852 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1853 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1854 end do
1856 case ('asoa3j' )
1857 ispcc=p_asoa3j-1
1858 be % v12(ispcc) % name = trim(adjustl(variable))
1859 read (be_unit) be12_evec_glo(ispcc,:,:)
1860 read (be_unit) be12_eval_glo(ispcc,:)
1861 read (be_unit) evec_loc
1862 read (be_unit) eval_loc
1863 do j=1,nj
1864 b = bin2d(1,j)
1865 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1866 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1867 end do
1869 case ('asoa3i' )
1870 ispcc=p_asoa3i-1
1871 be % v12(ispcc) % name = trim(adjustl(variable))
1872 read (be_unit) be12_evec_glo(ispcc,:,:)
1873 read (be_unit) be12_eval_glo(ispcc,:)
1874 read (be_unit) evec_loc
1875 read (be_unit) eval_loc
1876 do j=1,nj
1877 b = bin2d(1,j)
1878 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1879 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1880 end do
1882 case ('asoa4j' )
1883 ispcc=p_asoa4j-1
1884 be % v12(ispcc) % name = trim(adjustl(variable))
1885 read (be_unit) be12_evec_glo(ispcc,:,:)
1886 read (be_unit) be12_eval_glo(ispcc,:)
1887 read (be_unit) evec_loc
1888 read (be_unit) eval_loc
1889 do j=1,nj
1890 b = bin2d(1,j)
1891 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1892 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1893 end do
1895 case ('asoa4i' )
1896 ispcc=p_asoa4i-1
1897 be % v12(ispcc) % name = trim(adjustl(variable))
1898 read (be_unit) be12_evec_glo(ispcc,:,:)
1899 read (be_unit) be12_eval_glo(ispcc,:)
1900 read (be_unit) evec_loc
1901 read (be_unit) eval_loc
1902 do j=1,nj
1903 b = bin2d(1,j)
1904 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1905 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1906 end do
1908 case ('bsoa1j' )
1909 ispcc=p_bsoa1j-1
1910 be % v12(ispcc) % name = trim(adjustl(variable))
1911 read (be_unit) be12_evec_glo(ispcc,:,:)
1912 read (be_unit) be12_eval_glo(ispcc,:)
1913 read (be_unit) evec_loc
1914 read (be_unit) eval_loc
1915 do j=1,nj
1916 b = bin2d(1,j)
1917 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1918 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1919 end do
1921 case ('bsoa1i' )
1922 ispcc=p_bsoa1i-1
1923 be % v12(ispcc) % name = trim(adjustl(variable))
1924 read (be_unit) be12_evec_glo(ispcc,:,:)
1925 read (be_unit) be12_eval_glo(ispcc,:)
1926 read (be_unit) evec_loc
1927 read (be_unit) eval_loc
1928 do j=1,nj
1929 b = bin2d(1,j)
1930 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1931 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1932 end do
1934 case ('bsoa2j' )
1935 ispcc=p_bsoa2j-1
1936 be % v12(ispcc) % name = trim(adjustl(variable))
1937 read (be_unit) be12_evec_glo(ispcc,:,:)
1938 read (be_unit) be12_eval_glo(ispcc,:)
1939 read (be_unit) evec_loc
1940 read (be_unit) eval_loc
1941 do j=1,nj
1942 b = bin2d(1,j)
1943 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1944 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1945 end do
1947 case ('bsoa2i' )
1948 ispcc=p_bsoa2i-1
1949 be % v12(ispcc) % name = trim(adjustl(variable))
1950 read (be_unit) be12_evec_glo(ispcc,:,:)
1951 read (be_unit) be12_eval_glo(ispcc,:)
1952 read (be_unit) evec_loc
1953 read (be_unit) eval_loc
1954 do j=1,nj
1955 b = bin2d(1,j)
1956 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1957 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1958 end do
1960 case ('bsoa3j' )
1961 ispcc=p_bsoa3j-1
1962 be % v12(ispcc) % name = trim(adjustl(variable))
1963 read (be_unit) be12_evec_glo(ispcc,:,:)
1964 read (be_unit) be12_eval_glo(ispcc,:)
1965 read (be_unit) evec_loc
1966 read (be_unit) eval_loc
1967 do j=1,nj
1968 b = bin2d(1,j)
1969 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1970 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1971 end do
1973 case ('bsoa3i' )
1974 ispcc=p_bsoa3i-1
1975 be % v12(ispcc) % name = trim(adjustl(variable))
1976 read (be_unit) be12_evec_glo(ispcc,:,:)
1977 read (be_unit) be12_eval_glo(ispcc,:)
1978 read (be_unit) evec_loc
1979 read (be_unit) eval_loc
1980 do j=1,nj
1981 b = bin2d(1,j)
1982 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1983 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1984 end do
1986 case ('bsoa4j' )
1987 ispcc=p_bsoa4j-1
1988 be % v12(ispcc) % name = trim(adjustl(variable))
1989 read (be_unit) be12_evec_glo(ispcc,:,:)
1990 read (be_unit) be12_eval_glo(ispcc,:)
1991 read (be_unit) evec_loc
1992 read (be_unit) eval_loc
1993 do j=1,nj
1994 b = bin2d(1,j)
1995 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
1996 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
1997 end do
1999 case ('bsoa4i' )
2000 ispcc=p_bsoa4i-1
2001 be % v12(ispcc) % name = trim(adjustl(variable))
2002 read (be_unit) be12_evec_glo(ispcc,:,:)
2003 read (be_unit) be12_eval_glo(ispcc,:)
2004 read (be_unit) evec_loc
2005 read (be_unit) eval_loc
2006 do j=1,nj
2007 b = bin2d(1,j)
2008 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2009 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2010 end do
2012 case ('orgpaj' )
2013 ispcc=p_orgpaj-1
2014 be % v12(ispcc) % name = trim(adjustl(variable))
2015 read (be_unit) be12_evec_glo(ispcc,:,:)
2016 read (be_unit) be12_eval_glo(ispcc,:)
2017 read (be_unit) evec_loc
2018 read (be_unit) eval_loc
2019 do j=1,nj
2020 b = bin2d(1,j)
2021 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2022 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2023 end do
2025 case ('orgpai' )
2026 ispcc=p_orgpai-1
2027 be % v12(ispcc) % name = trim(adjustl(variable))
2028 read (be_unit) be12_evec_glo(ispcc,:,:)
2029 read (be_unit) be12_eval_glo(ispcc,:)
2030 read (be_unit) evec_loc
2031 read (be_unit) eval_loc
2032 do j=1,nj
2033 b = bin2d(1,j)
2034 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2035 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2036 end do
2038 case ('ecj' )
2039 ispcc=p_ecj-1
2040 be % v12(ispcc) % name = trim(adjustl(variable))
2041 read (be_unit) be12_evec_glo(ispcc,:,:)
2042 read (be_unit) be12_eval_glo(ispcc,:)
2043 read (be_unit) evec_loc
2044 read (be_unit) eval_loc
2045 do j=1,nj
2046 b = bin2d(1,j)
2047 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2048 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2049 end do
2051 case ('eci' )
2052 ispcc=p_eci-1
2053 be % v12(ispcc) % name = trim(adjustl(variable))
2054 read (be_unit) be12_evec_glo(ispcc,:,:)
2055 read (be_unit) be12_eval_glo(ispcc,:)
2056 read (be_unit) evec_loc
2057 read (be_unit) eval_loc
2058 do j=1,nj
2059 b = bin2d(1,j)
2060 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2061 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2062 end do
2064 case ('p25j' )
2065 ispcc=p_p25j-1
2066 be % v12(ispcc) % name = trim(adjustl(variable))
2067 read (be_unit) be12_evec_glo(ispcc,:,:)
2068 read (be_unit) be12_eval_glo(ispcc,:)
2069 read (be_unit) evec_loc
2070 read (be_unit) eval_loc
2071 do j=1,nj
2072 b = bin2d(1,j)
2073 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2074 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2075 end do
2077 case ('p25i' )
2078 ispcc=p_p25i-1
2079 be % v12(ispcc) % name = trim(adjustl(variable))
2080 read (be_unit) be12_evec_glo(ispcc,:,:)
2081 read (be_unit) be12_eval_glo(ispcc,:)
2082 read (be_unit) evec_loc
2083 read (be_unit) eval_loc
2084 do j=1,nj
2085 b = bin2d(1,j)
2086 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2087 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2088 end do
2090 case ('antha' )
2091 ispcc=p_antha-1
2092 be % v12(ispcc) % name = trim(adjustl(variable))
2093 read (be_unit) be12_evec_glo(ispcc,:,:)
2094 read (be_unit) be12_eval_glo(ispcc,:)
2095 read (be_unit) evec_loc
2096 read (be_unit) eval_loc
2097 do j=1,nj
2098 b = bin2d(1,j)
2099 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2100 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2101 end do
2103 case ('seas' )
2104 ispcc=p_seas-1
2105 be % v12(ispcc) % name = trim(adjustl(variable))
2106 read (be_unit) be12_evec_glo(ispcc,:,:)
2107 read (be_unit) be12_eval_glo(ispcc,:)
2108 read (be_unit) evec_loc
2109 read (be_unit) eval_loc
2110 do j=1,nj
2111 b = bin2d(1,j)
2112 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2113 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2114 end do
2116 case ('soila' )
2117 ispcc=p_soila-1
2118 be % v12(ispcc) % name = trim(adjustl(variable))
2119 read (be_unit) be12_evec_glo(ispcc,:,:)
2120 read (be_unit) be12_eval_glo(ispcc,:)
2121 read (be_unit) evec_loc
2122 read (be_unit) eval_loc
2123 do j=1,nj
2124 b = bin2d(1,j)
2125 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2126 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2127 end do
2129 case ('so2' )
2130 ispcc=p_so2-1
2131 be % v12(ispcc) % name = trim(adjustl(variable))
2132 read (be_unit) be12_evec_glo(ispcc,:,:)
2133 read (be_unit) be12_eval_glo(ispcc,:)
2134 read (be_unit) evec_loc
2135 read (be_unit) eval_loc
2136 do j=1,nj
2137 b = bin2d(1,j)
2138 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2139 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2140 end do
2142 case ('no2' )
2143 ispcc=p_no2-1
2144 be % v12(ispcc) % name = trim(adjustl(variable))
2145 read (be_unit) be12_evec_glo(ispcc,:,:)
2146 read (be_unit) be12_eval_glo(ispcc,:)
2147 read (be_unit) evec_loc
2148 read (be_unit) eval_loc
2149 do j=1,nj
2150 b = bin2d(1,j)
2151 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2152 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2153 end do
2155 case ('o3' )
2156 ispcc=p_o3-1
2157 be % v12(ispcc) % name = trim(adjustl(variable))
2158 read (be_unit) be12_evec_glo(ispcc,:,:)
2159 read (be_unit) be12_eval_glo(ispcc,:)
2160 read (be_unit) evec_loc
2161 read (be_unit) eval_loc
2162 do j=1,nj
2163 b = bin2d(1,j)
2164 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2165 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2166 end do
2168 case ('co' )
2169 ispcc=p_co-1
2170 be % v12(ispcc) % name = trim(adjustl(variable))
2171 read (be_unit) be12_evec_glo(ispcc,:,:)
2172 read (be_unit) be12_eval_glo(ispcc,:)
2173 read (be_unit) evec_loc
2174 read (be_unit) eval_loc
2175 do j=1,nj
2176 b = bin2d(1,j)
2177 be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2178 be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b)
2179 end do
2181 case ('w' )
2182 if ( use_cv_w ) then
2183 be % v11 % name = trim(adjustl(variable))
2184 read (be_unit) be11_evec_glo
2185 read (be_unit) be11_eval_glo
2186 read (be_unit) evec_loc
2187 read (be_unit) eval_loc
2188 do j=1,nj
2189 b = bin2d(1,j)
2190 be11_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
2191 be11_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
2192 end do
2193 end if
2195 case default;
2196 message(1)=' Read problem in eigen vectors/values in BE file '
2197 write (unit=message(2),fmt='(A,A)') ' Trying to read Eigenvectors for variable: ',trim(adjustl(variable))
2198 write (unit=message(3),fmt='(A)') ' Make sure you are using the correct be.dat file for your cv_options setting!'
2199 call da_error(__FILE__,__LINE__,message(1:3))
2200 end select
2202 !print '(A6,3(A,I4),I5)', trim(adjustl(variable)), ' i =', i, ' , ispcc =', ispcc,' , nk =', nk, num_bins2d
2203 end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra-1
2204 print*,"reading variables for chem_cv_options=108 (racm_soa_vbs_da) is done"
2205 end if ! ( chem_cv_options == 108 ) then ! racm_soa_vbs_da
2206 #endif
2208 deallocate (evec_loc)
2209 deallocate (eval_loc)
2211 ! 2.2.5 Read in 2D Control variable ps_u
2213 read (be_unit) variable
2214 read (be_unit) nk_2d, num_bins2d
2215 print *, trim(adjustl(variable)), nk_2d, num_bins2d
2217 !hcl-why
2218 !#ifdef CLOUD_CV
2219 ! nk_2d=1
2220 !#endif
2221 allocate (evec_loc(1:nk_2d,1:nk_2d,1:num_bins2d))
2222 allocate (eval_loc(1:nk_2d, 1:num_bins2d))
2223 read (be_unit) be5_evec_glo
2224 read (be_unit) be5_eval_glo
2225 read (be_unit) evec_loc
2226 read (be_unit) eval_loc
2227 if( (trim(adjustl(variable)) /= 'ps_u') .and. (trim(adjustl(variable)) /= 'ps') ) then
2228 message(1)=' Variable mismatch while transfering eigen values from BE file '
2229 if (cv_options == 7) then
2230 write (unit=message(2),fmt='(A,A)') ' Expected ps but got ',trim(adjustl(variable))
2231 else
2232 write (unit=message(2),fmt='(A,A)') ' Expected ps_u but got ',trim(adjustl(variable))
2233 endif
2234 call da_error(__FILE__,__LINE__,message(1:2))
2235 end if
2236 be % v5 % name = trim(adjustl(variable))
2237 do j=1,nj
2238 b = bin2d(1,j)
2239 be5_evec_loc(j,1:1,1:1) = evec_loc(1:1,1:1,b)
2240 be5_eval_loc(j,1:1 ) = eval_loc(1:1,b)
2241 end do
2243 deallocate (evec_loc)
2244 deallocate (eval_loc)
2246 if(use_radarobs .and. use_radar_rf .or. use_rad .and. crtm_cloud) then
2247 if ( cloud_cv_options == 1 ) be % v4 % name = 'qt '
2248 end if
2250 message(1:5) = ''
2251 write (unit=message(1),fmt='(3x,A,i2)') 'cloud_cv_options = ', cloud_cv_options
2252 #if (WRF_CHEM == 1)
2253 write (unit=message(1),fmt='(3x,A,i3)') 'chem_cv_options = ', chem_cv_options
2254 #endif
2255 write (unit=message(2),fmt='(3x,A,7A)') 'WRFDA dry control variables are: ', &
2256 trim(be % v1 % name), ', ', trim(be % v2 % name), ', ', &
2257 trim(be % v3 % name), ' and ', trim(be % v5 % name)
2258 write (unit=message(3),fmt='(3x,A,A)') &
2259 'WRFDA Humidity control variable is ', trim(be % v4 % name)
2260 if ( cloud_cv_options >= 2 ) then
2261 write (unit=message(4),fmt='(3x,A,10A)') 'WRFDA hydrometeor control variables are: ', &
2262 trim(be % v6 % name), ', ', trim(be % v7 % name), ', ', &
2263 trim(be % v8 % name), ', ', trim(be % v9 % name), ', and ', &
2264 trim(be % v10 % name)
2265 end if
2266 !#if (WRF_CHEM == 1)
2267 ! if ( chem_cv_options >= 10 ) then
2268 ! write (unit=message(4),fmt='(3x,A,10A)') 'WRFDA chem control variables are: ', &
2269 ! be % v12(:) % name
2270 ! end if
2271 !#endif
2272 if ( use_cv_w ) then
2273 write (unit=message(5),fmt='(3x,A)') &
2274 'w control variable is activated'
2275 end if
2276 call da_message(message(1:5))
2278 if (use_rf) then
2280 ! 3.0 Load the scale lengths
2282 ! 3.1 allocate the array for scale lengths
2283 allocate (rfls1_be(1:nk))
2284 allocate (rfls2_be(1:nk))
2285 allocate (rfls3_be(1:nk))
2286 allocate (rfls4_be(1:nk))
2288 allocate (be1_rf_lengthscale(1:kz))
2289 allocate (be2_rf_lengthscale(1:kz))
2290 allocate (be3_rf_lengthscale(1:kz))
2291 allocate (be4_rf_lengthscale(1:kz))
2292 allocate (be5_rf_lengthscale(1:1))
2295 ! ===================================================================
2296 ! This is used for do_normalize = .T., I am not sure if it is working
2297 ! properly (???), more testing needed. (YRG, 01/12/2012)
2298 !
2299 b=0 ! pointer to variable*mode
2300 allocate(nij(0:0,0:1,0:0))
2302 kzs=(/kz,kz,kz,kz,1/)
2303 be%v1%mz=1; be%v2%mz=1;be%v3%mz=1;be%v4%mz=1;be%v5%mz=1;
2304 mzs=(/be%v1%mz,be%v2%mz,be%v3%mz,be%v4%mz,be%v5%mz/)
2305 ! ====================================================================
2307 ! 3.2 read in the scale lengths
2309 print *, '----- read lengthscale --------'
2310 do i = 1 , num_cv_3d_basic
2311 read (be_unit) variable
2312 print *, trim(adjustl(variable))
2313 select case( trim(adjustl(variable)) )
2314 case ('psi', 'u')
2315 read(be_unit) rfls1_be
2316 if (.not.interpolate_stats) be1_rf_lengthscale = rfls1_be
2317 case ('chi_u', 'v')
2318 read(be_unit) rfls2_be
2319 if (.not.interpolate_stats) be2_rf_lengthscale = rfls2_be
2320 case ('t_u', 't')
2321 read(be_unit) rfls3_be
2322 if (.not.interpolate_stats) be3_rf_lengthscale = rfls3_be
2323 case ('rh_u' , 'rh')
2324 read(be_unit) rfls4_be
2325 if (.not.interpolate_stats) be4_rf_lengthscale = rfls4_be
2326 case default;
2327 message(1)='Read problem in lengthscales in be.dat'
2328 write(message(2),'("Trying to read lengthscales for variable ",I0,": ",A)')i,trim(adjustl(variable))
2329 call da_error(__FILE__,__LINE__,message(1:2))
2330 end select
2331 end do ! num_cv_3d_basic
2333 if ( cloud_cv_options == 2 ) then
2334 do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra
2335 read (be_unit) variable
2336 print *, trim(adjustl(variable))
2337 select case( trim(adjustl(variable)) )
2338 case ('qcloud')
2339 read(be_unit) be6_rf_lengthscale
2340 case ('qrain')
2341 read(be_unit) be7_rf_lengthscale
2342 case ('qice')
2343 read(be_unit) be8_rf_lengthscale
2344 case ('qsnow')
2345 read(be_unit) be9_rf_lengthscale
2346 case ('qgraup')
2347 read(be_unit) be10_rf_lengthscale
2348 case ('w')
2349 !if ( use_cv_w ) then
2350 ! if ( trim(adjustl(variable)) == 'w' ) then
2351 read(be_unit) be11_rf_lengthscale
2352 ! end if
2353 !end if
2354 case default;
2355 message(1)='Read problem in lengthscales in be.dat'
2356 write(message(2),'("Trying to read lengthscales for variable ",I0,": ",A)')i,trim(adjustl(variable))
2357 call da_error(__FILE__,__LINE__,message(1:2))
2358 end select
2359 end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra
2360 end if
2362 #if (WRF_CHEM == 1)
2363 if ( chem_cv_options == 10 ) then
2364 do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra
2365 read (be_unit) variable
2366 print *, trim(adjustl(variable))
2367 select case( trim(adjustl(variable)) )
2368 case ('dust_1')
2369 ispcc=p_dust_1-1
2370 read(be_unit) be12_rf_lengthscale(ispcc,:)
2371 case ('dust_2')
2372 ispcc=p_dust_2-1
2373 read(be_unit) be12_rf_lengthscale(ispcc,:)
2374 case ('dust_3')
2375 ispcc=p_dust_3-1
2376 read(be_unit) be12_rf_lengthscale(ispcc,:)
2377 case ('dust_4')
2378 ispcc=p_dust_4-1
2379 read(be_unit) be12_rf_lengthscale(ispcc,:)
2381 case ('seas_1')
2382 ispcc=p_seas_1-1
2383 read(be_unit) be12_rf_lengthscale(ispcc,:)
2384 case ('seas_2')
2385 ispcc=p_seas_2-1
2386 read(be_unit) be12_rf_lengthscale(ispcc,:)
2387 case ('seas_3')
2388 ispcc=p_seas_3-1
2389 read(be_unit) be12_rf_lengthscale(ispcc,:)
2390 case ('seas_4')
2391 ispcc=p_seas_4-1
2392 read(be_unit) be12_rf_lengthscale(ispcc,:)
2394 case ('oc1')
2395 ispcc=p_oc1-1
2396 read(be_unit) be12_rf_lengthscale(ispcc,:)
2397 case ('oc2')
2398 ispcc=p_oc2-1
2399 read(be_unit) be12_rf_lengthscale(ispcc,:)
2400 case ('bc1')
2401 ispcc=p_bc1-1
2402 read(be_unit) be12_rf_lengthscale(ispcc,:)
2403 case ('bc2')
2404 ispcc=p_bc2-1
2405 read(be_unit) be12_rf_lengthscale(ispcc,:)
2407 case ('so2')
2408 ispcc=p_so2-1
2409 read(be_unit) be12_rf_lengthscale(ispcc,:)
2410 case ('no2')
2411 ispcc=p_no2-1
2412 read(be_unit) be12_rf_lengthscale(ispcc,:)
2413 case ('o3')
2414 ispcc=p_o3-1
2415 read(be_unit) be12_rf_lengthscale(ispcc,:)
2416 case ('co')
2417 ispcc=p_co-1
2418 read(be_unit) be12_rf_lengthscale(ispcc,:)
2419 case ('sulf')
2420 ispcc=p_sulf-1
2421 read(be_unit) be12_rf_lengthscale(ispcc,:)
2422 case ('p25')
2423 ispcc=p_p25-1
2424 read(be_unit) be12_rf_lengthscale(ispcc,:)
2425 case ('p10')
2426 ispcc=p_p10-1
2427 read(be_unit) be12_rf_lengthscale(ispcc,:)
2429 case ('w')
2430 !if ( use_cv_w ) then
2431 ! if ( trim(adjustl(variable)) == 'w' ) then
2432 read(be_unit) be11_rf_lengthscale
2433 ! end if
2434 !end if
2435 case default;
2436 message(1)='Read problem in lengthscales in be.dat'
2437 write(message(2),'("Trying to read lengthscales for variable ",I0,": ",A)')i,trim(adjustl(variable))
2438 call da_error(__FILE__,__LINE__,message(1:2))
2439 end select
2440 end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra
2441 end if
2443 if ( chem_cv_options == 20 ) then
2444 do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra
2445 read (be_unit) variable
2446 print *, trim(adjustl(variable))
2447 select case( trim(adjustl(variable)) )
2448 case ('bc_1')
2449 ispcc=p_bc_a01-1
2450 read(be_unit) be12_rf_lengthscale(ispcc,:)
2451 case ('bc_2')
2452 ispcc=p_bc_a02-1
2453 read(be_unit) be12_rf_lengthscale(ispcc,:)
2454 case ('bc_3')
2455 ispcc=p_bc_a03-1
2456 read(be_unit) be12_rf_lengthscale(ispcc,:)
2457 case ('bc_4')
2458 ispcc=p_bc_a04-1
2459 read(be_unit) be12_rf_lengthscale(ispcc,:)
2461 case ('oc_1')
2462 ispcc=p_oc_a01-1
2463 read(be_unit) be12_rf_lengthscale(ispcc,:)
2464 case ('oc_2')
2465 ispcc=p_oc_a02-1
2466 read(be_unit) be12_rf_lengthscale(ispcc,:)
2467 case ('oc_3')
2468 ispcc=p_oc_a03-1
2469 read(be_unit) be12_rf_lengthscale(ispcc,:)
2470 case ('oc_4')
2471 ispcc=p_oc_a04-1
2472 read(be_unit) be12_rf_lengthscale(ispcc,:)
2474 case ('so4_1')
2475 ispcc=p_so4_a01-1
2476 read(be_unit) be12_rf_lengthscale(ispcc,:)
2477 case ('so4_2')
2478 ispcc=p_so4_a02-1
2479 read(be_unit) be12_rf_lengthscale(ispcc,:)
2480 case ('so4_3')
2481 ispcc=p_so4_a03-1
2482 read(be_unit) be12_rf_lengthscale(ispcc,:)
2483 case ('so4_4')
2484 ispcc=p_so4_a04-1
2485 read(be_unit) be12_rf_lengthscale(ispcc,:)
2487 case ('no3_1')
2488 ispcc=p_no3_a01-1
2489 read(be_unit) be12_rf_lengthscale(ispcc,:)
2490 case ('no3_2')
2491 ispcc=p_no3_a02-1
2492 read(be_unit) be12_rf_lengthscale(ispcc,:)
2493 case ('no3_3')
2494 ispcc=p_no3_a03-1
2495 read(be_unit) be12_rf_lengthscale(ispcc,:)
2496 case ('no3_4')
2497 ispcc=p_no3_a04-1
2498 read(be_unit) be12_rf_lengthscale(ispcc,:)
2500 case ('nh4_1')
2501 ispcc=p_nh4_a01-1
2502 read(be_unit) be12_rf_lengthscale(ispcc,:)
2503 case ('nh4_2')
2504 ispcc=p_nh4_a02-1
2505 read(be_unit) be12_rf_lengthscale(ispcc,:)
2506 case ('nh4_3')
2507 ispcc=p_nh4_a03-1
2508 read(be_unit) be12_rf_lengthscale(ispcc,:)
2509 case ('nh4_4')
2510 ispcc=p_nh4_a04-1
2511 read(be_unit) be12_rf_lengthscale(ispcc,:)
2513 case ('cl_1')
2514 ispcc=p_cl_a01-1
2515 read(be_unit) be12_rf_lengthscale(ispcc,:)
2516 case ('cl_2')
2517 ispcc=p_cl_a02-1
2518 read(be_unit) be12_rf_lengthscale(ispcc,:)
2519 case ('cl_3')
2520 ispcc=p_cl_a03-1
2521 read(be_unit) be12_rf_lengthscale(ispcc,:)
2522 case ('cl_4')
2523 ispcc=p_cl_a04-1
2524 read(be_unit) be12_rf_lengthscale(ispcc,:)
2526 case ('na_1')
2527 ispcc=p_na_a01-1
2528 read(be_unit) be12_rf_lengthscale(ispcc,:)
2529 case ('na_2')
2530 ispcc=p_na_a02-1
2531 read(be_unit) be12_rf_lengthscale(ispcc,:)
2532 case ('na_3')
2533 ispcc=p_na_a03-1
2534 read(be_unit) be12_rf_lengthscale(ispcc,:)
2535 case ('na_4')
2536 ispcc=p_na_a04-1
2537 read(be_unit) be12_rf_lengthscale(ispcc,:)
2539 case ('oin_1')
2540 ispcc=p_oin_a01-1
2541 read(be_unit) be12_rf_lengthscale(ispcc,:)
2542 case ('oin_2')
2543 ispcc=p_oin_a02-1
2544 read(be_unit) be12_rf_lengthscale(ispcc,:)
2545 case ('oin_3')
2546 ispcc=p_oin_a03-1
2547 read(be_unit) be12_rf_lengthscale(ispcc,:)
2548 case ('oin_4')
2549 ispcc=p_oin_a04-1
2550 read(be_unit) be12_rf_lengthscale(ispcc,:)
2552 case ('so2')
2553 ispcc=p_so2-1
2554 read(be_unit) be12_rf_lengthscale(ispcc,:)
2555 case ('no2')
2556 ispcc=p_no2-1
2557 read(be_unit) be12_rf_lengthscale(ispcc,:)
2558 case ('o3')
2559 ispcc=p_o3-1
2560 read(be_unit) be12_rf_lengthscale(ispcc,:)
2561 case ('co')
2562 ispcc=p_co-1
2563 read(be_unit) be12_rf_lengthscale(ispcc,:)
2565 case ('w')
2566 if ( use_cv_w ) then
2567 read(be_unit) be11_rf_lengthscale
2568 end if
2569 case default;
2570 message(1)='Read problem in lengthscales in be.dat'
2571 write(message(2),'("Trying to read lengthscales for variable ",I0,": ",A)')i,trim(adjustl(variable))
2572 call da_error(__FILE__,__LINE__,message(1:2))
2573 end select
2574 end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra
2575 end if ! if ( chem_cv_options == 20 ) then
2577 if ( chem_cv_options == 108 ) then ! racm_soa_vbs_da
2578 do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra
2579 read (be_unit) variable
2580 print *, trim(adjustl(variable))
2581 select case( trim(adjustl(variable)) )
2582 case ('so4aj')
2583 ispcc=p_so4aj-1
2584 read(be_unit) be12_rf_lengthscale(ispcc,:)
2585 case ('so4ai' )
2586 ispcc=p_so4ai-1
2587 read(be_unit) be12_rf_lengthscale(ispcc,:)
2588 case ('nh4aj' )
2589 ispcc=p_nh4aj-1
2590 read(be_unit) be12_rf_lengthscale(ispcc,:)
2591 case ('nh4ai' )
2592 ispcc=p_nh4ai-1
2593 read(be_unit) be12_rf_lengthscale(ispcc,:)
2594 case ('no3aj' )
2595 ispcc=p_no3aj-1
2596 read(be_unit) be12_rf_lengthscale(ispcc,:)
2597 case ('no3ai' )
2598 ispcc=p_no3ai-1
2599 read(be_unit) be12_rf_lengthscale(ispcc,:)
2600 case ('naaj' )
2601 ispcc=p_naaj-1
2602 read(be_unit) be12_rf_lengthscale(ispcc,:)
2603 case ('naai' )
2604 ispcc=p_naai-1
2605 read(be_unit) be12_rf_lengthscale(ispcc,:)
2606 case ('claj' )
2607 ispcc=p_claj-1
2608 read(be_unit) be12_rf_lengthscale(ispcc,:)
2609 case ('clai' )
2610 ispcc=p_clai-1
2611 read(be_unit) be12_rf_lengthscale(ispcc,:)
2612 case ('asoa1j' )
2613 ispcc=p_asoa1j-1
2614 read(be_unit) be12_rf_lengthscale(ispcc,:)
2615 case ('asoa1i' )
2616 ispcc=p_asoa1i-1
2617 read(be_unit) be12_rf_lengthscale(ispcc,:)
2618 case ('asoa2j' )
2619 ispcc=p_asoa2j-1
2620 read(be_unit) be12_rf_lengthscale(ispcc,:)
2621 case ('asoa2i' )
2622 ispcc=p_asoa2i-1
2623 read(be_unit) be12_rf_lengthscale(ispcc,:)
2624 case ('asoa3j' )
2625 ispcc=p_asoa3j-1
2626 read(be_unit) be12_rf_lengthscale(ispcc,:)
2627 case ('asoa3i' )
2628 ispcc=p_asoa3i-1
2629 read(be_unit) be12_rf_lengthscale(ispcc,:)
2630 case ('asoa4j' )
2631 ispcc=p_asoa4j-1
2632 read(be_unit) be12_rf_lengthscale(ispcc,:)
2633 case ('asoa4i' )
2634 ispcc=p_asoa4i-1
2635 read(be_unit) be12_rf_lengthscale(ispcc,:)
2636 case ('bsoa1j' )
2637 ispcc=p_bsoa1j-1
2638 read(be_unit) be12_rf_lengthscale(ispcc,:)
2639 case ('bsoa1i' )
2640 ispcc=p_bsoa1i-1
2641 read(be_unit) be12_rf_lengthscale(ispcc,:)
2642 case ('bsoa2j' )
2643 ispcc=p_bsoa2j-1
2644 read(be_unit) be12_rf_lengthscale(ispcc,:)
2645 case ('bsoa2i' )
2646 ispcc=p_bsoa2i-1
2647 read(be_unit) be12_rf_lengthscale(ispcc,:)
2648 case ('bsoa3j' )
2649 ispcc=p_bsoa3j-1
2650 read(be_unit) be12_rf_lengthscale(ispcc,:)
2651 case ('bsoa3i' )
2652 ispcc=p_bsoa3i-1
2653 read(be_unit) be12_rf_lengthscale(ispcc,:)
2654 case ('bsoa4j' )
2655 ispcc=p_bsoa4j-1
2656 read(be_unit) be12_rf_lengthscale(ispcc,:)
2657 case ('bsoa4i' )
2658 ispcc=p_bsoa4i-1
2659 read(be_unit) be12_rf_lengthscale(ispcc,:)
2660 case ('orgpaj' )
2661 ispcc=p_orgpaj-1
2662 read(be_unit) be12_rf_lengthscale(ispcc,:)
2663 case ('orgpai' )
2664 ispcc=p_orgpai-1
2665 read(be_unit) be12_rf_lengthscale(ispcc,:)
2666 case ('ecj' )
2667 ispcc=p_ecj-1
2668 read(be_unit) be12_rf_lengthscale(ispcc,:)
2669 case ('eci' )
2670 ispcc=p_eci-1
2671 read(be_unit) be12_rf_lengthscale(ispcc,:)
2672 case ('p25j' )
2673 ispcc=p_p25j-1
2674 read(be_unit) be12_rf_lengthscale(ispcc,:)
2675 case ('p25i' )
2676 ispcc=p_p25i-1
2677 read(be_unit) be12_rf_lengthscale(ispcc,:)
2678 case ('antha' )
2679 ispcc=p_antha-1
2680 read(be_unit) be12_rf_lengthscale(ispcc,:)
2681 case ('seas' )
2682 ispcc=p_seas-1
2683 read(be_unit) be12_rf_lengthscale(ispcc,:)
2684 case ('soila' )
2685 ispcc=p_soila-1
2686 read(be_unit) be12_rf_lengthscale(ispcc,:)
2687 case ('so2' )
2688 ispcc=p_so2-1
2689 read(be_unit) be12_rf_lengthscale(ispcc,:)
2690 case ('no2' )
2691 ispcc=p_no2-1
2692 read(be_unit) be12_rf_lengthscale(ispcc,:)
2693 case ('o3' )
2694 ispcc=p_o3-1
2695 read(be_unit) be12_rf_lengthscale(ispcc,:)
2696 case ('co' )
2697 ispcc=p_co-1
2698 read(be_unit) be12_rf_lengthscale(ispcc,:)
2699 case ('w')
2700 if ( use_cv_w ) then
2701 read(be_unit) be11_rf_lengthscale
2702 end if
2703 case default;
2704 message(1)='Read problem in lengthscales in be.dat'
2705 write(message(2),'("Trying to read lengthscales for variable ",I0,": ",A)')i,trim(adjustl(variable))
2706 call da_error(__FILE__,__LINE__,message(1:2))
2707 end select
2708 end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra
2709 endif !( chem_cv_options == 108 ) then ! racm_soa_vbs_da
2710 #endif
2712 ! Read in lengthscale of 2D Control variable ps_u
2714 read (be_unit) variable
2715 print *, trim(adjustl(variable))
2716 if ( trim(adjustl(variable)) /= 'ps_u' .and. &
2717 trim(adjustl(variable)) /= 'ps' ) then
2718 message(1)='Read problem in lengthscales in be.dat'
2719 write(message(2),'("Trying to read lengthscales for variable : ",A)') trim(adjustl(variable))
2720 call da_error(__FILE__,__LINE__,message(1:2))
2721 end if
2722 read (be_unit) be5_rf_lengthscale(1:1)
2724 if( do_normalize )then
2725 read(be_unit)do_normalize1
2726 if( do_normalize.neqv.do_normalize1 ) &
2727 call da_error(__FILE__,__LINE__,(/namv(i)//": do_normalize.neqv.do_normalize1"/))
2728 read(be_unit)nij
2729 if( i==1 )allocate(be%sd(nij(0,1,0),nij(0,0,0),sum(mzs)))
2730 do k=1,mzs(i)
2731 read(be_unit)be%sd(:,:,b+k)
2732 enddo
2733 do k=mzs(i)+1,kzs(i)! read and discard truncated modes:
2734 read(be_unit)(junk,j=1,nij(0,1,0)*nij(0,0,0))
2735 enddo
2736 write(*,'(A,": |sd[",A,"]|=",es9.3)')__FILE__,namv(i),sqrt(sum(be%sd(:,:,b+1:b+mzs(i))**2))
2737 b=b+mzs(i) ! point to next variable.
2738 endif
2740 if (print_detail_be) then
2741 ! Lengthscale variables (except 5!) are arrays of length "kz", so need to account for this in fmt statement
2742 write(write_fmt,'(A,I0,A)')'(A,',kz,'F10.5)'
2743 write(unit=message(1),fmt='(A)') 'BE Length scales:'
2744 write(unit=message(2),fmt=write_fmt) 'be1= ',be1_rf_lengthscale
2745 write(unit=message(3),fmt=write_fmt) 'be2= ',be2_rf_lengthscale
2746 write(unit=message(4),fmt=write_fmt) 'be3= ',be3_rf_lengthscale
2747 write(unit=message(5),fmt=write_fmt) 'be4= ',be4_rf_lengthscale
2748 write(unit=message(6),fmt='(A,F10.5)') 'be5= ',be5_rf_lengthscale
2749 call da_message(message(1:6))
2751 if ( cloud_cv_options >= 2 ) then
2752 write(unit=message(7),fmt=write_fmt) 'be6= ',be6_rf_lengthscale
2753 write(unit=message(8),fmt=write_fmt) 'be7= ',be7_rf_lengthscale
2754 write(unit=message(9),fmt=write_fmt) 'be8= ',be8_rf_lengthscale
2755 write(unit=message(10),fmt=write_fmt) 'be9= ',be9_rf_lengthscale
2756 write(unit=message(11),fmt=write_fmt) 'be10= ',be10_rf_lengthscale
2757 call da_message(message(7:11))
2758 end if
2759 !#if (WRF_CHEM == 1)
2760 ! if ( chem_cv_options >= 10 ) then
2761 ! do i=1,num_cv_3d_chem
2762 ! write(unit=message(i+6),fmt=write_fmt) 'be12 for ichem ',(i+5),'= ',be12_rf_lengthscale(i,:)
2763 ! end do
2764 ! call da_message(message(i+6:num_cv_3d_chem+6))
2765 ! end if
2766 !#endif
2767 if ( use_cv_w ) then
2768 write(unit=message(12),fmt=write_fmt) 'be11= ',be11_rf_lengthscale
2769 call da_message(message(12:12))
2770 end if
2771 end if ! print_detail_be
2773 end if ! use_rf
2774 !================== for interpolating CV5 =======================================
2775 ! 3.3 Interpolate the CV5 BE to the model resolution
2777 if (interpolate_stats) then
2779 ! 3.3.1 Compute the Vertical covariance matrix arrays for input CV5 BE:
2781 ! 3.3.1.1 Allocate the covariance arrays for input CV5 BE
2783 allocate (covm1_be(1:nk,1:nk))
2784 allocate (covm2_be(1:nk,1:nk))
2785 allocate (covm3_be(1:nk,1:nk))
2786 allocate (covm4_be(1:nk,1:nk))
2787 ! 3.3.1.2 Compute the vertical covariance matrix from Input CV5 BE
2788 ! eigenvector/eigenvalue:
2789 call da_Eigen_to_covmatrix(nk, be1_evec_glo, be1_eval_glo, covm1_be)
2790 call da_Eigen_to_covmatrix(nk, be2_evec_glo, be2_eval_glo, covm2_be)
2791 call da_Eigen_to_covmatrix(nk, be3_evec_glo, be3_eval_glo, covm3_be)
2792 call da_Eigen_to_covmatrix(nk, be4_evec_glo, be4_eval_glo, covm4_be)
2793 ! Deallocate the input CV5 BE Global eigenvector/eigenvalue:
2794 DEALLOCATE ( be1_eval_glo )
2795 DEALLOCATE ( be2_eval_glo )
2796 DEALLOCATE ( be3_eval_glo )
2797 DEALLOCATE ( be4_eval_glo )
2799 DEALLOCATE ( be1_evec_glo )
2800 DEALLOCATE ( be2_evec_glo )
2801 DEALLOCATE ( be3_evec_glo )
2802 DEALLOCATE ( be4_evec_glo )
2804 ! Allocate the model CV5 BE Global eigenvector/eigenvalue:
2805 ALLOCATE ( be1_eval_glo(1:kz) )
2806 ALLOCATE ( be2_eval_glo(1:kz) )
2807 ALLOCATE ( be3_eval_glo(1:kz) )
2808 ALLOCATE ( be4_eval_glo(1:kz) )
2810 ALLOCATE ( be1_evec_glo(1:kz,1:kz) )
2811 ALLOCATE ( be2_evec_glo(1:kz,1:kz) )
2812 ALLOCATE ( be3_evec_glo(1:kz,1:kz) )
2813 ALLOCATE ( be4_evec_glo(1:kz,1:kz) )
2815 ! 3.3.2 Convert the Vertical resolution from input CV5 BE to Model:
2817 ! 3.3.2.1 Allocate the Covariance matrix arrays for model
2819 allocate (covm1(1:kz,1:kz))
2820 allocate (covm2(1:kz,1:kz))
2821 allocate (covm3(1:kz,1:kz))
2822 allocate (covm4(1:kz,1:kz))
2824 ! 3.3.2.2 Allocate the regression coeff. arrays for model BE
2826 allocate (reg_psi_ps (1:kz))
2827 allocate (reg_psi_chi(1:kz))
2828 allocate (reg_psi_t (1:kz,1:kz))
2829 ! 3.3.2 Vertical resolution conversion from input BE(nk) to model(kz) BE:
2830 call da_chg_be_vres(kz, nk, xb%sigmah, eta_be,&
2831 reg_psi_chi_be,reg_psi_ps_be,reg_psi_t_be, &
2832 reg_psi_chi ,reg_psi_ps ,reg_psi_t , &
2833 covm1_be, covm2_be, covm3_be, covm4_be, &
2834 covm1 , covm2 , covm3 , covm4 , &
2835 rfls1_be, rfls2_be, rfls3_be, rfls4_be, &
2836 be1_rf_lengthscale, be2_rf_lengthscale, &
2837 be3_rf_lengthscale, be4_rf_lengthscale)
2839 ! 3.3.3 Deallocate the arrays for input CV5 BE:
2840 deallocate (reg_psi_ps_be )
2841 deallocate (reg_psi_chi_be)
2842 deallocate (reg_psi_t_be )
2844 ! 3.3.4 Assign the Regression coefficients for model BE:
2846 ! 3.3.4.1 Ps and Chi regression coefficients
2847 do k=1,kz
2848 do j=1,jy
2849 be%reg_psi_ps(j,k) = reg_psi_ps(k)
2850 be%reg_psi_chi (j,k)= reg_psi_chi(k)
2851 enddo
2852 enddo
2853 ! write(6,*)'be%reg_psi_ps(10,10)=',be%reg_psi_ps(10,10)
2854 ! write(6,*)'be%reg_psi_chi(10,10)=',be%reg_psi_chi(10,10)
2855 ! 3.3.4.2 Temp. regression coefficients
2856 do m=1,kz
2857 do k=1,kz
2858 do j=1,jy
2859 be%reg_psi_t(j,k,m) = reg_psi_t(k,m)
2860 enddo
2861 enddo
2862 enddo
2863 ! write(6,*)'be%reg_psi_t(10,10,10)=',be%reg_psi_t(10,10,10)
2864 ! 3.3.4.3 RWite out the Interpolated BE: Regression coefficients
2865 if (rootproc) then
2866 write(be_out_unit) reg_psi_chi
2867 write(be_out_unit) reg_psi_ps
2868 write(be_out_unit) reg_psi_t
2869 endif
2870 ! 3.3.4.4 Deallocate the domain-averaged Reg. Coeff. arrays:
2872 deallocate (reg_psi_ps )
2873 deallocate (reg_psi_chi)
2874 deallocate (reg_psi_t )
2875 ! 3.3.5 Re-compute the Eigenvector/eigenvalue from the model covariance matrix
2876 call da_gen_eigen(kz, covm1, be1_evec_glo, be1_eval_glo)
2877 call da_gen_eigen(kz, covm2, be2_evec_glo, be2_eval_glo)
2878 call da_gen_eigen(kz, covm3, be3_evec_glo, be3_eval_glo)
2879 call da_gen_eigen(kz, covm4, be4_evec_glo, be4_eval_glo)
2881 deallocate (covm1)
2882 deallocate (covm2)
2883 deallocate (covm3)
2884 deallocate (covm4)
2885 deallocate (covm1_be)
2886 deallocate (covm2_be)
2887 deallocate (covm3_be)
2888 deallocate (covm4_be)
2890 ! 3.3.6 Write out the interpolated CV5 BE
2892 ! 3.3.6.1 Eigenvector/eigenvalue:
2893 variable = be % v1 % name
2894 if (rootproc) then
2895 write(be_out_unit) variable
2896 write(be_out_unit) kz, num_bins2d_out
2897 write(be_out_unit) be1_evec_glo
2898 write(be_out_unit) be1_eval_glo
2899 ! - Because num_bins2d = 1, be1_evec_loc = be1_evec_glo:
2900 write(be_out_unit) be1_evec_glo
2901 write(be_out_unit) be1_eval_glo
2902 endif
2903 variable = be % v2 % name
2904 if (rootproc) then
2905 write(be_out_unit) variable
2906 write(be_out_unit) kz, num_bins2d_out
2907 write(be_out_unit) be2_evec_glo
2908 write(be_out_unit) be2_eval_glo
2909 ! - Because num_bins2d = 1, be2_evec_loc = be2_evec_glo:
2910 write(be_out_unit) be2_evec_glo
2911 write(be_out_unit) be2_eval_glo
2912 endif
2913 variable = be % v3 % name
2914 if (rootproc) then
2915 write(be_out_unit) variable
2916 write(be_out_unit) kz, num_bins2d_out
2917 write(be_out_unit) be3_evec_glo
2918 write(be_out_unit) be3_eval_glo
2919 ! - Because num_bins2d = 1, be3_evec_loc = be3_evec_glo:
2920 write(be_out_unit) be3_evec_glo
2921 write(be_out_unit) be3_eval_glo
2922 endif
2923 variable = be % v4 % name
2924 if (rootproc) then
2925 write(be_out_unit) variable
2926 write(be_out_unit) kz, num_bins2d_out
2927 write(be_out_unit) be4_evec_glo
2928 write(be_out_unit) be4_eval_glo
2929 ! - Because num_bins2d = 1, be4_evec_loc = be4_evec_glo:
2930 write(be_out_unit) be4_evec_glo
2931 write(be_out_unit) be4_eval_glo
2932 endif
2933 variable = be % v5 % name
2934 if (rootproc) then
2935 write(be_out_unit) variable
2936 write(be_out_unit) nk_2d, num_bins2d_out
2937 write(be_out_unit) be5_evec_glo
2938 write(be_out_unit) be5_eval_glo
2939 ! - Because num_bins2d = 1, be5_evec_loc = be5_evec_glo:
2940 write(be_out_unit) be5_evec_glo
2941 write(be_out_unit) be5_eval_glo
2942 endif
2944 ! 3.3.6.2 Scale-length
2945 variable = be % v1 % name
2946 if (rootproc) then
2947 write(be_out_unit) variable
2948 write(be_out_unit) be1_rf_lengthscale
2949 endif
2950 variable = be % v2 % name
2951 if (rootproc) then
2952 write(be_out_unit) variable
2953 write(be_out_unit) be2_rf_lengthscale
2954 endif
2955 variable = be % v3 % name
2956 if (rootproc) then
2957 write(be_out_unit) variable
2958 write(be_out_unit) be3_rf_lengthscale
2959 endif
2960 variable = be % v4 % name
2961 if (rootproc) then
2962 write(be_out_unit) variable
2963 write(be_out_unit) be4_rf_lengthscale
2964 endif
2965 variable = be % v5 % name
2966 if (rootproc) then
2967 write(be_out_unit) variable
2968 write(be_out_unit) be5_rf_lengthscale
2969 endif
2970 ! - Close the output unit:
2972 close(be_out_unit)
2973 ! 3.3.7 Deallocate the arrays for reading the input CV5 BE:
2975 DEALLOCATE ( be1_eval_loc )
2976 DEALLOCATE ( be2_eval_loc )
2977 DEALLOCATE ( be3_eval_loc )
2978 DEALLOCATE ( be4_eval_loc )
2979 DEALLOCATE ( be5_eval_loc )
2980 DEALLOCATE ( be1_evec_loc )
2981 DEALLOCATE ( be2_evec_loc )
2982 DEALLOCATE ( be3_evec_loc )
2983 DEALLOCATE ( be4_evec_loc )
2984 DEALLOCATE ( be5_evec_loc )
2986 ! Allocate the arrays for Model local BE:
2987 ALLOCATE ( be1_eval_loc (1:jy,1:kz) )
2988 ALLOCATE ( be2_eval_loc (1:jy,1:kz) )
2989 ALLOCATE ( be3_eval_loc (1:jy,1:kz) )
2990 ALLOCATE ( be4_eval_loc (1:jy,1:kz) )
2991 ALLOCATE ( be5_eval_loc (1:jy,1:1) )
2992 ALLOCATE ( be1_evec_loc(1:jy,1:kz,1:kz) )
2993 ALLOCATE ( be2_evec_loc(1:jy,1:kz,1:kz) )
2994 ALLOCATE ( be3_evec_loc(1:jy,1:kz,1:kz) )
2995 ALLOCATE ( be4_evec_loc(1:jy,1:kz,1:kz) )
2996 ALLOCATE ( be5_evec_loc(1:jy,1:1,1:1) )
2998 endif ! if interpolate_stats
3000 if ( use_rf ) then
3001 deallocate ( rfls1_be )
3002 deallocate ( rfls2_be )
3003 deallocate ( rfls3_be )
3004 deallocate ( rfls4_be )
3005 end if
3007 ! 4.0 Check and get the truncated number of the vertical modes
3008 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3010 if (vert_corr == vert_corr_2) then
3012 ! 4.1 Perform checks on eigenvectors:
3014 if (test_statistics) then
3015 call da_check_eof_decomposition(be1_eval_glo(:), be1_evec_glo(:,:), be % v1 % name)
3016 call da_check_eof_decomposition(be2_eval_glo(:), be2_evec_glo(:,:), be % v2 % name)
3017 call da_check_eof_decomposition(be3_eval_glo(:), be3_evec_glo(:,:), be % v3 % name)
3018 call da_check_eof_decomposition(be4_eval_glo(:), be4_evec_glo(:,:), be % v4 % name)
3020 if ( cloud_cv_options == 2 ) then
3021 call da_check_eof_decomposition(be6_eval_glo(:), be6_evec_glo(:,:), be % v6 % name)
3022 call da_check_eof_decomposition(be7_eval_glo(:), be7_evec_glo(:,:), be % v7 % name)
3023 call da_check_eof_decomposition(be8_eval_glo(:), be8_evec_glo(:,:), be % v8 % name)
3024 call da_check_eof_decomposition(be9_eval_glo(:), be9_evec_glo(:,:), be % v9 % name)
3025 call da_check_eof_decomposition(be10_eval_glo(:), be10_evec_glo(:,:), be % v10 % name)
3026 end if
3028 #if (WRF_CHEM == 1)
3029 if ( chem_cv_options >= 10 ) then
3030 do i=1,num_cv_3d_chem
3031 call da_check_eof_decomposition(be12_eval_glo(i,:), be12_evec_glo(i,:,:), be % v12(i) % name)
3032 end do
3033 end if
3034 #endif
3036 if ( use_cv_w ) then
3037 call da_check_eof_decomposition(be11_eval_glo(:), be11_evec_glo(:,:), be % v11 % name)
3038 end if
3040 end if
3042 ! 4.2 Truncate in vertical:
3044 call da_get_vertical_truncation(max_vert_var1, be1_eval_glo(:), be % v1)
3045 call da_get_vertical_truncation(max_vert_var2, be2_eval_glo(:), be % v2)
3046 call da_get_vertical_truncation(max_vert_var3, be3_eval_glo(:), be % v3)
3047 call da_get_vertical_truncation(max_vert_var4, be4_eval_glo(:), be % v4)
3049 if ( cloud_cv_options == 2 ) then
3050 call da_get_vertical_truncation(max_vert_var6, be6_eval_glo(:), be % v6)
3051 call da_get_vertical_truncation(max_vert_var7, be7_eval_glo(:), be % v7)
3052 call da_get_vertical_truncation(max_vert_var8, be8_eval_glo(:), be % v8)
3053 call da_get_vertical_truncation(max_vert_var9, be9_eval_glo(:), be % v9)
3054 call da_get_vertical_truncation(max_vert_var10,be10_eval_glo(:),be % v10)
3055 else
3056 if ( jb_factor > 0.0 ) then
3057 be % v6 % mz = xb % mkz
3058 be % v7 % mz = xb % mkz
3059 be % v8 % mz = xb % mkz
3060 be % v9 % mz = xb % mkz
3061 be % v10 % mz = xb % mkz
3062 else
3063 be % v6 % mz = 0
3064 be % v7 % mz = 0
3065 be % v8 % mz = 0
3066 be % v9 % mz = 0
3067 be % v10 % mz = 0
3068 end if
3069 end if
3071 #if (WRF_CHEM == 1)
3072 if ( chem_cv_options >= 10 ) then
3073 do i=1,num_cv_3d_chem
3074 call da_get_vertical_truncation(max_vert_var12(i), be12_eval_glo(i,:), be % v12(i))
3075 end do
3076 else
3077 if ( jb_factor > 0.0 ) then
3078 be % v12(:) % mz = xb % mkz
3079 else
3080 be % v12(:) % mz = 0
3081 end if
3082 end if
3083 #endif
3085 if ( use_cv_w ) then
3086 call da_get_vertical_truncation(max_vert_var11,be11_eval_glo(:),be % v11)
3087 else
3088 if ( jb_factor > 0.0 ) then
3089 be % v11 % mz = xb % mkz
3090 else
3091 be % v11 % mz = 0
3092 end if
3093 end if
3096 if (max_vert_var5 == 0.0) then
3097 be % v5 % mz = 0
3098 else
3099 be % v5 % mz = 1
3100 end if
3102 write (unit=stdout,fmt=*) ' '
3104 else
3106 ! 4.3 no truncated
3108 be % v1 % mz = xb % mkz
3109 be % v2 % mz = xb % mkz
3110 be % v3 % mz = xb % mkz
3111 be % v4 % mz = xb % mkz
3112 be % v5 % mz = xb % mkz
3114 be % v6 % mz = xb % mkz
3115 be % v7 % mz = xb % mkz
3116 be % v8 % mz = xb % mkz
3117 be % v9 % mz = xb % mkz
3118 be % v10% mz = xb % mkz
3119 be % v11% mz = xb % mkz
3121 #if (WRF_CHEM == 1)
3122 if ( chem_cv_options >= 10 ) then
3123 be % v12(:) % mz = xb % mkz
3124 end if
3125 #endif
3127 end if
3129 ! 5.0 Initialise control variable space components of header:
3130 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3132 ! 5.1 Compute the size of the control variables
3134 be % mix = ix
3135 be % mjy = jy
3137 ! 5.2 Transfer errors to error structure:
3139 call da_allocate_background_errors(jy, kz, be1_eval_glo, be1_evec_glo, be1_eval_loc, &
3140 be1_evec_loc, be % v1)
3141 call da_allocate_background_errors(jy, kz, be2_eval_glo, be2_evec_glo, be2_eval_loc, &
3142 be2_evec_loc, be % v2)
3143 call da_allocate_background_errors(jy, kz, be3_eval_glo, be3_evec_glo, be3_eval_loc, &
3144 be3_evec_loc, be % v3)
3145 call da_allocate_background_errors(jy, kz, be4_eval_glo, be4_evec_glo, be4_eval_loc, &
3146 be4_evec_loc, be % v4)
3148 ! 5.2.1 transfer the ps_u variance to be % v5:
3150 call da_allocate_background_errors(jy, 1, be5_eval_glo, be5_evec_glo, be5_eval_loc, &
3151 be5_evec_loc, be % v5)
3153 if ( cloud_cv_options >=2 ) then
3154 call da_allocate_background_errors(jy, kz, be6_eval_glo, be6_evec_glo, be6_eval_loc, &
3155 be6_evec_loc, be % v6)
3156 call da_allocate_background_errors(jy, kz, be7_eval_glo, be7_evec_glo, be7_eval_loc, &
3157 be7_evec_loc, be % v7)
3158 call da_allocate_background_errors(jy, kz, be8_eval_glo, be8_evec_glo, be8_eval_loc, &
3159 be8_evec_loc, be % v8)
3160 call da_allocate_background_errors(jy, kz, be9_eval_glo, be9_evec_glo, be9_eval_loc, &
3161 be9_evec_loc, be % v9)
3162 call da_allocate_background_errors(jy, kz, be10_eval_glo, be10_evec_glo, be10_eval_loc, &
3163 be10_evec_loc, be % v10)
3164 end if
3165 #if (WRF_CHEM == 1)
3166 if ( chem_cv_options >=10 ) then
3167 do i=1,num_cv_3d_chem
3168 call da_allocate_background_errors(jy, kz, be12_eval_glo(i,:), be12_evec_glo(i,:,:), be12_eval_loc(i,:,:), &
3169 be12_evec_loc(i,:,:,:), be % v12(i))
3170 end do
3171 end if
3172 #endif
3173 if ( use_cv_w ) then
3174 call da_allocate_background_errors(jy, kz, be11_eval_glo, be11_evec_glo, be11_eval_loc, &
3175 be11_evec_loc, be % v11)
3176 end if
3178 if (print_detail_be) then
3179 write (unit=stdout,fmt='(3x,a,i10)') "b) Finished eigenvector processing!"
3180 end if
3182 if( use_rf )then ! use recursive filters:
3184 ! 5.3 Convert the scale lengths in the real distance (meter)
3186 be1_rf_lengthscale(1:kz) = be1_rf_lengthscale(1:kz) * xb%ds
3187 be2_rf_lengthscale(1:kz) = be2_rf_lengthscale(1:kz) * xb%ds
3188 be3_rf_lengthscale(1:kz) = be3_rf_lengthscale(1:kz) * xb%ds
3189 be4_rf_lengthscale(1:kz) = be4_rf_lengthscale(1:kz) * xb%ds
3190 be5_rf_lengthscale(1:1) = be5_rf_lengthscale(1:1) * xb%ds
3192 if ( cloud_cv_options >= 2 ) then
3193 be6_rf_lengthscale(1:kz) = be6_rf_lengthscale(1:kz) * xb%ds
3194 be7_rf_lengthscale(1:kz) = be7_rf_lengthscale(1:kz) * xb%ds
3195 be8_rf_lengthscale(1:kz) = be8_rf_lengthscale(1:kz) * xb%ds
3196 be9_rf_lengthscale(1:kz) = be9_rf_lengthscale(1:kz) * xb%ds
3197 be10_rf_lengthscale(1:kz)= be10_rf_lengthscale(1:kz) * xb%ds
3198 end if
3199 #if (WRF_CHEM == 1)
3200 if ( chem_cv_options >= 10 ) then
3201 do i=1,num_cv_3d_chem
3202 be12_rf_lengthscale(i,1:kz) = be12_rf_lengthscale(i,1:kz) * xb%ds
3203 end do
3204 end if
3205 #endif
3206 if ( use_cv_w ) then
3207 be11_rf_lengthscale(1:kz)= be11_rf_lengthscale(1:kz) * xb%ds
3208 end if
3210 else ! use wavelets:
3211 ! kzs and mzs are used for wavelets
3212 kzs=(/kz,kz,kz,kz,1/)
3213 mzs=(/be%v1%mz,be%v2%mz,be%v3%mz,be%v4%mz,be%v5%mz/)
3215 read(be_unit)k ! trimmed string length
3216 if( k>len(variable) )then
3217 write(message(1),'(i0,">",i0," for use_rf=.false.")')k,len(variable)
3218 call da_error(__FILE__,__LINE__,message(1:1))
3219 endif
3220 read(be_unit)variable(1:k)
3221 if( variable(1:k) /= 'wavelet' )then
3222 write(message(1),'(A,"/=''wavelet''")')variable(1:k)
3223 call da_error(__FILE__,__LINE__,message(1:1))
3224 endif
3225 b = 0 ! pointer to variable*mode
3226 do m=1,5 ! variables loop:
3227 read(be_unit)k,nk
3228 if( k>len(variable) )then
3229 write(message(1),'(i0,">",i0," for ",A)')k,len(variable),namv(m)
3230 call da_error(__FILE__,__LINE__,message(1:1))
3231 elseif( nk /= kzs(m) )then
3232 write(message(1),'(i0,"=nk/=kzs(",i0,")=",i0," for ",A)')nk,m,kzs(m),namv(m)
3233 call da_error(__FILE__,__LINE__,message(1:1))
3234 endif
3235 read(be_unit)variable(1:k)
3236 if( variable(1:k) /= trim(namv(m)) )then
3237 write(message(1),'(A,"/=",A)')variable(1:k),trim(namv(m))
3238 call da_error(__FILE__,__LINE__,message(1:1))
3239 endif
3240 if( m==1 )then
3241 ! Possibly reassign namelist do_normalize value:
3242 read(be_unit)do_normalize,namw,lf,nb
3243 allocate(nij(0:nb,0:1,0:2))
3244 read(be_unit)nij ! wavelet indexes
3245 write(*,'(A,": ")',advance="no")__FILE__
3246 do i=0,nb ! wavelet-band loop:
3247 write(*,'(i2,"{",2(i3,","),i3,";",2(i3,","),i3,"}")',advance="no")&
3248 i,transpose(nij(i,:,:))
3249 enddo
3250 allocate(be%wsd(nij(0,1,2),nij(0,0,2),sum(mzs)), &
3251 !____________________max of {,i}dwtai wavelet scratch workspace sizes:
3252 ws(max(maxval(nij(0,:,2)), &
3253 2*floor(.5*(real(maxval(nij(0,:,0)))+lf))+lf-2)))
3254 if( do_normalize )allocate(be%sd(nij(0,1,0),nij(0,0,0),sum(mzs)))
3255 endif ! if( m==1 )
3256 do k=1,mzs(m)
3257 !___________mode-k field & wavelet-coefficient std. devs.:
3258 read(be_unit)be%wsd(:,:,b+k)
3259 if( do_normalize )read(be_unit)be%sd(:,:,b+k)
3260 enddo ! mode-k loop.
3261 do k=mzs(m)+1,kzs(m) ! read and discard truncated modes:
3262 read(be_unit)(junk,j=1,nij(0,1,2)*nij(0,0,2))
3263 if( do_normalize )read(be_unit)(junk,j=1,nij(0,1,0)*nij(0,0,0))
3264 enddo ! mode-k loop.
3265 write(*,'(A,": |wsd[",A,"]|=",es9.3)')__FILE__,namv(m),sqrt(sum(be%wsd(:,:,b+1:b+mzs(m))**2))
3266 if( do_normalize )write(*,'(A,": | sd[",A,"]|=",es9.3)')__FILE__,namv(m),sqrt(sum(be%sd(:,:,b+1:b+mzs(m))**2))
3267 b = b+mzs(m) ! point to next variable.
3268 enddo ! variables loop
3269 endif ! if( use_rf )
3271 ! 6.0 Perform checks on eigenvectors with be data structure:
3272 if (jb_factor > 0.0 .and. test_statistics) then
3273 call da_check_eof_decomposition(be%v1%val_g(:), be%v1%evec_g(:,:),&
3274 be%v1%name)
3275 call da_check_eof_decomposition(be%v2%val_g(:), be%v2%evec_g(:,:),&
3276 be%v2%name)
3277 call da_check_eof_decomposition(be%v3%val_g(:), be%v3%evec_g(:,:),&
3278 be%v3%name)
3279 call da_check_eof_decomposition(be%v4%val_g(:), be%v4%evec_g(:,:),&
3280 be%v4%name)
3282 if ( cloud_cv_options >= 2 ) then
3283 call da_check_eof_decomposition(be%v6%val_g(:), be%v6%evec_g(:,:),&
3284 be%v6%name)
3285 call da_check_eof_decomposition(be%v7%val_g(:), be%v7%evec_g(:,:),&
3286 be%v7%name)
3287 call da_check_eof_decomposition(be%v8%val_g(:), be%v8%evec_g(:,:),&
3288 be%v8%name)
3289 call da_check_eof_decomposition(be%v9%val_g(:), be%v9%evec_g(:,:),&
3290 be%v9%name)
3291 call da_check_eof_decomposition(be%v10%val_g(:), be%v10%evec_g(:,:),&
3292 be%v10%name)
3293 end if
3294 #if (WRF_CHEM == 1)
3295 if ( chem_cv_options >= 10 ) then
3296 do i=1,num_cv_3d_chem
3297 call da_check_eof_decomposition(be%v12(i)%val_g(:), be%v12(i)%evec_g(:,:),&
3298 be%v12(i)%name)
3299 end do
3300 end if
3301 #endif
3302 if ( use_cv_w ) then
3303 call da_check_eof_decomposition(be%v11%val_g(:), be%v11%evec_g(:,:),&
3304 be%v11%name)
3305 end if
3307 end if
3309 ! 6.1 Close the be unit
3311 close(be_unit)
3312 call da_free_unit(be_unit)
3314 if( use_rf )then
3315 ! 6.2 Keep the original be % v1, be % v2,...., and lengthscale in the first loop
3316 ! for the rescaling in the later loops:
3318 it = 1
3320 if (max_ext_its > 1 .and. jb_factor > 0.0) then
3322 if ( rootproc ) then
3323 write(unit=message(1),fmt='(A,I4)') '>>> Save the variances and scale-lengths in outer-loop', it
3324 call da_message(message(1:1))
3325 write(be_rf_unit) kz, jy, ix, be % v1 % mz, be % v2 % mz, be% v3 % mz, &
3326 be % v4 % mz, be % v5 % mz, xb % ds
3327 write(be_rf_unit) be % v1 % val, be % v2 % val, be% v3 % val, &
3328 be % v4 % val, be % v5 % val, &
3329 be1_rf_lengthscale, be2_rf_lengthscale, be3_rf_lengthscale, &
3330 be4_rf_lengthscale, be5_rf_lengthscale
3331 if ( cloud_cv_options >= 2 ) then
3332 write(be_rf_unit) be%v6%val, be%v7%val, be%v8%val, be%v9%val, be%v10%val
3333 write(be_rf_unit) be6_rf_lengthscale, be7_rf_lengthscale, be8_rf_lengthscale, &
3334 be9_rf_lengthscale, be10_rf_lengthscale
3335 end if
3336 #if (WRF_CHEM == 1)
3337 if ( chem_cv_options >= 10 ) then
3338 do i=1, num_cv_3d_chem
3339 write(be_rf_unit) be % v12(i) % val
3340 end do
3341 do i=1, num_cv_3d_chem
3342 write(be_rf_unit) be12_rf_lengthscale(i,:)
3343 end do
3344 end if
3345 #endif
3346 if ( use_cv_w ) then
3347 write(be_rf_unit) be%v11%val
3348 write(be_rf_unit) be11_rf_lengthscale
3349 end if
3350 end if ! rootproc
3352 if (print_detail_be ) then
3353 write(be_print_unit,'("it=",i2,2x,"kz=",i3,2x,"jy=",i4,2x,"ix=",i4,2x,"ds=",e12.5)') &
3354 it, kz, jy, ix, xb % ds
3355 write(be_print_unit,'("Original val and rf, and mz:",5i5)') &
3356 be % v1 % mz, be % v2 % mz, be% v3 % mz, be % v4 % mz, be % v5 % mz
3357 write(be_print_unit,'("mz=",i3,2x,"be%v1%val:"/(10e12.5))') be%v1%mz, be%v1%val(1,:)
3358 write(be_print_unit,'("mz=",i3,2x,"be%v2%val:"/(10e12.5))') be%v2%mz, be%v2%val(1,:)
3359 write(be_print_unit,'("mz=",i3,2x,"be%v3%val:"/(10e12.5))') be%v3%mz, be%v3%val(1,:)
3360 write(be_print_unit,'("mz=",i3,2x,"be%v4%val:"/(10e12.5))') be%v4%mz, be%v4%val(1,:)
3361 write(be_print_unit,'("mz=",i3,2x,"be%v5%val:"/(10e12.5))') be%v5%mz, be%v5%val(1,:)
3362 write(be_print_unit,'(/"scale-length: kz=",i3)') kz
3363 do i = 1,kz
3364 if (i == 1) then
3365 write(be_print_unit,'(i3,2x,5e15.5)') i,be1_rf_lengthscale(i), &
3366 be2_rf_lengthscale(i), be3_rf_lengthscale(i), be4_rf_lengthscale(i), &
3367 be5_rf_lengthscale(i)
3368 else
3369 write(be_print_unit,'(i3,2x,5e15.5)') i,be1_rf_lengthscale(i), &
3370 be2_rf_lengthscale(i), be3_rf_lengthscale(i), be4_rf_lengthscale(i)
3371 endif
3372 enddo
3374 endif
3376 endif
3378 ! 7.0 Apply empirical and recursive filter rescaling factor:
3379 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3381 call da_rescale_background_errors(var_scaling1(1), len_scaling1(1), &
3382 xb % ds, be1_rf_lengthscale, be % v1)
3383 call da_rescale_background_errors(var_scaling2(1), len_scaling2(1), &
3384 xb % ds, be2_rf_lengthscale, be % v2)
3385 call da_rescale_background_errors(var_scaling3(1), len_scaling3(1), &
3386 xb % ds, be3_rf_lengthscale, be % v3)
3387 call da_rescale_background_errors(var_scaling4(1), len_scaling4(1), &
3388 xb % ds, be4_rf_lengthscale, be % v4)
3389 call da_rescale_background_errors(var_scaling5(1), len_scaling5(1), &
3390 xb % ds, be5_rf_lengthscale, be % v5)
3391 if ( cloud_cv_options >= 2 ) then
3392 call da_rescale_background_errors(var_scaling6(1), len_scaling6(1), &
3393 xb % ds, be6_rf_lengthscale, be % v6)
3394 call da_rescale_background_errors(var_scaling7(1), len_scaling7(1), &
3395 xb % ds, be7_rf_lengthscale, be % v7)
3396 call da_rescale_background_errors(var_scaling8(1), len_scaling8(1), &
3397 xb % ds, be8_rf_lengthscale, be % v8)
3398 call da_rescale_background_errors(var_scaling9(1), len_scaling9(1), &
3399 xb % ds, be9_rf_lengthscale, be % v9)
3400 call da_rescale_background_errors(var_scaling10(1), len_scaling10(1), &
3401 xb % ds, be10_rf_lengthscale, be % v10)
3402 end if
3403 #if (WRF_CHEM == 1)
3404 if ( chem_cv_options >= 10 ) then
3405 do i=1,num_cv_3d_chem
3406 call da_rescale_background_errors(var_scaling12(i), len_scaling12(i), &
3407 xb % ds, be12_rf_lengthscale(i,:), be % v12(i))
3408 end do
3409 end if
3410 #endif
3411 if ( use_cv_w ) then
3412 call da_rescale_background_errors(var_scaling11(1), len_scaling11(1), &
3413 xb % ds, be11_rf_lengthscale, be % v11)
3414 end if
3416 end if ! use_rf
3418 if (print_detail_be .and. jb_factor > 0.0) then
3419 write(be_print_unit,'(/"============================================================")')
3420 write(be_print_unit,'("For outer loop ",i2)') it
3421 write(be_print_unit,'("it=",i2,2x,"kz=",i3,2x,"jy=",i4,2x,"ix=",i4,2x,"ds=",e12.5)') &
3422 it, kz, jy, ix, xb % ds
3423 write(be_print_unit,'("Namelist options specified for this iteration:")')
3424 write(be_print_unit,'("var_scaling1(it) = ",e12.5,2x,"len_scaling1(it) = "e12.5)')var_scaling1(it),len_scaling1(it)
3425 write(be_print_unit,'("var_scaling2(it) = ",e12.5,2x,"len_scaling2(it) = "e12.5)')var_scaling2(it),len_scaling2(it)
3426 write(be_print_unit,'("var_scaling3(it) = ",e12.5,2x,"len_scaling3(it) = "e12.5)')var_scaling3(it),len_scaling3(it)
3427 write(be_print_unit,'("var_scaling4(it) = ",e12.5,2x,"len_scaling4(it) = "e12.5)')var_scaling4(it),len_scaling4(it)
3428 write(be_print_unit,'("var_scaling5(it) = ",e12.5,2x,"len_scaling5(it) = "e12.5)')var_scaling5(it),len_scaling5(it)
3429 write(be_print_unit,'("Background error statistics for this iteration:")')
3430 write(be_print_unit,'("mz=",i3,2x,"be%v1%val:"/(10e12.5))') be%v1%mz, be%v1%val(1,:)
3431 write(be_print_unit,'("mz=",i3,2x,"be%v2%val:"/(10e12.5))') be%v2%mz, be%v2%val(1,:)
3432 write(be_print_unit,'("mz=",i3,2x,"be%v3%val:"/(10e12.5))') be%v3%mz, be%v3%val(1,:)
3433 write(be_print_unit,'("mz=",i3,2x,"be%v4%val:"/(10e12.5))') be%v4%mz, be%v4%val(1,:)
3434 write(be_print_unit,'("mz=",i3,2x,"be%v5%val:"/(10e12.5))') be%v5%mz, be%v5%val(1,:)
3436 if ( cloud_cv_options >= 2 ) then
3437 write(be_print_unit,'("mz=",i3,2x,"be%v6%val:"/(10e12.5))') be%v6%mz, be%v6%val(1,:)
3438 write(be_print_unit,'("mz=",i3,2x,"be%v7%val:"/(10e12.5))') be%v7%mz, be%v7%val(1,:)
3439 write(be_print_unit,'("mz=",i3,2x,"be%v8%val:"/(10e12.5))') be%v8%mz, be%v8%val(1,:)
3440 write(be_print_unit,'("mz=",i3,2x,"be%v9%val:"/(10e12.5))') be%v9%mz, be%v9%val(1,:)
3441 write(be_print_unit,'("mz=",i3,2x,"be%v10%val:"/(10e12.5))') be%v10%mz, be%v10%val(1,:)
3442 end if
3443 #if (WRF_CHEM == 1)
3444 if ( chem_cv_options >= 10 ) then
3445 do i=1,num_cv_3d_chem
3446 write(be_print_unit,'("mz=",i3,2x,"be%v12%val:"/(10e12.5))') be%v12(i)%mz, be%v12(i)%val(1,:)
3447 end do
3448 end if
3449 #endif
3450 if ( use_cv_w ) then
3451 write(be_print_unit,'("mz=",i3,2x,"be%v11%val:"/(10e12.5))') be%v11%mz, be%v11%val(1,:)
3452 end if
3454 write(be_print_unit,'(/"scale-length: kz=",i3)') kz
3455 write(be_print_unit,'("be%v1%rf_alpha:"/(10e12.5))') be % v1 % rf_alpha(:)
3456 write(be_print_unit,'("be%v2%rf_alpha:"/(10e12.5))') be % v2 % rf_alpha(:)
3457 write(be_print_unit,'("be%v3%rf_alpha:"/(10e12.5))') be % v3 % rf_alpha(:)
3458 write(be_print_unit,'("be%v4%rf_alpha:"/(10e12.5))') be % v4 % rf_alpha(:)
3459 write(be_print_unit,'("be%v5%rf_alpha:"/(10e12.5))') be % v5 % rf_alpha(:)
3461 if ( cloud_cv_options >= 2 ) then
3462 write(be_print_unit,'("be%v6%rf_alpha:"/(10e12.5))') be % v6 % rf_alpha(:)
3463 write(be_print_unit,'("be%v7%rf_alpha:"/(10e12.5))') be % v7 % rf_alpha(:)
3464 write(be_print_unit,'("be%v8%rf_alpha:"/(10e12.5))') be % v8 % rf_alpha(:)
3465 write(be_print_unit,'("be%v9%rf_alpha:"/(10e12.5))') be % v9 % rf_alpha(:)
3466 write(be_print_unit,'("be%v10%rf_alpha:"/(10e12.5))') be % v10 % rf_alpha(:)
3467 end if
3468 #if (WRF_CHEM == 1)
3469 if ( chem_cv_options >= 10 ) then
3470 do i=1,num_cv_3d_chem
3471 write(be_print_unit,'("be%v12%rf_alpha:"/(10e12.5))') be % v12(i) % rf_alpha(:)
3472 end do
3473 end if
3474 #endif
3475 if ( use_cv_w ) then
3476 write(be_print_unit,'("be%v11%rf_alpha:"/(10e12.5))') be % v11 % rf_alpha(:)
3477 end if
3479 write(be_print_unit,'(/"scale-length: kz=",i3)') kz
3480 do i = 1,kz
3481 if (i == 1) then
3482 write(be_print_unit,'(i3,2x,5e15.5)') i, be1_rf_lengthscale(i), be2_rf_lengthscale(i), &
3483 be3_rf_lengthscale(i), be4_rf_lengthscale(i), be5_rf_lengthscale(i)
3484 else
3485 write(be_print_unit,'(i3,2x,4e15.5)') i, be1_rf_lengthscale(i), be2_rf_lengthscale(i), &
3486 be3_rf_lengthscale(i), be4_rf_lengthscale(i)
3487 end if
3488 enddo
3489 if ( cloud_cv_options >= 2 ) then
3490 do i = 1,kz
3491 write(be_print_unit,'(i3,2x,5e15.5)') i, be6_rf_lengthscale(i), be7_rf_lengthscale(i), &
3492 be8_rf_lengthscale(i), be9_rf_lengthscale(i), be10_rf_lengthscale(i)
3493 enddo
3494 end if
3495 !#if (WRF_CHEM == 1)
3496 ! if ( chem_cv_options >= 10 ) then
3497 ! do i = 1,kz
3498 ! write(be_print_unit,'(i3,2x,22e15.5)') i, be12_rf_lengthscale(ic,i)
3499 ! enddo
3500 ! end if
3501 !#endif
3503 endif ! print_detail_be
3505 close(be_rf_unit)
3506 close(be_print_unit)
3507 call da_free_unit(be_rf_unit)
3508 call da_free_unit(be_print_unit)
3510 ! 8.0 deallocate input model state:
3511 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3513 deallocate (be1_eval_loc)
3514 deallocate (be2_eval_loc)
3515 deallocate (be3_eval_loc)
3516 deallocate (be4_eval_loc)
3517 deallocate (be5_eval_loc)
3519 if ( cloud_cv_options >= 2 ) then
3520 deallocate (be6_eval_loc)
3521 deallocate (be7_eval_loc)
3522 deallocate (be8_eval_loc)
3523 deallocate (be9_eval_loc)
3524 deallocate (be10_eval_loc)
3525 end if
3526 #if (WRF_CHEM == 1)
3527 if ( chem_cv_options >= 10 ) then
3528 deallocate (be12_eval_loc)
3529 end if
3530 #endif
3531 if ( use_cv_w ) then
3532 deallocate (be11_eval_loc)
3533 end if
3535 if( use_rf )then
3536 deallocate (be1_rf_lengthscale)
3537 deallocate (be2_rf_lengthscale)
3538 deallocate (be3_rf_lengthscale)
3539 deallocate (be4_rf_lengthscale)
3540 deallocate (be5_rf_lengthscale)
3542 if ( cloud_cv_options >= 2 ) then
3543 deallocate (be6_rf_lengthscale)
3544 deallocate (be7_rf_lengthscale)
3545 deallocate (be8_rf_lengthscale)
3546 deallocate (be9_rf_lengthscale)
3547 deallocate (be10_rf_lengthscale)
3548 end if
3549 #if (WRF_CHEM == 1)
3550 if ( chem_cv_options >= 10 ) then
3551 deallocate (be12_rf_lengthscale)
3552 end if
3553 #endif
3554 if ( use_cv_w ) then
3555 deallocate (be11_rf_lengthscale)
3556 end if
3557 endif
3559 if (vert_corr == vert_corr_2) then
3560 deallocate (be1_eval_glo)
3561 deallocate (be2_eval_glo)
3562 deallocate (be3_eval_glo)
3563 deallocate (be4_eval_glo)
3564 deallocate (be5_eval_glo)
3566 if ( cloud_cv_options >= 2 ) then
3567 deallocate (be6_eval_glo)
3568 deallocate (be7_eval_glo)
3569 deallocate (be8_eval_glo)
3570 deallocate (be9_eval_glo)
3571 deallocate (be10_eval_glo)
3572 end if
3573 #if (WRF_CHEM == 1)
3574 if ( chem_cv_options >= 10 ) then
3575 deallocate (be12_eval_glo)
3576 end if
3577 #endif
3578 if ( use_cv_w ) then
3579 deallocate (be11_eval_glo)
3580 end if
3582 deallocate (be1_evec_loc)
3583 deallocate (be2_evec_loc)
3584 deallocate (be3_evec_loc)
3585 deallocate (be4_evec_loc)
3586 deallocate (be5_evec_loc)
3588 if ( cloud_cv_options >= 2 ) then
3589 deallocate (be6_evec_loc)
3590 deallocate (be7_evec_loc)
3591 deallocate (be8_evec_loc)
3592 deallocate (be9_evec_loc)
3593 deallocate (be10_evec_loc)
3594 end if
3595 #if (WRF_CHEM == 1)
3596 if ( chem_cv_options >= 10 ) then
3597 deallocate (be12_evec_loc)
3598 end if
3599 #endif
3600 if ( use_cv_w ) then
3601 deallocate (be11_evec_loc)
3602 end if
3604 deallocate (be1_evec_glo)
3605 deallocate (be2_evec_glo)
3606 deallocate (be3_evec_glo)
3607 deallocate (be4_evec_glo)
3608 deallocate (be5_evec_glo)
3610 if ( cloud_cv_options >= 2 ) then
3611 deallocate (be6_evec_glo)
3612 deallocate (be7_evec_glo)
3613 deallocate (be8_evec_glo)
3614 deallocate (be9_evec_glo)
3615 deallocate (be10_evec_glo)
3616 end if
3617 #if (WRF_CHEM == 1)
3618 if ( chem_cv_options >= 10 ) then
3619 deallocate (be12_evec_glo)
3620 end if
3621 #endif
3622 if ( use_cv_w ) then
3623 deallocate (be11_evec_glo)
3624 end if
3626 end if
3628 deallocate (bin)
3629 deallocate (bin2d)
3631 else ! jb_factor <= 0.0
3632 ni = ix
3633 nj = jy
3634 nk = kz
3635 be % v1 % mz = 0
3636 be % v2 % mz = 0
3637 be % v3 % mz = 0
3638 be % v4 % mz = 0
3639 be % v5 % mz = 0
3640 be % v6 % mz = 0
3641 be % v7 % mz = 0
3642 be % v8 % mz = 0
3643 be % v9 % mz = 0
3644 be % v10 % mz = 0
3645 be % v11 % mz = 0
3646 end if ! jb_factor > 0.0
3648 if (be % ne > 0) then
3649 be % alpha % name = 'alpha'
3650 allocate (alpha_val(1:nk)) ! Not using jy dimension yet.
3651 allocate (alpha_evec(1:nk,1:nk)) ! Not using jy dimension yet.
3653 if ( alpha_vertloc_opt > 0 ) then ! Use vertical localization:
3654 if ( rootproc ) then
3655 call da_get_alpha_vertloc(xb, alpha_val, alpha_evec)
3656 end if
3657 call wrf_dm_bcast_real(alpha_val, nk)
3658 call wrf_dm_bcast_real(alpha_evec, nk*nk)
3659 be % alpha % mz = nk
3660 call da_get_vertical_truncation(max_vert_var_alpha, alpha_val, be % alpha)
3661 else
3662 be % alpha % mz = 1 ! No vertical localization.
3663 alpha_val(1) = 1.0
3664 alpha_val(2:kz) = 0.0
3665 alpha_evec(:,:) = 1.0
3666 end if
3667 mz = be % alpha % mz
3669 ! Alpha eigenvalues and eigenvectors:
3670 allocate (be % alpha % val(1:jy,1:mz)) ! Not using jy dimension but here for consistency.
3671 allocate (be % alpha % evec(1:nj,1:nk,1:mz))
3673 if ( anal_type_hybrid_dual_res ) then
3674 deallocate(be % alpha % val, be % alpha % evec)
3675 allocate (be % alpha % val(jds_int:jde_int,1:mz))
3676 allocate (be % alpha % evec(jds_int:jde_int,1:nk,1:mz))
3677 endif
3680 do m = 1, mz
3681 be % alpha % val(:,m) = sigma_alpha * alpha_val(m)
3682 do k = 1, nk
3683 be % alpha % evec(:,k,m) = alpha_evec(k,m)
3684 end do
3685 end do
3687 ! Alpha RF lengthscales and variance scaling factors:
3688 allocate (alpha_rf_lengthscale(1:mz))
3689 allocate (be % alpha % rf_alpha(1:mz))
3690 allocate (alpha_rf_scale_factor(1:mz))
3693 if ( anal_type_hybrid_dual_res ) then
3694 alpha_rf_lengthscale(1:mz) = 1000.0 * alpha_corr_scale / (grid%intermediate_grid % dx )
3695 else
3696 alpha_rf_lengthscale(1:mz) = 1000.0 * alpha_corr_scale / xb % ds ! Convert km to grid spacings.
3697 endif
3700 call da_calculate_rf_factors( alpha_rf_lengthscale(:), be % alpha % rf_alpha(:), &
3701 alpha_rf_scale_factor(:) )
3702 do m = 1, mz
3703 be % alpha % val(:,m) = alpha_rf_scale_factor(m) * be % alpha % val(:,m)
3704 end do
3706 if( .not. use_rf ) then
3707 allocate(be%alpha%wsd(nij(0,1,2),nij(0,0,2),mz))
3708 call random_number(be%alpha%wsd) ! need to parallelize
3709 if( do_normalize )then
3710 allocate(be%alpha%sd(nij(0,1,0),nij(0,0,0),mz))
3711 call random_number(be%alpha%sd) ! need to parallelize
3712 endif
3713 endif
3715 deallocate(alpha_val)
3716 deallocate(alpha_evec)
3717 deallocate(alpha_rf_lengthscale)
3718 deallocate(alpha_rf_scale_factor)
3720 else
3721 be % alpha % mz = 0
3722 end if
3724 be%cv_mz(1) = be%v1%mz
3725 be%cv_mz(2) = be%v2%mz
3726 be%cv_mz(3) = be%v3%mz
3727 be%cv_mz(4) = be%v4%mz
3728 be%cv_mz(5) = be%v5%mz
3729 if ( cloud_cv_options >= 2 ) then
3730 be%cv_mz(6) = be%v6%mz
3731 be%cv_mz(7) = be%v7%mz
3732 be%cv_mz(8) = be%v8%mz
3733 be%cv_mz(9) = be%v9%mz
3734 be%cv_mz(10) = be%v10%mz
3735 end if
3736 #if (WRF_CHEM == 1)
3737 if ( chem_cv_options >= 10 ) then
3738 be % ncv_mz_chemic = num_cv_3d_chem
3739 allocate ( be%cv_mz_chemic(1:be%ncv_mz_chemic) )
3740 be%cv_mz_chemic(:) = 0 ! initialize
3742 do i=1,num_cv_3d_chem
3743 be%cv_mz_chemic(i) = be%v12(i)%mz
3744 end do
3745 end if
3746 #endif
3747 if ( use_cv_w ) then
3748 if (be%ncv_mz ==8 .or. be%ncv_mz ==13 ) then
3749 be%cv_mz(be%ncv_mz-2) = be%v11%mz
3750 else
3751 be%cv_mz(num_cv_3d_basic+num_cv_3d_extra+num_cv_2d) = be%v11%mz
3752 end if
3753 !be%cv_mz(num_cv_3d_basic+num_cv_3d_extra+num_cv_2d) = be%v11%mz ! This is a bug. (Soyoung on Aug-2021)
3754 end if
3755 be%cv_mz(be%ncv_mz-1) = be%alpha%mz
3756 be%cv_mz(be%ncv_mz) = be%ne
3758 if (trace_use) call da_trace_exit("da_setup_be_regional")
3760 end subroutine da_setup_be_regional