MODULE sed_MUSTANG_MARS

#if defined key_sedim_MUSTANG && defined key_MARS

   !&E==========================================================================
   !&E                   ***  MODULE  sed_MUSTANG_MARS  ***
   !&E
   !&E
   !&E ** Purpose : concerns subroutines related to sediment dynamics link to hydrodynamic model
   !&E              to be used in MARS system
   !&E 
   !&E ** Description :
   !&E     subroutine sedinit_fromfile  ! reads filrepsed where all results of 
   !&E                                   ! sediment dyn. are stored (depends on hydro model)
   !&E     subroutine sed_gradvit         ! calcul gradient de vitesse, u*
   !&E     subroutine sed_obc_corflu      ! echange MPI des corrections de flux horizontaux pour sables
   !&E     subroutine sed_skinstress     ! computes the skin stress
   !&E     subroutine sed_exchange_w2s_MARS  ! echange MPI de flux verticaux a l'interface WS
   !&E     subroutine sed_exchange_s2w_MARS  ! echange MPI de flux verticaux a l'interface WS
   !&E     subroutine sed_outsaverestart ! save related field for future runs
   !&E
   !&E                    *********************************************************************
   !&E                    *****         ATTENTION                                           ***
   !&E                    *****                                                             ***
   !&E                    *****     OMP : fast all subroutines are into parallel section    ***
   !&E                    *****                                                             ***
   !&E                    *********************************************************************
   !&E ** History :
   !&E     !  2009-03  (Pierre Le Hir, Florence Cayocca)
   !&E     !  2010-06  (Florence Cayocca) corrections
   !&E     !  2010-06  (F. Cayocca, V. Garnier) MPI
   !&E     !  2010-07  (B. Thouvenin) taking into account of non constitutive particulate variables
   !&E     !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E     !  2011-02  (B. Thouvenin) * unique dynamical allocation of raphbx, raphby, frofonx, frofony
   !&E                                * get restart file at a different date
   !&E                                * evolution of diffusion/consolidation
   !&E                                * re-estimate of ndt_part according to the real settling velocity (key_parsub_newdt)
   !&E     !  2012-11 (F. Grasso) add a boolean to use "sandbottomcell" (need to add a line in variable.dat)
   !&E     !  2012-12 (F. Grasso) modif from Y. Kervella for Non Uniform Bed Coverage
   !&E     !  2012-12 (F. Grasso) to add a layer of fluid mud (creme de vase)
   !&E     !  2013-01 (F. Grasso) change the fill_value for sed conc, sal and temp (-rg_valmanq_io)
   !&E     !  2013-01 (F. Grasso) limitation of tenfon tot (<6 N.m-2)
   !&E     !  2014    (B. Thouvenin)  Evolutions, introduction of key_sand2D, introduction of initial adjustment of the sediment height 
   !&E     !  2014-05  (F. Grasso)    Modifications on consolidation subroutines based on "Grasso et al. (Ocean Dynynamics, 2014)"
   !&E                                * subroutine sed_consolid_mixsed
   !&E                                * subroutine sed_constitutivrel
   !&E     !  2014-05  (F. Cayocca) Update and modifs some names
   !&E                                * sand/mud erosion laws
   !&E                                * morpho module
   !&E     ! 2014-05  (R. Verney)  : Modification of settling velocities laws and formulations
   !&E     ! 2014-05  (F. Cayocca)   beginning of update morphodynamic (not operationnal - TO SEE LATER) -  
   !&E     ! 2014-05  (F. Cayocca) Update and modifs some names
   !&E     ! 2014-12  (B.Thouvenin) : parallelization OMP
   !&E     ! 2015-01  (B.Thouvenin) : corrections skinstress..
   !&E     ! 2015-02  (B.Thouvenin) : add key_nosubs_Dbio_insed if using sedimento only for constitutive particules
   !&E                                (with key_biolo and key_benthos for example) 
   !&E                                (idem as key_nosedflux of Katerina Kombiadou)
   !&E     ! 2015-04  (B.Thouvenin)  settling velocity of sand variables introduce in sed_init_ws_sand  
   !&E     ! 2015-07  (F. Grasso) Mise en place dragages et contrainte de fond (vague+courant) selon Soulsby (1995)
   !&E     ! 2015-04  (B.Thouvenin  P. Le Hir) : gestion of dissolved variables in interstitial (pore) waters  
   !E      ! 2015-04  (F.Cayocca, P. Le Hir, B. Thouvenin) corrections for morphodynamic and introduction of lateral slip
   !&E     ! 2015-04  (B.Thouvenin  P. Le Hir) : gestion of dissolved variables in interstitial (pore) waters  
   !&E     ! 2015-07  (B.Thouvenin  P. Le Hir) : correction lateral erosion of dry cell  
   !&E     ! 2015-07  (F. Grasso) Mise en place dragages et contrainte de fond (vague+courant) selon Soulsby (1995)
   !&E     ! 2015-09  (B.Thouvenin  P. Le Hir) : correction consolidation+diffusion and add bioturbation  (for key_mixsed) 
   !&E     ! 2015-12  (B.Thouvenin) : extraction of routines  for reorganization of module SEDIMARS to MUSTANG
   !&E
   !&E==========================================================================

