trunk/filtrez/filtreg.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écalculées
    ! pour l'opérateur 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ée : champ à filtrer, en sortie : champ filtré

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

    logical, intent(in):: intensive
    ! champ intensif ou extensif (pondéré par les aires)

    ! Local:
    LOGICAL griscal
    INTEGER nlat ! nombre de latitudes à filtrer
    integer nbniv ! nombre de niveaux verticaux à 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

    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

  END SUBROUTINE filtreg

end module filtreg_m
1	module filtreg_m
2
3	IMPLICIT NONE
4
5	contains
6
7	SUBROUTINE filtreg(champ, direct, intensive)
8
9	! From filtrez/filtreg.F, version 1.1.1.1, 2004/05/19 12:53:09
10	! 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	END SUBROUTINE filtreg
174
175	end module filtreg_m