/[lmdze]/trunk/Sources/filtrez/inifilr.f
ViewVC logotype

Diff of /trunk/Sources/filtrez/inifilr.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 139 by guez, Tue May 26 17:46:03 2015 UTC revision 140 by guez, Fri Jun 5 18:58:06 2015 UTC
# Line 3  module inifilr_m Line 3  module inifilr_m
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv    INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
6      ! jfiltn index of the last scalar line filtered in NH
7      ! jfilts index of the first line filtered in SH
8    
9    ! North:    ! North:
10    real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)    real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
# Line 24  contains Line 26  contains
26      ! H. Upadhyaya, O. Sharma      ! H. Upadhyaya, O. Sharma
27    
28      ! This routine computes the eigenfunctions of the laplacian on the      ! This routine computes the eigenfunctions of the laplacian on the
29      ! stretched grid, and the filtering coefficients.      ! stretched grid, and the filtering coefficients. The modes are
30      ! We designate:      ! filtered from modfrst to iim.
     ! eignfn eigenfunctions of the discrete laplacian  
     ! eigenvl eigenvalues  
     ! jfiltn index of the last scalar line filtered in NH  
     ! jfilts index of the first line filtered in SH  
     ! modfrst index of the mode from where modes are filtered  
     ! modemax maximum number of modes (im)  
     ! coefil filtering coefficients (lamda_max * cos(rlat) / lamda)  
     ! sdd SQRT(dx)  
31    
     ! The modes are filtered from modfrst to modemax.  
   
     USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &  
          eignfnv, modfrstu, modfrstv  
32      USE dimens_m, ONLY : iim, jjm      USE dimens_m, ONLY : iim, jjm
33      USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx      USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
34      use inifgn_m, only: inifgn      use inifgn_m, only: inifgn, eignfnu, eignfnv
35      use nr_util, only: pi      use nr_util, only: pi
36    
37      ! Local:      ! Local:
38      REAL dlatu(jjm)      REAL dlatu(jjm)
39      REAL rlamda(2: iim), eignvl(iim)      REAL rlamda(2: iim)
40        real eignvl(iim) ! eigenvalues
41      REAL lamdamax, cof      REAL lamdamax, cof
42      INTEGER i, j, modemax, imx, k, kf      INTEGER i, j, k, kf
43      REAL dymin, colat0      REAL dymin, colat0
44      REAL eignft(iim, iim), coff      REAL eignft(iim, iim), coff
45    
46        ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
47        real coefilu(iim, jjm), coefilv(iim, jjm)
48        real coefilu2(iim, jjm), coefilv2(iim, jjm)
49    
50        integer modfrstu(jjm), modfrstv(jjm)
51        ! index of the mode from where modes are filtered
52    
53      !-----------------------------------------------------------      !-----------------------------------------------------------
54    
55      print *, "Call sequence information: inifilr"      print *, "Call sequence information: inifilr"
# Line 73  contains Line 70  contains
70      ! Calcul de colat0      ! Calcul de colat0
71    
72      DO j = 1, jjm      DO j = 1, jjm
73         dlatu(j) = rlatu(j) - rlatu(j+1)         dlatu(j) = rlatu(j) - rlatu(j + 1)
74      END DO      END DO
75    
76      dymin = dlatu(1)      dymin = dlatu(1)
# Line 99  contains Line 96  contains
96    
97      ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv      ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
98    
99      modemax = iim      PRINT *, 'TRUNCATION AT ', iim
     imx = iim  
   
     PRINT *, 'TRUNCATION AT ', imx  
