Context Navigation

filtreg.F @ 575

Last change on this file since 575 was 495, checked in by lmdzadmin, 21 years ago
Optimisation de differentes routines, IM, MAF, FH LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 10.7 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
[495]	51	INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil(2),jffil(2)
[2]	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)
[495]	58	cIM REAL eignq(iim), sdd1(iim),sdd2(iim)
	59	REAL eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
[2]	60	LOGICAL griscal
	61	INTEGER hemisph, iaire
	62	c
	63
	64	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
	65	* STOP'Pas de transformee simple dans cette version'
	66
	67	IF( iter.EQ. 2 ) THEN
	68	PRINT *,' Pas d iteration du filtre dans cette version !'
	69	* , ' Utiliser old_filtreg et repasser !'
	70	STOP
	71	ENDIF
	72
	73	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
	74	PRINT *,' Cette routine ne calcule le filtre inverse que ',
	75	* ' sur la grille des scalaires !'
	76	STOP
	77	ENDIF
	78
	79	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
	80	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
	81	*,' corriger et repasser !'
	82	STOP
	83	ENDIF
	84	c
	85
	86	iim2 = iim * iim
	87	immjm = iim * jjm
	88	c
	89	c
	90	IF( griscal ) THEN
	91	IF( nlat. NE. jjp1 ) THEN
	92	PRINT 1111
	93	STOP
	94	ELSE
	95	c
	96	IF( iaire.EQ.1 ) THEN
	97	CALL SCOPY( iim, sddv, 1, sdd1, 1 )
	98	CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
	99	ELSE
	100	CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
	101	CALL SCOPY( iim, sddv, 1, sdd2, 1 )
	102	END IF
	103	c
	104	jdfil1 = 2
	105	jffil1 = jfiltnu
	106	jdfil2 = jfiltsu
	107	jffil2 = jjm
	108	END IF
	109	ELSE
	110	IF( nlat.NE.jjm ) THEN
	111	PRINT 2222
	112	STOP
	113	ELSE
	114	c
	115	IF( iaire.EQ.1 ) THEN
	116	CALL SCOPY( iim, sddu, 1, sdd1, 1 )
	117	CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
	118	ELSE
	119	CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
	120	CALL SCOPY( iim, sddu, 1, sdd2, 1 )
	121	END IF
	122	c
	123	jdfil1 = 1
	124	jffil1 = jfiltnv
	125	jdfil2 = jfiltsv
	126	jffil2 = jjm
	127	END IF
	128	END IF
	129	c
[495]	130	jdfil(1) = jdfil1
	131	jffil(1) = jffil1
	132	jdfil(2) = jdfil2
	133	jffil(2) = jffil2
[2]	134	c
	135	DO 100 hemisph = 1, 2
	136	c
[495]	137	c IF ( hemisph.EQ.1 ) THEN
	138	c jdfil = jdfil1
	139	c jffil = jffil1
	140	c ELSE
	141	c jdfil = jdfil2
	142	c jffil = jffil2
	143	c END IF
[2]	144
	145
	146	DO 50 l = 1, nbniv
[495]	147	DO 30 j = jdfil(hemisph),jffil(hemisph)
[2]	148
	149
	150	DO 5 i = 1, iim
	151	champ(i,j,l) = champ(i,j,l) * sdd1(i)
	152	5 CONTINUE
	153	c
[495]	154	30 CONTINUE
	155	50 CONTINUE
[2]	156
	157	IF( hemisph. EQ. 1 ) THEN
	158
	159	IF( ifiltre. EQ. -2 ) THEN
	160	#ifdef CRAY
[495]	161	DO l = 1, nbniv
	162	DO j = jdfil(hemisph),jffil(hemisph)
[2]	163	CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
	164	* 1, iim, iim )
[495]	165	ENDDO
	166	ENDDO
[2]	167	#else
	168	#ifdef BLAS
[495]	169	DO l = 1, nbniv
	170	DO j = jdfil(hemisph),jffil(hemisph)
[2]	171	CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
	172	. champ(1,j,l), 1, 0.0, eignq, 1)
[495]	173	ENDDO
	174	ENDDO
[2]	175	#else
[495]	176	DO l = 1, nbniv
	177	DO j = jdfil(hemisph),jffil(hemisph)
