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