1 | !> \file sliding-Heino_mod.f90 |
---|
2 | !! Module qui definit la loi de glissement dans les experiences Heino. |
---|
3 | !< |
---|
4 | |
---|
5 | !> \namespace sliding_dragging_heino |
---|
6 | !! Definition de la loi de glissement dans les experiences Heino. |
---|
7 | !! \author ... |
---|
8 | !! \date ... |
---|
9 | !! @note Used module |
---|
10 | !! @note - use module3D_phy |
---|
11 | !< |
---|
12 | |
---|
13 | module sliding_dragging_heino |
---|
14 | |
---|
15 | ! Défini la loi de glissement dans les expériences Heino |
---|
16 | |
---|
17 | USE module3D_phy |
---|
18 | |
---|
19 | |
---|
20 | implicit none |
---|
21 | |
---|
22 | |
---|
23 | integer,dimension(nx,ny) :: mksedim !< masque des regions sediment-hard rock |
---|
24 | integer,dimension(nx,ny) :: mkxsedim !< décalle sur les demi mailles |
---|
25 | integer,dimension(nx,ny) :: mkysedim !< décalle sur les demi mailles |
---|
26 | real :: Cr !< coefficient glissement pour hard rock |
---|
27 | real :: Cs !< coefficient glissement pour sediment |
---|
28 | real :: coefmax |
---|
29 | real :: coefslid |
---|
30 | real :: coefdrag |
---|
31 | !logical,dimension(nx,ny) :: gzmx_heino !< pour transporter la valeur dans flottab |
---|
32 | !logical,dimension(nx,ny) :: gzmy_heino |
---|
33 | real :: longslope |
---|
34 | real :: seuil_gliss !< seuil sur hw_mx pour avoir du glissement |
---|
35 | real, dimension(nx,ny) :: hw_mx |
---|
36 | real, dimension(nx,ny) :: hw_my |
---|
37 | |
---|
38 | character(len=80) :: filin |
---|
39 | |
---|
40 | |
---|
41 | contains |
---|
42 | !------------------------------------------------------------------------------- |
---|
43 | !> SUBROUTINE: init_sliding |
---|
44 | !! Cette routine fait l'initialisation du glissement et est appelée juste apres |
---|
45 | !! la lecture du masque |
---|
46 | !< |
---|
47 | subroutine init_sliding |
---|
48 | |
---|
49 | if (itracebug.eq.1) call tracebug(' Heino: entree dans routine init_sliding') |
---|
50 | |
---|
51 | ! Cette routine fait l'initialisation du glissement et est appelée juste apres |
---|
52 | ! la lecture du masque |
---|
53 | |
---|
54 | ! Cette routine fait l'initialisation du glissement |
---|
55 | |
---|
56 | loigliss=4 ! 4 -> que glissement Heino (pas de mode dragging) |
---|
57 | if (loigliss.eq.4) then |
---|
58 | flgzmx(:,:)=.false. |
---|
59 | flgzmy(:,:)=.false. |
---|
60 | endif |
---|
61 | coefmax=hwatermax |
---|
62 | seuil_gliss=0.0001 |
---|
63 | ! 5 -> glissement Heino + dragging |
---|
64 | |
---|
65 | ! ecriture dans le fichier parametres |
---|
66 | write(num_rep_42,*)'glissement Heino - initialisation masque' |
---|
67 | write(num_rep_42,*)'loi glissement=',loigliss |
---|
68 | ! write(num_rep_42,*)'coef=1, gliss si les 2 voisins temperes' |
---|
69 | write(num_rep_42,*)'glissement si eau basale (6 voisins) >=seuil_gliss=',seuil_gliss |
---|
70 | write(num_rep_42,*)'coef dependant de la hauteur d eau coefmax=',coefmax |
---|
71 | write(num_rep_42,*)' moyenne de la hauteur d eau a partir de 6 voisins' |
---|
72 | write(num_rep_42,*)'----------------------------------------------------------' |
---|
73 | |
---|
74 | |
---|
75 | ! fichier mask fourni par l'intercomparaison Heino |
---|
76 | ! on refait la lecture |
---|
77 | |
---|
78 | filin='mask.dat' |
---|
79 | filin = TRIM(DIRNAMEINP)//TRIM(filin) |
---|
80 | |
---|
81 | write(num_rep_42,*) 'type de socle ' |
---|
82 | write(num_rep_42,*) 'mksedim donne par ',filin |
---|
83 | |
---|
84 | open(num_coupe,file=filin) |
---|
85 | do i=1,6 |
---|
86 | read(num_coupe,*) ! 6 lignes de commentaires |
---|
87 | end do |
---|
88 | |
---|
89 | ! lecture proprement dite de la carte sediment |
---|
90 | do j=ny, 1, -1 |
---|
91 | read(num_coupe,'(81i1)') (mksedim(i,j), i=1,nx) |
---|
92 | end do |
---|
93 | |
---|
94 | close(num_coupe) |
---|
95 | |
---|
96 | |
---|
97 | ! test pour une calotte axi-symetrique avec glissement Cr partout |
---|
98 | !write(num_rep_42,*) 'masque sediment mis a 1 ' |
---|
99 | !mksedim(:,:)=min(mksedim(:,:),1) |
---|
100 | |
---|
101 | ! a la frontiere le demi noeud est sediment. |
---|
102 | ! calcul du masque sur les demi mailles : frontiere -> sedim |
---|
103 | mkxsedim(:,:)=max(mksedim(:,:),eoshift(mksedim(:,:),shift=-1,boundary=0,dim=1)) |
---|
104 | mkysedim(:,:)=max(mksedim(:,:),eoshift(mksedim(:,:),shift=-1,boundary=0,dim=2)) |
---|
105 | write(num_rep_42,*) 'sur les bords zone sediment, masque exterieur ' |
---|
106 | |
---|
107 | |
---|
108 | !mkxsedim(:,:)=min(mksedim(:,:),eoshift(mksedim(:,:),shift=-1,boundary=0,dim=1)) |
---|
109 | !mkysedim(:,:)=min(mksedim(:,:),eoshift(mksedim(:,:),shift=-1,boundary=0,dim=2)) |
---|
110 | |
---|
111 | !write(num_rep_42,*) 'sur les bords zone sediment, demi-maille sediment si les 2 mailles le sont' |
---|
112 | |
---|
113 | mslid_mx(:,:)=mkxsedim(:,:) |
---|
114 | mslid_my(:,:)=mkysedim(:,:) |
---|
115 | |
---|
116 | Cr=1.e5 ! en a-1 coefficient hard rock (loi en puissance 3) |
---|
117 | Cs=500. ! en a-1 coefficient sediment (lineaire) |
---|
118 | |
---|
119 | |
---|
120 | write(num_rep_42,*)'coefficients: Cr=',Cr,' Cs=',Cs,' (a-1)' |
---|
121 | ! write(num_rep_42,*)'glissement seulement si les 2 voisins sont au point de fusion' |
---|
122 | |
---|
123 | |
---|
124 | |
---|
125 | return |
---|
126 | end subroutine init_sliding |
---|
127 | |
---|
128 | !________________________________________________________________________________ |
---|
129 | |
---|
130 | !> SUBROUTINE: init_dragging |
---|
131 | !! Initialisation du basal dragging |
---|
132 | !< |
---|
133 | subroutine init_dragging |
---|
134 | |
---|
135 | |
---|
136 | |
---|
137 | if ((loigliss.eq.5).and.(Cs.gt.1.)) then |
---|
138 | coefdrag=rog/Cs |
---|
139 | Cs=0. |
---|
140 | write(num_rep_42,*)'dragging Heino' |
---|
141 | write(num_rep_42,*)'passage en stream si zone sediment + fusion' |
---|
142 | write(num_rep_42,*)'coefficient de frottement dans la zone sediment (Pa/m)=',coefdrag |
---|
143 | endif |
---|
144 | |
---|
145 | end subroutine init_dragging |
---|
146 | |
---|
147 | !________________________________________________________________________________ |
---|
148 | |
---|
149 | !> SUBROUTINE: sliding |
---|
150 | !! Cette routine calcule le terme ddbx et ddby dans le cas Heino |
---|
151 | !< |
---|
152 | subroutine sliding |
---|
153 | |
---|
154 | |
---|
155 | ! cette routine calcule le terme ddbx et ddby dans le cas Heino |
---|
156 | ! ------------------------------------------------- |
---|
157 | ! GLISSEMENT |
---|
158 | ! DDBX et DDBY termes dus au glissement |
---|
159 | ! relation avec la vitesse de glissement UXB et UYB |
---|
160 | ! UXB=-DDBX*SDX et UYB=-DDBY*SDY |
---|
161 | ! ------------------------------------------------- |
---|
162 | |
---|
163 | implicit none |
---|
164 | |
---|
165 | integer :: i_moins1,j_moins1,i_plus1,j_plus1 |
---|
166 | |
---|
167 | !!$real, dimension(nx,ny) :: hw_mx |
---|
168 | !!$real, dimension(nx,ny) :: hw_my |
---|
169 | |
---|
170 | if (itracebug.eq.1) call tracebug(' Heino: entree dans routine sliding') |
---|
171 | |
---|
172 | |
---|
173 | call moy_mxmy(nx,ny,Hwater,hw_mx,hw_my) |
---|
174 | |
---|
175 | |
---|
176 | |
---|
177 | slidy: do j=2,NY |
---|
178 | do I=2,NX-1 |
---|
179 | |
---|
180 | pose_y: if (.not.flgzmy(i,j)) then |
---|
181 | |
---|
182 | if (HMY(I,J).lt.1.) then ! glace tres fine |
---|
183 | ddby(i,j)=0. |
---|
184 | mslid_my(i,j)=0 |
---|
185 | |
---|
186 | ! else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).and.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then |
---|
187 | ! ddby(i,j)=0. ! pas de glissement seulement |
---|
188 | ! ! si les 2 voisins sont froids |
---|
189 | |
---|
190 | ! else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then |
---|
191 | |
---|
192 | ! else if ((ibase(i,j).eq.1).or.(ibase(i,j-1).eq.1)) then |
---|
193 | ! else if ((ibase(i,j).eq.1).and.(ibase(i,j-1).eq.1)) then |
---|
194 | |
---|
195 | else if (hw_my(i,j).lt.seuil_gliss) then |
---|
196 | ddby(i,j)=0. ! pas de glissement si tous les voisins |
---|
197 | ! sont a sec |
---|
198 | mslid_my(i,j)=0 |
---|
199 | flgzmy(i,j)=.false. |
---|
200 | |
---|
201 | else ! glissement |
---|
202 | ! coefslid=(hwater(i,j)+hwater(i,j-1))*0.5 |
---|
203 | coefslid=hw_my(i,j) |
---|
204 | coefslid=max(coefslid,0.) |
---|
205 | coefslid=min(coefslid/coefmax,1.) |
---|
206 | ! coefslid=1. |
---|
207 | |
---|
208 | if (mkysedim(i,j).eq.1) then ! loi hard rock |
---|
209 | ddby(i,j)=Cr*slope2my(i,j)*Hmy(i,j)*coefslid |
---|
210 | |
---|
211 | mslid_my(i,j)=1 |
---|
212 | |
---|
213 | else if ((mkysedim(i,j).eq.2).and.(loigliss.eq.4)) then ! loi sediment |
---|
214 | ddby(i,j)=Cs*Hmy(i,j)*coefslid |
---|
215 | |
---|
216 | mslid_my(i,j)=2 |
---|
217 | |
---|
218 | ! if (abs(sdy(i,j)).gt.1.e-8) then |
---|
219 | ! j_1=max(j-1,1) |
---|
220 | ! j1=min(j+1,ny) |
---|
221 | ! longslope=(sdy(i,j_1)+sdy(i,j)+sdy(i,j1))/3. |
---|
222 | ! ddby(i,j)=ddby(i,j)/sdy(i,j)*longslope |
---|
223 | ! endif |
---|
224 | |
---|
225 | else |
---|
226 | ddby(i,j)=0. |
---|
227 | mslid_my(i,j)=0 |
---|
228 | |
---|
229 | endif |
---|
230 | |
---|
231 | endif |
---|
232 | |
---|
233 | |
---|
234 | endif pose_y |
---|
235 | end do |
---|
236 | end do slidy |
---|
237 | |
---|
238 | |
---|
239 | slidx: do j=2,NY-1 |
---|
240 | do I=2,NX |
---|
241 | pose_x: if (.not.flgzmx(i,j)) then |
---|
242 | |
---|
243 | if (HMx(I,J).lt.1.) then ! glace tres fine |
---|
244 | ddbx(i,j)=0. |
---|
245 | mslid_mx(i,j)=0 |
---|
246 | |
---|
247 | ! else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).and.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then |
---|
248 | ! ddbx(i,j)=0. ! pas de glissement seulement |
---|
249 | ! ! si les 2 voisins sont froids |
---|
250 | |
---|
251 | ! else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then |
---|
252 | ! else if ((ibase(i,j).eq.1).or.(ibase(i-1,j).eq.1)) then |
---|
253 | |
---|
254 | ! else if ((ibase(i,j).eq.1).and.(ibase(i-1,j).eq.1)) then |
---|
255 | else if (hw_mx(i,j).lt.seuil_gliss) then |
---|
256 | |
---|
257 | ddbx(i,j)=0. ! pas de glissement si au moins un seul |
---|
258 | ! des 2 voisins est froid |
---|
259 | mslid_mx(i,j)=0 |
---|
260 | flgzmx(i,j)=.false. |
---|
261 | |
---|
262 | |
---|
263 | else ! glissement |
---|
264 | |
---|
265 | ! coefslid=(hwater(i,j)+hwater(i-1,j))*0.5 |
---|
266 | coefslid=hw_mx(i,j) |
---|
267 | coefslid=max(coefslid,0.) |
---|
268 | coefslid=min(coefslid/coefmax,1.) |
---|
269 | ! coefslid=1. |
---|
270 | |
---|
271 | if (mkxsedim(i,j).eq.1) then ! loi hard rock |
---|
272 | ddbx(i,j)=Cr*slope2mx(i,j)*Hmx(i,j)*coefslid |
---|
273 | mslid_mx(i,j)=1 |
---|
274 | |
---|
275 | else if ((mkxsedim(i,j).eq.2).and.(loigliss.eq.4)) then ! loi sediment |
---|
276 | ddbx(i,j)=Cs*Hmx(i,j)*coefslid |
---|
277 | mslid_mx(i,j)=2 |
---|
278 | |
---|
279 | ! if (abs(sdx(i,j)).gt.1.e-8) then |
---|
280 | ! i_1=max(i-1,1) |
---|
281 | ! i1=min(i+1,nx) |
---|
282 | ! longslope=(sdx(i_1,j)+sdx(i,j)+sdx(i1,j))/3. |
---|
283 | ! ddbx(i,j)=ddbx(i,j)/sdx(i,j)*longslope |
---|
284 | ! endif |
---|
285 | |
---|
286 | |
---|
287 | else |
---|
288 | ddbx(i,j)=0. |
---|
289 | mslid_mx(i,j)=0 |
---|
290 | |
---|
291 | endif |
---|
292 | |
---|
293 | endif |
---|
294 | |
---|
295 | endif pose_x |
---|
296 | end do |
---|
297 | end do slidx |
---|
298 | |
---|
299 | return |
---|
300 | end subroutine sliding |
---|
301 | |
---|
302 | !------------------------------------------------------------------------- |
---|
303 | |
---|
304 | !> SUBROUTINE: dragging |
---|
305 | !! Definit le basal drag |
---|
306 | !< |
---|
307 | |
---|
308 | subroutine dragging |
---|
309 | |
---|
310 | ! dans la zone sediment : appele si loigliss=5, |
---|
311 | ! defini le basal drag |
---|
312 | |
---|
313 | |
---|
314 | flgzmx(:,:)=.false. |
---|
315 | flgzmy(:,:)=.false. |
---|
316 | |
---|
317 | betamx(:,:)=1.e5 |
---|
318 | betamy(:,:)=1.e5 |
---|
319 | |
---|
320 | shelfy=.false. |
---|
321 | |
---|
322 | if (loigliss.eq.5) then ! dragging dans les regions mkysedim:2 |
---|
323 | |
---|
324 | call moy_mxmy(nx,ny,Hwater,hw_mx,hw_my) |
---|
325 | |
---|
326 | froty: do j=2,ny |
---|
327 | do i=2,nx-1 |
---|
328 | |
---|
329 | if (mkysedim(i,j).eq.2) then ! loi sediment |
---|
330 | |
---|
331 | ! seulement si au dessus du point de fusion |
---|
332 | ! if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then |
---|
333 | if (hw_my(i,j).lt.seuil_gliss) then |
---|
334 | ! base froide |
---|
335 | else |
---|
336 | coefslid=hw_my(i,j) |
---|
337 | coefslid=max(coefslid,seuil_gliss) |
---|
338 | coefslid=min(coefslid/coefmax,1.) |
---|
339 | |
---|
340 | betamy(i,j)=min(coefdrag/coefslid,1.e5) |
---|
341 | |
---|
342 | ! betamy(i,j)=coefdrag |
---|
343 | |
---|
344 | |
---|
345 | flgzmy(i,j)=.true. |
---|
346 | ddby(i,j)=0. |
---|
347 | mslid_my(i,j)=5 |
---|
348 | shelfy=.true. |
---|
349 | endif |
---|
350 | endif |
---|
351 | |
---|
352 | end do |
---|
353 | end do froty |
---|
354 | |
---|
355 | |
---|
356 | frotx: do j=2,ny-1 |
---|
357 | do i=2,nx |
---|
358 | |
---|
359 | if (mkxsedim(i,j).eq.2) then ! loi sediment |
---|
360 | |
---|
361 | ! seulement si au dessus du point de fusion |
---|
362 | ! if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then |
---|
363 | if (hw_mx(i,j).lt.seuil_gliss) then |
---|
364 | ! base froide |
---|
365 | else |
---|
366 | coefslid=hw_mx(i,j) |
---|
367 | coefslid=max(coefslid,seuil_gliss) |
---|
368 | coefslid=min(coefslid/coefmax,1.) |
---|
369 | betamx(i,j)=min(coefdrag/coefslid,1.e5) |
---|
370 | |
---|
371 | |
---|
372 | ! betamx(i,j)=coefdrag |
---|
373 | |
---|
374 | flgzmx(i,j)=.true. |
---|
375 | ddbx(i,j)=0. |
---|
376 | mslid_mx(i,j)=5 |
---|
377 | shelfy=.true. |
---|
378 | endif |
---|
379 | |
---|
380 | endif |
---|
381 | end do |
---|
382 | end do frotx |
---|
383 | |
---|
384 | gzmx(:,:)=flgzmx(:,:) |
---|
385 | gzmy(:,:)=flgzmy(:,:) |
---|
386 | ! gzmx_heino(:,:)=gzmx(:,:) |
---|
387 | ! gzmy_heino(:,:)=gzmy(:,:) |
---|
388 | |
---|
389 | |
---|
390 | |
---|
391 | end if |
---|
392 | |
---|
393 | return |
---|
394 | end subroutine dragging |
---|
395 | |
---|
396 | |
---|
397 | END MODULE sliding_dragging_heino |
---|
398 | |
---|
399 | |
---|