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

Diff of /trunk/filtrez/inifilr.f90

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

revision 140 by guez, Fri Jun 5 18:58:06 2015 UTC revision 154 by guez, Tue Jul 7 17:49:23 2015 UTC
# Line 31  contains Line 31  contains
31    
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, eignfnu, eignfnv      use inifgn_m, only: inifgn
35        use jumble, only: new_unit
36      use nr_util, only: pi      use nr_util, only: pi
37    
38      ! Local:      ! Local:
39      REAL dlatu(jjm)      REAL dlatu(jjm)
40      REAL rlamda(2: iim)      REAL rlamda(2: iim)
41      real eignvl(iim) ! eigenvalues      real eignvl(iim) ! eigenvalues sorted in descending order
42      REAL lamdamax, cof      REAL cof
43      INTEGER i, j, 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):      ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
51      real coefilu(iim, jjm), coefilv(iim, jjm)      real coefilu(iim, jjm), coefilv(iim, jjm)
52      real coefilu2(iim, jjm), coefilv2(iim, jjm)      real coefilu2(iim, jjm), coefilv2(iim, jjm)
53    
54      integer modfrstu(jjm), modfrstv(jjm)      ! Index of the mode from where modes are filtered:
55      ! index of the mode from where modes are filtered      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 68  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)))  
78    
79      DO j = 1, jjm      ! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
        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  
   
     PRINT *, 'TRUNCATION AT ', iim  
   
     DO j = 2, jjm / 2 + 1  
        IF (cos(rlatu(j)) / colat0 < 1. &  
             .and. rlamda(iim) * cos(rlatu(j)) < 1.) jfiltnu = j  
