/[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.f revision 82 by guez, Wed Mar 5 14:57:53 2014 UTC trunk/Sources/filtrez/inifilr.f revision 143 by guez, Tue Jun 9 14:32:46 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 jumble, only: new_unit
36      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  
37    
38      ! Local:      ! Local:
39      REAL dlonu(iim), dlatu(jjm)      REAL dlatu(jjm)
40      REAL rlamda(2: iim), eignvl(iim)      REAL rlamda(2: iim)
41        real eignvl(iim) ! eigenvalues
42      REAL lamdamax, cof      REAL cof
43      INTEGER i, j, modemax, imx, k, kf      INTEGER i, j, k, kf, unit
44      REAL dymin, dxmin, colat0      REAL colat0
45      REAL eignft(iim, iim), coff      REAL eignft(iim, iim), coff
46      EXTERNAL inifgn  
47        ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
48        real coefilu(iim, jjm), coefilv(iim, jjm)
49        real coefilu2(iim, jjm), coefilv2(iim, jjm)
50    
51        integer modfrstu(jjm), modfrstv(jjm)
52        ! index of the mode from where modes are filtered
53    
54      !-----------------------------------------------------------      !-----------------------------------------------------------
55    
56      print *, "Call sequence information: inifilr"      print *, "Call sequence information: inifilr"
57    
     DO i = 1, iim  
        dlonu(i) = xprimu(i)  
     END DO  
   
58      CALL inifgn(eignvl)      CALL inifgn(eignvl)
59    
60      PRINT *, 'EIGNVL '      call new_unit(unit)
61      PRINT "(1X, 5E13.6)", eignvl      open(unit, file = "eignvl.txt", status = "replace", action = "write")
62        write(unit, fmt = *) EIGNVL
63        close(unit)
64    
65      ! compute eigenvalues and eigenfunctions      ! compute eigenvalues and eigenfunctions
66      ! compute the filtering coefficients for scalar lines and      ! compute the filtering coefficients for scalar lines and
# Line 74  contains Line 71  contains
71      ! is set equal to zero for the regular grid case      ! is set equal to zero for the regular grid case
72    
73      ! Calcul de colat0      ! Calcul de colat0
74        forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
75      DO j = 1, jjm      colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
        dlatu(j) = rlatu(j) - rlatu(j+1)  
     END DO  
   
     dxmin = dlonu(1)  
     DO i = 2, iim  
        dxmin = min(dxmin, dlonu(i))  
     END DO  
     dymin = dlatu(1)  
     DO j = 2, jjm  
        dymin = min(dymin, dlatu(j))  
     END DO  
   
     colat0 = min(0.5, dymin/dxmin)  
   
     IF (.NOT. fxyhypb .AND. ysinus) THEN  
        colat0 = 0.6  
        ! À revoir pour ysinus  
        alphax = 0.  
     END IF  
   
76      PRINT *, 'colat0 = ', colat0      PRINT *, 'colat0 = ', colat0
     PRINT *, 'alphax = ', alphax  
   
     IF (alphax == 1.) THEN  
        PRINT *, 'alphax doit etre < a 1. Corriger '  
        STOP 1  
     END IF  
   
     lamdamax = iim / (pi * colat0 * (1. - alphax))  
     rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))  
77    
78      DO j = 1, jjm      rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
79         DO i = 1, iim      coefilu = 0.
80            coefilu(i, j) = 0.      coefilv = 0.
81            coefilv(i, j) = 0.      coefilu2 = 0.
82            coefilu2(i, j) = 0.      coefilv2 = 0.
           coefilv2(i, j) = 0.  
        end DO  
     END DO  
83    
84      ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv      ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv
85    
     modemax = iim  
     imx = iim  
   
     PRINT *, 'TRUNCATION AT ', imx  
   
