Index: /LMDZ4/branches/LMDZ4-dev/libf/dyn3d/groupeun.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/dyn3d/groupeun.F	(revision 1086)
+++ /LMDZ4/branches/LMDZ4-dev/libf/dyn3d/groupeun.F	(revision 1087)
@@ -2,6 +2,6 @@
 ! $Header$
 !
-      subroutine groupeun(jjmax,llmax,q)
-      implicit none
+      SUBROUTINE groupeun(jjmax,llmax,q)
+      IMPLICIT NONE
 
 #include "dimensions.h"
@@ -10,51 +10,150 @@
 #include "comgeom2.h"
 
-      integer jjmax,llmax
-      real q(iip1,jjmax,llmax)
+      INTEGER jjmax,llmax
+      REAL q(iip1,jjmax,llmax)
 
-      integer ngroup
-      parameter (ngroup=3)
+      INTEGER ngroup
+      PARAMETER (ngroup=3)
 
-      real airen,airecn,qn
-      real aires,airecs,qs
+      REAL airecn,qn
+      REAL airecs,qs
 
-      integer i,j,l,ig,j1,j2,i0,jd
+      INTEGER i,j,l,ig,j1,j2,i0,jd
 
-Champs 3D
+c--------------------------------------------------------------------c 
+c Strategie d'optimisation                                           c
+c stocker les valeurs systematiquement recalculees                   c
+c et identiques d'un pas de temps sur l'autre. Il s'agit des         c
+c aires des cellules qui sont sommees. S'il n'y a pas de changement  c
+c de grille au cours de la simulation tout devrait bien se passer.   c
+c Autre optimisation : determination des bornes entre lesquelles "j" c
+c varie, au lieu de faire un test à chaque fois...
+c--------------------------------------------------------------------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.
+
+      IF (first) THEN
+         CALL INIT_GROUPEUN(airen_tab, aires_tab)
+         first = .FALSE.
+      ENDIF
+
+c Champs 3D
       jd=jjp1-jjmax
-      do l=1,llm
-      j1=1+jd
-      j2=2
-      do ig=1,ngroup
-         do j=j1-jd,j2-jd
-c           print*,'groupe ',ig,'  j= ',j,2**(ngroup-ig+1),'pts groupes'
-            do i0=1,iim,2**(ngroup-ig+1)
-               airen=0.
-               airecn=0.
-               qn=0.
-               aires=0.
-               airecs=0.
-               qs=0.
-               do i=i0,i0+2**(ngroup-ig+1)-1
-                  airen=airen+aire(i,j)
-                  aires=aires+aire(i,jjp1-j+1)
-                  qn=qn+q(i,j,l)
-                  qs=qs+q(i,jjp1-j+1-jd,l)
-               enddo
-               airecn=0.
-               airecs=0.
-               do i=i0,i0+2**(ngroup-ig+1)-1
-                  q(i,j,l)=qn*aire(i,j)/airen
-                  q(i,jjp1-j+1-jd,l)=qs*aire(i,jjp1-j+1)/aires
-               enddo
-            enddo
-            q(iip1,j,l)=q(1,j,l)
-            q(iip1,jjp1-j+1-jd,l)=q(1,jjp1-j+1-jd,l)
-         enddo
-         j1=j2+1
-         j2=j2+2**ig
-      enddo
-      enddo
 
-      return
-      end
+      DO l=1,llm
+         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)
+                  qn=0.
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     qn=qn+q(i,j,l)
+                  ENDDO
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     q(i,j,l)=qn*airen_tab(i,j,jd)
+                  ENDDO
+               ENDDO
+               q(iip1,j,l)=q(1,j,l)
+            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)
+                  qs=0.
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     qs=qs+q(i,jjp1-j+1-jd,l)
+                  ENDDO
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     q(i,jjp1-j+1-jd,l)=qs*aires_tab(i,jjp1-j+1,jd)
+                  ENDDO
+               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
+
+      RETURN
+      END
+
+
+
+      SUBROUTINE INIT_GROUPEUN(airen_tab, aires_tab)
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.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
Index: /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/groupeun_p.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/groupeun_p.F	(revision 1086)
+++ /LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/groupeun_p.F	(revision 1087)
@@ -1,5 +1,5 @@
-      subroutine groupeun_p(jjmax,llmax,jjb,jje,q)
+      SUBROUTINE groupeun_p(jjmax,llmax,jjb,jje,q)
       USE parallel
-      implicit none
+      IMPLICIT NONE
 
 #include "dimensions.h"
@@ -8,70 +8,155 @@
 #include "comgeom2.h"
 
-      integer jjmax,llmax,jjb,jje
-      real q(iip1,jjmax,llmax)
+      INTEGER jjmax,llmax,jjb,jje
+      REAL q(iip1,jjmax,llmax)
 
-      integer ngroup
-      parameter (ngroup=3)
+      INTEGER ngroup
+      PARAMETER (ngroup=3)
 
