Changeset 357 for trunk

Timestamp:

06/10/22 14:12:38 (2 years ago)

Author:

aquiquet

Message:

Corrections for Coulomb friction with no hydrology, we now follow PISM formulation for the drag coefficient

File:

• trunk/SOURCES/dragging_coulomb_friction_simplhydro_mod.f90 (modified) (9 diffs)

trunk/SOURCES/dragging_coulomb_friction_simplhydro_mod.f90

@@ -55 +55 @@
 
 real :: till_c0        !  A friction coefficient
-real :: TFA_coeff !  Till friction angle parameter, could be spatially variable
+real :: TFA_min !  Till friction angle parameter, min value (below bed_crit)
+real :: TFA_max !  Till friction angle parameter, max value (above bed_crit)
+real :: bed_crit  !  Critical bed elevation for change in TFA
 real :: Ob_coeff  !  A friction coefficient, could be spatially variable
 real :: m_nolin   !  Schoof non-linear exponent
@@ -94 +96 @@
 integer i,j
 
-namelist/drag_coulomb_friction_simplhydro/till_c0, TFA_coeff, Ob_coeff,m_nolin,niter_nolin,betamax,betamin,bool_sedim,file_sedim,seuil_sedim,coef_sedim
+namelist/drag_coulomb_friction_simplhydro/till_c0, TFA_min, TFA_max, bed_crit, Ob_coeff,m_nolin,niter_nolin,betamax,betamin,bool_sedim,file_sedim,seuil_sedim,coef_sedim
 
 
@@ -113 +115 @@
 write(num_rep_42,*)
 write(num_rep_42,*) 'till_c0               = ', till_c0
-write(num_rep_42,*) 'TFA_coeff             = ', TFA_coeff
+write(num_rep_42,*) 'TFA_min               = ', TFA_min
+write(num_rep_42,*) 'TFA_max               = ', TFA_max
+write(num_rep_42,*) 'bed_crit              = ', bed_crit
 write(num_rep_42,*) 'Ob_coeff              = ', Ob_coeff
 write(num_rep_42,*) 'm_nolin               = ', m_nolin
@@ -125 +129 @@
 write(num_rep_42,*)'/'                            
 write(num_rep_42,428) '! till_c0: sediment cohesion (usually = 0)'
-write(num_rep_42,428) '! TFA_coeff: till friction angle'
+write(num_rep_42,428) '! TFA_min/max: till friction angle'
+write(num_rep_42,428) '! bed_crit: to change the value of TFA with depth (<0 for below pi sea level'
 write(num_rep_42,428) '! Ob_coeff: a friction coefficient'
 write(num_rep_42,428) '! m_nolin: non-linear exponent, from 1 to infinity (put -1 for infinity)'
@@ -138 +143 @@
 mstream_mx(:,:)=1
 mstream_my(:,:)=1
-
 
 ! coefficient permettant de modifier le basal drag.
@@ -200 +204 @@
 implicit none
 
-! Mohr-Coulomb friction law without hydrology as described in Pattyn et al. 2017.  
+! Mohr-Coulomb friction law without hydrology as described in PISM Martin et al. TC (2011).  
+real, dimension(nx,ny) :: TFA    ! till friction angle 
+real, dimension(nx,ny) :: frac_porewater  ! water pressure as a fraction of the total overburden pressure
 real, dimension(nx,ny) :: effective_pressure    ! local simple Neff
 real, dimension(nx,ny) :: effective_pressure_mx ! local simple Neff
@@ -211 +217 @@
 real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq
 
+real,dimension(nx,ny) :: Bsea ! temporary variable
+
 real,parameter        :: u0 = 100d0 ! threshold sliding speed
 
@@ -216 +224 @@
 
 !$OMP WORKSHARE
+Bsea(:,:) = Bsoc(:,:)-sealevel_2d(:,:)
 hwatmx(:,:)=0.
 hwatmy(:,:)=0.
-effective_pressure(:,:)= max (  RO*G*H(:,:)-max(-ROFRESH*G*(B(:,:)-sealevel_2d(:,:)), 0.) , 1. ) !limited to 1Pa
+where (Bsea(:,:).le.bed_crit)
+   TFA(:,:) = TFA_min
+elsewhere (Bsea(:,:).le.0.)
+   TFA(:,:) = TFA_min + TFA_max * (1-Bsea(:,:)/bed_crit)
+elsewhere
+   TFA(:,:) = TFA_max
+endwhere
+
+! [0,1] linear elevation weighing:
+frac_porewater(:,:) = 1. - min ( max ( Bsea(:,:), 0. ) , 1000. ) / 1000.
+
+
+effective_pressure(:,:) = max ( RO*G*H(:,:) * (1.-0.96 * frac_porewater(:,:)) , 1. ) !limited to 1Pa
 !$OMP END WORKSHARE
+
 !$OMP DO
 do i=2,nx
@@ -238 +260 @@
 !$OMP WORKSHARE
 
-tauc_mx(:,:) = till_c0 + Ob_coeff*TAND(TFA_coeff) * effective_pressure_mx(:,:) * coef_sedim_mx(:,:)
-tauc_my(:,:) = till_c0 + Ob_coeff*TAND(TFA_coeff) * effective_pressure_my(:,:) * coef_sedim_my(:,:)
+tauc_mx(:,:) = till_c0 + Ob_coeff*TAND(TFA) * effective_pressure_mx(:,:) * coef_sedim_mx(:,:)
+tauc_my(:,:) = till_c0 + Ob_coeff*TAND(TFA) * effective_pressure_my(:,:) * coef_sedim_my(:,:)
 
 ! ux/uy(:,:,nz) should be used but only uxbar/uybar are updated by diagno_L2