Context Navigation

mod_filtreg_p.F90 @ 5456

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

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

Download in other formats:

Original Format