Context Navigation

mod_filtreg_p.F @ 5437

Last change on this file since 5437 was 2125, checked in by Ehouarn Millour, 10 years ago
Bug fix: wrong (implicit) bounds for matmul operations, in the non-FFT filteringcase. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 13.7 KB

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

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

Download in other formats:

Original Format