da_setup_be_regional.inc
References to this file elsewhere.
1 subroutine da_setup_be_regional(xb, be)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Define and allocate components of background errors
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (xb_type), intent(in) :: xb ! First guess structure.
10 type (be_type), intent(out) :: be ! Back. errors structure.
11
12 integer :: i, j, k, m ! Loop counters.
13 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point
14 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point
15 integer :: bin_type ! Type of bin to average over.
16 integer :: num_bins ! Number of bins (3D fields).
17 integer :: num_bins2d ! Number of bins (3D fields).
18 real :: lat_min, lat_max, binwidth_lat ! Used if bin_type = 2 (degrees)..
19 real :: hgt_min, hgt_max, binwidth_hgt ! Used if bin_type = 2 (m). .
20
21 real, allocatable :: be1_eval_loc(:,:) ! Temp arrays.
22 real, allocatable :: be2_eval_loc(:,:) ! Temp arrays.
23 real, allocatable :: be3_eval_loc(:,:) ! Temp arrays.
24 real, allocatable :: be4_eval_loc(:,:) ! Temp arrays.
25 real, allocatable :: be5_eval_loc(:,:) ! Temp arrays.
26
27 real, allocatable :: be1_eval_glo(:) ! Global Eigenvalues.
28 real, allocatable :: be2_eval_glo(:) ! Global Eigenvalues.
29 real, allocatable :: be3_eval_glo(:) ! Global Eigenvalues.
30 real, allocatable :: be4_eval_glo(:) ! Global Eigenvalues.
31 real, allocatable :: be5_eval_glo(:) ! Global Eigenvalues.
32
33 real, allocatable :: be1_evec_loc(:,:,:) ! Local Eigenvectors.
34 real, allocatable :: be2_evec_loc(:,:,:) ! Local Eigenvectors.
35 real, allocatable :: be3_evec_loc(:,:,:) ! Local Eigenvectors.
36 real, allocatable :: be4_evec_loc(:,:,:) ! Local Eigenvectors.
37 real, allocatable :: be5_evec_loc(:,:,:) ! Local Eigenvectors.
38
39 real, allocatable :: be1_evec_glo(:,:) ! Global Eigenvectors.
40 real, allocatable :: be2_evec_glo(:,:) ! Global Eigenvectors.
41 real, allocatable :: be3_evec_glo(:,:) ! Global Eigenvectors.
42 real, allocatable :: be4_evec_glo(:,:) ! Global Eigenvectors.
43 real, allocatable :: be5_evec_glo(:,:) ! Global Eigenvectors.
44
45 real, allocatable :: be1_rf_lengthscale(:) ! RF lengthscale.
46 real, allocatable :: be2_rf_lengthscale(:) ! RF lengthscale.
47 real, allocatable :: be3_rf_lengthscale(:) ! RF lengthscale.
48 real, allocatable :: be4_rf_lengthscale(:) ! RF lengthscale.
49 real, allocatable :: be5_rf_lengthscale(:)
50 real, allocatable :: alpha_rf_lengthscale(:)
51
52 real, allocatable :: evec_loc(:,:,:) ! Latitudinally varying eigenvectors.
53 real, allocatable :: eval_loc(:,:) ! Latitudinally varying eigenvalues.
54
55 character*10 :: variable
56 integer :: ni, nj, nk, nk_2d, b
57 integer :: ix, jy, kz
58 real, allocatable :: regcoeff1(:)
59 real, allocatable :: regcoeff2(:,:)
60 real, allocatable :: regcoeff3(:,:,:)
61 real :: avg,avg2,avg3
62 integer :: be_unit
63
64 if (trace_use) call da_trace_entry("da_setup_be_regional")
65
66 write (unit=message(1),fmt='(A)') &
67 'Set up background errors for regional application'
68
69 be % v1 % name = 'psi ' ! Streamfunction
70 be % v2 % name = 'chi_u' ! Uncorrelated velocity potential.
71 be % v3 % name = 't_u' ! Unbalanced temperature.
72 be % v4 % name = 'q/qsg'
73 be % v5 % name = 'psfc' ! surface pressure
74 !clq
75 if(use_radarobs .and. use_radar_rf) then
76 be % v4 % name = 'qt '
77 end if
78 !clq end
79
80 write (unit=message(2),fmt='(3x,A,7A)') 'WRF-Var dry control variables are:', &
81 trim(be % v1 % name), ', ', trim(be % v2 % name), ', ', &
82 trim(be % v3 % name), ' and ', trim(be % v5 % name)
83
84 write (unit=message(3),fmt='(3x,A,A)') &
85 'Humidity control variable is ', trim(be % v4 % name)
86
87 ix = xb % mix
88 jy = xb % mjy
89 kz = xb % mkz
90
91 call da_get_unit(be_unit)
92 open(unit=be_unit,file="be.dat", status="old",form="unformatted")
93
94 rewind (be_unit)
95 read (be_unit) ni, nj, nk
96
97 allocate (bin(1:ni,1:nj,1:nk))
98 allocate (bin2d(1:ni,1:nj))
99
100 read (be_unit)bin_type
101 read (be_unit)lat_min, lat_max, binwidth_lat
102 read (be_unit)hgt_min, hgt_max, binwidth_hgt
103 read (be_unit)num_bins, num_bins2d
104 read (be_unit)bin(1:ni,1:nj,1:nk)
105 read (be_unit)bin2d(1:ni,1:nj)
106
107 ! 1.1 Read in regression coefficients
108
109 allocate (regcoeff1(1:num_bins))
110 allocate (regcoeff2(1:nk,1:num_bins2d))
111 allocate (regcoeff3(1:nk,1:nk,1:num_bins2d))
112
113 read (be_unit) regcoeff1
114 read (be_unit) regcoeff2
115 read (be_unit) regcoeff3
116
117 ! 1.2 Fill regression coeff. array
118
119 allocate (be%reg_chi(1:nj,1:nk))
120 allocate (be%reg_ps (1:nj,1:nk))
121 allocate (be%reg_t (1:nj,1:nk,1:nk))
122
123 do k=1,nk
124 do j =1, nj
125 b = bin(1,j,k)
126 be%reg_chi(j,k) = regcoeff1(b)
127 end do
128 end do
129
130 do j=1,nj
131 b = bin2d(1,j)
132 do k=1,nk
133 be%reg_ps(j,k) = regcoeff2(k,b)
134 end do
135 end do
136
137 do j=1,nj
138 b = bin2d(1,j)
139 do i=1,nk
140 do k=1,nk
141 be%reg_t(j,i,k) = regcoeff3(i,k,b)
142 end do
143 end do
144 end do
145
146 deallocate (regcoeff1)
147 deallocate (regcoeff2)
148 deallocate (regcoeff3)
149
150 ! 1.3 Domain_averaged regression coefficients
151
152 if (.not.lat_stats_option) then
153 write (unit=message(4), fmt='(3x, a)') &
154 'Using the averaged regression coefficients for unbalanced part'
155
156 do k=1,nk
157 avg= 0.0
158 avg2=0.0
159 do j=1,num_bins2d
160 avg= avg + be%reg_ps (j,k)/float(num_bins2d)
161 avg2= avg2 + be%reg_chi (j,k)/float(num_bins2d)
162 end do
163
164 do j=1,num_bins2d
165 be%reg_ps(j,k)=avg
166 be%reg_chi (j,k)=avg2
167 end do
168 end do
169
170 do m=1,nk
171 do k=1,nk
172 avg3= 0.0
173
174 do j=1,num_bins2d
175 avg3= avg3 + be%reg_t (j,k,m)/float(num_bins2d)
176 end do
177
178 do j=1,num_bins2d
179 be%reg_t(j,k,m)=avg3
180 end do
181 end do
182 end do
183 else
184 write (unit=message(4), fmt='(3x, a)') &
185 'Using the latitude-dependent regression coefficients for unbalanced part'
186 end if
187
188 call da_message(message(1:4))
189
190 ! 2.0 Load the eigenvector and eigenvalue
191
192 allocate (be1_eval_loc (1:jy,1:kz))
193 allocate (be2_eval_loc (1:jy,1:kz))
194 allocate (be3_eval_loc (1:jy,1:kz))
195 allocate (be4_eval_loc (1:jy,1:kz))
196 allocate (be5_eval_loc (1:jy,1:1))
197
198 if (vert_corr == vert_corr_2) then
199
200 allocate (be1_eval_glo(1:kz))
201 allocate (be2_eval_glo(1:kz))
202 allocate (be3_eval_glo(1:kz))
203 allocate (be4_eval_glo(1:kz))
204 allocate (be5_eval_glo(1:1))
205
206 allocate (be1_evec_loc(1:jy,1:kz,1:kz))
207 allocate (be2_evec_loc(1:jy,1:kz,1:kz))
208 allocate (be3_evec_loc(1:jy,1:kz,1:kz))
209 allocate (be4_evec_loc(1:jy,1:kz,1:kz))
210 allocate (be5_evec_loc(1:jy,1: 1,1: 1))
211
212 allocate (be1_evec_glo(1:kz,1:kz))
213 allocate (be2_evec_glo(1:kz,1:kz))
214 allocate (be3_evec_glo(1:kz,1:kz))
215 allocate (be4_evec_glo(1:kz,1:kz))
216 allocate (be5_evec_glo(1:1,1:1))
217 end if
218
219 ! 2.2 Read in the eigenvector and eigenvalue
220
221 ! 2.2.1 Control variable 1 (psi)
222
223 read (be_unit) variable
224 read (be_unit) nk, num_bins2d
225 read (be_unit) be1_evec_glo
226 read (be_unit) be1_eval_glo
227 allocate (evec_loc(1:nk,1:nk,1:num_bins2d))
228 allocate (eval_loc(1:nk, 1:num_bins2d))
229 read (be_unit) evec_loc
230 read (be_unit) eval_loc
231 if (variable(1:3) /= 'psi') then
232 message(1)=' Variable mismatch while transfering eigen values from BE file '
233 write (unit=message(2),fmt='(A,A)') ' Expected psi but got ',variable
234 call da_error(__FILE__,__LINE__,message(1:2))
235 end if
236 be % v1 % name = variable
237 do j=1,jy
238 b = bin2d(1,j)
239 be1_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
240 be1_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
241 end do
242
243 ! 2.2.2 Control variable 2 (chi_u)
244
245 read (be_unit) variable
246 read (be_unit) nk, num_bins2d
247 read (be_unit) be2_evec_glo
248 read (be_unit) be2_eval_glo
249 read (be_unit) evec_loc
250 read (be_unit) eval_loc
251 if (variable(1:5) /= 'chi_u') then
252 message(1)=' Variable mismatch while transfering eigen values from BE file '
253 write (unit=message(2),fmt='(A,A)') ' Expected chi_u but got ',variable
254 call da_error(__FILE__,__LINE__,message(1:2))
255 end if
256 be % v2 % name = variable
257 do j=1,jy
258 b = bin2d(1,j)
259 be2_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
260 be2_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
261 end do
262
263 ! 2.2.3 Control variable 3 (t_u)
264 read (be_unit) variable
265 read (be_unit) nk, num_bins2d
266 read (be_unit) be3_evec_glo
267 read (be_unit) be3_eval_glo
268 read (be_unit) evec_loc
269 read (be_unit) eval_loc
270 if (variable(1:3) /= 't_u') then
271 message(1)=' Variable mismatch while transfering eigen values from BE file '
272 write (unit=message(2),fmt='(A,A)') ' Expected t_u but got ',variable
273 call da_error(__FILE__,__LINE__,message(1:2))
274 end if
275 be % v3 % name = variable
276 do j=1,jy
277 b = bin2d(1,j)
278 be3_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
279 be3_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
280 end do
281
282 ! 2.2.4 Control variable 4 (q/qsg)
283
284 read (be_unit) variable
285 read (be_unit) nk, num_bins2d
286 read (be_unit) be4_evec_glo
287 read (be_unit) be4_eval_glo
288 read (be_unit) evec_loc
289 read (be_unit) eval_loc
290 if (variable(1:3) /= 'rh') then
291 message(1)=' Variable mismatch while transfering eigen values from BE file '
292 write (unit=message(2),fmt='(A,A)') ' Expected rh but got ',variable
293 call da_error(__FILE__,__LINE__,message(1:2))
294 end if
295 be % v4 % name = variable
296 do j=1,jy
297 b = bin2d(1,j)
298 be4_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b)
299 be4_eval_loc(j,1:nk ) = eval_loc(1:nk,b)
300 end do
301
302 deallocate (evec_loc)
303 deallocate (eval_loc)
304
305 ! 2.2.5 Control variable ps_u
306
307 read (be_unit) variable
308 read (be_unit) nk_2d, num_bins2d
309 allocate (evec_loc(1:nk_2d,1:nk_2d,1:num_bins2d))
310 allocate (eval_loc(1:nk_2d, 1:num_bins2d))
311 read (be_unit) be5_evec_glo
312 read (be_unit) be5_eval_glo
313 read (be_unit) evec_loc
314 read (be_unit) eval_loc
315 if (variable(1:4) /= 'ps_u') then
316 message(1)=' Variable mismatch while transfering eigen values from BE file '
317 write (unit=message(2),fmt='(A,A)') ' Expected ps_u but got ',variable
318 call da_error(__FILE__,__LINE__,message(1:2))
319 end if
320 be % v5 % name = variable
321 do j=1,jy
322 b = bin2d(1,j)
323 be5_evec_loc(j,1:1,1:1) = evec_loc(1:1,1:1,b)
324 be5_eval_loc(j,1:1 ) = eval_loc(1:1,b)
325 end do
326
327 deallocate (evec_loc)
328 deallocate (eval_loc)
329
330 ! 3.0 Check and get the truncated number of the vertical modes
331 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
332
333 if (vert_corr == vert_corr_2) then
334
335 ! 3.1 Perform checks on eigenvectors:
336
337 if (test_statistics) then
338 call da_check_eof_decomposition(be1_eval_glo(:), be1_evec_glo(:,:), be % v1 % name)
339 call da_check_eof_decomposition(be2_eval_glo(:), be2_evec_glo(:,:), be % v2 % name)
340 call da_check_eof_decomposition(be3_eval_glo(:), be3_evec_glo(:,:), be % v3 % name)
341 call da_check_eof_decomposition(be4_eval_glo(:), be4_evec_glo(:,:), be % v4 % name)
342 end if
343
344 ! 3.2 Truncate in vertical:
345
346 call da_get_vertical_truncation(max_vert_var1, be1_eval_glo(:), be % v1)
347 call da_get_vertical_truncation(max_vert_var2, be2_eval_glo(:), be % v2)
348 call da_get_vertical_truncation(max_vert_var3, be3_eval_glo(:), be % v3)
349 call da_get_vertical_truncation(max_vert_var4, be4_eval_glo(:), be % v4)
350
351 be % v5 % mz = 1
352
353 write (unit=stdout,fmt=*) ' '
354
355 else
356
357 ! 3.3 no truncated
358
359 be % v1 % mz = xb % mkz
360 be % v2 % mz = xb % mkz
361 be % v3 % mz = xb % mkz
362 be % v4 % mz = xb % mkz
363 be % v5 % mz = xb % mkz
364
365 end if
366
367 ! 4.0 Initialise control variable space components of header:
368 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
369
370 ! 4.1 Compute the size of the control variables
371
372 be % mix = ix
373 be % mjy = jy
374
375 ! 4.2 Transfer errors to error structure:
376
377 call da_allocate_background_errors(jy, kz, be1_eval_glo, be1_evec_glo, be1_eval_loc, &
378 be1_evec_loc, be % v1)
379 call da_allocate_background_errors(jy, kz, be2_eval_glo, be2_evec_glo, be2_eval_loc, &
380 be2_evec_loc, be % v2)
381 call da_allocate_background_errors(jy, kz, be3_eval_glo, be3_evec_glo, be3_eval_loc, &
382 be3_evec_loc, be % v3)
383 call da_allocate_background_errors(jy, kz, be4_eval_glo, be4_evec_glo, be4_eval_loc, &
384 be4_evec_loc, be % v4)
385
386 ! 4.2.1 transfer the ps_u variance to be % v5:
387
388 call da_allocate_background_errors(jy, 1, be5_eval_glo, be5_evec_glo, be5_eval_loc, &
389 be5_evec_loc, be % v5)
390
391 if (print_detail_be) then
392 write (unit=stdout,fmt='(3x,a,i10)') "b) Finished eigenvector processing!"
393 end if
394
395 ! 5.0 Load the scale lengths
396 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
397
398 ! 5.1 allocate the array for scale lengths
399
400 allocate (be1_rf_lengthscale(1:nk))
401 allocate (be2_rf_lengthscale(1:nk))
402 allocate (be3_rf_lengthscale(1:nk))
403 allocate (be4_rf_lengthscale(1:nk))
404 allocate (be5_rf_lengthscale(1:nk))
405
406 ! 5.2 read in the scale lengths
407
408 read (be_unit) variable
409 read (be_unit) be1_rf_lengthscale
410
411 read (be_unit) variable
412 read (be_unit) be2_rf_lengthscale
413
414 read (be_unit) variable
415 read (be_unit) be3_rf_lengthscale
416
417 read (be_unit) variable
418 read (be_unit) be4_rf_lengthscale
419
420 read (be_unit) variable
421 read (be_unit) be5_rf_lengthscale(1:1)
422 be%v5%name = variable
423
424 ! 5.3 Convert the scale lengths in the real distance (meter)
425
426 be1_rf_lengthscale(1:nk) = be1_rf_lengthscale(1:nk) * xb%ds
427 be2_rf_lengthscale(1:nk) = be2_rf_lengthscale(1:nk) * xb%ds
428 be3_rf_lengthscale(1:nk) = be3_rf_lengthscale(1:nk) * xb%ds
429 be4_rf_lengthscale(1:nk) = be4_rf_lengthscale(1:nk) * xb%ds
430 be5_rf_lengthscale(1:1) = be5_rf_lengthscale(1:1) * xb%ds
431
432 ! 6.0 Perform checks on eigenvectors with be data structure:
433 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
434
435 if (test_statistics) then
436 call da_check_eof_decomposition(be%v1%val_g(:), be%v1%evec_g(:,:),&
437 be%v1%name)
438 call da_check_eof_decomposition(be%v2%val_g(:), be%v2%evec_g(:,:),&
439 be%v2%name)
440 call da_check_eof_decomposition(be%v3%val_g(:), be%v3%evec_g(:,:),&
441 be%v3%name)
442 call da_check_eof_decomposition(be%v4%val_g(:), be%v4%evec_g(:,:),&
443 be%v4%name)
444 end if
445
446 ! 6.2 Close the be unit
447
448 close(be_unit)
449 call da_free_unit(be_unit)
450
451 ! 7.0 Apply empirical and recursive filter rescaling factor:
452 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
453
454 call da_rescale_background_errors(var_scaling1, len_scaling1, &
455 xb % ds, be1_rf_lengthscale, be % v1)
456 call da_rescale_background_errors(var_scaling2, len_scaling2, &
457 xb % ds, be2_rf_lengthscale, be % v2)
458 call da_rescale_background_errors(var_scaling3, len_scaling3, &
459 xb % ds, be3_rf_lengthscale, be % v3)
460 call da_rescale_background_errors(var_scaling4, len_scaling4, &
461 xb % ds, be4_rf_lengthscale, be % v4)
462 call da_rescale_background_errors(var_scaling5, len_scaling5, &
463 xb % ds, be5_rf_lengthscale, be % v5)
464
465 ! 8.0 deallocate input model state:
466 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
467
468 deallocate (be1_eval_loc)
469 deallocate (be2_eval_loc)
470 deallocate (be3_eval_loc)
471 deallocate (be4_eval_loc)
472 deallocate (be5_eval_loc)
473
474 if (vert_corr == vert_corr_2) then
475 deallocate (be1_eval_glo)
476 deallocate (be2_eval_glo)
477 deallocate (be3_eval_glo)
478 deallocate (be4_eval_glo)
479 deallocate (be5_eval_glo)
480
481 deallocate (be1_evec_loc)
482 deallocate (be2_evec_loc)
483 deallocate (be3_evec_loc)
484 deallocate (be4_evec_loc)
485 deallocate (be5_evec_loc)
486
487 deallocate (be1_evec_glo)
488 deallocate (be2_evec_glo)
489 deallocate (be3_evec_glo)
490 deallocate (be4_evec_glo)
491 deallocate (be5_evec_glo)
492
493 end if
494
495 deallocate (bin)
496 deallocate (bin2d)
497
498 if (be % ne > 0) then
499 be % alpha % mz = be % ne
500 be % alpha % name = 'alpha'
501
502 allocate (be5_eval_loc (1:1,1:num_bins2d))
503 allocate (be5_eval_glo(1:1))
504 allocate (be5_evec_loc(1:1,1:1,1:num_bins2d))
505 allocate (be5_evec_glo(1:1,1:1))
506
507 ! Transfer the alpha standard deviation to be % alpha:
508 allocate (be % alpha % val(1:jy,1:be % ne))
509 allocate (be % alpha % evec(1:jy,1:kz,1:be % ne))
510 allocate (be % alpha % rf_alpha(1:be % ne))
511
512 be % alpha % val(1:jy,1:be % ne) = alpha_std_dev
513 be % alpha % evec(1:jy,1:kz,1:be % ne) = 1.0
514 be % alpha % rf_alpha(1:be % ne) = 1.0
515
516 ! Include alpha lengthscale info:
517 allocate (alpha_rf_lengthscale(1:be % ne))
518 alpha_rf_lengthscale(:) = 1000.0 * alpha_corr_scale ! Convert km to m.
519
520 call da_rescale_background_errors(1.0, 1.0, &
521 xb % ds, alpha_rf_lengthscale, be % alpha)
522 else
523 be % alpha % mz = 0
524 end if
525
526 if (trace_use) call da_trace_exit("da_setup_be_regional")
527
528 end subroutine da_setup_be_regional
529
530