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