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