#include "toolcpp.h"
#include "coupleur_define_MUSTANG.h"

   !! * Modules used

   USE comvars2d,    ONLY : ig,id,jb,jh
   USE comsubstance, ONLY : nv_adv,nvp 
   USE comMUSTANG
   USE coupleur_MUSTANG
   
   IMPLICIT NONE

   !! * Accessibility

   ! functions & routines of this module, called outside :
   ! PUBLIC functions                                ! none ?
   PUBLIC sed_gradvit,sed_obc_corflu, sedinit_fromfile
   PUBLIC sed_skinstress,sed_outsaverestart,sed_exchange_s2w_MARS,sed_exchange_w2s_MARS

   PRIVATE

   
 CONTAINS
 
    !!==============================================================================
  SUBROUTINE sedinit_fromfile(ifirst,ilast,jfirst,jlast,rhoref,h0)
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sedinit_fromfile  ***
   !&E
   !&E ** Purpose : manages the fields to be re-initialized in case of the 
   !&E              continuation of a previous run
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  main 
   !&E
   !&E ** External calls : ionc4_openr, ionc4_read_time,
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (Pierre Le Hir) extracted from init.F90
   !&E       !  oct.2008 inserted into module_sedimento.F90 by F. Cayocca
   !&E       !  2009-02  modified for new sand/mud sediment (P. Le Hir & F. Cayocca)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   USE comsubstance, ONLY : nvp,nvpc,init_cv_name
   USE comvars2d,    ONLY : tsauv,t,saverestart_step,iscreenlog,l_restart_subs,iscreenlog
   USE parameters,   ONLY : riosh,grav,valmanq
   USE_MPI toolmpi,  ONLY : MASTER
   USE ionc4,        ONLY : ionc4_openr,ionc4_read_subxyt,ionc4_read_subzxyt,ionc4_close
   USE initMUSTANG,  ONLY : MUSTANG_init_ws_sand

   !! * Arguments
   INTEGER,INTENT(IN)                 ::   ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_m2p2),INTENT(IN)   :: h0
   REAL(KIND=rsh),INTENT(IN)          ::   rhoref

   !! * Local declarations
   CHARACTER*50     :: filesand,filemud
   CHARACTER*8      :: nomcv,nomcs
   INTEGER          :: i,j,k,iv,ivg  
   LOGICAL          :: l_readiv             
   REAL(KIND=rsh)                                             ::  somalp,dzsf,eptot,dzsmaxmin,h0max_def
   REAL(KIND=riosh),DIMENSION(ksdmax,limin:limax,ljmin:ljmax) ::  tzecrits
   REAL(KIND=riosh),DIMENSION(limin:limax,ljmin:ljmax)        ::  xind

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

     CALL MUSTANG_init_ws_sand(rhoref)

     CALL ionc4_openr(filrepsed,l_out_nc4par=.true.)

     CALL ionc4_read_subxyt(filrepsed,'ksmi',xind,limin,limax,ljmin,ljmax,1,0,0)
     DO j=ljmin,ljmax
        DO i=limin,limax
          ksmi(i,j)=INT(xind(i,j))
        END DO
     END DO

     CALL ionc4_read_subxyt(filrepsed,'ksma',xind,limin,limax,ljmin,ljmax,1,0,0)
     DO j=ljmin,ljmax 
        DO i=limin,limax 
          ksma(i,j)=INT(xind(i,j)) 
        END DO 
     END DO
     WHERE (h0(limin:limax,ljmin:ljmax) == -valmanq) ksmi(limin:limax,ljmin:ljmax)=1
     WHERE (h0(limin:limax,ljmin:ljmax) == -valmanq) ksma(limin:limax,ljmin:ljmax)=0

     IF(.NOT. l_dzsmaxuni) THEN 
       IF( dzsmaxuni==0.0_rsh) THEN ! dzsmax directement lu dans le fihcier initial
         CALL ionc4_read_subxyt(filrepsed,'dzsmax',xind,limin,limax,ljmin,ljmax,1,0,0)
         DO j=ljmin,ljmax 
           DO i=limin,limax 
             dzsmax(i,j)=INT(xind(i,j)) 
           END DO 
         END DO
       ELSE
         ! definition des dzsmax
         ! estime en fonction de la profondeur et de dzsmaxuni considere comme la valeur max de dzsmax
         ! ici valeur minimum de dzsmax en profondeur=dzsmaxuni/100
         ! profondeur max donnant la plus faible valeur de dzsma a informer, mise par defaut=50m
          dzsmaxmin=dzsmaxuni/100._rsh
          h0max_def=50.0_rsh
          write(iscreenlog,*)
          write(iscreenlog,*)'DZSMAX NON UNIFORM and not known (dzsmaxuni/=0) in file :',filrepsed
          write(iscreenlog,*)'by default linearly computed in MUSTANG_sedinitfromfile (sedim.F90)'
          write(iscreenlog,*)'from dzsmaxuni (',dzsmaxuni,'m) to dzsmaxuni/100 (',dzsmaxmin,'m) depending on water depth'
          write(iscreenlog,*)'dzsmax minimum for depth > ',h0max_def
          write(iscreenlog,*)'PROGRAM your own initialization of dzsmax (or READ a file) in MUSTANG_sedinitfromfile'
          dzsmax(limin:limax,ljmin:ljmax)=(h0max_def-h0(limin:limax,ljmin:ljmax))/h0max_def*dzsmaxuni+dzsmaxmin
          WHERE (dzsmax(limin:limax,ljmin:ljmax) > dzsmaxuni ) dzsmax(limin:limax,ljmin:ljmax)=dzsmaxuni
          WHERE (dzsmax(limin:limax,ljmin:ljmax) < dzsmaxmin ) dzsmax(limin:limax,ljmin:ljmax)=dzsmaxmin

         ! ou lus dans un fichier (a programmer)
       
       ENDIF
     ELSE
         dzsmax(limin:limax,ljmin:ljmax)=dzsmaxuni
     ENDIF
       

     CALL ionc4_read_subzxyt(filrepsed,'dzs',tzecrits,limin,limax,ljmin,ljmax,1,ksdmax,1,0,0,1)
     DO j=ljmin,ljmax 
      DO i=limin,limax 
        DO k = 1,ksdmax 
           dzs(k,i,j)=tzecrits(k,i,j)
        END DO 
      END DO 
     END DO

     DO iv =-1,nv_adv
       l_readiv=.false.
       IF (iv==-1) THEN
          nomcv='concTEM'
          nomcs='fxw2sTEM'
          l_readiv=.true.
       ELSE IF (iv==0) THEN
          nomcv='concSAL'
          nomcs='fxw2sSAL'
          l_readiv=.true.
       ELSE 
          WRITE(nomcv,10) iv
          WRITE(nomcs,20) iv
          IF (l_restart_subs .OR. iv <= nvpc .OR. (.NOT. l_restart_subs .AND. TRIM(init_cv_name(iv))/='none')) l_readiv=.true.
       ENDIF

       IF (l_readiv) THEN
         CALL ionc4_read_subzxyt(filrepsed,nomcv,tzecrits,                         &
                           limin,limax,ljmin,ljmax,1,ksdmax,1,0,0,1)
         DO j =ljmin,ljmax
           DO i =limin,limax
             DO k =1,ksdmax
               cv_sed(iv,k,i,j)=tzecrits(k,i,j)
             ENDDO
           ENDDO
         ENDDO

         CALL ionc4_read_subxyt(filrepsed,nomcs,xind,limin,limax,ljmin,ljmax,1,0,0)
         DO j = ljmin,ljmax
           DO i = limin,limax
             flx_w2s(iv,i,j) = xind(i,j)
           ENDDO
         ENDDO
      ELSE
         IF(l_initsed_vardiss) THEN
          !initialisation des variables dissoutes + salinite et temperature avec cv_wat du fond
          DO j = ljmin,ljmax
          DO i = limin,limax
           IF(iv==-1)THEN
             DO k=ksmi(i,j),ksma(i,j)
               cv_sed(iv,k,i,j)=temp_bottom_MUSTANG(i,j)
             ENDDO
           ELSE IF(iv==0) THEN
             DO k=ksmi(i,j),ksma(i,j)
               cv_sed(iv,k,i,j)=sal_bottom_MUSTANG(i,j)
             ENDDO
           ELSE IF(iv>nvp) THEN
             DO k=ksmi(i,j),ksma(i,j)
               cv_sed(iv,k,i,j)=cw_bottom_MUSTANG(iv,i,j)
             ENDDO
           ELSE
             cv_sed(iv,:,i,j)=cini_sed(iv)
           ENDIF
          ENDDO
          ENDDO         
         ELSE
           cv_sed(iv,:,:,:)=cini_sed(iv)
         ENDIF
 
      ENDIF

      ENDDO
      CALL ionc4_close(filrepsed)


     DO j=ljmin,ljmax
       DO i=MAX0(limin,ig(j)),MIN0(limax,id(j))
         z0sed(i,j)=z0seduni
         IF (h0(i,j) /= -valmanq) THEN
         DO k=ksmi(i,j),ksma(i,j)
           somalp=0.0_rsh
           c_sedtot(k,i,j)=0.0_rsh
           DO iv=1,nvpc
             c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
             somalp=somalp+cv_sed(iv,k,i,j)/ros(iv)
           END DO
           poro(k,i,j)=1.0_rsh-somalp
         END DO
        END IF
       END DO
     END DO
      
   10 FORMAT('concen',i2.2)
   20 FORMAT('flxw2s',i2.2)
      
  END SUBROUTINE sedinit_fromfile
!!===========================================================================================

  SUBROUTINE sed_gradvit(ifirst,ilast,jfirst,jlast,h0fond,rhoref)

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_gradvit  ***
   !&E
   !&E ** Purpose : calcul du taux de cisaillement  depend du code ?
   !&E                 ATTENTION ON EST EN ZONE PARALLEL OMP
   !&E
   !&E ** Description : to be programmed using hydrodynamic knowledge
   !&E           using htot, RESIDUAL_THICKNESS_WAT (h0fond), RHOREF, sig, epn, nz ..
   !&E     output : gradvit (in comMUSTANG)
   !&E
   !&E ** Called by :  step
   !&E
   !&E ** Reference : 
   !&E
   !&E ** History :
    !&E       !  2015-12  (B.Thouvenin) reorganization of module SEDIMARS
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
#ifdef key_MARS
   USE comvars3d,    ONLY : sig,nz,turb_nbeq, turb_2eq_option
   USE comturb,      ONLY : eps
#endif
   USE comMUSTANG,     ONLY: ustarbot,gradvit

   !! * Arguments (Private variables for OMP because locales variables)
   INTEGER, INTENT(IN)                        :: ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),INTENT(IN)                  :: h0fond,rhoref

   !! * Local declarations
   INTEGER                            :: i,j,k


!$OMP DO SCHEDULE(RUNTIME)
      DO j=jfirst,jlast
#ifdef key_MARS
      DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
#else
      DO i=ifirst,ilast
#endif

          IF(htot(i,j).GT.h0fond)  THEN
#ifdef key_MARS
            DO k=1,kmax
              IF (turb_nbeq == 2 .AND. turb_2eq_option == 2) THEN ! using k-epsilon closure scheme
                gradvit(k,i,j)=sqrt(eps(k,i,j)/nz(k,i,j))
              ELSE ! using nezu and nakagawa
                gradvit(k,i,j)=sqrt(ustarbot(i,j)**3._rsh/(0.4_rsh)*(-sig(k)/(1._rsh+sig(k))))
              ENDIF
            ENDDO
#else
           !  calcul specifique pour autre modele ??
#endif

          ENDIF
#ifdef key_MARS
        ENDIF
#endif
      ENDDO
      ENDDO
