Context Navigation

filtreg_p.F @ 1807

Last change on this file since 1807 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".

File size: 12.9 KB

Rev	Line
[1632]	1
	2
	3	SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
	4	& ifiltre, iaire, griscal ,iter)
	5	USE Parallel, only : OMP_CHUNK
	6	USE mod_filtre_fft
	7	USE timer_filtre
	8
	9	USE filtreg_mod
	10
	11	IMPLICIT NONE
	12
	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
	61
	62	REAL champ( iip1,nlat,nbniv)
	63
	64	LOGICAL griscal
	65	INTEGER hemisph, iaire
	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)
	72
	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
	85	IF (first) THEN
	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.
	93	ENDIF
	94
	95	c$OMP MASTER
	96	CALL start_timer
	97	c$OMP END MASTER
	98
	99	c-------------------------------------------------------c
	100
	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
	134	sdd1_type = type_sddv
	135	sdd2_type = type_unsddv
	136	ELSE
	137	sdd1_type = type_unsddv
	138	sdd2_type = type_sddv
	139	ENDIF
	140	c
	141	jdfil1 = 2
	142	jffil1 = jfiltnu
	143	jdfil2 = jfiltsu
	144	jffil2 = jjm
	145	ENDIF
	146	ELSE
	147	IF( nlat.NE.jjm ) THEN
	148	PRINT 2222
	149	STOP
	150	ELSE
	151	c
	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
	178
	179
	180	cccccccccccccccccccccccccccccccccccccccccccc
	181	c Utilisation du filtre classique
	182	cccccccccccccccccccccccccccccccccccccccccccc
	183
	184	IF (.NOT. use_filtre_fft) THEN
	185
	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
	203
	204	nbniv_loc = ll_nb
	205
	206	IF( hemisph.EQ.1 ) THEN
	207
	208	IF( ifiltre.EQ.-2 ) THEN
	209	DO j = jdfil,jffil
[1673]	210	#ifdef BLAS
[1632]	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)
[1673]	215	#else
[1680]	216	champ_fft(:iim,j-jdfil+1,:)
[1673]	217	& =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
	218	#endif
[1632]	219	ENDDO
	220
	221	ELSE IF ( griscal ) THEN
	222	DO j = jdfil,jffil
[1673]	223	#ifdef BLAS
[1632]	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)
[1673]	228	#else
[1680]	229	champ_fft(:iim,j-jdfil+1,:)
[1673]	230	& =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
	231	#endif
[1632]	232	ENDDO
	233
	234	ELSE
	235	DO j = jdfil,jffil
[1673]	236	#ifdef BLAS
[1632]	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)
[1673]	241	#else
[1680]	242	champ_fft(:iim,j-jdfil+1,:)
[1673]	243	& =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
	244	#endif
[1632]	245	ENDDO
	246
	247	ENDIF
	248
	249	ELSE
	250
	251	IF( ifiltre.EQ.-2 ) THEN
	252	DO j = jdfil,jffil
[1673]	253	#ifdef BLAS
[1632]	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)
[1673]	258	#else
[1680]	259	champ_fft(:iim,j-jdfil+1,:)
[1673]	260	& =matmul(matrinvs(:,:,j-jfiltsu+1),
	261	& champ_loc(:iim,j,:))
	262	#endif
[1632]	263	ENDDO
	264
	265	ELSE IF ( griscal ) THEN
	266
	267	DO j = jdfil,jffil
[1673]	268	#ifdef BLAS
[1632]	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)
[1673]	273	#else
[1680]	274	champ_fft(:iim,j-jdfil+1,:)
[1673]	275	& =matmul(matriceus(:,:,j-jfiltsu+1),
	276	& champ_loc(:iim,j,:))
	277	#endif
[1632]	278	ENDDO
	279
	280	ELSE
	281
	282	DO j = jdfil,jffil
[1673]	283	#ifdef BLAS
[1632]	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)
[1673]	288	#else
[1680]	289	champ_fft(:iim,j-jdfil+1,:)
[1673]	290	& =matmul(matricevs(:,:,j-jfiltsv+1),
	291	& champ_loc(:iim,j,:))
	292	#endif
[1632]	293	ENDDO
	294
	295	ENDIF
	296
	297	ENDIF
	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
	336	ENDIF
	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
	344	c$OMP END DO NOWAIT
	345
	346	ccccccccccccccccccccccccccccccccccccccccccccc
	347	c Utilisation du filtre FFT
	348	ccccccccccccccccccccccccccccccccccccccccccccc
	349
	350	ELSE
	351
	352	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	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
	360	ENDDO
	361	c$OMP END DO NOWAIT
	362
	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
	372
	373
	374	IF( ifiltre.EQ. 2 ) THEN
	375	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	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
	384	c$OMP END DO NOWAIT
	385	ELSE
	386
	387	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	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
	396	c$OMP END DO NOWAIT
	397	ENDIF
	398	c
	399	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	400	DO l=1,nbniv
	401	DO j=jdfil,jffil
	402	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
	403	champ( iip1,j,l ) = champ( 1,j,l )
	404	ENDDO
	405	ENDDO
	406	c$OMP END DO NOWAIT
	407	ENDIF
	408	c Fin de la zone de filtrage
	409
	410
	411	ENDDO
	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
	424	c
	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'/)
	429	c$OMP MASTER
	430	CALL stop_timer
	431	c$OMP END MASTER
	432	RETURN
	433	END

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

Download in other formats:

Original Format