gen_be_ensrf.f90
References to this file elsewhere.
1 program gen_be_ensrf
2 !
3 !----------------------------------------------------------------------
4 ! Purpose : Perform an Ensemble Square Root Filter (EnSRF) assimilation
5 ! using WRF forecast data.
6 !
7 ! Owner: Dale Barker (NCAR/MMM)
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,stderr,filename_len
17 use da_gen_be, only : da_stage0_initialize, da_get_field, da_get_trh
18 use da_tracing, only : da_trace_init,da_trace_entry,da_trace_exit, &
19 da_trace_report
20
21 implicit none
22
23 type info_type
24 character*40 :: id
25 character*40 :: platform
26 integer :: plat
27 integer :: variable
28 end type info_type
29
30 type loc_type
31 real :: x
32 real :: y
33 real :: z
34 end type loc_type
35
36 type ob_type
37 type (info_type):: info
38 type (loc_type) :: loc
39 real :: yo
40 real :: omb_mean
41 real :: oma_mean
42 real :: omb_rms
43 real :: oma_rms
44 real :: sigma_o
45 real :: sigma_b
46 end type ob_type
47
48 integer, parameter :: nv3d = 4 ! #3D variables (u, v, T, q).
49 integer, parameter :: nv2d = 1 ! #2D variables (ps).
50 integer, parameter :: unit = 100 ! Unit number.
51
52 character (len=filename_len) :: filestub ! General filename stub.
53 character (len=filename_len) :: input_file ! Input file.
54 character (len=filename_len) :: output_file ! Output file.
55 character (len=10) :: var ! Variable to search for.
56 character (len=3) :: ce ! Member index -> character.
57
58 integer :: num_members ! Ensemble size.
59 integer :: num_obs ! Number of observations.
60 integer :: ni ! 1st dimension size.
61 integer :: niu ! 1st dimension size (u)
62 integer :: nj ! 2nd dimension size.
63 integer :: njv ! 2nd dimension size (v)
64 integer :: nk ! 3rd dimension size.
65 integer :: nv ! Total number of 2D/3D variables.
66 integer :: nij, nijk, nijkv ! Dimensions.
67 integer :: o, member, i, j, k, v ! Loop counters.
68 integer :: imin, imax, jmin, jmax ! Min/max i,j for localization.
69 integer :: ijkv !
70 integer :: index !
71 integer :: xend
72 real :: num_members_inv
73 real :: num_members_inv1
74 real :: cov_inf_fac ! Covariance inflation factor.
75 real :: cov_loc_rad_m ! Covariance localization radius (m).
76 real :: cov_loc_rad ! Covariance localization radius (gridpts).
77 real :: ds ! Grid resolution (m).
78 real :: x1, y1 ! Grid location of ob.
79 real :: sigma_o2, sigma_b2
80 real :: sigma_hh_inv
81 real :: sigma_xh
82 real :: beta
83 real :: o_minus_b ! Observation minus forecast.
84 real :: kalman_gain
85 real :: x_mean
86 real :: h_mean
87 real :: hb_mean
88 real :: ha_mean
89 real :: sum_f, sum_a
90 real :: hf, ha
91
92 type (ob_type),pointer:: ob(:)
93
94 character(len=10),pointer:: varr(:) ! Variable name.
95 integer, pointer :: nkk(:) ! Vertical dimension of field.
96 integer, pointer :: xstart(:) ! Starting position of variable..
97
98 real, pointer :: x(:,:), xf(:,:)
99 ! real, pointer :: xf_cif(:,:)
100 real, pointer :: h(:), h_devn(:)
101 real, pointer :: x_devn(:)
102
103 real, pointer :: uc(:,:) ! u-wind (C grid).
104 real, pointer :: vc(:,:) ! v-wind (C grid).
105 real, pointer :: field(:,:) ! Input field.
106 real, pointer :: dummy(:,:) ! Dummy.
107 real, pointer :: corr(:,:) ! Correlation.
108 real, pointer :: xf_mean(:) ! Prior mean.
109 real, pointer :: xa_mean(:) ! Posterior mean.
110
111 namelist / gen_be_ensrf_nl / filestub, num_members, &
112 cov_inf_fac, cov_loc_rad_m
113
114 stderr = 0
115 stdout = 6
116
117 !---------------------------------------------------------------------------------------------
118 write(6,'(/a)')' [1] Initialize information.'
119 !---------------------------------------------------------------------------------------------
120
121 if (trace_use) call da_trace_init
122 if (trace_use) call da_trace_entry("gen_be_ensrf")
123
124 filestub = 'test'
125 num_members = 56
126 cov_inf_fac = 1.0
127 cov_loc_rad_m = 1500000.0
128
129 open(unit=unit, file='gen_be_ensrf_nl.nl', &
130 form='formatted', status='old', action='read')
131 read(unit, gen_be_ensrf_nl)
132 close(unit)
133
134 write(6,'(a,a)')' Filestub = ', trim(filestub)
135 write(6,'(a,i4)')' Number of ensemble members = ', num_members
136 write(6,'(a,f16.8)')' Covariance Inflation Factor = ', cov_inf_fac
137 write(6,'(a,f16.8)')' Covariance Localization Radius (m) = ', cov_loc_rad_m
138
139 num_members_inv = 1.0 / real(num_members)
140 num_members_inv1 = 1.0 / real(num_members-1)
141
142 !---------------------------------------------------------------------------------------------
143 write(6,'(/a)')' [2] Set up data dimensions and allocate arrays:'
144 !---------------------------------------------------------------------------------------------
145
146 ! Get grid dimensions from first T field:
147 input_file = trim(filestub)//'.e001'
148 var = 'T'
149 call da_stage0_initialize( input_file, var, ni, nj, nk, ds )
150 niu = ni+1 ! u i dimension is 1 larger.
151 njv = nj+1 ! v j dimension is 1 larger.
152 cov_loc_rad = cov_loc_rad_m / ds
153 nij = ni * nj
154 nijk = nij * nk
155 nijkv = nv3d * nijk + nv2d * nij ! #3D + #2D variables.
156 nv = nv3d + nv2d
157 print *, ni, nj, nk, nv, nij, nijk, nv3d, nv2d, nijkv
158 ! Allocate arrays:
159 allocate ( varr(1:nv) )
160 allocate ( nkk(1:nv) )
161 allocate ( xstart(1:nv) )
162 allocate( x(1:nijkv,1:num_members) )
163 allocate( xf(1:nijkv,1:num_members) )
164 ! allocate( xf_cif(1:nijkv,1:num_members) )
165 allocate( h(1:num_members) )
166 allocate( h_devn(1:num_members) )
167 allocate( x_devn(1:num_members) )
168 do v = 1, nv3d
169 nkk(v) = nk
170 end do
171 do v = nv3d + 1, nv
172 nkk(v) = 1
173 end do
174 allocate( field(1:ni,1:nj) )
175 allocate( uc(1:niu,1:nj) )
176 allocate( vc(1:ni,1:njv) )
177 allocate( dummy(1:ni,1:nj) )
178 allocate( corr(1:ni,1:nj) )
179
180 ! Hardwired:
181 varr(1) = 'U'
182 varr(2) = 'V'
183 varr(3) = 'T'
184 varr(4) = 'QVAPOR'
185 varr(5) = 'PSFC'
186
187 !---------------------------------------------------------------------------------------------
188 write(6,'(/a)')' [3] Read observations.'
189 !---------------------------------------------------------------------------------------------
190
191 input_file = 'observations'
192 open(unit, file = input_file, status='old')
193 !write(56,*)1
194 ! write(56,'(2a10,i6,5f16.8)')'test', 'dale', 3, &
195 ! 20.0, 20.0, 12.0, 250.0, 1.0
196 !stop
197 read(unit,*)num_obs
198 write(6,'(a,i10)')' Number of observations = ', num_obs
199 allocate( ob(1:num_obs) )
200
201 do o = 1, num_obs
202 read(unit,'(2a10,i6,5f16.8)')ob(o) % info % id, ob(o) % info % platform, &
203 ob(o) % info % variable, &
204 ob(o) % loc % x, ob(o) % loc % y ,ob(o) % loc % z, &
205 ob(o) % yo, ob(o) % sigma_o
206 end do
207
208 !---------------------------------------------------------------------------------------------
209 write(6,'(/a)')' [4] Extract necessary fields from WRF ensemble forecasts.'
210 !---------------------------------------------------------------------------------------------
211
212 xstart(1) = 1
213 do v = 2, nv
214 xstart(v) = xstart(v-1) + ( ni * nj * nkk(v-1) )
215 end do
216
217 do member = 1, num_members
218
219 write(UNIT=ce,FMT='(i3.3)')member
220 input_file = trim(filestub)//'.e'//ce
221
222 do v = 1, nv
223 do k = 1, nkk(v)
224
225 if ( varr(v) == 'U' ) then
226 call da_get_field( input_file, varr(v), 3, niu, nj, nk, k, uc )
227 do j = 1, nj
228 do i = 1, ni
229 field(i,j) = 0.5 * ( uc(i,j) + uc(i+1,j) )
230 end do
231 end do
232 end if
233
234 if ( varr(v) == 'V' ) then
235 call da_get_field( input_file, varr(v), 3, ni, njv, nk, k, vc )
236 do j = 1, nj
237 do i = 1, ni
238 field(i,j) = 0.5 * ( vc(i,j) + vc(i,j+1) )
239 end do
240 end do
241 end if
242
243 if ( varr(v) == "T" ) then
244 ! Read theta, and convert to temperature:
245 call da_get_trh( input_file, ni, nj, nk, k, field, dummy )
246 end if
247
248 ! Read mixing ratio, and convert to specific humidity:
249 if ( varr(v) == 'QVAPOR' ) then
250 call da_get_field( input_file, varr(v), 3, ni, nj, nk, k, field )
251 field(:,:) = field(:,:) / ( 1.0 + field(:,:) )
252 end if
253
254 ! Read surface pressure:
255 if ( varr(v) == 'PSFC' ) then
256 call da_get_field( input_file, varr(v), 2, ni, nj, nk, k, field )
257 end if
258
259 ! Fill 4D array:
260 index = xstart(v) + (k-1) * nij
261 do j = 1, nj
262 do i = 1, ni
263 xf(index,member) = field(i,j)
264 index = index + 1
265 end do
266 end do
267 end do ! k
268 end do ! v
269 end do !member
270
271 ! Initialize analysis as first guess:
272 x(1:nijkv,1:num_members) = xf(1:nijkv,1:num_members)
273
274 ! Perform initial covariance inflation (to boost input be):
275 call da_cov_inflation( nijkv, num_members, cov_inf_fac, x )
276 ! xf_cif(1:nijkv,1:num_members) = x(1:nijkv,1:num_members) ! Diagnostic
277
278 deallocate( uc )
279 deallocate( vc )
280 deallocate( field )
281 deallocate( dummy )
282
283 !---------------------------------------------------------------------------------------------
284 write(6,'(/a)')' [5] Perform EnSRF:'
285 !---------------------------------------------------------------------------------------------
286
287 do o = 1, num_obs
288 write(6,'(2(a,i8))')' Assimilating observation ', o, ' of ', num_obs
289 x1 = ob(o) % loc % x
290 y1 = ob(o) % loc % y
291 sigma_o2 = ob(o) % sigma_o * ob(o) % sigma_o
292
293 ! Perform observation operator:
294
295 do member = 1, num_members
296 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), x(:,member), h(member) )
297 end do
298
299 ! Mean and deviations of H(x):
300 h_mean = sum(h(1:num_members)) * num_members_inv
301 h_devn(1:num_members) = h(1:num_members) - h_mean
302
303 ! Covariance H Pf H^T:
304 sigma_b2 = sum(h_devn(1:num_members) * h_devn(1:num_members)) * num_members_inv1
305
306 ob(o) % sigma_b = sqrt(sigma_b2)
307
308 ! Store 1/(H Pf HT + R):
309 sigma_hh_inv = 1.0 / ( sigma_b2 + sigma_o2 )
310
311 ! Calculate EnSRF deviation Kalman Gain scaling:
312 beta = 1.0 / ( 1.0 + sqrt( sigma_o2 * sigma_hh_inv ) )
313
314 ! Ensemble mean innovation vector:
315 !!!DALE Uncomment for pseudo observation
316 ! print *, 'Running pseudo observation O-B=1.0:'
317 ! ob(o) % yo = h_mean + 1.0
318 !!!DALE
319 o_minus_b = ob(o) % yo - h_mean
320 write(6,'(a,2f15.5)')' Observation increment = ', o_minus_b
321 write(6,'(a,2f15.5)')' Observation/forecast error variance = ', sigma_o2, sigma_b2
322
323 !-----------------------------------------------------------------------------
324 ! Calculate local support:
325 !-----------------------------------------------------------------------------
326
327 imin = max(nint(ob(o) % loc % x - cov_loc_rad),1)
328 imax = min(nint(ob(o) % loc % x + cov_loc_rad),ni)
329 jmin = max(nint(ob(o) % loc % y - cov_loc_rad),1)
330 jmax = min(nint(ob(o) % loc % y + cov_loc_rad),nj)
331
332 ! write(6,'(a,2f8.2)')' Ob location(x/y) = ', ob(o) % loc % x, ob(o) % loc % y
333 ! write(6,'(a,4i5)')' Min/max i/j = ', imin, imax, jmin, jmax
334
335 do j = jmin, jmax
336 do i = imin, imax
337 call compact_support( x1, real(i), y1, real(j), cov_loc_rad, corr(i,j) )
338 end do
339 end do
340
341 do v = 1, nv
342 do k = 1, nkk(v)
343 do j = jmin, jmax
344 do i = imin, imax
345 ijkv = (v-1) * nijk + (k-1) * nij + (j-1) * ni + i
346
347 x_mean = sum(x(ijkv,1:num_members)) * num_members_inv
348
349 x_devn(1:num_members) = x(ijkv,1:num_members) - x_mean
350
351 ! Calculate PfHT:
352 sigma_xh = sum( x_devn(1:num_members) * h_devn(1:num_members) ) * &
353 num_members_inv1
354
355 ! Apply covariance localization:
356 sigma_xh = sigma_xh * corr(i,j)
357
358 kalman_gain = sigma_xh * sigma_hh_inv
359
360 x_mean = x_mean + kalman_gain * o_minus_b
361
362 x_devn(1:num_members) = x_devn(1:num_members) - &
363 beta * kalman_gain * h_devn(1:num_members)
364
365 x(ijkv,1:num_members) = x_mean + x_devn(1:num_members)
366
367 end do
368 end do
369 end do
370 end do
371 end do
372
373 !---------------------------------------------------------------------------------------------
374 write(6,'(/a)')' [4] Calculate diagnostics:'
375 !---------------------------------------------------------------------------------------------
376
377 ! Calculate mean analysis over all grid-points (diagnostic only):
378
379 allocate( xf_mean(1:nijkv) )
380 allocate( xa_mean(1:nijkv) )
381
382 do ijkv = 1, nijkv
383 xf_mean(ijkv) = sum(xf(ijkv,1:num_members)) * num_members_inv
384 xa_mean(ijkv) = sum(x(ijkv,1:num_members)) * num_members_inv
385 end do
386
387 ! Calculate grid-point diagnostics:
388
389 ! call da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf_cif, x )
390 call da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf, x )
391
392 ! Calculate observation diagnostics:
393
394 write(53,*)num_obs
395 do o = 1, num_obs
396
397 ! Ensemble mean-observation diagnostics:
398 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf_mean, hb_mean )
399 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xa_mean, ha_mean )
400
401 ob(o) % omb_mean = ob(o) % yo - hb_mean
402 ob(o) % oma_mean = ob(o) % yo - ha_mean
403
404 ! Ensemble spread-observation diagnostics:
405 sum_f = 0.0
406 sum_a = 0.0
407 do member = 1, num_members
408 ! call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf_cif(:,member), hf )
409 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf(:,member), hf )
410 call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), x(:,member), ha )
411 sum_f = sum_f + ( ob(o) % yo - hf )**2
412 sum_a = sum_a + ( ob(o) % yo - ha )**2
413 end do
414
415 ob(o) % omb_rms = sqrt( sum_f * num_members_inv )
416 ob(o) % oma_rms = sqrt( sum_a * num_members_inv )
417
418 write(53,'(2a10,i8,10f16.8)')ob(o) % info % id, ob(o) % info % platform, &
419 ob(o) % info % variable, &
420 ob(o) % loc % x, ob(o) % loc % y ,ob(o) % loc % z, &
421 ob(o) % yo, ob(o) % omb_mean, ob(o) % oma_mean, &
422 ob(o) % omb_rms, ob(o) % oma_rms, &
423 ob(o) % sigma_o, ob(o) % sigma_b
424 end do
425
426 !---------------------------------------------------------------------------------------------
427 write(6,'(/a)')' [4] Output EnSRF analysis ensemble:'
428 !---------------------------------------------------------------------------------------------
429
430 do member = 1, num_members
431 write(6,'(a,i4)')' Writing ensemble member increments for member', member
432 write(UNIT=ce,FMT='(i3.3)')member
433
434 do v = 1, nv
435
436 ! Output prior ensemble forecasts:
437 output_file = trim(varr(v))//'/'//trim(varr(v))//'.prior.e'//trim(ce)
438 open(unit, file = output_file, form='unformatted')
439 write(unit)ni, nj, nkk(v)
440 xend = xstart(v) + nij * nkk(v) - 1
441 write(unit)xf(xstart(v):xend,member)
442 close(unit)
443
444 ! Output posterior ensemble forecasts:
445 output_file = trim(varr(v))//'/'//trim(varr(v))//'.posterior.e'//trim(ce)
446 open(unit, file = output_file, form='unformatted')
447 write(unit)ni, nj, nkk(v)
448 xend = xstart(v) + nij * nkk(v) - 1
449 write(unit)x(xstart(v):xend,member)
450 close(unit)
451 end do
452 end do
453
454 do v = 1, nv
455
456 ! Output prior ensemble forecasts:
457 output_file = trim(varr(v))//'/'//trim(varr(v))//'.prior.mean'
458 open(unit, file = output_file, form='unformatted')
459 write(unit)ni, nj, nkk(v)
460 xend = xstart(v) + nij * nkk(v) - 1
461 write(unit)xf_mean(xstart(v):xend)
462 close(unit)
463
464 ! Output posterior ensemble forecasts:
465 output_file = trim(varr(v))//'/'//trim(varr(v))//'.posterior.mean'
466 open(unit, file = output_file, form='unformatted')
467 write(unit)ni, nj, nkk(v)
468 xend = xstart(v) + nij * nkk(v) - 1
469 write(unit)xa_mean(xstart(v):xend)
470 close(unit)
471 end do
472
473 if (trace_use) call da_trace_exit("gen_be_ensrf")
474 if (trace_use) call da_trace_report
475
476 deallocate( ob )
477
478 contains
479
480 !-----------------------------------------------------------------------------
481
482 subroutine da_obs_operator( nv, ni, nj, nkk, xstart, ob, x, h )
483
484 implicit none
485
486 integer, intent(in) :: nv, ni, nj
487 integer, intent(in) :: nkk(1:nv)
488 integer, intent(in) :: xstart(1:nv)
489 type (ob_type), intent(in) :: ob
490 real, intent(in), target :: x(:)
491 real, intent(out) :: h
492
493 integer :: v, nk
494 integer :: start
495 real :: ob_x, ob_y, ob_z ! Observation location.
496
497 v = ob % info % variable
498 nk = nkk(v)
499 start = xstart( v )
500
501 ! Perform 3D interpolation:
502
503 ob_x = ob % loc % x
504 ob_y = ob % loc % y
505 ob_z = ob % loc % z
506
507 call da_interpolate_3d( ni, nj, nk, ob % loc % x, ob % loc % y, ob % loc % z, &
508 x(start:start + ni * nj * nk ), h )
509
510 end subroutine da_obs_operator
511
512 subroutine da_interpolate_3d( ni, nj, nk, x, y, z, field, h )
513
514 implicit none
515 integer, intent(in) :: ni, nj, nk
516 real, intent(in) :: x, y, z ! Grid location of point.
517 real, intent(in) :: field(1:ni,1:nj,1:nk) ! Field to interpolate.
518 real, intent(out) :: h ! Interpolated value.
519
520 integer :: ii, jj, kk ! Grid locators.
521 real :: dx, dy, dz ! Interpolation weights.
522 real :: h1, h2 ! Interpolation values.
523 real :: h_below ! 2D interp values.
524 real :: h_above ! 2D interp values.
525
526 ii = int(x)
527 jj = int(y)
528 kk = int(z)
529
530 dx = x - real(ii)
531 dy = y - real(jj)
532 dz = z - real(kk)
533
534 ! 2D interpolation on layer below:
535 h1 = ( 1.0 - dx ) * field(ii,jj,kk) + dx * field(ii+1,jj,kk)
536 h2 = ( 1.0 - dx ) * field(ii,jj+1,kk) + dx * field(ii+1,jj+1,kk)
537 h_below = ( 1.0 - dy ) * h1 + dy * h2
538
539 ! 2D interpolation on layer above:
540 h1 = ( 1.0 - dx ) * field(ii,jj,kk+1) + dx * field(ii+1,jj,kk+1)
541 h2 = ( 1.0 - dx ) * field(ii,jj+1,kk+1) + dx * field(ii+1,jj+1,kk+1)
542 h_above = ( 1.0 - dy ) * h1 + dy * h2
543
544 ! Interpolation in vertical:
545 h = ( 1.0 - dz ) * h_below + dz * h_above
546
547 end subroutine da_interpolate_3d
548
549 subroutine compact_support( x1, x2, y1, y2, cov_loc_rad, corr )
550
551 ! Compact support according to (4.10) of Gaspari and Cohn (1999).
552
553 ! Assumes r >=0.
554
555 implicit none
556
557 real, intent(in) :: x1, x2 ! x values of two points.
558 real, intent(in) :: y1, y2 ! y values of two points.
559 real, intent(in) :: cov_loc_rad ! Cut-off correlation > 2c.
560 real, intent(out) :: corr ! Compactly-supported correlation.
561
562 real :: z ! Distance between points.
563 real :: r ! Ratio.
564 real :: r2, r3, r4, r5 ! Powers of r.
565
566 z = sqrt( ( x2 - x1 )**2 + ( y2 - y1)**2 )
567 r = z / cov_loc_rad
568
569 r2 = r * r
570 r3 = r2 * r
571 r4 = r3 * r
572 r5 = r4 * r
573
574 if ( r <= 1.0 ) then
575 corr = -0.25 * r5 + 0.5 * r4 + 0.625 * r3 - 5.0 * r2 / 3.0 + 1
576 else if ( r > 1.0 .and. r < 2.0 ) then
577 corr = r5 / 12.0 - 0.5 * r4 + 0.625 * r3 + 5.0 * r2 / 3.0 - &
578 5.0 * r + 4.0 - 2.0 / ( 3.0 * r )
579
580 else if ( r >= 2.0 ) then
581 corr = 0.0
582 end if
583
584 end subroutine compact_support
585
586 subroutine da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf, xa )
587
588 implicit none
589
590 integer, intent(in) :: num_members
591 integer, intent(in) :: ni
592 integer, intent(in) :: nj
593 integer, intent(in) :: nk
594 integer, intent(in) :: nv
595 integer, intent(in) :: nijkv
596 integer, intent(in) :: nkk(1:nv)
597 character*10,intent(in):: varr(1:nv)
598 real, intent(in) :: xf(1:nijkv,1:num_members)
599 real, intent(in) :: xa(1:nijkv,1:num_members)
600
601 character*10 :: name
602 integer :: v, k, j, i, ijkv
603 real :: num_members_inv, nij_inv
604 real :: mnsq_f, mnsq_a
605 real :: domain_mean_f, domain_mean_a
606 real :: domain_stdv_f, domain_stdv_a
607 real :: mean_f(1:ni,1:nj), mean_a(1:ni,1:nj)
608 real :: stdv_f(1:ni,1:nj), stdv_a(1:ni,1:nj)
609
610 num_members_inv = 1.0 / num_members
611 nij_inv = 1.0 / real ( ni * nj )
612
613 write(6,'(a)')' Variable, Level, Mean_F, Mean_A, StDv_F, StdV_A'
614
615 ijkv = 0
616 do v = 1, nv
617 name = varr(v)
618 do k = 1, nkk(v)
619 do j = 1, nj
620 do i = 1, ni
621 ijkv = ijkv + 1
622 mean_f(i,j) = sum(xf(ijkv,1:num_members)) * num_members_inv
623 mean_a(i,j) = sum(xa(ijkv,1:num_members)) * num_members_inv
624 mnsq_f = sum(xf(ijkv,1:num_members)**2) * num_members_inv
625 mnsq_a = sum(xa(ijkv,1:num_members)**2) * num_members_inv
626 stdv_f(i,j) = mnsq_f - mean_f(i,j) * mean_f(i,j)
627 stdv_a(i,j) = mnsq_a - mean_a(i,j) * mean_a(i,j)
628 stdv_f(i,j) = 0.0
629 if ( stdv_f(i,j) > 0.0 ) stdv_f(i,j) = sqrt(stdv_f(i,j))
630 stdv_a(i,j) = 0.0
631 if ( stdv_a(i,j) > 0.0 ) stdv_a(i,j) = sqrt(stdv_a(i,j))
632 end do
633 end do
634
635 domain_mean_f = sum(mean_f(:,:)) * nij_inv
636 domain_mean_a = sum(mean_a(:,:)) * nij_inv
637 domain_stdv_f = sum(stdv_f(:,:)) * nij_inv
638 domain_stdv_a = sum(stdv_a(:,:)) * nij_inv
639
640 write(6,'(a,i6,4f16.8)')trim(name), k, domain_mean_f, domain_mean_a, &
641 domain_stdv_f, domain_stdv_a
642 end do
643 end do
644
645 end subroutine da_get_grid_stats
646
647 subroutine da_cov_inflation( nijkv, num_members, cov_inf_fac, x )
648
649 implicit none
650
651 integer, intent(in) :: nijkv
652 integer, intent(in) :: num_members
653 real, intent(in) :: cov_inf_fac
654 real, intent(inout) :: x(1:nijkv,1:num_members)
655
656 integer :: ijkv
657 real :: num_members_inv
658 real :: x_devn(1:num_members)
659
660 num_members_inv = 1.0 / real(num_members)
661
662 do ijkv = 1, nijkv
663 x_mean = sum(x(ijkv,1:num_members)) * num_members_inv
664 x_devn(1:num_members) = cov_inf_fac * ( x(ijkv,1:num_members) - x_mean )
665 x(ijkv,1:num_members) = x_mean + x_devn(1:num_members)
666 end do
667
668 end subroutine da_cov_inflation
669
670 end program gen_be_ensrf
671