/[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

trunk/filtrez/inifilr.f90 revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/Sources/filtrez/inifilr.f revision 140 by guez, Fri Jun 5 18:58:06 2015 UTC
# Line 1  Line 1 
1  module inifilr_m  module inifilr_m
2    
   use dimens_m, only: iim  
   
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv    INTEGER jfiltnu, jfiltsu, jfiltnv, jfiltsv
6    INTEGER, PARAMETER:: nfilun=3, nfilus=2, nfilvn=2, nfilvs=2    ! jfiltn index of the last scalar line filtered in NH
7      ! jfilts index of the first line filtered in SH
8    
9      ! North:
10      real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :)
11      ! (iim, iim, 2:jfiltnu)
12    
13      real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
14    
15    real matriceun(iim,iim,nfilun), matriceus(iim,iim,nfilus)    ! South:
16    real matricevn(iim,iim,nfilvn), matricevs(iim,iim,nfilvs)    real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
17    real matrinvn(iim,iim,nfilun), matrinvs(iim,iim,nfilus)    ! (iim, iim, jfiltsu:jjm)
18    
19    private iim, nfilun, nfilus, nfilvn, nfilvs    real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
20    
21  contains  contains
22    
# Line 21  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)  
   
     ! The modes are filtered from modfrst to modemax.  
31    
32      USE dimens_m, ONLY : iim, jjm      USE dimens_m, ONLY : iim, jjm
33      use conf_gcm_m, ONLY : fxyhypb, ysinus      USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
34      USE comgeom, ONLY : rlatu, rlatv, xprimu      use inifgn_m, only: inifgn, eignfnu, eignfnv
35      use nr_util, only: pi      use nr_util, only: pi
     USE serre, ONLY : alphax  
     USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &  
          eignfnv, modfrstu, modfrstv  
36    
37      ! Local:      ! Local:
38      REAL dlonu(iim), 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, dxmin, colat0      REAL dymin, colat0
44      REAL eignft(iim, iim), coff      REAL eignft(iim, iim), coff
45      EXTERNAL inifgn  
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"
56    
     DO i = 1, iim  
        dlonu(i) = xprimu(i)  
     END DO  
   
57      CALL inifgn(eignvl)      CALL inifgn(eignvl)
58    
59      PRINT *, 'EIGNVL '      PRINT *, 'EIGNVL '
# Line 76  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    
     dxmin = dlonu(1)  
     DO i = 2, iim  
        dxmin = min(dxmin, dlonu(i))  
     END DO  
76      dymin = dlatu(1)      dymin = dlatu(1)
77      DO j = 2, jjm      DO j = 2, jjm
78         dymin = min(dymin, dlatu(j))         dymin = min(dymin, dlatu(j))
79      END DO      END DO
80    
81      colat0 = min(0.5, dymin/dxmin)      colat0 = min(0.5, dymin / minval(xprimu(:iim)))
   
     IF (.NOT. fxyhypb .AND. ysinus) THEN  
        colat0 = 0.6  
        ! À revoir pour ysinus  
        alphax = 0.  
     END IF  
82    
83      PRINT *, 'colat0 = ', colat0      PRINT *, 'colat0 = ', colat0
     PRINT *, 'alphax = ', alphax  
84    
85      IF (alphax == 1.) THEN      lamdamax = iim / (pi * colat0 / grossismx)
        PRINT *, 'alphax doit etre < a 1. Corriger '  
        STOP 1  
     END IF  
   
     lamdamax = iim / (pi * colat0 * (1. - alphax))  
