Context Navigation

mod_filtreg_p.F @ 2117

Last change on this file since 2117 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

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

Download in other formats:

Original Format