Context Navigation

filtreg.F90 @ 5348

Last change on this file since 5348 was 5297, checked in by abarral, 4 weeks ago
Turn gradsdef.h coefils.h into a module
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 8.2 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire, &
5	griscal ,iter)
6
7	USE filtreg_mod
8
9	USE dimensions_mod, ONLY: iim, jjm, llm, ndm
10	USE paramet_mod_h
11	USE coefils_mod_h
12	IMPLICIT NONE
13	!=======================================================================
14	!
15	! Auteur: P. Le Van 07/10/97
16	! ------
17	!
18	! Objet: filtre matriciel longitudinal ,avec les matrices precalculees
19	! pour l'operateur Filtre .
20	! ------
21	!
22	! Arguments:
23	! ----------
24	!
25	! nblat nombre de latitudes a filtrer
26	! nbniv nombre de niveaux verticaux a filtrer
27	! champ(iip1,nblat,nbniv) en entree : champ a filtrer
28	! en sortie : champ filtre
29	! ifiltre +1 Transformee directe
30	! -1 Transformee inverse
31	! +2 Filtre directe
32	! -2 Filtre inverse
33	!
34	! iaire 1 si champ intensif
35	! 2 si champ extensif (pondere par les aires)
36	!
37	! iter 1 filtre simple
38	!
39	!=======================================================================
40	!
41	!
42	! Variable Intensive
43	! ifiltre = 1 filtre directe
44	! ifiltre =-1 filtre inverse
45	!
46	! Variable Extensive
47	! ifiltre = 2 filtre directe
48	! ifiltre =-2 filtre inverse
49	!
50
51	INTEGER :: nlat,nbniv,ifiltre,iter
52	INTEGER :: i,j,l,k
53	INTEGER :: iim2,immjm
54	INTEGER :: jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
55
56	REAL :: champ( iip1,nlat,nbniv)
57
58	REAL :: eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
59	LOGICAL :: griscal
60	INTEGER :: hemisph, iaire
61
62	LOGICAL,SAVE :: first=.TRUE.
63
64	REAL, SAVE :: sdd12(iim,4)
65
66	INTEGER, PARAMETER :: type_sddu=1
67	INTEGER, PARAMETER :: type_sddv=2
68	INTEGER, PARAMETER :: type_unsddu=3
69	INTEGER, PARAMETER :: type_unsddv=4
70
71	INTEGER :: sdd1_type, sdd2_type
72
73	if ( iim == 1 ) return ! no filtre in 2D y-z
74
75	IF (first) THEN
76	sdd12(1:iim,type_sddu) = sddu(1:iim)
77	sdd12(1:iim,type_sddv) = sddv(1:iim)
78	sdd12(1:iim,type_unsddu) = unsddu(1:iim)
79	sdd12(1:iim,type_unsddv) = unsddv(1:iim)
80
81	first=.FALSE.
82	ENDIF
83
84	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) &
85	STOP 'Pas de transformee simple dans cette version'
86
87	IF( iter.EQ. 2 ) THEN
88	PRINT *,' Pas d iteration du filtre dans cette version !'&
89	& , ' Utiliser old_filtreg et repasser !'
90	STOP
91	ENDIF
92
93	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
94	PRINT *,' Cette routine ne calcule le filtre inverse que ' &
95	, ' sur la grille des scalaires !'
96	STOP
97	ENDIF
98
99	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
100	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
101	, ' corriger et repasser !'
102	STOP
103	ENDIF
104
105	iim2 = iim * iim
106	immjm = iim * jjm
107
108	IF( griscal ) THEN
109	IF( nlat.NE. jjp1 ) THEN
110	PRINT 1111
111	STOP
112	ELSE
113
114	IF( iaire.EQ.1 ) THEN
115	sdd1_type = type_sddv
116	sdd2_type = type_unsddv
117	ELSE
118	sdd1_type = type_unsddv
119	sdd2_type = type_sddv
120	ENDIF
121
122	! IF( iaire.EQ.1 ) THEN
123	! CALL SCOPY( iim, sddv, 1, sdd1, 1 )
124	! CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
125	! ELSE
126	! CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
127	! CALL SCOPY( iim, sddv, 1, sdd2, 1 )
128	! END IF
129
130	jdfil1 = 2
131	jffil1 = jfiltnu
132	jdfil2 = jfiltsu
133	jffil2 = jjm
134	END IF
135	ELSE
136	IF( nlat.NE.jjm ) THEN
137	PRINT 2222
138	STOP
139	ELSE
140
141	IF( iaire.EQ.1 ) THEN
142	sdd1_type = type_sddu
143	sdd2_type = type_unsddu
144	ELSE
145	sdd1_type = type_unsddu
146	sdd2_type = type_sddu
147	ENDIF
148
149	! IF( iaire.EQ.1 ) THEN
150	! CALL SCOPY( iim, sddu, 1, sdd1, 1 )
151	! CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
152	! ELSE
153	! CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
154	! CALL SCOPY( iim, sddu, 1, sdd2, 1 )
155	! END IF
156
157	jdfil1 = 1
158	jffil1 = jfiltnv
159	jdfil2 = jfiltsv
160	jffil2 = jjm
161	END IF
162	END IF
163
164	DO hemisph = 1, 2
165
166	IF ( hemisph.EQ.1 ) THEN
167	jdfil = jdfil1
168	jffil = jffil1
169	ELSE
170	jdfil = jdfil2
171	jffil = jffil2
172	END IF
173
174	DO l = 1, nbniv
175	DO j = jdfil,jffil
176	DO i = 1, iim
177	champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i)
178	END DO
179	END DO
180	END DO
181
182	IF( hemisph.EQ. 1 ) THEN
183
184	IF( ifiltre.EQ. -2 ) THEN
185
186	DO j = jdfil,jffil
187	#ifdef BLAS
188	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
189	matrinvn(1,1,j), &
190	iim, champ(1,j,1), iip1*nlat, 0.0, &
191	eignq(1,j-jdfil+1,1), iim*nlat)
192	#else
193	eignq(:,j-jdfil+1,:) &
194	= matmul(matrinvn(:,:,j), champ(:iim,j,:))
195	#endif
196	END DO
197
198	ELSE IF ( griscal ) THEN
199
200	DO j = jdfil,jffil
201	#ifdef BLAS
202	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
203	matriceun(1,1,j), &
204	iim, champ(1,j,1), iip1*nlat, 0.0, &
205	eignq(1,j-jdfil+1,1), iim*nlat)
206	#else
207	eignq(:,j-jdfil+1,:) &
208	= matmul(matriceun(:,:,j), champ(:iim,j,:))
209	#endif
210	END DO
211
212	ELSE
213
214	DO j = jdfil,jffil
215	#ifdef BLAS
216	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
217	matricevn(1,1,j), &
218	iim, champ(1,j,1), iip1*nlat, 0.0, &
219	eignq(1,j-jdfil+1,1), iim*nlat)
220	#else
221	eignq(:,j-jdfil+1,:) &
222	= matmul(matricevn(:,:,j), champ(:iim,j,:))
223	#endif
224	END DO
225
226	ENDIF
227
228	ELSE
229
230	IF( ifiltre.EQ. -2 ) THEN
231
232	DO j = jdfil,jffil
233	#ifdef BLAS
234	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
235	matrinvs(1,1,j-jfiltsu+1), &
236	iim, champ(1,j,1), iip1*nlat, 0.0, &
237	eignq(1,j-jdfil+1,1), iim*nlat)
238	#else
239	eignq(:,j-jdfil+1,:) &
240	= matmul(matrinvs(:,:,j-jfiltsu+1), &
241	champ(:iim,j,:))
242	#endif
243	END DO
244
245
246	ELSE IF ( griscal ) THEN
247
248	DO j = jdfil,jffil
249	#ifdef BLAS
250	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
251	matriceus(1,1,j-jfiltsu+1), &
252	iim, champ(1,j,1), iip1*nlat, 0.0, &
253	eignq(1,j-jdfil+1,1), iim*nlat)
254	#else
255	eignq(:,j-jdfil+1,:) &
256	= matmul(matriceus(:,:,j-jfiltsu+1), &
257	champ(:iim,j,:))
258	#endif
259	END DO
260
261	ELSE
262
263	DO j = jdfil,jffil
264	#ifdef BLAS
265	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
266	matricevs(1,1,j-jfiltsv+1), &
267	iim, champ(1,j,1), iip1*nlat, 0.0, &
268	eignq(1,j-jdfil+1,1), iim*nlat)
269	#else
270	eignq(:,j-jdfil+1,:) &
271	= matmul(matricevs(:,:,j-jfiltsv+1), &
272	champ(:iim,j,:))
273	#endif
274	END DO
275
276	ENDIF
277
278	ENDIF
279
280	IF( ifiltre.EQ. 2 ) THEN
281
282	DO l = 1, nbniv
283	DO j = jdfil,jffil
284	DO i = 1, iim
285	champ( i,j,l ) = &
286	(champ(i,j,l) + eignq(i,j-jdfil+1,l)) &
287	* sdd12(i,sdd2_type) ! sdd2(i)
288	END DO
289	END DO
290	END DO
291
292	ELSE
293
294	DO l = 1, nbniv
295	DO j = jdfil,jffil
296	DO i = 1, iim
297	champ( i,j,l ) = &
298	(champ(i,j,l) - eignq(i,j-jdfil+1,l)) &
299	* sdd12(i,sdd2_type) ! sdd2(i)
300	END DO
301	END DO
302	END DO
303
304	ENDIF
305
306	DO l = 1, nbniv
307	DO j = jdfil,jffil
308	champ( iip1,j,l ) = champ( 1,j,l )
309	END DO
310	END DO
311
312
313	ENDDO
314
315	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a&
316	& filtrer, sur la grille des scalaires'/)
317	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi&
318	& ltrer, sur la grille de V ou de Z'/)
319	RETURN
320	END SUBROUTINE filtreg

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format