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
18 use da_gen_be
19 use da_tracing
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 if (trace_use) call da_trace_init
75 if (trace_use) call da_trace_entry("gen_be_ep2")
76
77 call da_get_unit(gen_be_iunit)
78 call da_get_unit(gen_be_ounit)
79
80 remove_mean = .true.
81
82 numarg = iargc()
83 if ( numarg /= 3 )then
84 write(UNIT=6,FMT='(a)') &
85 "Usage: gen_be_ep2 date ne <wrf_file_stub> Stop"
86 stop
87 end if
88
89 ! Initialse to stop Cray compiler complaining
90 date=""
91 cne=""
92 filestub=""
93
94 call getarg( 1, date )
95 call getarg( 2, cne )
96 read(cne,'(i3)')ne
97 call getarg( 3, filestub )
98
99 if ( remove_mean ) then
100 write(6,'(a,a)')' Computing gen_be ensemble perturbation files for date ', date
101 else
102 write(6,'(a,a)')' Computing gen_be ensemble forecast files for date ', date
103 end if
104 write(6,'(a)')' Perturbations are in MODEL SPACE (u, v, t, q, ps)'
105 write(6,'(a,i4)')' Ensemble Size = ', ne
106 write(6,'(a,a)')' Input filestub = ', trim(filestub)
107
108 !---------------------------------------------------------------------------------------------
109 write(6,'(/a)')' [2] Set up data dimensions and allocate arrays:'
110 !---------------------------------------------------------------------------------------------
111
112 ! Get grid dimensions from first T field:
113 var = "T"
114 input_file = trim(filestub)//'.e001'
115 call da_stage0_initialize( input_file, var, dim1, dim2, dim3, ds )
116 dim1s = dim1+1 ! u i dimension is 1 larger.
117 dim2s = dim2+1 ! v j dimension is 1 larger.
118
119 ! Allocate arrays in output fields:
120 allocate( u(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to mass pts for output.
121 allocate( v(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to mass pts for output.
122 allocate( temp(1:dim1,1:dim2,1:dim3) )
123 allocate( q(1:dim1,1:dim2,1:dim3) )
124 allocate( psfc(1:dim1,1:dim2) )
125 allocate( u_mean(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to chi pts for output.
126 allocate( v_mean(1:dim1,1:dim2,1:dim3) )
127 allocate( temp_mean(1:dim1,1:dim2,1:dim3) )
128 allocate( q_mean(1:dim1,1:dim2,1:dim3) )
129 allocate( psfc_mean(1:dim1,1:dim2) )
130 allocate( u_mnsq(1:dim1,1:dim2,1:dim3) ) ! Note - interpolated to chi pts for output.
131 allocate( v_mnsq(1:dim1,1:dim2,1:dim3) )
132 allocate( temp_mnsq(1:dim1,1:dim2,1:dim3) )
133 allocate( q_mnsq(1:dim1,1:dim2,1:dim3) )
134 allocate( psfc_mnsq(1:dim1,1:dim2) )
135 u_mean = 0.0
136 v_mean = 0.0
137 temp_mean = 0.0
138 q_mean = 0.0
139 psfc_mean = 0.0
140 u_mnsq = 0.0
141 v_mnsq = 0.0
142 temp_mnsq = 0.0
143 q_mnsq = 0.0
144 psfc_mnsq = 0.0
145
146 ! Temporary arrays:
147 allocate( utmp(1:dim1s,1:dim2) ) ! u on Arakawa C-grid.
148 allocate( vtmp(1:dim1,1:dim2s) ) ! v on Arakawa C-grid.
149 allocate( ttmp(1:dim1,1:dim2) )
150 allocate( dummy(1:dim1,1:dim2) )
151
152 !---------------------------------------------------------------------------------------------
153 write(6,'(/a)')' [3] Extract necessary fields from input NETCDF files and output.'
154 !---------------------------------------------------------------------------------------------
155
156 do member = 1, ne
157
158 write(UNIT=ce,FMT='(i3.3)')member
159 input_file = trim(filestub)//'.e'//ce
160
161 do k = 1, dim3
162
163 ! Read u, v:
164 var = "U"
165 call da_get_field( input_file, var, 3, dim1s, dim2, dim3, k, utmp )
166 var = "V"
167 call da_get_field( input_file, var, 3, dim1, dim2s, dim3, k, vtmp )
168
169 ! Interpolate u to mass pts:
170 do j = 1, dim2
171 do i = 1, dim1
172 u(i,j,k) = 0.5 * ( utmp(i,j) + utmp(i+1,j) )
173 v(i,j,k) = 0.5 * ( vtmp(i,j) + vtmp(i,j+1) )
174 end do
175 end do
176
177 ! Read theta, and convert to temperature:
178 call da_get_trh( input_file, dim1, dim2, dim3, k, ttmp, dummy )
179 temp(:,:,k) = ttmp(:,:)
180
181 ! Read mixing ratio, and convert to specific humidity:
182 var = "QVAPOR"
183 call da_get_field( input_file, var, 3, dim1, dim2, dim3, k, dummy )
184 q(:,:,k) = dummy(:,:) / ( 1.0 + dummy(:,:) )
185
186 end do
187
188 ! Finally, extract surface pressure:
189 var = "PSFC"
190 call da_get_field( input_file, var, 2, dim1, dim2, dim3, 1, psfc )
191
192 ! Write out ensemble forecasts for this member:
193 output_file = 'tmp.e'//ce
194 open (gen_be_ounit, file = output_file, form='unformatted')
195 write(gen_be_ounit)date, dim1, dim2, dim3
196 write(gen_be_ounit)u
197 write(gen_be_ounit)v
198 write(gen_be_ounit)temp
199 write(gen_be_ounit)q
200 write(gen_be_ounit)psfc
201 close(gen_be_ounit)
202
203 ! Calculate accumulating mean and mean square:
204 member_inv = 1.0 / real(member)
205 u_mean = ( real( member-1 ) * u_mean + u ) * member_inv
206 v_mean = ( real( member-1 ) * v_mean + v ) * member_inv
207 temp_mean = ( real( member-1 ) * temp_mean + temp ) * member_inv
208 q_mean = ( real( member-1 ) * q_mean + q ) * member_inv
209 psfc_mean = ( real( member-1 ) * psfc_mean + psfc ) * member_inv
210 u_mnsq = ( real( member-1 ) * u_mnsq + u * u ) * member_inv
211 v_mnsq = ( real( member-1 ) * v_mnsq + v * v ) * member_inv
212 temp_mnsq = ( real( member-1 ) * temp_mnsq + temp * temp ) * member_inv
213 q_mnsq = ( real( member-1 ) * q_mnsq + q * q ) * member_inv
214 psfc_mnsq = ( real( member-1 ) * psfc_mnsq + psfc * psfc ) * member_inv
215
216 end do
217
218 deallocate( utmp )
219 deallocate( vtmp )
220 deallocate( ttmp )
221 deallocate( dummy )
222
223 !---------------------------------------------------------------------------------------------
224 write(6,'(/a)')' [4] Compute perturbations and output'
225 !---------------------------------------------------------------------------------------------
226
227 if ( remove_mean ) then
228 write(6,'(a)') " Calculate ensemble perturbations"
229 else
230 write(6,'(a)') " WARNING: Not removing ensemble mean (outputs are full-fields)"
231 end if
232
233 do member = 1, ne
234 write(UNIT=ce,FMT='(i3.3)')member
235
236 ! Re-read ensemble member standard fields:
237 input_file = 'tmp.e'//ce
238 open (gen_be_iunit, file = input_file, form='unformatted')
239 read(gen_be_iunit)date, dim1, dim2, dim3
240 read(gen_be_iunit)u
241 read(gen_be_iunit)v
242 read(gen_be_iunit)temp
243 read(gen_be_iunit)q
244 read(gen_be_iunit)psfc
245 close(gen_be_iunit)
246
247 if ( remove_mean ) then
248 u = u - u_mean
249 v = v - v_mean
250 temp = temp - temp_mean
251 q = q - q_mean
252 psfc = psfc - psfc_mean
253 end if
254
255 ! Write out perturbations for this member:
256
257 output_file = 'u/'//date(1:10)//'.u.e'//trim(ce) ! Output u.
258 open (gen_be_ounit, file = output_file, form='unformatted')
259 write(gen_be_ounit)dim1, dim2, dim3
260 write(gen_be_ounit)u
261 close(gen_be_ounit)
262
263 output_file = 'v/'//date(1:10)//'.v.e'//trim(ce) ! Output v.
264 open (gen_be_ounit, file = output_file, form='unformatted')
265 write(gen_be_ounit)dim1, dim2, dim3
266 write(gen_be_ounit)v
267 close(gen_be_ounit)
268
269 output_file = 't/'//date(1:10)//'.t.e'//trim(ce) ! Output t.
270 open (gen_be_ounit, file = output_file, form='unformatted')
271 write(gen_be_ounit)dim1, dim2, dim3
272 write(gen_be_ounit)temp
273 close(gen_be_ounit)
274
275 output_file = 'q/'//date(1:10)//'.q.e'//trim(ce) ! Output q.
276 open (gen_be_ounit, file = output_file, form='unformatted')
277 write(gen_be_ounit)dim1, dim2, dim3
278 write(gen_be_ounit)q
279 close(gen_be_ounit)
280
281 output_file = 'ps/'//date(1:10)//'.ps.e'//trim(ce) ! Output ps.
282 open (gen_be_ounit, file = output_file, form='unformatted')
283 write(gen_be_ounit)dim1, dim2, dim3
284 write(gen_be_ounit).false., .false.
285 write(gen_be_ounit)psfc
286 close(gen_be_ounit)
287
288 end do
289
290 ! Write out mean/stdv fields (stdv stored in mnsq arrays):
291 u_mnsq = sqrt( u_mnsq - u_mean * u_mean )
292 v_mnsq = sqrt( v_mnsq - v_mean * v_mean )
293 temp_mnsq = sqrt( temp_mnsq - temp_mean * temp_mean )
294 q_mnsq = sqrt( q_mnsq - q_mean * q_mean )
295 psfc_mnsq = sqrt( psfc_mnsq - psfc_mean * psfc_mean )
296
297 output_file = 'u/'//date(1:10)//'.u.mean' ! Output u.
298 open (gen_be_ounit, file = output_file, form='unformatted')
299 write(gen_be_ounit)dim1, dim2, dim3
300 write(gen_be_ounit)u_mean
301 close(gen_be_ounit)
302
303 output_file = 'u/'//date(1:10)//'.u.stdv' ! Output u.
304 open (gen_be_ounit, file = output_file, form='unformatted')
305 write(gen_be_ounit)dim1, dim2, dim3
306 write(gen_be_ounit)u_mnsq
307 close(gen_be_ounit)
308
309 output_file = 'v/'//date(1:10)//'.v.mean' ! Output v.
310 open (gen_be_ounit, file = output_file, form='unformatted')
311 write(gen_be_ounit)dim1, dim2, dim3
312 write(gen_be_ounit)v_mean
313 close(gen_be_ounit)
314
315 output_file = 'v/'//date(1:10)//'.v.stdv' ! Output v.
316 open (gen_be_ounit, file = output_file, form='unformatted')
317 write(gen_be_ounit)dim1, dim2, dim3
318 write(gen_be_ounit)v_mnsq
319 close(gen_be_ounit)
320
321 output_file = 't/'//date(1:10)//'.t.mean' ! Output t.
322 open (gen_be_ounit, file = output_file, form='unformatted')
323 write(gen_be_ounit)dim1, dim2, dim3
324 write(gen_be_ounit)temp_mean
325 close(gen_be_ounit)
326
327 output_file = 't/'//date(1:10)//'.t.stdv' ! Output t.
328 open (gen_be_ounit, file = output_file, form='unformatted')
329 write(gen_be_ounit)dim1, dim2, dim3
330 write(gen_be_ounit)temp_mnsq
331 close(gen_be_ounit)
332
333 output_file = 'q/'//date(1:10)//'.q.mean' ! Output q.
334 open (gen_be_ounit, file = output_file, form='unformatted')
335 write(gen_be_ounit)dim1, dim2, dim3
336 write(gen_be_ounit)q_mean
337 close(gen_be_ounit)
338
339 output_file = 'q/'//date(1:10)//'.q.stdv' ! Output q.
340 open (gen_be_ounit, file = output_file, form='unformatted')
341 write(gen_be_ounit)dim1, dim2, dim3
342 write(gen_be_ounit)q_mnsq
343 close(gen_be_ounit)
344
345 output_file = 'ps/'//date(1:10)//'.ps.mean' ! Output ps.
346 open (gen_be_ounit, file = output_file, form='unformatted')
347 write(gen_be_ounit)dim1, dim2, dim3
348 write(gen_be_ounit)psfc_mean
349 close(gen_be_ounit)
350
351 output_file = 'ps/'//date(1:10)//'.ps.stdv' ! Output ps.
352 open (gen_be_ounit, file = output_file, form='unformatted')
353 write(gen_be_ounit)dim1, dim2, dim3
354 write(gen_be_ounit)psfc_mnsq
355 close(gen_be_ounit)
356
357 call da_free_unit(gen_be_iunit)
358 call da_free_unit(gen_be_ounit)
359
360 if (trace_use) call da_trace_exit("gen_be_ep2")
361 if (trace_use) call da_trace_report
362
363 #ifdef crayx1
364 contains
365
366 subroutine getarg(i, harg)
367 implicit none
368 character(len=*) :: harg
369 integer :: ierr, ilen, i
370
371 call pxfgetarg(i, harg, ilen, ierr)
372 return
373 end subroutine getarg
374 #endif
375
376 end program gen_be_ep2
377