source: trunk/SOURCES/dragging_coulomb_friction_mod.f90 @ 334

Last change on this file since 334 was 334, checked in by aquiquet, 16 months ago

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

File size: 10.1 KB
Line 
1!> \file dragging_coulomb_friction_mod.f90
2!! Coulomb friction for a given value of the non-linearity
3!<
4
5!> \namespace  dragging_coulomb_friction
6!! Coulomb friction as defined in Pattyn,TC,2017 and reference therein
7!! \author Aurelien Quiquet
8!! \date 28/01/2021
9!! @note Non-linearity is provided in the param file.
10!! @note   * m should vary from 1 to infinity
11!! @note   * m=1 (q=1) -> linear friction as in Quiquet et al. (2018)
12!! @note   * m=infinity (q=0) -> -1 in param_list
13!! @note   * m=3 (q=1/3) is a classic (Weertman, 1974; Schoof, 2007)
14!! @note Used module
15!! @note   - use module3D_phy
16!! @note   - use param_phy_mod
17!! @note   - use fake_beta_iter_vitbil_mod
18!<
19
20module dragging_coulomb_friction
21
22!______________________________________________________________________________
23!
24! I took the formulation of Pattyn in his The Cryos. 2017 paper. It looks like:
25!
26!    taub = [ tauc / ( |ub|^(1-q) uO^q) ] ub
27
28! with:   -tauc is our cf x N (expressed differently in Pattyn).
29!         u0=100 m/yr
30!
31! What is inside [] is what we used as -Beta in Quiquet et al. (2018).
32!
33! For q=1, we get back our standard formulation (except that cf'=cf/100)
34!
35! In the param_list we will give:
36!    - cf as before
37!    - m=1/q (if m=-1 is provided -> q=0)
38!    - number of iterations (as it is non-linear)
39!    - betamin and betamax
40!
41! I also add the possibility to tune locally the cf coefficient to account for
42! the presence of sediments (could also be Bsoc_relax<0 for Antarctica)
43!______________________________________________________________________________
44
45
46use module3d_phy, only: nx,ny
47use fake_beta_iter_vitbil_mod
48
49implicit none
50
51real :: betamin   ! betamin : (Pa) frottement mini sous les streams
52
53real :: cf        !  A friction coefficient
54real :: m_nolin   !  Schoof non-linear exponent
55real :: q_nolin   !  q=1/m
56
57real, dimension(nx,ny) :: coef_sedim_mx, coef_sedim_my ! in case we use sediments
58
59real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat
60real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat
61real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat
62real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat
63logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta
64
65contains
66!-------------------------------------------------------------------------------
67
68!> SUBROUTINE: init_dragging
69!! Cette routine fait l'initialisation du dragging.
70!>
71subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
72
73use module3d_phy, only: niter_nolin,betamax,betamax_2d,inv_beta,mstream_mx,mstream_my,drag_mx,drag_my,num_rep_42
74use runparam,     only: itracebug
75
76use interface_input  !for sediments
77use io_netcdf_grisli !for sediments
78
79implicit none
80
81logical :: bool_sedim              ! sediments: yes or no
82real :: seuil_sedim                ! sediment thickness threshold for drag reduction
83real :: coef_sedim                 ! drag reduction factor
84real,dimension(nx,ny) :: h_sedim   ! sediment thickness (or rebounded bedrock)
85real,dimension(nx,ny) :: h_sedimmx,h_sedimmy   ! temporary arrays
86character(len=150) :: file_sedim   ! file for sediment thickness for HN or rebounded bsoc for Antar
87character(len=150) :: file_ncdf    ! fichier netcdf issue des fichiers .dat
88integer i,j
89
90namelist/drag_coulomb_friction/cf,m_nolin,niter_nolin,betamax,betamin,bool_sedim,file_sedim,seuil_sedim,coef_sedim
91
92if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope')
93
94
95! formats pour les ecritures dans 42
96428 format(A)
97
98! lecture des parametres du run                      block drag neff
99!--------------------------------------------------------------------
100
101rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
102read(num_param,drag_coulomb_friction)
103
104write(num_rep_42,428)'!___________________________________________________________' 
105write(num_rep_42,428) '&drag_coulomb_friction    ! nom du bloc dragging coulomb friction'
106write(num_rep_42,*)
107write(num_rep_42,*) 'cf                    = ', cf
108write(num_rep_42,*) 'm_nolin               = ', m_nolin
109write(num_rep_42,*) 'niter_nolin           = ', niter_nolin
110write(num_rep_42,*) 'betamax               = ', betamax
111write(num_rep_42,*) 'betamin               = ', betamin
112write(num_rep_42,*) 'bool_sedim            = ', bool_sedim
113write(num_rep_42,'(A,A,A)') 'file_sedim    = "',trim(file_sedim),'"'
114write(num_rep_42,*) 'seuil_sedim           = ', seuil_sedim
115write(num_rep_42,*) 'coef_sedim            = ', coef_sedim
116write(num_rep_42,*)'/'                           
117write(num_rep_42,428) '! cf: a friction coefficient (to be tuned)'
118write(num_rep_42,428) '! m_nolin: non-linear exponent, from 1 to infinity (put -1 for infinity)'
119write(num_rep_42,428) '!          m_nolin=1/q in Pattyn TC 2017, q in [0:1]'
120write(num_rep_42,428) '! niter_nolin: number of iterations to solve the non-linearity (expensive!)'
121write(num_rep_42,428) '! Possibility to add sediment tuning of cf, if false, the file is not read'
122write(num_rep_42,428) '! If sediments: where h_sedim > seuil_sedim, beta*coef_sedim'
123
124inv_beta=0
125!-------------------------------------------------------------------
126! masque stream
127mstream_mx(:,:)=1
128mstream_my(:,:)=1
129
130
131! coefficient permettant de modifier le basal drag.
132drag_mx(:,:)=1.
133drag_my(:,:)=1.
134     
135betamax_2d(:,:) = betamax
136
137if (m_nolin.lt.0) then
138   q_nolin = 0
139else
140   q_nolin = 1 / m_nolin
141endif
142
143if (bool_sedim) then
144   
145   file_sedim=trim(dirnameinp)//trim(file_sedim)
146   call lect_input(1,'h_sedim',1,h_sedim(:,:),file_sedim,file_ncdf)
147
148   h_sedimmx(1,:)=h_sedim(1,:)
149   h_sedimmy(:,1)=h_sedim(:,1)
150   do i=2,nx
151      h_sedimmx(i,:)=(h_sedim(i-1,:)+h_sedim(i,:))/2.
152   enddo
153   do j=2,ny
154      h_sedimmy(:,j)=(h_sedim(:,j-1)+h_sedim(:,j))/2.
155   enddo
156
157   where (h_sedimmx(:,:).gt.seuil_sedim)
158      coef_sedim_mx(:,:) = coef_sedim
159   elsewhere
160      coef_sedim_mx(:,:) = 1.
161   endwhere
162   where (h_sedimmy(:,:).gt.seuil_sedim)
163      coef_sedim_my(:,:) = coef_sedim
164   elsewhere
165      coef_sedim_my(:,:) = 1.
166   endwhere
167
168else
169   
170   coef_sedim_mx(:,:) = 1.
171   coef_sedim_my(:,:) = 1.
172
173endif
174
175return
176end subroutine init_dragging
177!________________________________________________________________________________
178
179!> SUBROUTINE: dragging
180!! Calcul le frottement basal
181!>
182!-------------------------------------------------------------------------
183subroutine dragging   ! defini la localisation des streams et le frottement basal
184!$ USE OMP_LIB
185
186use module3d_phy, only: nx,ny,betamax,beta_centre,betamx,betamy,neffmx,neffmy,hwater,flot,ux,uy
187
188implicit none
189
190real,dimension(nx,ny) :: tauc_mx   ! Yield stress, x direction
191real,dimension(nx,ny) :: tauc_my   ! Yield stress, y direction
192
193real,dimension(nx,ny) :: neff   ! pression effective noeuds majeurs
194real,dimension(nx,ny) :: hwatmx ! hauteur d'eau staggered grille - afq
195real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq
196
197real,parameter        :: u0 = 100d0 ! threshold sliding speed
198
199!$OMP PARALLEL
200
201! calcul de neff (pression effective sur noeuds majeurs)
202if (sum(neffmx(:,:)).le.0.) neffmx(:,:) =1.e8
203if (sum(neffmy(:,:)).le.0.) neffmy(:,:) =1.e8
204
205!$OMP DO
206do j=1,ny-1
207   do i=1,nx-1
208        neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4
209   enddo
210enddo
211!$OMP END DO
212!aurel, for the borders:
213!$OMP WORKSHARE
214neff(:,ny)=1.e5
215neff(nx,:)=1.e5
216! calcul de hwat (staggered)
217hwatmx(:,:)=0.
218hwatmy(:,:)=0.
219!$OMP END WORKSHARE
220!$OMP DO
221do i=2,nx
222   hwatmx(i,:)=(hwater(i-1,:)+hwater(i,:))/2.
223enddo
224!$OMP END DO
225!$OMP DO
226do j=2,ny
227   hwatmy(:,j)=(hwater(:,j-1)+hwater(:,j))/2.
228enddo
229!$OMP END DO
230
231!_____________________________
232! Coulomb friction computation:
233
234!$OMP WORKSHARE
235
236tauc_mx(:,:)= cf*neffmx(:,:)*coef_sedim_mx(:,:)
237tauc_my(:,:)= cf*neffmy(:,:)*coef_sedim_my(:,:)
238
239where (abs(ux(:,:,nz)).gt.1)
240   betamx(:,:) = tauc_mx(:,:) * ( abs(ux(:,:,nz))**(q_nolin-1.) ) / ( u0 ** q_nolin )
241elsewhere
242   betamx(:,:) = tauc_mx(:,:) / ( u0 ** q_nolin )
243endwhere
244where (abs(uy(:,:,nz)).gt.1)
245   betamy(:,:) = tauc_my(:,:) / ( abs(uy(:,:,nz))**(q_nolin-1.) ) / ( u0 ** q_nolin )
246elsewhere
247   betamy(:,:) = tauc_my(:,:) / ( u0 ** q_nolin )
248endwhere
249
250betamx(:,:)=max(betamx(:,:),betamin)
251betamy(:,:)=max(betamy(:,:),betamin)
252betamx(:,:)=min(betamx(:,:),betamax)
253betamy(:,:)=min(betamy(:,:),betamax)
254
255where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax
256where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax
257
258!$OMP END WORKSHARE
259
260
261! calcul de gzmx
262
263! points flottants
264!$OMP DO
265do j=2,ny
266   do i=2,nx
267      if (flot(i,j).and.(flot(i-1,j))) then
268         betamx(i,j)=0.
269      end if
270      if (flot(i,j).and.(flot(i,j-1))) then
271         betamy(i,j)=0.
272      end if
273   end do
274end do
275!$OMP END DO
276
277!$OMP DO
278do j=2,ny-1
279   do i=2,nx-1
280      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
281          + (betamy(i,j)+betamy(i,j+1)))*0.25
282   end do
283end do
284!$OMP END DO
285!$OMP END PARALLEL
286
287return
288end subroutine dragging
289
290
291
292subroutine mstream_dragging   ! defini la localisation des streams
293
294 
295use module3d_phy, only: nx,ny,betamax,betamx,betamy,fleuvemx,fleuvemy,gzmx,gzmy,flgzmx,flgzmy,flot,flotmx,flotmy
296
297implicit none
298
299!$OMP PARALLEL
300
301!$OMP WORKSHARE
302fleuvemx(:,:)=.false.
303fleuvemy(:,:)=.false.
304gzmx(:,:)=.true.
305gzmy(:,:)=.true.
306flgzmx(:,:)=.false.
307flgzmy(:,:)=.false.
308!$OMP END WORKSHARE
309
310
311! calcul de gzmx
312
313! points flottants : flgzmx mais pas gzmx
314!$OMP DO
315do j=2,ny
316   do i=2,nx
317      if (flot(i,j).and.(flot(i-1,j))) then
318         flotmx(i,j)=.true.
319         flgzmx(i,j)=.true.
320         gzmx(i,j)=.false.
321         betamx(i,j)=0.
322      end if
323      if (flot(i,j).and.(flot(i,j-1))) then
324         flotmy(i,j)=.true.
325         flgzmy(i,j)=.true.
326         gzmy(i,j)=.false.
327         betamy(i,j)=0.
328      end if
329   end do
330end do
331!$OMP END DO
332
333!--------- autres criteres
334
335! points poses
336!$OMP WORKSHARE
337gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.betamax)   !  Pas de calcul pour les points
338gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.betamax)   !  au fort beta
339
340flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
341flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
342
343fleuvemx(:,:)=gzmx(:,:)
344fleuvemy(:,:)=gzmy(:,:)
345!$OMP END WORKSHARE
346
347!$OMP END PARALLEL
348
349return
350end subroutine mstream_dragging
351
352
353
354
355
356
357end module dragging_coulomb_friction
358
Note: See TracBrowser for help on using the repository browser.