MODULE sedimorpho

#if defined key_sedim_MUSTANG && defined key_MARS

   !&E==========================================================================
   !&E                   ***  MODULE  sedimorpho  ***
   !&E
   !&E
   !&E ** Purpose : concerns subroutines related to sediment morphodynamic
   !&E              to be used in MARSiAM system
   !&E 
   !&E ** Description :
   !&E
   !&E     subroutine sed_morphoinit     !
   !&E     subroutine sed_morpho         !
   !&E
   !&E ** History :
   !&E     !  2009-03  (Pierre Le Hir, Florence Cayocca)
   !&E     ! 2014-05  (F. Cayocca)   beginning of update morphodynamic (not operationnal - TO SEE LATER) -  
   !&E     ! 2015-12  (B.Thouvenin) : extraction of routines 1DV for reorganization of module SEDIMARS
   !&E
   !&E==========================================================================

#include "toolcpp.h"

   !! * Modules used

   USE comMUSTANG

   IMPLICIT NONE

   !! * Accessibility

   ! functions & routines of this module, called outside :
   ! PUBLIC functions                                
   PUBLIC sed_morphoinit,sed_morpho,sed_alloc_morpho
#if defined key_toseelater
   ! creation des masques "mur" (10 mai 2012)
   ! cf. les routines morphoinit et morpho
   LOGICAL, ALLOCATABLE,DIMENSION(:,:)      :: mask_wallx,mask_wally 
#endif
   !---------------------------------------------------------

 
   !! * Private variables



 CONTAINS
 
    !!===========================================================================
 
  SUBROUTINE sed_alloc_morpho
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_alloc_morpho  ***
   !&E
   !&E ** Purpose : 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  main 
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  march 2009 (P. Le Hir) enables morphodynamics computation
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
   !!===========================================================================
    !!--------------------------------------------------------------------------  
   !! * Executable part

      PRINT_DBG*, 'DANS SED_ALLOC_MORPHO'
      
    ALLOCATE(morphox(liminm1:limax,ljmin:ljmax))
    ALLOCATE(morphoy(limin:limax,ljminm1:ljmax))
    ALLOCATE(hxi(liminm1:limax,ljmin:ljmax))
    ALLOCATE(hyi(limin:limax,ljminm1:ljmax))
    ALLOCATE(h0i(limin:limaxp1,ljmin:ljmaxp1))
    hxi(liminm1:limax,ljmin:ljmax)=0.0_rsh
    hyi(limin:limax,ljminm1:ljmax)=0.0_rsh

     PRINT_DBG*, 'FIN SED_ALLOC_MORPHO'

  END SUBROUTINE sed_alloc_morpho     

   !!===========================================================================
  
  SUBROUTINE sed_morphoinit
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_morphoinit  ***
   !&E
   !&E ** Purpose : 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  main 
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  march 2009 (P. Le Hir) enables morphodynamics computation
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,    ONLY : ig,id,jb,jh,hx,hy,h0,h0fond,l_obc_cycl_x,l_obc_cycl_y
   USE comvarp2d,    ONLY : ssh
   USE parameters,   ONLY : valmanq
#ifdef key_castest
   USE parameters,   ONLY : num_testcase
#endif

   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_j_rsh

   !! * Local declarations
   INTEGER                  :: i,j,k          
   REAL(KIND=rsh)           :: hultra

   !!--------------------------------------------------------------------------  
   !! * Executable part

      PRINT_DBG*, 'DANS SED_MORPHO_INIT'

#ifdef key_toseelater
      ! mask_wallx(i,j) & mask_wally(i,j) indicate whether morphodynamic coupling  
      ! is effective for hx(i,j) or hy(i,j) - coupling only if mask i set to .TRUE.
      mask_wallx(:,:)=.FALSE.
      mask_wally(:,:)=.FALSE.
      DO j=ljmin,ljmax
      DO i=limin,limax            ! modification des boucles MPI
         IF(hx(i,j)/=-valmanq-h0fond) mask_wallx(i,j)=.TRUE.
         IF(hy(i,j)/=-valmanq-h0fond) mask_wally(i,j)=.TRUE.
      ENDDO
      ENDDO
#endif

      morphox(liminm1:limax,ljmin  :ljmax)=1.0_rsh
      morphoy(limin  :limax,ljminm1:ljmax)=1.0_rsh

      DO j=ljminp1,ljmaxm1
      DO i=MAX0(liminp1,ig(j)+1),MIN0(limaxm1,id(j)-1)
        hsed(i,j)=0.0_rsh
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
          DO k=ksmi(i,j),ksma(i,j)
            hsed(i,j)=hsed(i,j)+dzs(k,i,j)
          ENDDO
        ENDIF
      ENDDO
      ENDDO

      ! hypothese gradient nul aux frontieres
      IF_MPI (ljmax==jmax) THEN
        DO i=MAX0(limin,ig(jmax)+1),MIN0(limax,id(jmax)-1)
            hsed(i,jmax)=hsed(i,jmax-1)
        ENDDO
      ENDIF_MPI
      IF_MPI (ljmin==jmin) THEN
        DO i=MAX0(limin,ig(jmin)+1),MIN0(limax,id(jmin)-1)
            hsed(i,jmin+1)=hsed(i,jmin+2)
        ENDDO
      ENDIF_MPI
      IF_MPI (limax == imax) THEN
