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

Diff of /trunk/filtrez/inifilr.f

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

revision 136 by guez, Thu Apr 30 18:35:49 2015 UTC revision 151 by guez, Tue Jun 23 15:14:20 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.
31      ! 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.  
   
     USE coefils, ONLY : coefilu, coefilu2, coefilv, coefilv2, eignfnu, &  
          eignfnv, modfrstu, modfrstv  
     USE comgeom, ONLY : rlatu, rlatv, xprimu  
32      USE dimens_m, ONLY : iim, jjm      USE dimens_m, ONLY : iim, jjm
33      use inifgn_m, only: inifgn      USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
34        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 : grossismx  
37    
38      ! Local:      ! Local:
39      REAL dlatu(jjm)      REAL dlatu(jjm)
40      REAL rlamda(2: iim), eignvl(iim)      REAL rlamda(2: iim)
41        real eignvl(iim) ! eigenvalues sorted in descending order
42      REAL lamdamax, cof      REAL cof
43      INTEGER i, j, modemax, imx, k, kf      INTEGER i, j, k, unit
44      REAL dymin, colat0      REAL colat0 ! > 0
45      REAL eignft(iim, iim), coff      REAL eignft(iim, iim), coff
46    
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        ! Index of the mode from where modes are filtered:
52        integer, allocatable:: modfrstnu(:), modfrstsu(:)
53        integer, allocatable:: modfrstnv(:), modfrstsv(:)
54    
55      !-----------------------------------------------------------      !-----------------------------------------------------------
56    
57      print *, "Call sequence information: inifilr"      print *, "Call sequence information: inifilr"
58    
59      CALL inifgn(eignvl)      CALL inifgn(eignvl)
60    
     PRINT *, 'EIGNVL '  
     PRINT "(1X, 5E13.6)", eignvl  
   
61      ! compute eigenvalues and eigenfunctions      ! compute eigenvalues and eigenfunctions
62      ! compute the filtering coefficients for scalar lines and      ! compute the filtering coefficients for scalar lines and
63      ! meridional wind v-lines      ! meridional wind v-lines
# Line 72  contains Line 67  contains
67      ! is set equal to zero for the regular grid case      ! is set equal to zero for the regular grid case
68    
69      ! Calcul de colat0      ! Calcul de colat0
70        forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
71      DO j = 1, jjm      colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
        dlatu(j) = rlatu(j) - rlatu(j+1)  
     END DO  
   
     dymin = dlatu(1)  
     DO j = 2, jjm  
        dymin = min(dymin, dlatu(j))  
     END DO  
   
     colat0 = min(0.5, dymin / minval(xprimu(:iim)))  
   
72      PRINT *, 'colat0 = ', colat0      PRINT *, 'colat0 = ', colat0
73    
74      lamdamax = iim / (pi * colat0 / grossismx)      rlamda = iim / (pi * colat0 / grossismx) / sqrt(abs(eignvl(2: iim)))
     rlamda = lamdamax / sqrt(abs(eignvl(2: iim)))  
   
     DO j = 1, jjm  
        DO i = 1, iim  
           coefilu(i, j) = 0.  
           coefilv(i, j) = 0.  
           coefilu2(i, j) = 0.  
           coefilv2(i, j) = 0.  
        end DO  
     END DO  
   
     ! Determination de jfiltnu, jfiltnv, jfiltsu, jfiltsv  
   
     modemax = iim  
     imx = iim  
   
     PRINT *, 'TRUNCATION AT ', imx  
   
     DO j = 2, jjm / 2 + 1  
        IF (cos(rlatu(j)) / colat0 < 1. &  
             .and. rlamda(imx) * cos(rlatu(j)) < 1.) jfiltnu = j  
   
        IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &  
             .and. rlamda(imx) * cos(rlatu(jjm - j + 2)) < 1.) &  
             jfiltsu = jjm - j + 2  
     END DO  
   
     DO j = 1, jjm/2  
        cof = cos(rlatv(j))/colat0  
        IF (cof < 1.) THEN  
           IF (rlamda(imx)*cos(rlatv(j)) < 1.) jfiltnv = j  
        END IF  
   
        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  
     END DO  
   
     IF (jfiltnu <= 0) jfiltnu = 1  
     IF (jfiltnu > jjm/2+1) THEN  
        PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu  
        STOP 1  
     END IF  