!$OMP END DO

  END SUBROUTINE sed_gradvit

   !!==============================================================================


  SUBROUTINE sed_obc_corflu
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_obc_corflu ***
   !&E
   !&E ** Purpose : 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : sed_MUSTANG_first
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (Pierre Le Hir) extracted from traceur.F90
   !&E       !  2008-10  (F. Cayocca)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   USE comvars2d,      ONLY : l_obc_cycl_x,l_obc_cycl_y
   USE parameters,     ONLY : liminm2,liminm1,limaxp1,limaxp2, &
                              ljminm2,ljminm1,ljmaxp1,ljmaxp2, &
                              imin,imax,jmin,jmax
   USE comsubstance,   ONLY : isand1,isand2
   USE_MPI toolmpi,    ONLY : ex_i_rsh,ex_j_rsh


   !! * Arguments

   !! * Local declarations
    INTEGER                              :: i,j,ivp


   !!--------------------------------------------------------------------------
   !! * Executable part
   
!$OMP SINGLE
     ! hypothese gradient nul aux frontieres 
     IF_MPI (ljmax == jmax) THEN
#ifdef key_MARS 
        DO i=MAX0(limin,ig(jmax)+1),MIN0(limax,id(jmax)-1)
#else
        DO i=limin,limax
#endif
          DO ivp=isand1,isand2
             corflux(ivp,i,jmax)=corflux(ivp,i,jmax-1)
             corfluy(ivp,i,jmax)=corfluy(ivp,i,jmax-1)
          ENDDO
        ENDDO
     ENDIF_MPI
     IF_MPI (ljmin == jmin) THEN
#ifdef key_MARS 
        DO i=MAX0(limin,ig(jmin)+1),MIN0(limax,id(jmin)-1)
#else
        DO i=limin,limax
#endif
          DO ivp=isand1,isand2
             corflux(ivp,i,jmin+1)=corflux(ivp,i,jmin+2)
             corfluy(ivp,i,jmin+1)=corfluy(ivp,i,jmin+2)
          ENDDO
        ENDDO
     ENDIF_MPI
     IF_MPI (limin==imax) THEN
#ifdef key_MARS 
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
#else
        DO j=ljmin,ljmax
#endif
          DO ivp=isand1,isand2
             corflux(ivp,imax,j)=corflux(ivp,imax-1,j)
             corfluy(ivp,imax,j)=corfluy(ivp,imax-1,j)
          ENDDO
        ENDDO
     ENDIF_MPI
     IF_MPI (limin == imin) THEN
#ifdef key_MARS 
        DO j=MAX0(ljmin,jb(imin)+1),MIN0(ljmax,jh(imin)-1)
#else
        DO j=ljmin,ljmax
#endif
          DO ivp=isand1,isand2
             corflux(ivp,imin+1,j)=corflux(ivp,imin+2,j)
             corfluy(ivp,imin+1,j)=corfluy(ivp,imin+2,j)
          ENDDO
        ENDDO
     ENDIF_MPI

     ! les quatre coins
     IF_MPI (limax == imax .AND. ljmax == jmax) THEN
          DO ivp=isand1,isand2
             corflux(ivp,imax,jmax)=corflux(ivp,imax-1,jmax-1)
             corfluy(ivp,imax,jmax)=corfluy(ivp,imax-1,jmax-1)
          ENDDO
     ENDIF_MPI
     IF_MPI (limax==imax .AND. ljmin==jmin) THEN
          DO ivp=isand1,isand2
             corflux(ivp,imax,jmin+1)=corflux(ivp,imax-1,jmin+2)
             corfluy(ivp,imax,jmin+1)=corfluy(ivp,imax-1,jmin+2)
          ENDDO
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmin==jmin) THEN
          DO ivp=isand1,isand2
             corflux(ivp,imin+1,jmin+1)=corflux(ivp,imin+2,jmin+2)
             corfluy(ivp,imin+1,jmin+1)=corfluy(ivp,imin+2,jmin+2)
          ENDDO
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
          DO ivp=isand1,isand2
             corflux(ivp,imin+1,jmax)=corflux(ivp,imin+2,jmax-1)
             corfluy(ivp,imin+1,jmax)=corfluy(ivp,imin+2,jmax-1)
          ENDDO
     ENDIF_MPI

#ifdef key_toseelater
!  PROBLEME : corflux, corfluy se 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
          DO ivp=isand1,isand2
             corflux(ivp,imin+1,j)=corflux(ivp,imax-1,j)
             corfluy(ivp,imin+1,j)=corfluy(ivp,imax-1,j)
             corflux(ivp,imax,j)=corflux(ivp,imin+2,j)
             corfluy(ivp,imax,j)=corfluy(ivp,imin+2,j)
          ENDDO
        ENDDO
     ENDIF
     ! condition de periodicite dans la direction y
     IF(l_obc_cycl_y) THEN
        DO i=liminp1,limax
          DO ivp=isand1,isand2
             corflux(ivp,i,jmin+1)=corflux(ivp,i,jmax-1)
             corfluy(ivp,i,jmin+1)=corfluy(ivp,i,jmax-1)
             corflux(ivp,i,jmax)=corflux(ivp,i,jmin+2)
             corfluy(ivp,i,jmax)=corfluy(ivp,i,jmin+2)
          ENDDO
        ENDDO
     ENDIF
#endif

!$OMP END SINGLE

OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(0,1,nv_adv,limin,limaxp1,ljmin,ljmax  ,corflux(:,limin:limaxp1,ljmin:ljmax  ))
   CALL_MPI ex_j_rsh(0,1,nv_adv,limin,limax  ,ljmin,ljmaxp1,corfluy(:,limin:limax  ,ljmin:ljmaxp1))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(corflux,corfluy)

!$OMP SINGLE

      ! corflux are interpolated on mesh edges (in u & v)
      ! -------------------------------------------------
    DO j=ljminp1,ljmaxm1
#ifdef key_MARS 
        DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)
#else
        DO i=limin,limax
#endif
          DO ivp=isand1,isand2
            corflux(ivp,i,j)=0.5_rsh*(corflux(ivp,i,j)+corflux(ivp,i+1,j))
            corfluy(ivp,i,j)=0.5_rsh*(corfluy(ivp,i,j)+corfluy(ivp,i,j+1))
           ENDDO
     ENDDO
   ENDDO

     ! hypothese gradient nul aux frontieres 
     ! corflux et corfluy sont ainsi calcules pour liminp1 -> limax et ljminp1 ->ljmax 
    IF_MPI (ljmax == jmax) THEN
#ifdef key_MARS 
        DO i=MAX0(limin,ig(jmax)+1),MIN0(limax,id(jmax)-1)
#else
        DO i=limin,limax
#endif
          DO ivp=isand1,isand2
             corflux(ivp,i,jmax)=corflux(ivp,i,jmax-1)
             corfluy(ivp,i,jmax)=corfluy(ivp,i,jmax-1)
          ENDDO
        ENDDO
     ENDIF_MPI
     IF_MPI (ljmin == jmin) THEN
#ifdef key_MARS 
        DO i=MAX0(limin,ig(jmin)+1),MIN0(limax,id(jmin)-1)
#else
        DO i=limin,limax
#endif
          DO ivp=isand1,isand2
             corflux(ivp,i,jmin+1)=corflux(ivp,i,jmin+2)
             corfluy(ivp,i,jmin+1)=corfluy(ivp,i,jmin+2)
          ENDDO
        ENDDO
     ENDIF_MPI
     IF_MPI (limax==imax) THEN
        DO j=limin,ljmax
          DO ivp=isand1,isand2
             corflux(ivp,imax,j)=corflux(ivp,imax-1,j)
             corfluy(ivp,imax,j)=corfluy(ivp,imax-1,j)
          ENDDO
        ENDDO
     ENDIF_MPI
     IF_MPI (limin == imin) THEN
        DO j=ljmin,ljmax
          DO ivp=isand1,isand2
             corflux(ivp,imin+1,j)=corflux(ivp,imin+2,j)
             corfluy(ivp,imin+1,j)=corfluy(ivp,imin+2,j)
          ENDDO
        ENDDO
     ENDIF_MPI

     ! les quatre coins 
     IF_MPI (limax == imax .AND. ljmax == jmax) THEN
          DO ivp=isand1,isand2
             corflux(ivp,imax,jmax)=corflux(ivp,imax-1,jmax-1)
             corfluy(ivp,imax,jmax)=corfluy(ivp,imax-1,jmax-1)
          ENDDO
     ENDIF_MPI
     IF_MPI (limax==imax .AND. ljmin==jmin) THEN
          DO ivp=isand1,isand2
             corflux(ivp,imax,jmin+1)=corflux(ivp,imax-1,jmin+2)
             corfluy(ivp,imax,jmin+1)=corfluy(ivp,imax-1,jmin+2)
          ENDDO
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmin==jmin) THEN
          DO ivp=isand1,isand2
             corflux(ivp,imin+1,jmin+1)=corflux(ivp,imin+2,jmin+2)
             corfluy(ivp,imin+1,jmin+1)=corfluy(ivp,imin+2,jmin+2)
          ENDDO
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
          DO ivp=isand1,isand2
             corflux(ivp,imin+1,jmax)=corflux(ivp,imin+2,jmax-1)
             corfluy(ivp,imin+1,jmax)=corfluy(ivp,imin+2,jmax-1)
          ENDDO
     ENDIF_MPI

