Context Navigation

filtreg.F90 @ 5272

Last change on this file since 5272 was 5272, checked in by abarral, 8 weeks ago
Turn paramet.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.4 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, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
11	ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
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
52
53	INCLUDE "coefils.h"
54
55	INTEGER :: nlat,nbniv,ifiltre,iter
56	INTEGER :: i,j,l,k
57	INTEGER :: iim2,immjm
58	INTEGER :: jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
59
60	REAL :: champ( iip1,nlat,nbniv)
61
62	REAL :: eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
63	LOGICAL :: griscal
64	INTEGER :: hemisph, iaire
65
66	LOGICAL,SAVE :: first=.TRUE.
67
68	REAL, SAVE :: sdd12(iim,4)
69
70	INTEGER, PARAMETER :: type_sddu=1
71	INTEGER, PARAMETER :: type_sddv=2
72	INTEGER, PARAMETER :: type_unsddu=3
73	INTEGER, PARAMETER :: type_unsddv=4
74
75	INTEGER :: sdd1_type, sdd2_type
76
77	if ( iim == 1 ) return ! no filtre in 2D y-z
78
79	IF (first) THEN
80	sdd12(1:iim,type_sddu) = sddu(1:iim)
81	sdd12(1:iim,type_sddv) = sddv(1:iim)
82	sdd12(1:iim,type_unsddu) = unsddu(1:iim)
83	sdd12(1:iim,type_unsddv) = unsddv(1:iim)
84
85	first=.FALSE.
86	ENDIF
87
88	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) &
89	STOP 'Pas de transformee simple dans cette version'
90
91	IF( iter.EQ. 2 ) THEN
92	PRINT *,' Pas d iteration du filtre dans cette version !'&
93	& , ' Utiliser old_filtreg et repasser !'
94	STOP
95	ENDIF
96
97	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
98	PRINT *,' Cette routine ne calcule le filtre inverse que ' &
99	, ' sur la grille des scalaires !'
100	STOP
101	ENDIF
102
103	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
104	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
105	, ' corriger et repasser !'
106	STOP
107	ENDIF
108
109	iim2 = iim * iim
110	immjm = iim * jjm
111
112	IF( griscal ) THEN
113	IF( nlat.NE. jjp1 ) THEN
114	PRINT 1111
115	STOP
116	ELSE
117
118	IF( iaire.EQ.1 ) THEN
119	sdd1_type = type_sddv
120	sdd2_type = type_unsddv
121	ELSE
122	sdd1_type = type_unsddv
123	sdd2_type = type_sddv
124	ENDIF
125
126	! IF( iaire.EQ.1 ) THEN
127	! CALL SCOPY( iim, sddv, 1, sdd1, 1 )
128	! CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
129	! ELSE
130	! CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
131	! CALL SCOPY( iim, sddv, 1, sdd2, 1 )
132	! END IF
133
134	jdfil1 = 2
135	jffil1 = jfiltnu
136	jdfil2 = jfiltsu
137	jffil2 = jjm
138	END IF
139	ELSE
140	IF( nlat.NE.jjm ) THEN
141	PRINT 2222
142	STOP
143	ELSE
144
145	IF( iaire.EQ.1 ) THEN
146	sdd1_type = type_sddu
147	sdd2_type = type_unsddu
148	ELSE
149	sdd1_type = type_unsddu
150	sdd2_type = type_sddu
151	ENDIF
152
153	! IF( iaire.EQ.1 ) THEN
154	! CALL SCOPY( iim, sddu, 1, sdd1, 1 )
155	! CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
156	! ELSE
157	! CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
158	! CALL SCOPY( iim, sddu, 1, sdd2, 1 )
159	! END IF
160
161	jdfil1 = 1
162	jffil1 = jfiltnv
163	jdfil2 = jfiltsv
164	jffil2 = jjm
165	END IF
166	END IF
167
168	DO hemisph = 1, 2
169
170	IF ( hemisph.EQ.1 ) THEN
171	jdfil = jdfil1
172	jffil = jffil1
173	ELSE
174	jdfil = jdfil2
175	jffil = jffil2
176	END IF
177
178	DO l = 1, nbniv
179	DO j = jdfil,jffil
180	DO i = 1, iim
181	champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i)
182	END DO
183	END DO
184	END DO
185
186	IF( hemisph.EQ. 1 ) THEN
187
188	IF( ifiltre.EQ. -2 ) THEN
189
190	DO j = jdfil,jffil
191	#ifdef BLAS
192	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
193	matrinvn(1,1,j), &
194	iim, champ(1,j,1), iip1*nlat, 0.0, &
195	eignq(1,j-jdfil+1,1), iim*nlat)
196	#else
197	eignq(:,j-jdfil+1,:) &
198	= matmul(matrinvn(:,:,j), champ(:iim,j,:))
199	#endif
200	END DO
201
202	ELSE IF ( griscal ) THEN
203
204	DO j = jdfil,jffil
205	#ifdef BLAS
206	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
207	matriceun(1,1,j), &
208	iim, champ(1,j,1), iip1*nlat, 0.0, &
209	eignq(1,j-jdfil+1,1), iim*nlat)
210	#else
211	eignq(:,j-jdfil+1,:) &
212	= matmul(matriceun(:,:,j), champ(:iim,j,:))
213	#endif
214	END DO
215
216	ELSE
217
218	DO j = jdfil,jffil
219	#ifdef BLAS
220	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
221	matricevn(1,1,j), &
222	iim, champ(1,j,1), iip1*nlat, 0.0, &
223	eignq(1,j-jdfil+1,1), iim*nlat)
224	#else
225	eignq(:,j-jdfil+1,:) &
226	= matmul(matricevn(:,:,j), champ(:iim,j,:))
227	#endif
228	END DO
229
230	ENDIF
231
232	ELSE
233
234	IF( ifiltre.EQ. -2 ) THEN
235
236	DO j = jdfil,jffil
237	#ifdef BLAS
238	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
239	matrinvs(1,1,j-jfiltsu+1), &
240	iim, champ(1,j,1), iip1*nlat, 0.0, &
241	eignq(1,j-jdfil+1,1), iim*nlat)
242	#else
243	eignq(:,j-jdfil+1,:) &
244	= matmul(matrinvs(:,:,j-jfiltsu+1), &
245	champ(:iim,j,:))
246	#endif
247	END DO
248
249
250	ELSE IF ( griscal ) THEN
251
252	DO j = jdfil,jffil
253	#ifdef BLAS
254	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
255	matriceus(1,1,j-jfiltsu+1), &
256	iim, champ(1,j,1), iip1*nlat, 0.0, &
257	eignq(1,j-jdfil+1,1), iim*nlat)
258	#else
259	eignq(:,j-jdfil+1,:) &
260	= matmul(matriceus(:,:,j-jfiltsu+1), &
261	champ(:iim,j,:))
262	#endif
263	END DO
264
265	ELSE
266
267	DO j = jdfil,jffil
268	#ifdef BLAS
269	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
270	matricevs(1,1,j-jfiltsv+1), &
271	iim, champ(1,j,1), iip1*nlat, 0.0, &
272	eignq(1,j-jdfil+1,1), iim*nlat)
273	#else
274	eignq(:,j-jdfil+1,:) &
275	= matmul(matricevs(:,:,j-jfiltsv+1), &
276	champ(:iim,j,:))
277	#endif
278	END DO
279
280	ENDIF
281
282	ENDIF
283
284	IF( ifiltre.EQ. 2 ) THEN
285
286	DO l = 1, nbniv
287	DO j = jdfil,jffil
288	DO i = 1, iim
289	champ( i,j,l ) = &
290	(champ(i,j,l) + eignq(i,j-jdfil+1,l)) &
291	* sdd12(i,sdd2_type) ! sdd2(i)
292	END DO
293	END DO
294	END DO
295
296	ELSE
297
298	DO l = 1, nbniv
299	DO j = jdfil,jffil
300	DO i = 1, iim
301	champ( i,j,l ) = &
302	(champ(i,j,l) - eignq(i,j-jdfil+1,l)) &
303	* sdd12(i,sdd2_type) ! sdd2(i)
304	END DO
305	END DO
306	END DO
307
308	ENDIF
309
310	DO l = 1, nbniv
311	DO j = jdfil,jffil
312	champ( iip1,j,l ) = champ( 1,j,l )
313	END DO
314	END DO
315
316
317	ENDDO
318
319	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a&
320	& filtrer, sur la grille des scalaires'/)
321	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi&
322	& ltrer, sur la grille de V ou de Z'/)
323	RETURN
324	END SUBROUTINE filtreg

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

Download in other formats:

Original Format