#ifdef key_castest
      IF (num_testcase==41) THEN   ! testcase estuar
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
           hsed(imax,j)=0.0_rsh
        ENDDO
      ELSE
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
           hsed(imax,j)=hsed(imax-1,j)
        ENDDO
      ENDIF
#else
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
           hsed(imax,j)=hsed(imax-1,j)
        ENDDO
#endif
      ENDIF_MPI

      IF_MPI(limin==imin) THEN
        DO j=MAX0(ljmin,jb(imin)+1),MIN0(ljmax,jh(imin)-1)
           hsed(imin+1,j)=hsed(imin+2,j)
        ENDDO
      ENDIF_MPI

      ! les quatre coins 
      IF_MPI (limax == imax .AND. ljmax == jmax) THEN
         hsed(limax,ljmax)=hsed(limax-1,ljmax-1)
      ENDIF_MPI

     IF_MPI (limax==imax .AND. ljmin==jmin) THEN
        hsed(limax,ljmin+1)=hsed(limax-1,ljmin+2)
     ENDIF_MPI

     IF_MPI (limin==imin .AND. ljmin==jmin) THEN
        hsed(limin+1,ljmin+1)=hsed(limin+2,ljmin+2)
     ENDIF_MPI

     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
        hsed(limin+1,ljmax)=hsed(limin+2,ljmax-1)
     ENDIF_MPI

#ifdef key_toseelater
!  PROBLEME : hsed n est pas alloue de imin a imax et de jmin a jmax
!             il faut les echanger avec les processeurs extremes et les ranger dans m2 ou p2
!             voir comme dans tsobc3dapply par exemple
     ! condition de periodicite sur hsed
     IF(l_obc_cycl_x) THEN
       DO j=ljminp1,ljmax
          hsed(imin+1,j)=hsed(imax-1,j)
          hsed(imax,j)=hsed(imin+2,j)
       ENDDO
     ENDIF
     ! condition de periodicite sur hsed
     IF(l_obc_cycl_y) THEN
       DO i=liminp1,limax
          hsed(i,jmin+1)=hsed(i,jmax-1)
          hsed(i,jmax)=hsed(i,jmin+2)
       ENDDO
     ENDIF
#endif

   CALL_MPI ex_i_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,hsed(limin:limaxp1,ljmin:ljmaxp1))
   CALL_MPI ex_j_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,hsed(limin:limaxp1,ljmin:ljmaxp1))