#ifdef key_toseelater
!  PROBLEME : corflux et corfluy 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
          DO ivp=isand1,isand2
             corflux(ivp,imin+1,j)=corflux(ivp,imax-1,j)
             corfluy(ivp,imin+1,j)=corfluy(ivp,imax-1,j)
             corflux(ivp,imax,j)=corflux(ivp,imin+2,j)
             corfluy(ivp,imax,j)=corfluy(ivp,imin+2,j)
          ENDDO
        ENDDO
     ENDIF
     ! condition de periodicite dans la direction y 
     IF(l_obc_cycl_y) THEN
        DO i=liminp1,limax
          DO ivp=isand1,isand2
             corflux(ivp,i,jmin+1)=corflux(ivp,i,jmax-1)
             corfluy(ivp,i,jmin+1)=corfluy(ivp,i,jmax-1)
             corflux(ivp,i,jmax)=corflux(ivp,i,jmin+2)
             corfluy(ivp,i,jmax)=corfluy(ivp,i,jmin+2)
          ENDDO
        ENDDO
     ENDIF
#endif

!$OMP END SINGLE

   PRINT_DBG*, 'FIN SED_OBC_corflu'

   
   END SUBROUTINE sed_obc_corflu

   !!==============================================================================
   SUBROUTINE sed_skinstress
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_skinstress  ***
   !&E
   !&E ** Purpose : computes current related bottom shear stress
   !&E
   !&E ** Description : to be programmed from hydrodynamic variables
   !&E                   rhoref, boz, hx, hy, ssh, h0fond, u,uz, v, vz, hminfrot,hmaxfrot,
   !&E         using some parameters and variables of MUSTANG :
   !&E                    z0sed,  htncrit
   !&E               and Wave coupling
   !&E
   !&E
   !&E ** Called by :  sed_MUSTANG_first
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (Pierre Le Hir) extracted from init.F90
   !&E       !  2008-10  (F. Cayocca) inserted into module_sedimento.F90
   !&E       ! 2011-11 (J. Vareilles) Change in the expression of the bottom stress to be compatible with the 2D/3D
   !&E                              distinctions as is done in the hydrodynamic part
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
   USE parameters,     ONLY : grav,pi,exp1,epsilon,valmanq
   USE comvars2d,      ONLY : hx,hy,l_obc_cycl_x,l_obc_cycl_y,hex,hey,htx,hty,hm,hminfrot,hmaxfrot,nbouc,h0fond,h0
   USE comvarp2d,      ONLY : ssh,u,v
   USE comvarp3d,      ONLY : uz,vz,boz
   USE comvars3d,      ONLY : rhoref,sig,dsigu,roswat_bot
#ifdef key_wave
   USE wave,           ONLY : wave_ubot, wave_tp,l_wave_rdir,wave_dir
#endif
#if defined key_siggen || defined key_gencoord
   USE toolgeom,       ONLY : f_hzu_x,f_hzu_y
   USE comsiggen,      ONLY : hzutx,hzuty
#endif
#if defined key_1DV_force_towfile
   USE readforcingk1DV, ONLY : mytauw_out      
#endif

   !! * Arguments

      
   USE_MPI toolmpi,    ONLY : ex_i_rsh,ex_j_rsh

   !! * Local declarations
   INTEGER i,j,k
   REAL(KIND=rsh)                                   :: sigfrot,vit,z,thoule,hhoule,fws2ij
   REAL(KIND=rsh)                                   :: alofwscd,tentot,xwa,coef_a,coef_b,       &
                                                       coef_m,coef_n,coef_p,coef_q,htn,alpha,   &
                                                       beta,cosamb,sinamb
      
   !REAL(KIND=rsh),PARAMETER                         :: zfrot=0.010_rsh
   REAL(KIND=rsh),PARAMETER                         :: zfrot=0.50_rsh
   
   REAL(KIND=rsh),PARAMETER                         :: a1=-0.07_rsh, a2=1.87_rsh, a3=-0.34_rsh, a4=-0.12_rsh,   &
                                                       b1=0.27_rsh, b2=0.51_rsh, b3=-0.1_rsh, b4=-0.24_rsh,     &
                                                       cm1=0.72_rsh, cm2=-0.33_rsh, cm3=0.08_rsh, cm4=0.34_rsh, &
                                                       cn1=0.78_rsh, cn2=-0.23_rsh, cn3=0.12_rsh, cn4=-0.12_rsh,&
                                                       p1=-0.75_rsh, p2=0.13_rsh, p3=0.12_rsh, p4=0.02_rsh,     &
                                                       q1=0.89_rsh, q2=0.4_rsh, q3=0.5_rsh, q4=-0.28_rsh,       &     
                                                       hminwave=0.05_rsh

   REAL(KIND=rsh),PARAMETER                         :: a1pa2=a1+a2,a3pa4=a3+a4,                 &
                                                       b1pb2=b1+b2,b3pb4=b3+b4,                 &
                                                       cm1pm2=cm1+cm2,cm3pm4=cm3+cm4,           &
                                                       cn1pn2=cn1+cn2,cn3pn4=cn3+cn4,           &
                                                       p1pp2=p1+p2,p3pp4=p3+p4,                 &
                                                       q1pq2=q1+q2,q3pq4=q3+q4                 
#ifdef key_wave
   REAL(KIND=rsh)                                   :: courant
#endif
   !!--------------------------------------------------------------------------
   !! * Executable part

!ATTENTION deja dans une region parallel OMP ouverte dans step.F90
!!!$OMP DO SCHEDULE (DYNAMIC,5)
!   DO j=ljmin,ljmax
!   DO i=limin,limax
!      rho_bott(i,j)=rhoref*(1-boz(1,i,j)/grav)
!!   cdcurrent=0.16_rsh/(LOG(1.0_rsh/z0seduni))**2
!      raphbx(i,j)=0.0_rsh
!      raphby(i,j)=0.0_rsh
!      frofonx(i,j)=0.0_rsh
!      frofony(i,j)=0.0_rsh
!   ENDDO
!   ENDDO
!!!$OMP END DO
!$OMP SINGLE
      !rho_bott(limin:limax,ljmin:ljmax)=rhoref*(1-boz(1,limin:limax,ljmin:ljmax)/grav)
      raphbx(limin:limax,ljmin:ljmax)=0.0_rsh
      raphby(limin:limax,ljmin:ljmax)=0.0_rsh
      frofonx(limin:limax,ljmin:ljmax)=0.0_rsh
      frofony(limin:limax,ljmin:ljmax)=0.0_rsh
!$OMP END SINGLE
! correction janvier 2015
!htx,hty n ont pas ete reactualises apres le calcul de ssh au temps n+1
! on les reactualise ici
#ifdef key_MPI_2D
!$OMP DO SCHEDULE(RUNTIME) PRIVATE(i,j)
#else
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(i,j)
#endif
   DO j=ljmin,ljmaxm1
     DO i=limin,limaxm1
       htx(i,j)=hx(i,j)+0.5_rsh*(ssh(i,j)+ssh(i+1,j))
       hty(i,j)=hy(i,j)+0.5_rsh*(ssh(i,j)+ssh(i,j+1))
#if defined key_siggen || defined key_gencoord
       hzutx(1,i,j)=f_hzu_x(hx(i,j),0.5_rsh*(ssh(i,j)+ssh(i+1,j)),1,i,j)
       hzuty(1,i,j)=f_hzu_y(hy(i,j),0.5_rsh*(ssh(i,j)+ssh(i,j+1)),1,i,j)
#endif
       IF ( htx(i,j)<=0.0_rsh .AND. ( hx(i,j)+ssh(i,j) > 0.0_rsh .OR. hx(i,j)+ssh(i+1,j) > 0.0_rsh) ) THEN
         htx(i,j)=hx(i,j)+MAX(ssh(i,j),ssh(i+1,j))
#if defined key_siggen || defined key_gencoord
         hzutx(1,i,j)=f_hzu_x(hx(i,j),MAX(ssh(i,j),ssh(i+1,j)),1,i,j)
