da_setup_be_global.inc
References to this file elsewhere.
1 subroutine da_setup_be_global (be)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Define and allocate components of background errors
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (be_type), intent(out) :: be ! Back. errors structure
10
11 integer :: nrec
12 integer :: max_wave_in ! Dimension of input power
13 integer :: i, j, k, n ! Loop counters
14 ! real, allocatable :: height(:,:,:) ! Height field.
15 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point
16 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point
17 integer :: bin_type ! Type of bin to average over.
18 integer :: num_bins ! Number of bins (3D fields).
19 integer :: num_bins2d ! Number of bins (3D fields).
20 logical :: dummy
21
22 real :: binwidth_lat ! Used if bin_type = 2 (degrees)
23 real :: binwidth_hgt ! Used if bin_type = 2 (m)
24 real :: hgt_min, hgt_max ! Used if bin_type = 2 (m)
25 real :: lat_min, lat_max ! Used if bin_type = 2 (degrees)
26
27 character*10 :: variable
28 integer :: ni, nj, nk, b, be_unit
29 integer, allocatable :: max_wave(:)
30 real, allocatable :: evec_g(:,:)
31 real, allocatable :: eval_g(:)
32 real, allocatable :: evec_loc(:,:,:)
33 real, allocatable :: eval_loc(:,:)
34 real, allocatable :: regcoeff1(:)
35 real, allocatable :: regcoeff2(:,:)
36 real, allocatable :: regcoeff3(:,:,:)
37 real, allocatable :: power(:) ! Temporary power spectrum.
38 real, allocatable :: power2d(:,:) ! Temporary power spectrum.
39
40 if (trace_use) call da_trace_entry("da_setup_be_global")
41
42 call da_message((/"[3.0] Set up background errors (be) for global WRF-Var"/))
43
44 be % max_wave = ide / 2 - 1
45
46 !---------------------------------------------------------------------
47 ! [1] Read in be metadata:
48 !---------------------------------------------------------------------
49
50 call da_get_unit(be_unit)
51 open (unit=be_unit,file="be.dat",status="old",form="unformatted")
52
53 read (be_unit, end= 99, err = 100) ni, nj, nk
54 read (be_unit, err = 100) bin_type
55 read (be_unit, err = 100) lat_min, lat_max, binwidth_lat
56 read (be_unit, err = 100) hgt_min, hgt_max, binwidth_hgt
57 read (be_unit, err = 100) num_bins, num_bins2d
58
59 allocate (bin(1:ni,1:nj,1:nk))
60 allocate (bin2d(1:ni,1:nj))
61 read (be_unit, err = 100) bin(1:ni,1:nj,1:nk)
62 read (be_unit, err = 100) bin2d(1:ni,1:nj)
63
64 if (ni /= ide .or. nj /= jde .or. nk /= kde) then
65 call da_error(__FILE__,__LINE__, &
66 (/"Cannot generate BE at this resolution"/))
67 end if
68
69 !---------------------------------------------------------------------
70 ! [2] Read in be regression coefficients:
71 !---------------------------------------------------------------------
72
73 allocate (regcoeff1(1:num_bins))
74 allocate (regcoeff2(1:nk,1:num_bins2d))
75 allocate (regcoeff3(1:nk,1:nk,1:num_bins2d))
76
77 read (be_unit,end= 99, err = 100) regcoeff1
78 read (be_unit,end= 99, err = 100) regcoeff2
79 read (be_unit,end= 99, err = 100) regcoeff3
80
81 allocate (be%reg_chi(1:jde,1:kde))
82 allocate (be%reg_ps (1:jde,1:kde))
83 allocate (be%reg_t (1:jde,1:kde,1:kde))
84
85 ! Fill regression coeff. array
86 do k=1,kde
87 do j =1, jde
88 b = bin(1,j,k) ! Assumes bins averaged over i direction.
89 be%reg_chi(j,k) = regcoeff1(b)
90 end do
91 end do
92
93 do j=1,jde
94 b = bin2d(1,j)
95 do k=1,kde
96 be%reg_ps(j,k) = regcoeff2(k,b)
97 end do
98 end do
99
100 do j=1,jde
101 b = bin2d(1,j)
102 do i=1,kde
103 do k=1,kde
104 be%reg_t(j,i,k) = regcoeff3(i,k,b)
105 end do
106 end do
107 end do
108
109 deallocate(regcoeff1)
110 deallocate(regcoeff2)
111 deallocate(regcoeff3)
112
113 !---------------------------------------------------------------------
114 ! [3] Read in be vertical eigenmodes:
115 !---------------------------------------------------------------------
116
117 do nrec = 1, 4
118 read (be_unit,end= 99, err = 100) variable
119 read (be_unit,end= 99, err = 100) nk, num_bins2d
120
121 allocate (evec_g(1:nk,1:nk))
122 allocate (eval_g(1:nk))
123 allocate (evec_loc(1:nk,1:nk,num_bins2d))
124 allocate (eval_loc(1:nk,num_bins2d))
125
126 read (be_unit,end= 99, err = 100) evec_g
127 read (be_unit,end= 99, err = 100) eval_g
128 read (be_unit,end= 99, err = 100) evec_loc
129 read (be_unit,end= 99, err = 100) eval_loc
130
131 if (nrec == 1) then
132 be % v1 % name = variable
133 call da_get_bins_info(nj, nk, bin2d, evec_g, eval_g, &
134 evec_loc, eval_loc, max_vert_var1, var_scaling1, be%v1)
135
136 else if (nrec == 2) then
137 be % v2 % name = variable
138 call da_get_bins_info(nj, nk, bin2d, evec_g, eval_g, &
139 evec_loc, eval_loc, max_vert_var2, var_scaling2, be%v2)
140
141 else if (nrec == 3) then
142 be % v3 % name = variable
143 call da_get_bins_info(nj, nk, bin2d, evec_g, eval_g, &
144 evec_loc, eval_loc, max_vert_var3, var_scaling3, be%v3)
145
146 else if (nrec == 4) then
147 be % v4 % name = variable
148 call da_get_bins_info(nj, nk, bin2d, evec_g, eval_g, &
149 evec_loc, eval_loc, max_vert_var4, var_scaling4, be%v4)
150 end if
151
152 deallocate (evec_g)
153 deallocate (eval_g)
154 deallocate (evec_loc)
155 deallocate (eval_loc)
156
157 end do ! loop nrec
158
159 deallocate (bin)
160 deallocate (bin2d)
161
162 !---------------------------------------------------------------------
163 ! [4] Read in be power spectra:
164 !---------------------------------------------------------------------
165
166 allocate(max_wave(1:nk))
167
168 do k = 1, nk
169 read (be_unit) variable
170 read (be_unit) max_wave_in, nrec
171 read (be_unit) dummy ! use to preserve file format
172 if (k == 1) then
173 allocate (power(0:max_wave_in)) ! Temporary.
174 allocate (power2d(0:max_wave_in,1:nk)) ! Temporary.
175 end if
176 read (be_unit) power(0:max_wave_in)
177 power2d(:,k) = power(:)
178
179 ! Truncate power spectra:
180 call da_truncate_spectra(be % max_wave, max_wave_in, power_truncation, &
181 power, max_wave(k))
182 end do
183
184 be % v1 % max_wave = maxval(max_wave(1:nk))
185 write (unit=stdout,fmt='(/3x,3a,i6)') &
186 'Horizontal truncation for ', be % v1 % name, ' = ', be % v1 % max_wave
187 allocate (be % v1 % power(0:be % v1 % max_wave,1:nk))
188 be % v1 % power(0:be % v1 % max_wave,1:nk) = power2d(0:be % v1 % max_wave,1:nk)
189 be % v1 % power(0,1:nk) = len_scaling1 * be % v1 % power(0,1:nk)
190
191 do k = 1, nk
192 read (be_unit) variable
193 read (be_unit) max_wave_in, nrec
194 read (be_unit) dummy ! use to preserve file format
195 read (be_unit) power(0:max_wave_in)
196 power2d(:,k) = power(:)
197
198 ! Truncate power spectra:
199 call da_truncate_spectra (be % max_wave, max_wave_in, power_truncation, &
200 power, max_wave(k))
201 end do
202
203 be % v2 % max_wave = maxval(max_wave(1:nk))
204 write (unit=stdout,fmt='(3x,3a,i6)') &
205 'Horizontal truncation for ', be % v2 % name, ' = ', be % v2 % max_wave
206 allocate (be % v2 % power(0:be % v2 % max_wave,1:nk))
207 be % v2 % power(0:be % v2 % max_wave,1:nk) = power2d(0:be % v2 % max_wave,1:nk)
208 be % v2 % power(0,1:nk) = len_scaling2 * be % v2 % power(0,1:nk)
209
210 do k = 1, nk
211 read (be_unit) variable
212 read (be_unit) max_wave_in, nrec
213 read (be_unit) dummy ! use to preserve file format
214 read (be_unit) power(0:max_wave_in)
215 power2d(:,k) = power(:)
216
217 ! Truncate power spectra:
218 call da_truncate_spectra (be % max_wave, max_wave_in, power_truncation, &
219 power, max_wave(k))
220 end do
221
222 be % v3 % max_wave = maxval(max_wave(1:nk))
223 write(unit=stdout,fmt='(3x,3a,i6)') &
224 'Horizontal truncation for ', be % v3 % name, ' = ', be % v3 % max_wave
225 allocate (be % v3 % power(0:be % v3 % max_wave,1:nk))
226 be % v3 % power(0:be % v3 % max_wave,1:nk) = power2d(0:be % v3 % max_wave,1:nk)
227 be % v3 % power(0,1:nk) = len_scaling3 * be % v3 % power(0,1:nk)
228
229 do k = 1, nk
230 read (be_unit) variable
231 read (be_unit) max_wave_in, nrec
232 read (be_unit) dummy ! use to preserve file format
233 read (be_unit) power(0:max_wave_in)
234 power2d(:,k) = power(:)
235
236 ! Truncate power spectra:
237 call da_truncate_spectra (be % max_wave, max_wave_in, power_truncation, &
238 power, max_wave(k))
239 end do
240
241 be % v4 % max_wave = maxval(max_wave(1:nk))
242 write (unit=stdout,fmt='(3x,3a,i6)') &
243 'Horizontal truncation for ', be % v4 % name, ' = ', be % v4 % max_wave
244 allocate (be % v4 % power(0:be % v4 % max_wave,1:nk))
245 be % v4 % power(0:be % v4 % max_wave,1:nk) = power2d(0:be % v4 % max_wave,1:nk)
246 be % v4 % power(0,1:nk) = len_scaling4 * be % v4 % power(0,1:nk)
247
248 ! ps_u:
249 read (be_unit) variable
250 be % v5 % name = variable
251 be % v5 % mz = 1
252 if (max_vert_var5 <= 0.0) be % v5 % mz = 0
253 read (be_unit) max_wave_in, nrec
254 read (be_unit) dummy ! use to preserve file format
255 read (be_unit) power(0:max_wave_in)
256
257 ! Truncate power spectra:
258 call da_truncate_spectra (be % max_wave, max_wave_in, power_truncation, &
259 power, be % v5 % max_wave)
260
261 write (unit=stdout,fmt='(3x,3a,i6)') &
262 'Horizontal truncation for ', be % v5 % name, ' = ', be % v5 % max_wave
263 allocate (be % v5 % power(0:be % v5 % max_wave,1))
264 be % v5 % power(0:be % v5 % max_wave,1) = power(0:be % v5 % max_wave)
265 be % v5 % power(0,1) = len_scaling5 * be%v5%power(0,1)
266
267 deallocate(power)
268
269 !--------------------------------------------------------------
270 ! [5] Perform checks on eigenvectors:
271 !--------------------------------------------------------------
272
273 if (test_statistics) then
274 call da_check_eof_decomposition(be%v1%val_g(:), be%v1%evec_g(:,:),&
275 be%v1%name)
276 call da_check_eof_decomposition(be%v2%val_g(:), be%v2%evec_g(:,:),&
277 be%v2%name)
278 call da_check_eof_decomposition(be%v3%val_g(:), be%v3%evec_g(:,:),&
279 be%v3%name)
280 call da_check_eof_decomposition(be%v4%val_g(:), be%v4%evec_g(:,:),&
281 be%v4%name)
282 end if
283
284 !--------------------------------------------------------------
285 ! [6] Set up alpha (flow-dependent) control variable "errors":
286 !--------------------------------------------------------------
287
288 if (be % ne > 0) then
289 be % alpha % mz = be % ne
290 be % alpha % name = 'alpha'
291 be % alpha % max_wave = alpha_truncation
292 if (print_detail_be) then
293 write(unit=stdout,fmt='(3x,3a,i6)') &
294 'Horizontal truncation for ', be % alpha % name, ' = ', &
295 be % alpha % max_wave
296 end if
297 allocate (power(0:be % alpha % max_wave))
298 call da_calc_power_spectrum(be % alpha % max_wave, power)
299
300 allocate (be % alpha % power(0:be % alpha % max_wave, be % ne))
301 do n = 1, be % ne
302 be % alpha % power(0:be % alpha % max_wave,n) = power(0:be % alpha % max_wave)
303 end do
304 end if
305
306 deallocate (max_wave)
307
308 write (unit=stdout,fmt='(A)') " "
309
310 close (be_unit)
311 call da_free_unit(be_unit)
312
313 if (trace_use) call da_trace_exit("da_setup_be_global")
314
315 return
316
317 99 write (unit=message(1),fmt='(a, i5)')' Unexpected end on BE-unit = ',be_unit
318 call da_error(__FILE__,__LINE__,message(1:1))
319
320 100 write (unit=message(1),fmt='(a, i5)')' Read error on BE-unit = ',be_unit
321 call da_error(__FILE__,__LINE__,message(1:1))
322
323 ! FIX? daft having these here
324
325 deallocate(power)
326 deallocate(power2d)
327
328 end subroutine da_setup_be_global
329
330