1 | !> \file dragging_heino_mod.f90 |
---|
2 | !! Module qui defini la loi de glissement dans les experiences Heino |
---|
3 | !< |
---|
4 | |
---|
5 | !> \namespace sliding_heino |
---|
6 | !! Definition de la loi de glissement dans les experiences Heino |
---|
7 | !! \author Catherine |
---|
8 | !! \date ... |
---|
9 | !! @note Used module |
---|
10 | !! @note - use module3D_phy |
---|
11 | !< |
---|
12 | |
---|
13 | module sliding_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 | |
---|
35 | |
---|
36 | |
---|
37 | contains |
---|
38 | !------------------------------------------------------------------------------- |
---|
39 | !> SUBROUTINE: init_sliding |
---|
40 | !! Cette routine fait l'initialisation du glissement et est appelée juste apres |
---|
41 | !! la lecture du masque |
---|
42 | !< |
---|
43 | subroutine init_sliding |
---|
44 | |
---|
45 | |
---|
46 | if (itracebug.eq.1) call tracebug(' Heino: entree dans routine init_sliding') |
---|
47 | |
---|
48 | ! Cette routine fait l'initialisation du glissement et est appelée juste apres |
---|
49 | ! la lecture du masque |
---|
50 | |
---|
51 | |
---|
52 | ! ecriture dans le fichier parametres |
---|
53 | write(num_rep_42,*)'glissement Heino - initialisation masque' |
---|
54 | write(num_rep_42,*)'loi glissement=',loigliss |
---|
55 | write(num_rep_42,*)'----------------------------------------------------------' |
---|
56 | |
---|
57 | |
---|
58 | ! Le masque est égal au mk0 lu du début, |
---|
59 | |
---|
60 | mksedim(:,:)=mk0(:,:) |
---|
61 | where (mk0>1) mk0=1 |
---|
62 | |
---|
63 | !do j=ny, 1, -1 |
---|
64 | ! write(6,'(81i1)') (mk0(i,j), i=1,nx) |
---|
65 | !end do |
---|
66 | !write(6,*) mksedim(41,41) |
---|
67 | !do j=ny, 1, -1 |
---|
68 | ! write(6,'(81i1)') (mksedim(i,j), i=1,nx) |
---|
69 | !end do |
---|
70 | |
---|
71 | ! calcul du masque sur les demi mailles : frontiere -> sedim |
---|
72 | |
---|
73 | do j=2,ny-1 |
---|
74 | do i=2,nx-1 |
---|
75 | mkxsedim(i,j)=max(mksedim(i,j),mksedim(i-1,j)) |
---|
76 | mkysedim(i,j)=max(mksedim(i,j),mksedim(i,j-1)) |
---|
77 | end do |
---|
78 | end do |
---|
79 | |
---|
80 | |
---|
81 | Cr=1.e5 ! en a-1 coefficient hard rock (loi en puissance 3) |
---|
82 | Cs=500. ! en a-1 coefficient sediment (lineaire) |
---|
83 | |
---|
84 | write(num_rep_42,*)'coefficients: Cr=',Cr,' Cs=',Cs,' (a-1)' |
---|
85 | write(num_rep_42,*)'glissement seulement si les 2 voisins sont au point de fusion' |
---|
86 | |
---|
87 | if ((loigliss.eq.5).and.(Cs.gt.1.)) then |
---|
88 | coefdrag=rog/Cs |
---|
89 | Cs=0. |
---|
90 | write(num_rep_42,*)'passage en stream si zone sediment + fusion' |
---|
91 | write(num_rep_42,*)'coefficient de frottement dans la zone sediment=',coefdrag |
---|
92 | endif |
---|
93 | |
---|
94 | return |
---|
95 | end subroutine init_sliding |
---|
96 | |
---|
97 | !________________________________________________________________________________ |
---|
98 | |
---|
99 | !>SUBROUTINE: sliding |
---|
100 | !! Cette routine calcule le terme ddbx et ddby dans le cas Heino |
---|
101 | !< |
---|
102 | subroutine sliding |
---|
103 | |
---|
104 | ! cette routine calcule le terme ddbx et ddby dans le cas Heino |
---|
105 | ! ------------------------------------------------- |
---|
106 | ! GLISSEMENT |
---|
107 | ! DDBX et DDBY termes dus au glissement |
---|
108 | ! relation avec la vitesse de glissement UXB et UYB |
---|
109 | ! UXB=-DDBX*SDX et UYB=-DDBY*SDY |
---|
110 | ! ------------------------------------------------- |
---|
111 | |
---|
112 | |
---|
113 | if (itracebug.eq.1) call tracebug(' Heino: entree dans routine sliding') |
---|
114 | |
---|
115 | coefmax=1. |
---|
116 | |
---|
117 | slidy: do j=2,NY |
---|
118 | do I=2,NX-1 |
---|
119 | pose_y: if (.not.flgzmy(i,j)) then |
---|
120 | |
---|
121 | if (HMY(I,J).lt.1.) then ! glace tres fine |
---|
122 | ddby(i,j)=0. |
---|
123 | |
---|
124 | ! else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).and.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then |
---|
125 | ! ddby(i,j)=0. ! pas de glissement seulement |
---|
126 | ! ! si les 2 voisins sont froids |
---|
127 | |
---|
128 | else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then |
---|
129 | ddby(i,j)=0. ! pas de glissement si au moins un seul |
---|
130 | ! des 2 voisins est froid |
---|
131 | |
---|
132 | else ! glissement |
---|
133 | coefslid=(hwater(i,j)+hwater(i,j-1))*0.5 |
---|
134 | |
---|
135 | coefslid=max(coefslid,0.) |
---|
136 | coefslid=min(coefslid/coefmax,1.) |
---|
137 | coefslid=1. |
---|
138 | |
---|
139 | if (mkysedim(i,j).eq.1) then ! loi hard rock |
---|
140 | ddby(i,j)=Cr*slope2my(i,j)*Hmy(i,j)*coefslid |
---|
141 | |
---|
142 | else if (mkysedim(i,j).eq.2) then ! loi sediment |
---|
143 | ddby(i,j)=Cs*Hmy(i,j)*coefslid |
---|
144 | |
---|
145 | ! if (abs(sdy(i,j)).gt.1.e-8) then |
---|
146 | ! j_1=max(j-1,1) |
---|
147 | ! j1=min(j+1,ny) |
---|
148 | ! longslope=(sdy(i,j_1)+sdy(i,j)+sdy(i,j1))/3. |
---|
149 | ! ddby(i,j)=ddby(i,j)/sdy(i,j)*longslope |
---|
150 | ! endif |
---|
151 | |
---|
152 | else |
---|
153 | ddby(i,j)=0. |
---|
154 | endif |
---|
155 | |
---|
156 | endif |
---|
157 | |
---|
158 | |
---|
159 | endif pose_y |
---|
160 | end do |
---|
161 | end do slidy |
---|
162 | |
---|
163 | |
---|
164 | slidx: do j=2,NY-1 |
---|
165 | do I=2,NX |
---|
166 | pose_x: if (.not.flgzmx(i,j)) then |
---|
167 | |
---|
168 | if (HMx(I,J).lt.1.) then ! glace tres fine |
---|
169 | ddbx(i,j)=0. |
---|
170 | |
---|
171 | ! else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).and.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then |
---|
172 | ! ddbx(i,j)=0. ! pas de glissement seulement |
---|
173 | ! ! si les 2 voisins sont froids |
---|
174 | |
---|
175 | else if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then |
---|
176 | ddbx(i,j)=0. ! pas de glissement si au moins un seul |
---|
177 | ! des 2 voisins est froid |
---|
178 | |
---|
179 | |
---|
180 | |
---|
181 | else ! glissement |
---|
182 | |
---|
183 | coefslid=(hwater(i,j)+hwater(i-1,j))*0.5 |
---|
184 | |
---|
185 | coefslid=max(coefslid,0.) |
---|
186 | coefslid=min(coefslid/coefmax,1.) |
---|
187 | coefslid=1. |
---|
188 | |
---|
189 | if (mkxsedim(i,j).eq.1) then ! loi hard rock |
---|
190 | ddbx(i,j)=Cr*slope2mx(i,j)*Hmx(i,j)*coefslid |
---|
191 | |
---|
192 | else if (mkxsedim(i,j).eq.2) then ! loi sediment |
---|
193 | ddbx(i,j)=Cs*Hmx(i,j)*coefslid |
---|
194 | |
---|
195 | ! if (abs(sdx(i,j)).gt.1.e-8) then |
---|
196 | ! i_1=max(i-1,1) |
---|
197 | ! i1=min(i+1,nx) |
---|
198 | ! longslope=(sdx(i_1,j)+sdx(i,j)+sdx(i1,j))/3. |
---|
199 | ! ddbx(i,j)=ddbx(i,j)/sdx(i,j)*longslope |
---|
200 | ! endif |
---|
201 | |
---|
202 | |
---|
203 | else |
---|
204 | ddbx(i,j)=0. |
---|
205 | endif |
---|
206 | |
---|
207 | endif |
---|
208 | |
---|
209 | endif pose_x |
---|
210 | end do |
---|
211 | end do slidx |
---|
212 | |
---|
213 | return |
---|
214 | end subroutine sliding |
---|
215 | |
---|
216 | !------------------------------------------------------------------------- |
---|
217 | |
---|
218 | !> SUBROUTINE: dragging |
---|
219 | !! Defini le basal dragging |
---|
220 | !< |
---|
221 | |
---|
222 | subroutine dragging |
---|
223 | |
---|
224 | ! dans la zone sediment : appele si loigliss=5, |
---|
225 | ! defini le basal drag |
---|
226 | |
---|
227 | flgzmx(:,:)=.false. |
---|
228 | flgzmy(:,:)=.false. |
---|
229 | |
---|
230 | betamx(:,:)=1.e5 |
---|
231 | betamy(:,:)=1.e5 |
---|
232 | |
---|
233 | froty: do j=2,ny |
---|
234 | do i=2,nx-1 |
---|
235 | |
---|
236 | if (mkysedim(i,j).eq.2) then ! loi sediment |
---|
237 | |
---|
238 | ! seulement si au dessus du point de fusion |
---|
239 | if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i,j-1,nz).lt.Tpmp(i,j-1,nz))) then |
---|
240 | ! base froide |
---|
241 | else |
---|
242 | betamy(:,:)=coefdrag |
---|
243 | flgzmy(i,j)=.true. |
---|
244 | endif |
---|
245 | endif |
---|
246 | |
---|
247 | end do |
---|
248 | end do froty |
---|
249 | |
---|
250 | |
---|
251 | frotx: do j=2,ny-1 |
---|
252 | do i=2,nx |
---|
253 | |
---|
254 | if (mkxsedim(i,j).eq.2) then ! loi sediment |
---|
255 | |
---|
256 | ! seulement si au dessus du point de fusion |
---|
257 | if ((T(i,j,nz).lt.Tpmp(i,j,nz)).or.(T(i-1,j,nz).lt.Tpmp(i-1,j,nz))) then |
---|
258 | ! base froide |
---|
259 | else |
---|
260 | betamx(:,:)=coefdrag |
---|
261 | flgzmx(i,j)=.true. |
---|
262 | endif |
---|
263 | |
---|
264 | endif |
---|
265 | end do |
---|
266 | end do frotx |
---|
267 | |
---|
268 | gzmx(:,:)=flgzmx(:,:) |
---|
269 | gzmy(:,:)=flgzmy(:,:) |
---|
270 | gzmx_heino(:,:)=gzmx(:,:) |
---|
271 | gzmy_heino(:,:)=gzmy(:,:) |
---|
272 | |
---|
273 | return |
---|
274 | end subroutine dragging |
---|
275 | |
---|
276 | |
---|
277 | |
---|
278 | END MODULE sliding_heino |
---|
279 | |
---|