#ifdef key_toseelater
      DO j=ljminp1,ljmaxm1
      DO i=MAX0(liminp1,ig(j)+1),MIN0(limaxm1,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
          IF(mask_wallx(i,j)) THEN
             hxi(i,j)=hx(i,j)+0.5_rsh*(hsed(i,j)+hsed(i+1,j))*morphox(i,j)
          ELSE
             hxi(i,j)=hx(i,j)
          ENDIF
          IF(mask_wally(i,j)) THEN
             hyi(i,j)=hy(i,j)+0.5_rsh*(hsed(i,j)+hsed(i,j+1))*morphoy(i,j)
          ELSE
             hyi(i,j)=hy(i,j)
          ENDIF
        ENDIF
      ENDDO
      ENDDO
#endif

      DO j=ljmin,ljmax
      DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
             hxi(i,j)=hx(i,j)+0.5_rsh*(hsed(i,j)+hsed(i+1,j))*morphox(i,j)
             hyi(i,j)=hy(i,j)+0.5_rsh*(hsed(i,j)+hsed(i,j+1))*morphoy(i,j)
        ENDIF
      ENDDO
      ENDDO

     ! limites
      IF_MPI (ljmax==jmax) THEN
        DO i=MAX0(limin,ig(jmax)+1),MIN0(limax,id(jmax)-1)
          hxi(i,ljmax)=hxi(i,ljmax-1)
          hyi(i,ljmax)=hyi(i,ljmax-1)
        ENDDO
      ENDIF_MPI
      IF_MPI (ljmin==jmin) THEN
        DO i=MAX0(limin,ig(jmin)+1),MIN0(limax,id(jmin)-1)
         hxi(i,ljmin+1)=hxi(i,ljmin+2)
         hyi(i,ljmin+1)=hyi(i,ljmin+2)
        ENDDO
      ENDIF_MPI
      IF_MPI (limax == imax) THEN
#ifdef key_castest
      IF (num_testcase==41) THEN
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
         hxi(limax,j)=hx(limax,j)
         hyi(limax,j)=hy(limax,j)
        ENDDO
      ELSE
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
         hxi(limax,j)=hxi(limax-1,j)
         hyi(limax,j)=hyi(limax-1,j)
        ENDDO
      ENDIF
#else
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
         hxi(limax,j)=hxi(limax-1,j)
         hyi(limax,j)=hyi(limax-1,j)
        ENDDO
#endif
      ENDIF_MPI
      IF_MPI(limin==imin) THEN
        DO j=MAX0(ljmin,jb(imin)+1),MIN0(ljmax,jh(imin)-1)
          hxi(limin+1,j)=hxi(limin+2,j)
          hyi(limin+1,j)=hyi(limin+2,j)
        ENDDO
      ENDIF_MPI

      ! les quatre coins
      IF_MPI (limax == imax .AND. ljmax == jmax) THEN
        hxi(limax,ljmax)=hxi(limax-1,ljmax-1)
        hyi(limax,ljmax)=hyi(limax-1,ljmax-1)
      ENDIF_MPI
      IF_MPI (limax==imax .AND. ljmin==jmin) THEN
        hxi(limax,jmin+1)=hxi(limax-1,ljmin+2)
        hyi(limax,jmin+1)=hxi(limax-1,ljmin+2)
      ENDIF_MPI
      IF_MPI (limin==imin .AND. ljmin==jmin) THEN
        hxi(limin+1,jmin+1)=hxi(limin+2,jmin+2)
        hyi(limin+1,jmin+1)=hyi(limin+2,jmin+2)
      ENDIF_MPI
      IF_MPI (limin==imin .AND. ljmax==jmax) THEN
        hxi(limin+1,jmax)=hxi(limin+2,jmax-1)
        hyi(limin+1,jmax)=hyi(limin+2,jmax-1)
      ENDIF_MPI

#ifdef key_toseelater
!  PROBLEME : hxi,hyi ne sont pas alloues de imin a imax et de jmin a jmax
!             il faut les echanger avec les processeurs extremes et les ranger dans m2 ou p2
!             voir comme dans tsobc3dapply par exemple

      ! condition de periodicite dans la direction x 
      IF(l_obc_cycl_x) THEN
        DO j=ljminp1,ljmax 
          IF_MPI (limin==imin) THEN 
          hxi(imin+1,j)=hxi(imax-1,j)
          hxi(imax,j)=hxi(imin+2,j)
          hyi(imin+1,j)=hyi(imax-1,j)
          hyi(imax,j)=hyi(imin+2,j)
        ENDDO
      ENDIF

      ! condition de periodicite dans la direction y
      IF(l_obc_cycl_y) THEN
        DO i=liminp1,limax
          hxi(i,jmin+1)=hxi(i,jmax-1)
          hxi(i,jmax)=hxi(i,jmin+2)
          hyi(i,jmin+1)=hyi(i,jmax-1)
          hyi(i,jmax)=hyi(i,jmin+2)
        ENDDO
      ENDIF
#endif

      ! evaluation de la cote des mailles
      ! ---------------------------------
      
      DO j=ljmin,ljmax
      DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
          h0(i,j)=MAX(hx(i,j),hy(i,j),hx(i-1,j),hy(i,j-1))+h0fond
          ssh(i,j)=MAX(ssh(i,j),-h0(i,j))
        ENDIF      
      ENDDO
      ENDDO

     ! modification des conditions aux frontieres
     ! hypothese gradient nul aux frontieres
#ifdef key_castest
     IF_MPI (ljmax == jmax) THEN
        DO i=MAX0(limin,ig(jmax)+1),MIN0(limax,id(jmax)-1)
          IF (.NOT. num_testcase==41) h0(i,jmax)=h0(i,jmax-1)
          ssh(i,jmax)=ssh(i,jmax-1)
        ENDDO
     ENDIF_MPI
     IF_MPI (ljmin == jmin) THEN
        DO i=MAX0(limin,ig(jmin)+1),MIN0(limax,id(jmin)-1)
          IF (.NOT. num_testcase==41) h0(i,jmin+1)=h0(i,jmin+2)
          ssh(i,jmin+1)=ssh(i,jmin+2)
        ENDDO
     ENDIF_MPI
     IF_MPI (limax==imax) THEN
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
           IF (.NOT. num_testcase==41) h0(imax,j)=h0(imax-1,j)
           ssh(imax,j)=ssh(imax-1,j)
        ENDDO
     ENDIF_MPI
     IF_MPI (limin == imin) THEN
        DO j=MAX0(ljmin,jb(imin)+1),MIN0(ljmax,jh(imin)-1)
          IF (.NOT. num_testcase==41) h0(imin+1,j)=h0(imin+2,j)
          ssh(imin+1,j)=ssh(imin+2,j)
        ENDDO
     ENDIF_MPI
#else
     IF_MPI (ljmax == jmax) THEN
        DO i=MAX0(limin,ig(jmax)+1),MIN0(limax,id(jmax)-1)
          h0(i,jmax)=h0(i,jmax-1)
          ssh(i,jmax)=ssh(i,jmax-1)
        ENDDO
     ENDIF_MPI
     IF_MPI (ljmin == jmin) THEN
        DO i=MAX0(limin,ig(jmin)+1),MIN0(limax,id(jmin)-1)
          h0(i,jmin+1)=h0(i,jmin+2)
          ssh(i,jmin+1)=ssh(i,jmin+2)
        ENDDO
     ENDIF_MPI
     IF_MPI (limax==imax) THEN
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
           h0(imax,j)=h0(imax-1,j)
           ssh(imax,j)=ssh(imax-1,j)
        ENDDO
     ENDIF_MPI
     IF_MPI (limin == imin) THEN
        DO j=MAX0(ljmin,jb(imin)+1),MIN0(ljmax,jh(imin)-1)
          h0(imin+1,j)=h0(imin+2,j)
          ssh(imin+1,j)=ssh(imin+2,j)
        ENDDO
     ENDIF_MPI
#endif

     ! les quatre coins 
     IF_MPI (limax == imax .AND. ljmax == jmax) THEN
        h0(imax,jmax)=h0(imax-1,jmax-1)
        ssh(imax,jmax)=ssh(imax-1,jmax-1)
     ENDIF_MPI
     IF_MPI (limax==imax .AND. ljmin==jmin) THEN
        h0(imax,jmin+1)=h0(imax-1,jmin+2)
        ssh(imax,jmin+1)=ssh(imax-1,jmin+2)
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmin==jmin) THEN
        h0(imin+1,jmin+1)=h0(imin+2,jmin+2)
        ssh(imin+1,jmin+1)=ssh(imin+2,jmin+2)
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
        h0(imin+1,jmax)=h0(imin+2,jmax-1)
        ssh(imin+1,jmax)=ssh(imin+2,jmax-1)
     ENDIF_MPI
