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