main.F90
11.3 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
program enkf_eco
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Program which uses the ensemble Kalman Filter to assimilate !
! in the 1D version of the ecosystem model. !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use mod_dimensions
use mod_states
use mod_observations
use mod_variances
use m_cholesky
use m_physpar
use m_mkmxlayer
use m_reference
use m_dump
use m_load
use m_tecfld
use m_pseudo
use m_load_ens
use m_dump_ens
use m_measure
use m_integrate
use m_ens_stat
use m_random
use m_analysis
use m_set_random_seed2
implicit none
type(states), allocatable :: mem(:) ! Contains the ensemble of model states
type(states) ref(ndim) ! Holds the full reference solution
type(states) bguess ! Holds the initial values (best guess solution)
! used when generating the init. ensemble
type(states) estimate(ndim)! The solution estimate
type(states) variance(ndim)! The variance estimate
type(states) forw ! Temporary model state
type(local_states) bnd0,bnd0w ! Boundary conditions
type(local_states) bndH,bndHw ! Boundary conditions
real workA(kdim,antvar)
real workB(antvar)
real cov(kdim,kdim), chov(kdim,kdim) ! Cov matrix and its Sqrt
type(states) ave ! Temporary model state
type(states) var ! Temporary model state
type(states) Mean ! Ensemble mean
real, allocatable :: cov1(:,:), chov1(:,:)
type(variances) vars
real, allocatable :: S(:,:) ! S(3*kdim,nrsamp)
real MM(ndim) ! Mixed layer depth
type(observations) d(nrobs,nrmes_t) , d_temp
type(lstates) meas_action
type(local_states) mse
real, allocatable :: U0(:,:),sig0(:),VT0(:,:)
real, allocatable :: U(:,:), sig(:),VT1(:,:),A0(:,:)
real, allocatable, dimension(:) :: work
real scale ! 1/nrsamp
! ME
type(local_states) mseres
! ME end
type(local_states) totrms
type(local_states) totsd
real,parameter :: range = 50. ! 50 m decorrelation length (in z)
real tfin ! final integration time
logical iniref ! genereate new ref sol or not
logical iniens
real dt,dz
real totdepth
integer mode, i_xp
integer, parameter :: n_xp = 10 ! repeat 10 times experiment
integer i,j,k,m,n,iens,na,nb,nstep,nrsamp
logical ex
character(len=1) tag1
character(len=3) tag3
real t0,t1
integer ierr,lwork
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Reading the parameters from the file "infile"
call physpar(tfin,dt,totdepth,dz,iniref,iniens,vars,mode,meas_action,nrsamp)
print *,'physpar done!'
allocate (mem(nrsamp))
allocate (cov1(nrsamp,nrsamp))
allocate (chov1(nrsamp,nrsamp))
allocate (S(3*kdim,nrsamp))
call mkmxlayer(MM,dt)
print *,'mxlayer done!'
call set_random_seed2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Sqrt spatial covariance matrix
print *, 'initial cov ...'
do i=1,kdim
do j=1,i
cov(i,j)=gaussian(float(abs(i-j))*dz/range)
cov(j,i)=cov(i,j)
end do
end do
print *, 'beginning cholesky(cov) ...'
chov=cov
call cholesky(chov,kdim)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! We need this matrix to be of dimension (nrsamp,nrsamp) to be able to use it flater when computing othogonal matrix.
! Sqrt spatial covariance matrix
print *, 'initial cov ...'
do i=1,nrsamp
do j=1,i
cov1(i,j)=gaussian(float(abs(i-j))*dz/range)
cov1(j,i)=cov1(i,j)
end do
end do
print *, 'beginning cholesky(cov1) ...'
chov1=cov1
call cholesky(chov1,nrsamp)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i_xp=1, n_xp
! Set initial and boundary conditions
! Boundary cond (no flux at surface)
bnd0=0.0
! Boundary cond (specified value at bottom)
bndH%N=10.0
bndH%P=0.0
bndH%H=0.0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Reference solution
print *,'iniref ',iniref
inquire(file='refsol.uf',exist=ex)
if ((iniref).or.(.not.ex)) then
print *,'Computing reference solution'
call reference(ref,bguess,MM,dz,totdepth,dt,tfin,vars,bnd0,bndH,chov)
call dump('refsol.uf',ref)
else
print *,'Loading reference solution from file'
call load('refsol.uf',ref)
!print *, ' Init Ref N(kdim) ' , ref(1)%N(kdim)
endif
print *,'done!'
totsd=0. ; totrms=0.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initial ensemble
print *,'iniens= ',iniens
inquire(file='iniens.uf',exist=ex)
if ((iniens).or.(.not.ex)) then
print *,'Computing initial ensemble'
do iens=1,nrsamp
call pseudo(worka,kdim,antvar,chov)
mem(iens)%N=max(0.0, bguess%N+sqrt(vars%ini%N)*workA(:,1) )
mem(iens)%P=max(0.0, bguess%P+sqrt(vars%ini%P)*workA(:,2) )
mem(iens)%H=max(0.0, bguess%H+sqrt(vars%ini%H)*workA(:,3) )
enddo
call dump_ens('iniens.uf',mem,nrsamp)
else
print *,'Loading initial ensemble from file'
call load_ens('iniens.uf',mem,nrsamp)
endif
print *,'done!'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!NEW PART WITH GENERATION OF START ENSEMBLE OF "NRSAMP_INI" MEMBERS AND THEN RESAMPLE "NRSAMP" MEMBERS .
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Compute ensemble mean
! scale = 1.0/float(nrsamp_ini)
! Mean=0.0
! var=0.0
! do i=1,nrsamp_ini
! Mean = Mean + mem(i)
! var=var+mem(i)*mem(i)
! enddo
! Mean=scale*Mean
! var=var*scale
! Substract mean from ensemble members
! do i=1,nrsamp_ini
! mem(i)%N=mem(i)%N-Mean%N
! mem(i)%P=mem(i)%P-Mean%P
! mem(i)%H=mem(i)%H-Mean%H
! var%N = var%N + mem(i)%N**2
! var%P = var%P + mem(i)%P**2
! var%H = var%H + mem(i)%H**2
! mem(i)%N=var%N*mem(i)%N
! mem(i)%P=var%P*mem(i)%P
!mem(i)%H=var%H*mem(i)%H
!enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Make an othogonal VT1 used as linear combination for final ensemble
! allocate (A0(nrsamp,nrsamp), U(nrsamp,nrsamp) , sig(nrsamp) , VT1(nrsamp,nrsamp))
! lwork=2*max(3*nrsamp+3*kdim,5*nrsamp)
! allocate(work(lwork))
!call pseudo(A0,nrsamp,nrsamp,chov1)
!call dgesvd('N', 'S', nrsamp, nrsamp, A0, nrsamp, sig, U, nrsamp, VT1, nrsamp, work, lwork, ierr)
! if(ierr/=0) print *, 'ierr',ierr
! deallocate(A0, sig, U, work)
! Compute SVD of oversized ensemble
! allocate (U0(3*kdim,nrsamp_ini))
!allocate (sig0(nrsamp_ini))
!allocate(VT0(1,1))
!lwork=2*max(3*nrsamp_ini+3*kdim,5*nrsamp_ini)
!allocate(work(lwork))
!sig0=0.0
!call dgesvd('S', 'N', 3*kdim, nrsamp_ini, mem, 3*kdim,sig0, U0, 3*kdim, VT0, nrsamp_ini , work, lwork, ierr)
!deallocate(work,VT0)
! if (ierr /= 0) then
! print *,'analysis: ierr from call dgesvd 0= ',ierr; stop
!endif
! do i=1,nrsamp_ini
! mem(i)%N=mem(i)%N/var%N
! mem(i)%P=mem(i)%P/var%P
! mem(i)%H=mem(i)%H/var%H
! enddo
! Generate first nrsamp members and add mean
! do i=1,nrsamp
! do j=1,nrsamp
! mem(i)%N=(U0(1:kdim,i)*sig0(i)/sqrt(float(nrsamp))*VT1(j,i))+Mean%N
! mem(i)%P=(U0(kdim+1:2*kdim,i)*sig0(i)/sqrt(float(nrsamp))*VT1(j,i))+Mean%P
! mem(i)%H=(U0(2*kdim+1:3*kdim,i)*sig0(i)/sqrt(float(nrsamp))*VT1(j,i))+Mean%H
!enddo
!enddo
!deallocate(U0,sig0,VT1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!END OF NEW PART
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Generate measurements by measuring the reference solution
call measure(d,ref,vars,meas_action)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Time stepping
! -compute statistics of initial ensemble
call ens_stat(mem,ave,var,nrsamp)
estimate(1)=ave
variance(1)=var
do m=1,nrmes_t+1
na=(m-1)*nres+1
nb=m*nres+1
print *,'na,nb,nres,nrmes_t ',na,nb,nres,nrmes_t
do nstep=na+1,nb
t0=dt*float(nstep-2)
t1=dt*float(nstep-1)
print '(a,2I5,2f8.2)','m,nstep,t0,t1:',m,nstep,t0,t1
do iens=1,nrsamp
bnd0w=bnd0
call random(workb,3)
bndHw%N=bndH%N+sqrt(vars%bH%N)*workb(1)
bndHw%P=bndH%P+sqrt(vars%bH%P)*workb(2)
bndHw%H=bndH%H+sqrt(vars%bH%H)*workb(3)
forw=mem(iens)
call integrate(forw,MM,bnd0w,bndHw,totdepth,t0,t1,dt)
call pseudo(worka,kdim,antvar,chov)
mem(iens)%N=max(0.0, forw%N+sqrt(dt)*sqrt(vars%dyn%N)*worka(:,1))
mem(iens)%P=max(0.0, forw%P+sqrt(dt)*sqrt(vars%dyn%P)*worka(:,2))
mem(iens)%H=max(0.0, forw%H+sqrt(dt)*sqrt(vars%dyn%H)*worka(:,3))
enddo
call ens_stat(mem,ave,var,nrsamp)
estimate(nstep)=ave
variance(nstep)=var
enddo
if (nstep < ndim) then
if (mode == 1) then
print *,'calling analysis',m
print *,'Meas. at t1 ',t1
call prepare_and_analysis(mem,d(:,m),nrsamp)
endif
endif
enddo
write(tag1,'(i1.1)')mode
write(tag3,'(i3.3)')nrsamp
call tecfld('output_'//tag1,ref,estimate,variance,dz,dt)
open(10,file='mse_'//tag3//'.dat')
do i=1,ndim
mse%N=sum(variance(i)%N(:))/float(kdim)
mse%P=sum(variance(i)%P(:))/float(kdim)
mse%H=sum(variance(i)%H(:))/float(kdim)
write(10,'(f12.4,4f12.4)')float(i-1)*dt,sqrt(mse%N), sqrt(mse%P), &
sqrt(mse%H), sqrt(mse%N+mse%P+mse%H)
totsd%N = totsd%N + mse%N
totsd%P = totsd%P + mse%P
totsd%H = totsd%H + mse%H
enddo
close(10)
! ME for the residuals:
open(10,file='mseres_'//tag3//'.dat')
do i=1,ndim
mseres%N=sum((estimate(i)%N(:)-ref(i)%N(:))**2)/float(kdim)
mseres%P=sum((estimate(i)%P(:)-ref(i)%P(:))**2)/float(kdim)
mseres%H=sum((estimate(i)%H(:)-ref(i)%H(:))**2)/float(kdim)
write(10,'(f12.4,4f12.4)')float(i-1)*dt,sqrt(mseres%N), &
sqrt(mseres%P), sqrt(mseres%H), &
sqrt(mseres%N+mseres%P+mseres%H)
totrms%N = totrms%N + mseres%N
totrms%P = totrms%P + mseres%P
totrms%H = totrms%H + mseres%H
enddo
close(10)
! ME res. end.
enddo
totsd%N = sqrt(totsd%N / float(ndim*n_xp) )
totsd%P = sqrt(totsd%P / float(ndim*n_xp) )
totsd%H = sqrt(totsd%H / float(ndim*n_xp) )
totrms%N = sqrt(totrms%N / float(ndim*n_xp) )
totrms%P = sqrt(totrms%P / float(ndim*n_xp) )
totrms%H = sqrt(totrms%H / float(ndim*n_xp) )
open(11, file='totsd.dat',position='append')
write(11,'(i5.3,4g15.5)') nrsamp, totsd%N, totsd%P, totsd%H, totsd%N+totsd%P+totsd%H
close(11)
open(12, file='totrms.dat',position='append')
write(12,'(i5.3,4g15.5)') nrsamp, totrms%N, totrms%P, totrms%H, totrms%N+totrms%P+totrms%H
close(12)
deallocate (mem)
deallocate (cov1)
deallocate (chov1)
deallocate (S)
end program enkf_eco