gen_be_ep2.f90
References to this file elsewhere.
1 program gen_be_ep2
2 !
3 !----------------------------------------------------------------------
4 ! Purpose : To convert WRF ensemble to format required for use as
5 ! flow-dependent perturbations in WRF-Var (alpha control variable,
6 ! alphacv_method = 2).
7 !
8 ! Owner: Dale Barker (NCAR/MMM)
9 ! Please acknowledge author/institute in work that uses this code.
10 !
11 !----------------------------------------------------------------------
12
13 #ifdef crayx1
14 #define iargc ipxfargc
15 #endif
16
17 use da_control, only : stderr, stdout, filename_len
18 use da_tools_serial, only : da_get_unit, da_free_unit
19 use da_gen_be, only : da_stage0_initialize, da_get_field, da_get_trh
20
21 implicit none
22
23 character (len=filename_len) :: filestub ! General filename stub.
24 character (len=filename_len) :: input_file ! Input file.
25 character (len=filename_len) :: output_file ! Output file.
26 character (len=10) :: date ! Character date.
27 character (len=10) :: var ! Variable to search for.
28 character (len=3) :: cne ! Ensemble size.
29 character (len=3) :: ce ! Member index -> character.
30
31 integer, external :: iargc
32 integer :: numarg
33 integer :: ne ! Ensemble size.
34 integer :: i, j, k, member ! Loop counters.
35 integer :: dim1 ! Dimensions of grid (T points).
36 integer :: dim1s ! Dimensions of grid (vor/psi pts).
37 integer :: dim2 ! Dimensions of grid (T points).
38 integer :: dim2s ! Dimensions of grid (vor/psi pts).
39 integer :: dim3 ! Dimensions of grid (T points).
40 real :: member_inv ! 1 / member.
41 real :: ds ! Grid resolution.
42 logical :: remove_mean ! Remove mean from standard fields.
43
44 real, allocatable :: u(:,:,:) ! u-wind.
45 real, allocatable :: v(:,:,:) ! v-wind.
46 real, allocatable :: temp(:,:,:) ! Temperature.
47 real, allocatable :: q(:,:,:) ! Specific humidity.
48 real, allocatable :: psfc(:,:) ! Surface pressure.
49 real, allocatable :: u_mean(:,:,:) ! u-wind.
50 real, allocatable :: v_mean(:,:,:) ! v-wind.
51 real, allocatable :: temp_mean(:,:,:) ! Temperature.
52 real, allocatable :: q_mean(:,:,:) ! Specific humidity.
53 real, allocatable :: psfc_mean(:,:) ! Surface pressure.
54 real, allocatable :: u_mnsq(:,:,:) ! u-wind.
55 real, allocatable :: v_mnsq(:,:,:) ! v-wind.
56 real, allocatable :: temp_mnsq(:,:,:) ! Temperature.
57 real, allocatable :: q_mnsq(:,:,:) ! Specific humidity.
58 real, allocatable :: psfc_mnsq(:,:) ! Surface pressure.
59
60 real, allocatable :: utmp(:,:) ! u-wind.
61 real, allocatable :: vtmp(:,:) ! v-wind.
62 real, allocatable :: ttmp(:,:) ! temperature.
63 real, allocatable :: dummy(:,:) ! dummy.
64
65 integer :: gen_be_iunit, gen_be_ounit
66
67 stderr = 0
68 stdout = 6
69
70 !---------------------------------------------------------------------------------------------
71 write(6,'(/a)')' [1] Initialize information.'
72 !---------------------------------------------------------------------------------------------
73
74 call da_get_unit(gen_be_iunit)
75 call da_get_unit(gen_be_ounit)
76
77 remove_mean = .true.
78
79 numarg = iargc()
80 if ( numarg /= 3 )then
81 write(UNIT=6,FMT='(a)') &
82 "Usage: gen_be_ep2 date ne <wrf_file_stub> Stop"
83 stop
84 end if
85
86 ! Initialse to stop Cray compiler complaining
87 date=""
88 cne=""
89 filestub=""
90
91 call getarg( 1, date )
92 call getarg( 2, cne )
93 read(cne,'(i3)')ne
94 call getarg( 3, filestub )
95
96 if ( remove_mean ) then
97 write(6,'(a,a)')' Computing gen_be ensemble perturbation files for date ', date
98 else
99 write(6,'(a,a)')' Computing gen_be ensemble forecast files for date ', date
100 end if
101 write(6,'(a)')' Perturbations are in MODEL SPACE (u, v, t, q, ps)'
102 write(6,'(a,i4)')' Ensemble Size = ', ne
103 write(6,'(a,a)')' Input filestub = ', trim(filestub)
104
105 !---------------------------------------------------------------------------------------------
106 write(6,'(/a)')' [2] Set up data dimensions and allocate arrays:'
107 !---------------------------------------------------------------------------------------------
108
109 ! Get grid dimensions from first T field:
110 var = "T"
111 input_file = trim(filestub)//'.e001'
112 call da_stage0_initialize( input_file, var, dim1, dim2, dim3, ds )
113 dim1s = dim1+1 ! u i dimension is 1 larger.
114 dim2s = dim2+1 ! v j dimension is 1 larger.
115
116 ! Allocate arrays in output fields:
117 allocate( u(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to mass pts for output.
118 allocate( v(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to mass pts for output.
119 allocate( temp(1:dim1,1:dim2,1:dim3) )
120 allocate( q(1:dim1,1:dim2,1:dim3) )
121 allocate( psfc(1:dim1,1:dim2) )
122 allocate( u_mean(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to chi pts for output.
123 allocate( v_mean(1:dim1,1:dim2,1:dim3) )
124 allocate( temp_mean(1:dim1,1:dim2,1:dim3) )
125 allocate( q_mean(1:dim1,1:dim2,1:dim3) )
126 allocate( psfc_mean(1:dim1,1:dim2) )
127 allocate( u_mnsq(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to chi pts for output.
128 allocate( v_mnsq(1:dim1,1:dim2,1:dim3) )
129 allocate( temp_mnsq(1:dim1,1:dim2,1:dim3) )
130 allocate( q_mnsq(1:dim1,1:dim2,1:dim3) )
131 allocate( psfc_mnsq(1:dim1,1:dim2) )
132 u_mean = 0.0
133 v_mean = 0.0
134 temp_mean = 0.0
135 q_mean = 0.0
136 psfc_mean = 0.0
137 u_mnsq = 0.0
138 v_mnsq = 0.0
139 temp_mnsq = 0.0
140 q_mnsq = 0.0
141 psfc_mnsq = 0.0
142
143 ! Temporary arrays:
144 allocate( utmp(1:dim1s,1:dim2) ) ! u on Arakawa C-grid.
145 allocate( vtmp(1:dim1,1:dim2s) ) ! v on Arakawa C-grid.
146 allocate( ttmp(1:dim1,1:dim2) )
147 allocate( dummy(1:dim1,1:dim2) )
148
149 !---------------------------------------------------------------------------------------------
150 write(6,'(/a)')' [3] Extract necessary fields from input NETCDF files and output.'
151 !---------------------------------------------------------------------------------------------
152
153 do member = 1, ne
154
155 write(UNIT=ce,FMT='(i3.3)')member
156 input_file = trim(filestub)//'.e'//ce
157
158 do k = 1, dim3
159
160 ! Read u, v:
161 var = "U"
162 call da_get_field( input_file, var, 3, dim1s, dim2, dim3, k, utmp )
163 var = "V"
164 call da_get_field( input_file, var, 3, dim1, dim2s, dim3, k, vtmp )
165
166 ! Interpolate u to mass pts:
167 do j = 1, dim2
168 do i = 1, dim1
169 u(i,j,k) = 0.5 * ( utmp(i,j) + utmp(i+1,j) )
170 v(i,j,k) = 0.5 * ( vtmp(i,j) + vtmp(i,j+1) )
171 end do
172 end do
173
174 ! Read theta, and convert to temperature:
175 call da_get_trh( input_file, dim1, dim2, dim3, k, ttmp, dummy )
176 temp(:,:,k) = ttmp(:,:)
177
178 ! Read mixing ratio, and convert to specific humidity:
179 var = "QVAPOR"
180 call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, dummy )
181 q(:,:,k) = dummy(:,:) / ( 1.0 + dummy(:,:) )
182
183 end do
184
185 ! Finally, extract surface pressure:
186 var = "PSFC"
187 call da_get_field( input_file, var, 2, dim1, dim2, dim3, 1, psfc )
188
189 ! Write out ensemble forecasts for this member:
190 output_file = 'tmp.e'//ce
191 open (gen_be_ounit, file = output_file, form='unformatted')
192 write(gen_be_ounit)date, dim1, dim2, dim3
193 write(gen_be_ounit)u
194 write(gen_be_ounit)v
195 write(gen_be_ounit)temp
196 write(gen_be_ounit)q
197 write(gen_be_ounit)psfc
198 close(gen_be_ounit)
199
200 ! Calculate accumulating mean and mean square:
201 member_inv = 1.0 / real(member)
202 u_mean = ( real( member-1 ) * u_mean + u ) * member_inv
203 v_mean = ( real( member-1 ) * v_mean + v ) * member_inv
204 temp_mean = ( real( member-1 ) * temp_mean + temp ) * member_inv
205 q_mean = ( real( member-1 ) * q_mean + q ) * member_inv
206 psfc_mean = ( real( member-1 ) * psfc_mean + psfc ) * member_inv
207 u_mnsq = ( real( member-1 ) * u_mnsq + u * u ) * member_inv
208 v_mnsq = ( real( member-1 ) * v_mnsq + v * v ) * member_inv
209 temp_mnsq = ( real( member-1 ) * temp_mnsq + temp * temp ) * member_inv
210 q_mnsq = ( real( member-1 ) * q_mnsq + q * q ) * member_inv
211 psfc_mnsq = ( real( member-1 ) * psfc_mnsq + psfc * psfc ) * member_inv
212
213 end do
214
215 deallocate( utmp )
216 deallocate( vtmp )
217 deallocate( ttmp )
218 deallocate( dummy )
219
220 !---------------------------------------------------------------------------------------------
221 write(6,'(/a)')' [4] Compute perturbations and output'
222 !---------------------------------------------------------------------------------------------
223
224 if ( remove_mean ) then
225 write(6,'(a)') " Calculate ensemble perturbations"
226 else
227 write(6,'(a)') " WARNING: Not removing ensemble mean (outputs are full-fields)"
228 end if
229
230 do member = 1, ne
231 write(UNIT=ce,FMT='(i3.3)')member
232
233 ! Re-read ensemble member standard fields:
234 input_file = 'tmp.e'//ce
235 open (gen_be_iunit, file = input_file, form='unformatted')
236 read(gen_be_iunit)date, dim1, dim2, dim3
237 read(gen_be_iunit)u
238 read(gen_be_iunit)v
239 read(gen_be_iunit)temp
240 read(gen_be_iunit)q
241 read(gen_be_iunit)psfc
242 close(gen_be_iunit)
243
244 if ( remove_mean ) then
245 u = u - u_mean
246 v = v - v_mean
247 temp = temp - temp_mean
248 q = q - q_mean
249 psfc = psfc - psfc_mean
250 end if
251
252 ! Write out perturbations for this member:
253
254 output_file = 'u/'//date(1:10)//'.u.e'//trim(ce) ! Output u.
255 open (gen_be_ounit, file = output_file, form='unformatted')
256 write(gen_be_ounit)dim1, dim2, dim3
257 write(gen_be_ounit)u
258 close(gen_be_ounit)
259
260 output_file = 'v/'//date(1:10)//'.v.e'//trim(ce) ! Output v.
261 open (gen_be_ounit, file = output_file, form='unformatted')
262 write(gen_be_ounit)dim1, dim2, dim3
263 write(gen_be_ounit)v
264 close(gen_be_ounit)
265
266 output_file = 't/'//date(1:10)//'.t.e'//trim(ce) ! Output t.
267 open (gen_be_ounit, file = output_file, form='unformatted')
268 write(gen_be_ounit)dim1, dim2, dim3
269 write(gen_be_ounit)temp
270 close(gen_be_ounit)
271
272 output_file = 'q/'//date(1:10)//'.q.e'//trim(ce) ! Output q.
273 open (gen_be_ounit, file = output_file, form='unformatted')
274 write(gen_be_ounit)dim1, dim2, dim3
275 write(gen_be_ounit)q
276 close(gen_be_ounit)
277
278 output_file = 'ps/'//date(1:10)//'.ps.e'//trim(ce) ! Output ps.
279 open (gen_be_ounit, file = output_file, form='unformatted')
280 write(gen_be_ounit)dim1, dim2, dim3
281 write(gen_be_ounit).false., .false.
282 write(gen_be_ounit)psfc
283 close(gen_be_ounit)
284
285 end do
286
287 ! Write out mean/stdv fields (stdv stored in mnsq arrays):
288 u_mnsq = sqrt( u_mnsq - u_mean * u_mean )
289 v_mnsq = sqrt( v_mnsq - v_mean * v_mean )
290 temp_mnsq = sqrt( temp_mnsq - temp_mean * temp_mean )
291 q_mnsq = sqrt( q_mnsq - q_mean * q_mean )
292 psfc_mnsq = sqrt( psfc_mnsq - psfc_mean * psfc_mean )
293
294 output_file = 'u/'//date(1:10)//'.u.mean' ! Output u.
295 open (gen_be_ounit, file = output_file, form='unformatted')
296 write(gen_be_ounit)dim1, dim2, dim3
297 write(gen_be_ounit)u_mean
298 close(gen_be_ounit)
299
300 output_file = 'u/'//date(1:10)//'.u.stdv' ! Output u.
301 open (gen_be_ounit, file = output_file, form='unformatted')
302 write(gen_be_ounit)dim1, dim2, dim3
303 write(gen_be_ounit)u_mnsq
304 close(gen_be_ounit)
305
306 output_file = 'v/'//date(1:10)//'.v.mean' ! Output v.
307 open (gen_be_ounit, file = output_file, form='unformatted')
308 write(gen_be_ounit)dim1, dim2, dim3
309 write(gen_be_ounit)v_mean
310 close(gen_be_ounit)
311
312 output_file = 'v/'//date(1:10)//'.v.stdv' ! Output v.
313 open (gen_be_ounit, file = output_file, form='unformatted')
314 write(gen_be_ounit)dim1, dim2, dim3
315 write(gen_be_ounit)v_mnsq
316 close(gen_be_ounit)
317
318 output_file = 't/'//date(1:10)//'.t.mean' ! Output t.
319 open (gen_be_ounit, file = output_file, form='unformatted')
320 write(gen_be_ounit)dim1, dim2, dim3
321 write(gen_be_ounit)temp_mean
322 close(gen_be_ounit)
323
324 output_file = 't/'//date(1:10)//'.t.stdv' ! Output t.
325 open (gen_be_ounit, file = output_file, form='unformatted')
326 write(gen_be_ounit)dim1, dim2, dim3
327 write(gen_be_ounit)temp_mnsq
328 close(gen_be_ounit)
329
330 output_file = 'q/'//date(1:10)//'.q.mean' ! Output q.
331 open (gen_be_ounit, file = output_file, form='unformatted')
332 write(gen_be_ounit)dim1, dim2, dim3
333 write(gen_be_ounit)q_mean
334 close(gen_be_ounit)
335
336 output_file = 'q/'//date(1:10)//'.q.stdv' ! Output q.
337 open (gen_be_ounit, file = output_file, form='unformatted')
338 write(gen_be_ounit)dim1, dim2, dim3
339 write(gen_be_ounit)q_mnsq
340 close(gen_be_ounit)
341
342 output_file = 'ps/'//date(1:10)//'.ps.mean' ! Output ps.
343 open (gen_be_ounit, file = output_file, form='unformatted')
344 write(gen_be_ounit)dim1, dim2, dim3
345 write(gen_be_ounit)psfc_mean
346 close(gen_be_ounit)
347
348 output_file = 'ps/'//date(1:10)//'.ps.stdv' ! Output ps.
349 open (gen_be_ounit, file = output_file, form='unformatted')
350 write(gen_be_ounit)dim1, dim2, dim3
351 write(gen_be_ounit)psfc_mnsq
352 close(gen_be_ounit)
353
354 call da_free_unit(gen_be_iunit)
355 call da_free_unit(gen_be_ounit)
356
357 #ifdef crayx1
358 contains
359
360 subroutine getarg(i, harg)
361 implicit none
362 character(len=*) :: harg
363 integer :: ierr, ilen, i
364
365 call pxfgetarg(i, harg, ilen, ierr)
366 return
367 end subroutine getarg
368 #endif
369
370 end program gen_be_ep2
371