Context Navigation

filtreg_p.F @ 1601

Last change on this file since 1601 was 1601, checked in by Ehouarn Millour, 13 years ago

Updates and upgrades of filter:

add minor corrections (r1591 of trunk) in filtreg_mod.F90: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong in extreme cases.
adapt filtreg_p.F (r1597 of trunk) so that BLAS routine DGEMM is only used if 'BLAS' preprocessing flag is set.
update FFT filter routines to match current 'trunk' versions (mostly cosmetic changes, exept for the implementation of use of FFTW in mod_fft_fftw.F90).
update "arch-PW6_VARGAS.fcm" to enable use of FFTW

NB: implemented FFTs assume working precision to be double precision ; trying to use them when working precision is single precision will clearly end in tragedy.
EM

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

File size: 12.9 KB

Rev	Line
[985]	1
	2
[763]	3	SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
[1146]	4	& ifiltre, iaire, griscal ,iter)
	5	USE Parallel, only : OMP_CHUNK
[985]	6	USE mod_filtre_fft
	7	USE timer_filtre
[1146]	8
	9	USE filtreg_mod
	10
[763]	11	IMPLICIT NONE
[1146]	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
[1146]	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
[1146]	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
[1146]	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
[1146]	99	c-------------------------------------------------------c
	100
[763]	101	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
[1146]	102	& STOP'Pas de transformee simple dans cette version'
	103
[763]	104	IF( iter.EQ. 2 ) THEN
[1146]	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
[1146]	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
[1146]	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
[1146]	129	PRINT 1111
	130	STOP
[763]	131	ELSE
[1146]	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
[1146]	141	jdfil1 = 2
	142	jffil1 = jfiltnu
	143	jdfil2 = jfiltsu
	144	jffil2 = jjm
	145	ENDIF
[763]	146	ELSE
[1146]	147	IF( nlat.NE.jjm ) THEN
	148	PRINT 2222
	149	STOP
	150	ELSE
[763]	151	c
[1146]	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
[1146]	184	IF (.NOT. use_filtre_fft) THEN
[985]	185
[1146]	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
[1146]	204	nbniv_loc = ll_nb
[763]	205
[1146]	206	IF( hemisph.EQ.1 ) THEN
	207
	208	IF( ifiltre.EQ.-2 ) THEN
	209	DO j = jdfil,jffil
[1601]	210	#ifdef BLAS
[1146]	211	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	212	& matrinvn(1,1,j), iim,
	213	& champ_loc(1,j,1), iip1*nlat, 0.0,
	214	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1601]	215	#else
	216	champ_fft(:,j-jdfil+1,:)
	217	& =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
	218	#endif
[1146]	219	ENDDO
	220
	221	ELSE IF ( griscal ) THEN
	222	DO j = jdfil,jffil
[1601]	223	#ifdef BLAS
[1146]	224	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	225	& matriceun(1,1,j), iim,
	226	& champ_loc(1,j,1), iip1*nlat, 0.0,
	227	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1601]	228	#else
	229	champ_fft(:,j-jdfil+1,:)
	230	& =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
	231	#endif
[1146]	232	ENDDO
	233
	234	ELSE
	235	DO j = jdfil,jffil
[1601]	236	#ifdef BLAS
[1146]	237	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	238	& matricevn(1,1,j), iim,
	239	& champ_loc(1,j,1), iip1*nlat, 0.0,
	240	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1601]	241	#else
	242	champ_fft(:,j-jdfil+1,:)
	243	& =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
	244	#endif
[1146]	245	ENDDO
	246
	247	ENDIF
	248
	249	ELSE
	250
	251	IF( ifiltre.EQ.-2 ) THEN
	252	DO j = jdfil,jffil
[1601]	253	#ifdef BLAS
[1146]	254	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	255	& matrinvs(1,1,j-jfiltsu+1), iim,
	256	& champ_loc(1,j,1), iip1*nlat, 0.0,
	257	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1601]	258	#else
	259	champ_fft(:,j-jdfil+1,:)
	260	& =matmul(matrinvs(:,:,j-jfiltsu+1),
	261	& champ_loc(:iim,j,:))
	262	#endif
[1146]	263	ENDDO
	264
	265	ELSE IF ( griscal ) THEN
	266
	267	DO j = jdfil,jffil
[1601]	268	#ifdef BLAS
[1146]	269	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	270	& matriceus(1,1,j-jfiltsu+1), iim,
	271	& champ_loc(1,j,1), iip1*nlat, 0.0,
	272	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1601]	273	#else
	274	champ_fft(:,j-jdfil+1,:)
	275	& =matmul(matriceus(:,:,j-jfiltsu+1),
	276	& champ_loc(:iim,j,:))
	277	#endif
[1146]	278	ENDDO
	279
	280	ELSE
	281
	282	DO j = jdfil,jffil
