Context Navigation

filtreg.F90 @ 5105

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

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

Download in other formats:

Original Format