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 IMPLICIT NONE c======================================================================= c c Auteur: P. Le Van 07/10/97 c ------ c c Objet: filtre matriciel longitudinal ,avec les matrices precalculees c pour l'operateur Filtre . c ------ c c Arguments: c ---------- c c c ibeg..iend lattitude a filtrer c nlat nombre de latitudes du champ c nbniv nombre de niveaux verticaux a filtrer c champ(iip1,nblat,nbniv) en entree : champ a filtrer c en sortie : champ filtre c ifiltre +1 Transformee directe c -1 Transformee inverse c +2 Filtre directe c -2 Filtre inverse c c iaire 1 si champ intensif c 2 si champ extensif (pondere par les aires) c c iter 1 filtre simple c c======================================================================= c c c Variable Intensive c ifiltre = 1 filtre directe c ifiltre =-1 filtre inverse c c Variable Extensive c ifiltre = 2 filtre directe c ifiltre =-2 filtre inverse c c INCLUDE "dimensions.h" INCLUDE "paramet.h" INCLUDE "coefils.h" c 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. c$OMP THREADPRIVATE(first) REAL, DIMENSION(iip1,jjb:jje,nbniv) :: champ_loc INTEGER :: ll_nb, nbniv_loc REAL, SAVE :: sdd12(iim,4) c$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 c$OMP MASTER CALL start_timer c$OMP END MASTER c-------------------------------------------------------c IF(ifiltre==1.or.ifiltre==-1) & CALL abort_gcm("mod_filtreg_p",'Pas de transformee &simple dans cette version',1) IF( iter== 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== -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/=2 .AND.ifiltre/= - 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 c iim2 = iim * iim immjm = iim * jjm c c IF( griscal ) THEN IF( nlat. NE. jjp1 ) THEN CALL abort_gcm("mod_filtreg_p"," nlat. NE. jjp1",1) ELSE c IF( iaire==1 ) THEN sdd1_type = type_sddv sdd2_type = type_unsddv ELSE sdd1_type = type_unsddv sdd2_type = type_sddv ENDIF c jdfil1 = 2 jffil1 = jfiltnu jdfil2 = jfiltsu jffil2 = jjm ENDIF ELSE IF( nlat/=jjm ) THEN CALL abort_gcm("mod_filtreg_p"," nlat. NE. jjm",1) ELSE c IF( iaire==1 ) THEN sdd1_type = type_sddu sdd2_type = type_unsddu ELSE sdd1_type = type_unsddu sdd2_type = type_sddu ENDIF c jdfil1 = 1 jffil1 = jfiltnv jdfil2 = jfiltsv jffil2 = jjm ENDIF ENDIF c DO hemisph = 1, 2 c IF ( hemisph==1 ) THEN cym jdfil = max(jdfil1,ibeg) jffil = min(jffil1,iend) ELSE cym jdfil = max(jdfil2,ibeg) jffil = min(jffil2,iend) ENDIF cccccccccccccccccccccccccccccccccccccccccccc c Utilisation du filtre classique cccccccccccccccccccccccccccccccccccccccccccc IF (.NOT. use_filtre_fft) THEN c !---------------------------------! c ! Agregation des niveau verticaux ! c ! uniquement necessaire pour une ! c ! execution OpenMP ! c !---------------------------------! ll_nb = 0 c$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 c$OMP END DO NOWAIT nbniv_loc = ll_nb IF( hemisph==1 ) THEN IF( ifiltre==-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==-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==2 ) THEN c !-------------------------------------! c ! Dés-agregation des niveau verticaux ! c ! uniquement necessaire pour une ! c ! execution OpenMP ! c !-------------------------------------! ll_nb = 0 c$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 c$OMP END DO NOWAIT ELSE ll_nb = 0 c$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 c$OMP END DO NOWAIT ENDIF c$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 c$OMP END DO NOWAIT ccccccccccccccccccccccccccccccccccccccccccccc c Utilisation du filtre FFT ccccccccccccccccccccccccccccccccccccccccccccc ELSE c$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 c$OMP END DO NOWAIT IF (jdfil<=jffil) THEN IF( ifiltre == -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== 2 ) THEN c$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 c$OMP END DO NOWAIT ELSE c$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 c$OMP END DO NOWAIT ENDIF c c$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 c$OMP END DO NOWAIT ENDIF c 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) c 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'/) c$OMP MASTER CALL stop_timer c$OMP END MASTER RETURN END SUBROUTINE filtreg_p END MODULE mod_filtreg_p