Blame view

m_analysis.F90 10.7 KB
26362238   Dany Dumont   premier depot
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

module m_analysis
contains
subroutine prepare_and_analysis(mem,obs,nrsamp)
!  This subroutine computes the EnKF analysis. The improved analysed 
!  estimate psia(iens), when psif(iens) is the forcast, is given by
!  the following equation: 
!  psia(iens)=psif(iens)+{[H*(A-Mean)*(A-Mean)^T]^T*scale}*P^(-1)*h
!  H is the observation matrix, A contains the ensemble, Mean contains
!  the ensemble mean, h is the residual between the observations and 
!  the ensemble members at the position of the observations.

   use mod_dimensions
   use mod_states
   use mod_observations
   use m_random
   implicit none

! Interface
   integer,            intent(in) :: nrsamp
   type(observations), intent(in) :: obs(nrobs)
   type(states),       intent(inout) :: mem(nrsamp)

! Locals
   integer iens, i, j, k, m, nrsigma, obstot
   type(observations), allocatable :: d(:)      ! The active observations
   type(observations), allocatable :: tmpobs(:) ! Observations with noise

   type(states) Mean                            ! ensemble mean

   type(states), allocatable :: tmpens(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   real,    allocatable  :: h(:,:), E(:,:), d_vec(:)
   real,    allocatable  :: tmp(:,:)
   real,    allocatable  :: P(:,:)
 
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! For eigen value decomposition
   real,    allocatable  :: ZZ(:,:)  !contains the eigenvectors of P
   real,    allocatable  :: w(:)     !dsyevx workspace
   real,    allocatable  :: ss(:)    !contains the eigenvalues of P
   integer, allocatable  :: iw(:)    !dsyevx workspace
   integer, allocatable  :: ipvt(:)  !dsyevx workspace
   integer                  info     !completion code for dsyevx 
   real                     abstol   !The abs error tolerance for the eigenvalues
   real,   allocatable   :: ap(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   real,    allocatable  :: lambda(:)
   integer                  nndim      !parameter in the dgemm subroutine
   type(states), allocatable :: rtd(:) !Damped Transpose of Representer matrix
   integer                  ix,iy,iz,m2,m3,idim !temporary counters
   integer                  toteig   !total number of eigenvectors found from subroutine dsyevx
   real                     VL, VU
   integer                  IL, IU
   integer mp_my_threadnum, mp_numthreads
   real                     summ, summ2
   real,    allocatable  :: L(:,:)
   real                     scale      !1/nrsamp
   character(len=4) tag4
   character(len=3) tag3
   character(len=7) tag7

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute ensemble mean
   scale = 1.0/float(nrsamp)
   Mean=0.0
   do iens=1,nrsamp
      Mean = Mean + mem(iens)
   enddo
   Mean=scale*Mean

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Pick out active observations
   obstot=0
   do m=1,nrobs
      if (obs(m)%action) obstot=obstot+1
   enddo
   print *,'obstot= ',obstot
   allocate(d(obstot))
   obstot=0
   do m=1,nrobs
      if (obs(m)%action) then
         obstot=obstot+1
         d(obstot)=obs(m)
      endif
   enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute the residual between the observations and the value of the
! ensemble members at the position of the observations.
!        "h(:,iens)=d(:)%d-H*mem(iens)"
! and the matrix  "tmp = H*(A(iens) - Mean)" which are used in the analysis.

   allocate(h(obstot,nrsamp))
   allocate(E(obstot,nrsamp))
   allocate(tmp(nrsamp,obstot))  !GE used instead of allocate(tmp(obstot,nrsamp))
   allocate(tmpobs(obstot))
   allocate(d_vec(obstot))

! Compute non-perturbed innovations
   do m=1, obstot
      if (trim(d(m)%ch)=='N') then
         d_vec(m) = d(m)%d - Mean%N(d(m)%k)
      elseif (trim(d(m)%ch)=='P') then
         d_vec(m) = d(m)%d - Mean%P(d(m)%k)
      elseif (trim(d(m)%ch)=='H') then
         d_vec(m) = d(m)%d - Mean%H(d(m)%k)
      else
         stop 'analysis: error computing d_vec'
      endif
   enddo

   do iens=1,nrsamp

! Add noise to the observations 
      call random(tmpobs(:)%d,obstot) 
      call random(E(:,iens),obstot) 

      do m=1,obstot
         tmpobs(m)%ch=d(m)%ch
         tmpobs(m)%d = d(m)%d + sqrt(d(m)%var)*tmpobs(m)%d
         tmpobs(m)%d=max(0.0,tmpobs(m)%d)  ! ME: Pga neg. obs.
         !E(m,iens) = tmpobs(m)%d - d(m)%d
         E(m,iens) =  E(m,iens) * sqrt(d(m)%var)
      enddo

      do m=1,obstot
! h(:,iens) = obs(:)%d - H*mem(iens)
         if (trim(d(m)%ch)=='N') then
            h(m,iens)=tmpobs(m)%d-mem(iens)%N(d(m)%k)  
         elseif (trim(d(m)%ch)=='P') then
            h(m,iens)=tmpobs(m)%d-mem(iens)%P(d(m)%k)  
         elseif (trim(d(m)%ch)=='H') then
            h(m,iens)=tmpobs(m)%d-mem(iens)%H(d(m)%k)  
         else
            stop 'analysis: errror computing h'
         endif

! tmp = H*(mem(iens) - Mean)
         if (trim(tmpobs(m)%ch)=='N') then
            tmp(iens,m)=mem(iens)%N(d(m)%k)  - Mean%N(d(m)%k)
         elseif (trim(tmpobs(m)%ch)=='P') then
            tmp(iens,m)=mem(iens)%P(d(m)%k)  - Mean%P(d(m)%k)
         elseif (trim(tmpobs(m)%ch)=='H') then
            tmp(iens,m)=mem(iens)%H(d(m)%k)  - Mean%H(d(m)%k)
         else
            stop 'analysis: error computing tmp'
         endif

       enddo
   enddo

   print *,'tmp= done'

   ! CAll standard analysis
    !call analysis(mem,h,E,transpose(tmp),3*kdim,nrsamp,obstot,.false.)
    !return


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The observations have uncorrelated data errors
   allocate(P(obstot,obstot))
   P=0.0
   do m=1,obstot
      P(m,m)=d(m)%var
!      print '(4f9.4)', P(m,:)
   enddo

   call analysis2(mem,h,P,transpose(tmp),3*kdim,nrsamp,obstot,.false.)
   return

   !call analysis4(mem,P,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
   !return
  
   !call analysis4c(mem,P,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
   !return

  ! call analysis5(mem,P,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
  ! return
 
    !call analysis5c(mem,P,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
    !return

   !call analysis6(mem,E,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
   !return

   !call analysis6c(mem,E,transpose(tmp),d_vec,3*kdim,nrsamp,obstot,.false.)
   !return
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Construct the P matrix "P + scale*H*P^f*H"
! defined as:  P=scale*matmul(tmp,transpose(tmp))+P
! Here a rewriting is used to save compiler induced memory allocation
! of temporary variables.  Note that tmp has already been transposed
! to access memory with stride one.  The alternative is depressing!

   do j=1,obstot
   do i=j,obstot
      P(i,j)=scale*dot_product(tmp(:,i),tmp(:,j))+P(i,j)
    
      P(j,i)  = P(i,j)
   enddo
   enddo
   print *,'P ok'
   print '(a,i8)','Upper left of P'
 !  print '(10g11.2)',P(1:min(10,obstot),1:min(10,obstot))
   print '(4g11.2)',P(1:4,1:4)
   print *


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Calculate the eigendecomposition of P using eispack. P = ZZ*ss*ZZ

   allocate(ZZ(obstot,obstot))
   allocate(w(8*obstot))
   allocate(ss(obstot))
   allocate(iw(5*obstot))
   allocate(ipvt(obstot))
   allocate(lambda(obstot))

   ZZ=0.0; w=0.0; ss=0.0; iw=0; ipvt=0; lambda=0.0



! Eigen value factorization using rsm:
!  call rsm(obstot,obstot,P,ss,obstot,ZZ,w,ipvt,info)


   allocate (ap(obstot*(obstot+1)/2) )
   k=0
   do j=1,obstot
   do i=1,j
      k=k+1
    ap(k)=P(i,j) 
    enddo
    enddo



   call dspev(21,ap,ss,ZZ,obstot,obstot,w,2*obstot)
   deallocate(ap)



   print *,'Eigenvalues of P:'
   print '(20g10.2)',ss(obstot:1:-1)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Write the eigenvalues to a file eig.dat
   summ=sum(ss)
   print *,'sum(eigenvalues)=',summ

   open(10,file='eig.dat')
      do i=obstot,1,-1
!         write(10,'(i4,g13.5)')i,ss(i)/summ
         write(10,'(i4,g13.5)')obstot-i+1,ss(i)/summ
      enddo
   close(10)
   print *,'eig.dat done'


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Find the number of significant eigenvalues 
!   To avoid a singular invers of P=ZZ*ss*ZZ and a noisy solution, we
!   truncate the spectrum of P. lambda contains the truncated inverse of
!   the matrix diag(ss)

   print *,'computing the number of significant eigen values'
   summ2=0.0
   do i=obstot,1,-1
      summ2=summ2+ss(i)
      if (summ2/summ>0.99) then
         nrsigma=obstot-i+1
         exit
      endif
   enddo
   lambda=0.0
   do k=obstot,nrsigma,-1
      lambda(k)=1.0/ss(k)
   enddo
   print *,'Number of significant eigenvalues are: ',nrsigma


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The analysis increment is given as 
!     scale*(A-Mean)*[H*(A-Mean)]^T*{ZZ*lambda*ZZ^T*h}
! which also can be written as
!     scale*(A-Mean)*tmp^T*{ZZ*lambda*ZZ^T*h}
! or
!     scale*(A-Mean)*L

! ZZ^T*h  (rotated measurement residuals, or residuals represented in the ZZ space)
   h=matmul(transpose(ZZ),h)
   print *,'made ZZ^T*h'

! lambda*(ZZ^T*h)
   do j=1,nrsamp
      h(:,j)=lambda(:)*h(:,j)
   enddo
   print *,'made lambda*(ZZ^T*h)='
   print '(15E10.2)',h(:,2)

! ZZ*(lambda*(ZZ^T*h))
   h=matmul(ZZ,h)
   print *,'made ZZ*(lambda*(ZZ^T*h))'

! tmp^T*(ZZ*lambda*ZZ^T*h)            ! L: Not in use now!
   allocate(L(nrsamp,nrsamp))
   L=matmul(tmp,h)

   print *,'made coefficients L=tmp^T*(ZZ*lambda*ZZ^T*h)'


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Final analysis (A-Mean)*L. 
! ndim = nx*ny*nz*(number of 3D fields)+nx*ny*(number of 2D fields)
   nndim=kdim*antvar

   allocate(tmpens(nrsamp));  print *,'allocated tmpens'

!$OMP PARALLEL DO PRIVATE(iens)
   do iens=1,nrsamp
      tmpens(iens)=mem(iens)-Mean
   enddo
   

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Analysis using representer approach

! - compute represeters, ``rtd = scale*tmpens*tmp''.
   allocate (rtd(obstot))

   do i = 1, obstot
      rtd(i)=0.0 ! Initialize representers
   enddo
   print *,'Computing the representer matrix rtd ...'
   call dgemm('n','n',nndim,obstot,nrsamp,scale,tmpens,nndim,tmp,nrsamp,1.0,rtd,nndim)
   print *,'... done representers'

   deallocate (tmpens)
   deallocate (tmp)

! - compute analysis ``A = rtd*h + A''.

   print *,'Computing the new analysis ...'
   call dgemm('n','n',nndim,nrsamp,obstot,1.0,rtd,nndim,h,obstot,1.0,mem,nndim)
   print *,'... done'
   
   deallocate (rtd)

   do iens=1,nrsamp
     mem(iens)%N=max(0.0,mem(iens)%N) 
     mem(iens)%P=max(0.0,mem(iens)%P)
     mem(iens)%H=max(0.0,mem(iens)%H)
   enddo

end subroutine prepare_and_analysis

end module m_analysis