Context Navigation

filtreg_p.F @ 1793

Last change on this file since 1793 was 1680, checked in by Ehouarn Millour, 12 years ago

bug correction in filtreg_p.F on array boundaries (was a problem only in cases where BLAS was not used).
modified "makelmdz_fcm" script: added '-mem' option to trigger using "dyn3dmem" dynamics instead of (default) "dyn3dpar".

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
[1597]	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)
[1597]	215	#else
[1680]	216	champ_fft(:iim,j-jdfil+1,:)
[1597]	217	& =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
	218	#endif
[1146]	219	ENDDO
	220
	221	ELSE IF ( griscal ) THEN
	222	DO j = jdfil,jffil
[1597]	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)
[1597]	228	#else
[1680]	229	champ_fft(:iim,j-jdfil+1,:)
[1597]	230	& =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
	231	#endif
[1146]	232	ENDDO
	233
	234	ELSE
	235	DO j = jdfil,jffil
[1597]	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)
[1597]	241	#else
[1680]	242	champ_fft(:iim,j-jdfil+1,:)
[1597]	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
[1597]	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)
[1597]	258	#else
[1680]	259	champ_fft(:iim,j-jdfil+1,:)
[1597]	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
[1597]	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)
[1597]	273	#else
[1680]	274	champ_fft(:iim,j-jdfil+1,:)
[1597]	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
[1597]	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)
[1597]	288	#else
[1680]	289	champ_fft(:iim,j-jdfil+1,:)
[1597]	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