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