75    
76      IF (jfiltsu <= 0) jfiltsu = 1      ! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
     IF (jfiltsu > jjm + 1) THEN  
        PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu  
        STOP 1  
     END IF  
77    
78      IF (jfiltnv <= 0) jfiltnv = 1      jfiltnu = (jjm + 1) / 2
79      IF (jfiltnv > jjm/2) THEN      do while (cos(rlatu(jfiltnu)) >= colat0 &
80         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv           .or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
81           jfiltnu = jfiltnu - 1
82        end do
83    
84        jfiltsu = jjm / 2 + 2
85        do while (cos(rlatu(jfiltsu)) >= colat0 &
86             .or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
87           jfiltsu = jfiltsu + 1
88        end do
89    
90        jfiltnv = jjm / 2
91        do while ((cos(rlatv(jfiltnv)) >= colat0 &
92             .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
93           jfiltnv = jfiltnv - 1
94        end do
95    
96        if (cos(rlatv(jfiltnv)) >= colat0 &
97             .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
98           ! {jfiltnv == 1}
99           PRINT *, 'Could not find jfiltnv.'
100         STOP 1         STOP 1
101      END IF      END IF
102    
103      IF (jfiltsv <= 0) jfiltsv = 1      jfiltsv = (jjm + 1)/ 2 + 1
104      IF (jfiltsv > jjm) THEN      do while ((cos(rlatv(jfiltsv)) >= colat0 &
105         PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv           .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
106           jfiltsv = jfiltsv + 1
107        end do
108    
109        IF (cos(rlatv(jfiltsv)) >= colat0 &
110             .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
111           ! {jfiltsv == jjm}
112           PRINT *, 'Could not find jfiltsv.'
113         STOP 1         STOP 1
114      END IF      END IF
115    
116      PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &      PRINT *, 'jfiltnu =', jfiltnu
117           jfiltsu      PRINT *, 'jfiltsu =', jfiltsu
118        PRINT *, 'jfiltnv =', jfiltnv
119      ! Determination de coefilu, coefilv, n=modfrstu, modfrstv      PRINT *, 'jfiltsv =', jfiltsv
120    
121      DO j = 1, jjm      ! Determination de coefilu, coefilv, modfrst[ns][uv]:
122         modfrstu(j) = iim  
123         modfrstv(j) = iim      allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
124      END DO      allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
125        coefilu = 0.
126        coefilv = 0.
127        coefilu2 = 0.
128        coefilv2 = 0.
129    
130      DO j = 2, jfiltnu      DO j = 2, jfiltnu
131         DO k = 2, modemax         modfrstnu(j) = 2
132            cof = rlamda(k) * cos(rlatu(j))         do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
133            IF (cof < 1.) exit              .and. modfrstnu(j) <= iim - 1)
134         end DO            modfrstnu(j) = modfrstnu(j) + 1
135         if (k == modemax + 1) cycle         end do
136         modfrstu(j) = k  
137           if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
138         kf = modfrstu(j)            DO k = modfrstnu(j), iim
139         DO k = kf, modemax               cof = rlamda(k) * cos(rlatu(j))
140            cof = rlamda(k)*cos(rlatu(j))               coefilu(k, j) = cof - 1.
141            coefilu(k, j) = cof - 1.               coefilu2(k, j) = cof**2 - 1.
142            coefilu2(k, j) = cof*cof - 1.            end DO
143         end DO         end if
144      END DO      END DO
145    
146      DO j = 1, jfiltnv      DO j = 1, jfiltnv
147         DO k = 2, modemax         modfrstnv(j) = 2
148            cof = rlamda(k)*cos(rlatv(j))         do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
149            IF (cof < 1.) exit              .and. modfrstnv(j) <= iim - 1)
150         end DO            modfrstnv(j) = modfrstnv(j) + 1
151         if (k == modemax + 1) cycle         end do
152         modfrstv(j) = k  
153           if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
154         kf = modfrstv(j)            DO k = modfrstnv(j), iim
155         DO k = kf, modemax               cof = rlamda(k) * cos(rlatv(j))
156            cof = rlamda(k)*cos(rlatv(j))               coefilv(k, j) = cof - 1.
157            coefilv(k, j) = cof - 1.               coefilv2(k, j) = cof**2 - 1.
158            coefilv2(k, j) = cof*cof - 1.            end DO
159         end DO         end if
160      end DO      end DO
161    
162      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
163         DO k = 2, modemax         modfrstsu(j) = 2
164            cof = rlamda(k)*cos(rlatu(j))         do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
165            IF (cof < 1.) exit              .and. modfrstsu(j) <= iim - 1)
166         end DO            modfrstsu(j) = modfrstsu(j) + 1
167         if (k == modemax + 1) cycle         end do
168         modfrstu(j) = k  
169           if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
170         kf = modfrstu(j)            DO k = modfrstsu(j), iim
171         DO k = kf, modemax               cof = rlamda(k) * cos(rlatu(j))
172            cof = rlamda(k)*cos(rlatu(j))               coefilu(k, j) = cof - 1.
173            coefilu(k, j) = cof - 1.               coefilu2(k, j) = cof**2 - 1.
174            coefilu2(k, j) = cof*cof - 1.            end DO
175         end DO         end if
176      end DO      end DO
177    
178      DO j = jfiltsv, jjm      DO j = jfiltsv, jjm
179         DO k = 2, modemax         modfrstsv(j) = 2
180            cof = rlamda(k)*cos(rlatv(j))         do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
181            IF (cof < 1.) exit              .and. modfrstsv(j) <= iim - 1)
182         end DO            modfrstsv(j) = modfrstsv(j) + 1
183         if (k == modemax + 1) cycle         end do
184         modfrstv(j) = k  
185           if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
186         kf = modfrstv(j)            DO k = modfrstsv(j), iim
187         DO k = kf, modemax               cof = rlamda(k) * cos(rlatv(j))
188            cof = rlamda(k)*cos(rlatv(j))               coefilv(k, j) = cof - 1.
189            coefilv(k, j) = cof - 1.               coefilv2(k, j) = cof**2 - 1.
190            coefilv2(k, j) = cof*cof - 1.            end DO
191         end DO         end if
192      END DO      END DO
193    
194      IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN      call new_unit(unit)
195         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv      open(unit, file = "inifilr_out.txt", status = "replace", action = "write")
196         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu      write(unit, fmt = *) '"EIGNVL"', eignvl
197        write(unit, fmt = *) '"modfrstnu"', modfrstnu
198         PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &      write(unit, fmt = *) '"modfrstsu"', modfrstsu
199              jfiltsu      write(unit, fmt = *) '"modfrstnv"', modfrstnv
200      END IF      write(unit, fmt = *) '"modfrstsv"', modfrstsv
201        close(unit)
     PRINT *, 'Modes premiers v '  
     PRINT 334, modfrstv  
     PRINT *, 'Modes premiers u '  
     PRINT 334, modfrstu  
202    
203      allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))      allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
204      allocate(matricevn(iim, iim, jfiltnv))      allocate(matricevn(iim, iim, jfiltnv))
# Line 247  contains Line 210  contains
210    
211      DO j = 2, jfiltnu      DO j = 2, jfiltnu
212         DO i = 1, iim         DO i = 1, iim
213            IF (i < modfrstu(j)) then            IF (i < modfrstnu(j)) then
214               coff = 0.               coff = 0.
215            else            else
216               coff = coefilu(i, j)               coff = coefilu(i, j)
217            end IF            end IF
218            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
219         END DO         END DO
220         matriceun(:, :, j) = matmul(eignfnv, eignft)         matriceun(:, :, j) = matmul(eignfnv, eignft)
221      END DO      END DO
222    
223      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
224         DO i = 1, iim         DO i = 1, iim
225            IF (i < modfrstu(j)) then            IF (i < modfrstsu(j)) then
226               coff = 0.               coff = 0.
227            else            else
228               coff = coefilu(i, j)               coff = coefilu(i, j)
# Line 274  contains Line 237  contains
237    
238      DO j = 1, jfiltnv      DO j = 1, jfiltnv
239         DO i = 1, iim         DO i = 1, iim
240            IF (i < modfrstv(j)) then            IF (i < modfrstnv(j)) then
241               coff = 0.               coff = 0.
242            else            else
243               coff = coefilv(i, j)               coff = coefilv(i, j)
244            end IF            end IF
245            eignft(i, :) = eignfnu(:, i)*coff            eignft(i, :) = eignfnu(:, i) * coff
246         END DO         END DO
247         matricevn(:, :, j) = matmul(eignfnu, eignft)         matricevn(:, :, j) = matmul(eignfnu, eignft)
248      END DO      END DO
249    
250      DO j = jfiltsv, jjm      DO j = jfiltsv, jjm
251         DO i = 1, iim         DO i = 1, iim
252            IF (i < modfrstv(j)) then            IF (i < modfrstsv(j)) then
253               coff = 0.               coff = 0.
254            else            else
255               coff = coefilv(i, j)               coff = coefilv(i, j)
256            end IF            end IF
257            eignft(i, :) = eignfnu(:, i)*coff            eignft(i, :) = eignfnu(:, i) * coff
258         END DO         END DO
259         matricevs(:, :, j) = matmul(eignfnu, eignft)         matricevs(:, :, j) = matmul(eignfnu, eignft)
260      END DO      END DO
# Line 301  contains Line 264  contains
264    
265      DO j = 2, jfiltnu      DO j = 2, jfiltnu
266         DO i = 1, iim         DO i = 1, iim
267            IF (i < modfrstu(j)) then            IF (i < modfrstnu(j)) then
268               coff = 0.               coff = 0.
269            else            else
270               coff = coefilu(i, j)/(1.+coefilu(i, j))               coff = coefilu(i, j) / (1. + coefilu(i, j))
271            end IF            end IF
272            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
273         END DO         END DO
274         matrinvn(:, :, j) = matmul(eignfnv, eignft)         matrinvn(:, :, j) = matmul(eignfnv, eignft)
275      END DO      END DO
276    
277      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
278         DO i = 1, iim         DO i = 1, iim
279            IF (i < modfrstu(j)) then            IF (i < modfrstsu(j)) then
280               coff = 0.               coff = 0.
281            else            else
282               coff = coefilu(i, j)/(1.+coefilu(i, j))               coff = coefilu(i, j) / (1. + coefilu(i, j))
283            end IF            end IF
284            eignft(i, :) = eignfnv(:, i)*coff            eignft(i, :) = eignfnv(:, i) * coff
285         END DO         END DO
286         matrinvs(:, :, j) = matmul(eignfnv, eignft)         matrinvs(:, :, j) = matmul(eignfnv, eignft)
287      END DO      END DO
288    
 334 FORMAT (1X, 24I3)  
   
289    END SUBROUTINE inifilr    END SUBROUTINE inifilr
290    
291  end module inifilr_m  end module inifilr_m
Legend:
Removed from v.136  
changed lines
  Added in v.151
  ViewVC Help
Powered by ViewVC 1.1.21