gen_be_ep1.f90
References to this file elsewhere.
1 program gen_be_ep1
2
3 use da_control, only : stderr,stdout, &
4 filename_len, da_advance_cymdh
5 use da_gen_be, only : da_filter_regcoeffs
6 use da_tools1, only : da_get_unit, da_free_unit
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*3 :: ce ! Member index -> character.
15 integer :: ni, nj, nk, nkdum ! Grid dimensions.
16 integer :: i, j, k, member ! Loop counters.
17 integer :: b ! Bin marker.
18 integer :: sdate, cdate, edate ! Starting, current ending dates.
19 integer :: interval ! Period between dates (hours).
20 integer :: ne ! Number of ensemble members.
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 integer :: num_passes ! Recursive filter passes.
25 integer :: count ! Counter.
26 real :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
27 real :: binwidth_lat ! Used if bin_type = 2 (degrees).
28 real :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
29 real :: binwidth_hgt ! Used if bin_type = 2 (m).
30 real :: rf_scale ! Recursive filter scale.
31 real :: count_inv ! 1 / count.
32
33 real, allocatable :: psi(:,:,:) ! psi.
34 real, allocatable :: chi(:,:,:) ! chi.
35 real, allocatable :: temp(:,:,:) ! Temperature.
36 real, allocatable :: rh(:,:,:) ! Relative humidity.
37 real, allocatable :: ps(:,:) ! Surface pressure.
38 real, allocatable :: chi_u(:,:,:) ! chi.
39 real, allocatable :: temp_u(:,:,:) ! Temperature.
40 real, allocatable :: ps_u(:,:) ! Surface pressure.
41
42 real, allocatable :: psi_mnsq(:,:,:) ! psi.
43 real, allocatable :: chi_mnsq(:,:,:) ! chi.
44 real, allocatable :: temp_mnsq(:,:,:) ! Temperature.
45 real, allocatable :: rh_mnsq(:,:,:) ! Relative humidity.
46 real, allocatable :: ps_mnsq(:,:) ! Surface pressure.
47 real, allocatable :: chi_u_mnsq(:,:,:) ! chi.
48 real, allocatable :: temp_u_mnsq(:,:,:) ! Temperature.
49 real, allocatable :: ps_u_mnsq(:,:) ! Surface pressure.
50
51 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
52 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
53
54 real, allocatable :: regcoeff1(:) ! psi/chi regression cooefficient.
55 real, allocatable :: regcoeff2(:,:) ! psi/ps regression cooefficient.
56 real, allocatable :: regcoeff3(:,:,:) ! psi/T regression cooefficient.
57
58 namelist / gen_be_stage2a_nl / start_date, end_date, interval, &
59 ne, num_passes, rf_scale
60
61 integer :: ounit, gen_be_ounit, namelist_unit, iunit
62
63 stderr = 0
64 stdout = 6
65
66 !---------------------------------------------------------------------------------------------
67 write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
68 !---------------------------------------------------------------------------------------------
69
70 call da_get_unit(ounit)
71 call da_get_unit(iunit)
72 call da_get_unit(gen_be_ounit)
73 call da_get_unit(namelist_unit)
74
75 start_date = '2004030312'
76 end_date = '2004033112'
77 interval = 24
78 ne = 1
79 num_passes = 0
80 rf_scale = 1.0
81
82 open(unit=namelist_unit, file='gen_be_stage2a_nl.nl', &
83 form='formatted', status='old', action='read')
84 read(namelist_unit, gen_be_stage2a_nl)
85 close(namelist_unit)
86
87 read(start_date(1:10), fmt='(i10)')sdate
88 read(end_date(1:10), fmt='(i10)')edate
89 write(6,'(4a)')' Computing control variable fields'
90 write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
91 write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
92 write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
93
94 date = start_date
95 cdate = sdate
96
97 !---------------------------------------------------------------------------------------------
98 write(6,'(2a)')' [2] Read regression coefficients and bin information:'
99 !---------------------------------------------------------------------------------------------
100
101 filename = 'be.dat'
102 open (iunit, file = filename, form='unformatted')
103
104 read(iunit)ni, nj, nk
105 read(iunit)bin_type
106 read(iunit)lat_min, lat_max, binwidth_lat
107 read(iunit)hgt_min, hgt_max, binwidth_hgt
108 read(iunit)num_bins, num_bins2d
109
110 allocate( bin(1:ni,1:nj,1:nk) )
111 allocate( bin2d(1:ni,1:nj) )
112 allocate( regcoeff1(1:num_bins) )
113 allocate( regcoeff2(1:nk,1:num_bins2d) )
114 allocate( regcoeff3(1:nk,1:nk,1:num_bins2d) )
115
116 read(iunit)bin(1:ni,1:nj,1:nk)
117 read(iunit)bin2d(1:ni,1:nj)
118 read(iunit)regcoeff1
119 read(iunit)regcoeff2
120 read(iunit)regcoeff3
121
122 close(iunit)
123
124 allocate( psi(1:ni,1:nj,1:nk) )
125 allocate( chi(1:ni,1:nj,1:nk) )
126 allocate( temp(1:ni,1:nj,1:nk) )
127 allocate( rh(1:ni,1:nj,1:nk) )
128 allocate( ps(1:ni,1:nj) )
129 allocate( chi_u(1:ni,1:nj,1:nk) )
130 allocate( temp_u(1:ni,1:nj,1:nk) )
131 allocate( ps_u(1:ni,1:nj) )
132
133 if ( num_passes > 0 ) then
134
135 write(6,'(a,i4,a)')' [3] Apply ', num_passes, ' pass recursive filter to regression coefficients:'
136 call da_filter_regcoeffs( ni, nj, nk, num_bins, num_bins2d, num_passes, rf_scale, bin, &
137 regcoeff1, regcoeff2, regcoeff3 )
138 else
139 write(6,'(a)')' [3] num_passes = 0. Bypassing recursive filtering.'
140 end if
141
142 !---------------------------------------------------------------------------------------------
143 write(6,'(a)')' [4] Read standard fields, and compute unbalanced control variable fields:'
144 !---------------------------------------------------------------------------------------------
145
146 date = start_date
147 cdate = sdate
148
149 do while ( cdate <= edate )
150 write(6,'(a,a)')' Calculating unbalanced fields for date ', date
151
152 do member = 1, ne
153
154 write(ce,'(i3.3)')member
155
156 ! Read psi predictor:
157 variable = 'psi'
158 filename = trim(variable)//'/'//date(1:10)
159 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
160 open (iunit, file = filename, form='unformatted')
161 read(iunit)ni, nj, nk
162 read(iunit)psi
163 close(iunit)
164
165 ! Calculate unbalanced chi:
166 variable = 'chi'
167 filename = trim(variable)//'/'//date(1:10)
168 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
169 open (iunit, file = filename, form='unformatted')
170 read(iunit)ni, nj, nk
171 read(iunit)chi
172 close(iunit)
173
174 do k = 1, nk
175 do j = 1, nj
176 do i = 1, ni
177 b = bin(i,j,k)
178 chi_u(i,j,k) = chi(i,j,k) - regcoeff1(b) * psi(i,j,k)
179 end do
180 end do
181 end do
182
183 variable = 'chi_u'
184 filename = trim(variable)//'/'//date(1:10)
185 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
186 open (ounit, file = filename, form='unformatted')
187 write(ounit)ni, nj, nk
188 write(ounit)chi_u
189 close(ounit)
190
191 ! Calculate unbalanced T:
192 variable = 't'
193 filename = trim(variable)//'/'//date(1:10)
194 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
195 open (iunit, file = filename, form='unformatted')
196 read(iunit)ni, nj, nk
197 read(iunit)temp
198 close(iunit)
199
200 do j = 1, nj
201 do i = 1, ni
202 b = bin2d(i,j)
203 do k = 1, nk
204 temp_u(i,j,k) = temp(i,j,k) - SUM(regcoeff3(k,1:nk,b) * psi(i,j,1:nk))
205 end do
206 end do
207 end do
208
209 variable = 't_u'
210 filename = trim(variable)//'/'//date(1:10)
211 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
212 open (ounit, file = filename, form='unformatted')
213 write(ounit)ni, nj, nk
214 write(ounit)temp_u
215 close(ounit)
216
217 ! Calculate unbalanced ps:
218 variable = 'ps'
219 filename = trim(variable)//'/'//date(1:10)
220 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
221 open (iunit, file = filename, form='unformatted')
222 read(iunit)ni, nj, nkdum
223 read(iunit)ps
224 close(iunit)
225
226 do j = 1, nj
227 do i = 1, ni
228 b = bin2d(i,j)
229 ps_u(i,j) = ps(i,j) - SUM(regcoeff2(1:nk,b) * psi(i,j,1:nk))
230 end do
231 end do
232
233 variable = 'ps_u'
234 filename = trim(variable)//'/'//date(1:10)
235 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
236 open (ounit, file = filename, form='unformatted')
237 write(ounit)ni, nj, 1
238 write(ounit)ps_u
239 close(ounit)
240
241 end do ! End loop over ensemble members.
242
243 !--------------------------------------------------------------------------------------
244 ! Calculate mean control variables (diagnostics):
245 !--------------------------------------------------------------------------------------
246
247 ! Read psi predictor:
248 variable = 'psi'
249 filename = trim(variable)//'/'//date(1:10)
250 filename = trim(filename)//'.'//trim(variable)//'.mean'
251 open (iunit, file = filename, form='unformatted')
252 read(iunit)ni, nj, nk
253 read(iunit)psi
254 close(iunit)
255
256 ! Calculate unbalanced chi:
257 variable = 'chi'
258 filename = trim(variable)//'/'//date(1:10)
259 filename = trim(filename)//'.'//trim(variable)//'.mean'
260 open (iunit, file = filename, form='unformatted')
261 read(iunit)ni, nj, nk
262 read(iunit)chi
263 close(iunit)
264
265 do k = 1, nk
266 do j = 1, nj
267 do i = 1, ni
268 b = bin(i,j,k)
269 chi_u(i,j,k) = chi(i,j,k) - regcoeff1(b) * psi(i,j,k)
270 end do
271 end do
272 end do
273
274 variable = 'chi_u'
275 filename = trim(variable)//'/'//date(1:10)
276 filename = trim(filename)//'.'//trim(variable)//'.mean'
277 open (ounit, file = filename, form='unformatted')
278 write(ounit)ni, nj, nk
279 write(ounit)chi_u
280 close(ounit)
281
282 ! Calculate unbalanced T:
283 variable = 't'
284 filename = trim(variable)//'/'//date(1:10)
285 filename = trim(filename)//'.'//trim(variable)//'.mean'
286 open (iunit, file = filename, form='unformatted')
287 read(iunit)ni, nj, nk
288 read(iunit)temp
289 close(iunit)
290
291 do j = 1, nj
292 do i = 1, ni
293 b = bin2d(i,j)
294 do k = 1, nk
295 temp_u(i,j,k) = temp(i,j,k) - SUM(regcoeff3(k,1:nk,b) * psi(i,j,1:nk))
296 end do
297 end do
298 end do
299
300 variable = 't_u'
301 filename = trim(variable)//'/'//date(1:10)
302 filename = trim(filename)//'.'//trim(variable)//'.mean'
303 open (ounit, file = filename, form='unformatted')
304 write(ounit)ni, nj, nk
305 write(ounit)temp_u
306 close(ounit)
307
308 ! Calculate unbalanced ps:
309 variable = 'ps'
310 filename = trim(variable)//'/'//date(1:10)
311 filename = trim(filename)//'.'//trim(variable)//'.mean'
312 open (iunit, file = filename, form='unformatted')
313 read(iunit)ni, nj, nkdum
314 read(iunit)ps
315 close(iunit)
316
317 do j = 1, nj
318 do i = 1, ni
319 b = bin2d(i,j)
320 ps_u(i,j) = ps(i,j) - SUM(regcoeff2(1:nk,b) * psi(i,j,1:nk))
321 end do
322 end do
323
324 variable = 'ps_u'
325 filename = trim(variable)//'/'//date(1:10)
326 filename = trim(filename)//'.'//trim(variable)//'.mean'
327 open (ounit, file = filename, form='unformatted')
328 write(ounit)ni, nj, 1
329 write(ounit)ps_u
330 close(ounit)
331
332 ! Calculate next date:
333 call da_advance_cymdh( date, interval, new_date )
334 date = new_date
335 read(date(1:10), fmt='(i10)')cdate
336 end do ! End loop over times.
337
338 deallocate( bin )
339 deallocate( bin2d )
340 deallocate( regcoeff1 )
341 deallocate( regcoeff2 )
342 deallocate( regcoeff3 )
343
344 !---------------------------------------------------------------------------------------------
345 write(6,'(a)')' [5] Compute mean square statistics:'
346 !---------------------------------------------------------------------------------------------
347
348 allocate( psi_mnsq(1:ni,1:nj,1:nk) )
349 allocate( chi_mnsq(1:ni,1:nj,1:nk) )
350 allocate( temp_mnsq(1:ni,1:nj,1:nk) )
351 allocate( rh_mnsq(1:ni,1:nj,1:nk) )
352 allocate( ps_mnsq(1:ni,1:nj) )
353 allocate( chi_u_mnsq(1:ni,1:nj,1:nk) )
354 allocate( temp_u_mnsq(1:ni,1:nj,1:nk) )
355 allocate( ps_u_mnsq(1:ni,1:nj) )
356
357 date = start_date
358 cdate = sdate
359
360 do while ( cdate <= edate )
361 count = 0
362
363 psi_mnsq = 0.0
364 chi_mnsq = 0.0
365 temp_mnsq = 0.0
366 rh_mnsq = 0.0
367 ps_mnsq = 0.0
368 chi_u_mnsq = 0.0
369 temp_u_mnsq = 0.0
370 ps_u_mnsq = 0.0
371
372 do member = 1, ne
373 write(ce,'(i3.3)')member
374 count = count + 1
375 count_inv = 1.0 / real(count)
376
377 variable = 'psi'
378 filename = trim(variable)//'/'//date(1:10)
379 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
380 open (iunit, file = filename, form='unformatted')
381 read(iunit)ni, nj, nk
382 read(iunit)psi
383 close(iunit)
384
385 variable = 'chi'
386 filename = trim(variable)//'/'//date(1:10)
387 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
388 open (iunit, file = filename, form='unformatted')
389 read(iunit)ni, nj, nk
390 read(iunit)chi
391 close(iunit)
392
393 variable = 't'
394 filename = trim(variable)//'/'//date(1:10)
395 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
396 open (iunit, file = filename, form='unformatted')
397 read(iunit)ni, nj, nk
398 read(iunit)temp
399 close(iunit)
400
401 variable = 'rh'
402 filename = trim(variable)//'/'//date(1:10)
403 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
404 open (iunit, file = filename, form='unformatted')
405 read(iunit)ni, nj, nk
406 read(iunit)rh
407 close(iunit)
408
409 variable = 'ps'
410 filename = trim(variable)//'/'//date(1:10)
411 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
412 open (iunit, file = filename, form='unformatted')
413 read(iunit)ni, nj, nkdum
414 read(iunit)ps
415 close(iunit)
416
417 variable = 'chi_u'
418 filename = trim(variable)//'/'//date(1:10)
419 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
420 open (iunit, file = filename, form='unformatted')
421 read(iunit)ni, nj, nk
422 read(iunit)chi_u
423 close(iunit)
424
425 variable = 't_u'
426 filename = trim(variable)//'/'//date(1:10)
427 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
428 open (iunit, file = filename, form='unformatted')
429 read(iunit)ni, nj, nk
430 read(iunit)temp_u
431 close(iunit)
432
433 variable = 'ps_u'
434 filename = trim(variable)//'/'//date(1:10)
435 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
436 open (iunit, file = filename, form='unformatted')
437 read(iunit)ni, nj, nkdum
438 read(iunit)ps
439 close(iunit)
440
441 ! Calculate accumulating mean:
442 psi_mnsq = ( real( count-1 ) * psi_mnsq + psi * psi ) * count_inv
443 chi_mnsq = ( real( count-1 ) * chi_mnsq + chi * chi ) * count_inv
444 temp_mnsq = ( real( count-1 ) * temp_mnsq + temp * temp ) * count_inv
445 rh_mnsq = ( real( count-1 ) * rh_mnsq + rh * rh ) * count_inv
446 ps_mnsq = ( real( count-1 ) * ps_mnsq + ps * ps ) * count_inv
447 chi_u_mnsq = ( real( count-1 ) * chi_u_mnsq + chi_u * chi_u ) * count_inv
448 temp_u_mnsq = ( real( count-1 ) * temp_u_mnsq + temp_u * temp_u ) * count_inv
449 ps_u_mnsq = ( real( count-1 ) * ps_u_mnsq + ps_u * ps_u ) * count_inv
450 end do ! End loop over ensemble members.
451
452 psi_mnsq = sqrt(psi_mnsq) ! Convert mnsq to stdv (mean=0)
453 chi_mnsq = sqrt(chi_mnsq) ! Convert mnsq to stdv (mean=0)
454 temp_mnsq = sqrt(temp_mnsq) ! Convert mnsq to stdv (mean=0)
455 rh_mnsq = sqrt(rh_mnsq) ! Convert mnsq to stdv (mean=0)
456 ps_mnsq = sqrt(ps_mnsq) ! Convert mnsq to stdv (mean=0)
457 chi_u_mnsq = sqrt(chi_u_mnsq) ! Convert mnsq to stdv (mean=0)
458 temp_u_mnsq = sqrt(temp_u_mnsq) ! Convert mnsq to stdv (mean=0)
459 ps_u_mnsq = sqrt(ps_u_mnsq) ! Convert mnsq to stdv (mean=0)
460
461 ! Write mean square statistics:
462 filename = 'psi/'//date(1:10)//'.psi.stdv'
463 open (gen_be_ounit, file = filename, form='unformatted')
464 write(gen_be_ounit)ni, nj, nk
465 write(gen_be_ounit)psi_mnsq
466 close(gen_be_ounit)
467
468 filename = 'chi/'//date(1:10)//'.chi.stdv'
469 open (gen_be_ounit, file = filename, form='unformatted')
470 write(gen_be_ounit)ni, nj, nk
471 write(gen_be_ounit)chi_mnsq
472 close(gen_be_ounit)
473
474 filename = 't/'//date(1:10)//'.t.stdv'
475 open (gen_be_ounit, file = filename, form='unformatted')
476 write(gen_be_ounit)ni, nj, nk
477 write(gen_be_ounit)temp_mnsq
478 close(gen_be_ounit)
479
480 filename = 'rh/'//date(1:10)//'.rh.stdv'
481 open (gen_be_ounit, file = filename, form='unformatted')
482 write(gen_be_ounit)ni, nj, nk
483 write(gen_be_ounit)rh_mnsq
484 close(gen_be_ounit)
485
486 filename = 'ps/'//date(1:10)//'.ps.stdv'
487 open (gen_be_ounit, file = filename, form='unformatted')
488 write(gen_be_ounit)ni, nj, nk
489 write(gen_be_ounit)ps_mnsq
490 close(gen_be_ounit)
491
492 filename = 'chi_u/'//date(1:10)//'.chi_u.stdv'
493 open (gen_be_ounit, file = filename, form='unformatted')
494 write(gen_be_ounit)ni, nj, nk
495 write(gen_be_ounit)chi_u_mnsq
496 close(gen_be_ounit)
497
498 filename = 't_u/'//date(1:10)//'.t_u.stdv'
499 open (gen_be_ounit, file = filename, form='unformatted')
500 write(gen_be_ounit)ni, nj, nk
501 write(gen_be_ounit)temp_u_mnsq
502 close(gen_be_ounit)
503
504 filename = 'ps_u/'//date(1:10)//'.ps_u.stdv'
505 open (gen_be_ounit, file = filename, form='unformatted')
506 write(gen_be_ounit)ni, nj, nk
507 write(gen_be_ounit)ps_u_mnsq
508 close(gen_be_ounit)
509
510 ! Calculate next date:
511 call da_advance_cymdh( date, interval, new_date )
512 date = new_date
513 read(date(1:10), fmt='(i10)')cdate
514 end do ! End loop over times.
515
516 deallocate( psi )
517 deallocate( chi )
518 deallocate( temp )
519 deallocate( rh )
520 deallocate( ps )
521 deallocate( chi_u )
522 deallocate( temp_u )
523 deallocate( ps_u )
524 deallocate( psi_mnsq )
525 deallocate( chi_mnsq )
526 deallocate( temp_mnsq )
527 deallocate( rh_mnsq )
528 deallocate( ps_mnsq )
529 deallocate( chi_u_mnsq )
530 deallocate( temp_u_mnsq )
531 deallocate( ps_u_mnsq )
532
533 call da_free_unit(ounit)
534 call da_free_unit(iunit)
535 call da_free_unit(gen_be_ounit)
536 call da_free_unit(namelist_unit)
537
538 end program gen_be_ep1
539