-      real airen,airecn,qn
-      real aires,airecs,qs
+      REAL airecn,qn
+      REAL airecs,qs
 
-      integer i,j,l,ig,j1,j2,i0,jd
+      INTEGER i,j,l,ig,j1,j2,i0,jd
 
-Champs 3D
+c--------------------------------------------------------------------c 
+c Strategie d'optimisation                                           c
+c stocker les valeurs systematiquement recalculees                   c
+c et identiques d'un pas de temps sur l'autre. Il s'agit des         c
+c aires des cellules qui sont sommees. S'il n'y a pas de changement  c
+c de grille au cours de la simulation tout devrait bien se passer.   c
+c Autre optimisation : determination des bornes entre lesquelles "j" c
+c varie, au lieu de faire un test à chaque fois...
+c--------------------------------------------------------------------c 
+
+      INTEGER j_start, j_finish
+
+      REAL, SAVE :: airen_tab(iip1,jjp1,0:1)
+      REAL, SAVE :: aires_tab(iip1,jjp1,0:1)
+!$OMP THREADPRIVATE(airen_tab, aires_tab)
+
+      LOGICAL, SAVE :: first = .TRUE.
+!$OMP THREADPRIVATE(first)
+
+      IF (first) THEN
+         CALL INIT_GROUPEUN_P(airen_tab, aires_tab)
+         first = .FALSE.
+      ENDIF
+
+c Champs 3D
       jd=jjp1-jjmax
 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-      do l=1,llm
-      j1=1+jd
-      j2=2
-      do ig=1,ngroup
-         do j=j1-jd,j2-jd
-c           print*,'groupe ',ig,'  j= ',j,2**(ngroup-ig+1),'pts groupes'
-            if ( j >= jjb .AND. j <= jje) THEN
-             
-              do i0=1,iim,2**(ngroup-ig+1)
-                 
-                 airen=0.
-                 airecn=0.
-                 qn=0.
-                 
-                 do i=i0,i0+2**(ngroup-ig+1)-1
-                    airen=airen+aire(i,j)
-                    qn=qn+q(i,j,l)
-                 enddo
-                 airecn=0.
-                 do i=i0,i0+2**(ngroup-ig+1)-1
-                   q(i,j,l)=qn*aire(i,j)/airen
-                 enddo
-              enddo
-              q(iip1,j,l)=q(1,j,l)
-              
-            endif
+      DO l=1,llm
+         j1=1+jd
+         j2=2
+         DO ig=1,ngroup
+
+c     Concerne le pole nord
+            j_start  = MAX(jjb, j1-jd)
+            j_finish = MIN(jje, j2-jd)
+            DO j=j_start, j_finish
+               DO i0=1,iim,2**(ngroup-ig+1)
+                  qn=0.
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     qn=qn+q(i,j,l)
+                  ENDDO
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     q(i,j,l)=qn*airen_tab(i,j,jd)
+                  ENDDO
+               ENDDO
+               q(iip1,j,l)=q(1,j,l)
+            ENDDO
+        
+!c     Concerne le pole sud
+            j_start  = MAX(1+jjp1-jje-jd, j1-jd)
+            j_finish = MIN(1+jjp1-jjb-jd, j2-jd)
+            DO j=j_start, j_finish
+               DO i0=1,iim,2**(ngroup-ig+1)
+                  qs=0.
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     qs=qs+q(i,jjp1-j+1-jd,l)
+                  ENDDO
+                  DO i=i0,i0+2**(ngroup-ig+1)-1
+                     q(i,jjp1-j+1-jd,l)=qs*aires_tab(i,jjp1-j+1,jd)
+                  ENDDO
+               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 INIT_GROUPEUN_P(airen_tab, aires_tab)
+
+      USE parallel
+      IMPLICIT NONE
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comconst.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
             
-            if ( jjp1-j+1-jd >= jjb .AND. jjp1-j+1-jd <= jje) THEN
-             
-              do i0=1,iim,2**(ngroup-ig+1)
-                 aires=0.
-                 airecs=0.
-                 qs=0.
-                 do i=i0,i0+2**(ngroup-ig+1)-1
-                    aires=aires+aire(i,jjp1-j+1)
-                    qs=qs+q(i,jjp1-j+1-jd,l)
-                 enddo
-                 airecs=0.
-                 do i=i0,i0+2**(ngroup-ig+1)-1
-                   q(i,jjp1-j+1-jd,l)=qs*aire(i,jjp1-j+1)/aires
-                 enddo
-              enddo
-              q(iip1,jjp1-j+1-jd,l)=q(1,jjp1-j+1-jd,l)
-           
-            endif
-         enddo
-             
-           j1=j2+1
-           j2=j2+2**ig
-      enddo
-      enddo
-c$OMP END DO NOWAIT
-      return
-      end
+!     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