100    
101      DO j = 2, jjm / 2 + 1      DO j = 2, jjm / 2 + 1
102         IF (cos(rlatu(j)) / colat0 < 1. &         IF (cos(rlatu(j)) / colat0 < 1. &
103              .and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j              .and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j
104    
105         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
106              .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) &              .and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &
107              jfiltsu = jjm - j + 2              jfiltsu = jjm - j + 2
108      END DO      END DO
109    
110      DO j = 1, jjm/2      DO j = 1, jjm / 2
111         cof = cos(rlatv(j))/colat0         IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &
112         IF (cof < 1.) THEN              jfiltnv = j
113            IF (rlamda(imx)*cos(rlatv(j)) < 1.) jfiltnv = j  
114         END IF         IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &
115                * cos(rlatv(jjm - j + 1)) < 1.) jfiltsv = jjm - j + 1
        cof = cos(rlatv(jjm-j+1))/colat0  
        IF (cof < 1.) THEN  
           IF (rlamda(imx)*cos(rlatv(jjm-j+1)) < 1.) jfiltsv = jjm - j + 1  
        END IF  
116      END DO      END DO
117    
118      IF (jfiltnu <= 0) jfiltnu = 1      IF (jfiltnu <= 0) jfiltnu = 1
119      IF (jfiltnu > jjm/2+1) THEN      IF (jfiltnu > jjm / 2 + 1) THEN
120         PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu         PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
121         STOP 1         STOP 1
122      END IF      END IF
# Line 138  contains Line 128  contains
128      END IF      END IF
129    
130      IF (jfiltnv <= 0) jfiltnv = 1      IF (jfiltnv <= 0) jfiltnv = 1
131      IF (jfiltnv > jjm/2) THEN      IF (jfiltnv > jjm / 2) THEN
132         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
133         STOP 1         STOP 1
134      END IF      END IF
# Line 160  contains Line 150  contains
150      END DO      END DO
151    
152      DO j = 2, jfiltnu      DO j = 2, jfiltnu
153         DO k = 2, modemax         DO k = 2, iim
154            cof = rlamda(k) * cos(rlatu(j))            IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
           IF (cof < 1.) exit  
155         end DO         end DO
156         if (k == modemax + 1) cycle         if (k == iim + 1) cycle
157         modfrstu(j) = k         modfrstu(j) = k
158    
159         kf = modfrstu(j)         kf = modfrstu(j)
160         DO k = kf, modemax         DO k = kf, iim
161            cof = rlamda(k)*cos(rlatu(j))            cof = rlamda(k) * cos(rlatu(j))
162            coefilu(k, j) = cof - 1.            coefilu(k, j) = cof - 1.
163            coefilu2(k, j) = cof*cof - 1.            coefilu2(k, j) = cof**2 - 1.
164         end DO         end DO
165      END DO      END DO
166    
167      DO j = 1, jfiltnv      DO j = 1, jfiltnv
168         DO k = 2, modemax         DO k = 2, iim
169            cof = rlamda(k)*cos(rlatv(j))            IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
           IF (cof < 1.) exit  
170         end DO         end DO
171         if (k == modemax + 1) cycle         if (k == iim + 1) cycle
172         modfrstv(j) = k         modfrstv(j) = k
173    
174         kf = modfrstv(j)         kf = modfrstv(j)
175         DO k = kf, modemax         DO k = kf, iim
176            cof = rlamda(k)*cos(rlatv(j))            cof = rlamda(k) * cos(rlatv(j))
177            coefilv(k, j) = cof - 1.            coefilv(k, j) = cof - 1.
178            coefilv2(k, j) = cof*cof - 1.            coefilv2(k, j) = cof**2 - 1.
179         end DO         end DO
180      end DO      end DO
181    
182      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
183         DO k = 2, modemax         DO k = 2, iim
184            cof = rlamda(k)*cos(rlatu(j))            IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
           IF (cof < 1.) exit  