!
!  initial updating of  hx and hy:
     CALL_MPI ex_i_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,h0(limin:limaxp1,ljmin:ljmaxp1))
     CALL_MPI ex_j_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,h0(limin:limaxp1,ljmin:ljmaxp1))
      DO j=ljmin,MIN0(ljmax,jmax-1)   !(avant : jusqu a ljmax)
      DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-2)
            hx(i,j)=MIN(h0(i,j),h0(i+1,j))-h0fond
            hy(i,j)=MIN(h0(i,j),h0(i,j+1))-h0fond
     enddo
     enddo
    
#ifdef key_toseelater
!  PROBLEME : h0 , ssh ne sont pas alloues de imin a imac et de jmin a jmax
!             il faut les echanger avec les processeurs extremes et les ranger dans m2 ou p2
!             voir comme dans tsobc3dapply par exemple
     ! condition de periodicite dans la direction x
     IF(l_obc_cycl_x) THEN
        DO j=ljminp1,ljmax
          h0(imin+1,j)=h0(imax-1,j)
          h0(imax,j)=h0(imin+2,j)
          ssh(imin+1,j)=ssh(imax-1,j)
          ssh(imax,j)=ssh(imin+2,j)
        ENDDO
     ENDIF

     ! condition de periodicite dans la direction y
     IF(l_obc_cycl_y) THEN
        DO i=liminp1,limax
          h0(i,jmin+1)=h0(i,jmax-1)
          h0(i,jmax)=h0(i,jmin+2)
          ssh(i,jmin+1)=ssh(i,jmax-1)
          ssh(i,jmax)=ssh(i,jmin+2)
        ENDDO
     ENDIF
#endif

#ifdef key_toseelater
!ATTENTION : hupper et dhcell non utilises, a integrer a l avenir pour bien gerer 
!            le decouvrement (coherence avec fwet de SIAM) - voir PLH
      hultra=10000.0_rsh
      DO j=ljmin,ljmax
      DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN

          hupper(i,j)=MIN(hx(i,j)*morphox(i,j)+(1.0_rsh-morphox(i,j))*hultra,         &
                          hx(i-1,j)*morphox(i-1,j)+(1.0_rsh-morphox(i-1,j))*hultra,   &
                          hy(i,j)*morphoy(i,j)+(1.0_rsh-morphoy(i,j))*hultra,         &
                          hy(i,j-1)*morphoy(i,j-1)+(1.0_rsh-morphoy(i,j-1))*hultra)
          IF(hupper(i,j).GE.hultra)hupper(i,j)=-99.99_rsh
          dhcell(i,j)=h0(i,j)-hupper(i,j)
        ENDIF
      ENDDO  
      ENDDO        
#endif

