Context Navigation

filtreg_p.F @ 1208

Last change on this file since 1208 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: 11.9 KB

Rev	Line
[985]	1
	2
[763]	3	SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
[1086]	4	& ifiltre, iaire, griscal ,iter)
	5	USE Parallel, only : OMP_CHUNK
[985]	6	USE mod_filtre_fft
	7	USE timer_filtre
[1086]	8
[1108]	9	USE filtreg_mod
[1086]	10
[763]	11	IMPLICIT NONE
[1086]	12
[763]	13	c=======================================================================
	14	c
	15	c Auteur: P. Le Van 07/10/97
	16	c ------
	17	c
	18	c Objet: filtre matriciel longitudinal ,avec les matrices precalculees
	19	c pour l'operateur Filtre .
	20	c ------
	21	c
	22	c Arguments:
	23	c ----------
	24	c
	25	c
	26	c ibeg..iend lattitude a filtrer
	27	c nlat nombre de latitudes du champ
	28	c nbniv nombre de niveaux verticaux a filtrer
	29	c champ(iip1,nblat,nbniv) en entree : champ a filtrer
	30	c en sortie : champ filtre
	31	c ifiltre +1 Transformee directe
	32	c -1 Transformee inverse
	33	c +2 Filtre directe
	34	c -2 Filtre inverse
	35	c
	36	c iaire 1 si champ intensif
	37	c 2 si champ extensif (pondere par les aires)
	38	c
	39	c iter 1 filtre simple
	40	c
	41	c=======================================================================
	42	c
	43	c
	44	c Variable Intensive
	45	c ifiltre = 1 filtre directe
	46	c ifiltre =-1 filtre inverse
	47	c
	48	c Variable Extensive
	49	c ifiltre = 2 filtre directe
	50	c ifiltre =-2 filtre inverse
	51	c
	52	c
	53	#include "dimensions.h"
	54	#include "paramet.h"
	55	#include "coefils.h"
	56	c
	57	INTEGER ibeg,iend,nlat,nbniv,ifiltre,iter
	58	INTEGER i,j,l,k
	59	INTEGER iim2,immjm
	60	INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
[1086]	61
[763]	62	REAL champ( iip1,nlat,nbniv)
[985]	63
[763]	64	LOGICAL griscal
	65	INTEGER hemisph, iaire
[985]	66
	67	REAL :: champ_fft(iip1,nlat,nbniv)
	68	REAL :: champ_in(iip1,nlat,nbniv)
	69
	70	LOGICAL,SAVE :: first=.TRUE.
	71	c$OMP THREADPRIVATE(first)
[763]	72
[1086]	73	REAL, DIMENSION(iip1,nlat,nbniv) :: champ_loc
	74	INTEGER :: ll_nb, nbniv_loc
	75	REAL, SAVE :: sdd12(iim,4)
	76	c$OMP THREADPRIVATE(sdd12)
	77
	78	INTEGER, PARAMETER :: type_sddu=1
	79	INTEGER, PARAMETER :: type_sddv=2
	80	INTEGER, PARAMETER :: type_unsddu=3
	81	INTEGER, PARAMETER :: type_unsddv=4
	82
	83	INTEGER :: sdd1_type, sdd2_type
	84
[985]	85	IF (first) THEN
[1086]	86	sdd12(1:iim,type_sddu) = sddu(1:iim)
	87	sdd12(1:iim,type_sddv) = sddv(1:iim)
	88	sdd12(1:iim,type_unsddu) = unsddu(1:iim)
	89	sdd12(1:iim,type_unsddv) = unsddv(1:iim)
	90
	91	CALL Init_timer
	92	first=.FALSE.
[985]	93	ENDIF
	94
	95	c$OMP MASTER
	96	CALL start_timer
	97	c$OMP END MASTER
	98
[1086]	99	c-------------------------------------------------------c
	100
[763]	101	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
[1086]	102	& STOP'Pas de transformee simple dans cette version'
	103
[763]	104	IF( iter.EQ. 2 ) THEN
[1086]	105	PRINT *,' Pas d iteration du filtre dans cette version !'
	106	& , ' Utiliser old_filtreg et repasser !'
	107	STOP
[763]	108	ENDIF
	109
	110	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
[1086]	111	PRINT *,' Cette routine ne calcule le filtre inverse que '
	112	& , ' sur la grille des scalaires !'
	113	STOP
[763]	114	ENDIF
	115
	116	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
[1086]	117	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
	118	& , ' corriger et repasser !'
	119	STOP
