!> \file dragging_hudson_mod.f90
!! Module qui definie la loi de glissement dans les expériences Heino
!<

!> \namespace  dragging_hudson
!! Definie la loi de glissement dans les expériences Heino
!! \author ...
!! \date ...
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use sedim_declar
!<



module dragging_hudson

! Defini la loi de glissement dans les expériences Heino
! modification par rapport a la version de Catherine : ice-shelves + streams autorises

USE module3D_phy
USE sedim_declar

implicit none

real :: Cfs                      !< coefficient glissement pour sediment
real :: coefmax
real :: coefslid
real :: coefdrag
real :: seuil_sedim              !< seuil sur hw_mx pour avoir du glissement
real, dimension(nx,ny) :: hw_mx
real, dimension(nx,ny) :: hw_my

!character(len=80) :: filin


contains

!> SUBROUTINE: init_dragging
!! Initialisation du dragging basal 
!<

  subroutine init_dragging

    implicit none

! formats pour les ecritures dans 42
428 format(A)

    namelist/drag_hudson/seuil_sedim,coefmax,Cfs 

! lecture des parametres du run                block drag_hudson
!--------------------------------------------------------------------

    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
    read(num_param,drag_hudson)
    write(num_rep_42,428)'!___________________________________________________________' 
    write(num_rep_42,428) '&drag_hudson                 ! nom du bloc drag_hudson'
    write(num_rep_42,*)
    write(num_rep_42,*) 'seuil_sedim    = ', seuil_sedim 
    write(num_rep_42,*) 'coefmax        = ', coefmax
    write(num_rep_42,*) 'Cfs            = ', Cfs
    write(num_rep_42,*)'/'
    write(num_rep_42,428) '! seuil_sedim : seuil pour passer en stream si zone sediment'
    write(num_rep_42,428) '! coefmax : hauteur d eau max dans le sediment (=hwatermax)'
    write(num_rep_42,428) '! Cfs : coefficient de la loi de frottement sediment'


!coefmax=hwatermax
!seuil_sedim=0.0001

! fichier mask fourni par l'intercomparaison Heino 
! on refait la lecture 

!      filin='mask-shelf-40.dat'

!      filin = TRIM(DIRNAMEINP)//TRIM(filin)

!       write(42,*) 'type de socle '
!       write(42,*) 'mksedim donne par ',filin

!open(11,file=filin)
!do i=1,6
!      read(11,*)  ! 6 lignes de commentaires
!end do


! lecture proprement dite de la carte sediment
!do j=ny, 1, -1
!    read(11,'(87i1)') (mksedim(i,j), i=1,nx)
!end do


!close(11)

! test pour une calotte axi-symetrique avec glissement Cr partout
!write(num_rep_42,*) 'masque sediment mis a 1 '
!mksedim(:,:)=min(mksedim(:,:),1)

! a la frontiere le demi noeud est sediment.
! calcul du masque sur les demi mailles : frontiere -> sedim
!mkxsedim(:,:)=max(mksedim(:,:),eoshift(mksedim(:,:),shift=-1,boundary=0,dim=1))
!mkysedim(:,:)=max(mksedim(:,:),eoshift(mksedim(:,:),shift=-1,boundary=0,dim=2))
! write(42,*) 'sur les bords zone sediment, masque exterieur '

!if ((loigliss.eq.5).and.(Cs.gt.1.)) then
    coefdrag=rog/Cfs
    Cfs=0.
    write(42,*)'dragging Heino'
    write(42,*)'passage en stream si zone sediment + fusion'
    write(42,*)'coefficient de frottement dans la zone sediment (Pa/m)=',coefdrag
!endif

  end subroutine init_dragging

!> SUBROUTINE: dragging
!! Defini le basal drag
!<

  subroutine dragging   

! dans la zone sediment : appele si loigliss=5,
! defini le basal drag


!flgzmx(:,:)=.false.
!flgzmy(:,:)=.false.

    betamx(:,:)=1.e5
    betamy(:,:)=1.e5

    shelfy=.false.



    call moy_mxmy(nx,ny,Hwater,hw_mx,hw_my)

    froty: do j=2,ny
       do i=2,nx-1

          if (mkysedim(i,j).eq.2) then  ! loi sediment

            ! seulement si au dessus du point de fusion 
!            if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then 
             if (hw_my(i,j).ge.seuil_sedim) then
                if (.not.flot(i,j)) then
                   gzmy(i,j)=.true.
                   coefslid=hw_my(i,j)
                   coefslid=max(coefslid,seuil_sedim)
                   coefslid=min(coefslid/coefmax,1.)
                   betamy(i,j)=min(coefdrag/coefslid,1.e5)

!               betamy(i,j)=coefdrag


                   flgzmy(i,j)=.true.
                   ddby(i,j)=0.
                   mslid_my(i,j)=5
                   shelfy=.true.
                endif
             endif
          endif

       end do
    end do froty


    frotx: do j=2,ny-1
       do i=2,nx

          if (mkxsedim(i,j).eq.2) then  ! loi sediment

            ! seulement si au dessus du point de fusion 
!            if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then 
             if (hw_mx(i,j).ge.seuil_sedim) then
                if (.not.flot(i,j)) then
                   gzmx(i,j)=.true.
                   coefslid=hw_mx(i,j)
                   coefslid=max(coefslid,seuil_sedim)
                   coefslid=min(coefslid/coefmax,1.)
                   betamx(i,j)=min(coefdrag/coefslid,1.e5)

!               betamx(i,j)=coefdrag

                   flgzmx(i,j)=.true.
                   ddbx(i,j)=0.
                   mslid_mx(i,j)=5
                   shelfy=.true.
                endif
             endif
            
          endif
       end do
    end do frotx

!   gzmx(:,:)=flgzmx(:,:)
!   gzmy(:,:)=flgzmy(:,:)
!   gzmx_heino(:,:)=gzmx(:,:)
!   gzmy_heino(:,:)=gzmy(:,:)


    return
  end subroutine dragging


END MODULE dragging_hudson