[1601]	283	#ifdef BLAS
[1146]	284	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	285	& matricevs(1,1,j-jfiltsv+1), iim,
	286	& champ_loc(1,j,1), iip1*nlat, 0.0,
	287	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1601]	288	#else
	289	champ_fft(:,j-jdfil+1,:)
	290	& =matmul(matricevs(:,:,j-jfiltsv+1),
	291	& champ_loc(:iim,j,:))
	292	#endif
[1146]	293	ENDDO
	294
	295	ENDIF
	296
[985]	297	ENDIF
[1146]	298	! c
	299	IF( ifiltre.EQ.2 ) THEN
	300
	301	c !-------------------------------------!
	302	c ! Dés-agregation des niveau verticaux !
	303	c ! uniquement necessaire pour une !
	304	c ! execution OpenMP !
	305	c !-------------------------------------!
	306	ll_nb = 0
	307	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	308	DO l = 1, nbniv
	309	ll_nb = ll_nb + 1
	310	DO j = jdfil,jffil
	311	DO i = 1, iim
	312	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
	313	& + champ_fft(i,j-jdfil+1,ll_nb))
	314	& * sdd12(i,sdd2_type)
	315	ENDDO
	316	ENDDO
	317	ENDDO
	318	c$OMP END DO NOWAIT
	319
	320	ELSE
	321
	322	ll_nb = 0
	323	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	324	DO l = 1, nbniv_loc
	325	ll_nb = ll_nb + 1
	326	DO j = jdfil,jffil
	327	DO i = 1, iim
	328	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
	329	& - champ_fft(i,j-jdfil+1,ll_nb))
	330	& * sdd12(i,sdd2_type)
	331	ENDDO
	332	ENDDO
	333	ENDDO
	334	c$OMP END DO NOWAIT
	335
[985]	336	ENDIF
[1146]	337
	338	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	339	DO l = 1, nbniv
	340	DO j = jdfil,jffil
	341	champ( iip1,j,l ) = champ( 1,j,l )
	342	ENDDO
	343	ENDDO
[763]	344	c$OMP END DO NOWAIT
[1146]	345
[985]	346	ccccccccccccccccccccccccccccccccccccccccccccc
	347	c Utilisation du filtre FFT
	348	ccccccccccccccccccccccccccccccccccccccccccccc
	349
[1146]	350	ELSE
[985]	351
	352	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1146]	353	DO l=1,nbniv
	354	DO j=jdfil,jffil
	355	DO i = 1, iim
	356	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
	357	champ_fft( i,j,l) = champ(i,j,l)
	358	ENDDO
	359	ENDDO
[985]	360	ENDDO
	361	c$OMP END DO NOWAIT
	362
[1146]	363	IF (jdfil<=jffil) THEN
	364	IF( ifiltre. EQ. -2 ) THEN
	365	CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	366	ELSE IF ( griscal ) THEN
	367	CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	368	ELSE
	369	CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	370	ENDIF
	371	ENDIF
[985]	372
	373
[1146]	374	IF( ifiltre.EQ. 2 ) THEN
[985]	375	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1146]	376	DO l=1,nbniv
	377	DO j=jdfil,jffil
	378	DO i = 1, iim
	379	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
	380	& *sdd12(i,sdd2_type)
	381	ENDDO
	382	ENDDO
	383	ENDDO
[985]	384	c$OMP END DO NOWAIT
[1146]	385	ELSE
[985]	386
	387	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1146]	388	DO l=1,nbniv
	389	DO j=jdfil,jffil
	390	DO i = 1, iim
	391	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
	392	& *sdd12(i,sdd2_type)
	393	ENDDO
	394	ENDDO
	395	ENDDO
[985]	396	c$OMP END DO NOWAIT
[1146]	397	ENDIF
[985]	398	c
	399	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
[1146]	400	DO l=1,nbniv
	401	DO j=jdfil,jffil
[985]	402	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
[1146]	403	champ( iip1,j,l ) = champ( 1,j,l )
	404	ENDDO
	405	ENDDO
[985]	406	c$OMP END DO NOWAIT
[1146]	407	ENDIF
[985]	408	c Fin de la zone de filtrage
	409
	410
[1146]	411	ENDDO
[985]	412
	413	! DO j=1,nlat
	414	!
	415	! PRINT *,"check FFT ----> Delta(",j,")=",
	416	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
	417	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
	418	! ENDDO
	419
	420	! PRINT *,"check FFT ----> Delta(",j,")=",
	421	! & sum(champ-champ_fft)/sum(champ)
	422	!
	423
[763]	424	c
[1146]	425	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
	426	& filtrer, sur la grille des scalaires'/)
	427	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
	428	& ltrer, sur la grille de V ou de Z'/)
[985]	429	c$OMP MASTER
	430	CALL stop_timer
	431	c$OMP END MASTER
[763]	432	RETURN
	433	END

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

Download in other formats:

Original Format