gen_be_stage1.f90
References to this file elsewhere.
1 program gen_be_stage1
2 !
3 !----------------------------------------------------------------------
4 ! Purpose : to remove the binned mean from the perturbation fields.
5 !
6 ! Input : binary files: "pert.ccyymmddhh.e"ce for ENS or
7 ! "pert.ccyymmddhh.e001" for NMC.
8 !
9 ! Output : binary files for use of the gen_be_stage2:
10 !
11 !----------------------------------------------------------------------
12 !
13 use da_control, only : stderr, stdout, filename_len
14 use da_tools_serial, only : da_get_unit,da_advance_cymdh
15 use da_gen_be, only : da_create_bins
16
17 implicit none
18
19 character*10 :: start_date, end_date ! Starting and ending dates (ccyymmddhh).
20 character*10 :: date, new_date ! Current date (ccyymmddhh).
21 character*10 :: variable ! Variable name.
22 character*3 :: be_method ! Be method (NMC, or ENS)
23 character*3 :: ce ! Ensemble member index.
24 character(len=filename_len) :: dat_dir ! Input data directory.
25 character(len=filename_len) :: filename ! Input filename.
26 integer :: count ! Counter.
27 integer :: ni, nj, nk ! Dimensions read in.
28 integer :: member ! Loop counter
29 integer :: sdate, cdate, edate ! Starting, current ending dates.
30 integer :: interval ! Interval between file times (hours).
31 integer :: ne ! Number of ensemble members.
32 integer :: bin_type ! Type of bin to average over.
33 integer :: num_bins ! Number of bins (3D fields).
34 integer :: num_bins2d ! Number of bins (2D fields).
35 real :: count_inv ! 1 / count.
36 real :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
37 real :: binwidth_lat ! Used if bin_type = 2 (degrees).
38 real :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
39 real :: binwidth_hgt ! Used if bin_type = 2 (m).
40
41 real, allocatable :: ps_prime(:,:) ! Surface pressure perturbation.
42 real, allocatable :: t_prime(:,:,:) ! Temperature perturbation.
43 real, allocatable :: psi_prime(:,:,:) ! Streamfunction perturbation.
44 real, allocatable :: chi_prime(:,:,:) ! Velocity Potential perturbation.
45 real, allocatable :: rh_prime(:,:,:) ! Relative Humidity Perturbation.
46 real, allocatable :: height(:,:,:) ! Geopotential height.
47 real, allocatable :: latitude(:,:) ! Latitude (radians)
48 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
49 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
50 real, allocatable :: psi_mean(:,:,:) ! Mean field.
51 real, allocatable :: chi_mean(:,:,:) ! Mean field.
52 real, allocatable :: t_mean(:,:,:) ! Mean field.
53 real, allocatable :: rh_mean(:,:,:) ! Mean field.
54 real, allocatable :: ps_mean(:,:) ! Mean field.
55
56 namelist / gen_be_stage1_nl / start_date, end_date, interval, &
57 be_method, ne, bin_type, &
58 lat_min, lat_max, binwidth_lat, &
59 hgt_min, hgt_max, binwidth_hgt, dat_dir
60
61 integer :: ounit,iunit,namelist_unit
62
63 stderr = 0
64 stdout = 6
65
66 write(unit=stdout,fmt='(a)') &
67 ' [1] Initialize namelist variables and other scalars.'
68
69 call da_get_unit(ounit)
70 call da_get_unit(iunit)
71 call da_get_unit(namelist_unit)
72
73 start_date = '2004030312'
74 end_date = '2004033112'
75 interval = 24
76 be_method = 'NMC'
77 ne = 1
78 bin_type = 5 ! 0 = Every pt, 1 = x direction, 2 = latitude, ....
79 lat_min = -90.0
80 lat_max = 90.0
81 binwidth_lat = 10.0
82 hgt_min = 0.0
83 hgt_max = 20000.0
84 binwidth_hgt = 1000.0
85 dat_dir = '/data2/hcshin/youn/DIFF63'
86
87 open(unit=namelist_unit, file='gen_be_stage1_nl.nl', &
88 form='formatted', status='old', action='read')
89 read(namelist_unit, gen_be_stage1_nl)
90 close(namelist_unit)
91
92 if ( be_method /= "ENS" ) ne = 1
93
94 read(start_date(1:10), fmt='(i10)')sdate
95 read(end_date(1:10), fmt='(i10)')edate
96 write(6,'(4a)')' Computing statistics for dates ', start_date, ' to ', end_date
97 write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
98 write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
99
100 date = start_date
101 cdate = sdate
102 count = 0
103
104 !---------------------------------------------------------------------------------------------
105 write(6,'(a)')' [2] Read fields from standard files, and calculate mean fields'
106 !---------------------------------------------------------------------------------------------
107
108 do while ( cdate <= edate )
109 do member = 1, ne
110 count = count + 1
111 count_inv = 1.0 / real(count)
112
113 write(6,'(a,a)')' Processing data for date ', date
114
115 if ( be_method == 'NMC' ) then
116 filename = trim(dat_dir)//'/pert.'//date(1:10)//'.e001'
117 else
118 write(UNIT=ce,FMT='(i3.3)')member
119 filename = trim(dat_dir)//'/pert.'//date(1:10)//'.e'//trim(ce)
120 endif
121
122 open (iunit, file = trim(filename), form='unformatted')
123 read(iunit)date, ni, nj, nk
124
125 if ( count == 1 ) then
126 write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk
127 allocate( ps_prime(1:ni,1:nj) )
128 allocate( t_prime(1:ni,1:nj,1:nk) )
129 allocate( psi_prime(1:ni,1:nj,1:nk) )
130 allocate( chi_prime(1:ni,1:nj,1:nk) )
131 allocate( rh_prime(1:ni,1:nj,1:nk) )
132 allocate( height(1:ni,1:nj,1:nk) )
133 allocate( latitude(1:ni,1:nj) )
134 allocate( psi_mean(1:ni,1:nj,1:nk) )
135 allocate( chi_mean(1:ni,1:nj,1:nk) )
136 allocate( t_mean(1:ni,1:nj,1:nk) )
137 allocate( rh_mean(1:ni,1:nj,1:nk) )
138 allocate( ps_mean(1:ni,1:nj) )
139 allocate( bin(1:ni,1:nj,1:nk) )
140 allocate( bin2d(1:ni,1:nj) )
141 psi_mean(:,:,:) = 0.0
142 chi_mean(:,:,:) = 0.0
143 t_mean(:,:,:) = 0.0
144 rh_mean(:,:,:) = 0.0
145 ps_mean(:,:) = 0.0
146
147 end if
148
149 read(iunit)psi_prime
150 read(iunit)chi_prime
151 read(iunit)t_prime
152 read(iunit)rh_prime
153 read(iunit)ps_prime
154 read(iunit)height
155 read(iunit)latitude
156 close(iunit)
157
158
159 if ( count == 1 ) then
160 call da_create_bins( ni, nj, nk, bin_type, num_bins, num_bins2d, bin, bin2d, &
161 lat_min, lat_max, binwidth_lat, &
162 hgt_min, hgt_max, binwidth_hgt, latitude, height )
163 end if
164
165 !---------------------------------------------------------------------------------------------
166 ! write(6,(2a)) [2] Calculate time/ensemble mean.
167 !---------------------------------------------------------------------------------------------
168
169 psi_mean = ( real( count-1 ) * psi_mean + psi_prime ) * count_inv
170 chi_mean = ( real( count-1 ) * chi_mean + chi_prime ) * count_inv
171 t_mean = ( real( count-1 ) * t_mean + t_prime ) * count_inv
172 rh_mean = ( real( count-1 ) * rh_mean + rh_prime ) * count_inv
173 ps_mean = ( real( count-1 ) * ps_mean + ps_prime ) * count_inv
174
175 end do ! End loop over ensemble members.
176
177 ! Calculate next date:
178 call da_advance_cymdh( date, interval, new_date )
179 date = new_date
180 read(date(1:10), fmt='(i10)')cdate
181 end do ! End loop over times.
182
183 !---------------------------------------------------------------------------------------------
184 write(6,'(a)')' [2] Read fields again, and remove time/ensemble/area mean'
185 !---------------------------------------------------------------------------------------------
186
187 date = start_date
188 cdate = sdate
189
190 do while ( cdate <= edate )
191 do member = 1, ne
192
193 write(6,'(a,a)')' Removing mean for date ', date
194
195 if ( be_method == 'NMC' ) then
196 filename = trim(dat_dir)//'/pert.'//date(1:10)//'.e001'
197 else
198 write(UNIT=ce,FMT='(i3.3)')member
199 filename = trim(dat_dir)//'/pert.'//date(1:10)//'.e'//trim(ce)
200 endif
201
202 open (iunit, file = trim(filename), form='unformatted')
203 read(iunit)date, ni, nj, nk
204 read(iunit)psi_prime
205 read(iunit)chi_prime
206 read(iunit)t_prime
207 read(iunit)rh_prime
208 read(iunit)ps_prime
209 read(iunit)height
210 read(iunit)latitude
211 close(iunit)
212
213 !---------------------------------------------------------------------------------------------
214 ! write(6,(2a)) [2] Remove mean.
215 !---------------------------------------------------------------------------------------------
216
217 psi_prime = psi_prime - psi_mean
218 chi_prime = chi_prime - chi_mean
219 t_prime = t_prime - t_mean
220 rh_prime = rh_prime - rh_mean
221 ps_prime = ps_prime - ps_mean
222
223 !---------------------------------------------------------------------------------------------
224 ! write(6,(2a)) [2] Write fields to file for further processing.
225 !---------------------------------------------------------------------------------------------
226
227 write(ce,'(i3.3)')member
228
229 ! Write necessary full-fields:
230 variable = 'fullflds'
231 filename = trim(variable)//'/'//date(1:10)
232 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
233 open (ounit, file = filename, form='unformatted')
234 write(ounit)ni, nj, nk
235 write(ounit)latitude
236 write(ounit)height
237 close(ounit)
238
239 ! Write psi:
240 variable = 'psi'
241 filename = trim(variable)//'/'//date(1:10)
242 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
243 open (ounit, file = filename, form='unformatted')
244 write(ounit)ni, nj, nk
245 write(ounit)psi_prime
246 close(ounit)
247
248 ! Write chi:
249 variable = 'chi'
250 filename = trim(variable)//'/'//date(1:10)
251 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
252 open (ounit, file = filename, form='unformatted')
253 write(ounit)ni, nj, nk
254 write(ounit)chi_prime
255 close(ounit)
256
257 ! Write T:
258 variable = 't'
259 filename = trim(variable)//'/'//date(1:10)
260 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
261 open (ounit, file = filename, form='unformatted')
262 write(ounit)ni, nj, nk
263 write(ounit)t_prime
264 close(ounit)
265
266 ! Write RH:
267 variable = 'rh'
268 filename = trim(variable)//'/'//date(1:10)
269 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
270 open (ounit, file = filename, form='unformatted')
271 write(ounit)ni, nj, nk
272 write(ounit)rh_prime
273 close(ounit)
274
275 ! Write ps:
276 variable = 'ps' ! 2D field
277 filename = trim(variable)//'/'//date(1:10)
278 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
279 open (ounit, file = filename, form='unformatted')
280 write(ounit)ni, nj, 1
281 write(ounit)ps_prime
282 close(ounit)
283
284 end do ! End loop over ensemble members.
285
286 ! Calculate next date:
287 call da_advance_cymdh( date, interval, new_date )
288 date = new_date
289 read(date(1:10), fmt='(i10)')cdate
290 end do ! End loop over times.
291
292 ! Finally, write bin info:
293
294 filename = 'bin.data'
295 open (ounit, file = filename, form='unformatted')
296 write(ounit)bin_type
297 write(ounit)lat_min, lat_max, binwidth_lat
298 write(ounit)hgt_min, hgt_max, binwidth_hgt
299 write(ounit)num_bins, num_bins2d
300 write(ounit)bin(1:ni,1:nj,1:nk)
301 write(ounit)bin2d(1:ni,1:nj)
302 close(ounit)
303
304 end program gen_be_stage1