source: branches/GRISLIv3/SOURCES/MISMIP3D_files/dragging_mismip3d_mod.f90 @ 364

!> \file dragging_mismip3d_mod.f90
!<

!> \namespace dragging_mismip3
!!  Module qui fait le dragging de MISMIP avec un beta non lineraire
!!  tau_b = C  | U_b |**(m-1)  U_b

!! \author Cat
!! \date 11 mars 2012
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!<

module dragging_mismip3


  use module3d_phy
  use param_phy_mod
  use interface_input

  implicit none

  real, dimension(nx,ny)   :: Vcol_x              !< uniquement pour compatibilite dans le spinup
  real, dimension(nx,ny)   :: Vcol_y              !< uniquement pour compatibilite
  real, dimension(nx,ny)   :: Vsl_x               !<
  real, dimension(nx,ny)   :: Vsl_y               !<
  logical :: corr_def = .False.                   !< for deformation correction (compatibility)


real :: beta_limgz          !< when beta gt beta_limgz -> not gzmx

!  equation glissement        tau_b = C  | U_b |**(m-1)  U_b

real :: Coef_C              ! coefficient du glissement dans equation (C)
real :: exp_m               ! exposant du glissement 
real :: m_1                 ! exposant dans l'equation
real :: lim_unorm           ! vitesse correspondante a beta_limgz
integer :: i_dome           ! indice du noeud milieu en x
integer :: j_dome           ! indice du noeud milieu en y
 
contains
!-------------------------------------------------------------------------------

!> subroutine: init_dragging
!! Cette routine fait l'initialisation du dragging.
!<

  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.

    implicit none

    if (itracebug.eq.1)  call tracebug(' entree dans init_dragging de MISMIP')
! les coefficients sont codes en dur

    Coef_C = 1.e7                         ! Coefficient glissement en Pa m-1/3 s-1/3
    exp_m  = 1./3.                        ! exposant du glissement (m)
    Coef_C = secyear**(-exp_m) * Coef_C   ! en Pa m-1/3 a-1/3
    m_1    = exp_m - 1                    ! exposant dans le calcul de beta


