MODULE mod_filtreg_p CONTAINS SUBROUTINE filtreg_p ( champ,jjb,jje, ibeg, iend, nlat, nbniv, & ifiltre, iaire, griscal ,iter) USE parallel_lmdz, only : OMP_CHUNK USE mod_filtre_fft_loc, ONLY: use_filtre_fft, filtre_u_fft, & filtre_v_fft, filtre_inv_fft USE timer_filtre, ONLY: init_timer, start_timer, stop_timer USE filtreg_mod, ONLY: matrinvn, matrinvs, matriceun, matriceus, & matricevn, matricevs USE dimensions_mod, ONLY: iim, jjm, llm, ndm USE paramet_mod_h USE coefils_mod_h IMPLICIT NONE !======================================================================= ! ! Auteur: P. Le Van 07/10/97 ! ------ ! ! Objet: filtre matriciel longitudinal ,avec les matrices precalculees ! pour l'operateur Filtre . ! ------ ! ! Arguments: ! ---------- ! ! ! ibeg..iend lattitude a filtrer ! nlat nombre de latitudes du champ ! nbniv nombre de niveaux verticaux a filtrer ! champ(iip1,nblat,nbniv) en entree : champ a filtrer ! en sortie : champ filtre ! ifiltre +1 Transformee directe ! -1 Transformee inverse ! +2 Filtre directe ! -2 Filtre inverse ! ! iaire 1 si champ intensif ! 2 si champ extensif (pondere par les aires) ! ! iter 1 filtre simple ! !======================================================================= ! ! ! Variable Intensive ! ifiltre = 1 filtre directe ! ifiltre =-1 filtre inverse ! ! Variable Extensive ! ifiltre = 2 filtre directe ! ifiltre =-2 filtre inverse ! ! INTEGER,INTENT(IN) :: jjb,jje,ibeg,iend,nlat,nbniv,ifiltre,iter INTEGER,INTENT(IN) :: iaire LOGICAL,INTENT(IN) :: griscal REAL,INTENT(INOUT) :: champ( iip1,jjb:jje,nbniv) INTEGER :: i,j,l,k INTEGER :: iim2,immjm INTEGER :: jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil INTEGER :: hemisph REAL :: champ_fft(iip1,jjb:jje,nbniv) ! REAL :: champ_in(iip1,jjb:jje,nbniv) LOGICAL,SAVE :: first=.TRUE. !$OMP THREADPRIVATE(first) REAL, DIMENSION(iip1,jjb:jje,nbniv) :: champ_loc INTEGER :: ll_nb, nbniv_loc REAL, SAVE :: sdd12(iim,4) !$OMP THREADPRIVATE(sdd12) INTEGER, PARAMETER :: type_sddu=1 INTEGER, PARAMETER :: type_sddv=2 INTEGER, PARAMETER :: type_unsddu=3 INTEGER, PARAMETER :: type_unsddv=4 INTEGER :: sdd1_type, sdd2_type CHARACTER (LEN=132) :: abort_message IF (first) THEN sdd12(1:iim,type_sddu) = sddu(1:iim) sdd12(1:iim,type_sddv) = sddv(1:iim) sdd12(1:iim,type_unsddu) = unsddu(1:iim) sdd12(1:iim,type_unsddv) = unsddv(1:iim) CALL Init_timer first=.FALSE. ENDIF !$OMP MASTER CALL start_timer !$OMP END MASTER !-------------------------------------------------------c IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) & CALL abort_gcm("mod_filtreg_p",'Pas de transformee& &simple dans cette version',1) IF( iter.EQ. 2 ) THEN PRINT *,' Pas d iteration du filtre dans cette version !'& & , ' Utiliser old_filtreg et repasser !' CALL abort_gcm("mod_filtreg_p","stopped",1) ENDIF IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN PRINT *,' Cette routine ne calcule le filtre inverse que ' & , ' sur la grille des scalaires !' CALL abort_gcm("mod_filtreg_p","stopped",1) ENDIF IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' & , ' corriger et repasser !' CALL abort_gcm("mod_filtreg_p","stopped",1) ENDIF ! iim2 = iim * iim immjm = iim * jjm ! ! IF( griscal ) THEN IF( nlat.NE. jjp1 ) THEN CALL abort_gcm("mod_filtreg_p"," nlat.NE. jjp1",1) ELSE ! IF( iaire.EQ.1 ) THEN sdd1_type = type_sddv sdd2_type = type_unsddv ELSE sdd1_type = type_unsddv sdd2_type = type_sddv ENDIF ! jdfil1 = 2 jffil1 = jfiltnu jdfil2 = jfiltsu jffil2 = jjm ENDIF ELSE IF( nlat.NE.jjm ) THEN CALL abort_gcm("mod_filtreg_p"," nlat.NE. jjm",1) ELSE ! IF( iaire.EQ.1 ) THEN sdd1_type = type_sddu sdd2_type = type_unsddu ELSE sdd1_type = type_unsddu sdd2_type = type_sddu ENDIF ! jdfil1 = 1 jffil1 = jfiltnv jdfil2 = jfiltsv jffil2 = jjm ENDIF ENDIF ! DO hemisph = 1, 2 ! IF ( hemisph.EQ.1 ) THEN !ym jdfil = max(jdfil1,ibeg) jffil = min(jffil1,iend) ELSE !ym jdfil = max(jdfil2,ibeg) jffil = min(jffil2,iend) ENDIF !ccccccccccccccccccccccccccccccccccccccccccc ! Utilisation du filtre classique !ccccccccccccccccccccccccccccccccccccccccccc IF (.NOT. use_filtre_fft) THEN ! !---------------------------------! ! ! Agregation des niveau verticaux ! ! ! uniquement necessaire pour une ! ! ! execution OpenMP ! ! !---------------------------------! ll_nb = 0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv ll_nb = ll_nb+1 DO j = jdfil,jffil DO i = 1, iim champ_loc(i,j,ll_nb) = & champ(i,j,l) * sdd12(i,sdd1_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT nbniv_loc = ll_nb IF( hemisph.EQ.1 ) THEN IF( ifiltre.EQ.-2 ) THEN DO j = jdfil,jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matrinvn(1,1,j), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim,j,1:nbniv_loc)= & matmul(matrinvn(1:iim,1:iim,j), & champ_loc(1:iim,j,1:nbniv_loc)) #endif ENDDO ELSE IF ( griscal ) THEN DO j = jdfil,jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matriceun(1,1,j), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim,j,1:nbniv_loc)= & matmul(matriceun(1:iim,1:iim,j), & champ_loc(1:iim,j,1:nbniv_loc)) #endif ENDDO ELSE DO j = jdfil,jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matricevn(1,1,j), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim,j,1:nbniv_loc)= & matmul(matricevn(1:iim,1:iim,j), & champ_loc(1:iim,j,1:nbniv_loc)) #endif ENDDO ENDIF ELSE IF( ifiltre.EQ.-2 ) THEN DO j = jdfil,jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matrinvs(1,1,j-jfiltsu+1), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim,j,1:nbniv_loc)= & matmul(matrinvs(1:iim,1:iim,j-jfiltsu+1), & champ_loc(1:iim,j,1:nbniv_loc)) #endif ENDDO ELSE IF ( griscal ) THEN DO j = jdfil,jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matriceus(1,1,j-jfiltsu+1), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim,j,1:nbniv_loc)= & matmul(matriceus(1:iim,1:iim,j-jfiltsu+1), & champ_loc(1:iim,j,1:nbniv_loc)) #endif ENDDO ELSE DO j = jdfil,jffil #ifdef BLAS CALL SGEMM("N", "N", iim, nbniv_loc, iim, 1.0, & matricevs(1,1,j-jfiltsv+1), iim, & champ_loc(1,j,1), iip1*(jje-jjb+1), 0.0, & champ_fft(1,j,1), iip1*(jje-jjb+1)) #else champ_fft(1:iim,j,1:nbniv_loc)= & matmul(matricevs(1:iim,1:iim,j-jfiltsv+1), & champ_loc(1:iim,j,1:nbniv_loc)) #endif ENDDO ENDIF ENDIF ! c IF( ifiltre.EQ.2 ) THEN ! !-------------------------------------! ! ! Dés-agregation des niveau verticaux ! ! ! uniquement necessaire pour une ! ! ! execution OpenMP ! ! !-------------------------------------! ll_nb = 0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv ll_nb = ll_nb + 1 DO j = jdfil,jffil DO i = 1, iim champ( i,j,l ) = (champ_loc(i,j,ll_nb) & + champ_fft(i,j,ll_nb)) & * sdd12(i,sdd2_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT ELSE ll_nb = 0 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv ll_nb = ll_nb + 1 DO j = jdfil,jffil DO i = 1, iim champ( i,j,l ) = (champ_loc(i,j,ll_nb) & - champ_fft(i,j,ll_nb)) & * sdd12(i,sdd2_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT ENDIF !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l = 1, nbniv DO j = jdfil,jffil ! ! add redundant longitude champ( iip1,j,l ) = champ( 1,j,l ) ENDDO ENDDO !$OMP END DO NOWAIT !cccccccccccccccccccccccccccccccccccccccccccc ! Utilisation du filtre FFT !cccccccccccccccccccccccccccccccccccccccccccc ELSE !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,nbniv DO j=jdfil,jffil DO i = 1, iim champ( i,j,l)= champ(i,j,l)*sdd12(i,sdd1_type) champ_fft( i,j,l) = champ(i,j,l) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT IF (jdfil<=jffil) THEN IF( ifiltre.EQ. -2 ) THEN CALL Filtre_inv_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv) ELSE IF ( griscal ) THEN CALL Filtre_u_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv) ELSE CALL Filtre_v_fft(champ_fft,jjb,jje,jdfil,jffil,nbniv) ENDIF ENDIF IF( ifiltre.EQ. 2 ) THEN !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,nbniv DO j=jdfil,jffil DO i = 1, iim champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l)) & *sdd12(i,sdd2_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT ELSE !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,nbniv DO j=jdfil,jffil DO i = 1, iim champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l)) & *sdd12(i,sdd2_type) ENDDO ENDDO ENDDO !$OMP END DO NOWAIT ENDIF ! !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO l=1,nbniv DO j=jdfil,jffil ! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l ) ! ! add redundant longitude champ( iip1,j,l ) = champ( 1,j,l ) ENDDO ENDDO !$OMP END DO NOWAIT ENDIF ! Fin de la zone de filtrage ENDDO ! DO j=1,nlat ! ! PRINT *,"check FFT ----> Delta(",j,")=", ! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)), ! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:)) ! ENDDO ! PRINT *,"check FFT ----> Delta(",j,")=", ! & sum(champ-champ_fft)/sum(champ) ! ! 1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a& & filtrer, sur la grille des scalaires'/) 2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi& & ltrer, sur la grille de V ou de Z'/) !$OMP MASTER CALL stop_timer !$OMP END MASTER RETURN END SUBROUTINE filtreg_p END MODULE mod_filtreg_p