Context Navigation

filtreg_p.F @ 883

Last change on this file since 883 was 763, checked in by Laurent Fairhead, 17 years ago
Merge entre la version V3_conv et le HEAD Nouvelles routines parallélisées YM, JG, LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 8.1 KB

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

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

Download in other formats:

Original Format