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