da_solve_etkf.inc
References to this file elsewhere.
1 subroutine da_solve_etkf( ndim,nanals,nobs,ens,ens_ob,oberrvar,obs,nout,&
2 naccumt1,naccumt2,nstartaccum1,nstartaccum2,tainflatinput,rhoinput)
3
4 !-----------------------------------------------------------------------
5 ! Purpose: ETKF perturbation updates
6 ! Xuguang Wang, Jan. 2006
7 ! Dale Barker. January 2007. Implicit none, Convert to f90, and enclose within gen_be_etkf.f90 driver within WRF.
8 ! Dale Barker. Also modify inflation factor method (nstartaccum1,nstartaccum2 redundant).
9 !
10 ! references:
11 ! Bishop et al 2001 MWR,
12 ! Wang and Bishop 2003 MWR,
13 ! Wang et al. 2004, MWR
14 ! Wang et al. 2006, MWR
15 !
16 !1) nanals, ensemble size
17 !2) ndim, dimension of the perturbation vector that needs to be updated
18 !3) nobs, number of observations assmiilated
19 !4) ens, array of perturbations before (Xf') and after (Xa') ETKF update
20 !5) ens_ob, array of HXf
21 !6) oberrvar, observation error variance, listed in the same sequence as HXf
22 !7) obs, observations assmiilated
23 !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
24 !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
25 !10) nstartaccum1, the cycle from which the accumulation of previous naccumt1 cycles in 8) starts
26 !11) nstartaccum2, the cycle from which the accumulation of previous naccumt2 in 9) starts
27 !12) tainflatinput, pre-specified inflation, if not using adaptive inflation
28 !13) rhoinput, pre-specified rho factor, if not using adaptively determined rho factor
29 !14) nout, record number for output squared innovations and the forecast error variance
30 !-----------------------------------------------------------------------projected onto ensemble subspace, which is related to 8) and 9)
31
32 implicit none
33
34 integer, intent(in) :: nanals,ndim,nobs
35 real, intent(inout), dimension(ndim,nanals) :: ens
36 real, intent(inout), dimension(nobs,nanals) :: ens_ob
37 real, intent(in), dimension(nobs) :: oberrvar
38 real, intent(in), dimension(nobs) :: obs
39 real, intent(in) :: tainflatinput,rhoinput
40 integer, intent(in) :: nout,naccumt1,naccumt2,nstartaccum1,nstartaccum2
41
42 integer :: n ! Loop counters.
43 integer :: nmin ! Minimum nout value to use.
44 real :: nanals_inv ! 1 / nanals.
45 real :: ainflat_mean ! Rolling mean inflation factor.
46
47 real, dimension(nobs) :: ensmean_ob
48 real, dimension(nobs) :: obsinc
49 integer :: ij, nanal, i, j ! Loop counters.
50 integer :: info, lwork
51 real, allocatable, dimension(:) :: work
52 real, dimension(1) :: workt
53 real, dimension(nanals) :: eignv, eignv1
54 real, dimension(nanals, nanals) :: hzthz, C, TR
55 real, dimension(nanals-1, nanals) :: CT
56 real, dimension(nanals, nanals-1) :: T
57 ! real, dimension(nanals, nanals-1) :: cgamma
58 ! real, dimension(nobs,nanals-1) :: E
59 real, dimension(nobs,nanals) :: enspert_ob_tmp
60 ! file for output eigen values filename_eignv
61 ! file for output factors filename_factor
62 character (len=150) :: filename_factor
63 ! character (len=150) :: filename_eignv
64 ! file for squared innovation vectors
65 ! file for squared projection of fcst err variance onto ensemble subspace
66 ! character (len=150) :: filename_inno2
67 ! character (len=150) :: filename_proj2
68 ! real, dimension(nanals-1) :: proj
69 real :: tracehpfht
70 ! real :: ta, tb
71 ! real :: sum_eignv
72 ! real :: aftersum_eignv
73 real :: ainflat
74 ! real :: tainflat, rho
75 ! real :: proj2, proj2sum
76 ! real :: proj2sum1
77 ! real :: squareinnosum, squareinnosum1
78 real :: squareinno
79
80 if (trace_use) call da_trace_entry("da_solve_etkf")
81
82 filename_factor = "inflation_factor.dat"
83 open(109,form="unformatted",file=filename_factor,access="direct",recl=8)
84
85 !------------------------------------------------------------------------------
86 ! [1] Compute mean(H(xf)) and H(xf)'.
87 !------------------------------------------------------------------------------
88
89 nanals_inv = 1.0 / real(nanals)
90 do ij = 1, nobs
91 ensmean_ob(ij) = sum(ens_ob(ij,:)) * nanals_inv
92 end do
93
94 do nanal = 1, nanals
95 ens_ob(:,nanal) = ens_ob(:,nanal) - ensmean_ob(:)
96 end do
97
98 !------------------------------------------------------------------------------
99 ! [2] Calculate HZTHZ in Bishop et al. 2001
100 !------------------------------------------------------------------------------
101
102 do i = 1, nobs
103 enspert_ob_tmp(i,:) = ens_ob(i,:)/sqrt(oberrvar(i))
104 end do
105
106 call da_innerprod(enspert_ob_tmp,hzthz,nobs,nanals)
107 hzthz = hzthz/float(nanals-1)
108
109 !------------------------------------------------------------------------------
110 ! [3] Calculate C and Gamma in Bishop et al. 2001
111 !------------------------------------------------------------------------------
112 ! in outputs, hzthz contains C, eignv contains gamma
113 call dsyev('V', 'L', nanals , hzthz, nanals, eignv1, workt, -1, info)
114 lwork = int(workt(1))
115 allocate (work(lwork))
116 call dsyev('V', 'L', nanals , hzthz, nanals, eignv1, work, lwork, info)
117 deallocate(work)
118
119 ! note eignv1 output from SSYEV is in ascending order !!!
120 ! re-order both eigenvalues and eigenvectors first
121
122 do i = 1, nanals
123 eignv(i) = eignv1(nanals-i+1)
124 end do
125
126 ! note eigenvectors output from SSYEV are stored in columns !!!
127 do i = 1, nanals
128 C(:,i) = hzthz(:,nanals-i+1)
129 end do
130
131 !------------------------------------------------------------------------------
132 ! [4] Calculate inflation factor:
133 !------------------------------------------------------------------------------
134
135 if ( naccumt1 > 0 ) then
136
137 tracehpfht = 0.0
138 do i = 1, nanals-1
139 tracehpfht = tracehpfht + eignv(i)
140 end do
141
142 obsinc(:) = ensmean_ob(:) - obs(:)
143 squareinno = sum( obsinc(:) * obsinc(:) / oberrvar(:) )
144 ainflat = ( squareinno - real(nobs) ) / tracehpfht
145 write(109,rec=nout) ainflat
146
147 ! Calculate running mean inflation factor:
148 nmin = max( 1, nout - naccumt1 + 1 )
149 ainflat_mean = 0.0
150 do n = nmin, nout
151 read(109,rec=n) ainflat
152 ainflat_mean = ainflat_mean + ainflat
153 end do
154 ainflat_mean = ainflat_mean / real( nout - nmin + 1 )
155 write (unit=stdout,fmt='(/a,f15.5)') " Current Inflation factor = ", ainflat
156 write (unit=stdout,fmt='(a,f15.5)') " Rolling mean inflation factor = ", ainflat_mean
157 else
158 ! This is for pre-specified inflation factor
159 ainflat = tainflatinput
160 end if
161
162 !------------------------------------------------------------------------------
163 ! [5] Calculate the grand transformation matrix:
164 !------------------------------------------------------------------------------
165
166 do i = 1, nanals
167 do j = 1, nanals-1
168 T(i,j) = C(i,j)*sqrt(1.0/(eignv(j)+1.0))
169 end do
170 end do
171
172 do i = 1, nanals-1
173 do j = 1, nanals
174 CT(i,j) = C(j,i)
175 end do
176 end do
177
178 call da_matmulti(T,CT,TR,nanals,nanals,nanals-1)
179
180 ! Apply inflation:
181 tr = sqrt(ainflat_mean) * tr
182
183 !------------------------------------------------------------------------------
184 ! [4] Calculate the rescaled ETKF perturbations
185 !------------------------------------------------------------------------------
186
187 call da_matmultiover(ens, tr, ndim, nanals)
188
189 if (trace_use) call da_trace_exit("da_solve_etkf")
190
191 end subroutine da_solve_etkf
192