#endif
       END IF
       IF ( hty(i,j)<=0.0_rsh .AND. ( hy(i,j)+ssh(i,j) > 0.0_rsh .OR. hy(i,j)+ssh(i,j+1) > 0.0_rsh )) THEN
         hty(i,j)=hy(i,j)+MAX(ssh(i,j),ssh(i,j+1))
#if defined key_siggen || defined key_gencoord
         hzuty(1,i,j)=f_hzu_y(hy(i,j),MAX(ssh(i,j),ssh(i,j+1)),1,i,j)
#endif
       ENDIF
     END DO

     IF_MPI (limax==imax) THEN
     IF (hy(imax,j) /= -valmanq-h0fond) THEN
       hty(imax,j)=hy(imax,j)+0.5_rsh*(ssh(imax,j)+ssh(imax,j+1))
#if defined key_siggen || defined key_gencoord
       hzuty(1,imax,j)=f_hzu_y(hy(imax,j),0.5_rsh*(ssh(imax,j)+ssh(imax,j+1)),1,imax,j)
#endif
       IF ( hty(imax,j)<=0.0_rsh .AND. ( hy(imax,j)+ssh(imax,j) > 0.0_rsh .OR. hy(imax,j)+ssh(imax,j+1) > 0.0_rsh )) THEN
           hty(imax,j)=hy(imax,j)+MAX(ssh(imax,j),ssh(imax,j+1))
#if defined key_siggen || defined key_gencoord
           hzuty(1,imax,j)=f_hzu_y(hy(imax,j),MAX(ssh(imax,j),ssh(imax,j+1)),1,imax,j)
#endif
       ENDIF
     END IF
     ENDIF_MPI
   END DO ! end do for j loop
!$OMP END DO
!$OMP SINGLE
  IF_MPI (ljmax==jmax) THEN
     DO i=limin,limaxm1
       IF (hx(i,jmax) /= -valmanq-h0fond) THEN
         htx(i,jmax)=hx(i,jmax)+0.5_rsh*(ssh(i,jmax)+ssh(i+1,jmax))
#if defined key_siggen || defined key_gencoord
         hzutx(1,i,jmax)=f_hzu_x(hx(i,jmax),0.5_rsh*(ssh(i,jmax)+ssh(i+1,jmax)),1,i,jmax)
#endif
         IF ( htx(i,jmax)<=0.0_rsh .AND. ( hx(i,jmax)+ssh(i,jmax) > 0.0_rsh .OR. hx(i,jmax)+ssh(i+1,jmax) > 0.0_rsh) ) THEN
             htx(i,jmax)=hx(i,jmax)+MAX(ssh(i,jmax),ssh(i+1,jmax))
#if defined key_siggen || defined key_gencoord
             hzutx(1,i,jmax)=f_hzu_x(hx(i,jmax),MAX(ssh(i,jmax),ssh(i+1,jmax)),1,i,jmax)
#endif
         END IF
       END IF
     ENDDO
   ENDIF_MPI
!$OMP END SINGLE

! modification de PLH pour le cas estuar, mais pourquoi ?
!IF(l_morphocoupl) THEN
!!$OMP DO SCHEDULE(RUNTIME)
!   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
!       IF (htx(i,j).GE.htncrit_eros)raphbx(i,j)=1.0_rsh
!       IF (hty(i,j).GE.htncrit_eros)raphby(i,j)=1.0_rsh
!      ENDIF
!     ENDDO
!   ENDDO
!!!$OMP END DO

!ELSE

!$OMP DO SCHEDULE(RUNTIME) PRIVATE(i,j)
   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
        IF (htx(i,j).LE.hm) THEN
           raphbx(i,j)=ABS(u(i,j))/(ABS(u(i,j))+epsilon)
        ELSE
           raphbx(i,j)=ABS(uz(1,i,j))/(ABS(uz(1,i,j))+epsilon)
        ENDIF

        IF(hty(i,j).LE.hm) THEN
           raphby(i,j)=ABS(v(i,j))/(ABS(v(i,j))+epsilon)
        ELSE
           raphby(i,j)=ABS(vz(1,i,j))/(ABS(vz(1,i,j))+epsilon)
        ENDIF
      ENDIF
     ENDDO
   ENDDO
!$OMP END DO
!ENDIF

   ! the bottom stress needs to be computed at a fairly "uniform" distance from the bed
   ! throughout the domain. This distance is called zfrot (=0.5), and transcribed into the
   ! corresponding sigma layer
   ! Otherwise, the bottom stress is simply computed in the first layer 

!FLO : rajouter un booleen dans parasedim pour choisir l option
   
!$OMP DO SCHEDULE(RUNTIME) PRIVATE(i,j,vit,k,z)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       IF(htx(i,j).LE.hm) THEN            !under hm, hydrodynamic 2D calculation
          vit=SQRT(u(i,j)**2+0.0625_rsh*(v(i,j)+v(i+1,j)            &
              +v(i,j-1)+v(i+1,j-1))**2)*u(i,j)
          frofonx(i,j)=0.16_rsh*(LOG(MAX(MIN(htx(i,j),hmaxfrot),hminfrot)/(z0sed(i,j)*exp1)))**(-2)*vit
       ELSE                               !above hm, hydrodynamic 3D calculation
          k=1
#if defined key_siggen || defined key_gencoord
          z=0.5_rsh*hzutx(k,i,j)*dsigu(k)
#else
          z=(1.0_rsh+sig(k))*htx(i,j)
#endif
          vit=SQRT(uz(k,i,j)**2+0.0625_rsh*(vz(k,i,j)+vz(k,i+1,j)            &
               +vz(k,i,j-1)+vz(k,i+1,j-1))**2)*uz(k,i,j)
          frofonx(i,j)=0.16_rsh*(LOG(z/z0sed(i,j)+1.0_rsh))**(-2)*vit
       ENDIF
       
       IF(hty(i,j).LE.hm) THEN            !under hm, hydrodynamic 2D calculation
          vit=SQRT(0.0625_rsh*(u(i,j)+u(i,j+1)+u(i-1,j+1)         &
                  +u(i-1,j))**2+v(i,j)**2)*v(i,j)
          frofony(i,j)=0.16_rsh*(LOG(MAX(MIN(hty(i,j),hmaxfrot),hminfrot)/(z0sed(i,j)*exp1)))**(-2)*vit
       ELSE                               !above hm, hydrodynamic 3D calculation
          k=1
#if defined key_siggen || defined key_gencoord
          z=0.5_rsh*hzuty(k,i,j)*dsigu(k)
#else
          z=(1.0_rsh+sig(k))*hty(i,j)
#endif
          vit=SQRT(0.0625_rsh*(uz(k,i,j)+uz(k,i,j+1)+uz(k,i-1,j+1)         &
                           +uz(k,i-1,j))**2+vz(k,i,j)**2)*vz(k,i,j)
          frofony(i,j)=0.16_rsh*(LOG(z/z0sed(i,j)+1.0_rsh))**(-2)*vit
       ENDIF
     ENDDO
   ENDDO
!$OMP END DO


OMPMPI barrier
OMPMPI master
   CALL_MPI ex_i_rsh(-1,0,1,liminm1,limax,ljmin  ,ljmax,raphbx (liminm1:limax,ljmin  :ljmax))
   CALL_MPI ex_i_rsh(-1,0,1,liminm1,limax,ljmin  ,ljmax,frofonx(liminm1:limax,ljmin  :ljmax))
   CALL_MPI ex_j_rsh(-1,0,1,limin  ,limax,ljminm1,ljmax,raphby (limin  :limax,ljminm1:ljmax))
   CALL_MPI ex_j_rsh(-1,0,1,limin  ,limax,ljminm1,ljmax,frofony(limin  :limax,ljminm1:ljmax))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(raphbx,frofonx,raphby,frofony)



   ! calcul de la tension (en N/m2)
    ! ------------------------------
!$OMP SINGLE
   tenfonc(limin:limax,ljmin:ljmax)=0.0_rsh
   tenfonw(limin:limax,ljmin:ljmax)=0.0_rsh
   tenfon(limin:limax,ljmin:ljmax)=0.0_rsh
!$OMP END SINGLE

