source: trunk/SOURCES/dragging_coulomb_friction_mod.f90 @ 334

!> \file dragging_coulomb_friction_mod.f90
!! Coulomb friction for a given value of the non-linearity
!<

!> \namespace  dragging_coulomb_friction
!! Coulomb friction as defined in Pattyn,TC,2017 and reference therein
!! \author Aurelien Quiquet
!! \date 28/01/2021
!! @note Non-linearity is provided in the param file.
!! @note   * m should vary from 1 to infinity
!! @note   * m=1 (q=1) -> linear friction as in Quiquet et al. (2018)
!! @note   * m=infinity (q=0) -> -1 in param_list
!! @note   * m=3 (q=1/3) is a classic (Weertman, 1974; Schoof, 2007)
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!! @note   - use fake_beta_iter_vitbil_mod
!<

module dragging_coulomb_friction

!______________________________________________________________________________
!
! I took the formulation of Pattyn in his The Cryos. 2017 paper. It looks like:
! 
!    taub = [ tauc / ( |ub|^(1-q) uO^q) ] ub
!  
! with:   -tauc is our cf x N (expressed differently in Pattyn).
!         u0=100 m/yr 
!
! What is inside [] is what we used as -Beta in Quiquet et al. (2018).
!
! For q=1, we get back our standard formulation (except that cf'=cf/100)
!
! In the param_list we will give:
!    - cf as before
!    - m=1/q (if m=-1 is provided -> q=0)
!    - number of iterations (as it is non-linear)
!    - betamin and betamax
!
! I also add the possibility to tune locally the cf coefficient to account for
! the presence of sediments (could also be Bsoc_relax<0 for Antarctica)
!______________________________________________________________________________


use module3d_phy, only: nx,ny
use fake_beta_iter_vitbil_mod

implicit none

real :: betamin   ! betamin : (Pa) frottement mini sous les streams

real :: cf        !  A friction coefficient
real :: m_nolin   !  Schoof non-linear exponent
real :: q_nolin   !  q=1/m

real, dimension(nx,ny) :: coef_sedim_mx, coef_sedim_my ! in case we use sediments

real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta

contains
!-------------------------------------------------------------------------------

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

use module3d_phy, only: niter_nolin,betamax,betamax_2d,inv_beta,mstream_mx,mstream_my,drag_mx,drag_my,num_rep_42
use runparam,     only: itracebug

use interface_input  !for sediments
use io_netcdf_grisli !for sediments

implicit none

logical :: bool_sedim              ! sediments: yes or no
real :: seuil_sedim                ! sediment thickness threshold for drag reduction
real :: coef_sedim                 ! drag reduction factor
real,dimension(nx,ny) :: h_sedim   ! sediment thickness (or rebounded bedrock)
real,dimension(nx,ny) :: h_sedimmx,h_sedimmy   ! temporary arrays
character(len=150) :: file_sedim   ! file for sediment thickness for HN or rebounded bsoc for Antar
character(len=150) :: file_ncdf    ! fichier netcdf issue des fichiers .dat
integer i,j

namelist/drag_coulomb_friction/cf,m_nolin,niter_nolin,betamax,betamin,bool_sedim,file_sedim,seuil_sedim,coef_sedim

if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope')


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

! lecture des parametres du run                      block drag neff
!--------------------------------------------------------------------

rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
read(num_param,drag_coulomb_friction)

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&drag_coulomb_friction    ! nom du bloc dragging coulomb friction'
write(num_rep_42,*)
write(num_rep_42,*) 'cf                    = ', cf
write(num_rep_42,*) 'm_nolin               = ', m_nolin
write(num_rep_42,*) 'niter_nolin           = ', niter_nolin
write(num_rep_42,*) 'betamax               = ', betamax
write(num_rep_42,*) 'betamin               = ', betamin
write(num_rep_42,*) 'bool_sedim            = ', bool_sedim
write(num_rep_42,'(A,A,A)') 'file_sedim    = "',trim(file_sedim),'"'
write(num_rep_42,*) 'seuil_sedim           = ', seuil_sedim
write(num_rep_42,*) 'coef_sedim            = ', coef_sedim
write(num_rep_42,*)'/'                            
write(num_rep_42,428) '! cf: a friction coefficient (to be tuned)'
write(num_rep_42,428) '! m_nolin: non-linear exponent, from 1 to infinity (put -1 for infinity)'
write(num_rep_42,428) '!          m_nolin=1/q in Pattyn TC 2017, q in [0:1]'
write(num_rep_42,428) '! niter_nolin: number of iterations to solve the non-linearity (expensive!)'
write(num_rep_42,428) '! Possibility to add sediment tuning of cf, if false, the file is not read'
write(num_rep_42,428) '! If sediments: where h_sedim > seuil_sedim, beta*coef_sedim'

inv_beta=0
!-------------------------------------------------------------------
! masque stream
mstream_mx(:,:)=1
mstream_my(:,:)=1


! coefficient permettant de modifier le basal drag.
drag_mx(:,:)=1.
drag_my(:,:)=1.
      
betamax_2d(:,:) = betamax

if (m_nolin.lt.0) then
   q_nolin = 0
else
   q_nolin = 1 / m_nolin
endif

if (bool_sedim) then
   
   file_sedim=trim(dirnameinp)//trim(file_sedim)
   call lect_input(1,'h_sedim',1,h_sedim(:,:),file_sedim,file_ncdf)

   h_sedimmx(1,:)=h_sedim(1,:)
   h_sedimmy(:,1)=h_sedim(:,1)
   do i=2,nx
      h_sedimmx(i,:)=(h_sedim(i-1,:)+h_sedim(i,:))/2.
   enddo
   do j=2,ny
      h_sedimmy(:,j)=(h_sedim(:,j-1)+h_sedim(:,j))/2.
   enddo

   where (h_sedimmx(:,:).gt.seuil_sedim)
      coef_sedim_mx(:,:) = coef_sedim
   elsewhere
      coef_sedim_mx(:,:) = 1.
   endwhere
   where (h_sedimmy(:,:).gt.seuil_sedim)
      coef_sedim_my(:,:) = coef_sedim
   elsewhere
      coef_sedim_my(:,:) = 1.
   endwhere

else
   
   coef_sedim_mx(:,:) = 1.
   coef_sedim_my(:,:) = 1.

endif

return
end subroutine init_dragging
!________________________________________________________________________________

!> SUBROUTINE: dragging
!! Calcul le frottement basal
!>
!-------------------------------------------------------------------------
subroutine dragging   ! defini la localisation des streams et le frottement basal
!$ USE OMP_LIB

use module3d_phy, only: nx,ny,betamax,beta_centre,betamx,betamy,neffmx,neffmy,hwater,flot,ux,uy

implicit none

real,dimension(nx,ny) :: tauc_mx   ! Yield stress, x direction
real,dimension(nx,ny) :: tauc_my   ! Yield stress, y direction

real,dimension(nx,ny) :: neff   ! pression effective noeuds majeurs
real,dimension(nx,ny) :: hwatmx ! hauteur d'eau staggered grille - afq
real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq

real,parameter        :: u0 = 100d0 ! threshold sliding speed

!$OMP PARALLEL

! calcul de neff (pression effective sur noeuds majeurs)
if (sum(neffmx(:,:)).le.0.) neffmx(:,:) =1.e8
if (sum(neffmy(:,:)).le.0.) neffmy(:,:) =1.e8

!$OMP DO
do j=1,ny-1
   do i=1,nx-1
        neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
   enddo
enddo
!$OMP END DO
!aurel, for the borders:
!$OMP WORKSHARE
neff(:,ny)=1.e5
neff(nx,:)=1.e5
! calcul de hwat (staggered)
hwatmx(:,:)=0.
hwatmy(:,:)=0.
!$OMP END WORKSHARE
!$OMP DO
do i=2,nx
   hwatmx(i,:)=(hwater(i-1,:)+hwater(i,:))/2.
enddo
!$OMP END DO
!$OMP DO
do j=2,ny
   hwatmy(:,j)=(hwater(:,j-1)+hwater(:,j))/2.
enddo
!$OMP END DO

!_____________________________
! Coulomb friction computation:

!$OMP WORKSHARE

tauc_mx(:,:)= cf*neffmx(:,:)*coef_sedim_mx(:,:)
tauc_my(:,:)= cf*neffmy(:,:)*coef_sedim_my(:,:)

where (abs(ux(:,:,nz)).gt.1)
   betamx(:,:) = tauc_mx(:,:) * ( abs(ux(:,:,nz))**(q_nolin-1.) ) / ( u0 ** q_nolin )
elsewhere
   betamx(:,:) = tauc_mx(:,:) / ( u0 ** q_nolin )
endwhere
where (abs(uy(:,:,nz)).gt.1)
   betamy(:,:) = tauc_my(:,:) / ( abs(uy(:,:,nz))**(q_nolin-1.) ) / ( u0 ** q_nolin )
elsewhere
   betamy(:,:) = tauc_my(:,:) / ( u0 ** q_nolin )
endwhere

betamx(:,:)=max(betamx(:,:),betamin)
betamy(:,:)=max(betamy(:,:),betamin)
betamx(:,:)=min(betamx(:,:),betamax)
betamy(:,:)=min(betamy(:,:),betamax)

where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax
where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax

!$OMP END WORKSHARE


! calcul de gzmx 

! points flottants
!$OMP DO
do j=2,ny
   do i=2,nx
      if (flot(i,j).and.(flot(i-1,j))) then
         betamx(i,j)=0.
      end if
      if (flot(i,j).and.(flot(i,j-1))) then
         betamy(i,j)=0.
      end if
   end do
end do
!$OMP END DO

!$OMP DO
do j=2,ny-1
   do i=2,nx-1
      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
          + (betamy(i,j)+betamy(i,j+1)))*0.25
   end do
