gen_be_stage3.f90
References to this file elsewhere.
1 program gen_be_stage3
2
3 use da_control, only : stderr, stdout, filename_len, vertical_ip
4 use da_gen_be, only : da_eof_decomposition_test, da_eof_decomposition, &
5 da_transform_vptovv, da_create_bins
6 use da_tools_serial, only : da_get_unit, da_advance_cymdh
7
8 implicit none
9
10 character*10 :: start_date, end_date ! Starting and ending dates.
11 character*10 :: date, new_date ! Current date (ccyymmddhh).
12 character*10 :: variable ! Variable name
13 character(len=filename_len) :: filename ! Input filename.
14 character*2 :: ck ! Level index -> character.
15 character*3 :: ce ! Member index -> character.
16 integer :: ni, nj, nk ! Dimensions read in.
17 integer :: sdate, cdate, edate ! Starting, current ending dates.
18 integer :: interval ! Period between dates (hours).
19 integer :: ne ! Number of ensemble members.
20 integer :: i, j, k, k1, k2, b, member ! Loop counters.
21 integer :: bin_type ! Type of bin to average over.
22 integer :: num_bins ! Number of bins (3D fields).
23 integer :: num_bins2d ! Number of bins (2D fields).
24 real :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
25 real :: binwidth_lat ! Used if bin_type = 2 (degrees).
26 real :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
27 real :: binwidth_hgt ! Used if bin_type = 2 (m).
28 real :: inv_nij ! 1 / (ni*nj).
29 real :: mean_field ! Mean field.
30 real :: coeffa, coeffb ! Accumulating mean coefficients.
31 logical :: first_time ! True if first file.
32 logical :: testing_eofs ! True if testing EOF decomposition.
33 logical :: use_global_eofs ! True if projected data uses global EOFs.
34 logical :: data_on_levels ! True if output level data (diagnostic).
35 logical :: twod_field ! True if 2D field.
36
37 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
38 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
39 integer, allocatable:: bin_pts2d(:) ! Number of points in bin (2D fields).
40 real, allocatable :: latitude(:,:) ! Latitude (degrees, from south).
41 real, allocatable :: height(:,:,:) ! Height field.
42 real, allocatable :: field(:,:,:) ! Input field.
43 real, allocatable :: field_out(:,:,:) ! Field projected into EOF space.
44 real, allocatable :: vertical_wgt(:,:,:) ! Inner product for vertical transform.
45 real, allocatable :: bv(:,:,:) ! Vertical BE for this time.
46 real, allocatable :: work(:,:) ! EOF work array.
47 real, allocatable :: e_vec_loc(:,:,:) ! Latitudinally varying eigenvectors.
48 real, allocatable :: e_val_loc(:,:) ! Latitudinally varying eigenvalues.
49 real, allocatable :: e_vec(:,:) ! Domain-averaged eigenvectors.
50 real, allocatable :: e_val(:) ! Domain-averaged eigenvalues.
51 real, allocatable :: evec(:,:,:) ! Gridpoint eigenvectors.
52 real, allocatable :: eval(:,:) ! Gridpoint sqrt(eigenvalues).
53
54 namelist / gen_be_stage3_nl / start_date, end_date, interval, variable, &
55 ne, bin_type, &
56 lat_min, lat_max, binwidth_lat, &
57 hgt_min, hgt_max, binwidth_hgt, &
58 testing_eofs, use_global_eofs, data_on_levels
59
60 integer :: ounit,iunit,namelist_unit
61
62 stderr = 0
63 stdout = 6
64
65 write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
66
67 call da_get_unit(ounit)
68 call da_get_unit(iunit)
69 call da_get_unit(namelist_unit)
70
71
72 vertical_ip = 0
73 twod_field = .false.
74
75 start_date = '2004030312'
76 end_date = '2004033112'
77 interval = 24
78 variable = 'psi'
79 ne = 1
80 bin_type = 5
81 lat_min = -90.0
82 lat_max = 90.0
83 binwidth_lat = 10.0
84 hgt_min = 0.0
85 hgt_max = 20000.0
86 binwidth_hgt = 1000.0
87 testing_eofs = .true.
88 use_global_eofs = .true.
89 data_on_levels = .false.
90
91 open(unit=namelist_unit, file='gen_be_stage3_nl.nl', &
92 form='formatted', status='old', action='read')
93 read(namelist_unit, gen_be_stage3_nl)
94 close(namelist_unit)
95
96 read(start_date(1:10), fmt='(i10)')sdate
97 read(end_date(1:10), fmt='(i10)')edate
98 write(6,'(4a)')' Computing vertical error statistics for dates ', start_date, ' to ', end_date
99 write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
100 write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
101
102 date = start_date
103 cdate = sdate
104 first_time = .true.
105
106 if ( .not. data_on_levels ) then ! Can bypass [2] and [3] if outputing on levels.
107
108 write(6,'(2a)')' [2] Process fields for variable ', variable
109
110 do while ( cdate <= edate )
111 do member = 1, ne
112 write(6,'(5a,i4)')' Processing data for date ', date, ', variable ', trim(variable), &
113 ' and member ', member
114
115 write(ce,'(i3.3)')member
116
117 ! Read Full-fields:
118 filename = 'fullflds'//'/'//date(1:10)//'.'//'fullflds'//'.e'//ce
119 open (iunit, file = filename, form='unformatted')
120
121 read(iunit)ni, nj, nk
122 if ( first_time ) then
123 write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk
124 allocate( latitude(1:ni,1:nj) )
125 allocate( height(1:ni,1:nj,1:nk) )
126 allocate( bin(1:ni,1:nj,1:nk) )
127 allocate( bin2d(1:ni,1:nj) )
128 end if
129 read(iunit)latitude
130 read(iunit)height
131
132 close(iunit)
133
134 if ( first_time ) then
135 ! Create and sort into bins:
136 call da_create_bins( ni, nj, nk, bin_type, num_bins, num_bins2d, bin, bin2d, &
137 lat_min, lat_max, binwidth_lat, &
138 hgt_min, hgt_max, binwidth_hgt, latitude, height )
139
140 allocate( bin_pts2d(1:num_bins2d) )
141 bin_pts2d(:) = 0
142 end if
143
144 filename = trim(variable)//'/'//date(1:10)
145 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
146 if ( trim(variable) == 'ps_u' .or. trim(variable) == 'ps' ) then ! 2D field
147 twod_field = .true.
148 filename = trim(filename)//'.01'
149 end if
150
151 open (iunit, file = filename, form='unformatted')
152 read(iunit)ni, nj, nk
153
154 if ( first_time ) then
155 inv_nij = 1.0 / real(ni*nj)
156 allocate( field(1:ni,1:nj,1:nk) )
157 allocate( bv(1:nk,1:nk,1:num_bins2d) )
158 bv(:,:,:) = 0.0
159 first_time = .false.
160 end if
161
162 read(iunit)field
163 close(iunit)
164
165 ! Remove mean field:
166 do k = 1, nk
167 mean_field = sum(field(1:ni,1:nj,k)) * inv_nij
168 field(1:ni,1:nj,k) = field(1:ni,1:nj,k) - mean_field
169 end do
170
171 do j = 1, nj
172 do i = 1, ni
173 b = bin2d(i,j)
174 bin_pts2d(b) = bin_pts2d(b) + 1
175 coeffa = 1.0 / real(bin_pts2d(b))
176 coeffb = real(bin_pts2d(b)-1) * coeffa
177 do k1 = 1, nk
178 do k2 = 1, k1
179 bv(k1,k2,b) = coeffb * bv(k1,k2,b) + coeffa * field(i,j,k1) * field(i,j,k2)
180 end do
181 end do
182 end do
183 end do
184 end do ! End loop over ensemble members.
185
186 ! Calculate next date:
187 call da_advance_cymdh( date, interval, new_date )
188 date = new_date
189 read(date(1:10), fmt='(i10)')cdate
190 end do ! End loop over times.
191
192 ! Fill in upper-right part of BE matrix by symmetry:
193
194 do b = 1, num_bins2d
195 do k1 = 1, nk
196 do k2 = k1+1, nk ! Symmetry.
197 bv(k1,k2,b) = bv(k2,k1,b)
198 end do
199 end do
200 end do
201
202 write(6,'(2a)')' [3] Calculate eigenvectors and eigenvalues for variable ', variable
203
204 allocate( work(1:nk,1:nk) )
205 allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) )
206 allocate( e_val_loc(1:nk,1:num_bins2d) )
207 allocate( e_vec(1:nk,1:nk) )
208 allocate( e_val(1:nk) )
209
210 ! Latitudinally varying BE decomposition:
211 do b = 1, num_bins2d
212 write(6,'(2(a,i6))')' Calculate eigenvectors and eigenvalues for bin ', b, &
213 ' of ', num_bins2d
214
215 work(1:nk,1:nk) = bv(1:nk,1:nk,b)
216 call da_eof_decomposition( nk, work, e_vec, e_val )
217 e_vec_loc(1:nk,1:nk,b) = e_vec(1:nk,1:nk)
218 e_val_loc(1:nk,b) = e_val(1:nk)
219 end do
220
221 ! Domain-averaged BE decomposition:
222 work(1:nk,1:nk) = 0.0
223 do b = 1, num_bins2d
224 work(1:nk,1:nk) = work(1:nk,1:nk) + bv(1:nk,1:nk,b)
225 end do
226 work(1:nk,1:nk) = work(1:nk,1:nk) / real( num_bins2d )
227
228 call da_eof_decomposition( nk, work, e_vec, e_val )
229
230 if ( testing_eofs ) then
231 call da_eof_decomposition_test( nk, work, e_vec, e_val )
232 end if
233
234 ! Output eigenvectors, eigenvalues for use in WRF_Var:
235 filename = 'gen_be_stage3.'//trim(variable)//'.dat'
236 open (ounit, file = filename, form='unformatted')
237 write(ounit)variable
238 write(ounit)nk, num_bins2d
239 write(ounit)e_vec
240 write(ounit)e_val
241 write(ounit)e_vec_loc
242 write(ounit)e_val_loc
243 close(ounit)
244
245 ! Decide on local or domain-averaged EOFs for horizontal correlations:
246 if ( use_global_eofs ) then
247 do b = 1, num_bins2d
248 e_vec_loc(1:nk,1:nk,b) = e_vec(1:nk,1:nk)
249 e_val_loc(1:nk,b) = e_val(1:nk)
250 end do
251 end if
252
253 ! Map binned eigenvectors to x, y grid, and take sqrt(this is used in WRF-Var):
254 allocate( evec(1:nj,1:nk,1:nk) )
255 allocate( eval(1:nj,1:nk) )
256
257 do j = 1, nj
258 do i = 1, ni
259 b = bin2d(i,j)
260 evec(j,1:nk,1:nk) = e_vec_loc(1:nk,1:nk,b)
261 eval(j,1:nk) = sqrt(e_val_loc(1:nk,b))
262 end do
263 end do
264
265 end if ! End bypass [2] and [3] if outputing on levels.
266
267 write(6,'(2a)')' [4] Transform perturbations (or not), and output.'
268
269 date = start_date
270 cdate = sdate
271 first_time = .true.
272
273 if ( .not. twod_field ) then
274 do while ( cdate <= edate )
275 do member = 1, ne
276 write(6,'(5a,i4)')' Date = ', date, ', variable ', trim(variable), &
277 ' and member ', member
278
279 write(ce,'(i3.3)')member
280
281 filename = trim(variable)//'/'//date(1:10)
282 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
283
284 open (iunit, file = filename, form='unformatted')
285 read(iunit)ni, nj, nk
286
287 if ( first_time ) then
288 if ( data_on_levels) allocate( latitude(1:ni,1:nj) ) ! Not allocated earlier.
289 if ( data_on_levels) allocate( field(1:ni,1:nj,1:nk) ) ! Not allocated earlier.
290 allocate( field_out(1:ni,1:nj,1:nk) )
291 allocate( vertical_wgt(1:ni,1:nj,1:nk) )
292 vertical_wgt(1:ni,1:nj,1:nk) = 1.0 ! vertical_ip = 0 hardwired.
293 first_time = .false.
294 end if
295
296 read(iunit)field
297 close(iunit)
298
299 if ( data_on_levels ) then
300 ! Keep data on vertical levels:
301 field_out(:,:,:) = field(:,:,:)
302 else
303 ! Project fields onto vertical modes:
304 call da_transform_vptovv( evec, eval, vertical_wgt, &
305 field, field_out, nk, &
306 1, nk, & ! WRF ids, ide etc.
307 1, ni, 1, nj, 1, nk, & ! WRF ims, ime etc.
308 1, ni, 1, nj, 1, nk ) ! WRF its, ite etc.
309 end if
310
311 ! Output fields (split into 2D files to allow parallel horizontal treatment):
312
313 do k = 1, nk
314 write(ck,'(i2)')k
315 if ( k < 10 ) ck = '0'//ck(2:2)
316 ! Assumes variable directory has been created by script:
317 filename = trim(variable)//'/'//date(1:10)//'.'//trim(variable)
318 filename = trim(filename)//'.e'//ce//'.'//ck
319 open (ounit, file = filename, form='unformatted')
320 write(ounit)ni, nj, k
321 write(ounit)field_out(1:ni,1:nj,k)
322 close(ounit)
323 end do
324 end do ! End loop over members.
325
326 ! Calculate next date:
327 call da_advance_cymdh( date, interval, new_date )
328 date = new_date
329 read(date(1:10), fmt='(i10)')cdate
330 end do
331 end if
332
333 end program gen_be_stage3
334