<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='DA_SOLVE_ETKF'><A href='../../html_code/etkf/da_solve_etkf.inc.html#DA_SOLVE_ETKF' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_solve_etkf( ndim,nanals,nobs,ens,ens_ob,oberrvar,obs,nout,&,6
naccumt1,naccumt2,nstartaccum1,nstartaccum2,tainflatinput,rhoinput, &
infl_fac_file,eigen_val_file,inno2_val_file,proj2_val_file,infl_fac_TRNK, &
infl_fac_WG03,infl_fac_WG07,infl_fac_BOWL)
!-----------------------------------------------------------------------
! Purpose: ETKF perturbation updates
! Xuguang Wang, January 2006
! Dale Barker. January 2007 Implicit none, Convert to f90, and
! enclose within gen_be_etkf.f90 driver within WRF.
! Also modify inflation factor method (nstartaccum1,
! nstartaccum2 redundant).
! Chris Snyder, October 2008
! Arthur P. Mizzi, February 2011 Extensive modifications: add comments, clean up code,
! correct inflation factor algorithms, add other
! inflation factor algorithms, etc.
! NOTE: gen_be_etkf.f90 is also modified to match changes herein.
!
! references:
! Bishop et al 2001 MWR,
! Wang and Bishop 2003 MWR,
! Wang et al. 2004, MWR
! Wang et al. 2006, MWR
!
!1) nanals, ensemble size
!2) ndim, dimension of the perturbation vector that needs to be updated
!3) nobs, number of observations assimilated
!4) ens, array of perturbations before (Xf') and after (Xa') ETKF update
!5) ens_ob, array of HXf
!6) oberrvar, observation error variance, listed in the same sequence as HXf
!7) obs, observations assimilated
!8) naccumt1, number of previous cycles immediately before the current cycle,
! which is needed for calculating adaptive inflation factor.
! naccumt1 < 0 for pre-specified inflation
!9) naccumt2, number of previous cycles immediately before the current cycle,
! which is needed for calculating the rho factor in the latest version of ETKF.
! naccumt2 < 0 for using the older version of the ETKF.
! naccumt2 = 0 for using pre-specified rho factor
!10) nstartaccum1, the cycle from which the accumulation of previous naccumt1 cycles in 8) starts
!11) nstartaccum2, the cycle from which the accumulation of previous naccumt2 in 9) starts
!12) tainflatinput, pre-specified inflation, if not using adaptive inflation
!13) rhoinput, pre-specified rho factor, if not using adaptively determined rho factor
!14) nout, record number for output squared innovations and the forecast error variance
!---------------projected onto ensemble subspace, which is related to 8) and 9)
implicit none
integer, intent(in) :: nanals,ndim,nobs
real, intent(inout), dimension(ndim,nanals) :: ens
real, intent(inout), dimension(nobs,nanals) :: ens_ob
real, intent(in), dimension(nobs) :: oberrvar
real, intent(in), dimension(nobs) :: obs
real, intent(in) :: tainflatinput,rhoinput
integer, intent(in) :: nout,naccumt1,naccumt2,nstartaccum1,nstartaccum2
integer :: n, nstart, ndivd ! Loop counters.
integer :: nmin ! Minimum nout value to use.
real :: nanals_inv ! 1 / nanals.
real :: nanals_m1_inv ! 1 / (nanals-1).
real :: ainflat_mean ! Rolling mean inflation factor.
real :: ainflat_mean_old ! Rolling mean inflation factor read in from file.
real :: ainflat_tmp ! Instantaneous ratio of sqd innovations minus
! obs error to forecast variance at obs locations.
real, dimension(nobs) :: ensmean_ob
real, dimension(nobs) :: obsinc
integer :: ij, nanal, i, j, k ! Loop counters.
integer :: info, lwork
real, allocatable, dimension(:) :: work
real, dimension(1) :: workt
real, dimension(nanals) :: eignv, eignv1
real, dimension(nanals, nanals) :: hzthz, hzthz_sv,C, TR
real, dimension(nanals-1, nanals) :: CT
real, dimension(nanals, nanals-1) :: T
real, dimension(nanals, nanals-1) :: cgamma
real, dimension(nobs,nanals-1) :: E
real, dimension(nobs,nanals) :: enspert_ob_tmp
character (len=150) :: filename_factor
character (len=150) :: filename_eigenv
character (len=150) :: filename_inno2
character (len=150) :: filename_proj2
real, dimension(nanals-1) :: proj
real :: tracehpfht, tracehpfht_old
real :: ta, tb
real :: sum_eignv,sumv
real :: aftersum_eignv, trm
real :: ainflat, ainflat_tmp_sav, ainflat_mn
real :: tainflat, rho
real :: proj2, proj2sum
real :: proj2sum1
real :: squareinnosum, squareinnosum1
real :: squareinno, obsinc_sq
real :: squareinno_mn, squareinno_old, squareinno_sav
character (len=150):: infl_fac_file,eigen_val_file,inno2_val_file,proj2_val_file
logical:: infl_fac_TRNK,infl_fac_WG03,infl_fac_WG07,infl_fac_BOWL
!
! BEGIN EXECUTABLE CODE
if (trace_use) call da_trace_entry
("da_solve_etkf")
!
! SWITCH TO CHOOSE INFLATION FACTOR ALGORITHM
print *, "Enter ETKF with nout = ",nout
print *, "infl_fac_TRNK ", infl_fac_TRNK
print *, "infl_fac_WG03 ", infl_fac_WG03
print *, "infl_fac_WG07 ", infl_fac_WG07
print *, "infl_fac_BOWL ", infl_fac_BOWL
rho=1.
!
! OPEN THE INFLATION FACTOR HISTORY FILE
!
filename_factor = trim(infl_fac_file)
open(109,form="unformatted",file=filename_factor,access="direct",recl=8)
!
filename_eigenv = trim(eigen_val_file)
open(119,form="unformatted",file=filename_eigenv,access="direct",recl=8)
!
filename_inno2 = trim(inno2_val_file)
open(99,form="unformatted",file=filename_inno2,access="direct",recl=8)
!
filename_proj2 = trim(proj2_val_file)
open(89,form="unformatted",file=filename_proj2,access="direct",recl=8)
filename_proj2 = "ETKF_PARAMS"
open(120,form="unformatted",file=filename_proj2)
!
!------------------------------------------------------------------------------
print *, "[1] Compute mean(H(xf)) and H(xf)'"
!------------------------------------------------------------------------------
!
! CALCULATE THE ENSEMBLE MEAN AND CREATE ENSEMBLE PERTURBATIONS
! FOR THE OB.ETKF.exxx DATA
nanals_inv = 1.0 / real(nanals)
nanals_m1_inv = 1.0 / real(nanals-1)
do ij = 1, nobs
ensmean_ob(ij) = sum(ens_ob(ij,:)) * nanals_inv
enddo
do nanal = 1, nanals
ens_ob(:,nanal) = ens_ob(:,nanal) - ensmean_ob(:)
enddo
!------------------------------------------------------------------------------
print *, "[2] Calculate Transform(HZ)(HZ) in Bishop et al. 2001 Eq.(11)"
!------------------------------------------------------------------------------
!
! NORMALIZE H(xf)' BY SQRT(R) TO GET Tilda(H) THEN CALCULATE
! (HZ)T(HZ) IN EQN(11) AND NORMALIZE TO GET THE ENSEMBLE FORECAST
! COVARIANCE MATRIX IN hzths(nanals,nanals)
do i = 1, nobs
enspert_ob_tmp(i,:) = ens_ob(i,:)/sqrt(oberrvar(i))
end do
call da_innerprod(enspert_ob_tmp,hzthz,nobs,nanals)
hzthz = hzthz*nanals_m1_inv
hzthz_sv = hzthz
!------------------------------------------------------------------------------
print *, "[3] Calculate C and Gamma in Bishop et al. 2001 Eq.(11)"
!------------------------------------------------------------------------------
!
! ON OUTPUT hzthz IS C AND CONTAINS THE EIGENVECTORS
! ON OUTPUT egnv1 IS GAMMA AND CONTAINS THE EIGENVALUES
! CALCULATE THE EIGENVALUES AND VECTORS
call dsyev('V', 'L', nanals , hzthz, nanals, eignv1, workt, -1, info)
lwork = int(workt(1))
allocate (work(lwork))
call dsyev('V', 'L', nanals , hzthz, nanals, eignv1, work, lwork, info)
deallocate(work)
!
! REORDER THE EIGENVALUES AND EIGENVECTORS BECAUSE THEY ARE IN ASCENDING ORDER
do i = 1, nanals
eignv(i) = eignv1(nanals-i+1)
print *, "Eigenvalue ",i,eignv(i)
end do
print *, " "
!
! THE EIGENVECTORS ARE STORED AS COLUMN VECTORS
do i = 1, nanals
C(:,i) = hzthz(:,nanals-i+1)
end do
write(120) nobs,nanals
write(120) eignv
write(120) C
!
! CHECK EIGENVALUES AND VECTORS
sumv=0.
do k=1,nanals
do i=1,nanals
do j=1,nanals
trm=0.
if(i.eq.j) trm=eignv(k)
sumv=sumv+(hzthz_sv(i,j)-trm)*C(j,k)
enddo
enddo
print *, "Eigenvalue-vector check",k, sumv
enddo
!------------------------------------------------------------------------------
print *, "[4] Calculate inflation factor"
!------------------------------------------------------------------------------
!
! SNYDER INFLATION FACTOR ALGORITHM
if (infl_fac_WG03) then
print *, "USING WG03 VERSION OF INFLATION FACTOR ALGORITHM"
!
! IF (naccumt1 <= 0) USE PRESPECIFIED INFLATION FACTOR
if ( naccumt1 > 0 ) then
!
! CALCULATE TRACE BY SUMMING EIGENVECTORS
tracehpfht = 0.0
do i = 1, nanals-1
tracehpfht = tracehpfht + eignv(i)
end do
!
! CALCULATE OF NEGATIVE OF ENSEMBLE MEAN INNOVATION
! THEN SUM THE SQUARES AND NORMALIZE BY THE OBSERVATION ERROR VARIANCE
! THIS (d_tilda)T(d_tilda) WANG AND BISHOP (2003) EQ.(16)
obsinc(:) = ensmean_ob(:) - obs(:)
obsinc_sq = sum( obsinc(:) * obsinc(:))
squareinno = sum( obsinc(:) * obsinc(:) / oberrvar(:) )
write(99,rec=nout) squareinno
write(120) obsinc_sq, squareinno
!
! CALCULATE RUNNING MEAN OF SQUARED INNOVATIONS
print *, "squareinno_now", squareinno
if ( nout .eq. 1 ) then
squareinno_mn = squareinno
else
nmin = max( 1, nout - naccumt1 + 1 )
squareinno_mn = 0.0
do n = nmin, nout-1
read(99,rec=n) squareinno_old
print *, "square inno ", n, squareinno_old
read(109,rec=n) ainflat
print *, "ainflat_old ", n, ainflat
squareinno_mn = squareinno_mn + squareinno_old
end do
squareinno_mn = (squareinno_mn + squareinno)/real(nout - nmin + 1)
endif
!
! RATIO OF SQUARED INNOVATIONS MINUS OBS ERROR TO ENSEMBLE VARIANCE AT OBS LOCATIONS
! (ainflat_tmp) SHOULD CONVERGE TO ONE
!
! CALCULATE INITIAL ESTIMATE OF INFLATION FACTOR WANG and BISHOP (2003) EQ.(18)
ainflat_tmp = ( squareinno_mn - real(nobs) ) / tracehpfht
print *, "squareinno_mn ",squareinno_mn
print *, "nobs ",real(nobs)
print *, "trace ",tracehpfht
print *, "ainflat ",ainflat_tmp
if (ainflat_tmp .le. 0) then
ainflat_tmp = 1.
endif
!
! UPDATE INITIAL ESTIMATE FOLLOWING WANG AND BISHOP (2003) EQ.(17)
if ( nout .eq. 1 ) then
ainflat = 1
else
read(109,rec=nout-1) ainflat
endif
ainflat_mean = ainflat * ainflat_tmp
print *, "ainflat_old ",ainflat
print *, "ainflat_tmp ",ainflat_tmp
print *, "ainflat_mean ",ainflat_mean
write(120) ainflat_mean
close(120)
!
! SAVE THE FINAL SMOOTHED INFLATION FACTOR
write(109,rec=nout) ainflat_mean
else
!
! THIS IS FOR PRESPECIFIED INFLATION FACTOR
print *, 'Using prespecified inflation factor: ',tainflatinput
ainflat_mean = tainflatinput
end if
end if
!
! USE THE TRUNK INFLATION FACTOR ALGORITHM
if(infl_fac_TRNK) then
print *, "USING TRNK VERSION OF INFLATION FACTOR ALGORITHM"
!
! IF (naccumt1 <= 0) USE PRESPECIFIED INFLATION FACTOR
if ( naccumt1 > 0 ) then
!
! CALCULATE TRACE BY SUMMING EIGENVECTORS
tracehpfht = 0.0
do i = 1, nanals-1
tracehpfht = tracehpfht + eignv(i)
end do
!
! CALCULATE OF NEGATIVE OF ENSEMBLE MEAN INNOVATION
! THEN SUM THE SQUARES AND NORMALIZE BY THE OBSERVATION ERROR VARIANCE
! THIS (d_tilda)T(d_tilda) WANG AND BISHOP (2003) EQ.(16)
obsinc(:) = ensmean_ob(:) - obs(:)
obsinc_sq = sum( obsinc(:) * obsinc(:))
squareinno = sum( obsinc(:) * obsinc(:) / oberrvar(:) )
write(120) obsinc_sq, squareinno
!
! CALCULATE INITIAL ESTIMATE OF INFLATION FACTOR WANG and BISHOP (2003) EQ.(18)
! AND SAVE RAW INFLATION FACTOR
ainflat_tmp = ( squareinno - real(nobs) ) / tracehpfht
print *, "squareinno ",squareinno
print *, "nobs ",real(nobs)
print *, "trace ",tracehpfht
print *, "ainflat ",ainflat_tmp
if(ainflat_tmp.le.0.) then
ainflat_tmp=1.
endif
!
! UPDATE INITIAL ESTIMATE FOLLOWING WANG AND BISHOP (2003) EQ.(17)
if ( nout .eq. 1 ) then
ainflat = 1
else
read(109,rec=nout-1) ainflat
endif
ainflat_mean = ainflat * ainflat_tmp
print *, "ainflat_old ",ainflat
print *, "ainflat_tmp ",ainflat_tmp
print *, "ainflat_mn ",ainflat_mean
write(120) ainflat_mean
close(120)
!
! CALCULATE RUNNING MEAN INFLATION FACTOR
nmin = max( 1, nout - naccumt1 + 1 )
ainflat_mn = 0.0
do n = nmin, nout-1
read(109,rec=n) ainflat
print *, "ainflat_old ", n, ainflat
ainflat_mn = ainflat_mn + ainflat
end do
ainflat_mean = (ainflat_mn+ainflat_mean) / real( nout - nmin + 1 )
print *, "ainflat_mean ",ainflat_mean
write(109,rec=nout) ainflat_mean
else
!
! THIS IS FOR PRESPECIFIED INFLATION FACTOR
ainflat_mean = tainflatinput
end if
end if
!
! USE THE BOWLER INFLATION FACTOR ALGORITHM
if(infl_fac_BOWL) then
print *, "USING BOWL INFLATION FACTOR ALGORITHM"
!
! IF (naccumt1 <= 0) USE PRESPECIFIED INFLATION FACTOR
if ( naccumt1 > 0 ) then
!
! CALCULATE TRACE BY SUMMING EIGENVECTORS
tracehpfht = 0.0
do i = 1, nanals-1
tracehpfht = tracehpfht + eignv(i)
end do
write(119,rec=nout) tracehpfht
!
! CALCULATE OF NEGATIVE OF ENSEMBLE MEAN INNOVATION
! THEN SUM THE SQUARES AND NORMALIZE BY THE OBSERVATION ERROR VARIANCE
! THIS (d_tilda)T(d_tilda) WANG AND BISHOP (2003) EQ.(16)
obsinc(:) = ensmean_ob(:) - obs(:)
obsinc_sq = sum( obsinc(:) * obsinc(:))
squareinno = sum( obsinc(:) * obsinc(:) / oberrvar(:) )
write(120) obsinc_sq, squareinno
!
! CALCULATE INITIAL ESTIMATE OF INFLATION FACTOR WANG and BISHOP (2003) EQ.(18)
! AND SAVE RAW INFLATION FACTOR
ainflat_tmp = ( squareinno - real(nobs) )
!
! CALCULATE BOWLER'S GEOMETRIC MEAN ADJUSTMENT
if(nout.eq.1) then
tracehpfht_old=ainflat_tmp
else
read(119,rec=nout-1) tracehpfht_old
print *, "trace_old ", tracehpfht_old
endif
ainflat_tmp = sqrt(ainflat_tmp * tracehpfht_old)/tracehpfht
print *, "squareinno ",squareinno
print *, "nobs ",real(nobs)
print *, "trace ",tracehpfht
print *, "ainflat ",ainflat_tmp
if(ainflat_tmp.le.0.) then
ainflat_tmp=1.
endif
!
! UPDATE INITIAL ESTIMATE FOLLOWING WANG AND BISHOP (2003) EQ.(17)
if ( nout .eq. 1 ) then
ainflat = 1
else
read(109,rec=nout-1) ainflat
endif
ainflat_mean = ainflat * ainflat_tmp
print *, "ainflat_old ",ainflat
print *, "ainflat_tmp ",ainflat_tmp
print *, "ainflat_mn ",ainflat_mean
write(120) ainflat_mean
close(120)
write(109,rec=nout) ainflat_mean
else
!
! THIS IS FOR PRESPECIFIED INFLATION FACTOR
ainflat_mean = tainflatinput
end if
end if
!
! USE THE WANG ET. AL. 2007 INFLATION FACTOR ALGORITHM INCLUDING THE RHO FACTOR
if(infl_fac_WG07) then
print *, "USING WG07 INFLATION FACTOR ALGORITHM"
!
! IF (naccumt1 <= 0) USE PRESPECIFIED INFLATION FACTOR
if ( naccumt1 > 0 ) then
!
! CALCULATE TRACE BY SUMMING EIGENVECTORS
tracehpfht = 0.0
do i = 1, nanals-1
tracehpfht = tracehpfht + eignv(i)
end do
!
! CALCULATE OF NEGATIVE OF ENSEMBLE MEAN INNOVATION
! THEN SUM THE SQUARES AND NORMALIZE BY THE OBSERVATION ERROR VARIANCE
! THIS (d_tilda)T(d_tilda) WANG AND BISHOP (2003) EQ.(16)
obsinc(:) = ensmean_ob(:) - obs(:)
obsinc_sq = sum( obsinc(:) * obsinc(:))
squareinno = sum( obsinc(:) * obsinc(:) / oberrvar(:) )
write(99,rec=nout) squareinno
write(120) obsinc_sq, squareinno
!
! CALCULATE RUNNING MEAN OF SQUARED INNOVATIONS
print *, "squareinno_now ", squareinno
if (nout .eq. 1) then
squareinno_mn = squareinno
else
nmin = max( 1, nout - naccumt1 + 1 )
squareinno_mn = 0.0
do n = nmin, nout-1
read(99,rec=n) squareinno_old
print *, "square inno_old ", n, squareinno_old
read(109,rec=n) ainflat
print *, "ainflat_old ", n, ainflat
squareinno_mn = squareinno_mn + squareinno_old
enddo
squareinno_mn = (squareinno_mn + squareinno)/real(nout - nmin + 1)
endif
!
! BOUND PRE-INFLATION FACTOR BY ZERO
ainflat_tmp = (squareinno_mn - real(nobs)) / tracehpfht
print *, "squareinno_mn ",squareinno_mn
print *, "nobs ",real(nobs)
print *, "trace ",tracehpfht
print *, "ainflat ",ainflat_tmp
if(ainflat_tmp .le. 0.) then
ainflat_tmp=1.0
endif
!
! CALCULATE INFLATION FACTOR
if(nout .eq. 1) then
ainflat = 1.
else
read(109, rec=nout-1) ainflat
endif
ainflat_mean = ainflat * ainflat_tmp
print *, "ainflat_old ",ainflat
print *, "ainflat_tmp ",ainflat_tmp
print *, "ainflat_mean ",ainflat_mean
write(109, rec=nout) ainflat_mean
write(120) ainflat_mean
close(120)
else
!
! THIS IS FOR PRESPECIFIED INFLATION FACTOR
ainflat_mean = tainflatinput
write(109, rec=nout) ainflat_mean
end if
!
! CALCULATE THE RHO FACTOR FOLLOWING WANG ET. AL. 2007
if (naccumt2 > 0) then
!
! CALCULATE THE ENSEMBLE SUBSPACE EIGENVECTORS E=R(-1/2)HZCgamma^(-1/2) (SEE WANG ET. AL. 2007)
do i = 1, nanals
do j = 1, nanals-1
cgamma(i,j) = C(i,j)*sqrt(1.0/eignv(j))
enddo
enddo
!
! CALCULATE R(-1/2)HZ
do i = 1, nobs
enspert_ob_tmp(i,:) = ens_ob(i,:)/sqrt(oberrvar(i))/sqrt(real(nanals-1))
enddo
call da_matmulti(enspert_ob_tmp,cgamma,E,nobs,nanals-1,nanals)
!
! PROJECT NORMALIZED (DIVIDED BY OBERRSTDEV) INNOVATION VECTOR ONTO E
proj = 0.0
do i = 1, nanals-1
do k = 1, nobs
proj(i) = proj(i) + obsinc(k)/sqrt(oberrvar(k))*E(k,i)
enddo
enddo
!
! Get rho = (sum(proj*proj)-dim. of proj))/(normalized innovation^2-nobs)
! since nanals is relatively small, need follow wang et al (2007) for accumulation
! relative error = sqrt(2/(dim of proj)) = e.g., sqrt(2/207)=10%
! e.g., 50mem 2wks = sqrt(2/(14*49)) = 5% = e.g. 3wks sqrt(2/(21*49))=4%
!
! AVERAGE PROJ2
proj2 = sum(proj*proj)
write(89,rec=nout) proj2
if (nout .eq. 1) then
proj2sum = proj2
else
proj2sum = 0.0
nmin = max( 1, nout - naccumt2 + 1 )
do n = nmin, nout
read(89,rec=n) proj2sum1
print *, "proj2_old ",n, proj2sum1
proj2sum = proj2sum + proj2sum1/real(nout - nmin + 1)
enddo
print *, "proj2_mn ", proj2sum
endif
!
! AVERAGE SQUAREINNO
if (nout .eq. 1) then
squareinno_mn = squareinno
else
nmin = max( 1, nout - naccumt2 + 1 )
squareinno_mn = 0.0
do n = nmin, nout
read(99,rec=n) squareinno
print *, "squareinno_old ",n, squareinno
squareinno_mn = squareinno_mn + squareinno/real(nout - nmin + 1)
enddo
endif
print *, "squareinno_mn ",squareinno_mn
!
! CALCULATE RHO
rho = (proj2sum-(real(nanals-1)))/(squareinno_mn-real(nobs))
print *, "proj_mn ",proj2sum
print *, "nanals-1 ",nanals-1
print *, "squareinno_mn ",squareinno_mn
print *, "nobs ",nobs
print *, "rho ",rho
if (rho <= 0) then
rho = 1.0
endif
else
!
! THIS IS FOR PRE_SPECIFIED RHO FACTOR
rho = rhoinput
endif
endif
!
!------------------------------------------------------------------------------
print *, "[5] Calculate the grand transformation matrix"
!------------------------------------------------------------------------------
!
! CALCULATE THE TRANSFORMATION MATRIX. THIS FORM OF TRANSFORM MATRIX IS FROM
! WANG, ET. AL., (2007) EQ.(8). THE FOLLOWING FORM
! IS FROM BISHOP, ET. AL., (2001)
do i = 1, nanals
do j = 1, nanals-1
T(i,j) = C(i,j)*sqrt(1.0/(rho*eignv(j)+1.0))
enddo
enddo
do i = 1, nanals-1
do j = 1, nanals
CT(i,j) = C(j,i)
end do
end do
call da_matmulti(T,CT,TR,nanals,nanals,nanals-1)
!
! APPLY THE SQUARE ROOT OF THE INFLATION FACTOR
TR = sqrt(ainflat_mean) * TR
if(infl_fac_WG07) then
print *, "ENS"
print *, "MAX ENS ", maxval(ens)
print *, "MIN ENS ", minval(ens)
print *, "EIGEN/(EIGEN*RHO+1)"
print *, eignv/(eignv*rho+1.0)
sum_eignv = 0.
aftersum_eignv = 0.
do i=1,nanals
sum_eignv = sum_eignv + eignv(i)
aftersum_eignv = aftersum_eignv + eignv(i)/(eignv(i)*rho+1.0)
enddo
print *, "SUM_EIGNV=", sum_eignv
print *, "SUM_AFTEREIGNV=", aftersum_eignv
print *, "AVG_ENS=", sum(ens)/real(ndim*nanals)
print *, "STD_ENS=", sum(ens*ens)/float(ndim*nanals)
endif
!
!------------------------------------------------------------------------------
print *, "[6] Calculate the rescaled ETKF perturbations"
!------------------------------------------------------------------------------
!
! CALCULATE THE UPDATED PERTURBATIONS
call da_matmultiover(ens, TR, ndim, nanals)
if(infl_fac_WG07) then
print *, "ENS"
print *, "MAX ENS ", maxval(ens)
print *, "MIN ENS ", minval(ens)
print *, "AFTER AVG_ENS=", sum(ens)/real(ndim*nanals)
print *, "AFTER STD_ENS=", sum(ens*ens)/float(ndim*nanals)
endif
if (trace_use) call da_trace_entry("da_solve_etkf")
end subroutine da_solve_etkf
!
<A NAME='DA_SOLVE_LETKF'><A href='../../html_code/etkf/da_solve_etkf.inc.html#DA_SOLVE_LETKF' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_solve_letkf( ndim,nanals,nobs,ens,ens_ob,oberrvar,obs,nout,&,4
naccumt1,naccumt2,nstartaccum1,nstartaccum2,tainflatinput,rhoinput, &
infl_fac_file,eigen_val_file,inno2_val_file,proj2_val_file,infl_fac_TRNK, &
infl_fac_WG03,infl_fac_WG07,infl_fac_BOWL,yo_lat,yo_lon,yo_prs,xf_lat,xf_lon, &
ijk_idx, apl_idx, apl_ndim, napl1, napl2, nxs, nys, nzs, nv,infl_let_file)
!-----------------------------------------------------------------------
! Purpose: LETKF perturbation updates
! Arthur P. Mizzi, February 2011 Modified the da_solve_etkf routine to become
! LETKF. Under development.
!
!1) nanals, ensemble size
!2) ndim, dimension of the perturbation vector that needs to be updated
!3) nobs, number of observations assimilated
!4) ens, array of perturbations before (Xf') and after (Xa') ETKF update
!5) ens_ob, array of HXf
!6) oberrvar, observation error variance, listed in the same sequence as HXf
!7) obs, observations assimilated
!8) naccumt1, number of previous cycles immediately before the current cycle,
! which is needed for calculating adaptive inflation factor.
! naccumt1 < 0 for pre-specified inflation
!9) naccumt2, number of previous cycles immediately before the current cycle,
! which is needed for calculating the rho factor in the latest version of ETKF.
! naccumt2 < 0 for using the older version of the ETKF.
! naccumt2 = 0 for using pre-specified rho factor
!10) nstartaccum1, the cycle from which the accumulation of previous naccumt1 cycles in 8) starts
!11) nstartaccum2, the cycle from which the accumulation of previous naccumt2 in 9) starts
!12) tainflatinput, pre-specified inflation, if not using adaptive inflation
!13) rhoinput, pre-specified rho factor, if not using adaptively determined rho factor
!14) nout, record number for output squared innovations and the forecast error variance
!---------------projected onto ensemble subspace, which is related to 8) and 9)
implicit none
integer, intent(in) :: nanals,ndim,nobs
real, intent(inout), dimension(ndim,nanals) :: ens
real, intent(inout), dimension(nobs,nanals) :: ens_ob
real, intent(in), dimension(nobs) :: oberrvar
real, intent(in), dimension(nobs) :: obs
real, intent(in) :: tainflatinput,rhoinput
integer, intent(in) :: nout,naccumt1,naccumt2,nstartaccum1,nstartaccum2
integer, intent(in) :: napl1, napl2, nxs, nys, nzs, nv
integer, intent(in), dimension(nv,nxs,nys,nzs) :: ijk_idx
integer, intent(in), dimension(napl1,napl2) :: apl_idx
integer, intent(in), dimension(napl1) :: apl_ndim
integer :: n, nstart, ndivd ! Loop counters.
integer :: nmin ! Minimum nout value to use.
real :: nanals_inv ! 1 / nanals.
real :: nanals_m1_inv ! 1 / (nanals-1).
real :: ainflat_mean ! Rolling mean inflation factor.
real :: ainflat_mean_old ! Rolling mean inflation factor read in from file.
real :: ainflat_tmp ! Instantaneous ratio of sqd innovations minus
! obs error to forecast variance at obs locations.
real, dimension(nobs) :: ensmean_ob
real, dimension(nobs) :: obsinc
integer :: ij, nanal, i, j, k ! Loop counters.
integer :: info, lwork
real, allocatable, dimension(:) :: work
real, dimension(1) :: workt
real, dimension(nanals) :: eignv, eignv1
real, dimension(nanals, nanals) :: hzthz, C, TR
real, dimension(nanals-1, nanals) :: CT
real, dimension(nanals, nanals-1) :: T
real, dimension(nanals, nanals-1) :: cgamma
real, dimension(nobs,nanals-1) :: E
real, dimension(nobs,nanals) :: enspert_ob_tmp
character (len=150) :: filename_letkf
character (len=150) :: filename_factor
character (len=150) :: filename_eigenv
character (len=150) :: filename_inno2
character (len=150) :: filename_proj2
real, dimension(nanals-1) :: proj
real :: tracehpfht
real :: ta, tb
real :: sum_eignv
real :: aftersum_eignv
real :: ainflat, ainflat_tmp_sav
real :: tainflat, rho
real :: proj2, proj2sum
real :: proj2sum1
real :: squareinnosum, squareinnosum1
real :: squareinno
character (len=150) :: infl_fac_file,eigen_val_file,inno2_val_file,proj2_val_file
character (len=150) :: infl_let_file
logical :: infl_fac_TRNK,infl_fac_WG03,infl_fac_WG07, infl_fac_BOWL
real, intent(in), dimension(nobs) :: yo_lat,yo_lon,yo_prs
real, intent(in), dimension(ndim) :: xf_lat,xf_lon
real, dimension(ndim,nanals) :: ens_p
real, dimension(ndim,naccumt1) :: ainflat_v, ainflat_u, squareinno_v
real :: z_rad, l_rad, ref_lat, ref_lon, ref_prs, squareinno_mn, squareinno_sav
real :: summ, rearth, pi, deg2rad,d_lat, d_lon, ainflat_tmp_mn, ainflat_snyd
real :: ainflat_apm,zinno_mn,zinno_vr,zinfl_apm_mn,zinfl_apm_vr,zinfl_snyd_mn
real :: zinfl_snyd_vr
integer :: ij_cnt, ijp, nout_p, nout_pp, iprd, icol, idx, idxx, idy, ii, jj
logical :: avg_wang, avg_snyd, avg_apm
!
! BEGIN EXECUTABLE CODE
if (trace_use) call da_trace_entry("da_solve_etkf")
!
! SWITCH TO CHOOSE INFLATION FACTOR ALGORITHM
print *, "Enter LETKF with nout = ",nout
print *, "infl_fac_TRNK ", infl_fac_TRNK
print *, "infl_fac_WG03 ", infl_fac_WG03
print *, "infl_fac_WG07 ", infl_fac_WG07
print *, "infl_fac_BOWL ", infl_fac_BOWL
rho=1.
infl_fac_WG03 = .true.
!
! PRINT INPUT OBS
! do i = 1, nobs
! print *, i,obs(i),oberrvar(i),ens_ob(i,1), &
! ens_ob(i,5),ens_ob(i,10)
! end do
!
! OPEN THE INFLATION FACTOR HISTORY FILE
!
filename_factor = trim(infl_fac_file)
open(109,form="unformatted",file=filename_factor,access="direct",recl=8)
!
filename_letkf = trim(infl_let_file)
open(110,form="unformatted",file=filename_letkf)
!
filename_eigenv = trim(eigen_val_file)
open(119,form="unformatted",file=filename_eigenv,access="direct",recl=8)
!
filename_inno2 = trim(inno2_val_file)
open(99,form="unformatted",file=filename_inno2)
!
filename_proj2 = trim(proj2_val_file)
open(89,form="unformatted",file=filename_proj2,access="direct",recl=8)
!
!------------------------------------------------------------------------------
print *, "[1] Compute mean(H(xf)) and H(xf)'"
!------------------------------------------------------------------------------
!
! CALCULATE THE ENSEMBLE MEAN AND CREATE ENSEMBLE PERTURBATIONS
! FOR THE OB.ETKF.exxx DATA
nanals_inv = 1.0 / real(nanals)
nanals_m1_inv = 1.0 / real(nanals-1)
do ij = 1, nobs
ensmean_ob(ij) = sum(ens_ob(ij,:)) * nanals_inv
enddo
do nanal = 1, nanals
ens_ob(:,nanal) = ens_ob(:,nanal) - ensmean_ob(:)
enddo
!
!-----------------------------------------------------------------------------
print *,"[2] Loop over locations and calculate Transform(HZ)(HZ)"
print *," in Bishop et al. 2001 Eq.(11)"
!------------------------------------------------------------------------------
!
z_rad = 5000000.
pi = 3.1415927
rearth = 6378388.
deg2rad = pi/180.
d_lat=2.*pi*rearth/360.
idy=0
do ii=1,nxs-1
do jj=1,nys-1
idy=idy+1
ij=ijk_idx(5,ii,jj,1)
ref_lat = xf_lat(ij)
ref_lon = xf_lon(ij)
if(ref_lon .lt. 0.) ref_lon=360.+ref_lon
!
! NORMALIZE H(xf)' BY SQRT(R) TO GET Tilda(H) THEN CALCULATE
! (HZ)T(HZ) IN EQN(11) AND NORMALIZE TO GET THE ENSEMBLE FORECAST
! COVARIANCE MATRIX IN hzths(nanals,nanals)
ij_cnt=0
obsinc(:)=0.
enspert_ob_tmp(:,:)=0.
! print *, 'APM: ref_lat, raf_lon ',ref_lat, ref_lon
do ijp = 1, nobs
d_lon=2.*pi*rearth*cos((ref_lat+yo_lat(ijp))/2.*deg2rad)/360.
l_rad = sqrt(((yo_lat(ijp)-ref_lat)*d_lat)**2 + &
((yo_lon(ijp)-ref_lon)*d_lon)**2)
! print *, 'APM: l_rad, z_rad ',l_rad, z_rad
if(l_rad .le. z_rad) then
ij_cnt=ij_cnt+1
enspert_ob_tmp(ij_cnt,:) = ens_ob(ijp,:)/sqrt(oberrvar(ijp))
obsinc(ij_cnt)=(ensmean_ob(ijp)-obs(ijp))/sqrt(oberrvar(ijp))
endif
enddo
print *, 'APM: local nobs ',ij_cnt
ij=idy
hzthz=0.
do i=1,nanals
do j=1,nanals
do ijp=1,ij_cnt
hzthz(i,j)=hzthz(i,j)+enspert_ob_tmp(ijp,i)*enspert_ob_tmp(ijp,j)
enddo
enddo
enddo
hzthz = hzthz*nanals_m1_inv
! do i=1,nanals
! print *,i,j,hzthz(i,1),hzthz(i,2),hzthz(i,3),hzthz(i,4),hzthz(i,5)
! enddo
!
!------------------------------------------------------------------------------
! print *, "[3] Calculate local C and Gamma in Bishop et al. 2001 Eq.(11)"
!------------------------------------------------------------------------------
!
! ON OUTPUT hzthz IS C AND CONTAINS THE EIGENVECTORS
! ON OUTPUT egnv1 IS GAMMA AND CONTAINS THE EIGENVALUES
! CALCULATE THE EIGENVALUES AND VECTORS
call dsyev('V', 'L', nanals , hzthz, nanals, eignv1, workt, -1, info)
lwork = int(workt(1))
allocate (work(lwork))
call dsyev('V', 'L', nanals , hzthz, nanals, eignv1, work, lwork, info)
deallocate(work)
!
! REORDER THE EIGENVALUES AND EIGENVECTORS BECAUSE THEY ARE IN ASCENDING ORDER
do i = 1, nanals
eignv(i) = eignv1(nanals-i+1)
enddo
! print *, "eigen values ",eignv
!
! THE EIGENVECTORS ARE STORED AS COLUMN VECTORS
do i = 1, nanals
C(:,i) = hzthz(:,nanals-i+1)
enddo
! do i=1,nanals
! print *,i,j,C(i,1),C(i,2),C(i,3),C(i,4),C(i,5)
! enddo
!
!------------------------------------------------------------------------------
print *, "[4] Calculate inflation factor"
!------------------------------------------------------------------------------
!
! WG03 INFLATION FACTOR ALGORITHM
if (infl_fac_WG03) then
!
! IF (naccumt1 <= 0) USE PRESPECIFIED INFLATION FACTOR
if ( naccumt1 > 0 ) then
!
! CALCULATE TRACE BY SUMMING EIGENVECTORS
tracehpfht = 0.0
do i = 1, nanals-1
tracehpfht = tracehpfht + eignv(i)
enddo
!
! SET FLAG FOR AVERAGING SCHEME
avg_wang=.false.
avg_apm=.false.
avg_snyd=.true.
!
! CALCULATE OF NEGATIVE OF ENSEMBLE MEAN INNOVATION
! THEN SUM THE SQUARES AND NORMALIZE BY THE OBSERVATION ERROR VARIANCE
! THIS (d_tilda)T(d_tilda) WANG AND BISHOP (2003) EQ.(16)
squareinno = sum(obsinc(:) * obsinc(:))
! print *, 'local squareinno ',squareinno
!
! Wang averaging
if (ii.eq.1 .and. jj.eq.1 .and. nout.ne.1) then
read(99) squareinno_v
endif
squareinno_v(ij,1) = squareinno
nout_p = nout
if (nout .gt. naccumt1) nout_p = naccumt1
squareinno_mn=0.
do n = 1, nout_p
squareinno_mn = squareinno_mn + squareinno_v(ij,n)/real(nout_p)
end do
squareinno_sav=squareinno
if(avg_wang) then
squareinno=squareinno_mn
endif
nout_pp=nout_p+1
if (nout_pp .gt.naccumt1) nout_pp = naccumt1
do n = nout_pp, 2, -1
squareinno_v(ij,n)=squareinno_v(ij,n-1)
enddo
squareinno_v(ij,1)=0.
if (ii.eq.nxs-1 .and. jj.eq.nys-1) then
rewind(99)
write(99) squareinno_v
endif
! print *, "squareinno_sv ",squareinno_sav
! print *, "squareinno_mn ",squareinno_mn
!
! BEGIN SNYDER CHANGES
! RATIO OF SQUARED INNOVATIONS MINUS OBS ERROR TO ENSEMBLE VARIANCE AT OBS LOCATIONS
! (ainflat_tmp) SHOULD CONVERGE TO ONE
!
! CALCULATE INITIAL ESTIMATE OF INFLATION FACTOR WANG and BISHOP (2003) EQ.(18)
ainflat_tmp = ( squareinno - real(ij_cnt) ) / tracehpfht
if (ainflat_tmp .le. 0.) then
ainflat_tmp = 1.
endif
!
! UPDATE INITIAL ESTIMATE FOLLOWING WANG AND BISHOP (2003) EQ.(17)
! ainflat(ij,i) is a inflation factor history storage stack: ij is
! the location index, i is the history index (i=1 is top of stack
! i=naccumt1 is bottom of stack).
if ( nout .eq. 1 ) then
ainflat_v(ij,2) = 1.
else if (ii.eq.1 .and. jj.eq.1) then
read(110) ainflat_v
read(110) ainflat_u
endif
!
! APM averaging
! ainflat_u stores the pre-infaltion factor
! ainflat_v stores the final inflation factor
!
ainflat_u(ij,1) = ainflat_tmp
ainflat_tmp_mn=0.
nout_p = nout
if (nout .gt. naccumt1) nout_p = naccumt1
do n = 1, nout_p
ainflat_tmp_mn = ainflat_tmp_mn + ainflat_u(ij,n)/real(nout_p)
end do
ainflat_apm = ainflat_v(ij,2) * ainflat_tmp_mn
!
! Snyder averaging
ainflat_v(ij,1) = ainflat_v(ij,2) * ainflat_tmp
ainflat_snyd = 0.
nout_p = nout
if (nout .gt. naccumt1) nout_p = naccumt1
do n = 1, nout_p
ainflat_snyd = ainflat_snyd + ainflat_v(ij,n)/real(nout_p)
end do
!
! SAVE THE FINAL SMOOTHED INFLATION FACTOR
print *, "pre-inflation factor ", ainflat_tmp
print *, "inflation factor Snyder AVG", ainflat_snyd
! print *, "inflation factor APM AVG", ainflat_apm
if(avg_wang) then
ainflat_mean=ainflat_apm
else if(avg_snyd) then
ainflat_mean=ainflat_snyd
endif
!
ainflat_v(ij,1) = ainflat_mean
nout_pp=nout_p+1
if (nout_pp .gt.naccumt1) nout_pp = naccumt1
do n = nout_pp, 2, -1
ainflat_v(ij,n)=ainflat_v(ij,n-1)
ainflat_u(ij,n)=ainflat_u(ij,n-1)
enddo
ainflat_v(ij,1)=0.
ainflat_u(ij,1)=0.
if (ii.eq.nxs-1 .and. jj.eq.nys-1) then
rewind(110)
write(110) ainflat_v
write(110) ainflat_u
endif
!
! END SNYDER CHANGES
else
!
! THIS IS FOR PRESPECIFIED INFLATION FACTOR
! print *, 'Using prespecified inflation factor: ',tainflatinput
ainflat_mean = tainflatinput
endif
endif
!
! USE THE TRNK INFLATION FACTOR ALGORITHM
if(infl_fac_TRNK) then
print *, "USING TRNK VERSION OF INFLATION FACTOR ALGORITHM"
print *, "LETKF - only works with WG03 inflation factor"
stop
!
! IF (naccumt1 <= 0) USE PRESPECIFIED INFLATION FACTOR
if ( naccumt1 > 0 ) then
!
! CALCULATE TRACE BY SUMMING EIGENVECTORS
tracehpfht = 0.0
do i = 1, nanals-1
tracehpfht = tracehpfht + eignv(i)
enddo
!
! CALCULATE OF NEGATIVE OF ENSEMBLE MEAN INNOVATION
! THEN SUM THE SQUARES AND NORMALIZE BY THE OBSERVATION ERROR VARIANCE
! THIS (d_tilda)T(d_tilda) WANG AND BISHOP (2003) EQ.(16)
obsinc(:) = ensmean_ob(:) - obs(:)
squareinno = sum( obsinc(:) * obsinc(:) / oberrvar(:) )
print *, 'EIGENVALUES ',eignv
print *, 'TRACE ',tracehpfht
print *, 'SUM NORMALIZED SQUARED INNOV ',squareinno
print *, 'NUMBER OBS ',nobs
print *, 'EIGENVECTORS '
do i=1,nanals
write(unit=stdout,fmt='(10(g10.4,1x))') (C(i,n),n=1,nanals)
enddo
print *, ' '
!
! CALCULATE INITIAL ESTIMATE OF INFLATION FACTOR WANG and BISHOP (2003) EQ.(18)
! AND SAVE RAW INFLATION FACTOR
ainflat = ( squareinno - real(nobs) ) / tracehpfht
print *, 'RAW INFL FAC ',ainflat
if(ainflat.le.0.) then
ainflat=1.
endif
print *, 'CURRENT INFL FAC ',ainflat
write(109,rec=nout) ainflat
!
! CALCULATE RUNNING MEAN INFLATION FACTOR
nmin = max( 1, nout - naccumt1 + 1 )
ainflat_mean = 0.0
do n = nmin, nout
read(109,rec=n) ainflat
ainflat_mean = ainflat_mean + ainflat
end do
ainflat_mean = ainflat_mean / real( nout - nmin + 1 )
write (unit=stdout,fmt='(/a,f15.5)') " Current Inflation factor = ", ainflat
write (unit=stdout,fmt='(a,f15.5)') " Rolling mean inflation factor = ", ainflat_mean
else
! THIS IS FOR PRESPECIFIED INFLATION FACTOR
ainflat_mean = tainflatinput
end if
end if
!
! USE THE WG07 INFLATION FACTOR ALGORITHM INCLUDING THE RHO FACTOR
if(infl_fac_WG07) then
print *, "USING WG07 INFLATION FACTOR ALGORITHM"
print *, "LETKF - only works with WG03 inflation factor"
stop
!
! CALCULATE TRACE BY SUMMING EIGENVECTORS
tracehpfht = 0.0
do i = 1, nanals-1
tracehpfht = tracehpfht + eignv(i)
end do
!
! CALCULATE OF NEGATIVE OF ENSEMBLE MEAN INNOVATION
! THEN SUM THE SQUARES AND NORMALIZE BY THE OBSERVATION ERROR VARIANCE
! THIS (d_tilda)T(d_tilda) WANG AND BISHOP (2003) EQ.(16)
obsinc(:) = ensmean_ob(:) - obs(:)
squareinno = sum( obsinc(:) * obsinc(:) / oberrvar(:) )
write(99,rec=nout) squareinno
print *, 'MAX(obsinc) ', maxval(obsinc)
print *, 'MIN(obsinc) ', minval(obsinc)
print *, 'SQUAREINNO WRITE, NOUT ',squareinno, nout
print *, 'EIGENVALUES ',eignv
print *, 'TRACE ',tracehpfht
print *, 'NUMBER OBS ',nobs
print *, 'EIGENVECTORS '
do i=1,nanals
write(unit=stdout,fmt='(10(g10.4,1x))') (C(i,n),n=1,nanals)
enddo
!
! IF (naccumt1 <= 0) USE PRESPECIFIED INFLATION FACTOR
if ( naccumt1 > 0 ) then
squareinnosum = squareinno
if(nout >= nstartaccum1) then
squareinnosum = 0.0
nstart=max(1,nout-naccumt1+1)
ndivd=nout-nstart+1
do n = nstart, nout
read(99,rec=n) squareinnosum1
print *, 'APM: squareinno archive ',n, squareinnosum1
squareinnosum = squareinnosum + squareinnosum1/real(ndivd)
enddo
endif
!
! CALCULATE INFLATION FACTOR AND BOUND BY ZERO
ainflat = (squareinnosum - real(nobs)) / tracehpfht
if(ainflat .le. 0.) then
ainflat=1.0
endif
!
! CALCULATE RUNNING MEAN INFLATION FACTOR
if(nout <= nstartaccum1) then
ainflat_mean = ainflat
else
read(109, rec=nout-1) ainflat_mean
ainflat_mean=ainflat_mean*ainflat
endif
write(109, rec=nout) ainflat_mean
print *, 'CURRENT INFLATION FACTOR ',ainflat
print *, 'FINAL INFLATION FACTOR ',ainflat_mean
else
!
! THIS IS FOR PRESPECIFIED INFLATION FACTOR
ainflat_mean = tainflatinput
write(109, rec=nout) ainflat_mean
end if
!
! CALCULATE THE RHO FACTOR FOLLOWING WANG ET. AL. 2007
if (naccumt2 > 0) then
!
! CALCULATE THE ENSEMBLE SUBSPACE EIGENVECTORS E=R(-1/2)HZCgamma^(-1/2) (SEE WANG ET. AL. 2007)
do i = 1, nanals
do j = 1, nanals-1
cgamma(i,j) = C(i,j)*sqrt(1.0/eignv(j))
enddo
enddo
!
! CALCULATE R(-1/2)HZ
do i = 1, nobs
enspert_ob_tmp(i,:) = ens_ob(i,:)/sqrt(oberrvar(i))/sqrt(real(nanals-1))
enddo
call da_matmulti(enspert_ob_tmp,cgamma,E,nobs,nanals-1,nanals)
!
! PROJECT NORMALIZED (DIVIDED BY OBERRSTDEV) INNOVATION VECTOR ONTO E
proj = 0.0
do i = 1, nanals-1
do k = 1, nobs
proj(i) = proj(i) + obsinc(k)/sqrt(oberrvar(k))*E(k,i)
enddo
enddo
!
! get rho = (sum(proj*proj)-dim. of proj))/(normalized innovation^2-nobs)
! since nanals is relatively small, need follow wang et al (2007) for accumulation
! relative error = sqrt(2/(dim of proj)) = e.g., sqrt(2/207)=10%
! e.g., 50mem 2wks = sqrt(2/(14*49)) = 5% = e.g. 3wks sqrt(2/(21*49))=4%
!
! AVERAGE PROJ2
proj2 = sum(proj*proj)
write(89,rec=nout) proj2
proj2sum = proj2
if (nout >= nstartaccum2) then
proj2sum = 0.0
nstart = max(1,nout-naccumt2+1)
ndivd = nout-nstart+1
do n = nstart, nout
read(89,rec=n) proj2sum1
proj2sum = proj2sum + proj2sum1/real(ndivd)
enddo
endif
!
! AVERAGE SQUAREINNO
squareinnosum = squareinno
if (nout >= nstartaccum2) then
squareinnosum = 0.0
nstart = max(1,nout-naccumt2+1)
ndivd = nout-nstart+1
print *, 'NSTART, NOUT, NACCUMT1 ',nstart,nout,naccumt1
do n = nstart, nout
read(99,rec=n) squareinnosum1
squareinnosum = squareinnosum + squareinnosum1/real(ndivd)
print *, "SQUAREINNOSUM, SQUAREINNOSUM1 = ", squareinnosum,squareinnosum1
enddo
endif
!
! CALCULATE RHO
rho = (proj2sum-(real(nanals-1)))/(squareinnosum-real(nobs))
if (rho <= 0) then
rho = 1.0
endif
endif
if (naccumt2 < 0) then
!
! THIS IS FOR WANG AND BISHOP 2003 ORIGINAL FORMULATION
rho = 1.0
endif
if (naccumt2 == 0) then
!
! THIS IS FOR PRE_SPECIFIED RHO FACTOR
rho = rhoinput
endif
print *, "rho = ", rho
print *, "proj2sum=", proj2sum
endif
!
!------------------------------------------------------------------------------
! print *, "[5] Calculate the grand transformation matrix"
!------------------------------------------------------------------------------
!
! CALCULATE THE TRANSFORMATION MATRIX. THIS FORM OF TRANSFORM MATRIX IS FROM
! WANG, ET. AL., (2007) EQ.(8). THE FOLLOWING FORM
! IS FROM BISHOP, ET. AL., (2001)
do i = 1, nanals
do j = 1, nanals-1
T(i,j) = C(i,j)*sqrt(1.0/(rho*eignv(j)+1.0))
enddo
enddo
do i = 1, nanals-1
do j = 1, nanals
CT(i,j) = C(j,i)
end do
end do
call da_matmulti(T,CT,TR,nanals,nanals,nanals-1)
!
! APPLY THE SQUARE ROOT OF THE INFLATION FACTOR
TR = sqrt(ainflat_mean) * TR
if(infl_fac_WG07) then
print *, "ENS"
print *, "MAX ENS ", maxval(ens)
print *, "MIN ENS ", minval(ens)
print *, "EIGEN/(EIGEN*RHO+1)"
print *, eignv/(eignv*rho+1.0)
sum_eignv = 0.
aftersum_eignv = 0.
do i=1,nanals
sum_eignv = sum_eignv + eignv(i)
aftersum_eignv = aftersum_eignv + eignv(i)/(eignv(i)*rho+1.0)
enddo
print *, "SUM_EIGNV=", sum_eignv
print *, "SUM_AFTEREIGNV=", aftersum_eignv
print *, "AVG_ENS=", sum(ens)/real(ndim*nanals)
print *, "STD_ENS=", sum(ens*ens)/float(ndim*nanals)
endif
!
!------------------------------------------------------------------------------
! print *, "[6] Calculate the rescaled ETKF perturbations for local row only"
!------------------------------------------------------------------------------
!
! CALCULATE THE UPDATED PERTURBATIONS AND SAVE IN TEMPORARY ARRAY
! FOR GREATER EFFICIENCY IN THE LETKF APPLY THE LOCALIZED TRANSFORM
! TO ALL LIKE LOCATIONS
!
do idx=1,apl_ndim(ij)
idxx=apl_idx(ij,idx)
do icol = 1,nanals
summ = 0.
do iprd = 1,nanals
summ = summ + ens(idxx,iprd)*TR(iprd,icol)
enddo
ens_p(idxx,icol) = summ
enddo
enddo
if(infl_fac_WG07) then
print *, "ENS"
print *, "MAX ENS ", maxval(ens)
print *, "MIN ENS ", minval(ens)
print *, "AFTER AVG_ENS=", sum(ens)/real(ndim*nanals)
print *, "AFTER STD_ENS=", sum(ens*ens)/float(ndim*nanals)
endif
enddo
enddo
!
! CALCULATE SPATIAL MEAN AND VARIANCE OF INFLATION FACGTOR PARAMETERS
zinno_mn=0.
zinno_vr=0.
zinfl_apm_mn=0.
zinfl_apm_vr=0.
zinfl_snyd_mn=0.
zinfl_snyd_vr=0.
ij=0
do ii=1, nxs-1
do jj=1, nys-1
ij=ij+1
zinno_mn = zinno_mn + squareinno_v(ij,2)/float((nxs-1)*(nys-1))
zinfl_apm_mn = zinfl_apm_mn + ainflat_u(ij,2)/float((nxs-1)*(nys-1))
zinfl_snyd_mn = zinfl_snyd_mn + ainflat_v(ij,2)/float((nxs-1)*(nys-1))
enddo
!
! SAVE THE SPATIAL DOMAIN AVERAGE INFLATION FACTOR FOR PLOTTING PURPOSES
write(109,rec=nout) zinfl_snyd_mn
!
enddo
ij=0
do ii=1, nxs-1
do jj=1, nys-1
ij=ij+1
zinno_vr = zinno_vr + (squareinno_v(ij,2)-zinno_mn)**2/ &
float((nxs-1)*(nys-1)-1)
zinfl_apm_vr = zinfl_apm_vr + (ainflat_u(ij,2)-zinfl_apm_mn)**2/ &
float((nxs-1)*(nys-1)-1)
zinfl_snyd_vr = zinfl_snyd_vr + (ainflat_v(ij,2)-zinfl_snyd_mn)**2/ &
float((nxs-1)*(nys-1)-1)
enddo
enddo
print *, 'squareinno: mn, vr ',zinno_mn,zinno_vr
print *, 'ainflat_apm: mn, vr ',zinfl_apm_mn,zinfl_apm_vr
print *, 'ainflat_snyd: mn, vr ',zinfl_snyd_mn,zinfl_snyd_vr
!
! PASS UPDATED PERTURBATIONS TO WORKING ARRAY
ens(:,:)=ens_p(:,:)
!
if (trace_use) call da_trace_entry("da_solve_etkf")
end subroutine da_solve_letkf
!
<A NAME='RAND_FILTER'><A href='../../html_code/etkf/da_solve_etkf.inc.html#RAND_FILTER' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine rand_filter(y_in,y_sig_in,yo_in,yo_obs_typ_in,yo_typ_in,yo_subtyp_in, &
yo_lat_in,yo_lon_in,yo_prs_in,yo_tim_in,y_ot,y_sig_ot,yo_ot, &
yo_obs_typ_ot,yo_typ_ot,yo_subtyp_ot,yo_lat_ot,yo_lon_ot, &
yo_prs_ot,yo_tim_ot,idx_in,nobs_in,nobs_ot,nobs_flt,num_mems,iseed)
implicit none
!
! DEFINE INCOMMING/OUTGOING VARIABLES
integer,intent(in) :: nobs_in,nobs_flt,nobs_ot,num_mems,iseed
integer,intent(in),dimension(nobs_in,num_mems) :: idx_in
real,intent(in),dimension(nobs_in,num_mems) :: y_in, y_sig_in, yo_in, &
yo_typ_in, yo_subtyp_in, yo_lat_in,yo_lon_in, yo_prs_in, yo_tim_in
real,intent(out),dimension(nobs_ot,num_mems) :: y_ot
real,intent(out),dimension(nobs_ot) :: y_sig_ot,yo_ot,yo_typ_ot, &
yo_subtyp_ot,yo_lat_ot,yo_lon_ot,yo_prs_ot, yo_tim_ot
character(len=10),intent(in),dimension(nobs_in,num_mems) :: yo_obs_typ_in
character(len=10),intent(out),dimension(nobs_ot) :: yo_obs_typ_ot
!
! DEFINE WORK VARIABLES
integer, allocatable, dimension(:) :: idx_flt,idx_rnd,iput,iget
integer :: isize,inum,istr,idx,i,ii
real :: znum,zscl
!
! SET PARAMETERS
zscl=1000000.
!
! SETUP MAPPING BETWEEN INCOMMING AND PRE_FILTERED INDEXES
idx=0
allocate(idx_flt(nobs_flt))
do i=1,nobs_in
do ii=1,num_mems
if(idx_in(i,ii) .ne. -999) then
idx=idx+1
idx_flt(idx)=i
go to 1010
endif
enddo
1010 continue
enddo
! print *, ' '
! do i=1,nobs_flt
! print *, 'idx flt ',i,idx_flt(i)
! enddo
!
! SET THE RANDOM NUMBER SEED
call random_seed(SIZE=isize)
allocate (iput(isize))
do i=1,isize
iput(i)=iseed + (i-1)
enddo
call random_seed(PUT=iput)
!
! GET NOBS_FLT RANDOM NUMBERS BETWEEN 1 AND NOBS_FLT
istr=0
allocate(idx_rnd(nobs_ot))
do i=1,nobs_ot
1000 continue
call random_number(znum)
inum=ifix(zscl*znum)
if(inum.ge.1 .and. inum.le.nobs_flt) then
if(istr.eq. 0) then
istr=istr+1
idx_rnd(i)=inum
else
do ii=1,istr
if(idx_rnd(ii) .ne. inum) then
if(ii .eq. istr) then
istr=istr+1
idx_rnd(i)=inum
endif
else
go to 1000
endif
enddo
endif
else
go to 1000
endif
enddo
! print *, ' '
! do i=1,nobs_ot
! print *, 'idx rnd ',i,idx_rnd(i),idx_flt(idx_rnd(i))
! enddo
!
! ASSIGN THE RANDOMLY FILTERED OBS TO THE OUTGOING ARRAYS
! print *, ' '
do i=1,nobs_ot
y_sig_ot(i)=y_sig_in(idx_flt(idx_rnd(i)),1)
yo_ot(i)=yo_in(idx_flt(idx_rnd(i)),1)
yo_obs_typ_ot(i)=yo_obs_typ_in(idx_flt(idx_rnd(i)),1)
yo_typ_ot(i)=yo_typ_in(idx_flt(idx_rnd(i)),1)
yo_subtyp_ot(i)=yo_subtyp_in(idx_flt(idx_rnd(i)),1)
yo_lat_ot(i)=yo_lat_in(idx_flt(idx_rnd(i)),1)
yo_lon_ot(i)=yo_lon_in(idx_flt(idx_rnd(i)),1)
yo_prs_ot(i)=yo_prs_in(idx_flt(idx_rnd(i)),1)
yo_tim_ot(i)=yo_tim_in(idx_flt(idx_rnd(i)),1)
do ii=1,num_mems
y_ot(i,ii)=y_in(idx_flt(idx_rnd(i)),ii)
enddo
enddo
return
end subroutine rand_filter