86      rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))      rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))
87    
88      DO j = 1, jjm      DO j = 1, jjm
# Line 118  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 157  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 179  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 255  contains Line 222  contains
222      PRINT *, 'Modes premiers u '      PRINT *, 'Modes premiers u '
223      PRINT 334, modfrstu      PRINT 334, modfrstu
224    
225      IF (nfilun < jfiltnu) THEN      allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
226         PRINT *, 'le parametre nfilun utilise pour la matrice ', &      allocate(matricevn(iim, iim, jfiltnv))
227              'matriceun est trop petit ! '      allocate(matricevs(iim, iim, jfiltsv:jjm))
228         PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnu      allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
        PRINT *, 'Pour information, nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', jfiltnu, &  
             jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
        STOP 1  
     END IF  
     IF (nfilun > jfiltnu+2) THEN  
        PRINT *, 'le parametre nfilun utilise pour la matrice ', &  
             'matriceun est trop grand ! Gachis de memoire ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnu  
        PRINT *, 'Pour information, nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
     END IF  
     IF (nfilus < jjm-jfiltsu+1) THEN  
        PRINT *, 'le parametre nfilus utilise pour la matrice ', &  
             'matriceus est trop petit ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', &  
             jjm - jfiltsu + 1  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
        STOP 1  
     END IF  
     IF (nfilus > jjm-jfiltsu+3) THEN  
        PRINT *, 'le parametre nfilus utilise pour la matrice ', &  
             'matriceus est trop grand ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', &  
             jjm - jfiltsu + 1  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
     END IF  
     IF (nfilvn < jfiltnv) THEN  
        PRINT *, 'le parametre nfilvn utilise pour la matrice ', &  
             'matricevn est trop petit ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnv  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
        STOP 1  
     END IF  
     IF (nfilvn > jfiltnv+2) THEN  
        PRINT *, 'le parametre nfilvn utilise pour la matrice ', &  
             'matricevn est trop grand ! Gachis de memoire ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', jfiltnv  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
     END IF  
     IF (nfilvs < jjm-jfiltsv+1) THEN  
        PRINT *, 'le parametre nfilvs utilise pour la matrice ', &  
             'matricevs est trop petit ! Le changer dans parafilt.h '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', &  
             jjm - jfiltsv + 1  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
        STOP 1  
     END IF  
     IF (nfilvs > jjm-jfiltsv+3) THEN  
        PRINT *, 'le parametre nfilvs utilise pour la matrice ', &  
             'matricevs est trop grand ! Gachis de memoire ! '  
        PRINT *, 'Le changer dans parafilt.h et le mettre a ', &  
             jjm - jfiltsv + 1  
        PRINT *, 'Pour information , nfilun, nfilus, nfilvn, nfilvs ', &  
             'doivent etre egaux successivement a ', &  
             jfiltnu, jjm - jfiltsu + 1, jfiltnv, jjm - jfiltsv + 1  
     END IF  
229    
230      ! Calcul de la matrice filtre 'matriceu' pour les champs situes      ! Calcul de la matrice filtre 'matriceu' pour les champs situes
231      ! sur la grille scalaire      ! sur la grille scalaire
# Line 338  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 352  contains Line 251  contains
251            end IF            end IF
252            eignft(i, :) = eignfnv(:, i) * coff            eignft(i, :) = eignfnv(:, i) * coff
253         END DO         END DO
254         matriceus(:, :, j - jfiltsu + 1) = matmul(eignfnv, eignft)         matriceus(:, :, j) = matmul(eignfnv, eignft)
255      END DO      END DO
256    
257      ! Calcul de la matrice filtre 'matricev' pour les champs situes      ! Calcul de la matrice filtre 'matricev' pour les champs situes
# Line 365  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 377  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-jfiltsv+1) = matmul(eignfnu, eignft)         matricevs(:, :, j) = matmul(eignfnu, eignft)
282      END DO      END DO
283    
284      ! Calcul de la matrice filtre 'matrinv' pour les champs situes      ! Calcul de la matrice filtre 'matrinv' pour les champs situes
# Line 390  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 402  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-jfiltsu+1) = matmul(eignfnv, eignft)         matrinvs(:, :, j) = matmul(eignfnv, eignft)
309      END DO      END DO
310    
311  334 FORMAT (1X, 24I3)  334 FORMAT (1X, 24I3)
Legend:
Removed from v.76  
changed lines
  Added in v.140
  ViewVC Help
Powered by ViewVC 1.1.21