! $Header$

SUBROUTINE groupeun(jjmax, llmax, q)

  USE comconst_mod, ONLY: ngroup

  IMPLICIT NONE

  include "dimensions.h"
  include "paramet.h"
  include "comgeom2.h"

  INTEGER :: jjmax, llmax
  REAL :: q(iip1, jjmax, llmax)

  ! INTEGER ngroup
  ! PARAMETER (ngroup=3)

  REAL :: airecn, qn
  REAL :: airecs, qs

  INTEGER :: i, j, l, ig, ig2, j1, j2, i0, jd

  !--------------------------------------------------------------------c
  ! Strategie d'optimisation                                           c
  ! stocker les valeurs systematiquement recalculees                   c
  ! et identiques d'un pas de temps sur l'autre. Il s'agit des         c
  ! aires des cellules qui sont sommees. S'il n'y a pas de changement  c
  ! de grille au cours de la simulation tout devrait bien se passer.   c
  ! Autre optimisation : determination des bornes entre lesquelles "j" c
  ! varie, au lieu de faire un test à chaque fois...
  !--------------------------------------------------------------------c

  INTEGER :: j_start, j_finish

  REAL, SAVE :: airen_tab(iip1, jjp1, 0:1)
  REAL, SAVE :: aires_tab(iip1, jjp1, 0:1)

  LOGICAL, SAVE :: first = .TRUE.
  ! INTEGER,SAVE :: i_index(iim,ngroup)
  INTEGER :: offset
  ! REAL         :: qsum(iim/ngroup)

  IF (first) THEN
    CALL INIT_GROUPEUN(airen_tab, aires_tab)
    first = .FALSE.
  ENDIF


  ! Champs 3D
  jd = jjp1 - jjmax
  !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
  DO l = 1, llm
    j1 = 1 + jd
    j2 = 2
    DO ig = 1, ngroup

      ! Concerne le pole nord
      j_start = j1 - jd
      j_finish = j2 - jd
      DO ig2 = 1, ngroup - ig + 1
        offset = 2**(ig2 - 1)
        DO j = j_start, j_finish
          !CDIR NODEP
          !CDIR ON_ADB(q)
          DO i0 = 1, iim, 2**ig2
            q(i0, j, l) = q(i0, j, l) + q(i0 + offset, j, l)
          ENDDO
        ENDDO
      ENDDO

      DO j = j_start, j_finish
        !CDIR NODEP
        !CDIR ON_ADB(q)
        DO i = 1, iim
          q(i, j, l) = q(i - MOD(i - 1, 2**(ngroup - ig + 1)), j, l)
        ENDDO
      ENDDO

      DO j = j_start, j_finish
        !CDIR ON_ADB(airen_tab)
        !CDIR ON_ADB(q)
        DO i = 1, iim
          q(i, j, l) = q(i, j, l) * airen_tab(i, j, jd)
        ENDDO
        q(iip1, j, l) = q(1, j, l)
      ENDDO

      !c     Concerne le pole sud
      j_start = j1 - jd
      j_finish = j2 - jd
      DO ig2 = 1, ngroup - ig + 1
        offset = 2**(ig2 - 1)
        DO j = j_start, j_finish
          !CDIR NODEP
          !CDIR ON_ADB(q)
          DO i0 = 1, iim, 2**ig2
            q(i0, jjp1 - j + 1 - jd, l) = q(i0, jjp1 - j + 1 - jd, l) &
                    + q(i0 + offset, jjp1 - j + 1 - jd, l)
          ENDDO
        ENDDO
      ENDDO

      DO j = j_start, j_finish
        !CDIR NODEP
        !CDIR ON_ADB(q)
        DO i = 1, iim
          q(i, jjp1 - j + 1 - jd, l) = q(i - MOD(i - 1, 2**(ngroup - ig + 1)), &
                  jjp1 - j + 1 - jd, l)
        ENDDO
      ENDDO

      DO j = j_start, j_finish
        !CDIR ON_ADB(aires_tab)
        !CDIR ON_ADB(q)
        DO i = 1, iim
          q(i, jjp1 - j + 1 - jd, l) = q(i, jjp1 - j + 1 - jd, l) * &
                  aires_tab(i, jjp1 - j + 1, jd)
        ENDDO
        q(iip1, jjp1 - j + 1 - jd, l) = q(1, jjp1 - j + 1 - jd, l)
      ENDDO

      j1 = j2 + 1
      j2 = j2 + 2**ig
    ENDDO
  ENDDO
  !$OMP END DO NOWAIT

  RETURN
END SUBROUTINE groupeun


SUBROUTINE INIT_GROUPEUN(airen_tab, aires_tab)

  USE comconst_mod, ONLY: ngroup

  IMPLICIT NONE

  include "dimensions.h"
  include "paramet.h"
  include "comgeom2.h"

  ! INTEGER ngroup
  ! PARAMETER (ngroup=3)

  REAL :: airen, airecn
  REAL :: aires, airecs

  INTEGER :: i, j, l, ig, j1, j2, i0, jd

  INTEGER :: j_start, j_finish

  REAL :: airen_tab(iip1, jjp1, 0:1)
  REAL :: aires_tab(iip1, jjp1, 0:1)

  DO jd = 0, 1
    j1 = 1 + jd
    j2 = 2
    DO ig = 1, ngroup

      ! c     Concerne le pole nord
      j_start = j1 - jd
      j_finish = j2 - jd
      DO j = j_start, j_finish
        DO i0 = 1, iim, 2**(ngroup - ig + 1)
          airen = 0.
          DO i = i0, i0 + 2**(ngroup - ig + 1) - 1
            airen = airen + aire(i, j)
          ENDDO
          DO i = i0, i0 + 2**(ngroup - ig + 1) - 1
            airen_tab(i, j, jd) = &
                    aire(i, j) / airen
          ENDDO
        ENDDO
      ENDDO

      ! c     Concerne le pole sud
      j_start = j1 - jd
      j_finish = j2 - jd
      DO j = j_start, j_finish
        DO i0 = 1, iim, 2**(ngroup - ig + 1)
          aires = 0.
          DO i = i0, i0 + 2**(ngroup - ig + 1) - 1
            aires = aires + aire(i, jjp1 - j + 1)
          ENDDO
          DO i = i0, i0 + 2**(ngroup - ig + 1) - 1
            aires_tab(i, jjp1 - j + 1, jd) = &
                    aire(i, jjp1 - j + 1) / aires
          ENDDO
        ENDDO
      ENDDO

      j1 = j2 + 1
      j2 = j2 + 2**ig
    ENDDO
  ENDDO

  RETURN
END SUBROUTINE INIT_GROUPEUN