[763]	120	ENDIF
	121	c
	122
	123	iim2 = iim * iim
	124	immjm = iim * jjm
	125	c
	126	c
	127	IF( griscal ) THEN
	128	IF( nlat. NE. jjp1 ) THEN
[1086]	129	PRINT 1111
	130	STOP
[763]	131	ELSE
[1086]	132	c
	133	IF( iaire.EQ.1 ) THEN
	134	sdd1_type = type_sddv
	135	sdd2_type = type_unsddv
	136	ELSE
	137	sdd1_type = type_unsddv
	138	sdd2_type = type_sddv
	139	ENDIF
[763]	140	c
[1086]	141	jdfil1 = 2
	142	jffil1 = jfiltnu
	143	jdfil2 = jfiltsu
	144	jffil2 = jjm
	145	ENDIF
[763]	146	ELSE
[1086]	147	IF( nlat.NE.jjm ) THEN
	148	PRINT 2222
	149	STOP
	150	ELSE
[763]	151	c
[1086]	152	IF( iaire.EQ.1 ) THEN
	153	sdd1_type = type_sddu
	154	sdd2_type = type_unsddu
	155	ELSE
	156	sdd1_type = type_unsddu
	157	sdd2_type = type_sddu
	158	ENDIF
	159	c
	160	jdfil1 = 1
	161	jffil1 = jfiltnv
	162	jdfil2 = jfiltsv
	163	jffil2 = jjm
	164	ENDIF
	165	ENDIF
	166	c
	167	DO hemisph = 1, 2
	168	c
	169	IF ( hemisph.EQ.1 ) THEN
	170	cym
	171	jdfil = max(jdfil1,ibeg)
	172	jffil = min(jffil1,iend)
	173	ELSE
	174	cym
	175	jdfil = max(jdfil2,ibeg)
	176	jffil = min(jffil2,iend)
	177	ENDIF
[985]	178
	179
	180	cccccccccccccccccccccccccccccccccccccccccccc
	181	c Utilisation du filtre classique
	182	cccccccccccccccccccccccccccccccccccccccccccc
	183
[1086]	184	IF (.NOT. use_filtre_fft) THEN
[985]	185
[1086]	186	c !---------------------------------!
	187	c ! Agregation des niveau verticaux !
	188	c ! uniquement necessaire pour une !
	189	c ! execution OpenMP !
	190	c !---------------------------------!
	191	ll_nb = 0
	192	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	193	DO l = 1, nbniv
	194	ll_nb = ll_nb+1
	195	DO j = jdfil,jffil
	196	DO i = 1, iim
	197	champ_loc(i,j,ll_nb) =
	198	& champ(i,j,l) * sdd12(i,sdd1_type)
	199	ENDDO
	200	ENDDO
	201	ENDDO
	202	c$OMP END DO NOWAIT
[763]	203
[1086]	204	nbniv_loc = ll_nb
[763]	205
[1086]	206	IF( hemisph.EQ.1 ) THEN
	207
	208	IF( ifiltre.EQ.-2 ) THEN
	209	DO j = jdfil,jffil
	210	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	211	& matrinvn(1,1,j), iim,
	212	& champ_loc(1,j,1), iip1*nlat, 0.0,
	213	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
	214	ENDDO
	215
	216	ELSE IF ( griscal ) THEN
	217	DO j = jdfil,jffil
	218	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	219	& matriceun(1,1,j), iim,
	220	& champ_loc(1,j,1), iip1*nlat, 0.0,
	221	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
	222	ENDDO
	223
	224	ELSE
	225	DO j = jdfil,jffil
	226	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	227	& matricevn(1,1,j), iim,
	228	& champ_loc(1,j,1), iip1*nlat, 0.0,
	229	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
	230	ENDDO
	231
	232	ENDIF
	233
	234	ELSE
	235
	236	IF( ifiltre.EQ.-2 ) THEN
	237	DO j = jdfil,jffil
	238	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	239	& matrinvs(1,1,j-jfiltsu+1), iim,
	240	& champ_loc(1,j,1), iip1*nlat, 0.0,
	241	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
	242	ENDDO
	243
	244	ELSE IF ( griscal ) THEN
	245
	246	DO j = jdfil,jffil
	247	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	248	& matriceus(1,1,j-jfiltsu+1), iim,
	249	& champ_loc(1,j,1), iip1*nlat, 0.0,
	250	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
	251	ENDDO
	252
	253	ELSE
	254
	255	DO j = jdfil,jffil
	256	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	257	& matricevs(1,1,j-jfiltsv+1), iim,
	258	& champ_loc(1,j,1), iip1*nlat, 0.0,
	259	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
	260	ENDDO
	261
	262	ENDIF
	263