!$OMP DO SCHEDULE(RUNTIME) PRIVATE(i,j,htn)
   DO j=ljminp2,ljmax
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j))
       IF (h0(i,j) /= -valmanq) THEN
         htn=h0(i,j)+ssh(i,j)

         IF(htn .LT.htncrit_eros+h0fond .AND. dry_cell(i,j)==0) THEN
           dry_cell(i,j)=1
         ENDIF
         IF(htn .GT.2.0_rsh*htncrit_eros+h0fond .AND. dry_cell(i,j)==1) THEN
           dry_cell(i,j)=0
         ENDIF

         IF(dry_cell(i,j)==0) THEN
           tenfonc(i,j)=SQRT(((frofonx(i,j)*raphbx(i,j)+frofonx(i-1,j)        &
                     *raphbx(i-1,j))/(raphbx(i,j)                             &
                     +raphbx(i-1,j)+epsilon))**2+((frofony(i,j)*raphby(i,j)   &
                     +frofony(i,j-1)*raphby(i,j-1))/(raphby(i,j)              &
                     +raphby(i,j-1)+epsilon))**2)*roswat_bot(i,j)
            tenfonc(i,j)=min(10.0_rsh,tenfonc(i,j))  !limitation sur tenfonc,FG(21/07/2015)

#ifdef key_wave
            fws2ij=fws2
            IF (l_fricwave .AND. wave_ubot(i,j)*wave_tp(i,j) > 0.0_rsh .AND.  &
                     wave_tp(i,j) > 0.001_rsh.AND.wave_ubot(i,j)>0.001_rsh) THEN
              fws2ij=0.5_rsh*1.39_rsh*(wave_ubot(i,j)*wave_tp(i,j)/REAL(2.0_rlg*pi*z0sed(i,j),rsh))**(-0.52_rsh)
            ENDIF

            ! calcul de tenfonw (vagues)
            ! --------------------------
!FLO 11jan2011           tenfonw(i,j)=(rho_bott(i,j)*fws2*wave_ubot(i,j)**2)
            tenfonw(i,j)=(1025._rsh*fws2ij*wave_ubot(i,j)**2)
            !tenfon(i,j)=SQRT(tenfonw(i,j)**2+tenfonc(i,j)**2)
            tenfonw(i,j)=min(15.0_rsh,tenfonw(i,j))  !limitation sur tenfonw, FG(21/07/2015)

            IF (l_wave_rdir ) THEN
               courant=SQRT(u(i,j)**2+v(i,j)**2)
               IF(tenfonc(i,j) > 0.0_rsh .AND. tenfonw(i,j) > 0.0_rsh .AND. courant> 0.0_rsh) THEN

                ! calcul de tenfon avec la formulation de Soulsby (1995)
                ! ======================================================
                ! calcul de tenfonc (courant) influence par les vagues
                ! ----------------------------------------------------
                  tenfonc(i,j)=tenfonc(i,j)*(1+(1.2*(tenfonw(i,j)/(tenfonw(i,j)+tenfonc(i,j)))**3.2))
                  !tenfonc(i,j)=min(10.0_rsh,tenfonc(i,j))  !limitation sur tenfonc, FG(21/07/2015)

                ! calcul de la difference d angle entre la direction des vagues et du courant
                ! ---------------------------------------------------------------------------
                ! calcul de l orientation du courant par rapport au nord
                  alpha=ACOS(v(i,j)/courant)   ! en radians
                ! calcul de l orientation des vagues par rapport au nord
                  beta=wave_dir(i,j)*pi/180   ! en radians
                ! calcul de cos(alpha-beta) et sin(alpha-beta)
                  cosamb=ABS(COS(alpha-beta))
                  sinamb=ABS(SIN(alpha-beta))

                ! calcul de tenfon (vagues + courant)
                ! -----------------------------------
                  tenfon(i,j)=SQRT((tenfonc(i,j)+tenfonw(i,j)*cosamb)**2 +(tenfonw(i,j)*sinamb)**2)

              ELSE
                 tenfon(i,j)=tenfonw(i,j)+tenfonc(i,j)
              ENDIF
            ELSE

            ! calcul de tenfon avec la moyenne arithmetique de tw,tc (maximise tenfon, max=25N/m2)
            ! =======================================================================
               tenfon(i,j)=tenfonw(i,j)+tenfonc(i,j)

            ! calcul de tenfon avec la moyenne quadratique de tw,tc (minimise tenfon, max=18N/m2)
            ! =======================================================================
            !   tenfon(i,j)=SQRT(tenfonw(i,j)**2+tenfonc(i,j)**2)
            ENDIF
#elif key_1DV_force_towfile
            tenfon(i,j)=mytauw_out
#else
            tenfon(i,j)=tenfonc(i,j)
#endif
         ENDIF
       ELSE
         tenfonc(i,j)=-valmanq
         tenfonw(i,j)=-valmanq
         tenfon(i,j)=-valmanq
       ENDIF
     ENDDO
   ENDDO
!$OMP END DO

!$OMP SINGLE
  ! modification des conditions aux frontieres (28 fevrier 2012)
     ! hypothese gradient nul aux frontieres 
     IF_MPI (ljmax == jmax) THEN
        DO i=MAX0(limin,ig(jmax)+1),MIN0(limax,id(jmax)-1)
          tenfon(i,jmax)=tenfon(i,jmax-1)
        ENDDO
     ENDIF_MPI
     IF_MPI (ljmin == jmin) THEN
        DO i=MAX0(limin,ig(jmin)+1),MIN0(limax,id(jmin)-1)
          tenfon(i,jmin+1)=tenfon(i,jmin+2)
        ENDDO
     ENDIF_MPI
     IF_MPI (limax==imax) THEN
        DO j=MAX0(ljmin,jb(imax)+1),MIN0(ljmax,jh(imax)-1)
           tenfon(imax,j)=tenfon(imax-1,j)
        ENDDO
     ENDIF_MPI
     IF_MPI (limin == imin) THEN
        DO j=MAX0(ljmin,jb(imin)+1),MIN0(ljmax,jh(imin)-1)
          tenfon(imin+1,j)=tenfon(imin+2,j)
        ENDDO
     ENDIF_MPI

     ! les quatre coins (28 fevrier 2012)
     IF_MPI (limax == imax .AND. ljmax == jmax) THEN
        tenfon(imax,jmax)=tenfon(imax-1,jmax-1)
     ENDIF_MPI
     IF_MPI (limax==imax .AND. ljmin==jmin) THEN
        tenfon(imax,jmin+1)=tenfon(imax-1,jmin+2)
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmin==jmin) THEN
        tenfon(imin+1,jmin+1)=tenfon(imin+2,jmin+2)
     ENDIF_MPI
     IF_MPI (limin==imin .AND. ljmax==jmax) THEN
        tenfon(imin+1,jmax)=tenfon(imin+2,jmax-1)
     ENDIF_MPI
     
#ifdef key_toseelater
!  PROBLEME : tenfon 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 dans la direction x (28 fevrier 2012)
     IF(l_obc_cycl_x) THEN
        DO j=ljminp1,ljmax
          tenfon(imin+1,j)=tenfon(imax-1,j)
          tenfon(imax,j)=tenfon(imin+2,j)
        ENDDO
     ENDIF
     ! condition de periodicite dans la direction y (9 mai 2012)
     IF(l_obc_cycl_y) THEN
        DO i=liminp1,limax
          tenfon(i,jmin+1)=tenfon(i,jmax-1)
          tenfon(i,jmax)=tenfon(i,jmin+2)
        ENDDO
     ENDIF
#endif
     
   PRINT_DBG*,'controle tenfon max: ',MAXVAL(tenfonc),MAXVAL(tenfonw),MAXVAL((tenfon))
   PRINT_DBG*,'control tenfon max: ',MAXVAL(tenfonc),MAXVAL(polation des corrections de flux pour les sables dans les coins de la grille, roswat_bot)
   PRINT_DBG*,'controle tenfon max: ',MAXLOC(tenfonc),MAXLOC(roswat_bot)
   PRINT_DBG*,'controle tenfon max: ',MAXVAL(tenfon)
   PRINT_DBG*, 'FIN SED_SKINSTRESS'  
!$OMP END SINGLE


  END SUBROUTINE sed_skinstress
 
  !!==============================================================================

  SUBROUTINE sed_exchange_w2s_MARS

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_exchange_MPI  ***
   !&E
   !&E ** Purpose :  
   !&E                 ATTENTION ON EST EN ZONE PARALLEL OMP
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  sed_MUSTANG_deposition
   !&E
   !&E ** Reference : 
   !&E
   !&E ** History :
    !&E       !  2015-12  (B.Thouvenin) reorganization of module SEDIMARS=MUSTANG
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_j_rsh
   USE parameters,   ONLY : liminm1,limin,limax,limaxp1,ljminm1,ljmin,ljmax,ljmaxp1
 
   !! * Arguments


  !! * Local declarations

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


! slip depo
! echange MPI
    
