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

Diff of /trunk/filtrez/filtreg_scal.f

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

trunk/libf/filtrez/filtreg.f90 revision 25 by guez, Fri Mar 5 16:43:45 2010 UTC trunk/filtrez/filtreg.f revision 107 by guez, Thu Sep 11 15:09:15 2014 UTC
# Line 1  Line 1 
1  SUBROUTINE filtreg(champ, nlat, nbniv, ifiltre, iaire, griscal, iter)  module filtreg_m
   
   ! From filtrez/filtreg.F, version 1.1.1.1 2004/05/19 12:53:09  
   
   ! Auteur: P. Le Van 07/10/97  
   ! Objet: filtre matriciel longitudinal, avec les matrices précalculées  
   ! pour l'operateur filtre.  
   
   USE dimens_m  
   USE paramet_m  
   USE parafilt  
   USE coefils  
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    INTEGER, intent(in):: nlat ! nombre de latitudes a filtrer  contains
6    integer, intent(in):: nbniv ! nombre de niveaux verticaux a filtrer  
7      SUBROUTINE filtreg(champ, direct, intensive)
8    
9    REAL, intent(inout):: champ(iip1, nlat, nbniv)      ! From filtrez/filtreg.F, version 1.1.1.1, 2004/05/19 12:53:09
10    ! en entree : champ a filtrer, en sortie : champ filtre      ! Author: P. Le Van
11        ! Objet : filtre matriciel longitudinal, avec les matrices précalculées
12        ! pour l'opérateur filtre.
13    
14        USE coefils, ONLY: sddu, sddv, unsddu, unsddv
15        USE dimens_m, ONLY: iim, jjm
16        use inifilr_m, only: jfiltnu, jfiltnv, jfiltsu, jfiltsv, matriceun, &
17             matriceus, matricevn, matricevs, matrinvn, matrinvs
18        use nr_util, only: assert
19    
20        REAL, intent(inout):: champ(:, :, :) ! (iim + 1, nlat, nbniv)
21        ! en entrée : champ à filtrer, en sortie : champ filtré
22    
23        logical, intent(in):: direct ! filtre direct ou inverse
24    
25        logical, intent(in):: intensive
26        ! champ intensif ou extensif (pondéré par les aires)
27    
28        ! Local:
29        LOGICAL griscal
30        INTEGER nlat ! nombre de latitudes à filtrer
31        integer nbniv ! nombre de niveaux verticaux à filtrer
32        INTEGER jdfil1, jdfil2, jffil1, jffil2, jdfil, jffil
33        INTEGER i, j, l, k
34        REAL eignq(iim), sdd1(iim), sdd2(iim)
35        INTEGER hemisph
36    
37        !-----------------------------------------------------------
38    
39        call assert(size(champ, 1) == iim + 1, "filtreg iim + 1")
40        nlat = size(champ, 2)
41        nbniv = size(champ, 3)
42        call assert(nlat == jjm .or. nlat == jjm + 1, "filtreg nlat")
43        griscal = nlat == jjm + 1
44    
45        IF (.not. direct .AND. nlat == jjm) THEN
46           PRINT *, 'filtreg: inverse filter on scalar grid only'
47           STOP 1
48        END IF
49    
50        IF (griscal) THEN
51           IF (intensive) THEN
52              sdd1 = sddv
53              sdd2 = unsddv
54           ELSE
55              sdd1 = unsddv
56              sdd2 = sddv
57           END IF
58    
59           jdfil1 = 2
60           jffil1 = jfiltnu
61           jdfil2 = jfiltsu
62           jffil2 = jjm
63        ELSE
64           IF (intensive) THEN
65              sdd1 = sddu
66              sdd2 = unsddu
67           ELSE
68              sdd1 = unsddu
69              sdd2 = sddu
70           END IF
71    
72           jdfil1 = 1
73           jffil1 = jfiltnv
74           jdfil2 = jfiltsv
75           jffil2 = jjm
76        END IF
77    
78        DO hemisph = 1, 2
79           IF (hemisph==1) THEN
80              jdfil = jdfil1
81              jffil = jffil1
82           ELSE
83              jdfil = jdfil2
84              jffil = jffil2
85           END IF
86    
87           loop_vertical: DO l = 1, nbniv
88              loop_latitude: DO j = jdfil, jffil
89                 DO i = 1, iim
90                    champ(i, j, l) = champ(i, j, l)*sdd1(i)
91                 END DO
92    
93                 IF (hemisph==1) THEN
94                    IF (.not. direct) THEN
95                       DO k = 1, iim
96                          eignq(k) = 0.
97                       END DO
98                       DO k = 1, iim
99                          DO i = 1, iim
100                             eignq(k) = eignq(k) + matrinvn(k, i, j)*champ(i, j, l)
101                          END DO
102                       END DO
103                    ELSE IF (griscal) THEN
104                       DO k = 1, iim
105                          eignq(k) = 0.
106                       END DO
107                       DO i = 1, iim
108                          DO k = 1, iim
109                             eignq(k) = eignq(k) + matriceun(k, i, j) &
110                                  * champ(i, j, l)
111                          END DO
112                       END DO
113                    ELSE
114                       DO k = 1, iim
115                          eignq(k) = 0.
116                       END DO
117                       DO i = 1, iim
118                          DO k = 1, iim
119                             eignq(k) = eignq(k) + matricevn(k, i, j) &
120                                  * champ(i, j, l)
121                          END DO
122                       END DO
123                    END IF
124                 ELSE
125                    IF (.not. direct) THEN
126                       DO k = 1, iim
127                          eignq(k) = 0.
128                       END DO
129                       DO i = 1, iim
130                          DO k = 1, iim
131                             eignq(k) = eignq(k) + matrinvs(k, i, j-jfiltsu+1) &
132                                  *champ(i, j, l)
133                          END DO
134                       END DO
135                    ELSE IF (griscal) THEN
136                       DO k = 1, iim
137                          eignq(k) = 0.
138                       END DO
139                       DO i = 1, iim
140                          DO k = 1, iim
141                             eignq(k) = eignq(k) + matriceus(k, i, j-jfiltsu+1) &
142                                  *champ(i, j , l)
143                          END DO
144                       END DO
145                    ELSE
146                       DO k = 1, iim
147                          eignq(k) = 0.
148                       END DO
149                       DO i = 1, iim
150                          DO k = 1, iim
151                             eignq(k) = eignq(k) + matricevs(k, i, j-jfiltsv+1) &
152                                  *champ(i, j , l)
153                          END DO
154                       END DO
155                    END IF
156                 END IF
157    
158                 IF (direct) THEN
159                    DO i = 1, iim
160                       champ(i, j, l) = (champ(i, j, l)+eignq(i))*sdd2(i)
161                    end DO
162                 ELSE
163                    DO i = 1, iim
164                       champ(i, j, l) = (champ(i, j, l)-eignq(i))*sdd2(i)
165                    end DO
166                 END IF
167    
168                 champ(iim + 1, j, l) = champ(1, j, l)
169              END DO loop_latitude
170           END DO loop_vertical
171        end DO
172    
173    integer, intent(in):: ifiltre    END SUBROUTINE filtreg
   !  +1 Transformee directe  
   ! -1 Transformee inverse  
   ! +2 Filtre directe  
   ! -2 Filtre inverse  
   ! Variable Intensive  
   ! ifiltre = 1 filtre directe  
   ! ifiltre =-1 filtre inverse  
   ! Variable Extensive  
   ! ifiltre = 2 filtre directe  
   ! ifiltre =-2 filtre inverse  
   
   integer, intent(in):: iaire  
   !  1 si champ intensif  
   ! 2 si champ extensif (pondere par les aires)  
   
   integer, intent(in):: iter  
   !  1 filtre simple  
   
   LOGICAL, intent(in):: griscal  
   
   ! Variables local to the procedure:  
   
   INTEGER jdfil1, jdfil2, jffil1, jffil2, jdfil, jffil  
   INTEGER i, j, l, k  
   REAL matriceun, matriceus, matricevn, matricevs, matrinvn, matrinvs  
   COMMON /matrfil/matriceun(iim, iim, nfilun), matriceus(iim, iim, nfilus), &  
        matricevn(iim, iim, nfilvn), matricevs(iim, iim, nfilvs), &  
        matrinvn(iim, iim, nfilun), matrinvs(iim, iim, nfilus)  
   REAL eignq(iim), sdd1(iim), sdd2(iim)  
   INTEGER hemisph  
   
   !-----------------------------------------------------------  
   
   IF (ifiltre==1 .OR. ifiltre==-1) STOP &  
        'Pas de transformee simple dans cette version'  
   
   IF (iter==2) THEN  
      PRINT *, ' Pas d iteration du filtre dans cette version !', &  
           ' Utiliser old_filtreg et repasser !'  
      STOP  
   END IF  
   
   IF (ifiltre==-2 .AND. .NOT. griscal) THEN  
      PRINT *, ' Cette routine ne calcule le filtre inverse que ', &  
           ' sur la grille des scalaires !'  
      STOP  
   END IF  
   
   IF (ifiltre/=2 .AND. ifiltre/=-2) THEN  
      PRINT *, ' Probleme dans filtreg car ifiltre NE 2 et NE -2', &  
           ' corriger et repasser !'  
      STOP  
   END IF  
   
   IF (griscal) THEN  
      IF (nlat/=jjp1) THEN  
         PRINT 1111  
         STOP  
      ELSE  
   
         IF (iaire==1) THEN  
            CALL scopy(iim, sddv, 1, sdd1, 1)  
            CALL scopy(iim, unsddv, 1, sdd2, 1)  
         ELSE  
            CALL scopy(iim, unsddv, 1, sdd1, 1)  
            CALL scopy(iim, sddv, 1, sdd2, 1)  
         END IF  
   
         jdfil1 = 2  
         jffil1 = jfiltnu  
         jdfil2 = jfiltsu  
         jffil2 = jjm  
      END IF  
   ELSE  
      IF (nlat/=jjm) THEN  
         PRINT 2222  
         STOP  
      ELSE  
   
         IF (iaire==1) THEN  
            CALL scopy(iim, sddu, 1, sdd1, 1)  
            CALL scopy(iim, unsddu, 1, sdd2, 1)  
         ELSE  
            CALL scopy(iim, unsddu, 1, sdd1, 1)  
            CALL scopy(iim, sddu, 1, sdd2, 1)  
         END IF  
   
         jdfil1 = 1  
         jffil1 = jfiltnv  
         jdfil2 = jfiltsv  
         jffil2 = jjm  
      END IF  
   END IF  
   
   
   DO hemisph = 1, 2  
   
      IF (hemisph==1) THEN  
         jdfil = jdfil1  
         jffil = jffil1  
      ELSE  
         jdfil = jdfil2  
         jffil = jffil2  
      END IF  
   
   
      DO l = 1, nbniv  
         DO j = jdfil, jffil  
   
   
            DO i = 1, iim  
               champ(i, j, l) = champ(i, j, l)*sdd1(i)  
            END DO  
   
   
            IF (hemisph==1) THEN  
   
               IF (ifiltre==-2) THEN  
                  DO k = 1, iim  
                     eignq(k) = 0.0  
                  END DO  
                  DO k = 1, iim  
                     DO i = 1, iim  
                        eignq(k) = eignq(k) + matrinvn(k, i, j)*champ(i, j, l)  
                     END DO  
                  END DO  
               ELSE IF (griscal) THEN  
                  DO k = 1, iim  
                     eignq(k) = 0.0  
                  END DO  
                  DO i = 1, iim  
                     DO k = 1, iim  
                        eignq(k) = eignq(k) + matriceun(k, i, j)*champ(i, j, l)  
                     END DO  
                  END DO  
               ELSE  
                  DO k = 1, iim  
                     eignq(k) = 0.0  
                  END DO  
                  DO i = 1, iim  
                     DO k = 1, iim  
                        eignq(k) = eignq(k) + matricevn(k, i, j)*champ(i, j, l)  
                     END DO  
                  END DO  
               END IF  
   
            ELSE  
   
               IF (ifiltre==-2) THEN  
                  DO k = 1, iim  
                     eignq(k) = 0.0  
                  END DO  
                  DO i = 1, iim  
                     DO k = 1, iim  
                        eignq(k) = eignq(k) + matrinvs(k, i, j-jfiltsu+1)*champ(i, j, &  
                             l)  
                     END DO  
                  END DO  
               ELSE IF (griscal) THEN  
                  DO k = 1, iim  
                     eignq(k) = 0.0  
                  END DO  
                  DO i = 1, iim  
                     DO k = 1, iim  
                        eignq(k) = eignq(k) + matriceus(k, i, j-jfiltsu+1)*champ(i, j &  
                             , l)  
                     END DO  
                  END DO  
               ELSE  
                  DO k = 1, iim  
                     eignq(k) = 0.0  
                  END DO  
                  DO i = 1, iim  
                     DO k = 1, iim  
                        eignq(k) = eignq(k) + matricevs(k, i, j-jfiltsv+1)*champ(i, j &  
                             , l)  
                     END DO  
                  END DO  
               END IF  
   
            END IF  
   
            IF (ifiltre==2) THEN  
               DO i = 1, iim  
                  champ(i, j, l) = (champ(i, j, l)+eignq(i))*sdd2(i)  
               end DO  
            ELSE  
               DO i = 1, iim  
                  champ(i, j, l) = (champ(i, j, l)-eignq(i))*sdd2(i)  
               end DO  
            END IF  
   
            champ(iip1, j, l) = champ(1, j, l)  
   
         END DO  
   
      END DO  
   
   end DO  
   
 1111 FORMAT (//20X, 'ERREUR dans le dimensionnement du tableau &  
        & CHAMP a filtrer, sur la grille des scalaires'/)  
 2222 FORMAT (//20X, 'ERREUR dans le dimensionnement du tableau &  
        &CHAMP a filtrer, sur la grille de V ou de Z'/)  
174    
175  END SUBROUTINE filtreg  end module filtreg_m
Legend:
Removed from v.25  
changed lines
  Added in v.107
  ViewVC Help
Powered by ViewVC 1.1.21