Context Navigation

filtreg_p.F @ 5097

Last change on this file since 5097 was 5091, checked in by abarral, 4 months ago
Move lmdz_netcdf_format.F90 -> lmdz_cppkeys_wrapper.F90 to handle other CPP keys Replace all (except wrapper) use of CPP_PHYS by fortran logical Refactor makelmdz_fcm (put blocks into functions, use modern bash)
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.0 KB

Rev	Line
[1632]	1
	2
	3	SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
	4	& ifiltre, iaire, griscal ,iter)
[1823]	5	USE parallel_lmdz, only : OMP_CHUNK
[1632]	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
[4593]	53	INCLUDE "dimensions.h"
	54	INCLUDE "paramet.h"
	55	INCLUDE "coefils.h"
[1632]	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
[5082]	101	IF(ifiltre==1.or.ifiltre==-1)
[4469]	102	& CALL abort_gcm("fitreg_p","Pas de transformee simple
	103	&dans cette version",1)
[1632]	104
[5082]	105	IF( iter== 2 ) THEN
[1632]	106	PRINT *,' Pas d iteration du filtre dans cette version !'
	107	& , ' Utiliser old_filtreg et repasser !'
[4469]	108	CALL abort_gcm("fitreg_p","stopped",1)
[1632]	109	ENDIF
	110
[5082]	111	IF( ifiltre== -2 .AND..NOT.griscal ) THEN
[1632]	112	PRINT *,' Cette routine ne calcule le filtre inverse que '
	113	& , ' sur la grille des scalaires !'
[4469]	114	CALL abort_gcm("fitreg_p","stopped",1)
[1632]	115	ENDIF
	116
[5082]	117	IF( ifiltre/=2 .AND.ifiltre/= - 2 ) THEN
[1632]	118	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
	119	& , ' corriger et repasser !'
[4469]	120	CALL abort_gcm("fitreg_p","stopped",1)
	121	ENDIF
[1632]	122	c
	123
	124	iim2 = iim * iim
	125	immjm = iim * jjm
	126	c
	127	c
	128	IF( griscal ) THEN
	129	IF( nlat. NE. jjp1 ) THEN
[4469]	130	CALL abort_gcm("fitreg_p","nlat. NE. jjp1",1)
[1632]	131	ELSE
	132	c
[5082]	133	IF( iaire==1 ) THEN
[1632]	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
[5082]	147	IF( nlat/=jjm ) THEN
[4469]	148	CALL abort_gcm("fitreg_p","nlat. NE. jjm",1)
[1632]	149	ELSE
	150	c
[5082]	151	IF( iaire==1 ) THEN
[1632]	152	sdd1_type = type_sddu
	153	sdd2_type = type_unsddu
	154	ELSE
	155	sdd1_type = type_unsddu
	156	sdd2_type = type_sddu
	157	ENDIF
	158	c
	159	jdfil1 = 1
	160	jffil1 = jfiltnv
	161	jdfil2 = jfiltsv
	162	jffil2 = jjm
	163	ENDIF
	164	ENDIF
	165	c
	166	DO hemisph = 1, 2
	167	c
[5082]	168	IF ( hemisph==1 ) THEN
[1632]	169	cym
	170	jdfil = max(jdfil1,ibeg)
	171	jffil = min(jffil1,iend)
	172	ELSE
	173	cym
	174	jdfil = max(jdfil2,ibeg)
	175	jffil = min(jffil2,iend)
	176	ENDIF
	177
	178
	179	cccccccccccccccccccccccccccccccccccccccccccc
	180	c Utilisation du filtre classique
	181	cccccccccccccccccccccccccccccccccccccccccccc
	182
	183	IF (.NOT. use_filtre_fft) THEN
	184
	185	c !---------------------------------!
	186	c ! Agregation des niveau verticaux !
	187	c ! uniquement necessaire pour une !
	188	c ! execution OpenMP !
	189	c !---------------------------------!
	190	ll_nb = 0
	191	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	192	DO l = 1, nbniv
	193	ll_nb = ll_nb+1
	194	DO j = jdfil,jffil
	195	DO i = 1, iim
	196	champ_loc(i,j,ll_nb) =
	197	& champ(i,j,l) * sdd12(i,sdd1_type)
	198	ENDDO
	199	ENDDO
	200	ENDDO
	201	c$OMP END DO NOWAIT
	202
	203	nbniv_loc = ll_nb
	204
[5082]	205	IF( hemisph==1 ) THEN
[1632]	206
[5082]	207	IF( ifiltre==-2 ) THEN
[1632]	208	DO j = jdfil,jffil
[1673]	209	#ifdef BLAS
[1632]	210	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	211	& matrinvn(1,1,j), iim,
	212	& champ_loc(1,j,1), iip1*nlat, 0.0,
	213	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1673]	214	#else
