SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv, . ifiltre, iaire, griscal ,iter) USE Parallel, only : OMP_CHUNK USE mod_filtre_fft USE timer_filtre 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 "parafilt.h" #include "coefils.h" c INTEGER ibeg,iend,nlat,nbniv,ifiltre,iter INTEGER i,j,l,k INTEGER iim2,immjm INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil REAL champ( iip1,nlat,nbniv) REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus) , , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs) , , matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus) cym REAL eignq(iim), sdd1(iim),sdd2(iim) REAL eignq(iim) REAL :: sdd1(iim),sdd2(iim) LOGICAL griscal INTEGER hemisph, iaire REAL :: champ_fft(iip1,nlat,nbniv) REAL :: champ_in(iip1,nlat,nbniv) REAL,SAVE,TARGET :: sddu_loc(iim) REAL,SAVE,TARGET :: sddv_loc(iim) REAL,SAVE,TARGET :: unsddu_loc(iim) REAL,SAVE,TARGET :: unsddv_loc(iim) c$OMP THREADPRIVATE(sddu_loc,sddv_loc,unsddu_loc,unsddv_loc) LOGICAL,SAVE :: first=.TRUE. c$OMP THREADPRIVATE(first) IF (first) THEN sddu_loc(1:iim)=sddu(1:iim) sddv_loc(1:iim)=sddv(1:iim) unsddu_loc(1:iim)=unsddu(1:iim) unsddv_loc(1:iim)=unsddv(1:iim) CALL Init_timer first=.FALSE. c PRINT *,"----> sddu_loc=",sddu_loc c PRINT *,"----> sddv_loc=",sddv_loc c PRINT *,"----> unsddu_loc=",unsddu_loc c PRINT *,"----> unsddv_loc=",unsddv_loc ENDIF c$OMP MASTER CALL start_timer c$OMP END MASTER IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) * STOP'Pas de transformee simple dans cette version' IF( iter.EQ. 2 ) THEN PRINT *,' Pas d iteration du filtre dans cette version !' * , ' Utiliser old_filtreg et repasser !' STOP ENDIF IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN PRINT *,' Cette routine ne calcule le filtre inverse que ', * ' sur la grille des scalaires !' STOP ENDIF IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2' *,' corriger et repasser !' STOP ENDIF c iim2 = iim * iim immjm = iim * jjm c c IF( griscal ) THEN IF( nlat. NE. jjp1 ) THEN PRINT 1111 STOP ELSE c IF( iaire.EQ.1 ) THEN cym CALL SCOPY( iim, sddv, 1, sdd1, 1 ) cym CALL SCOPY( iim, unsddv, 1, sdd2, 1 ) cym sdd1=>sddv_loc cym sdd2=>unsddv_loc sdd1(1:iim)=sddv_loc(1:iim) sdd2(1:iim)=unsddv_loc(1:iim) ELSE cym CALL SCOPY( iim, unsddv, 1, sdd1, 1 ) cym CALL SCOPY( iim, sddv, 1, sdd2, 1 ) sdd1(1:iim)=unsddv_loc(1:iim) sdd2(1:iim)=sddv_loc(1:iim) END IF c jdfil1 = 2 jffil1 = jfiltnu jdfil2 = jfiltsu jffil2 = jjm END IF ELSE IF( nlat.NE.jjm ) THEN PRINT 2222 STOP ELSE c IF( iaire.EQ.1 ) THEN cym CALL SCOPY( iim, sddu, 1, sdd1, 1 ) cym CALL SCOPY( iim, unsddu, 1, sdd2, 1 ) cym sdd1=>sddu_loc cym sdd2=>unsddu_loc sdd1(1:iim)=sddu_loc(1:iim) sdd2(1:iim)=unsddu_loc(1:iim) ELSE cym CALL SCOPY( iim, unsddu, 1, sdd1, 1 ) cym CALL SCOPY( iim, sddu, 1, sdd2, 1 ) cym sdd1=>unsddu_loc cym sdd2=>sddu_loc sdd1(1:iim)=unsddu_loc(1:iim) sdd2(1:iim)=sddu_loc(1:iim) END IF c jdfil1 = 1 jffil1 = jfiltnv jdfil2 = jfiltsv jffil2 = jjm END IF END IF c PRINT *,"APPEL a filtreg --> sdd1=",sdd1 c PRINT *,"APPEL a filtreg --> sdd2=",sdd2 c PRINT *,"----> sddu_loc=",sddu_loc c PRINT *,"----> sddv_loc=",sddv_loc c PRINT *,"----> unsddu_loc=",unsddu_loc c PRINT *,"----> unsddv_loc=",unsddv_loc c c DO 100 hemisph = 1, 2 c IF ( hemisph.EQ.1 ) THEN c ym jdfil = max(jdfil1,ibeg) jffil = min(jffil1,iend) ELSE c ym jdfil = max(jdfil2,ibeg) jffil = min(jffil2,iend) END IF cccccccccccccccccccccccccccccccccccccccccccc c Utilisation du filtre classique cccccccccccccccccccccccccccccccccccccccccccc IF (.NOT. use_filtre_fft) THEN c$OMP DO SCHEDULE(STATIC,OMP_CHUNK) DO 50 l = 1, nbniv DO 30 j = jdfil,jffil DO 5 i = 1, iim champ(i,j,l) = champ(i,j,l) * sdd1(i) 5 CONTINUE c IF( hemisph. EQ. 1 ) THEN IF( ifiltre. EQ. -2 ) THEN CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim, . champ(1,j,l), 1, 0.0, eignq, 1) ELSE IF ( griscal ) THEN CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim, . champ(1,j,l), 1, 0.0, eignq, 1) ELSE CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim, . champ(1,j,l), 1, 0.0, eignq, 1) ENDIF ELSE IF( ifiltre. EQ. -2 ) THEN CALL SGEMV("N",iim,iim,1.0, matrinvs(1,1,j-jfiltsu+1),iim, . champ(1,j,l), 1, 0.0, eignq, 1) ELSE IF ( griscal ) THEN CALL SGEMV("N",iim,iim,1.0,matriceus(1,1,j-jfiltsu+1),iim, . champ(1,j,l), 1, 0.0, eignq, 1) ELSE CALL SGEMV("N",iim,iim,1.0,matricevs(1,1,j-jfiltsv+1),iim, . champ(1,j,l), 1, 0.0, eignq, 1) ENDIF ENDIF c IF( ifiltre.EQ. 2 ) THEN DO 15 i = 1, iim champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i) 15 CONTINUE ELSE DO 16 i=1,iim champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i) 16 CONTINUE ENDIF c champ( iip1,j,l ) = champ( 1,j,l ) c 30 CONTINUE c 50 CONTINUE 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)*sdd1(i) champ_fft( i,j,l) = champ(i,j,l) ENDDO ENDDO ENDDO c$OMP END DO NOWAIT IF (jdfil<=jffil) THEN IF( ifiltre. EQ. -2 ) THEN CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv) ELSE IF ( griscal ) THEN CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv) ELSE CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv) ENDIF ENDIF IF( ifiltre.EQ. 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)) . *sdd2(i) 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)) . *sdd2(i) 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 ) champ( iip1,j,l ) = champ( 1,j,l ) ENDDO ENDDO c$OMP END DO NOWAIT ENDIF c Fin de la zone de filtrage 100 CONTINUE ! 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