!> \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(:,:) ! ux/uy(:,:,nz) should be used but only uxbar/uybar are updated by diagno_L2 ! anyway: ux/uy(:,:,nz) are uxbar/uybar (as it should be???) where (abs(uxbar(:,:)).gt.1) betamx(:,:) = tauc_mx(:,:) * ( abs(uxbar(:,:))**(q_nolin-1.) ) / ( u0 ** q_nolin ) elsewhere betamx(:,:) = tauc_mx(:,:) / ( u0 ** q_nolin ) endwhere where (abs(uybar(:,:)).gt.1) betamy(:,:) = tauc_my(:,:) * ( abs(uybar(:,:))**(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