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