80    
81         IF (cos(rlatu(jjm - j + 2)) / colat0 < 1. &      jfiltnu = (jjm + 1) / 2
82              .and. rlamda(iim) * cos(rlatu(jjm - j + 2)) < 1.) &      do while (cos(rlatu(jfiltnu)) >= colat0 &
83              jfiltsu = jjm - j + 2           .or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
84      END DO         jfiltnu = jfiltnu - 1
85        end do
86      DO j = 1, jjm / 2  
87         IF (cos(rlatv(j)) / colat0 < 1. .and. rlamda(iim) * cos(rlatv(j)) < 1.) &      jfiltsu = jjm / 2 + 2
88              jfiltnv = j      do while (cos(rlatu(jfiltsu)) >= colat0 &
89             .or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
90         IF (cos(rlatv(jjm - j + 1)) / colat0 < 1. .and. rlamda(iim) &         jfiltsu = jfiltsu + 1
91              * cos(rlatv(jjm - j + 1)) < 1.) jfiltsv = jjm - j + 1      end do
92      END DO  
93        jfiltnv = jjm / 2
94      IF (jfiltnu <= 0) jfiltnu = 1      do while ((cos(rlatv(jfiltnv)) >= colat0 &
95      IF (jfiltnu > jjm / 2 + 1) THEN           .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
96         PRINT *, 'jfiltnu en dehors des valeurs acceptables ', jfiltnu         jfiltnv = jfiltnv - 1
97         STOP 1      end do
98      END IF  
99        if (cos(rlatv(jfiltnv)) >= colat0 &
100      IF (jfiltsu <= 0) jfiltsu = 1           .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
101      IF (jfiltsu > jjm + 1) THEN         ! {jfiltnv == 1}
102         PRINT *, 'jfiltsu en dehors des valeurs acceptables ', jfiltsu         PRINT *, 'Could not find jfiltnv.'
103         STOP 1         STOP 1
104      END IF      END IF
105    
106      IF (jfiltnv <= 0) jfiltnv = 1      jfiltsv = (jjm + 1)/ 2 + 1
107      IF (jfiltnv > jjm / 2) THEN      do while ((cos(rlatv(jfiltsv)) >= colat0 &
108         PRINT *, 'jfiltnv en dehors des valeurs acceptables ', jfiltnv           .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      IF (jfiltsv <= 0) jfiltsv = 1      PRINT *, 'jfiltnu =', jfiltnu
120      IF (jfiltsv > jjm) THEN      PRINT *, 'jfiltsu =', jfiltsu
121         PRINT *, 'jfiltsv en dehors des valeurs acceptables ', jfiltsv      PRINT *, 'jfiltnv =', jfiltnv
122         STOP 1      PRINT *, 'jfiltsv =', jfiltsv
123      END IF  
124        ! Determination de coefilu, coefilv, modfrst[ns][uv]:
125      PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &  
126           jfiltsu      allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
127        allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
128      ! Determination de coefilu, coefilv, n=modfrstu, modfrstv      coefilu = 0.
129        coefilv = 0.
130      DO j = 1, jjm      coefilu2 = 0.
131         modfrstu(j) = iim      coefilv2 = 0.
        modfrstv(j) = iim  
     END DO  
132    
133      DO j = 2, jfiltnu      DO j = 2, jfiltnu
134         DO k = 2, iim         modfrstnu(j) = 2
135            IF (rlamda(k) * cos(rlatu(j)) < 1.) exit         do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
136         end DO              .and. modfrstnu(j) <= iim - 1)
137         if (k == iim + 1) cycle            modfrstnu(j) = modfrstnu(j) + 1
138         modfrstu(j) = k         end do
139    
140         kf = modfrstu(j)         if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
141         DO k = kf, iim            DO k = modfrstnu(j), iim
142            cof = rlamda(k) * cos(rlatu(j))               cof = rlamda(k) * cos(rlatu(j))
143            coefilu(k, j) = cof - 1.               coefilu(k, j) = cof - 1.
144            coefilu2(k, j) = cof**2 - 1.               coefilu2(k, j) = cof**2 - 1.
145         end DO            end DO
146           end if
147      END DO      END DO
148    
149      DO j = 1, jfiltnv      DO j = 1, jfiltnv
150         DO k = 2, iim         modfrstnv(j) = 2
151            IF (rlamda(k) * cos(rlatv(j)) < 1.) exit         do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
152         end DO              .and. modfrstnv(j) <= iim - 1)
153         if (k == iim + 1) cycle            modfrstnv(j) = modfrstnv(j) + 1
154         modfrstv(j) = k         end do
155    
156         kf = modfrstv(j)         if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
157         DO k = kf, iim            DO k = modfrstnv(j), iim
158            cof = rlamda(k) * cos(rlatv(j))               cof = rlamda(k) * cos(rlatv(j))
159            coefilv(k, j) = cof - 1.               coefilv(k, j) = cof - 1.
160            coefilv2(k, j) = cof**2 - 1.               coefilv2(k, j) = cof**2 - 1.
161         end DO            end DO
162           end if
163      end DO      end DO
164    
165      DO j = jfiltsu, jjm      DO j = jfiltsu, jjm
166         DO k = 2, iim         modfrstsu(j) = 2
167            IF (rlamda(k) * cos(rlatu(j)) < 1.) exit         do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
168         end DO              .and. modfrstsu(j) <= iim - 1)
169         if (k == iim + 1) cycle            modfrstsu(j) = modfrstsu(j) + 1
170         modfrstu(j) = k         end do
171    
172         kf = modfrstu(j)         if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
173         DO k = kf, iim            DO k = modfrstsu(j), iim
174            cof = rlamda(k) * cos(rlatu(j))               cof = rlamda(k) * cos(rlatu(j))
175            coefilu(k, j) = cof - 1.               coefilu(k, j) = cof - 1.
176            coefilu2(k, j) = cof**2 - 1.               coefilu2(k, j) = cof**2 - 1.
177         end DO            end DO
178           end if
179      end DO      end DO
180    
181      DO j = jfiltsv, jjm      DO j = jfiltsv, jjm
182         DO k = 2, iim         modfrstsv(j) = 2
183            IF (rlamda(k) * cos(rlatv(j)) < 1.) exit         do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
184         end DO              .and. modfrstsv(j) <= iim - 1)
185         if (k == iim + 1) cycle            modfrstsv(j) = modfrstsv(j) + 1
186         modfrstv(j) = k         end do
187    
188         kf = modfrstv(j)         if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
189         DO k = kf, iim            DO k = modfrstsv(j), iim
190            cof = rlamda(k) * cos(rlatv(j))               cof = rlamda(k) * cos(rlatv(j))
191            coefilv(k, j) = cof - 1.               coefilv(k, j) = cof - 1.
192            coefilv2(k, j) = cof**2 - 1.               coefilv2(k, j) = cof**2 - 1.
193         end DO            end DO
194      END DO         end if
195        END DO
196      IF (jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2) THEN  
197         IF (jfiltnv == jfiltsv) jfiltsv = 1 + jfiltnv      call new_unit(unit)
198         IF (jfiltnu == jfiltsu) jfiltsu = 1 + jfiltnu      open(unit, file = "inifilr_out.txt", status = "replace", action = "write")
199        write(unit, fmt = *) '"EIGNVL"', eignvl
200         PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &      write(unit, fmt = *) '"modfrstnu"', modfrstnu
201              jfiltsu      write(unit, fmt = *) '"modfrstsu"', modfrstsu
202      END IF      write(unit, fmt = *) '"modfrstnv"', modfrstnv
203        write(unit, fmt = *) '"modfrstsv"', modfrstsv
204      PRINT *, 'Modes premiers v '      close(unit)
     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 232  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)
# Line 244  contains Line 225  contains
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 259  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)
# Line 271  contains Line 252  contains
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)
# Line 286  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))
# Line 298  contains Line 279  contains
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))
# Line 308  contains Line 289  contains
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.140  
changed lines
  Added in v.154
  ViewVC Help
Powered by ViewVC 1.1.21