OMPMPI barrier
OMPMPI master
    CALL_MPI ex_i_rsh(0,1,nvp ,liminm1,limaxp1,ljminm1,ljmaxp1,flx_w2s_corim1(:,liminm1:limaxp1  ,ljminm1:ljmaxp1 ))
    CALL_MPI ex_i_rsh(-1,0,nvp,liminm1,limaxp1,ljminm1,ljmaxp1,flx_w2s_corip1(:,liminm1:limaxp1  ,ljminm1:ljmaxp1 ))
    CALL_MPI ex_j_rsh(0,1,nvp,liminm1,limaxp1,ljminm1,ljmaxp1,flx_w2s_corjm1(:,liminm1:limaxp1  ,ljminm1:ljmaxp1 ))
    CALL_MPI ex_j_rsh(-1,0,nvp,liminm1,limaxp1,ljminm1,ljmaxp1,flx_w2s_corjp1(:,liminm1:limaxp1  ,ljminm1:ljmaxp1 ))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(flx_w2s_corim1,flx_w2s_corip1,flx_w2s_corjm1,flx_w2s_corjp1)


  END SUBROUTINE sed_exchange_w2s_MARS

  !!==============================================================================

  SUBROUTINE sed_exchange_s2w_MARS

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_exchange_MPI  ***
   !&E
   !&E ** Purpose :  
   !&E                 ATTENTION ON EST EN ZONE PARALLEL OMP
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  sed_MUSTANG_first
   !&E
   !&E ** Reference : 
   !&E
   !&E ** History :
    !&E       !  2015-12  (B.Thouvenin) reorganization of module SEDIMARS=MUSTANG
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
   USE_MPI toolmpi,  ONLY : ex_i_rsh,ex_j_rsh
   USE parameters,   ONLY : liminm1,limin,limax,limaxp1,ljminm1,ljmin,ljmax,ljmaxp1
 
   !! * Arguments


  !! * Local declarations

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


! lateral erosion of dry cell
! echange MPI
    
OMPMPI barrier
OMPMPI master
    CALL_MPI ex_i_rsh(0,1,nv_adv+2 ,liminm1,limaxp1,ljminm1,ljmaxp1,flx_s2w_corim1(:,liminm1:limaxp1  ,ljminm1:ljmaxp1 ))
    CALL_MPI ex_i_rsh(-1,0,nv_adv+2,liminm1,limaxp1,ljminm1,ljmaxp1,flx_s2w_corip1(:,liminm1:limaxp1  ,ljminm1:ljmaxp1 ))
    CALL_MPI ex_j_rsh(0,1,nv_adv+2,liminm1,limaxp1,ljminm1,ljmaxp1,flx_s2w_corjm1(:,liminm1:limaxp1  ,ljminm1:ljmaxp1 ))
    CALL_MPI ex_j_rsh(-1,0,nv_adv+2,liminm1,limaxp1,ljminm1,ljmaxp1,flx_s2w_corjp1(:,liminm1:limaxp1  ,ljminm1:ljmaxp1 ))
