! Amaury: on a ce fichier + lmdz_filtreg_p ! C'est totalement incompréhensible ;_;
! A minima il faut un nom clair pour chaque

SUBROUTINE filtreg_p( champ, ibeg, iend, nlat, nbniv, &
        ifiltre, iaire, griscal ,iter)
  USE parallel_lmdz, ONLY: OMP_CHUNK
  USE mod_filtre_fft
  USE lmdz_timer_filtre
  USE lmdz_coefils, ONLY: jfiltnu, jfiltnv, jfiltsu, jfiltsv, sddu, sddv, unsddu, unsddv, modfrstv, modfrstu
  USE lmdz_filtreg, ONLY: matrinvn, matrinvs, matriceun, matriceus, matricevn, matricevs

  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
  !
  !
  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"

  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)

  LOGICAL :: griscal
  INTEGER :: hemisph, iaire

  REAL :: champ_fft(iip1,nlat,nbniv)
  REAL :: champ_in(iip1,nlat,nbniv)

  LOGICAL,SAVE     :: first=.TRUE.
!$OMP THREADPRIVATE(first)

  REAL, DIMENSION(iip1,nlat,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

  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==1.OR.ifiltre==-1) &
        CALL abort_gcm("fitreg_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("fitreg_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("fitreg_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("fitreg_p","stopped",1)
   ENDIF
  !

  iim2   = iim * iim
  immjm  = iim * jjm
  !
  !
  IF( griscal )   THEN
     IF( nlat /= jjp1 )  THEN
        CALL abort_gcm("fitreg_p","nlat. NE. jjp1",1)
     ELSE
  !
        IF( iaire==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/=jjm )  THEN
        CALL abort_gcm("fitreg_p","nlat. NE. jjm",1)
     ELSE
  !
        IF( iaire==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==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==1 )      THEN

           IF( ifiltre==-2 )   THEN
              DO j = jdfil,jffil
#ifdef BLAS
                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
                       matrinvn(1,1,j), iim, &
                       champ_loc(1,j,1), iip1*nlat, 0.0, &
                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
#else
                 champ_fft(:iim,j-jdfil+1,:) &
                       =matmul(matrinvn(:,:,j),champ_loc(:iim,j,:))
#endif
              ENDDO

           ELSE IF ( griscal )     THEN
              DO j = jdfil,jffil
#ifdef BLAS
                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
                       matriceun(1,1,j), iim, &
                       champ_loc(1,j,1), iip1*nlat, 0.0, &
                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
#else
                 champ_fft(:iim,j-jdfil+1,:) &
                       =matmul(matriceun(:,:,j),champ_loc(:iim,j,:))
#endif
              ENDDO

           ELSE
              DO j = jdfil,jffil
#ifdef BLAS
                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
                       matricevn(1,1,j), iim, &
                       champ_loc(1,j,1), iip1*nlat, 0.0, &
                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
#else
                 champ_fft(:iim,j-jdfil+1,:) &
                       =matmul(matricevn(:,:,j),champ_loc(:iim,j,:))
#endif
              ENDDO

           ENDIF

        ELSE

           IF( ifiltre==-2 )   THEN
              DO j = jdfil,jffil
#ifdef BLAS
                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
                       matrinvs(1,1,j-jfiltsu+1), iim, &
                       champ_loc(1,j,1), iip1*nlat, 0.0, &
                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
#else
                 champ_fft(:iim,j-jdfil+1,:) &
                       =matmul(matrinvs(:,:,j-jfiltsu+1), &
                       champ_loc(:iim,j,:))
#endif
              ENDDO

           ELSE IF ( griscal )     THEN

              DO j = jdfil,jffil
#ifdef BLAS
                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
                       matriceus(1,1,j-jfiltsu+1), iim, &
                       champ_loc(1,j,1), iip1*nlat, 0.0, &
                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
#else
                 champ_fft(:iim,j-jdfil+1,:) &
                       =matmul(matriceus(:,:,j-jfiltsu+1), &
                       champ_loc(:iim,j,:))
#endif
              ENDDO

           ELSE

              DO j = jdfil,jffil
#ifdef BLAS
                 CALL DGEMM("N", "N", iim, nbniv_loc, iim, 1.0, &
                       matricevs(1,1,j-jfiltsv+1), iim, &
                       champ_loc(1,j,1), iip1*nlat, 0.0, &
                       champ_fft(1,j-jdfil+1,1), iip1*nlat)
#else
                 champ_fft(:iim,j-jdfil+1,:) &
                       =matmul(matricevs(:,:,j-jfiltsv+1), &
                       champ_loc(:iim,j,:))
#endif
              ENDDO

           ENDIF

        ENDIF
  ! c
        IF( ifiltre==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-jdfil+1,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_loc
              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-jdfil+1,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
              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 == -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== 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 )
              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)

  !
!$OMP MASTER
  CALL stop_timer
!$OMP END MASTER

END SUBROUTINE filtreg_p