[2]	178	DO k = 1, iim
[495]	179	c eignq(k) = 0.0
	180	eignq(k,j,l) = 0.0
[2]	181	ENDDO
	182	DO k = 1, iim
	183	DO i = 1, iim
[495]	184	c eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
	185	eignq(k,j,l) = eignq(k,j,l) + matrinvn(k,i,j)*champ(i,j,l)
[2]	186	ENDDO
	187	ENDDO
[495]	188	ENDDO
	189	ENDDO
[2]	190	#endif
	191	#endif
	192	ELSE IF ( griscal ) THEN
	193	#ifdef CRAY
[495]	194	DO l = 1, nbniv
	195	DO j = jdfil(hemisph),jffil(hemisph)
[2]	196	CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
	197	* 1, iim, iim )
[495]	198	ENDDO
	199	ENDDO
[2]	200	#else
	201	#ifdef BLAS
[495]	202	DO l = 1, nbniv
	203	DO j = jdfil(hemisph),jffil(hemisph)
[2]	204	CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
	205	. champ(1,j,l), 1, 0.0, eignq, 1)
[495]	206	ENDDO
	207	ENDDO
[2]	208	#else
[495]	209	DO l = 1, nbniv
	210	DO j = jdfil(hemisph),jffil(hemisph)
[2]	211	DO k = 1, iim
[495]	212	c eignq(k) = 0.0
	213	eignq(k,j,l) = 0.0
[2]	214	ENDDO
	215	DO i = 1, iim
	216	DO k = 1, iim
[495]	217	c eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
	218	eignq(k,j,l) = eignq(k,j,l) + matriceun(k,i,j)*champ(i,j,l)
[2]	219	ENDDO
	220	ENDDO
[495]	221	ENDDO
	222	ENDDO
[2]	223	#endif
	224	#endif
	225	ELSE
	226	#ifdef CRAY
[495]	227	DO l = 1, nbniv
	228	DO j = jdfil(hemisph),jffil(hemisph)
[2]	229	CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
	230	* 1, iim, iim )
[495]	231	ENDDO
	232	ENDDO
[2]	233	#else
	234	#ifdef BLAS
[495]	235	DO l = 1, nbniv
	236	DO j = jdfil(hemisph),jffil(hemisph)
[2]	237	CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
	238	. champ(1,j,l), 1, 0.0, eignq, 1)
[495]	239	ENDDO
	240	ENDDO
[2]	241	#else
[495]	242	DO l = 1, nbniv
	243	DO j = jdfil(hemisph),jffil(hemisph)
[2]	244	DO k = 1, iim
[495]	245	c eignq(k) = 0.0
	246	eignq(k,j,l) = 0.0
[2]	247	ENDDO
	248	DO i = 1, iim
	249	DO k = 1, iim
[495]	250	c eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
	251	eignq(k,j,l) = eignq(k,j,l) + matricevn(k,i,j)*champ(i,j,l)
[2]	252	ENDDO
	253	ENDDO
[495]	254	ENDDO
	255	ENDDO
[2]	256	#endif
	257	#endif
	258	ENDIF
	259
	260	ELSE
	261
	262	IF( ifiltre. EQ. -2 ) THEN
	263	#ifdef CRAY
[495]	264	DO l = 1, nbniv
	265	DO j = jdfil(hemisph),jffil(hemisph)
[2]	266	CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
	267	* eignq, 1, iim, iim )
[495]	268	ENDDO
	269	ENDDO
[2]	270	#else
	271	#ifdef BLAS
[495]	272	DO l = 1, nbniv
	273	DO j = jdfil(hemisph),jffil(hemisph)
[2]	274	CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
	275	. champ(1,j,l), 1, 0.0, eignq, 1)
[495]	276	ENDDO
	277	ENDDO
[2]	278	#else
[495]	279	DO l = 1, nbniv
	280	DO j = jdfil(hemisph),jffil(hemisph)
[2]	281	DO k = 1, iim
[495]	282	c eignq(k) = 0.0
	283	eignq(k,j,l) = 0.0
[2]	284	ENDDO
	285	DO i = 1, iim
	286	DO k = 1, iim
[495]	287	c eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
	288	eignq(k,j,l) = eignq(k,j,l) +
	289	.matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