[985]	264	ENDIF
[1086]	265	! c
	266	IF( ifiltre.EQ.2 ) THEN
	267
	268	c !-------------------------------------!
	269	c ! Dés-agregation des niveau verticaux !
	270	c ! uniquement necessaire pour une !
	271	c ! execution OpenMP !
	272	c !-------------------------------------!
	273	ll_nb = 0
	274	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	275	DO l = 1, nbniv
	276	ll_nb = ll_nb + 1
	277	DO j = jdfil,jffil
	278	DO i = 1, iim
	279	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
	280	& + champ_fft(i,j-jdfil+1,ll_nb))
	281	& * sdd12(i,sdd2_type)
	282	ENDDO
	283	ENDDO
	284	ENDDO
	285	c$OMP END DO NOWAIT
	286
	287	ELSE
	288
	289	ll_nb = 0
	290	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	291	DO l = 1, nbniv_loc
	292	ll_nb = ll_nb + 1
	293	DO j = jdfil,jffil
	294	DO i = 1, iim
	295	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
	296	& - champ_fft(i,j-jdfil+1,ll_nb))
	297	& * sdd12(i,sdd2_type)
	298	ENDDO
	299	ENDDO
	300	ENDDO
	301	c$OMP END DO NOWAIT
	302
[985]	303	ENDIF
[1086]	304
	305	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	306	DO l = 1, nbniv
	307	DO j = jdfil,jffil
	308	champ( iip1,j,l ) = champ( 1,j,l )
	309	ENDDO
	310	ENDDO
[763]	311	c$OMP END DO NOWAIT
[1086]	312
[985]	313	ccccccccccccccccccccccccccccccccccccccccccccc
	314	c Utilisation du filtre FFT
	315	ccccccccccccccccccccccccccccccccccccccccccccc
	316
[1086]	317	ELSE
[985]	318
	319	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1086]	320	DO l=1,nbniv
	321	DO j=jdfil,jffil
	322	DO i = 1, iim
	323	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
	324	champ_fft( i,j,l) = champ(i,j,l)
	325	ENDDO
	326	ENDDO
[985]	327	ENDDO
	328	c$OMP END DO NOWAIT
	329
[1086]	330	IF (jdfil<=jffil) THEN
	331	IF( ifiltre. EQ. -2 ) THEN
	332	CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	333	ELSE IF ( griscal ) THEN
	334	CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	335	ELSE
	336	CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	337	ENDIF
	338	ENDIF
[985]	339
	340
[1086]	341	IF( ifiltre.EQ. 2 ) THEN
[985]	342	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1086]	343	DO l=1,nbniv
	344	DO j=jdfil,jffil
	345	DO i = 1, iim
	346	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
	347	& *sdd12(i,sdd2_type)
	348	ENDDO
	349	ENDDO
	350	ENDDO
[985]	351	c$OMP END DO NOWAIT
[1086]	352	ELSE
[985]	353
	354	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1086]	355	DO l=1,nbniv
	356	DO j=jdfil,jffil
	357	DO i = 1, iim
	358	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
	359	& *sdd12(i,sdd2_type)
	360	ENDDO
	361	ENDDO
	362	ENDDO
[985]	363	c$OMP END DO NOWAIT
[1086]	364	ENDIF
[985]	365	c
	366	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1086]	367	DO l=1,nbniv
	368	DO j=jdfil,jffil
[985]	369	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
[1086]	370	champ( iip1,j,l ) = champ( 1,j,l )
	371	ENDDO
	372	ENDDO
[985]	373	c$OMP END DO NOWAIT
[1086]	374	ENDIF
[985]	375	c Fin de la zone de filtrage
	376
	377
[1086]	378	ENDDO
[985]	379
	380	! DO j=1,nlat
	381	!
	382	! PRINT *,"check FFT ----> Delta(",j,")=",
	383	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
	384	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
	385	! ENDDO
	386
	387	! PRINT *,"check FFT ----> Delta(",j,")=",
	388	! & sum(champ-champ_fft)/sum(champ)
	389	!
	390
[763]	391	c
[1086]	392	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
	393	& filtrer, sur la grille des scalaires'/)
	394	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
	395	& ltrer, sur la grille de V ou de Z'/)
[985]	396	c$OMP MASTER
	397	CALL stop_timer
	398	c$OMP END MASTER
[763]	399	RETURN
	400	END

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

Download in other formats:

Original Format