Context Navigation

filtreg.F @ 1119

Last change on this file since 1119 was 1108, checked in by yann meurdesoif, 16 years ago

module filtre -> module filtreg_mod ; use filtre -> use filtreg_mod

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 9.2 KB

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

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

Download in other formats:

Original Format