Blame view

analysis2.F90 4.6 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
subroutine analysis2(A, D, R, S, ndim, nrens, nrobs, verbose)
! Computes the analysed ensemble for A 
   use m_multa
   implicit none
   integer, intent(in) :: ndim             ! dimension of model state
   integer, intent(in) :: nrens            ! number of ensemble members
   integer, intent(in) :: nrobs            ! number of observations
   
   real, intent(inout) :: A(ndim,nrens)   ! ensemble matrix
   real, intent(in)    :: D(nrobs,nrens)  ! matrix  holding observation innovations
   real, intent(in)    :: S(nrobs,nrens)  ! matrix holding HA' 
   real, intent(inout) :: R(nrobs,nrobs)  ! Error covariance matrix for observations
   logical, intent(in) :: verbose


   real, allocatable, dimension(:,:) :: X1,X2,U,X4,Reps,V
   real X3(nrobs,nrens)
   real, allocatable, dimension(:)   :: sig,work

   real sigsum,sigsum1,oneobs(1,1)
   integer ierr,nrsigma,i,j,lwork, m
   integer iblkmax
   character(len=2) tag2



   if (nrobs > 1) then
!      R=float(nrens)*R+matmul(S,transpose(S))
      call dgemm('n','t', nrobs, nrobs, nrens, &
                 1.0, S, nrobs, &
                      S, nrobs, &
        float(nrens), R, nrobs)



      allocate (U(nrobs,nrobs)  )
      allocate (V(nrobs,nrobs)  )
      allocate (sig(nrobs)  )
      lwork=2*max(3*nrobs+nrobs,5*nrobs)
      allocate(work(lwork))
      sig=0.0
      call dgesvd('A', 'A', nrobs, nrobs, R, nrobs, sig, U, nrobs, V, nrobs, work, lwork, ierr)
      deallocate(work)
      if (ierr /= 0) then
         print *,'ierr from call dgesvd= ',ierr
         stop
      endif

      open(10,file='sigma.dat')
         do i=1,nrobs
            write(10,'(i5,g12.3)')i,sig(i)
         enddo
      close(10)

      sigsum=sum( sig(1:nrobs) )
      sigsum1=0.0
   ! Significant eigenvalues.
      nrsigma=0
      do i=1,nrobs                 ! singular values are in descending order
         if (sigsum1/sigsum < 0.999) then
            nrsigma=nrsigma+1
            sigsum1=sigsum1+sig(i)
            sig(i) = 1.0/sig(i)
         else
            sig(i:nrobs)=0.0
            exit
         endif
      enddo

      if (verbose) then
         write(*,'(a,i5,g13.5)') ' dominant sing. values and share ',nrsigma,sigsum1/sigsum
         write(*,'(8g12.3)')1./sig(1:nrsigma), sig(nrsigma+1:nrobs)
      endif

      allocate (X1(nrobs,nrobs))
      do i=1,nrobs
      do j=1,nrobs
         X1(i,j)=sig(i)*U(j,i)
      enddo
      enddo
      deallocate(sig)

      allocate (X2(nrobs,nrens))
!     X2=matmul(X1,D)
      call dgemm('n','n',nrobs,nrens,nrobs,1.0,X1,nrobs,D ,nrobs,0.0,X2,nrobs)
      deallocate(X1) 

!     X3=matmul(V,X2)
      call dgemm('t','n',nrobs,nrens,nrobs,1.0,V ,nrobs,X2,nrobs,0.0,X3,nrobs)
      deallocate(V)
      deallocate(X2)

   else
      oneobs=matmul(S,transpose(S))+R
      print *,'oneobs: ',oneobs(1,1)
      X3=D/oneobs(1,1)
   endif

   if (2_8*ndim*nrobs < 1_8*nrens*(nrobs+ndim)) then
!    Code for few observations ( m<nN/(2n-N) )
      if (verbose) print * ,'analysis: Representer approach is used'
      allocate (Reps(ndim,nrobs))

!    Reps=matmul(A,transpose(S))
      call dgemm('n','t',ndim,nrobs,nrens,1.0,A,ndim,S,nrobs,0.0,Reps,ndim)
!      call printreps(Reps,nrobs)

!    A=A+matmul(Reps,X3)
      call dgemm('n','n',ndim,nrens,nrobs,1.0,Reps,ndim,X3,nrobs,1.0,A,ndim)
      deallocate(Reps)

      tag2(1:2)='X3'
      open(10,file='X.uf',form='unformatted')
         write(10)tag2,nrens,nrobs,X3,S
      close(10)

   else
      if (verbose) print * ,'analysis: X5 appraoch is used'
      allocate(X4(nrens,nrens))
!      X4=matmul(transpose(S),X3)
      call dgemm('t','n',nrens,nrens,nrobs,1.0,S,nrobs,X3,nrobs,0.0,X4,nrens)
      do i=1,nrens
         X4(i,i)=X4(i,i)+1.0  ! Addition with Af
      enddo

      iblkmax=min(ndim,200)
      call  multa(A, X4, ndim, nrens, iblkmax )

      tag2(1:2)='X5'
      open(10,file='X.uf',form='unformatted')
         write(10)tag2,nrens,X4
      close(10)

      open(10,file='X5.dat')
         write(10,*)'TITLE = "X5"'
         write(10,*)'VARIABLES = "iens_f" "iens_a" "X5"'
         write(10,'(2(a,i5),a)')' ZONE  F=BLOCK, I=',nrens,', J=',nrens,', K=1'
         write(10,'(30I5)')((i,i=1,nrens),j=1,nrens)
         write(10,'(30I5)')((j,i=1,nrens),j=1,nrens)
         write(10,'(10(1x,e12.5))')((X4(i,j),i=1,nrens),j=1,nrens)
      close(10)

      open(10,file='X5col.dat')
         write(10,'(a)')'These should all be ones'
         do j=1,nrens
            write(10,'(i5,f10.4)')j,sum(X4(:,j))  
         enddo
      close(10)

      open(10,file='X5row.dat')
         do j=1,nrens
            write(10,'(i5,f10.4)')j,sum(X4(j,:))/float(nrens)
          enddo
      close(10)

      deallocate(X4)
   endif 


end subroutine analysis2