Context Navigation

filtreg.F @ 5065

Last change on this file since 5065 was 4593, checked in by yann meurdesoif, 17 months ago

Replace #include (c preprocessor) by INCLUDE (fortran keyword)

in phylmd (except rrtm and ecrad) filtrez, dy3dmem and dyn3dcommon

Other directories will follow
YM

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: 9.3 KB

Rev	Line
[524]	1	!
	2	! $Header$
	3	!
	4	SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire,
[1146]	5	& griscal ,iter)
	6
	7	USE filtreg_mod
	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
[4593]	48	INCLUDE "dimensions.h"
	49	INCLUDE "paramet.h"
	50	INCLUDE "coefils.h"
[524]	51
[1146]	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
[1146]	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
[4440]	74	if ( iim == 1 ) return ! no filtre in 2D y-z
	75
[1146]	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
[524]	85	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
[1146]	86	& STOP'Pas de transformee simple dans cette version'
	87
[524]	88	IF( iter.EQ. 2 ) THEN
[1146]	89	PRINT *,' Pas d iteration du filtre dans cette version !'
	90	& , ' Utiliser old_filtreg et repasser !'
	91	STOP
[524]	92	ENDIF
[1146]	93
[524]	94	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
[1146]	95	PRINT *,' Cette routine ne calcule le filtre inverse que '
	96	& , ' sur la grille des scalaires !'
	97	STOP
[524]	98	ENDIF
[1146]	99
[524]	100	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
[1146]	101	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
	102	& , ' corriger et repasser !'
	103	STOP
[524]	104	ENDIF
[1146]	105
[524]	106	iim2 = iim * iim
	107	immjm = iim * jjm
[1146]	108
[524]	109	IF( griscal ) THEN
	110	IF( nlat. NE. jjp1 ) THEN
[1146]	111	PRINT 1111
	112	STOP
[524]	113	ELSE
[1146]	114
	115	IF( iaire.EQ.1 ) THEN
[1279]	116	sdd1_type = type_sddv
	117	sdd2_type = type_unsddv
[1146]	118	ELSE
[1279]	119	sdd1_type = type_unsddv
	120	sdd2_type = type_sddv
[1146]	121	ENDIF
	122
	123	c IF( iaire.EQ.1 ) THEN
	124	c CALL SCOPY( iim, sddv, 1, sdd1, 1 )
	125	c CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
	126	c ELSE
	127	c CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
	128	c CALL SCOPY( iim, sddv, 1, sdd2, 1 )
	129	c END IF
	130
	131	jdfil1 = 2
	132	jffil1 = jfiltnu
	133	jdfil2 = jfiltsu
	134	jffil2 = jjm
	135	END IF
[524]	136	ELSE
[1146]	137	IF( nlat.NE.jjm ) THEN
	138	PRINT 2222
	139	STOP
	140	ELSE
	141
	142	IF( iaire.EQ.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	c IF( iaire.EQ.1 ) THEN
	151	c CALL SCOPY( iim, sddu, 1, sdd1, 1 )
	152	c CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
	153	c ELSE
	154	c CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
	155	c CALL SCOPY( iim, sddu, 1, sdd2, 1 )
	156	c END IF
	157
	158	jdfil1 = 1
	159	jffil1 = jfiltnv
	160	jdfil2 = jfiltsv
	161	jffil2 = jjm
	162	END IF
[524]	163	END IF
[1146]	164
	165	DO hemisph = 1, 2
	166
	167	IF ( hemisph.EQ.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. EQ. 1 ) THEN
	184
	185	IF( ifiltre. EQ. -2 ) THEN
	186
	187	DO j = jdfil,jffil
[524]	188	#ifdef BLAS
[1705]	189	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]	190	& matrinvn(1,1,j),
	191	& iim, champ(1,j,1), iip1*nlat, 0.0,
	192	& eignq(1,j-jdfil+1,1), iim*nlat)
[524]	193	#else
[1146]	194	eignq(:,j-jdfil+1,:)
	195	$ = matmul(matrinvn(:,:,j), champ(:iim,j,:))
[524]	196	#endif
[1146]	197	END DO
	198
	199	ELSE IF ( griscal ) THEN
	200
	201	DO j = jdfil,jffil
[524]	202	#ifdef BLAS
[1705]	203	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]	204	& matriceun(1,1,j),
	205	& iim, champ(1,j,1), iip1*nlat, 0.0,
	206	& eignq(1,j-jdfil+1,1), iim*nlat)
[524]	207	#else
[1146]	208	eignq(:,j-jdfil+1,:)
	209	$ = matmul(matriceun(:,:,j), champ(:iim,j,:))
[524]	210	#endif
[1146]	211	END DO
	212
	213	ELSE
	214
	215	DO j = jdfil,jffil
[524]	216	#ifdef BLAS
[1705]	217	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]	218	& matricevn(1,1,j),
	219	& iim, champ(1,j,1), iip1*nlat, 0.0,
	220	& eignq(1,j-jdfil+1,1), iim*nlat)
[524]	221	#else
[1146]	222	eignq(:,j-jdfil+1,:)
	223	$ = matmul(matricevn(:,:,j), champ(:iim,j,:))
[524]	224	#endif
[1146]	225	END DO
	226
	227	ENDIF
	228
	229	ELSE
	230
	231	IF( ifiltre. EQ. -2 ) THEN
	232
	233	DO j = jdfil,jffil
[524]	234	#ifdef BLAS
[1705]	235	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]	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)
[524]	239	#else
[1146]	240	eignq(:,j-jdfil+1,:)
	241	$ = matmul(matrinvs(:,:,j-jfiltsu+1),
	242	$ champ(:iim,j,:))
[524]	243	#endif
[1146]	244	END DO
	245
	246
	247	ELSE IF ( griscal ) THEN
	248
	249	DO j = jdfil,jffil
[524]	250	#ifdef BLAS
[1705]	251	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]	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)
[524]	255	#else
[1146]	256	eignq(:,j-jdfil+1,:)
	257	$ = matmul(matriceus(:,:,j-jfiltsu+1),
	258	$ champ(:iim,j,:))
[524]	259	#endif
[1146]	260	END DO
	261
	262	ELSE
	263
	264	DO j = jdfil,jffil
[524]	265	#ifdef BLAS
[1705]	266	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
[1146]	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)
[524]	270	#else
[1146]	271	eignq(:,j-jdfil+1,:)
	272	$ = matmul(matricevs(:,:,j-jfiltsv+1),
	273	$ champ(:iim,j,:))
[524]	274	#endif
[1146]	275	END DO
	276
	277	ENDIF
	278
	279	ENDIF
	280
	281	IF( ifiltre.EQ. 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
[524]	292
[1146]	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
[524]	316	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
[1146]	317	& filtrer, sur la grille des scalaires'/)
[524]	318	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
[1146]	319	& ltrer, sur la grille de V ou de Z'/)
[524]	320	RETURN
	321	END

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

Download in other formats:

Original Format