! h0i : total depth invariable (while h0 and hsed may vary)
      do j=ljmin,ljmax
      do i=limin,limax
        h0i(i,j)=h0(i,j)+hsed(i,j)
      enddo
      enddo

    CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
    CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
    CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,h0(liminm2:limaxp2,ljminm2:ljmaxp2))
    CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,h0(liminm2:limaxp2,ljminm2:ljmaxp2))
   ! exchange hupper and dhcell in the future if required
   ! hxi, hyi are not exchanged because neighboring cells are not used

      PRINT_DBG*, 'FIN SED_MORPHO_INIT'

  END SUBROUTINE sed_morphoinit     

   !!===========================================================================
 
  SUBROUTINE sed_morpho
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_morpho  ***
   !&E
   !&E ** Purpose : 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  step, at the end
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-03  (P. Le Hir) enables morphodynamics computation
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,    ONLY : ig,id,jb,jh,hx,hy,h0,h0fond,l_obc_cycl_x,l_obc_cycl_y
   USE comvars2d,    ONLY : htx,hty
   USE comvarp2d,    ONLY : ssh
   USE parameters,   ONLY : valmanq
   USE comMUSTANG
   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_j_rsh
#ifdef key_castest
   USE parameters,   ONLY : num_testcase
#endif

   !! * Local declarations
   INTEGER                  :: i,j,k
   REAL(KIND=rsh)           :: hultra                
   
   !!--------------------------------------------------------------------------
   !! * Executable part

      PRINT_DBG*, 'DANS SED_MORPHO'

      DO j=ljmin,ljmax
        DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
         IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN

          hsed(i,j)=0.0_rsh
          DO k=ksmi(i,j),ksma(i,j)
            hsed(i,j)=hsed(i,j)+dzs(k,i,j)
          ENDDO
         ENDIF
        ENDDO
      ENDDO
      ! hypothese gradient nul aux frontieres

      IF_MPI (ljmax==jmax) THEN
        DO i=MAX0(limin,ig(ljmax)+1),MIN0(limax,id(ljmax)-1)
            hsed(i,ljmax)=hsed(i,ljmax-1)
        ENDDO
      ENDIF_MPI
      IF_MPI (ljmin==jmin) THEN
        DO i=MAX0(limin,ig(ljmin)+1),MIN0(limax,id(ljmin)-1)
            hsed(i,ljmin+1)=hsed(i,ljmin+2)
        ENDDO
      ENDIF_MPI
#ifdef key_castest
      IF (num_testcase.NE.41) THEN
#endif
        IF_MPI (limax == imax) THEN
          DO j=MAX0(ljmin,jb(limax)+1),MIN0(ljmax,jh(limax)-1)
             hsed(limax,j)=hsed(limax-1,j)
          ENDDO
        ENDIF_MPI
#ifdef key_castest
      ENDIF
#endif

      IF_MPI(limin==imin) THEN
        DO j=MAX0(ljmin,jb(limin)+1),MIN0(ljmax,jh(limin)-1)
           hsed(limin+1,j)=hsed(limin+2,j)
        ENDDO
      ENDIF_MPI

      ! les quatre coins

      IF_MPI (limax == imax .AND. ljmax == jmax) THEN
         hsed(limax,ljmax)=hsed(limax-1,ljmax-1)
      ENDIF_MPI

     IF_MPI (limax==imax .AND. ljmin==jmin) THEN
        hsed(limax,ljmin+1)=hsed(limax-1,ljmin+2)
     ENDIF_MPI

     IF_MPI (limin==imin .AND. ljmin==jmin) THEN
        hsed(limin+1,ljmin+1)=hsed(limin+2,ljmin+2)
     ENDIF_MPI

     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
        hsed(limin+1,ljmax)=hsed(limin+2,ljmax-1)
     ENDIF_MPI

#ifdef key_toseelater
!  PROBLEME : hsed n est pas alloue de imin a imax et de jmin a jmax
!             il faut les echanger avec les processeurs extremes et les ranger dans m2 ou p2
!             voir comme dans tsobc3dapply par exemple
     ! condition de periodicite sur hsed
     IF(l_obc_cycl_x) THEN
       DO j=ljminp1,ljmax
          hsed(imin+1,j)=hsed(imax-1,j)
          hsed(imax,j)=hsed(imin+2,j)
       ENDDO
     ENDIF
     ! condition de periodicite sur hsed
     IF(l_obc_cycl_y) THEN
       DO i=liminp1,limax
         hsed(i,jmin+1)=hsed(i,jmax-1)
         hsed(i,jmax)=hsed(i,jmin+2)
       ENDDO
     ENDIF
