gen_be_stage2a.f90
References to this file elsewhere.
1 program gen_be_stage2a
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_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 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 :: rf_scale ! Recursive filter scale.
29
30 real, allocatable :: psi(:,:,:) ! psi.
31 real, allocatable :: chi(:,:,:) ! chi.
32 real, allocatable :: temp(:,:,:) ! Temperature.
33 real, allocatable :: ps(:,:) ! Surface pressure.
34 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
35 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
36
37 real, allocatable :: regcoeff1(:) ! psi/chi regression cooefficient
38 real, allocatable :: regcoeff2(:,:) ! psi/ps regression cooefficient.
39 real, allocatable :: regcoeff3(:,:,:) ! psi/T regression cooefficient.
40
41 namelist / gen_be_stage2a_nl / start_date, end_date, interval, &
42 ne, num_passes, rf_scale
43
44 integer :: ounit,iunit,namelist_unit
45
46 stderr = 0
47 stdout = 6
48
49 !---------------------------------------------------------------------------------------------
50 write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
51 !---------------------------------------------------------------------------------------------
52
53 call da_get_unit(ounit)
54 call da_get_unit(iunit)
55 call da_get_unit(namelist_unit)
56
57 start_date = '2004030312'
58 end_date = '2004033112'
59 interval = 24
60 ne = 1
61 num_passes = 0
62 rf_scale = 1.0
63
64 open(unit=namelist_unit, file='gen_be_stage2a_nl.nl', &
65 form='formatted', status='old', action='read')
66 read(namelist_unit, gen_be_stage2a_nl)
67 close(namelist_unit)
68
69 read(start_date(1:10), fmt='(i10)')sdate
70 read(end_date(1:10), fmt='(i10)')edate
71 write(6,'(4a)')' Computing control variable fields'
72 write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
73 write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
74 write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
75
76 date = start_date
77 cdate = sdate
78
79 !---------------------------------------------------------------------------------------------
80 write(6,'(2a)')' [2] Read regression coefficients and bin information:'
81 !---------------------------------------------------------------------------------------------
82
83 filename = 'gen_be_stage2.dat'
84 open (iunit, file = filename, form='unformatted')
85 read(iunit)ni, nj, nk
86 read(iunit)num_bins, num_bins2d
87
88 allocate( bin(1:ni,1:nj,1:nk) )
89 allocate( bin2d(1:ni,1:nj) )
90 allocate( psi(1:ni,1:nj,1:nk) )
91 allocate( chi(1:ni,1:nj,1:nk) )
92 allocate( temp(1:ni,1:nj,1:nk) )
93 allocate( ps(1:ni,1:nj) )
94 allocate( regcoeff1(1:num_bins) )
95 allocate( regcoeff2(1:nk,1:num_bins2d) )
96 allocate( regcoeff3(1:nk,1:nk,1:num_bins2d) )
97
98 read(iunit)regcoeff1
99 read(iunit)regcoeff2
100 read(iunit)regcoeff3
101 close(iunit)
102
103 ! Read bin info:
104 filename = 'bin.data'
105 open (iunit, file = filename, form='unformatted')
106 read(iunit)bin_type
107 read(iunit)lat_min, lat_max, binwidth_lat
108 read(iunit)hgt_min, hgt_max, binwidth_hgt
109 read(iunit)num_bins, num_bins2d
110 read(iunit)bin(1:ni,1:nj,1:nk)
111 read(iunit)bin2d(1:ni,1:nj)
112 close(iunit)
113
114 if ( num_passes > 0 ) then
115
116 !---------------------------------------------------------------------------------------------
117 write(6,'(a,i4,a)')' [3] Apply ', num_passes, ' pass recursive filter to regression coefficients:'
118 !---------------------------------------------------------------------------------------------
119 call da_filter_regcoeffs( ni, nj, nk, num_bins, num_bins2d, num_passes, rf_scale, bin, &
120 regcoeff1, regcoeff2, regcoeff3 )
121 else
122 write(6,'(a)')' [3] num_passes = 0. Bypassing recursive filtering.'
123 end if
124
125 !---------------------------------------------------------------------------------------------
126 write(6,'(a)')' [4] Read standard fields, and compute control variable fields:'
127 !---------------------------------------------------------------------------------------------
128
129 date = start_date
130 cdate = sdate
131
132 do while ( cdate <= edate )
133 write(6,'(a,a)')' Calculating unbalanced fields for date ', date
134
135 do member = 1, ne
136
137 write(ce,'(i3.3)')member
138
139 ! Read psi predictor:
140 variable = 'psi'
141 filename = trim(variable)//'/'//date(1:10)
142 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
143 open (iunit, file = filename, form='unformatted')
144 read(iunit)ni, nj, nk
145 read(iunit)psi
146 close(iunit)
147
148 ! Calculate unbalanced chi:
149 variable = 'chi'
150 filename = trim(variable)//'/'//date(1:10)
151 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
152 open (iunit, file = filename, form='unformatted')
153 read(iunit)ni, nj, nk
154 read(iunit)chi
155 close(iunit)
156
157 do k = 1, nk
158 do j = 1, nj
159 do i = 1, ni
160 b = bin(i,j,k)
161 chi(i,j,k) = chi(i,j,k) - regcoeff1(b) * psi(i,j,k)
162 end do
163 end do
164 end do
165
166 variable = 'chi_u'
167 filename = trim(variable)//'/'//date(1:10)
168 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
169 open (ounit, file = filename, form='unformatted')
170 write(ounit)ni, nj, nk
171 write(ounit)chi
172 close(ounit)
173
174 ! Calculate unbalanced T:
175 variable = 't'
176 filename = trim(variable)//'/'//date(1:10)
177 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
178 open (iunit, file = filename, form='unformatted')
179 read(iunit)ni, nj, nk
180 read(iunit)temp
181 close(iunit)
182
183 do j = 1, nj
184 do i = 1, ni
185 b = bin2d(i,j)
186 do k = 1, nk
187 temp(i,j,k) = temp(i,j,k) - SUM(regcoeff3(k,1:nk,b) * psi(i,j,1:nk))
188 end do
189 end do
190 end do
191
192 variable = 't_u'
193 filename = trim(variable)//'/'//date(1:10)
194 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
195 open (ounit, file = filename, form='unformatted')
196 write(ounit)ni, nj, nk
197 write(ounit)temp
198 close(ounit)
199
200 ! Calculate unbalanced ps:
201 variable = 'ps'
202 filename = trim(variable)//'/'//date(1:10)
203 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
204 open (iunit, file = filename, form='unformatted')
205 read(iunit)ni, nj, nkdum
206 read(iunit)ps
207 close(iunit)
208
209 do j = 1, nj
210 do i = 1, ni
211 b = bin2d(i,j)
212 ps(i,j) = ps(i,j) - SUM(regcoeff2(1:nk,b) * psi(i,j,1:nk))
213 end do
214 end do
215
216 variable = 'ps_u'
217 filename = trim(variable)//'/'//date(1:10)
218 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
219 open (ounit, file = filename, form='unformatted')
220 write(ounit)ni, nj, 1
221 write(ounit)ps
222 close(ounit)
223
224 end do ! End loop over ensemble members.
225
226 ! Calculate next date:
227 call da_advance_cymdh( date, interval, new_date )
228 date = new_date
229 read(date(1:10), fmt='(i10)')cdate
230 end do ! End loop over times.
231
232 end program gen_be_stage2a
233