<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='GEN_BE_ENSRF'><A href='../../html_code/gen_be/gen_be_ensrf.f90.html#GEN_BE_ENSRF' TARGET='top_target'><IMG SRC="../../gif/bar_yellow.gif" border=0></A>
program gen_be_ensrf,9
!
!----------------------------------------------------------------------
! Purpose : Perform an Ensemble Square Root Filter (EnSRF) assimilation
! using WRF forecast data.
!
! Owner: Dale Barker (NCAR/MMM)
! Please acknowledge author/institute in work that uses this code.
!
!----------------------------------------------------------------------
#ifdef crayx1
#define iargc ipxfargc
#endif
use da_control
, only : stdout,stderr,filename_len
use da_gen_be, only : da_stage0_initialize, da_get_field, da_get_trh
implicit none
type info_type
character*40 :: id
character*40 :: platform
integer :: plat
integer :: variable
end type info_type
type loc_type
real :: x
real :: y
real :: z
end type loc_type
type ob_type
type (info_type):: info
type (loc_type) :: loc
real :: yo
real :: omb_mean
real :: oma_mean
real :: omb_rms
real :: oma_rms
real :: sigma_o
real :: sigma_b
end type ob_type
integer, parameter :: nv3d = 4 ! #3D variables (u, v, T, q).
integer, parameter :: nv2d = 1 ! #2D variables (ps).
integer, parameter :: unit = 100 ! Unit number.
character (len=filename_len) :: filestub ! General filename stub.
character (len=filename_len) :: input_file ! Input file.
character (len=filename_len) :: output_file ! Output file.
character (len=10) :: var ! Variable to search for.
character (len=3) :: ce ! Member index -> character.
integer :: num_members ! Ensemble size.
integer :: num_obs ! Number of observations.
integer :: ni ! 1st dimension size.
integer :: niu ! 1st dimension size (u)
integer :: nj ! 2nd dimension size.
integer :: njv ! 2nd dimension size (v)
integer :: nk ! 3rd dimension size.
integer :: nv ! Total number of 2D/3D variables.
integer :: nij, nijk, nijkv ! Dimensions.
integer :: o, member, i, j, k, v ! Loop counters.
integer :: imin, imax, jmin, jmax ! Min/max i,j for localization.
integer :: ijkv !
integer :: index !
integer :: xend
real :: num_members_inv
real :: num_members_inv1
real :: cov_inf_fac ! Covariance inflation factor.
real :: cov_loc_rad_m ! Covariance localization radius (m).
real :: cov_loc_rad ! Covariance localization radius (gridpts).
real :: ds ! Grid resolution (m).
real :: x1, y1 ! Grid location of ob.
real :: sigma_o2, sigma_b2
real :: sigma_hh_inv
real :: sigma_xh
real :: beta
real :: o_minus_b ! Observation minus forecast.
real :: kalman_gain
real :: x_mean
real :: h_mean
real :: hb_mean
real :: ha_mean
real :: sum_f, sum_a
real :: hf, ha
type (ob_type),pointer:: ob(:)
character(len=10),allocatable :: varr(:) ! Variable name.
integer, allocatable :: nkk(:) ! Vertical dimension of field.
integer, allocatable :: xstart(:) ! Starting position of variable..
real, allocatable :: x(:,:), xf(:,:)
! real, allocatable :: xf_cif(:,:)
real, allocatable :: h(:), h_devn(:)
real, allocatable :: x_devn(:)
real, allocatable :: uc(:,:) ! u-wind (C grid).
real, allocatable :: vc(:,:) ! v-wind (C grid).
real, allocatable :: field(:,:) ! Input field.
real, allocatable :: dummy(:,:) ! Dummy.
real, allocatable :: corr(:,:) ! Correlation.
real, allocatable :: xf_mean(:) ! Prior mean.
real, allocatable :: xa_mean(:) ! Posterior mean.
namelist / gen_be_ensrf_nl / filestub, num_members, &
cov_inf_fac, cov_loc_rad_m
stderr = 0
stdout = 6
!---------------------------------------------------------------------------------------------
write(6,'(/a)')' [1] Initialize information.'
!---------------------------------------------------------------------------------------------
filestub = 'test'
num_members = 56
cov_inf_fac = 1.0
cov_loc_rad_m = 1500000.0
open(unit=unit, file='gen_be_ensrf_nl.nl', &
form='formatted', status='old', action='read')
read(unit, gen_be_ensrf_nl)
close(unit)
write(6,'(a,a)')' Filestub = ', trim(filestub)
write(6,'(a,i4)')' Number of ensemble members = ', num_members
write(6,'(a,f16.8)')' Covariance Inflation Factor = ', cov_inf_fac
write(6,'(a,f16.8)')' Covariance Localization Radius (m) = ', cov_loc_rad_m
num_members_inv = 1.0 / real(num_members)
num_members_inv1 = 1.0 / real(num_members-1)
!---------------------------------------------------------------------------------------------
write(6,'(/a)')' [2] Set up data dimensions and allocate arrays:'
!---------------------------------------------------------------------------------------------
! Get grid dimensions from first T field:
input_file = trim(filestub)//'.e001'
var = 'T'
call da_stage0_initialize( input_file, var, ni, nj, nk, ds )
niu = ni+1 ! u i dimension is 1 larger.
njv = nj+1 ! v j dimension is 1 larger.
cov_loc_rad = cov_loc_rad_m / ds
nij = ni * nj
nijk = nij * nk
nijkv = nv3d * nijk + nv2d * nij ! #3D + #2D variables.
nv = nv3d + nv2d
print *, ni, nj, nk, nv, nij, nijk, nv3d, nv2d, nijkv
! Allocate arrays:
allocate ( varr(1:nv) )
allocate ( nkk(1:nv) )
allocate ( xstart(1:nv) )
allocate( x(1:nijkv,1:num_members) )
allocate( xf(1:nijkv,1:num_members) )
! allocate( xf_cif(1:nijkv,1:num_members) )
allocate( h(1:num_members) )
allocate( h_devn(1:num_members) )
allocate( x_devn(1:num_members) )
do v = 1, nv3d
nkk(v) = nk
end do
do v = nv3d + 1, nv
nkk(v) = 1
end do
allocate( field(1:ni,1:nj) )
allocate( uc(1:niu,1:nj) )
allocate( vc(1:ni,1:njv) )
allocate( dummy(1:ni,1:nj) )
allocate( corr(1:ni,1:nj) )
! Hardwired:
varr(1) = 'U'
varr(2) = 'V'
varr(3) = 'T'
varr(4) = 'QVAPOR'
varr(5) = 'PSFC'
!---------------------------------------------------------------------------------------------
write(6,'(/a)')' [3] Read observations.'
!---------------------------------------------------------------------------------------------
input_file = 'observations'
open(unit, file = input_file, status='old')
!write(56,*)1
! write(56,'(2a10,i6,5f16.8)')'test', 'dale', 3, &
! 20.0, 20.0, 12.0, 250.0, 1.0
!stop
read(unit,*)num_obs
write(6,'(a,i10)')' Number of observations = ', num_obs
allocate( ob(1:num_obs) )
do o = 1, num_obs
read(unit,'(2a10,i6,5f16.8)')ob(o) % info % id, ob(o) % info % platform, &
ob(o) % info % variable, &
ob(o) % loc % x, ob(o) % loc % y ,ob(o) % loc % z, &
ob(o) % yo, ob(o) % sigma_o
end do
!---------------------------------------------------------------------------------------------
write(6,'(/a)')' [4] Extract necessary fields from WRF ensemble forecasts.'
!---------------------------------------------------------------------------------------------
xstart(1) = 1
do v = 2, nv
xstart(v) = xstart(v-1) + ( ni * nj * nkk(v-1) )
end do
do member = 1, num_members
write(UNIT=ce,FMT='(i3.3)')member
input_file = trim(filestub)//'.e'//ce
do v = 1, nv
do k = 1, nkk(v)
if ( varr(v) == 'U' ) then
call da_get_field( input_file, varr(v), 3, niu, nj, nk, k, uc )
do j = 1, nj
do i = 1, ni
field(i,j) = 0.5 * ( uc(i,j) + uc(i+1,j) )
end do
end do
end if
if ( varr(v) == 'V' ) then
call da_get_field( input_file, varr(v), 3, ni, njv, nk, k, vc )
do j = 1, nj
do i = 1, ni
field(i,j) = 0.5 * ( vc(i,j) + vc(i,j+1) )
end do
end do
end if
if ( varr(v) == "T" ) then
! Read theta, and convert to temperature:
call da_get_trh( input_file, ni, nj, nk, k, field, dummy )
end if
! Read mixing ratio, and convert to specific humidity:
if ( varr(v) == 'QVAPOR' ) then
call da_get_field( input_file, varr(v), 3, ni, nj, nk, k, field )
field(:,:) = field(:,:) / ( 1.0 + field(:,:) )
end if
! Read surface pressure:
if ( varr(v) == 'PSFC' ) then
call da_get_field( input_file, varr(v), 2, ni, nj, nk, k, field )
end if
! Fill 4D array:
index = xstart(v) + (k-1) * nij
do j = 1, nj
do i = 1, ni
xf(index,member) = field(i,j)
index = index + 1
end do
end do
end do ! k
end do ! v
end do !member
! Initialize analysis as first guess:
x(1:nijkv,1:num_members) = xf(1:nijkv,1:num_members)
! Perform initial covariance inflation (to boost input be):
call da_cov_inflation( nijkv, num_members, cov_inf_fac, x )
! xf_cif(1:nijkv,1:num_members) = x(1:nijkv,1:num_members) ! Diagnostic
deallocate( uc )
deallocate( vc )
deallocate( field )
deallocate( dummy )
!---------------------------------------------------------------------------------------------
write(6,'(/a)')' [5] Perform EnSRF:'
!---------------------------------------------------------------------------------------------
do o = 1, num_obs
write(6,'(2(a,i8))')' Assimilating observation ', o, ' of ', num_obs
x1 = ob(o) % loc % x
y1 = ob(o) % loc % y
sigma_o2 = ob(o) % sigma_o * ob(o) % sigma_o
! Perform observation operator:
do member = 1, num_members
call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), x(:,member), h(member) )
end do
! Mean and deviations of H(x):
h_mean = sum(h(1:num_members)) * num_members_inv
h_devn(1:num_members) = h(1:num_members) - h_mean
! Covariance H Pf H^T:
sigma_b2 = sum(h_devn(1:num_members) * h_devn(1:num_members)) * num_members_inv1
ob(o) % sigma_b = sqrt(sigma_b2)
! Store 1/(H Pf HT + R):
sigma_hh_inv = 1.0 / ( sigma_b2 + sigma_o2 )
! Calculate EnSRF deviation Kalman Gain scaling:
beta = 1.0 / ( 1.0 + sqrt( sigma_o2 * sigma_hh_inv ) )
! Ensemble mean innovation vector:
!!!DALE Uncomment for pseudo observation
! print *, 'Running pseudo observation O-B=1.0:'
! ob(o) % yo = h_mean + 1.0
!!!DALE
o_minus_b = ob(o) % yo - h_mean
write(6,'(a,2f15.5)')' Observation increment = ', o_minus_b
write(6,'(a,2f15.5)')' Observation/forecast error variance = ', sigma_o2, sigma_b2
!-----------------------------------------------------------------------------
! Calculate local support:
!-----------------------------------------------------------------------------
imin = max(nint(ob(o) % loc % x - cov_loc_rad),1)
imax = min(nint(ob(o) % loc % x + cov_loc_rad),ni)
jmin = max(nint(ob(o) % loc % y - cov_loc_rad),1)
jmax = min(nint(ob(o) % loc % y + cov_loc_rad),nj)
! write(6,'(a,2f8.2)')' Ob location(x/y) = ', ob(o) % loc % x, ob(o) % loc % y
! write(6,'(a,4i5)')' Min/max i/j = ', imin, imax, jmin, jmax
do j = jmin, jmax
do i = imin, imax
call compact_support( x1, real(i), y1, real(j), cov_loc_rad, corr(i,j) )
end do
end do
do v = 1, nv
do k = 1, nkk(v)
do j = jmin, jmax
do i = imin, imax
ijkv = (v-1) * nijk + (k-1) * nij + (j-1) * ni + i
x_mean = sum(x(ijkv,1:num_members)) * num_members_inv
x_devn(1:num_members) = x(ijkv,1:num_members) - x_mean
! Calculate PfHT:
sigma_xh = sum( x_devn(1:num_members) * h_devn(1:num_members) ) * &
num_members_inv1
! Apply covariance localization:
sigma_xh = sigma_xh * corr(i,j)
kalman_gain = sigma_xh * sigma_hh_inv
x_mean = x_mean + kalman_gain * o_minus_b
x_devn(1:num_members) = x_devn(1:num_members) - &
beta * kalman_gain * h_devn(1:num_members)
x(ijkv,1:num_members) = x_mean + x_devn(1:num_members)
end do
end do
end do
end do
end do
!---------------------------------------------------------------------------------------------
write(6,'(/a)')' [4] Calculate diagnostics:'
!---------------------------------------------------------------------------------------------
! Calculate mean analysis over all grid-points (diagnostic only):
allocate( xf_mean(1:nijkv) )
allocate( xa_mean(1:nijkv) )
do ijkv = 1, nijkv
xf_mean(ijkv) = sum(xf(ijkv,1:num_members)) * num_members_inv
xa_mean(ijkv) = sum(x(ijkv,1:num_members)) * num_members_inv
end do
! Calculate grid-point diagnostics:
! call da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf_cif, x )
call da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf, x )
! Calculate observation diagnostics:
write(53,*)num_obs
do o = 1, num_obs
! Ensemble mean-observation diagnostics:
call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf_mean, hb_mean )
call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xa_mean, ha_mean )
ob(o) % omb_mean = ob(o) % yo - hb_mean
ob(o) % oma_mean = ob(o) % yo - ha_mean
! Ensemble spread-observation diagnostics:
sum_f = 0.0
sum_a = 0.0
do member = 1, num_members
! call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf_cif(:,member), hf )
call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf(:,member), hf )
call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), x(:,member), ha )
sum_f = sum_f + ( ob(o) % yo - hf )**2
sum_a = sum_a + ( ob(o) % yo - ha )**2
end do
ob(o) % omb_rms = sqrt( sum_f * num_members_inv )
ob(o) % oma_rms = sqrt( sum_a * num_members_inv )
write(53,'(2a10,i8,10f16.8)')ob(o) % info % id, ob(o) % info % platform, &
ob(o) % info % variable, &
ob(o) % loc % x, ob(o) % loc % y ,ob(o) % loc % z, &
ob(o) % yo, ob(o) % omb_mean, ob(o) % oma_mean, &
ob(o) % omb_rms, ob(o) % oma_rms, &
ob(o) % sigma_o, ob(o) % sigma_b
end do
!---------------------------------------------------------------------------------------------
write(6,'(/a)')' [4] Output EnSRF analysis ensemble:'
!---------------------------------------------------------------------------------------------
do member = 1, num_members
write(6,'(a,i4)')' Writing ensemble member increments for member', member
write(UNIT=ce,FMT='(i3.3)')member
do v = 1, nv
! Output prior ensemble forecasts:
output_file = trim(varr(v))//'/'//trim(varr(v))//'.prior.e'//trim(ce)
open(unit, file = output_file, form='unformatted')
write(unit)ni, nj, nkk(v)
xend = xstart(v) + nij * nkk(v) - 1
write(unit)xf(xstart(v):xend,member)
close(unit)
! Output posterior ensemble forecasts:
output_file = trim(varr(v))//'/'//trim(varr(v))//'.posterior.e'//trim(ce)
open(unit, file = output_file, form='unformatted')
write(unit)ni, nj, nkk(v)
xend = xstart(v) + nij * nkk(v) - 1
write(unit)x(xstart(v):xend,member)
close(unit)
end do
end do
do v = 1, nv
! Output prior ensemble forecasts:
output_file = trim(varr(v))//'/'//trim(varr(v))//'.prior.mean'
open(unit, file = output_file, form='unformatted')
write(unit)ni, nj, nkk(v)
xend = xstart(v) + nij * nkk(v) - 1
write(unit)xf_mean(xstart(v):xend)
close(unit)
! Output posterior ensemble forecasts:
output_file = trim(varr(v))//'/'//trim(varr(v))//'.posterior.mean'
open(unit, file = output_file, form='unformatted')
write(unit)ni, nj, nkk(v)
xend = xstart(v) + nij * nkk(v) - 1
write(unit)xa_mean(xstart(v):xend)
close(unit)
end do
deallocate( ob )
contains
!-----------------------------------------------------------------------------
<A NAME='DA_OBS_OPERATOR'><A href='../../html_code/gen_be/gen_be_ensrf.f90.html#DA_OBS_OPERATOR' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_obs_operator( nv, ni, nj, nkk, xstart, ob, x, h ) 5,1
implicit none
integer, intent(in) :: nv, ni, nj
integer, intent(in) :: nkk(1:nv)
integer, intent(in) :: xstart(1:nv)
type (ob_type), intent(in) :: ob
! WHY? was this a target?
! real, intent(in), target :: x(:)
real, intent(in) :: x(:)
real, intent(out) :: h
integer :: v, nk
integer :: start
real :: ob_x, ob_y, ob_z ! Observation location.
v = ob % info % variable
nk = nkk(v)
start = xstart( v )
! Perform 3D interpolation:
ob_x = ob % loc % x
ob_y = ob % loc % y
ob_z = ob % loc % z
call da_interpolate_3d( ni, nj, nk, ob % loc % x, ob % loc % y, ob % loc % z, &
x(start:start + ni * nj * nk ), h )
end subroutine da_obs_operator
<A NAME='DA_INTERPOLATE_3D'><A href='../../html_code/gen_be/gen_be_ensrf.f90.html#DA_INTERPOLATE_3D' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_interpolate_3d( ni, nj, nk, x, y, z, field, h ) 1
implicit none
integer, intent(in) :: ni, nj, nk
real, intent(in) :: x, y, z ! Grid location of point.
real, intent(in) :: field(1:ni,1:nj,1:nk) ! Field to interpolate.
real, intent(out) :: h ! Interpolated value.
integer :: ii, jj, kk ! Grid locators.
real :: dx, dy, dz ! Interpolation weights.
real :: h1, h2 ! Interpolation values.
real :: h_below ! 2D interp values.
real :: h_above ! 2D interp values.
ii = int(x)
jj = int(y)
kk = int(z)
dx = x - real(ii)
dy = y - real(jj)
dz = z - real(kk)
! 2D interpolation on layer below:
h1 = ( 1.0 - dx ) * field(ii,jj,kk) + dx * field(ii+1,jj,kk)
h2 = ( 1.0 - dx ) * field(ii,jj+1,kk) + dx * field(ii+1,jj+1,kk)
h_below = ( 1.0 - dy ) * h1 + dy * h2
! 2D interpolation on layer above:
h1 = ( 1.0 - dx ) * field(ii,jj,kk+1) + dx * field(ii+1,jj,kk+1)
h2 = ( 1.0 - dx ) * field(ii,jj+1,kk+1) + dx * field(ii+1,jj+1,kk+1)
h_above = ( 1.0 - dy ) * h1 + dy * h2
! Interpolation in vertical:
h = ( 1.0 - dz ) * h_below + dz * h_above
end subroutine da_interpolate_3d
<A NAME='COMPACT_SUPPORT'><A href='../../html_code/gen_be/gen_be_ensrf.f90.html#COMPACT_SUPPORT' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine compact_support( x1, x2, y1, y2, cov_loc_rad, corr ) 1
! Compact support according to (4.10) of Gaspari and Cohn (1999).
! Assumes r >=0.
implicit none
real, intent(in) :: x1, x2 ! x values of two points.
real, intent(in) :: y1, y2 ! y values of two points.
real, intent(in) :: cov_loc_rad ! Cut-off correlation > 2c.
real, intent(out) :: corr ! Compactly-supported correlation.
real :: z ! Distance between points.
real :: r ! Ratio.
real :: r2, r3, r4, r5 ! Powers of r.
z = sqrt( ( x2 - x1 )**2 + ( y2 - y1)**2 )
r = z / cov_loc_rad
r2 = r * r
r3 = r2 * r
r4 = r3 * r
r5 = r4 * r
if ( r <= 1.0 ) then
corr = -0.25 * r5 + 0.5 * r4 + 0.625 * r3 - 5.0 * r2 / 3.0 + 1
else if ( r > 1.0 .and. r < 2.0 ) then
corr = r5 / 12.0 - 0.5 * r4 + 0.625 * r3 + 5.0 * r2 / 3.0 - &
5.0 * r + 4.0 - 2.0 / ( 3.0 * r )
else if ( r >= 2.0 ) then
corr = 0.0
end if
end subroutine compact_support
<A NAME='DA_GET_GRID_STATS'><A href='../../html_code/gen_be/gen_be_ensrf.f90.html#DA_GET_GRID_STATS' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf, xa ) 1
implicit none
integer, intent(in) :: num_members
integer, intent(in) :: ni
integer, intent(in) :: nj
integer, intent(in) :: nk
integer, intent(in) :: nv
integer, intent(in) :: nijkv
integer, intent(in) :: nkk(1:nv)
character*10,intent(in):: varr(1:nv)
real, intent(in) :: xf(1:nijkv,1:num_members)
real, intent(in) :: xa(1:nijkv,1:num_members)
character*10 :: name
integer :: v, k, j, i, ijkv
real :: num_members_inv, nij_inv
real :: mnsq_f, mnsq_a
real :: domain_mean_f, domain_mean_a
real :: domain_stdv_f, domain_stdv_a
real :: mean_f(1:ni,1:nj), mean_a(1:ni,1:nj)
real :: stdv_f(1:ni,1:nj), stdv_a(1:ni,1:nj)
num_members_inv = 1.0 / num_members
nij_inv = 1.0 / real ( ni * nj )
write(6,'(a)')' Variable, Level, Mean_F, Mean_A, StDv_F, StdV_A'
ijkv = 0
do v = 1, nv
name = varr(v)
do k = 1, nkk(v)
do j = 1, nj
do i = 1, ni
ijkv = ijkv + 1
mean_f(i,j) = sum(xf(ijkv,1:num_members)) * num_members_inv
mean_a(i,j) = sum(xa(ijkv,1:num_members)) * num_members_inv
mnsq_f = sum(xf(ijkv,1:num_members)**2) * num_members_inv
mnsq_a = sum(xa(ijkv,1:num_members)**2) * num_members_inv
stdv_f(i,j) = mnsq_f - mean_f(i,j) * mean_f(i,j)
stdv_a(i,j) = mnsq_a - mean_a(i,j) * mean_a(i,j)
stdv_f(i,j) = 0.0
if ( stdv_f(i,j) > 0.0 ) stdv_f(i,j) = sqrt(stdv_f(i,j))
stdv_a(i,j) = 0.0
if ( stdv_a(i,j) > 0.0 ) stdv_a(i,j) = sqrt(stdv_a(i,j))
end do
end do
domain_mean_f = sum(mean_f(:,:)) * nij_inv
domain_mean_a = sum(mean_a(:,:)) * nij_inv
domain_stdv_f = sum(stdv_f(:,:)) * nij_inv
domain_stdv_a = sum(stdv_a(:,:)) * nij_inv
write(6,'(a,i6,4f16.8)')trim(name), k, domain_mean_f, domain_mean_a, &
domain_stdv_f, domain_stdv_a
end do
end do
end subroutine da_get_grid_stats
<A NAME='DA_COV_INFLATION'><A href='../../html_code/gen_be/gen_be_ensrf.f90.html#DA_COV_INFLATION' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_cov_inflation( nijkv, num_members, cov_inf_fac, x ) 1
implicit none
integer, intent(in) :: nijkv
integer, intent(in) :: num_members
real, intent(in) :: cov_inf_fac
real, intent(inout) :: x(1:nijkv,1:num_members)
integer :: ijkv
real :: num_members_inv
real :: x_devn(1:num_members)
num_members_inv = 1.0 / real(num_members)
do ijkv = 1, nijkv
x_mean = sum(x(ijkv,1:num_members)) * num_members_inv
x_devn(1:num_members) = cov_inf_fac * ( x(ijkv,1:num_members) - x_mean )
x(ijkv,1:num_members) = x_mean + x_devn(1:num_members)
end do
end subroutine da_cov_inflation
end program gen_be_ensrf