Context Navigation

filtreg.F @ 320

Last change on this file since 320 was 231, checked in by lmdzadmin, 23 years ago
Merge avec la branche principale LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 7.8 KB

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

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

Download in other formats:

Original Format