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