[1680]	215	champ_fft(:iim,j-jdfil+1,:)
[1673]	216	& =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
	217	#endif
[1632]	218	ENDDO
	219
	220	ELSE IF ( griscal ) THEN
	221	DO j = jdfil,jffil
[1673]	222	#ifdef BLAS
[1632]	223	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	224	& matriceun(1,1,j), iim,
	225	& champ_loc(1,j,1), iip1*nlat, 0.0,
	226	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1673]	227	#else
[1680]	228	champ_fft(:iim,j-jdfil+1,:)
[1673]	229	& =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
	230	#endif
[1632]	231	ENDDO
	232
	233	ELSE
	234	DO j = jdfil,jffil
[1673]	235	#ifdef BLAS
[1632]	236	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	237	& matricevn(1,1,j), iim,
	238	& champ_loc(1,j,1), iip1*nlat, 0.0,
	239	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1673]	240	#else
[1680]	241	champ_fft(:iim,j-jdfil+1,:)
[1673]	242	& =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
	243	#endif
[1632]	244	ENDDO
	245
	246	ENDIF
	247
	248	ELSE
	249
[5082]	250	IF( ifiltre==-2 ) THEN
[1632]	251	DO j = jdfil,jffil
[1673]	252	#ifdef BLAS
[1632]	253	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	254	& matrinvs(1,1,j-jfiltsu+1), iim,
	255	& champ_loc(1,j,1), iip1*nlat, 0.0,
	256	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1673]	257	#else
[1680]	258	champ_fft(:iim,j-jdfil+1,:)
[1673]	259	& =matmul(matrinvs(:,:,j-jfiltsu+1),
	260	& champ_loc(:iim,j,:))
	261	#endif
[1632]	262	ENDDO
	263
	264	ELSE IF ( griscal ) THEN
	265
	266	DO j = jdfil,jffil
[1673]	267	#ifdef BLAS
[1632]	268	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	269	& matriceus(1,1,j-jfiltsu+1), iim,
	270	& champ_loc(1,j,1), iip1*nlat, 0.0,
	271	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1673]	272	#else
[1680]	273	champ_fft(:iim,j-jdfil+1,:)
[1673]	274	& =matmul(matriceus(:,:,j-jfiltsu+1),
	275	& champ_loc(:iim,j,:))
	276	#endif
[1632]	277	ENDDO
	278
	279	ELSE
	280
	281	DO j = jdfil,jffil
[1673]	282	#ifdef BLAS
[1632]	283	CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0,
	284	& matricevs(1,1,j-jfiltsv+1), iim,
	285	& champ_loc(1,j,1), iip1*nlat, 0.0,
	286	& champ_fft(1,j-jdfil+1,1), iip1*nlat)
[1673]	287	#else
[1680]	288	champ_fft(:iim,j-jdfil+1,:)
[1673]	289	& =matmul(matricevs(:,:,j-jfiltsv+1),
	290	& champ_loc(:iim,j,:))
	291	#endif
[1632]	292	ENDDO
	293
	294	ENDIF
	295
	296	ENDIF
	297	! c