!   limites sur beta
    beta_limgz = 1.e6                     ! valeur maxi de beta <-->  U = 0
    betamax    = beta_limgz
    lim_unorm  = (beta_limgz/Coef_C)**(1./m_1)  ! vitesse correspondante a beta_limgz (de l'ordre de 5 mm/an)


! coefficient permettant de modifier le basal drag.
! drag est celui fixe au depart qui reflete les proprietes du socle
    drag_mx(:,:)=1.
    drag_my(:,:)=1.

! mstream est celui courant, eventuellement affecte par les conditions basales
    mstream_mx(:,:)=1.
    mstream_my(:,:)=1.

! il faut quand meme au moins un point immobile 
! Ca fait 4 points pour des raisons de symetrie

    i_dome = (nx-1)/2 +1
    j_dome = (ny-1)/2 +1

    drag_mx(i_dome:i_dome+1,j_dome         )  = 0
    drag_my(i_dome        ,j_dome:j_dome+1 )  = 0
    mstream_mx(i_dome:i_dome+1,j_dome         )  = 0
    mstream_my(i_dome        ,j_dome:j_dome+1 )  = 0


! initialisation de certains tableaux
    ilemx(:,:)  =  .false.
    ilemy(:,:)  =  .false.


! initialisation des vitesses de glissement : lues dans lect_mismip_mod

! initialisation du beta en fonction des vitesses
! 

    call drag_non_lin_stag (nx,ny,Uxflgz,Uyflgz,betamx,betamy)

    return
  end subroutine init_dragging


!________________________________________________________________________________

!> subroutine: dragging
!! Defini la localisation des streams et le frottement basal
!! est appele par flottab
!! en entrant dans cette routine les logiques
!!  float et flotmx/flotmy sont connus 
!<
!-------------------------------------------------------------------------
  subroutine dragging   ! defini la localisation des streams et le frottement basal


    if (itracebug.eq.1)  call tracebug(' entree dans dragging de MISMIP')

    call drag_non_lin_stag (nx,ny,Uxflgz,Uyflgz,betamx,betamy)


! pas d'iles
    ilemx(:,:)  =  .false.
    ilemy(:,:)  =  .false.


    where (mstream_mx(:,:).eq.1)    ! SSA  betamx garde la valeur
       gzmx(:,:) = .true.


!       betamx(:,:) = beta_limgz     ! attention ligne pour debug qui enleve le gissement

    elsewhere
       gzmx(:,:) = .false.
       betamx(:,:) = beta_limgz
       Uxflgz(:,:) = 0.
    end where

    where (flotmx(:,:))             ! pas de frottement sous les parties flottantes
       betamx(:,:) = 0.
    end where

    where (mstream_my(:,:).eq.1)    ! SSA  betamy garde la valeur
       gzmy(:,:) = .true.


!       betamy(:,:) = beta_limgz     ! attention ligne pour debug qui enleve le gissement



    elsewhere
       gzmy(:,:) = .false.
       betamy(:,:) = beta_limgz
       Uyflgz(:,:) = 0.
    end where

    where (flotmy(:,:))             ! pas de frottement sous les parties flottantes
       betamy(:,:) = 0.
    end where
    
  if (itracebug.eq.1)  call tracebug(' fin  dragging de MISMIP')

  end subroutine dragging

!________________________________________________________________________________________
!< calcule la valeur de beta dependant de la vitesse basale


  subroutine drag_non_lin_stag (nxx,nyy,Ux_mx,Uy_my,beta_x,beta_y)

    implicit none

                                                            !declarations variables dummy
    integer                               :: nxx,nyy        !< dimension des tableaux
    real, dimension(nxx,nyy),intent(in)   :: Ux_mx          !< vitesse de glissement x sur le noeud mineur x
    real, dimension(nxx,nyy),intent(in)   :: Uy_my          !< vitesse de glissement y sur le noeud mineur y
    real, dimension(nxx,nyy),intent(out)  :: beta_x         !< valeur de beta sur le noeud mineur x
    real, dimension(nxx,nyy),intent(out)  :: beta_y         !< valeur de beta sur le noeud mineur y

                                                          ! variables locales
    real, dimension(nxx,nyy)              :: Uxmy           ! vitesse selon x en my
    real, dimension(nxx,nyy)              :: Uymx           ! vitesse selon y en mx
    real, dimension(nxx,nyy)              :: Unorm_mx       !< norme vitesse de glissement sur le noeud mineur x
    real, dimension(nxx,nyy)              :: Unorm_my       !< norme vitesse de glissement sur le noeud mineur y


  if (itracebug.eq.1)  call tracebug(' entree dans drag_non_lin_stag MISMIP')


    do j=2,nyy-1
       do i=2,nxx
          Uymx(i,j)      = ((Uy_my(i,j)+Uy_my(i-1,j+1))+(Uy_my(i-1,j)+Uy_my(i,j+1)))/4.  ! vitesse selon y en mx
          Unorm_mx(i,j)  = (Ux_mx(i,j)**2+Uymx(i,j)**2)**0.5
        end do
    end do


    do j=2,nyy
       do i=2,nxx-1
          Uxmy(i,j)      = ((Ux_mx(i,j)+Ux_mx(i+1,j-1))+(Ux_mx(i,j-1)+Ux_mx(i+1,j)))/4.  ! vitesse selon y en mx
          Unorm_my(i,j)  = (Uy_my(i,j)**2+Uxmy(i,j)**2)**0.5
       end do
    end do


    do j=2,nyy-1
       do i=2,nxx

          if (Unorm_mx(i,j).gt.lim_unorm) then
             beta_x(i,j)    =  Coef_C * Unorm_mx(i,j)**m_1
          else
             beta_x(i,j)    =  beta_limgz
          end if
       end do
    end do


    do j=2,nyy
       do i=2,nxx-1

          if (Unorm_my(i,j).gt.lim_unorm) then
             beta_y(i,j)    =  Coef_C * Unorm_my(i,j)**m_1
          else
             beta_y(i,j)    =  beta_limgz
          end if
       end do
    end do

  end subroutine drag_non_lin_stag

  end module dragging_mismip3