gen_be_etkf.f90
References to this file elsewhere.
1 program gen_be_etkf
2
3 !----------------------------------------------------------------------
4 ! Purpose : Perform an Ensemble Transform Kalman Filter (ETKF) rescaling
5 ! of WRF ensemble forecast data.
6 !
7 ! Owner: Dale Barker (NCAR/MMM) - WRF wrappper, Xuguang Wang (NOAA) - ETKF algorithm.
8 ! Please acknowledge author/institute in work that uses this code.
9 !
10 !----------------------------------------------------------------------
11
12 #ifdef crayx1
13 #define iargc ipxfargc
14 #endif
15
16 use da_control, only : trace_use, stdout,filename_len
17 ! use da_gen_be
18 use da_etkf, only : da_solve_etkf
19 use da_reporting, only : da_error, message
20
21 implicit none
22
23 #include "netcdf.inc"
24
25 integer, parameter :: max_num_vars = 50 ! Maximum number of variables.
26 integer, parameter :: max_num_dims = 20 ! Maximum number of dimensions.
27 integer, parameter :: unit = 100 ! Unit number.
28
29 character (len=filename_len) :: input_file ! Input file.
30 character (len=filename_len) :: output_file ! Output file.
31 character (len=3) :: ce ! Member index -> character.
32
33 integer :: num_members ! Ensemble size.
34 integer :: nv ! Number of variables.
35 integer :: num_obs ! Number of observations.
36 integer :: naccumt1 ! Number of previous cycles.
37 integer :: naccumt2 ! Number of previous cycles.
38 integer :: nstartaccum1 ! Cycle from which naccumt1 cycle starts.
39 integer :: nstartaccum2 ! Cycle from which naccumt2 cycle starts.
40 integer :: nout ! Output record for inn. vec./ob. error var.
41 integer :: length ! Filename length.
42 integer :: rcode ! NETCDF return code.
43 integer :: cdfid ! NETCDF file IDs.
44 integer :: member, v, o, i, j, k, ijkv ! Loop counters.
45 integer :: ivtype ! 4=integer, 5=real, 6=d.p.
46 integer :: natts ! Number of field attributes.
47
48 integer :: index ! Array index.
49 integer :: nijkv ! Array counter.
50 integer :: iend ! End of array
51 real :: num_members_inv ! 1 / ensemble size.
52 real :: tainflatinput ! Pre-specified inflation, if not using adaptive inflation.
53 real :: rhoinput ! Pre-specified inflation, if not using adaptive rho factor.
54 real :: ds ! Grid resolution (m).
55
56 character(len=10) :: cv(1:max_num_vars) ! Default array of variable names.
57 integer :: id_var(1:max_num_vars) ! NETCDF variable ID.
58 integer :: ndims(1:max_num_vars) ! Number of field dimensions.
59 integer :: istart(1:max_num_vars) ! Start index.
60 integer :: dimids(1:max_num_dims) ! Array of dimension IDs.
61 integer :: one(1:max_num_dims) ! Array of dimension starts.
62 integer :: dims(1:max_num_vars,1:max_num_dims) ! Array of dimensions.
63 integer :: dim_prod(1:max_num_dims) ! Product of array dimensions.
64
65 real (kind=4), allocatable :: data_r(:,:,:) ! Data array.
66
67 real, pointer :: xf(:,:) ! Ensemble perturbations.
68 real, pointer :: xf_mean(:) ! Ensemble perturbation mean.
69 real, pointer :: xf_vari(:) ! Ensemble perturbation variance.
70 real, pointer :: y(:,:) ! H(xf).
71 real, pointer :: sigma_o2(:) ! Ob error variance.
72 real, pointer :: yo(:) ! Observation.
73 real, pointer :: ens_mean(:) ! Variable ensemble mean.
74 real, pointer :: ens_stdv_pert_prior(:) ! Variable prior perturbation std. dev.
75 real, pointer :: ens_stdv_pert_poste(:) ! Variable posterior perturbation std. dev.
76
77 namelist / gen_be_etkf_nl / num_members, nv, cv, &
78 naccumt1, naccumt2, nstartaccum1, nstartaccum2, &
79 nout, tainflatinput, rhoinput
80
81 !---------------------------------------------------------------------------------------------
82 write(stdout,'(/a)')' [1] Initialize information.'
83 !-----------------------------------------------------------------------------------------
84
85 num_members = 56
86 nv = 1
87 cv = "U"
88
89 naccumt1 = 0
90 naccumt2 = 0
91 nstartaccum1 = 0
92 nstartaccum2 = 0
93 nout = 1
94 tainflatinput = 0.0
95 rhoinput = 0.0
96
97 open(unit=unit, file='gen_be_etkf_nl.nl', &
98 form='formatted', status='old', action='read')
99 read(unit, gen_be_etkf_nl)
100 close(unit)
101
102 write(stdout,'(a,i4)')' Number of ensemble members = ', num_members
103 write(stdout,'(a,i4)')' Number of prognostic variables = ', nv
104 write(stdout,'(50a)')' List of prognostic variables = ', cv(1:nv)
105 write(stdout,'(a,i4)')' naccumt1 = ', naccumt1
106 write(stdout,'(a,i4)')' naccumt2 = ', naccumt2
107 write(stdout,'(a,i4)')' nstartaccum1 = ', nstartaccum1
108 write(stdout,'(a,i4)')' nstartaccum2 = ', nstartaccum2
109 write(stdout,'(a,i4)')' nout = ', nout
110 write(stdout,'(a,f15.5)')' tainflatinput = ', tainflatinput
111 write(stdout,'(a,f15.5)')' rhoinput = ', rhoinput
112
113 num_members_inv = 1.0 / real(num_members)
114
115 allocate( ens_mean(1:nv) )
116 allocate( ens_stdv_pert_prior(1:nv) )
117 allocate( ens_stdv_pert_poste(1:nv) )
118
119 !-----------------------------------------------------------------------------------------
120 write(stdout,'(/a)')' [2] Read observation information.'
121 !-----------------------------------------------------------------------------------------
122
123 do member = 1, num_members
124
125 write(unit=ce,FMT='(i3.3)')member
126 input_file = 'ob.etkf.e'//ce
127 open(unit, file = input_file, status='old')
128 read(unit,*)num_obs
129
130 if ( member == 1 ) then
131 write(stdout,'(a,i10)')' Number of observations = ', num_obs
132 allocate( y(1:num_obs,1:num_members) )
133 allocate( sigma_o2(1:num_obs) )
134 allocate( yo(1:num_obs) )
135 end if
136
137 do o = 1, num_obs
138 read(unit,'(3f17.7)')yo(o), y(o,member), sigma_o2(o)
139 ! Convert yo-H(xb) to H(xb):
140 y(o,member) = yo(o) - y(o,member)
141 ! Convert ob error standard deviation to variance:
142 sigma_o2(o) = sigma_o2(o) * sigma_o2(o)
143 end do
144 end do
145
146 !-----------------------------------------------------------------------------------------
147 write(stdout,'(/a)')' [3] Set up arrays using input ensemble mean forecast'
148 !-----------------------------------------------------------------------------------------
149
150 ! Open mean:
151 input_file = 'etkf_input'
152 length = len_trim(input_file)
153 rcode = nf_open( input_file(1:length), NF_NOWRITE, cdfid )
154
155 do v = 1, nv
156
157 ! Check variable is in file, and get variable id:
158 rcode = nf_inq_varid ( cdfid, cv(v), id_var(v) )
159 if ( rcode /= 0 ) then
160 write(UNIT=message(1),FMT='(A,A)') &
161 cv(v), ' variable is not in input file'
162 call da_error(__FILE__,__LINE__,message(1:1))
163 end if
164
165 ! Get number of dimensions, and check of real type:
166 dimids = 0
167 rcode = nf_inq_var( cdfid, id_var(v), cv(v), ivtype, ndims(v), dimids, natts )
168 if ( ivtype /= 5 ) then
169 write(UNIT=message(1),FMT='(A,A)') cv(v), ' variable is not real type'
170 call da_error(__FILE__,__LINE__,message(1:1))
171 end if
172
173 ! Get dimensions of field:
174 dims(v,:) = 0
175 do i = 1, ndims(v)
176 rcode = nf_inq_dimlen( cdfid, dimids(i), dims(v,i) )
177 end do
178 end do
179
180 one = 1
181 istart(1) = 1
182 do v = 2, nv
183 istart(v) = istart(v-1) + product(dims(v-1,1:ndims(v-1)-1))
184 end do
185 nijkv = istart(nv) + product(dims(nv,1:ndims(nv)-1)) - 1
186
187 allocate( xf_mean(1:nijkv) )
188 do v = 1, nv
189 allocate( data_r(dims(v,1),dims(v,2),dims(v,3)))
190
191 call ncvgt( cdfid, id_var(v), one, dims(v,:), data_r, rcode)
192
193 ! Fill 4D array:
194 index = istart(v)
195 do k = 1, dims(v,3)
196 do j = 1, dims(v,2)
197 do i = 1, dims(v,1)
198 xf_mean(index) = data_r(i,j,k)
199 index = index + 1
200 end do
201 end do
202 end do
203 deallocate( data_r )
204 end do ! v
205
206 rcode = nf_close( cdfid )
207
208 !-----------------------------------------------------------------------------------------
209 write(stdout,'(/a)')' [4] Extract necessary fields from WRF ensemble forecasts.'
210 !-----------------------------------------------------------------------------------------
211
212 allocate( xf(1:nijkv,1:num_members) )
213 allocate( xf_vari(1:nijkv) )
214
215 do member = 1, num_members
216 write(UNIT=ce,FMT='(i3.3)')member
217
218 ! Open file:
219 input_file = 'etkf_input.e'//ce
220 length = len_trim(input_file)
221 rcode = nf_open( input_file(1:length), NF_NOWRITE, cdfid )
222
223 do v = 1, nv
224 allocate( data_r(dims(v,1),dims(v,2),dims(v,3)))
225
226 call ncvgt( cdfid, id_var(v), one, dims(v,:), data_r, rcode)
227
228 ! Fill 4D array:
229 index = istart(v)
230 do k = 1, dims(v,3)
231 do j = 1, dims(v,2)
232 do i = 1, dims(v,1)
233 xf(index,member) = data_r(i,j,k)
234 index = index + 1
235 end do
236 end do
237 end do
238 deallocate( data_r )
239 end do ! v
240 rcode = nf_close( cdfid )
241 end do !member
242
243 ! Convert ensemble forecasts to perturbations:
244 do ijkv = 1, nijkv
245 xf(ijkv,1:num_members) = xf(ijkv,1:num_members) - xf_mean(ijkv)
246 xf_vari(ijkv) = sum(xf(ijkv,1:num_members)**2) * num_members_inv
247 end do
248
249 ! Print prior mean, ensemble standard deviation:
250 do v = 1, nv
251 iend = istart(v) + product(dims(v,1:ndims(v)-1)) - 1
252 ens_mean(v) = sum(xf_mean(istart(v):iend)) / real(iend - istart(v) + 1)
253 ens_stdv_pert_prior(v) =sqrt( sum(xf_vari(istart(v):iend)) / &
254 real(iend - istart(v) + 1) )
255 end do
256
257 !-----------------------------------------------------------------------------------------
258 write(stdout,'(/a)')' [5] Call ETKF:'
259 !-----------------------------------------------------------------------------------------
260
261 call da_solve_etkf( nijkv, num_members, num_obs, xf, y, sigma_o2, yo, nout, &
262 naccumt1, naccumt2, nstartaccum1, nstartaccum2, tainflatinput, &
263 rhoinput )
264
265 ! Calculate posterior ensemble standard deviation and output:
266 do ijkv = 1, nijkv
267 xf_vari(ijkv) = sum(xf(ijkv,1:num_members)**2) * num_members_inv
268 end do
269
270 ! Print posterior mean, ensemble standard deviation:
271 write(stdout,'(5a)')' v', ' Variable ', ' Ensemble Mean', &
272 ' Prior Pert StDv', ' Post. Pert. StDv'
273 do v = 1, nv
274 iend = istart(v) + product(dims(v,1:ndims(v)-1)) - 1
275 ens_stdv_pert_poste(v) =sqrt( sum(xf_vari(istart(v):iend)) / &
276 real(iend - istart(v) + 1) )
277
278 write(stdout,'(i4,1x,a10,3f17.7)')v, cv(v), &
279 ens_mean(v), ens_stdv_pert_prior(v), ens_stdv_pert_poste(v)
280 end do
281
282 !-----------------------------------------------------------------------------------------
283 write(stdout,'(/a)')' [6] Output ETKF analysis ensemble:'
284 !-----------------------------------------------------------------------------------------
285
286 do member = 1, num_members
287 write(UNIT=ce,FMT='(i3.3)')member
288
289 ! Open file:
290 input_file = 'etkf_output.e'//ce
291 length = len_trim(input_file)
292 rcode = nf_open( input_file(1:length), NF_WRITE, cdfid )
293 if ( rcode /= 0 ) then
294 write(UNIT=message(1),FMT='(A,A)') &
295 ' Error opening netcdf file ', input_file(1:length)
296 call da_error(__FILE__,__LINE__,message(1:1))
297 end if
298
299 do v = 1, nv
300 allocate( data_r(dims(v,1),dims(v,2),dims(v,3)))
301
302 ! Add updated perturbations back to ensemble mean and output
303 index = istart(v)
304 do k = 1, dims(v,3)
305 do j = 1, dims(v,2)
306 do i = 1, dims(v,1)
307 data_r(i,j,k) = xf_mean(index) + xf(index,member)
308 index = index + 1
309 end do
310 end do
311 end do
312 call ncvpt( cdfid, id_var(v), one, dims(v,:), data_r, rcode)
313 deallocate( data_r )
314 end do ! v
315 rcode = nf_close( cdfid )
316 end do !member
317
318 deallocate( ens_mean )
319 deallocate( ens_stdv_pert_prior )
320 deallocate( ens_stdv_pert_poste )
321
322 end program gen_be_etkf
323