[5082]	298	IF( ifiltre==2 ) THEN
[1632]	299
	300	c !-------------------------------------!
	301	c ! Dés-agregation des niveau verticaux !
	302	c ! uniquement necessaire pour une !
	303	c ! execution OpenMP !
	304	c !-------------------------------------!
	305	ll_nb = 0
	306	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	307	DO l = 1, nbniv
	308	ll_nb = ll_nb + 1
	309	DO j = jdfil,jffil
	310	DO i = 1, iim
	311	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
	312	& + champ_fft(i,j-jdfil+1,ll_nb))
	313	& * sdd12(i,sdd2_type)
	314	ENDDO
	315	ENDDO
	316	ENDDO
	317	c$OMP END DO NOWAIT
	318
	319	ELSE
	320
	321	ll_nb = 0
	322	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	323	DO l = 1, nbniv_loc
	324	ll_nb = ll_nb + 1
	325	DO j = jdfil,jffil
	326	DO i = 1, iim
	327	champ( i,j,l ) = (champ_loc(i,j,ll_nb)
	328	& - champ_fft(i,j-jdfil+1,ll_nb))
	329	& * sdd12(i,sdd2_type)
	330	ENDDO
	331	ENDDO
	332	ENDDO
	333	c$OMP END DO NOWAIT
	334
	335	ENDIF
	336
	337	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	338	DO l = 1, nbniv
	339	DO j = jdfil,jffil
	340	champ( iip1,j,l ) = champ( 1,j,l )
	341	ENDDO
	342	ENDDO
	343	c$OMP END DO NOWAIT
	344
	345	ccccccccccccccccccccccccccccccccccccccccccccc
	346	c Utilisation du filtre FFT
	347	ccccccccccccccccccccccccccccccccccccccccccccc
	348
	349	ELSE
	350
	351	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	352	DO l=1,nbniv
	353	DO j=jdfil,jffil
	354	DO i = 1, iim
	355	champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type)
	356	champ_fft( i,j,l) = champ(i,j,l)
	357	ENDDO
	358	ENDDO
	359	ENDDO
	360	c$OMP END DO NOWAIT
	361
	362	IF (jdfil<=jffil) THEN
[5091]	363	IF( ifiltre == -2 ) THEN
[1632]	364	CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	365	ELSE IF ( griscal ) THEN
	366	CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	367	ELSE
	368	CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
	369	ENDIF
	370	ENDIF
	371
	372
[5082]	373	IF( ifiltre== 2 ) THEN
[1632]	374	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	375	DO l=1,nbniv
	376	DO j=jdfil,jffil
	377	DO i = 1, iim
	378	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
	379	& *sdd12(i,sdd2_type)
	380	ENDDO
	381	ENDDO
	382	ENDDO
	383	c$OMP END DO NOWAIT
	384	ELSE
	385
	386	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	387	DO l=1,nbniv
	388	DO j=jdfil,jffil
	389	DO i = 1, iim
	390	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
	391	& *sdd12(i,sdd2_type)
	392	ENDDO
	393	ENDDO
	394	ENDDO
	395	c$OMP END DO NOWAIT
	396	ENDIF
	397	c
	398	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
	399	DO l=1,nbniv
	400	DO j=jdfil,jffil
	401	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
	402	champ( iip1,j,l ) = champ( 1,j,l )
	403	ENDDO
	404	ENDDO
	405	c$OMP END DO NOWAIT
	406	ENDIF
	407	c Fin de la zone de filtrage
	408
	409
	410	ENDDO
	411
	412	! DO j=1,nlat
	413	!
	414	! PRINT *,"check FFT ----> Delta(",j,")=",
	415	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
	416	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
	417	! ENDDO
	418
	419	! PRINT *,"check FFT ----> Delta(",j,")=",
	420	! & sum(champ-champ_fft)/sum(champ)
	421	!
	422
	423	c
	424	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
	425	& filtrer, sur la grille des scalaires'/)
	426	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
	427	& ltrer, sur la grille de V ou de Z'/)
	428	c$OMP MASTER
	429	CALL stop_timer
	430	c$OMP END MASTER
	431	RETURN
	432	END

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

Download in other formats:

Original Format