end do
!$OMP END DO
!$OMP END PARALLEL

return
end subroutine dragging



subroutine mstream_dragging   ! defini la localisation des streams

  
use module3d_phy, only: nx,ny,betamax,betamx,betamy,fleuvemx,fleuvemy,gzmx,gzmy,flgzmx,flgzmy,flot,flotmx,flotmy

implicit none

!$OMP PARALLEL

!$OMP WORKSHARE
fleuvemx(:,:)=.false.
fleuvemy(:,:)=.false.
gzmx(:,:)=.true.
gzmy(:,:)=.true.
flgzmx(:,:)=.false.
flgzmy(:,:)=.false.
!$OMP END WORKSHARE


! calcul de gzmx 

! points flottants : flgzmx mais pas gzmx
!$OMP DO
do j=2,ny
   do i=2,nx
      if (flot(i,j).and.(flot(i-1,j))) then
         flotmx(i,j)=.true.
         flgzmx(i,j)=.true.
         gzmx(i,j)=.false.
         betamx(i,j)=0.
      end if
      if (flot(i,j).and.(flot(i,j-1))) then
         flotmy(i,j)=.true.
         flgzmy(i,j)=.true.
         gzmy(i,j)=.false.
         betamy(i,j)=0.
      end if
   end do
end do
!$OMP END DO

!--------- autres criteres

! points poses
!$OMP WORKSHARE
gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.betamax)   !  Pas de calcul pour les points
gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.betamax)   !  au fort beta

flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)

fleuvemx(:,:)=gzmx(:,:)
fleuvemy(:,:)=gzmy(:,:)
!$OMP END WORKSHARE

!$OMP END PARALLEL

return
end subroutine mstream_dragging






end module dragging_coulomb_friction