source: trunk/SOURCES/Heino_files/dragging_heino_mod.f90

!> \file dragging_heino_mod.f90
!! Module qui defini la loi de glissement dans les experiences Heino
!<

!> \namespace sliding_heino
!! Definition de la loi de glissement dans les experiences Heino
!! \author Catherine
!! \date ...
!! @note Used module
!! @note   - use module3D_phy
!<

module sliding_heino

! Défini la loi de glissement dans les expériences Heino

USE module3D_phy


implicit none


integer,dimension(nx,ny) :: mksedim      !< masque des regions sediment-hard rock
integer,dimension(nx,ny) :: mkxsedim     !< décalle sur les demi mailles
integer,dimension(nx,ny) :: mkysedim     !< décalle sur les demi mailles
real :: Cr                            !< coefficient glissement pour hard rock
real :: Cs                            !< coefficient glissement pour sediment
real :: coefmax
real :: coefslid
real :: coefdrag
logical,dimension(nx,ny) :: gzmx_heino     !< pour transporter la valeur dans flottab
logical,dimension(nx,ny) :: gzmy_heino
real :: longslope



contains
!-------------------------------------------------------------------------------
!> SUBROUTINE: init_sliding
!! Cette routine fait l'initialisation du glissement et est appelée juste apres
!! la lecture du masque
!<
subroutine init_sliding


if (itracebug.eq.1)  call tracebug(' Heino: entree dans routine init_sliding')

! Cette routine fait l'initialisation du glissement et est appelée juste apres
! la lecture du masque


!    ecriture dans le fichier parametres
     write(num_rep_42,*)'glissement Heino -  initialisation masque'
     write(num_rep_42,*)'loi glissement=',loigliss
     write(num_rep_42,*)'----------------------------------------------------------'


! Le masque est égal au mk0 lu du début, 
     
mksedim(:,:)=mk0(:,:)
where (mk0>1) mk0=1

!do j=ny, 1, -1
!    write(6,'(81i1)') (mk0(i,j), i=1,nx)
!end do
!write(6,*) mksedim(41,41)
!do j=ny, 1, -1
!    write(6,'(81i1)') (mksedim(i,j), i=1,nx)
!end do

! calcul du masque sur les demi mailles : frontiere -> sedim

do j=2,ny-1
   do i=2,nx-1
      mkxsedim(i,j)=max(mksedim(i,j),mksedim(i-1,j))
      mkysedim(i,j)=max(mksedim(i,j),mksedim(i,j-1))
   end do
end do


Cr=1.e5   ! en a-1  coefficient hard rock (loi en puissance 3)
Cs=500.    ! en a-1  coefficient sediment (lineaire)

   write(num_rep_42,*)'coefficients:   Cr=',Cr,'    Cs=',Cs,' (a-1)'
   write(num_rep_42,*)'glissement seulement si les 2 voisins sont au point de fusion'

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

      return
      end subroutine init_sliding

!________________________________________________________________________________

!>SUBROUTINE: sliding
!! Cette routine calcule le terme ddbx et ddby dans le cas Heino
!<
subroutine sliding

! cette routine calcule le terme ddbx et ddby dans le cas Heino
!      -------------------------------------------------
!      GLISSEMENT
!      DDBX et DDBY termes dus au glissement 
!      relation avec la vitesse de glissement UXB et UYB 
!      UXB=-DDBX*SDX et UYB=-DDBY*SDY
!      -------------------------------------------------


if (itracebug.eq.1)  call tracebug(' Heino: entree dans routine sliding')

coefmax=1.

slidy: do j=2,NY
   do I=2,NX-1
pose_y:  if (.not.flgzmy(i,j)) then

   if (HMY(I,J).lt.1.) then        ! glace tres fine
      ddby(i,j)=0.

!   else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).and.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then 
!      ddby(i,j)=0.                 ! pas de glissement seulement
!                                   ! si les 2 voisins sont froids

   else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then 
      ddby(i,j)=0.                 ! pas de glissement si au moins un seul 
                                   ! des 2 voisins est froid

   else                            ! glissement
      coefslid=(hwater(i,j)+hwater(i,j-1))*0.5

      coefslid=max(coefslid,0.)
      coefslid=min(coefslid/coefmax,1.)
      coefslid=1.

      if (mkysedim(i,j).eq.1) then  ! loi hard rock
         ddby(i,j)=Cr*slope2my(i,j)*Hmy(i,j)*coefslid

      else if (mkysedim(i,j).eq.2) then  ! loi sediment
         ddby(i,j)=Cs*Hmy(i,j)*coefslid

!         if (abs(sdy(i,j)).gt.1.e-8) then
!            j_1=max(j-1,1)
!            j1=min(j+1,ny)
!            longslope=(sdy(i,j_1)+sdy(i,j)+sdy(i,j1))/3.
!            ddby(i,j)=ddby(i,j)/sdy(i,j)*longslope
!         endif

      else
         ddby(i,j)=0.
      endif

   endif


endif pose_y
end do
end do slidy


slidx: do j=2,NY-1
   do I=2,NX
pose_x:  if (.not.flgzmx(i,j)) then

   if (HMx(I,J).lt.1.) then        ! glace tres fine
      ddbx(i,j)=0.

!   else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).and.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then 
!      ddbx(i,j)=0.                 ! pas de glissement seulement
!                                   ! si les 2 voisins sont froids

   else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then 
      ddbx(i,j)=0.                 ! pas de glissement si au moins un seul 
                                   ! des 2 voisins est froid



   else                            ! glissement

      coefslid=(hwater(i,j)+hwater(i-1,j))*0.5

      coefslid=max(coefslid,0.)
      coefslid=min(coefslid/coefmax,1.)
      coefslid=1.

      if (mkxsedim(i,j).eq.1) then  ! loi hard rock
         ddbx(i,j)=Cr*slope2mx(i,j)*Hmx(i,j)*coefslid

      else if (mkxsedim(i,j).eq.2) then  ! loi sediment
         ddbx(i,j)=Cs*Hmx(i,j)*coefslid

!         if (abs(sdx(i,j)).gt.1.e-8) then
!            i_1=max(i-1,1)
!            i1=min(i+1,nx)
!            longslope=(sdx(i_1,j)+sdx(i,j)+sdx(i1,j))/3.
!            ddbx(i,j)=ddbx(i,j)/sdx(i,j)*longslope
!         endif


      else
         ddbx(i,j)=0.
      endif

   endif

endif pose_x
end do
end do slidx

return
end subroutine  sliding

!-------------------------------------------------------------------------

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

subroutine dragging   

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

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

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

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 
            ! base froide
         else
            betamy(:,:)=coefdrag
            flgzmy(i,j)=.true.
         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 
            ! base froide
         else
            betamx(:,:)=coefdrag
            flgzmx(i,j)=.true.
         endif

      endif
   end do
end do frotx

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

return
end subroutine dragging



END MODULE  sliding_heino