gen_be_stage2.f90
References to this file elsewhere.
1 program gen_be_stage2
2
3 use da_control, only : stdout, stderr, filename_len
4 use da_tools_serial, only : da_get_unit, da_advance_cymdh
5 use da_gen_be, only : da_eof_decomposition,da_eof_decomposition_test
6
7 implicit none
8
9 real, parameter :: variance_threshold = 1e-6 ! Percentage of <psi psi> variance discarded.
10
11 character*10 :: start_date, end_date ! Starting and ending dates.
12 character*10 :: date, new_date ! Current date (ccyymmddhh).
13 character*10 :: variable ! Variable name
14 character(len=filename_len) :: filename ! Input filename.
15 character*3 :: ce ! Member index -> character.
16 integer :: ni, nj, nk, nkdum ! Grid dimensions.
17 integer :: i, j, k, member, k2, k3, m ! Loop counters.
18 integer :: b ! Bin marker.
19 integer :: sdate, cdate, edate ! Starting, current ending dates.
20 integer :: interval ! Period between dates (hours).
21 integer :: ne ! Number of ensemble members.
22 integer :: mmax ! Maximum mode (after variance truncation)..
23 integer :: bin_type ! Type of bin to average over.
24 integer :: num_bins ! Number of bins (3D fields).
25 integer :: num_bins2d ! Number of bins (2D fields).
26 real :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
27 real :: binwidth_lat ! Used if bin_type = 2 (degrees).
28 real :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
29 real :: binwidth_hgt ! Used if bin_type = 2 (m).
30 real :: coeffa, coeffb ! Accumulating mean coefficients.
31 real :: total_variance ! Total variance of <psi psi> matrix.
32 real :: cumul_variance ! Cumulative variance of <psi psi> matrix.
33 real :: summ ! Summation dummy.
34 logical :: first_time ! True if first file.
35 logical :: testing_eofs ! True if testing EOF decomposition.
36
37 real, allocatable :: latitude(:,:) ! Latitude (degrees, from south).
38 real, allocatable :: height(:,:,:) ! Height field.
39 real, allocatable :: psi(:,:,:) ! psi.
40 real, allocatable :: chi(:,:,:) ! chi.
41 real, allocatable :: temp(:,:,:) ! Temperature.
42 real, allocatable :: ps(:,:) ! Surface pressure.
43 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
44 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
45 integer, allocatable:: bin_pts(:) ! Number of points in bin (3D fields).
46 integer, allocatable:: bin_pts2d(:) ! Number of points in bin (2D fields).
47 real, allocatable :: covar1(:) ! Covariance between input fields.
48 real, allocatable :: covar2(:,:) ! Covariance between input fields.
49 real, allocatable :: covar3(:,:,:) ! Covariance between input fields.
50 real, allocatable :: var1(:) ! Autocovariance of field.
51 real, allocatable :: var2(:,:,:) ! Autocovariance of field.
52 real, allocatable :: var2_inv(:,:,:) ! Inverse Autocovariance of field.
53
54 real, allocatable :: work(:,:) ! EOF work array.
55 real, allocatable :: evec(:,:) ! Gridpoint eigenvectors.
56 real, allocatable :: eval(:) ! Gridpoint sqrt(eigenvalues).
57 real, allocatable :: LamInvET(:,:) ! ET/sqrt(Eigenvalue).
58 real, allocatable :: regcoeff1(:) ! psi/chi regression cooefficient.
59 real, allocatable :: regcoeff2(:,:) ! psi/ps regression cooefficient.
60 real, allocatable :: regcoeff3(:,:,:) ! psi/T regression cooefficient.
61
62 namelist / gen_be_stage2_nl / start_date, end_date, interval, &
63 ne, testing_eofs
64
65 integer :: ounit,iunit,namelist_unit
66
67
68 stderr = 0
69 stdout = 6
70
71 !---------------------------------------------------------------------------------------------
72 write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
73 !---------------------------------------------------------------------------------------------
74
75 call da_get_unit(ounit)
76 call da_get_unit(iunit)
77 call da_get_unit(namelist_unit)
78
79
80 start_date = '2004030312'
81 end_date = '2004033112'
82 interval = 24
83 ne = 1
84 testing_eofs = .true.
85
86 open(unit=namelist_unit, file='gen_be_stage2_nl.nl', &
87 form='formatted', status='old', action='read')
88 read(namelist_unit, gen_be_stage2_nl)
89 close(namelist_unit)
90
91 read(start_date(1:10), fmt='(i10)')sdate
92 read(end_date(1:10), fmt='(i10)')edate
93 write(6,'(4a)')' Computing regression coefficients'
94 write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
95 write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
96 write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
97
98 date = start_date
99 cdate = sdate
100
101 !---------------------------------------------------------------------------------------------
102 write(6,'(2a)')' [2] Read fields, and calculate correlations'
103 !---------------------------------------------------------------------------------------------
104
105 first_time = .true.
106
107 do while ( cdate <= edate )
108 write(6,'(a,a)')' Processing data for date ', date
109
110 do member = 1, ne
111
112 write(ce,'(i3.3)')member
113
114 ! Read Full-fields:
115 variable = 'fullflds'
116 filename = trim(variable)//'/'//date(1:10)
117 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
118 open (iunit, file = filename, form='unformatted')
119
120 read(iunit)ni, nj, nk
121 if ( first_time ) then
122 write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk
123 allocate( latitude(1:ni,1:nj) )
124 allocate( height(1:ni,1:nj,1:nk) )
125 allocate( bin(1:ni,1:nj,1:nk) )
126 allocate( bin2d(1:ni,1:nj) )
127 allocate( psi(1:ni,1:nj,1:nk) )
128 allocate( chi(1:ni,1:nj,1:nk) )
129 allocate( temp(1:ni,1:nj,1:nk) )
130 allocate( ps(1:ni,1:nj) )
131 end if
132 read(iunit)latitude
133 read(iunit)height
134
135 close(iunit)
136
137 if ( first_time ) then
138
139 ! Read bin info:
140 filename = 'bin.data'
141 open (iunit, file = filename, form='unformatted')
142 read(iunit)bin_type
143 read(iunit)lat_min, lat_max, binwidth_lat
144 read(iunit)hgt_min, hgt_max, binwidth_hgt
145 read(iunit)num_bins, num_bins2d
146 read(iunit)bin(1:ni,1:nj,1:nk)
147 read(iunit)bin2d(1:ni,1:nj)
148 close(iunit)
149
150 allocate( bin_pts(1:num_bins) )
151 allocate( bin_pts2d(1:num_bins2d) )
152 allocate( covar1(1:num_bins) )
153 allocate( covar2(1:nk,1:num_bins2d) )
154 allocate( covar3(1:nk,1:nk,1:num_bins2d) )
155 allocate( var1(1:num_bins) )
156 allocate( var2(1:nk,1:nk,1:num_bins2d) )
157 bin_pts(:) = 0
158 bin_pts2d(:) = 0
159 covar1(:) = 0.0
160 covar2(:,:) = 0.0
161 covar3(:,:,:) = 0.0
162 var1(:) = 0.0
163 var2(:,:,:) = 0.0
164 first_time = .false.
165 end if
166
167 ! Read psi:
168 variable = 'psi'
169 filename = trim(variable)//'/'//date(1:10)
170 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
171 open (iunit, file = filename, form='unformatted')
172 read(iunit)ni, nj, nk
173 read(iunit)psi
174 close(iunit)
175
176 ! Read chi:
177 variable = 'chi'
178 filename = trim(variable)//'/'//date(1:10)
179 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
180 open (iunit, file = filename, form='unformatted')
181 read(iunit)ni, nj, nk
182 read(iunit)chi
183 close(iunit)
184
185 ! Read T:
186 variable = 't'
187 filename = trim(variable)//'/'//date(1:10)
188 filename = trim(filename)//'.'//trim(variable)//'.e'//ce
189 open (iunit, file = filename, form='unformatted')
190 read(iunit)ni, nj, nk
191 read(iunit)temp
192 close(iunit)
193
194 ! Read ps:
195 variable = 'ps'
196 filename = trim(variable)//'/'//date(1:10)
197 filename = trim(filename)//'.'//trim(variable)//'.e'//ce//'.01'
198 open (iunit, file = filename, form='unformatted')
199 read(iunit)ni, nj, nkdum
200 read(iunit)ps
201 close(iunit)
202
203 ! Calculate psi/chi covariances:
204
205 do k = 1, nk
206 do j = 1, nj
207 do i = 1, ni
208 b = bin(i,j,k)
209 bin_pts(b) = bin_pts(b) + 1
210 coeffa = 1.0 / real(bin_pts(b))
211 coeffb = real(bin_pts(b)-1) * coeffa
212 covar1(b) = coeffb * covar1(b) + coeffa * psi(i,j,k) * chi(i,j,k)
213 var1(b) = coeffb * var1(b) + coeffa * psi(i,j,k) * psi(i,j,k)
214 end do
215 end do
216 end do
217
218 ! Calculate psi/ps, and psi/T, and psi/psi covariances:
219
220 do j = 1, nj
221 do i = 1, ni
222 b = bin2d(i,j)
223 bin_pts2d(b) = bin_pts2d(b) + 1
224 coeffa = 1.0 / real(bin_pts2d(b))
225 coeffb = real(bin_pts2d(b)-1) * coeffa
226 do k = 1, nk
227 ! psi/ps:
228 covar2(k,b) = coeffb * covar2(k,b) + coeffa * ps(i,j) * psi(i,j,k)
229
230 ! psi/T:
231 do k2 = 1, nk
232 covar3(k,k2,b) = coeffb * covar3(k,k2,b) + &
233 coeffa * temp(i,j,k) * psi(i,j,k2)
234 end do
235
236 ! psi/psi (symmetric):
237 do k2 = 1, k
238 var2(k,k2,b) = coeffb * var2(k,k2,b) + &
239 coeffa * psi(i,j,k) * psi(i,j,k2)
240 end do
241 end do
242 end do
243 end do
244
245 end do ! End loop over ensemble members.
246
247 ! Calculate next date:
248 call da_advance_cymdh( date, interval, new_date )
249 date = new_date
250 read(date(1:10), fmt='(i10)')cdate
251 end do ! End loop over times.
252
253 deallocate ( latitude ) ! Not needed any more.
254 deallocate ( height ) ! Not needed any more.
255
256 ! Fill in psi/psi covariance by symmetry:
257 do b = 1, num_bins2d
258 do k = 1, nk
259 do k2 = k+1, nk ! Symmetry.
260 var2(k,k2,b) = var2(k2,k,b)
261 end do
262 end do
263 end do
264
265 !---------------------------------------------------------------------------------------------
266 write(6,'(2a)')' [3] Calculate eigenvectors, eigenvalues and inverse for psi/psi covariance '
267 !---------------------------------------------------------------------------------------------
268
269 allocate( work(1:nk,1:nk) )
270 allocate( evec(1:nk,1:nk) )
271 allocate( eval(1:nk) )
272 allocate( LamInvET(1:nk,1:nk) )
273 allocate( var2_inv(1:nk,1:nk,1:num_bins2d) )
274
275 do b = 1, num_bins2d
276 LamInvET(:,:) = 0.0
277 work(1:nk,1:nk) = var2(1:nk,1:nk,b)
278 call da_eof_decomposition( nk, work, evec, eval )
279
280 if ( testing_eofs ) then
281 call da_eof_decomposition_test( nk, work, evec, eval )
282 testing_eofs = .false.
283 end if
284
285 ! Truncate eigenvalues to ensure inverse is not dominated by rounding error:
286 summ = 0.0
287 do m = 1, nk
288 summ = summ + eval(m)
289 end do
290 total_variance = summ
291
292 cumul_variance = 0.0
293 mmax = nk
294 do m = 1, nk
295 cumul_variance = cumul_variance + eval(m) / total_variance
296 if ( cumul_variance > 1.0 - variance_threshold ) then
297 mmax = m - 1
298 exit
299 end if
300 end do
301 write(6,'(2(a,i6),2(a,1pe11.5))') ' Bin = ', b, ', <psipsi> truncation = ', mmax, &
302 ', Total Variance = ', total_variance, &
303 ', Condition number = ', eval(1) / eval(nk-1)
304
305 ! Lam{-1} . E^T:
306 do k = 1, nk
307 do m = 1, mmax
308 LamInvET(m,k) = evec(k,m) / eval(m)
309 end do
310 end do
311
312 ! <psi psi>^{-1} = E . Lam{-1} . E^T:
313
314 do k = 1, nk
315 do k2 = 1, k
316 summ = 0.0
317 do m = 1, nk
318 summ = summ + evec(k,m) * LamInvET(m,k2)
319 end do
320 var2_inv(k,k2,b) = summ
321 end do
322 end do
323
324 do k = 1, nk
325 do k2 = k+1, nk ! Symmetry.
326 var2_inv(k,k2,b) = var2_inv(k2,k,b)
327 end do
328 end do
329 end do
330 !---------------------------------------------------------------------------------------------
331 write(6,'(2a)')' [4] Calculate regression coefficients '
332 !---------------------------------------------------------------------------------------------
333
334 allocate( regcoeff1(1:num_bins) )
335 allocate( regcoeff2(1:nk,1:num_bins2d) )
336 allocate( regcoeff3(1:nk,1:nk,1:num_bins2d) )
337
338 ! psi/chi:
339 do b = 1, num_bins
340 regcoeff1(b) = covar1(b) / var1(b)
341 end do
342
343 ! psi/ps:
344 do b = 1, num_bins2d
345 do k = 1, nk
346 summ = 0.0
347 do k2 = 1, nk
348 summ = summ + covar2(k2,b) * var2_inv(k2,k,b)
349 end do
350 regcoeff2(k,b) = summ
351 end do
352 end do
353
354 ! psi/T:
355 do b = 1, num_bins2d
356 do k = 1, nk
357 do k2 = 1, nk
358 summ = 0.0
359 do k3 = 1, nk
360 summ = summ + covar3(k,k3,b) * var2_inv(k3,k2,b)
361 end do
362 regcoeff3(k,k2,b) = summ
363 end do
364 end do
365 end do
366
367 ! Output regression coefficients for use in 3/4D-Var:
368 filename = 'gen_be_stage2.dat'
369 open (ounit, file = filename, form='unformatted')
370 write(ounit)ni, nj, nk
371 write(ounit)num_bins, num_bins2d
372 write(ounit)regcoeff1
373 write(ounit)regcoeff2
374 write(ounit)regcoeff3
375 close(ounit)
376
377 do k = 1, nk
378 do j = 1, nj
379 b = bin2d(1,j)
380 write(61,'(3i6,1pe13.5)')k, j, b, covar2(k,b)
381 write(62,'(3i6,1pe13.5)')k, j, b, var2(k,k,b)
382 write(63,'(3i6,1pe13.5)')k, j, b, var2_inv(k,k,b)
383 end do
384 end do
385
386 end program gen_be_stage2