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