source: trunk/SOURCES/dragging_coulomb_friction_mod.f90 @ 348

Last change on this file since 348 was 348, checked in by aquiquet, 3 years ago

Bug correction for Coulomb friction law

File size: 10.2 KB
RevLine 
[333]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
[334]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)
[333]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
[334]57real, dimension(nx,ny) :: coef_sedim_mx, coef_sedim_my ! in case we use sediments
58
[333]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
[334]76use interface_input  !for sediments
77use io_netcdf_grisli !for sediments
78
[333]79implicit none
80
[334]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
[333]89
[334]90namelist/drag_coulomb_friction/cf,m_nolin,niter_nolin,betamax,betamin,bool_sedim,file_sedim,seuil_sedim,coef_sedim
91
[333]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
[334]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
[333]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!)'
[334]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'
[333]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
[334]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
[333]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
[334]236tauc_mx(:,:)= cf*neffmx(:,:)*coef_sedim_mx(:,:)
237tauc_my(:,:)= cf*neffmy(:,:)*coef_sedim_my(:,:)
[333]238
[348]239! ux/uy(:,:,nz) should be used but only uxbar/uybar are updated by diagno_L2
240! anyway: ux/uy(:,:,nz) are uxbar/uybar (as it should be???)
241where (abs(uxbar(:,:)).gt.1)
242   betamx(:,:) = tauc_mx(:,:) * ( abs(uxbar(:,:))**(q_nolin-1.) ) / ( u0 ** q_nolin )
[333]243elsewhere
244   betamx(:,:) = tauc_mx(:,:) / ( u0 ** q_nolin )
245endwhere
[348]246where (abs(uybar(:,:)).gt.1)
247   betamy(:,:) = tauc_my(:,:) * ( abs(uybar(:,:))**(q_nolin-1.) ) / ( u0 ** q_nolin )
[333]248elsewhere
249   betamy(:,:) = tauc_my(:,:) / ( u0 ** q_nolin )
250endwhere
251
252betamx(:,:)=max(betamx(:,:),betamin)
253betamy(:,:)=max(betamy(:,:),betamin)
254betamx(:,:)=min(betamx(:,:),betamax)
255betamy(:,:)=min(betamy(:,:),betamax)
256
257where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax
258where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax
259
260!$OMP END WORKSHARE
261
262
263! calcul de gzmx
264
265! points flottants
266!$OMP DO
267do j=2,ny
268   do i=2,nx
269      if (flot(i,j).and.(flot(i-1,j))) then
270         betamx(i,j)=0.
271      end if
272      if (flot(i,j).and.(flot(i,j-1))) then
273         betamy(i,j)=0.
274      end if
275   end do
276end do
277!$OMP END DO
278
279!$OMP DO
280do j=2,ny-1
281   do i=2,nx-1
282      beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) &
283          + (betamy(i,j)+betamy(i,j+1)))*0.25
284   end do
285end do
286!$OMP END DO
287!$OMP END PARALLEL
288
289return
290end subroutine dragging
291
292
293
294subroutine mstream_dragging   ! defini la localisation des streams
295
296 
297use module3d_phy, only: nx,ny,betamax,betamx,betamy,fleuvemx,fleuvemy,gzmx,gzmy,flgzmx,flgzmy,flot,flotmx,flotmy
298
299implicit none
300
301!$OMP PARALLEL
302
303!$OMP WORKSHARE
304fleuvemx(:,:)=.false.
305fleuvemy(:,:)=.false.
306gzmx(:,:)=.true.
307gzmy(:,:)=.true.
308flgzmx(:,:)=.false.
309flgzmy(:,:)=.false.
310!$OMP END WORKSHARE
311
312
313! calcul de gzmx
314
315! points flottants : flgzmx mais pas gzmx
316!$OMP DO
317do j=2,ny
318   do i=2,nx
319      if (flot(i,j).and.(flot(i-1,j))) then
320         flotmx(i,j)=.true.
321         flgzmx(i,j)=.true.
322         gzmx(i,j)=.false.
323         betamx(i,j)=0.
324      end if
325      if (flot(i,j).and.(flot(i,j-1))) then
326         flotmy(i,j)=.true.
327         flgzmy(i,j)=.true.
328         gzmy(i,j)=.false.
329         betamy(i,j)=0.
330      end if
331   end do
332end do
333!$OMP END DO
334
335!--------- autres criteres
336
337! points poses
338!$OMP WORKSHARE
339gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.betamax)   !  Pas de calcul pour les points
340gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.betamax)   !  au fort beta
341
342flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
343flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:)
344
345fleuvemx(:,:)=gzmx(:,:)
346fleuvemy(:,:)=gzmy(:,:)
347!$OMP END WORKSHARE
348
349!$OMP END PARALLEL
350
351return
352end subroutine mstream_dragging
353
354
355
356
357
358
359end module dragging_coulomb_friction
360
Note: See TracBrowser for help on using the repository browser.