86      DO j = 2, jjm / 2 + 1      DO j = 2, jjm / 2 + 1
87         IF (cos(rlatu(j)) / colat0 < 1. &         IF (cos(rlatu(j)) / colat0 < 1. &
88              .and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j              .and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j
89    
90         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &
91              .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) &              .and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &
92              jfiltsu = jjm - j + 2              jfiltsu = jjm - j + 2
93      END DO      END DO
94    
95      DO j = 1, jjm/2      DO j = 1, jjm / 2
96         cof = cos(rlatv(j))/colat0         IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &
97         IF (cof < 1.) THEN              jfiltnv = j
98            IF (rlamda(imx)*cos(rlatv(j)) < 1.) jfiltnv = j  
99         END IF         IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &
100                * 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  
101      END DO      END DO
102    
103      IF (jfiltnu <= 0) jfiltnu = 1      IF (jfiltnu <= 0) jfiltnu = 1
104      IF (jfiltnu > jjm/2+1) THEN      IF (jfiltnu > jjm / 2 + 1) THEN
105         PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu         PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu
106         STOP 1         STOP 1
107      END IF      END IF
# Line 157  contains Line 113  contains
113      END IF      END IF
114    
115      IF (jfiltnv <= 0) jfiltnv = 1      IF (jfiltnv <= 0) jfiltnv = 1
116      IF (jfiltnv > jjm/2) THEN      IF (jfiltnv > jjm / 2) THEN
117         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv
118         STOP 1         STOP 1
119      END IF      END IF
# Line 179  contains Line 135  contains
135      END DO      END DO
136    
137      DO j = 2, jfiltnu      DO j = 2, jfiltnu
138         DO k = 2, modemax         DO k = 2, iim
139            cof = rlamda(k) * cos(rlatu(j))            IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
           IF (cof < 1.) exit  
140         end DO         end DO
141         if (k == modemax + 1) cycle         if (k == iim + 1) cycle
142         modfrstu(j) = k         modfrstu(j) = k
143    
144         kf = modfrstu(j)         kf = modfrstu(j)
145         DO k = kf, modemax         DO k = kf, iim
146            cof = rlamda(k)*cos(rlatu(j))            cof = rlamda(k) * cos(rlatu(j))
147            coefilu(k, j) = cof - 1.            coefilu(k, j) = cof - 1.
148            coefilu2(k, j) = cof*cof - 1.            coefilu2(k, j) = cof**2 - 1.
149         end DO         end DO
150      END DO      END DO
151    
152      DO j = 1, jfiltnv      DO j = 1, jfiltnv
153         DO k = 2, modemax         DO k = 2, iim
154            cof = rlamda(k)*cos(rlatv(j))            IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
           IF (cof < 1.) exit  
155         end DO         end DO
156         if (k == modemax + 1) cycle         if (k == iim + 1) cycle
157         modfrstv(j) = k         modfrstv(j) = k
158    
159         kf = modfrstv(j)         kf = modfrstv(j)
160         DO k = kf, modemax         DO k = kf, iim
161            cof = rlamda(k)*cos(rlatv(j))            cof = rlamda(k) * cos(rlatv(j))
162            coefilv(k, j) = cof - 1.            coefilv(k, j) = cof - 1.
163            coefilv2(k, j) = cof*cof - 1.            coefilv2(k, j) = cof**2 - 1.
164         end DO         end DO
165      end DO      end DO
166    
167      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
168         DO k = 2, modemax         DO k = 2, iim
169            cof = rlamda(k)*cos(rlatu(j))            IF (rlamda(k) * cos(rlatu(j)) < 1.) exit
           IF (cof < 1.) exit  
170         end DO         end DO
171         if (k == modemax + 1) cycle         if (k == iim + 1) cycle
172         modfrstu(j) = k         modfrstu(j) = k
173    
174         kf = modfrstu(j)         kf = modfrstu(j)
175         DO k = kf, modemax         DO k = kf, iim
176            cof = rlamda(k)*cos(rlatu(j))            cof = rlamda(k) * cos(rlatu(j))
177            coefilu(k, j) = cof - 1.            coefilu(k, j) = cof - 1.
178            coefilu2(k, j) = cof*cof - 1.            coefilu2(k, j) = cof**2 - 1.
179         end DO         end DO
180      end DO      end DO
181    
182      DO j = jfiltsv, jjm      DO j = jfiltsv, jjm
183         DO k = 2, modemax         DO k = 2, iim
184            cof = rlamda(k)*cos(rlatv(j))            IF (rlamda(k) * cos(rlatv(j)) < 1.) exit
           IF (cof < 1.) exit  
185         end DO         end DO
186         if (k == modemax + 1) cycle         if (k == iim + 1) cycle
187         modfrstv(j) = k         modfrstv(j) = k
188    
189         kf = modfrstv(j)         kf = modfrstv(j)
190         DO k = kf, modemax         DO k = kf, iim
191            cof = rlamda(k)*cos(rlatv(j))            cof = rlamda(k) * cos(rlatv(j))
192            coefilv(k, j) = cof - 1.            coefilv(k, j) = cof - 1.
193            coefilv2(k, j) = cof*cof - 1.            coefilv2(k, j) = cof**2 - 1.
194         end DO         end DO
195      END DO      END DO
196    
197      IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN      IF (jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2) THEN
198         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv
199         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu
200    
# Line 255  contains Line 207  contains
207      PRINT *, 'Modes premiers u '      PRINT *, 'Modes premiers u '
208      PRINT 334, modfrstu      PRINT 334, modfrstu
209    
210      IF (nfilun < jfiltnu) THEN      allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
211         PRINT *, 'le parametre nfilun utilise pour la matrice ', &      allocate(matricevn(iim, iim, jfiltnv))
212              'matriceun est trop petit ! '      allocate(matricevs(iim, iim, jfiltsv:jjm))
213         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  
214    
215      ! Calcul de la matrice filtre 'matriceu' pour les champs situes      ! Calcul de la matrice filtre 'matriceu' pour les champs situes
216      ! sur la grille scalaire      ! sur la grille scalaire
# Line 338  contains Line 222  contains
222            else            else
223               coff = coefilu(i, j)               coff = coefilu(i, j)
224            end IF            end IF
225            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
226         END DO         END DO
227         matriceun(:, :, j) = matmul(eignfnv, eignft)         matriceun(:, :, j) = matmul(eignfnv, eignft)
228      END DO      END DO
# Line 352  contains Line 236  contains
236            end IF            end IF
237            eignft(i, :) = eignfnv(:, i) * coff            eignft(i, :) = eignfnv(:, i) * coff
238         END DO         END DO
239         matriceus(:, :, j - jfiltsu + 1) = matmul(eignfnv, eignft)         matriceus(:, :, j) = matmul(eignfnv, eignft)
240      END DO      END DO
241    
242      ! Calcul de la matrice filtre 'matricev' pour les champs situes      ! Calcul de la matrice filtre 'matricev' pour les champs situes
# Line 365  contains Line 249  contains
249            else            else
250               coff = coefilv(i, j)               coff = coefilv(i, j)
251            end IF            end IF
252            eignft(i, :) = eignfnu(:, i)*coff            eignft(i, :) = eignfnu(:, i) * coff
253         END DO         END DO
254         matricevn(:, :, j) = matmul(eignfnu, eignft)         matricevn(:, :, j) = matmul(eignfnu, eignft)
255      END DO      END DO
# Line 377  contains Line 261  contains
261            else            else
262               coff = coefilv(i, j)               coff = coefilv(i, j)
263            end IF            end IF
264            eignft(i, :) = eignfnu(:, i)*coff            eignft(i, :) = eignfnu(:, i) * coff
265         END DO         END DO
266         matricevs(:, :, j-jfiltsv+1) = matmul(eignfnu, eignft)         matricevs(:, :, j) = matmul(eignfnu, eignft)
267      END DO      END DO
268    
269      ! Calcul de la matrice filtre 'matrinv' pour les champs situes      ! Calcul de la matrice filtre 'matrinv' pour les champs situes
# Line 390  contains Line 274  contains
274            IF (i < modfrstu(j)) then            IF (i < modfrstu(j)) then
275               coff = 0.               coff = 0.
276            else            else
277               coff = coefilu(i, j)/(1.+coefilu(i, j))               coff = coefilu(i, j) / (1. + coefilu(i, j))
278            end IF            end IF
279            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
280         END DO         END DO
281         matrinvn(:, :, j) = matmul(eignfnv, eignft)         matrinvn(:, :, j) = matmul(eignfnv, eignft)
282      END DO      END DO
# Line 402  contains Line 286  contains
286            IF (i < modfrstu(j)) then            IF (i < modfrstu(j)) then
287               coff = 0.               coff = 0.
288            else            else
289               coff = coefilu(i, j)/(1.+coefilu(i, j))               coff = coefilu(i, j) / (1. + coefilu(i, j))
290            end IF            end IF
291            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
292         END DO         END DO
293         matrinvs(:, :, j-jfiltsu+1) = matmul(eignfnv, eignft)         matrinvs(:, :, j) = matmul(eignfnv, eignft)
294      END DO      END DO
295    
296  334 FORMAT (1X, 24I3)  334 FORMAT (1X, 24I3)
Legend:
Removed from v.82  
changed lines
  Added in v.143
  ViewVC Help
Powered by ViewVC 1.1.21