Context Navigation

filtreg_p.F @ 1102

Last change on this file since 1102 was 985, checked in by Laurent Fairhead, 16 years ago
Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 10.0 KB

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

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

Download in other formats:

Original Format