Context Navigation

filtreg_p.F @ 1673

Last change on this file since 1673 was 1673, checked in by Laurent Fairhead, 12 years ago

Fin du phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ5 r1671
Il reste quelques routines a verifier (en particulier ce qui touche a l'etude des cas academiques)
et la validation a effectuer

End of the phasing of the localised (low memory) parallel dynamics package with the
LMDZ5 trunk (r1671)
Some routines still need some checking (in particular the academic cases) and some
validation is still required

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
	216	champ_fft(:,j-jdfil+1,:)
	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
	229	champ_fft(:,j-jdfil+1,:)
	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
	242	champ_fft(:,j-jdfil+1,:)
	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
	259	champ_fft(:,j-jdfil+1,:)
	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
	274	champ_fft(:,j-jdfil+1,:)
	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
	289	champ_fft(:,j-jdfil+1,:)
	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