185         end DO         end DO
186         if (k == modemax + 1) cycle         if (k == iim + 1) cycle
187         modfrstu(j) = k         modfrstu(j) = k
188    
189         kf = modfrstu(j)         kf = modfrstu(j)
190         DO k = kf, modemax         DO k = kf, iim
191            cof = rlamda(k)*cos(rlatu(j))            cof = rlamda(k) * cos(rlatu(j))
192            coefilu(k, j) = cof - 1.            coefilu(k, j) = cof - 1.
193            coefilu2(k, j) = cof*cof - 1.            coefilu2(k, j) = cof**2 - 1.
194         end DO         end DO
195      end DO      end DO
196    
197      DO j = jfiltsv, jjm      DO j = jfiltsv, jjm
198         DO k = 2, modemax         DO k = 2, iim
199            cof = rlamda(k)*cos(rlatv(j))            IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
           IF (cof < 1.) exit  
200         end DO         end DO
201         if (k == modemax + 1) cycle         if (k == iim + 1) cycle
202         modfrstv(j) = k         modfrstv(j) = k
203    
204         kf = modfrstv(j)         kf = modfrstv(j)
205         DO k = kf, modemax         DO k = kf, iim
206            cof = rlamda(k)*cos(rlatv(j))            cof = rlamda(k) * cos(rlatv(j))
207            coefilv(k, j) = cof - 1.            coefilv(k, j) = cof - 1.
208            coefilv2(k, j) = cof*cof - 1.            coefilv2(k, j) = cof**2 - 1.
209         end DO         end DO
210      END DO      END DO
211    
212      IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN      IF (jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2) THEN
213         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
214         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
215    
# Line 251  contains Line 237  contains
237            else            else
238               coff = coefilu(i, j)               coff = coefilu(i, j)
239            end IF            end IF
240            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
241         END DO         END DO
242         matriceun(:, :, j) = matmul(eignfnv, eignft)         matriceun(:, :, j) = matmul(eignfnv, eignft)
243      END DO      END DO
# Line 278  contains Line 264  contains
264            else            else
265               coff = coefilv(i, j)               coff = coefilv(i, j)
266            end IF            end IF
267            eignft(i, :) = eignfnu(:, i)*coff            eignft(i, :) = eignfnu(:, i) * coff
268         END DO         END DO
269         matricevn(:, :, j) = matmul(eignfnu, eignft)         matricevn(:, :, j) = matmul(eignfnu, eignft)
270      END DO      END DO
# Line 290  contains Line 276  contains
276            else            else
277               coff = coefilv(i, j)               coff = coefilv(i, j)
278            end IF            end IF
279            eignft(i, :) = eignfnu(:, i)*coff            eignft(i, :) = eignfnu(:, i) * coff
280         END DO         END DO
281         matricevs(:, :, j) = matmul(eignfnu, eignft)         matricevs(:, :, j) = matmul(eignfnu, eignft)
282      END DO      END DO
# Line 303  contains Line 289  contains
289            IF (i < modfrstu(j)) then            IF (i < modfrstu(j)) then
290               coff = 0.               coff = 0.
291            else            else
292               coff = coefilu(i, j)/(1.+coefilu(i, j))               coff = coefilu(i, j) / (1. + coefilu(i, j))
293            end IF            end IF
294            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
295         END DO         END DO
296         matrinvn(:, :, j) = matmul(eignfnv, eignft)         matrinvn(:, :, j) = matmul(eignfnv, eignft)
297      END DO      END DO
# Line 315  contains Line 301  contains
301            IF (i < modfrstu(j)) then            IF (i < modfrstu(j)) then
302               coff = 0.               coff = 0.
303            else            else
304               coff = coefilu(i, j)/(1.+coefilu(i, j))               coff = coefilu(i, j) / (1. + coefilu(i, j))
305            end IF            end IF
306            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
307         END DO         END DO
308         matrinvs(:, :, j) = matmul(eignfnv, eignft)         matrinvs(:, :, j) = matmul(eignfnv, eignft)
309      END DO      END DO
Legend:
Removed from v.139  
changed lines
  Added in v.140
  ViewVC Help
Powered by ViewVC 1.1.21