OMPMPI end master
OMPMPI barrier
OMPMPI flush(flx_s2w_corim1,flx_s2w_corip1,flx_s2w_corjm1,flx_s2w_corjp1)


  END SUBROUTINE sed_exchange_s2w_MARS

 !!===========================================================================
 
  SUBROUTINE sed_outsaverestart

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_outsaverestart  ***
   !&E
   !&E ** Purpose : saves sediment-related field for future runs
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  step 
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E
   !&E       !  2008-10  (F. Cayocca) inserted into module_sedimento.F90
   !&E       !  2009-02  (P. Le Hir) modified declarations
   !
   !objet    :  gestion de la reprise de calcul a partir d un etat
         !           enregistre dans un fichier (netcdf ou binaire acces direct)
         !           lors d un run anterieur
   !
   !     appelant         : stepu
         !                :
   !     appeles          : ionc4_createfile, ionc4_torigine, 
   !     appeles          : ionc4_createvar, ionc4_createvar
   !     appeles          : ionc4_createvar, ionc4_write_time
   !     appeles          : ionc4_write_xy, ionc4_write_xyt, ionc4_write_zxyt
         !                :
   !     canal d ecriture :
         !                :
   !     canal de lecture : 26, filrep
         !                :
   !     cle cpp          : key_notra
         !                :
   !     variables in     : t, u, uz, xu, v, vz, xv, sal, temp, tz
   !     variables in     : ect, boz, nz, kz, wz, vish_xe, vish_phi
   !     variables in     : filsauvsed, lon2d,lat2d, ntra
   !     variables in     : cdate, h0,sig
         !                :
   !     variables out    : 
         !                :
   !     variables inout  : tsauv
         !                :
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
   !!   ------------
   USE parameters,   ONLY :  riosh,rg_valmanq_io,lchain
   USE comvars2d,    ONLY :  lon2d_g,lat2d_g,lon2d,lat2d,l_saverestart_bydate,date_ref, &
                             tsauv,cdate,saverestart_step,name_out_h0,name_out_xe,t
   USE_MPI toolmpi,  ONLY : ADD_TAG,MASTER
   USE output,       ONLY : l_out_nc4par
   USE ionc4,       ONLY : ionc4_openr,ionc4_read_subxyt,ionc4_read_subzxyt,ionc4_close,  &
                           ionc4_createfile,ionc4_gatt,ionc4_createvar,ionc4_write_sig,&
                           ionc4_vatt_fill,ionc4_torigine,ionc4_write_time,ionc4_write_xyt,&
                           ionc4_write_zxyt,ionc4_sync
                           
   !! * Local declarations
   !!   ------------------

   CHARACTER(LEN=8)      :: nomc
   CHARACTER(LEN=13)     :: prefsauv
   CHARACTER(LEN=lchain) :: chain,filsauvsed
   INTEGER               :: i,j,k,l,iv,ifois
   REAL(KIND=riosh),DIMENSION(limin:limax,ljmin:ljmax)          :: xind
   REAL(KIND=riosh),DIMENSION(ksdmin:ksdmax,limin:limax,ljmin:ljmax)   :: tlits
   REAL(KIND=riosh),DIMENSION(0:ksdmax+1)                       :: sigsed

   DATA ifois /0/
   
   !! * Executable part
   !!   ---------------

   PRINT_DBG*, 'enter sed_outsaverestart'

   ! initialisation of sigsed
   DO k=ksdmin,ksdmax
      sigsed(k)=k
   END DO
   sigsed(0)=0
   sigsed(ksdmax+1)=0
   
   ! creation
   ! ********

   IF (ifois.EQ.0 .AND. .NOT.(l_saverestart_bydate)) THEN
     ifois=1
     filsauvsed='./save_sedim.nc'
     IF_MPI(.NOT. l_out_nc4par) CALL ADD_TAG(filsauvsed,LEN(filsauvsed))
     IF ( l_out_nc4par ) THEN
        CALL ionc4_createfile(filsauvsed,lon2d_g,lat2d_g,imin,imax,1,          &
                              jmin,jmax,1,ksdmin,ksdmax,1,ksdmax,0,l_out_nc4par = .TRUE.)
     ELSE
        CALL ionc4_createfile(filsauvsed,lon2d,lat2d,limin,limax,1,          &
                              ljmin,ljmax,1,ksdmin,ksdmax,1,ksdmax,0)
     ENDIF
    write(*,*)filsauvsed

   ! save attribute of the spatial extension of the domain
   ! -----------------------------------------------------
     CALL ionc4_gatt(filsauvsed,'global_imin',imin)
     CALL ionc4_gatt(filsauvsed,'global_imax',imax)
     CALL ionc4_gatt(filsauvsed,'global_jmin',jmin)
     CALL ionc4_gatt(filsauvsed,'global_jmax',jmax)

   ! indispensable pour concatenation pas terrible a revoir
     CALL ionc4_createvar(filsauvsed,'SIG',' ','sigma variable',&
#if defined key_siggen || defined key_gencoord
                          standard_name= 'ocean_s_variable',&
#else
                          standard_name= 'ocean_sigma_variable',&
#endif
                          l_auxcoordinate=.false.,&
                          valid_min= -1.0_riosh,valid_max= 0.0_riosh,fill_value= rg_valmanq_io,&
                          dims= 'z', l_out_nc4par = l_out_nc4par)
     CALL ionc4_write_sig(filsauvsed,'SIG',sigsed,ksdmax)

     CALL ionc4_torigine(filsauvsed,date_ref)
     CALL ionc4_createvar(filsauvsed,'ksmi','s.o','ksmi',fill_value=rg_valmanq_io, &
                          dims='xyt',valid_min=0.0_riosh,valid_max=10.0_riosh,l_out_nc4par=l_out_nc4par)
     CALL ionc4_createvar(filsauvsed,'ksma','s.o','ksma',fill_value=rg_valmanq_io, &
                          dims='xyt',valid_min=0.0_riosh,valid_max=1000.0_riosh,l_out_nc4par=l_out_nc4par)
     CALL ionc4_createvar(filsauvsed,'dzsmax','s.o','ksma',fill_value=rg_valmanq_io, &
                          dims='xyt',valid_min=0.0_riosh,valid_max=1.0_riosh,l_out_nc4par=l_out_nc4par)
     CALL ionc4_createvar(filsauvsed,'dzs','m','dzs',fill_value=rg_valmanq_io, &
                          dims='zxyt',valid_min=0.0_riosh,valid_max=100.0_riosh,l_out_nc4par=l_out_nc4par)
     DO iv=-1,nv_adv     ! pour T et S
       IF (iv>0) WRITE(nomc,10) iv
       IF (iv==-1) nomc='concTEM'
       IF (iv==0) nomc='concSAL'
       CALL ionc4_createvar(filsauvsed,nomc,'u.s.i','concentration ',&
                            fill_value=rg_valmanq_io,dims='zxyt',l_out_nc4par=l_out_nc4par)
       IF (iv>0) WRITE(nomc,20) iv
       IF (iv==-1) nomc='fxw2sTEM'
       IF (iv==0) nomc='fxw2sSAL'
       CALL ionc4_createvar(filsauvsed,nomc,'u.s.i','flx_w2s ',&
                            fill_value=rg_valmanq_io,dims='xyt',l_out_nc4par=l_out_nc4par)
     END DO 
   ENDIF

   IF (l_saverestart_bydate.AND.t.GE.tsauv) THEN
     prefsauv='./save_sedim_'
     filsauvsed=prefsauv//cdate(7:10)//cdate(4:5)//  &
                    cdate(1:2)//cdate(12:13)//cdate(15:16)//cdate(18:19)//'.nc'
     IF_MPI(.NOT. l_out_nc4par) CALL ADD_TAG(filsauvsed,LEN(filsauvsed))
     IF ( l_out_nc4par ) THEN
        CALL ionc4_createfile(filsauvsed,lon2d_g,lat2d_g,imin,imax,1,          &
                              jmin,jmax,1,ksdmin,ksdmax,1,ksdmax,0,l_out_nc4par = .TRUE.)
     ELSE
        CALL ionc4_createfile(filsauvsed,lon2d,lat2d,limin,limax,1,          &
                              ljmin,ljmax,1,ksdmin,ksdmax,1,ksdmax,0)
     ENDIF

   ! save attribute of the spatial extension of the domain
   ! -----------------------------------------------------
     CALL ionc4_gatt(filsauvsed,'global_imin',imin)
     CALL ionc4_gatt(filsauvsed,'global_imax',imax)
     CALL ionc4_gatt(filsauvsed,'global_jmin',jmin)
     CALL ionc4_gatt(filsauvsed,'global_jmax',jmax)

   ! indispensable pour concatenation pas terrible a revoir
     CALL ionc4_createvar(filsauvsed,'SIG',' ','sigma variable',&
#if defined key_siggen || defined key_gencoord
                          standard_name= 'ocean_s_variable',&
#else
                          standard_name= 'ocean_sigma_variable',&
#endif
                          l_auxcoordinate=.false.,&
                          valid_min= -1.0_riosh,valid_max= 0.0_riosh,fill_value= rg_valmanq_io,&
                          dims= 'z', l_out_nc4par = l_out_nc4par)
     CALL ionc4_write_sig(filsauvsed,'SIG',sigsed,ksdmax)

     CALL ionc4_torigine(filsauvsed,date_ref)
     CALL ionc4_createvar(filsauvsed,'ksmi','s.o','ksmi',fill_value=rg_valmanq_io, &
                          dims='xyt',valid_min= 0.0_riosh,valid_max= 10.0_riosh,l_out_nc4par=l_out_nc4par)
     CALL ionc4_createvar(filsauvsed,'ksma','s.o','ksma',fill_value=rg_valmanq_io, &
                          dims='xyt',valid_min= 0.0_riosh,valid_max= 1000.0_riosh,l_out_nc4par=l_out_nc4par)
     CALL ionc4_createvar(filsauvsed,'dzsmax','s.o','ksma',fill_value=rg_valmanq_io, &
                          dims='xyt',valid_min=0.0_riosh,valid_max=1.0_riosh,l_out_nc4par=l_out_nc4par)
     CALL ionc4_createvar(filsauvsed,'dzs','m','dzs',fill_value=rg_valmanq_io, &
                          dims='zxyt',valid_min= 0.0_riosh,valid_max= 100.0_riosh,l_out_nc4par=l_out_nc4par)
     DO iv=-1,nv_adv     ! pour T et S
       IF (iv>0) WRITE(nomc,10) iv
       IF (iv==-1) nomc='concTEM'
       IF (iv==0) nomc='concSAL'
       CALL ionc4_createvar(filsauvsed,nomc,'u.s.i','concentration ',&
                            fill_value=rg_valmanq_io,dims='zxyt',l_out_nc4par=l_out_nc4par)
       IF (iv>0) WRITE(nomc,20) iv
       IF (iv==-1) nomc='fxw2sTEM'
       IF (iv==0) nomc='fxw2sSAL'
       CALL ionc4_createvar(filsauvsed,nomc,'u.s.i','flx_w2s ',&
                            fill_value=rg_valmanq_io,dims='xyt',l_out_nc4par=l_out_nc4par)
     END DO 
   ENDIF

   ! ecriture
   ! ********

   IF (t.GE.tsauv) THEN
     CALL ionc4_write_time(filsauvsed,1,t)
     DO j=ljmin,ljmax
       DO i=limin,limax
         xind(i,j)=REAL(ksmi(i,j))
       END DO
     END DO
     CALL ionc4_write_xyt(filsauvsed,'ksmi',xind,limin,limax,ljmin,ljmax,1,l_out_nc4par=l_out_nc4par)

     DO j=ljmin,ljmax
       DO i=limin,limax
         xind(i,j)=REAL(ksma(i,j))
       END DO
     END DO
     CALL ionc4_write_xyt(filsauvsed,'ksma',xind,limin,limax,ljmin,ljmax,1,l_out_nc4par=l_out_nc4par)
 
     DO j=ljmin,ljmax
       DO i=limin,limax
         xind(i,j)=REAL(dzsmax(i,j))
       END DO
     END DO
     CALL ionc4_write_xyt(filsauvsed,'dzsmax',xind,limin,limax,ljmin,ljmax,1,l_out_nc4par=l_out_nc4par)
 
     DO j = ljmin,ljmax
       DO i = limin,limax
         DO k = 1,ksdmax
           tlits(k,i,j) = dzs(k,i,j)
         END DO
       END DO
     END DO
     CALL ionc4_write_zxyt(filsauvsed,'dzs',tlits,limin,limax,ljmin,ljmax,ksdmin,ksdmax,1,l_out_nc4par=l_out_nc4par)

     DO iv=-1,nv_adv     ! pour T et S
       IF (iv>0) WRITE(nomc,10) iv
       IF (iv==-1) nomc='concTEM'
       IF (iv==0) nomc='concSAL'
       DO j = ljmin,ljmax
         DO i = limin,limax
           DO k = 1,ksdmax
             tlits(k,i,j) = cv_sed(iv,k,i,j)
           ENDDO
         ENDDO
       ENDDO
       CALL ionc4_write_zxyt(filsauvsed,nomc,tlits,limin,limax,ljmin,ljmax,ksdmin,ksdmax,1,l_out_nc4par=l_out_nc4par)
       IF (iv>0) WRITE(nomc,20) iv
       IF (iv==-1) nomc='fxw2sTEM'
       IF (iv==0) nomc='fxw2sSAL'
       DO j = ljmin,ljmax
         DO i = limin,limax
           xind(i,j) = flx_w2s(iv,i,j)
         ENDDO
       ENDDO
       CALL ionc4_write_xyt(filsauvsed,nomc,xind,limin,limax,ljmin,ljmax,1,l_out_nc4par=l_out_nc4par)       
     ENDDO
     tsauv=t+saverestart_step

     IF (l_saverestart_bydate) THEN
       CALL ionc4_close(filsauvsed)
     ELSE
       ! To write the data on the disk and not loose data in case of run crash
       CALL ionc4_sync(filsauvsed)
     ENDIF

   ENDIF

   10 FORMAT('concen',i2.2)
   20 FORMAT('flxw2s',i2.2)

   PRINT_DBG*, 'exit sed_outsaverestart'

  END SUBROUTINE sed_outsaverestart

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

END MODULE sed_MUSTANG_MARS