[2]	290	ENDDO
	291	ENDDO
[495]	292	ENDDO
	293	ENDDO
[2]	294	#endif
	295	#endif
	296	ELSE IF ( griscal ) THEN
	297	#ifdef CRAY
[495]	298	DO l = 1, nbniv
	299	DO j = jdfil(hemisph),jffil(hemisph)
[2]	300	CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
	301	* eignq, 1, iim, iim )
[495]	302	ENDDO
	303	ENDDO
[2]	304	#else
	305	#ifdef BLAS
[495]	306	DO l = 1, nbniv
	307	DO j = jdfil(hemisph),jffil(hemisph)
[2]	308	CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
	309	. champ(1,j,l), 1, 0.0, eignq, 1)
[495]	310	ENDDO
	311	ENDDO
[2]	312	#else
[495]	313	DO l = 1, nbniv
	314	DO j = jdfil(hemisph),jffil(hemisph)
[2]	315	DO k = 1, iim
[495]	316	c eignq(k) = 0.0
	317	eignq(k,j,l) = 0.0
[2]	318	ENDDO
	319	DO i = 1, iim
	320	DO k = 1, iim
[495]	321	c eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
	322	eignq(k,j,l) = eignq(k,j,l) +
	323	.matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
[2]	324	ENDDO
	325	ENDDO
[495]	326	ENDDO
	327	ENDDO
[2]	328	#endif
	329	#endif
	330	ELSE
	331	#ifdef CRAY
[495]	332	DO l = 1, nbniv
	333	DO j = jdfil(hemisph),jffil(hemisph)
[2]	334	CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
	335	* eignq, 1, iim, iim )
[495]	336	ENDDO
	337	ENDDO
[2]	338	#else
	339	#ifdef BLAS
[495]	340	DO l = 1, nbniv
	341	DO j = jdfil(hemisph),jffil(hemisph)
[2]	342	CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
	343	. champ(1,j,l), 1, 0.0, eignq, 1)
[495]	344	ENDDO
	345	ENDDO
[2]	346	#else
[495]	347	DO l = 1, nbniv
	348	DO j = jdfil(hemisph),jffil(hemisph)
[2]	349	DO k = 1, iim
[495]	350	c eignq(k) = 0.0
	351	eignq(k,j,l) = 0.0
[2]	352	ENDDO
	353	DO i = 1, iim
	354	DO k = 1, iim
[495]	355	c eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
	356	eignq(k,j,l) = eignq(k,j,l) +
	357	.matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
[2]	358	ENDDO
	359	ENDDO
[495]	360	ENDDO
	361	ENDDO
[2]	362	#endif
	363	#endif
	364	ENDIF
	365
	366	ENDIF
	367	c
	368	IF( ifiltre.EQ. 2 ) THEN
[495]	369	DO l = 1, nbniv
	370	DO j = jdfil(hemisph),jffil(hemisph)
[2]	371	DO 15 i = 1, iim
[495]	372	c champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
	373	champ( i,j,l ) = ( champ(i,j,l) + eignq(i,j,l) ) * sdd2(i)
[2]	374	15 CONTINUE
[495]	375	ENDDO
	376	ENDDO
[2]	377	ELSE
[495]	378	DO l = 1, nbniv
	379	DO j = jdfil(hemisph),jffil(hemisph)
[2]	380	DO 16 i=1,iim
[495]	381	c champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
	382	champ( i,j,l ) = ( champ(i,j,l) - eignq(i,j,l) ) * sdd2(i)
[2]	383	16 CONTINUE
[495]	384	ENDDO
	385	ENDDO
[2]	386	ENDIF
	387	c
[495]	388	DO l = 1, nbniv
	389	DO j = jdfil(hemisph),jffil(hemisph)
[2]	390	champ( iip1,j,l ) = champ( 1,j,l )
[495]	391	ENDDO
	392	ENDDO
[2]	393	c
[495]	394	c 30 CONTINUE
[2]	395	c
[495]	396	c 50 CONTINUE
[2]	397	c
	398	100 CONTINUE
	399	c
[231]	400	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
	401	*filtrer, sur la grille des scalaires'/)
	402	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
	403	*ltrer, sur la grille de V ou de Z'/)
[2]	404	RETURN
	405	END

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

Download in other formats:

Original Format