trunk/filtrez/filtreg_scal.f

module filtreg_m

  IMPLICIT NONE

contains

  SUBROUTINE filtreg(champ, direct, intensive)

    ! From filtrez/filtreg.F, version 1.1.1.1, 2004/05/19 12:53:09
    ! Author: P. Le Van
    ! Objet : filtre matriciel longitudinal, avec les matrices pr\'ecalcul\'ees
    ! pour l'op\'erateur filtre. 

    USE coefils, ONLY: sddu, sddv, unsddu, unsddv
    USE dimens_m, ONLY: iim, jjm
    use inifilr_m, only: jfiltnu, jfiltnv, jfiltsu, jfiltsv, matriceun, &
         matriceus, matricevn, matricevs, matrinvn, matrinvs
    use nr_util, only: assert

    REAL, intent(inout):: champ(:, :, :) ! (iim + 1, nlat, nbniv)
    ! en entr\'ee : champ \`a filtrer, en sortie : champ filtr\'e

    logical, intent(in):: direct ! filtre direct ou inverse 

    logical, intent(in):: intensive
    ! champ intensif ou extensif (pond\'er\'e par les aires)

    ! Local:
    LOGICAL griscal
    INTEGER nlat ! nombre de latitudes \`a filtrer
    integer nbniv ! nombre de niveaux verticaux \`a filtrer 
    INTEGER jdfil1, jdfil2, jffil1, jffil2, jdfil, jffil
    INTEGER i, j, l, k
    REAL eignq(iim), sdd1(iim), sdd2(iim)
    INTEGER hemisph

    !-----------------------------------------------------------

    call assert(size(champ, 1) == iim + 1, "filtreg iim + 1")
    nlat = size(champ, 2)
    nbniv = size(champ, 3)
    call assert(nlat == jjm .or. nlat == jjm + 1, "filtreg nlat")
    griscal = nlat == jjm + 1

    IF (.not. direct .AND. nlat == jjm) THEN
       PRINT *, 'filtreg: inverse filter on scalar grid only'
       STOP 1
    END IF

    IF (griscal) THEN
       IF (intensive) THEN
          sdd1 = sddv
          sdd2 = unsddv
       ELSE
          sdd1 = unsddv
          sdd2 = sddv
       END IF

       jdfil1 = 2
       jffil1 = jfiltnu
       jdfil2 = jfiltsu
       jffil2 = jjm
    ELSE
       IF (intensive) THEN
          sdd1 = sddu
          sdd2 = unsddu
       ELSE
          sdd1 = unsddu
          sdd2 = sddu
       END IF

       jdfil1 = 1
       jffil1 = jfiltnv
       jdfil2 = jfiltsv
       jffil2 = jjm
    END IF

    loop_hemisph: DO hemisph = 1, 2
       IF (hemisph==1) THEN
          jdfil = jdfil1
          jffil = jffil1
       ELSE
          jdfil = jdfil2
          jffil = jffil2
       END IF

       loop_vertical: DO l = 1, nbniv
          loop_latitude: 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 (.not. direct) THEN
                   DO k = 1, iim
                      eignq(k) = 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.
                   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.
                   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 (.not. direct) THEN
                   DO k = 1, iim
                      eignq(k) = 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.
                   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.
                   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 (direct) 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(iim + 1, j, l) = champ(1, j, l)
          END DO loop_latitude
       END DO loop_vertical
    end DO loop_hemisph

  END SUBROUTINE filtreg

end module filtreg_m
1	guez	27	module filtreg_m
2	guez	3
3	guez	27	IMPLICIT NONE
4	guez	3
5	guez	27	contains
6	guez	3
7	guez	107	SUBROUTINE filtreg(champ, direct, intensive)
8	guez	3
9	guez	27	! From filtrez/filtreg.F, version 1.1.1.1, 2004/05/19 12:53:09
10			! Author: P. Le Van
11	guez	133	! Objet : filtre matriciel longitudinal, avec les matrices pr\'ecalcul\'ees
12			! pour l'op\'erateur filtre.
13	guez	3
14	guez	71	USE coefils, ONLY: sddu, sddv, unsddu, unsddv
15	guez	64	USE dimens_m, ONLY: iim, jjm
16	guez	54	use inifilr_m, only: jfiltnu, jfiltnv, jfiltsu, jfiltsv, matriceun, &
17			matriceus, matricevn, matricevs, matrinvn, matrinvs
18	guez	107	use nr_util, only: assert
19	guez	3
20	guez	107	REAL, intent(inout):: champ(:, :, :) ! (iim + 1, nlat, nbniv)
21	guez	133	! en entr\'ee : champ \`a filtrer, en sortie : champ filtr\'e
22	guez	3
23	guez	107	logical, intent(in):: direct ! filtre direct ou inverse
24	guez	3
25	guez	107	logical, intent(in):: intensive
26	guez	133	! champ intensif ou extensif (pond\'er\'e par les aires)
27	guez	3
28	guez	107	! Local:
29			LOGICAL griscal
30	guez	133	INTEGER nlat ! nombre de latitudes \`a filtrer
31			integer nbniv ! nombre de niveaux verticaux \`a filtrer
32	guez	27	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	guez	25
37	guez	27	!-----------------------------------------------------------
38	guez	25
39	guez	107	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	guez	25
45	guez	107	IF (.not. direct .AND. nlat == jjm) THEN
46			PRINT *, 'filtreg: inverse filter on scalar grid only'
47	guez	54	STOP 1
48	guez	27	END IF
49	guez	25
50	guez	27	IF (griscal) THEN
51	guez	107	IF (intensive) THEN
52			sdd1 = sddv
53			sdd2 = unsddv
54	guez	27	ELSE
55	guez	107	sdd1 = unsddv
56			sdd2 = sddv
57			END IF
58	guez	25
59	guez	107	jdfil1 = 2
60			jffil1 = jfiltnu
61			jdfil2 = jfiltsu
62			jffil2 = jjm
63	guez	27	ELSE
64	guez	107	IF (intensive) THEN
65			sdd1 = sddu
66			sdd2 = unsddu
67	guez	27	ELSE
68	guez	107	sdd1 = unsddu
69			sdd2 = sddu
70			END IF
71	guez	25
72	guez	107	jdfil1 = 1
73			jffil1 = jfiltnv
74			jdfil2 = jfiltsv
75			jffil2 = jjm
76	guez	27	END IF
77	guez	25
78	guez	133	loop_hemisph: DO hemisph = 1, 2
79	guez	27	IF (hemisph==1) THEN
80			jdfil = jdfil1
81			jffil = jffil1
82			ELSE
83			jdfil = jdfil2
84			jffil = jffil2
85			END IF
86	guez	25
87	guez	57	loop_vertical: DO l = 1, nbniv
88			loop_latitude: DO j = jdfil, jffil
89	guez	27	DO i = 1, iim
90			champ(i, j, l) = champ(i, j, l)*sdd1(i)
91			END DO
92	guez	25
93	guez	27	IF (hemisph==1) THEN
94	guez	107	IF (.not. direct) THEN
95	guez	27	DO k = 1, iim
96	guez	54	eignq(k) = 0.
97	guez	27	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	guez	54	eignq(k) = 0.
106	guez	27	END DO
107			DO i = 1, iim
108			DO k = 1, iim
109	guez	32	eignq(k) = eignq(k) + matriceun(k, i, j) &
110			* champ(i, j, l)
111	guez	27	END DO
112			END DO
113			ELSE
114			DO k = 1, iim
115	guez	54	eignq(k) = 0.
116	guez	27	END DO
117			DO i = 1, iim
118			DO k = 1, iim
119	guez	32	eignq(k) = eignq(k) + matricevn(k, i, j) &
120			* champ(i, j, l)
121	guez	27	END DO
122			END DO
123			END IF
124			ELSE
125	guez	107	IF (.not. direct) THEN
126	guez	27	DO k = 1, iim
127	guez	54	eignq(k) = 0.
128	guez	27	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	guez	54	eignq(k) = 0.
138	guez	27	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	guez	54	eignq(k) = 0.
148	guez	27	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	guez	25
158	guez	107	IF (direct) THEN
159	guez	27	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	guez	25
168	guez	27	champ(iim + 1, j, l) = champ(1, j, l)
169	guez	57	END DO loop_latitude
170			END DO loop_vertical
171	guez	133	end DO loop_hemisph
172	guez	27
173			END SUBROUTINE filtreg
174
175			end module filtreg_m