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, nt, i, j, k ! 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) :: cgamma, T
55 real, dimension(nobs,nanals-1) :: E
56 real, dimension(nobs,nanals) :: enspert_ob_tmp
57 ! file for output eigen values filename_eignv
58 ! file for output factors filename_factor
59 character (len=150) :: filename_factor,filename_eignv
60 ! file for squared innovation vectors
61 ! file for squared projection of fcst err variance onto ensemble subspace
62 character (len=150) :: filename_inno2,filename_proj2
63 real, dimension(nanals-1) :: proj
64 real :: ta, tb, tracehpfht, aftersum_eignv, sum_eignv
65 real :: ainflat, tainflat, rho
66 real :: proj2, proj2sum, proj2sum1
67 real :: squareinnosum, squareinno, squareinnosum1
68
69 filename_factor = trim(directory)//"/inflation_factor.dat"
70 open(109,form="unformatted",file=filename_factor,access="direct",recl=8)
71
72 !------------------------------------------------------------------------------
73 ! [1] Compute mean(H(xf)) and H(xf)'.
74 !------------------------------------------------------------------------------
75
76 nanals_inv = 1.0 / real(nanals)
77 do ij = 1, nobs
78 ensmean_ob(ij) = sum(ens_ob(ij,:)) * nanals_inv
79 end do
80
81 do nanal = 1, nanals
82 ens_ob(:,nanal) = ens_ob(:,nanal) - ensmean_ob(:)
83 end do
84
85 !------------------------------------------------------------------------------
86 ! [2] Calculate HZTHZ in Bishop et al. 2001
87 !------------------------------------------------------------------------------
88
89 do i = 1, nobs
90 enspert_ob_tmp(i,:) = ens_ob(i,:)/sqrt(oberrvar(i))
91 end do
92
93 call da_innerprod(enspert_ob_tmp,hzthz,nobs,nanals)
94 hzthz = hzthz/float(nanals-1)
95
96 !------------------------------------------------------------------------------
97 ! [3] Calculate C and Gamma in Bishop et al. 2001
98 !------------------------------------------------------------------------------
99 ! in outputs, hzthz contains C, eignv contains gamma
100 call dsyev('V', 'L', nanals , hzthz, nanals, eignv1, workt, -1, info)
101 lwork = int(workt(1))
102 allocate (work(lwork))
103 call dsyev('V', 'L', nanals , hzthz, nanals, eignv1, work, lwork, info)
104 deallocate(work)
105
106 ! note eignv1 output from SSYEV is in ascending order !!!
107 ! re-order both eigenvalues and eigenvectors first
108
109 do i = 1, nanals
110 eignv(i) = eignv1(nanals-i+1)
111 end do
112
113 ! note eigenvectors output from SSYEV are stored in columns !!!
114 do i = 1, nanals
115 C(:,i) = hzthz(:,nanals-i+1)
116 end do
117
118 !------------------------------------------------------------------------------
119 ! [4] Calculate inflation factor:
120 !------------------------------------------------------------------------------
121
122 if ( naccumt1 > 0 ) then
123
124 tracehpfht = 0.0
125 do i = 1, nanals-1
126 tracehpfht = tracehpfht + eignv(i)
127 end do
128
129 obsinc(:) = ensmean_ob(:) - obs(:)
130 squareinno = sum( obsinc(:) * obsinc(:) / oberrvar(:) )
131 ainflat = ( squareinno - real(nobs) ) / tracehpfht
132 write(109,rec=nout) ainflat
133
134 ! Calculate running mean inflation factor:
135 nmin = max( 1, nout - naccumt1 + 1 )
136 ainflat_mean = 0.0
137 do n = nmin, nout
138 read(109,rec=n) ainflat
139 ainflat_mean = ainflat_mean + ainflat
140 end do
141 ainflat_mean = ainflat_mean / real( nout - nmin + 1 )
142 write (unit=stdout,fmt='(a,i6,2f15.5)') "nout, ainflat, ainflat_mean = ", &
143 nout, ainflat, ainflat_mean
144 else
145 ! This is for pre-specified inflation factor
146 ainflat = tainflatinput
147 end if
148
149 !------------------------------------------------------------------------------
150 ! [5] Calculate the grand transformation matrix:
151 !------------------------------------------------------------------------------
152
153 do i = 1, nanals
154 do j = 1, nanals-1
155 T(i,j) = C(i,j)*sqrt(1./(eignv(j)+1.))
156 end do
157 end do
158
159 do i = 1, nanals-1
160 do j = 1, nanals
161 CT(i,j) = C(j,i)
162 end do
163 end do
164
165 call da_matmulti(T,CT,TR,nanals,nanals,nanals-1)
166
167 ! Apply inflation:
168 tr = sqrt(ainflat_mean) * tr
169
170 !------------------------------------------------------------------------------
171 ! [4] Calculate the ETKF initial perturbations
172 !------------------------------------------------------------------------------
173
174 call da_matmultiover(ens, tr, ndim, nanals)
175
176 end subroutine da_solve_etkf
177