Changeset 334 for trunk/SOURCES/dragging_coulomb_friction_mod.f90

Timestamp:

01/30/21 10:50:17 (3 years ago)

Author:

aquiquet

Message:

Coulomb friction law: possibility to tune the drag coefficient given a sediment thickness map (or relaxed bedrock elevation)

File:

• trunk/SOURCES/dragging_coulomb_friction_mod.f90 (modified) (7 diffs)

trunk/SOURCES/dragging_coulomb_friction_mod.f90

@@ -38 +38 @@
 !    - 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)
 !______________________________________________________________________________
 
@@ -51 +54 @@
 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
@@ -69 +74 @@
 use runparam,     only: itracebug
 
+use interface_input  !for sediments
+use io_netcdf_grisli !for sediments
+
 implicit none
 
-namelist/drag_coulomb_friction/cf,m_nolin,niter_nolin,betamax,betamin
+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')
@@ -93 +110 @@
 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)'
@@ -98 +119 @@
 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
@@ -116 +139 @@
 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
 
@@ -179 +234 @@
 !$OMP WORKSHARE
 
-tauc_mx(:,:)= cf*neffmx(:,:)
-tauc_my(:,:)= cf*neffmy(:,:)
+tauc_mx(:,:)= cf*neffmx(:,:)*coef_sedim_mx(:,:)
+tauc_my(:,:)= cf*neffmy(:,:)*coef_sedim_my(:,:)
 
 where (abs(ux(:,:,nz)).gt.1)