#endif

     DO j=ljmin,ljmax
     DO i=limin,limax
       h0(i,j)=h0i(i,j)-hsed(i,j)
     ENDDO
     ENDDO

     CALL_MPI ex_i_rsh(-1,2,1,liminm1,limaxp2,ljminm1,ljmaxp2,hsed(liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL_MPI ex_j_rsh(-1,2,1,liminm1,limaxp2,ljminm1,ljmaxp2,hsed(liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL_MPI ex_i_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,h0(limin:limaxp1,ljmin:ljmaxp1))
     CALL_MPI ex_j_rsh(0,1,1,limin,limaxp1,ljmin,ljmaxp1,h0(limin:limaxp1,ljmin:ljmaxp1))

#ifdef key_toseelater  
     ! probleme sur les masques
      DO j=ljmin,ljmax
        DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
         IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
          ! METHODE M1 : modification des hx/hy par la moyenne des milieux de mailles de part et d autre 
           IF(mask_wallx(i,j)) THEN
             hx(i,j)=hxi(i,j)-0.5_rsh*(hsed(i,j)+hsed(i+1,j))*morphox(i,j)
          ENDIF
          IF(mask_wally(i,j)) THEN
             hy(i,j)=hyi(i,j)-0.5_rsh*(hsed(i,j)+hsed(i,j+1))*morphoy(i,j)
          ENDIF
          ! METHODE M2i : modification des hx/hy par la maille voisine la plus basse
          ! hx(i,j)=hxi(i,j)-MIN(hsed(i,j),hsed(i+1,j))*morphox(i,j)
          ! hy(i,j)=hyi(i,j)-MIN(hsed(i,j),hsed(i,j+1))*morphoy(i,j)
         ENDIF
        ENDDO
      ENDDO
      
     ! modification des conditions aux frontieres 
      IF_MPI (ljmax==jmax) THEN
        DO i=MAX0(limin,ig(ljmax)+1),MIN0(limax,id(ljmax)-1)
           IF(mask_wallx(i,ljmax)) hx(i,ljmax)=hx(i,ljmax-1)
           IF(mask_wally(i,ljmax)) hy(i,ljmax)=hy(i,ljmax-1)
        ENDDO
      ENDIF_MPI
      IF_MPI (ljmin==jmin) THEN
        DO i=MAX0(limin,ig(ljmin)+1),MIN0(limax,id(ljmin)-1)
            IF(mask_wallx(i,ljmin+1)) hx(i,ljmin+1)=hx(i,ljmin+2)
            IF(mask_wally(i,ljmin+1)) hy(i,ljmin+1)=hy(i,ljmin+2)
        ENDDO
      ENDIF_MPI

      IF_MPI (limax == imax) THEN
        DO j=MAX0(ljmin,jb(limax)+1),MIN0(ljmax,jh(limax)-1)
           IF(mask_wallx(limax,j)) hx(limax,j)=hx(limax-1,j)
           IF(mask_wally(limax,j)) hy(limax,j)=hy(limax-1,j)
        ENDDO
      ENDIF_MPI

      IF_MPI(limin==imin) THEN
        DO j=MAX0(ljmin,jb(limin)+1),MIN0(ljmax,jh(limin)-1)
           IF(mask_wallx(limin+1,j)) hx(limin+1,j)=hx(limin+2,j)
           IF(mask_wally(limin+1,j)) hy(limin+1,j)=hy(limin+2,j)
        ENDDO
      ENDIF_MPI

      ! les quatre coins 
      IF_MPI (limax == imax .AND. ljmax == jmax) THEN
         IF(mask_wallx(limax,ljmax)) hx(limax,ljmax)=hx(limax-1,ljmax-1)
         IF(mask_wally(limax,ljmax)) hy(limax,ljmax)=hy(limax-1,ljmax-1)
      ENDIF_MPI

     IF_MPI (limax==imax .AND. ljmin==jmin) THEN
        IF(mask_wallx(limax,ljmin+1)) hx(limax,ljmin+1)=hx(limax-1,ljmin+2)
        IF(mask_wally(limax,ljmin+1)) hy(limax,ljmin+1)=hy(limax-1,ljmin+2)
     ENDIF_MPI

     IF_MPI (limin==imin .AND. ljmin==jmin) THEN
        IF(mask_wallx(limin+1,ljmin+1)) hx(limin+1,ljmin+1)=hx(limin+2,ljmin+2)
        IF(mask_wally(limin+1,ljmin+1)) hy(limin+1,ljmin+1)=hy(limin+2,ljmin+2)
     ENDIF_MPI

     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
        IF(mask_wallx(limin+1,ljmax)) hx(limin+1,ljmax)=hx(limin+2,ljmax-1)
        IF(mask_wally(limin+1,ljmax)) hy(limin+1,ljmax)=hy(limin+2,ljmax-1)
     ENDIF_MPI

     ! actualisation des masques
     mask_wallx(:,:)=.FALSE.
     mask_wally(:,:)=.FALSE.
     DO j=ljmin,ljmax
     DO i=limin,limax
       IF(hx(i,j)/=-valmanq-h0fond) mask_wallx(i,j)=.TRUE.
       IF(hy(i,j)/=-valmanq-h0fond) mask_wally(i,j)=.TRUE.
     ENDDO
     ENDDO
#endif 

      IF (l_bathy_smoothing) THEN
        DO j=ljmin,ljmax
        DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-2)
          IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-2) THEN
            hx(i,j)=hxi(i,j)-0.3_rsh*(hsed(i,j)+hsed(i+1,j))                                  &
                          -0.1_rsh*(hsed(i-1,j)+hsed(i+2,j))                                &
                          -0.05_rsh*(hsed(i,j+1)+hsed(i,j-1)+hsed(i+1,j+1)+hsed(i+1,j-1))
            hy(i,j)=hyi(i,j)-0.3_rsh*(hsed(i,j)+hsed(i,j+1))                                  &
                          -0.1_rsh*(hsed(i,j-1)+hsed(i,j+2))                                &
                          -0.05_rsh*(hsed(i-1,j)+hsed(i+1,j)+hsed(i-1,j+1)+hsed(i+1,j+1))       
          ENDIF
        ENDDO
        ENDDO
      ELSE
        DO j=ljmin,ljmaxm1
        DO i=limin,limaxm1
          IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
            hx(i,j)=MIN(h0(i,j),h0(i+1,j))-h0fond
            hy(i,j)=MIN(h0(i,j),h0(i,j+1))-h0fond
          ENDIF
        ENDDO
        ENDDO
      ENDIF

      ! modification des conditions aux frontieres 
      IF_MPI (ljmax==jmax) THEN
        DO i=MAX0(limin,ig(ljmax)+1),MIN0(limax,id(ljmax)-1)
#ifdef key_castest_estuary
          hx(i,ljmax)=-7.0_rsh
          hy(i,ljmax)=-7.0_rsh
          hx(i,ljmax-1)=-7.0_rsh
          hy(i,ljmax-1)=-7.0_rsh
          hy(i,ljmax-2)=-7.0_rsh
#else
          hx(i,ljmax)=hx(i,ljmax-1)
          hy(i,ljmax)=hy(i,ljmax-1)
#endif
        ENDDO
      ENDIF_MPI
      IF_MPI (ljmin==jmin) THEN
        DO i=MAX0(limin,ig(ljmin)+1),MIN0(limax,id(ljmin)-1)
#ifdef key_castest_estuary
          hx(i,ljmin)=-7.0_rsh
          hy(i,ljmin)=-7.0_rsh
          hx(i,ljmin+1)=-7.0_rsh
          hy(i,ljmin+1)=-7.0_rsh
#else
          hx(i,ljmin+1)=hx(i,ljmin+2)
          hy(i,ljmin+1)=hy(i,ljmin+2)
#endif
        ENDDO
      ENDIF_MPI
      IF_MPI (limax == imax) THEN
        DO j=MAX0(ljmin,jb(limax)+1),MIN0(ljmax,jh(limax)-1)
#ifdef key_castest_estuary
          hx(limax,j)=-7.0_rsh
          hy(limax,j)=-7.0_rsh
          hx(limax-1,j)=-7.0_rsh
#else
          hx(limax,j)=hx(limax-1,j)
          hy(limax,j)=hy(limax-1,j)
#endif
        ENDDO
      ENDIF_MPI
      IF_MPI(limin==imin) THEN
        DO j=MAX0(ljmin,jb(limin)+1),MIN0(ljmax,jh(limin)-1)
          hx(limin+1,j)=hx(limin+2,j)
          hy(limin+1,j)=hy(limin+2,j)
        ENDDO
      ENDIF_MPI

      ! les quatre coins 
      IF_MPI (limax == imax .AND. ljmax == jmax) THEN
        hx(limax,ljmax)=hx(limax-1,ljmax-1)
        hy(limax,ljmax)=hy(limax-1,ljmax-1)
      ENDIF_MPI
      IF_MPI (limax==imax .AND. ljmin==jmin) THEN
        hx(limax,ljmin+1)=hx(limax-1,ljmin+2)
        hy(limax,ljmin+1)=hy(limax-1,ljmin+2)
      ENDIF_MPI
      IF_MPI (limin==imin .AND. ljmin==jmin) THEN
        hx(limin+1,ljmin+1)=hx(limin+2,ljmin+2)
        hy(limin+1,ljmin+1)=hy(limin+2,ljmin+2)
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
        hx(limin+1,ljmax)=hx(limin+2,ljmax-1)
        hy(limin+1,ljmax)=hy(limin+2,ljmax-1)
     ENDIF_MPI

#ifdef key_toseelater
!  PROBLEME : hx et hy ne sont pas alloues de imin a imax et de jmin a jmax
!             il faut les echanger avec les processeurs extremes et les ranger dans m2 ou p2
!             voir comme dans tsobc3dapply par exemple
     ! condition de periodicite dans la direction x
     IF(l_obc_cycl_x) THEN
       DO j=ljminp1,ljmax
          hx(imin+1,j)=hx(imax-1,j)
          hx(imax,j)=hx(imin+2,j)
          hy(imin+1,j)=hy(imax-1,j)
          hy(imax,j)=hy(imin+2,j)
       ENDDO
     ENDIF
     
     ! condition de periodicite dans la direction y
     IF(l_obc_cycl_y) THEN
       DO i=liminp1,limax
          hx(i,jmin+1)=hx(i,jmax-1)
          hx(i,jmax)=hx(i,jmin+2)
          hy(i,jmin+1)=hy(i,jmax-1)
          hy(i,jmax)=hy(i,jmin+2)
       ENDDO
     ENDIF
#endif

      ! evaluation de la cote des mailles
      !---------------------------------

      DO j=ljmin,ljmax
        DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
         IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
            ssh(i,j)=MAX(ssh(i,j),-h0(i,j))
         ENDIF
        ENDDO
      ENDDO

     ! hypothese gradient nul aux frontieres 
     IF_MPI (ljmax == jmax) THEN
        DO i=MAX0(limin,ig(jmax)+1),MIN0(limax,id(jmax)-1)
          h0(i,ljmax)=h0(i,ljmax-1)
          ssh(i,ljmax)=ssh(i,ljmax-1)
        ENDDO
     ENDIF_MPI

     IF_MPI (ljmin == jmin) THEN
        DO i=MAX0(limin,ig(jmin)+1),MIN0(limax,id(jmin)-1)
          h0(i,ljmin+1)=h0(i,ljmin+2)
          ssh(i,ljmin+1)=ssh(i,ljmin+2)
        ENDDO
     ENDIF_MPI

     IF_MPI (limax==imax) THEN

        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
#ifdef key_castest_estuary
          h0(limax,j)=-5.00_rsh
#else
          h0(limax,j)=h0(limax-1,j)
          ssh(limax,j)=ssh(limax-1,j)
#endif
        ENDDO
     ENDIF_MPI

     IF_MPI (limin == imin) THEN
        DO j=MAX0(ljmin,jb(imin)+1),MIN0(ljmax,jh(imin)-1)
          h0(limin+1,j)=h0(limin+2,j)
          ssh(limin+1,j)=ssh(limin+2,j)
        ENDDO
     ENDIF_MPI

     ! les quatre coins 
     IF_MPI (limax == imax .AND. ljmax == jmax) THEN
        h0(limax,ljmax)=h0(limax-1,ljmax-1)
        ssh(limax,ljmax)=ssh(limax-1,ljmax-1)
     ENDIF_MPI

     IF_MPI (limax==imax .AND. ljmin==jmin) THEN
        h0(limax,ljmin+1)=h0(limax-1,ljmin+2)
        ssh(limax,ljmin+1)=ssh(limax-1,ljmin+2)
     ENDIF_MPI

     IF_MPI (limin==imin .AND. ljmin==jmin) THEN
        h0(limin+1,ljmin+1)=h0(limin+2,ljmin+2)
        ssh(limin+1,ljmin+1)=ssh(limin+2,ljmin+2)
     ENDIF_MPI

     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
        h0(limin+1,ljmax)=h0(limin+2,ljmax-1)
        ssh(limin+1,ljmax)=ssh(limin+2,ljmax-1)
     ENDIF_MPI

#ifdef key_toseelater
!  PROBLEME : h0 , ssh ne sont pas alloues de imin a imxc et de jmin a jmax
!             il faut les echanger avec les processeurs extremes et les ranger dans m2 ou p2
!             voir comme dans tsobc3dapply par exemple
     ! condition de periodicite dans la direction x
     IF(l_obc_cycl_x) THEN
        DO j=ljminp1,ljmax
          h0(imin+1,j)=h0(imax-1,j)
          h0(imax,j)=h0(imin+2,j)
          ssh(imin+1,j)=ssh(imax-1,j)
          ssh(imax,j)=ssh(imin+2,j)
        ENDDO
     ENDIF
     ! condition de periodicite dans la direction y
     IF(l_obc_cycl_y) THEN
        DO i=liminp1,limax
          h0(i,jmin+1)=h0(i,jmax-1)
          h0(i,jmax)=h0(i,jmin+2)
          ssh(i,jmin+1)=ssh(i,jmax-1)
          ssh(i,jmax)=ssh(i,jmin+2)
        ENDDO
     ENDIF
#endif

!ATTENTION : hupper et dhcell non utilises, a integrer a l avenir pour bien gerer 
!            le decouvrement (coherence avec fwet de SIAM) - voir PLH

#ifdef key_toseelater
      hultra=10000.0_rsh
      DO j=ljmin,ljmax
        DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
         IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN

          hupper(i,j)=MIN(hx(i,j)*morphox(i,j)+(1.0_rsh-morphox(i,j))*hultra,         &
                          hx(i-1,j)*morphox(i-1,j)+(1.0_rsh-morphox(i-1,j))*hultra,   &
                          hy(i,j)*morphoy(i,j)+(1.0_rsh-morphoy(i,j))*hultra,         &
                          hy(i,j-1)*morphoy(i,j-1)+(1.0_rsh-morphoy(i,j-1))*hultra)
          IF(hupper(i,j).GE.hultra)hupper(i,j)=-99.99_rsh
          dhcell(i,j)=h0(i,j)-hupper(i,j)
         ENDIF
        ENDDO
      ENDDO
#endif

   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,h0(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,h0(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hx(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hx(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hy(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hy(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))

  END SUBROUTINE sed_morpho      

   !!==============================================================================
#endif

END MODULE sedimorpho
