1 | !> \file resol_adv_diff_2D.F90 |
---|
2 | !! Module pour la resolution de l'equation d advection diffusion 2D |
---|
3 | !< |
---|
4 | |
---|
5 | !> \namespace reso_adv_diff_2D |
---|
6 | !! Module pour la resolution de l'equation d advection diffusion 2D |
---|
7 | !! \author ... |
---|
8 | !! \date ... |
---|
9 | !! @note Used module |
---|
10 | !! @note - use module3D_phy |
---|
11 | !< |
---|
12 | |
---|
13 | module reso_adv_diff_2D |
---|
14 | |
---|
15 | use module3D_phy |
---|
16 | |
---|
17 | |
---|
18 | |
---|
19 | implicit none |
---|
20 | real :: omega ! parametre schema temporel de la resolution partie diffusion |
---|
21 | real :: mu_adv ! parametre schema temporel de la resolution advection |
---|
22 | |
---|
23 | real :: upwind ! schema spatial pour l'advection |
---|
24 | |
---|
25 | integer, dimension(2) :: ijmax ! position de maxval |
---|
26 | integer :: iFAIL |
---|
27 | |
---|
28 | contains |
---|
29 | |
---|
30 | subroutine init_reso_adv_diff_2D |
---|
31 | |
---|
32 | |
---|
33 | write(num_rep_42,*)'partie diffusive' |
---|
34 | write(num_rep_42,*)'le type de schema temporel diffusif depend de omega' |
---|
35 | write(num_rep_42,*)'0 -> explicite, 0.5 -> Crank-Nicolson' |
---|
36 | write(num_rep_42,*)'1 -> semi implicite, >1 -> over-implicite' |
---|
37 | |
---|
38 | omega = 2.5 |
---|
39 | |
---|
40 | write(num_rep_42,*)'omega = ',omega |
---|
41 | |
---|
42 | write(num_rep_42,*) |
---|
43 | write(num_rep_42,*)'partie advective' |
---|
44 | write(num_rep_42,*)' le schéma temporel advectif dépend de mu' |
---|
45 | write(num_rep_42,*)'0 -> explicite, 0.5 -> Crank-Nicolson' |
---|
46 | write(num_rep_42,*)'1 -> semi implicite, >1 -> over-implicite' |
---|
47 | |
---|
48 | mu_adv=1. |
---|
49 | write(num_rep_42,*)'mu_adv = ',mu_adv |
---|
50 | |
---|
51 | write(num_rep_42,*)' le schéma spatial dépend de upwind' |
---|
52 | write(num_rep_42,*)'1 -> schema upwind, 0.5 -> schema centre' |
---|
53 | |
---|
54 | upwind=1 |
---|
55 | write(num_rep_42,*)'upwind = ',upwind |
---|
56 | |
---|
57 | |
---|
58 | write(num_rep_42,*)'------------------------------------------------------' |
---|
59 | |
---|
60 | return |
---|
61 | end subroutine init_reso_adv_diff_2D |
---|
62 | !------------------------------------------------------------------ |
---|
63 | |
---|
64 | subroutine resol_adv_diff_2D (Dfx,Dfy,advx,advy,vieuxH) |
---|
65 | |
---|
66 | implicit none |
---|
67 | |
---|
68 | real,dimension(nx,ny), intent(in) :: Dfx ! terme diffusif selon x |
---|
69 | real,dimension(nx,ny), intent(in) :: Dfy ! terme diffusif selon y |
---|
70 | real,dimension(nx,ny), intent(in) :: Advx ! terme advectif selon x |
---|
71 | real,dimension(nx,ny), intent(in) :: Advy ! terme advectif selon x |
---|
72 | real,dimension(nx,ny), intent(in) :: vieuxH ! terme advectif selon x |
---|
73 | |
---|
74 | |
---|
75 | ! tableaux de travail. resolution M H = Frelax |
---|
76 | real,dimension(nx,ny) :: crelax ! diagnonale de M |
---|
77 | real,dimension(nx,ny) :: arelax ! sous diagonale selon x |
---|
78 | real,dimension(nx,ny) :: brelax ! sur diagonale selon x |
---|
79 | real,dimension(nx,ny) :: drelax ! sous diagonale selon y |
---|
80 | real,dimension(nx,ny) :: erelax ! sur diagonale selon y |
---|
81 | real,dimension(nx,ny) :: frelax ! vecteur |
---|
82 | real,dimension(nx,ny) :: c_west ! sur demi mailles Ux |
---|
83 | real,dimension(nx,ny) :: c_east ! sur demi mailles Ux |
---|
84 | real,dimension(nx,ny) :: c_north ! sur demi mailles Uy |
---|
85 | real,dimension(nx,ny) :: c_south !sur demi mailles Ux |
---|
86 | |
---|
87 | real,dimension(nx,ny) :: bdx ! pente socle |
---|
88 | real,dimension(nx,ny) :: bdy ! pente socle |
---|
89 | |
---|
90 | real,dimension(nx,ny) :: hdx ! pente epaisseur |
---|
91 | real,dimension(nx,ny) :: hdy ! pente epaisseur |
---|
92 | |
---|
93 | real :: frdx,frdy ! pour calcul frelax : termes diffusion |
---|
94 | real :: fraxw,fraxe,frays,frayn ! termes advection |
---|
95 | |
---|
96 | real,dimension(nx,ny) :: deltah ! dans calcul relax |
---|
97 | real :: delh ! dans calcul relax |
---|
98 | real :: testh ! dans calcul relax |
---|
99 | |
---|
100 | logical :: stopp |
---|
101 | integer :: ntour |
---|
102 | real :: reste |
---|
103 | |
---|
104 | |
---|
105 | integer :: it1,it2,jt1,jt2 ! pour des tests d'asymétrie |
---|
106 | |
---|
107 | ! attention H et bm sont passés par le module module3D_phy |
---|
108 | |
---|
109 | if (itracebug.eq.1) call tracebug(' Entree dans routine resolution_diffusion') |
---|
110 | |
---|
111 | ! calcul de bdx et hdx |
---|
112 | hdx(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0.,dim=1)) |
---|
113 | hdy(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0.,dim=2)) |
---|
114 | bdx(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0.,dim=1)) |
---|
115 | bdy(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0.,dim=2)) |
---|
116 | |
---|
117 | |
---|
118 | ! initialisations (qui feront aussi les conditions aux limites) |
---|
119 | Arelax(:,:)=0. |
---|
120 | Brelax(:,:)=0. |
---|
121 | Drelax(:,:)=0. |
---|
122 | Erelax(:,:)=0. |
---|
123 | Crelax(:,:)=1. |
---|
124 | Frelax(:,:)=0. |
---|
125 | DeltaH(:,:)=0. |
---|
126 | |
---|
127 | H(1,:)=0. ! bord gauche |
---|
128 | H(nx,:)=0. ! bord droit |
---|
129 | H(:,1)=0. ! bord bas |
---|
130 | H(:,ny)=0. ! bord haut |
---|
131 | |
---|
132 | |
---|
133 | ! schema spatial |
---|
134 | |
---|
135 | if (upwind.eq.1.) then !schema amont |
---|
136 | |
---|
137 | where (Advx(:,:).ge.0.) |
---|
138 | c_west(:,:)=1. |
---|
139 | c_east(:,:)=0. |
---|
140 | elsewhere |
---|
141 | c_west(:,:)=0. |
---|
142 | c_east(:,:)=1. |
---|
143 | end where |
---|
144 | |
---|
145 | where (Advy(:,:).ge.0.) |
---|
146 | c_south(:,:)=1. |
---|
147 | c_north(:,:)=0. |
---|
148 | elsewhere |
---|
149 | c_south(:,:)=0. |
---|
150 | c_north(:,:)=1. |
---|
151 | end where |
---|
152 | |
---|
153 | else if (upwind.lt.1.) then ! schema centre |
---|
154 | c_west(:,:)=0.5 |
---|
155 | c_east(:,:)=0.5 |
---|
156 | c_south(:,:)=0.5 |
---|
157 | c_north(:,:)=0.5 |
---|
158 | end if |
---|
159 | |
---|
160 | ! attribution des elements des diagonales |
---|
161 | do j=2,ny-1 |
---|
162 | do i=2,nx-1 |
---|
163 | |
---|
164 | ! sous diagonale en x |
---|
165 | arelax(i,j)=-omega*Dtdx2*Dfx(i,j) & ! partie diffusive en x |
---|
166 | -mu_adv*dtdx*c_west(i,j)*Advx(i,j) ! partie advective en x |
---|
167 | |
---|
168 | ! sur diagonale en x |
---|
169 | brelax(i,j)=-omega*Dtdx2*Dfx(i+1,j) & ! partie diffusive |
---|
170 | +mu_adv*dtdx*c_east(i+1,j)*Advx(i+1,j) ! partie advective |
---|
171 | |
---|
172 | ! sous diagonale en y |
---|
173 | drelax(i,j)=-omega*Dtdx2*Dfy(i,j) & ! partie diffusive en y |
---|
174 | -mu_adv*dtdx*c_south(i,j)*Advy(i,j) ! partie advective en y |
---|
175 | |
---|
176 | ! sur diagonale en y |
---|
177 | erelax(i,j)=-omega*Dtdx2*Dfy(i,j+1) & ! partie diffusive |
---|
178 | +mu_adv*dtdx*c_north(i,j+1)*Advy(i,j+1) ! partie advective |
---|
179 | |
---|
180 | |
---|
181 | |
---|
182 | ! diagonale |
---|
183 | crelax(i,j)=omega*Dtdx2*((Dfx(i+1,j)+Dfx(i,j)) & ! partie diffusive en x |
---|
184 | +(Dfy(i,j+1)+Dfy(i,j))) ! partie diffusive en y |
---|
185 | crelax(i,j)=crelax(i,j)+mu_adv*dtdx* & |
---|
186 | ( (c_west(i+1,j)*Advx(i+1,j)-c_east(i,j)*Advx(i,j)) & !partie advective en x |
---|
187 | +(c_south(i,j+1)*Advy(i,j+1)-c_north(i,j)*Advy(i,j))) !partie advective en y |
---|
188 | crelax(i,j)=1.+crelax(i,j) ! partie temporelle |
---|
189 | |
---|
190 | |
---|
191 | ! terme du vecteur |
---|
192 | |
---|
193 | frdx= -Dfx(i,j) * (Bdx(i,j) +(1.-omega)*Hdx(i,j)) & ! partie diffusive en x |
---|
194 | +Dfx(i+1,j)*(Bdx(i+1,j)+(1.-omega)*Hdx(i+1,j)) |
---|
195 | |
---|
196 | frdy= -Dfy(i,j) * (Bdy(i,j) +(1.-omega)*Hdy(i,j)) & ! partie diffusive en y |
---|
197 | +Dfy(i,j+1)*(Bdy(i,j+1)+(1.-omega)*Hdy(i,j+1)) |
---|
198 | |
---|
199 | fraxw= -c_west(i,j)* Advx(i,j) * vieuxH(i-1,j) & ! partie advective en x |
---|
200 | +c_west(i+1,j)*Advx(i+1,j)*vieuxH(i,j) ! venant de l'west |
---|
201 | |
---|
202 | fraxe= -c_east(i,j) * Advx(i,j) * vieuxH(i,j) & ! partie advective en x |
---|
203 | +c_east(i+1,j)*Advx(i+1,j)*vieuxH(i+1,j) ! venant de l'est |
---|
204 | |
---|
205 | frays= -c_south(i,j) * Advy(i,j) * vieuxH(i,j-1) & ! partie advective en y |
---|
206 | +c_south(i,j+1)*Advy(i,j+1)*vieuxH(i,j) ! venant du sud |
---|
207 | |
---|
208 | frayn= -c_north(i,j) * Advy(i,j) * vieuxH(i,j) & ! partie advective en y |
---|
209 | +c_north(i,j+1)*Advy(i,j+1)*vieuxH(i,j+1) ! venant du nord |
---|
210 | |
---|
211 | |
---|
212 | |
---|
213 | |
---|
214 | frelax(i,j)=vieuxH(i,j)+Dt*(bm(i,j)-bmelt(i,j))+dtdx*(frdx+frdy) & |
---|
215 | + (1.-mu_adv)*dtdx*((fraxw+fraxe)+(frays+frayn)) |
---|
216 | |
---|
217 | |
---|
218 | end do |
---|
219 | end do |
---|
220 | |
---|
221 | debug_3D(:,:,44)=bm(:,:)-bmelt(:,:) |
---|
222 | |
---|
223 | |
---|
224 | ! Boucle de relaxation : |
---|
225 | ! ---------------------- |
---|
226 | |
---|
227 | !!$write(166,*)'time,diagonale', time,dt,maxval(crelax(:,:)),maxloc(crelax) |
---|
228 | !!$write(166,*)'autres',maxval(abs(arelax(:,:))+abs(brelax(:,:))+abs(drelax(:,:))+abs(erelax(:,:))), & |
---|
229 | !!$maxloc(abs(arelax(:,:))+abs(brelax(:,:))+abs(drelax(:,:))+abs(erelax(:,:))) |
---|
230 | |
---|
231 | !!$ |
---|
232 | !!$if (time.gt.22100.) then |
---|
233 | !!$ijmax(:)=maxloc(abs(arelax(:,:))+abs(brelax(:,:))+abs(drelax(:,:))+abs(erelax(:,:))) |
---|
234 | !!$ |
---|
235 | !!$write(166,*)'arelax',arelax(ijmax(1),ijmax(2)),-omega*Dtdx2*dfx(ijmax(1),ijmax(2)),-mu_adv*dtdx*advx(ijmax(1),ijmax(2)) |
---|
236 | !!$write(166,*)'brelax',brelax(ijmax(1),ijmax(2)),-omega*Dtdx2*dfx(ijmax(1)+1,ijmax(2)),-mu_adv*dtdx*advx(ijmax(1)+1,ijmax(2)) |
---|
237 | !!$write(166,*)'drelax',drelax(ijmax(1),ijmax(2)),-omega*Dtdx2*dfy(ijmax(1),ijmax(2)),-mu_adv*dtdx*advy(ijmax(1),ijmax(2)) |
---|
238 | !!$write(166,*)'erelax',erelax(ijmax(1),ijmax(2)),-omega*Dtdx2*dfy(ijmax(1),ijmax(2)+1),-mu_adv*dtdx*advy(ijmax(1),ijmax(2)+1) |
---|
239 | !!$ |
---|
240 | !!$end if |
---|
241 | |
---|
242 | stopp = .false. |
---|
243 | ntour=0 |
---|
244 | |
---|
245 | !!$do j=2,ny-1 |
---|
246 | !!$ do i=2,nx-1 |
---|
247 | !!$ if (abs(crelax(i,j)).le.1.e-15) then |
---|
248 | !!$ write(166,*) 'time ',time,i,j,crelax(i,j) |
---|
249 | !!$ write(166,'(6(es12.5,1x))') advx(i-1,j),advx(i,j),advx(i+1,j),advy(i,j-1),advy(i,j),advy(i,j+1) |
---|
250 | !!$ write(166,'(6(es12.5,1x))') dfx(i-1,j),dfx(i,j),dfx(i+1,j),dfy(i,j-1),dfy(i,j),dfy(i,j+1) |
---|
251 | !!$ end if |
---|
252 | !!$ end do |
---|
253 | !!$end do |
---|
254 | |
---|
255 | !!$write(166,*) 'dtdx=',dtdx,'dt=',dt |
---|
256 | !!$write(166,*) 'arelax,brelax,crelax,drelax,erelax,frelax' |
---|
257 | !!$ |
---|
258 | !!$do j=43,45 |
---|
259 | !!$write(166,*) 'j=',j |
---|
260 | !!$ write(166,*) arelax(118:122,j) |
---|
261 | !!$ write(166,*) brelax(118:122,j) |
---|
262 | !!$ write(166,*) crelax(118:122,j) |
---|
263 | !!$ write(166,*) drelax(118:122,j) |
---|
264 | !!$ write(166,*) erelax(118:122,j) |
---|
265 | !!$ write(166,*) frelax(118:122,j) |
---|
266 | !!$ |
---|
267 | !!$end do |
---|
268 | |
---|
269 | |
---|
270 | relax_loop: do while(.not.stopp) |
---|
271 | ntour=ntour+1 |
---|
272 | !!$if (time.gt.22102.) then |
---|
273 | !!$ write(166,*) time, ntour,testh,dt |
---|
274 | !!$end if |
---|
275 | |
---|
276 | |
---|
277 | do j=2,ny-1 |
---|
278 | do i=2,nx-1 |
---|
279 | |
---|
280 | !!$ reste = (((arelax(i,j)*H(i-1,j) + brelax(i,j)*H(i+1,j)) & |
---|
281 | !!$ + (drelax(i,j)*H(i,j-1) + erelax(i,j)*H(i,j+1))) & |
---|
282 | !!$ + crelax(i,j)*H(i,j))- frelax(i,j) |
---|
283 | |
---|
284 | reste = (((arelax(i,j)*H(i-1,j) +drelax(i,j)*H(i,j-1)) & |
---|
285 | + (brelax(i,j)*H(i+1,j) + erelax(i,j)*H(i,j+1))) & |
---|
286 | + crelax(i,j)*H(i,j))- frelax(i,j) |
---|
287 | |
---|
288 | if (ntour.eq.1) debug_3D(i,j,49)=(((arelax(i,j)*H(i-1,j) +drelax(i,j)*H(i,j-1)) & |
---|
289 | + (brelax(i,j)*H(i+1,j) + erelax(i,j)*H(i,j+1))) & |
---|
290 | + crelax(i,j)*H(i,j)) |
---|
291 | |
---|
292 | |
---|
293 | deltaH(i,j) = reste/crelax(i,j) |
---|
294 | |
---|
295 | end do |
---|
296 | end do |
---|
297 | |
---|
298 | debug_3D(:,:,50)=arelax(:,:) |
---|
299 | debug_3D(:,:,51)=brelax(:,:) |
---|
300 | debug_3D(:,:,52)=crelax(:,:) |
---|
301 | debug_3D(:,:,53)=drelax(:,:) |
---|
302 | debug_3D(:,:,54)=erelax(:,:) |
---|
303 | debug_3D(:,:,55)=frelax(:,:) |
---|
304 | |
---|
305 | |
---|
306 | |
---|
307 | |
---|
308 | H(:,:)=H(:,:)-deltaH(:,:) |
---|
309 | |
---|
310 | !!$if (ntour.eq.1) then |
---|
311 | !!$write(166,*) 'deltaH' |
---|
312 | !!$ |
---|
313 | !!$do j=43,45 |
---|
314 | !!$write(166,*) 'j=',j |
---|
315 | !!$ write(166,*) deltaH(118:122,j) |
---|
316 | !!$end do |
---|
317 | !!$ |
---|
318 | !!$end if |
---|
319 | |
---|
320 | ! critere d'arret: |
---|
321 | ! ---------------- |
---|
322 | |
---|
323 | delh=0 |
---|
324 | |
---|
325 | |
---|
326 | do j=2,ny-1 |
---|
327 | do i=2,nx-1 |
---|
328 | delh=delh+deltaH(i,j)**2 |
---|
329 | end do |
---|
330 | end do |
---|
331 | |
---|
332 | !!$if (time.gt.22102.) then |
---|
333 | !!$ write(166,*) delh,maxval(deltaH(:,:)),maxloc(deltaH(:,:)) |
---|
334 | !!$end if |
---|
335 | |
---|
336 | if (delh.gt.0.) then |
---|
337 | testh=sqrt(delh)/((nx-2)*(ny-2)) |
---|
338 | else |
---|
339 | testh=0. |
---|
340 | endif |
---|
341 | stopp = (testh.lt.1.e-4).or.(ntour.gt.100) |
---|
342 | |
---|
343 | !write(166,*) 'sortie relax',time, ntour,testh,maxval(deltaH(:,:)) |
---|
344 | end do relax_loop |
---|
345 | |
---|
346 | !write(166,*) time,testh,ntour |
---|
347 | |
---|
348 | !!$! conditions aux limites |
---|
349 | !!$ |
---|
350 | |
---|
351 | ! tests d'asymétrie |
---|
352 | !!$it1=2 |
---|
353 | !!$it2=80 |
---|
354 | !!$jt1=41 |
---|
355 | !!$jt2=41 |
---|
356 | !!$if (H(it1,jt1).ne.H(it2,jt2)) then |
---|
357 | !!$ write(6,*)'asymetrie H time',time,'dans resol' |
---|
358 | !!$ write(6,*) 'bm',bm(it1,jt1),bm(it2,jt2) |
---|
359 | !!$ write(6,*) 'vieuxH',vieuxH(it1,jt1),vieuxH(it2,jt2) |
---|
360 | !!$ write(6,*) 'H',H(it1,jt1),H(it2,jt2) |
---|
361 | !!$ write(6,*) 'dfx',dfx(it1,jt1),dfx(it1+1,jt1),dfx(it2,jt2),dfx(it2+1,jt2) |
---|
362 | !!$ write(6,*) 'diffmx',diffmx(it1,jt1),diffmx(it1+1,jt1),diffmx(it2,jt2),diffmx(it2+1,jt2) |
---|
363 | !!$ write(6,*) 'uxbar',uxbar(it1,jt1),uxbar(it1+1,jt1),uxbar(it2,jt2),uxbar(it2+1,jt2) |
---|
364 | !!$ write(6,*) 'dfy',dfy(it1,jt1),dfy(it2,jt2),dfy(it1,jt1+1),dfy(it2,jt2+1) |
---|
365 | !!$ write(6,*) 'sdx',sdx(it1,jt1),sdx(it2,jt2),sdx(it1+1,jt1),sdx(it2-1,jt2) |
---|
366 | !!$ write(6,*) 'sdy',sdy(it1,jt1),sdy(it2,jt2),sdy(it1,jt1+1),sdy(it2,jt2+1) |
---|
367 | !!$ write(6,*) 'sdxmy',sdxmy(it1,jt1),sdxmy(it2,jt2),sdxmy(it1+1,jt1),sdxmy(it2+1,jt2) |
---|
368 | !!$ write(6,*) 'sdymx',sdymx(it1,jt1),sdymx(it2,jt2),sdymx(it1,jt1+1),sdymx(it2,jt2+1) |
---|
369 | !!$ write(6,*)'ddbx',ddbx(it1,jt1),ddby(it2,jt2) |
---|
370 | !!$ write(6,*)'ddx(i,j,1)',ddx(it1,it2,1),ddy(it2,it1,1) |
---|
371 | !!$ write(6,*) 'arelax',arelax(it1,jt1),arelax(it2,jt2),arelax(it1+1,jt1),arelax(it2-1,jt2) |
---|
372 | !!$ write(6,*) 'brelax',brelax(it1,jt1),brelax(it2,jt2),brelax(it1+1,jt1),brelax(it2-1,jt2) |
---|
373 | !!$ write(6,*) 'drelax',drelax(it1,jt1),drelax(it2,jt2) |
---|
374 | !!$ write(6,*) 'erelax',erelax(it1,jt1),erelax(it2,jt2) |
---|
375 | !!$ write(6,*) 'crelax',crelax(it1,jt1),crelax(it2,jt2) |
---|
376 | !!$ write(6,*) 'frelax',frelax(it1,jt1),frelax(it2,jt2) |
---|
377 | !!$ stop |
---|
378 | !!$endif |
---|
379 | |
---|
380 | |
---|
381 | |
---|
382 | ! write(166,*) '---------------------------------------------------------------' |
---|
383 | |
---|
384 | |
---|
385 | |
---|
386 | return |
---|
387 | |
---|
388 | end subroutine resol_adv_diff_2D |
---|
389 | !_______________________________________________________________________________ |
---|
390 | !> In this version, H can be prescribed on some nodes |
---|
391 | !! eventually changing with time |
---|
392 | !! value H_presc et test I_Hpresc |
---|
393 | subroutine resol_adv_diff_2D_Hpresc (Dfx,Dfy,advx,advy,vieuxH,H_presc,i_Hpresc) |
---|
394 | |
---|
395 | implicit none |
---|
396 | |
---|
397 | real,dimension(nx,ny), intent(in) :: Dfx ! terme diffusif selon x |
---|
398 | real,dimension(nx,ny), intent(in) :: Dfy ! terme diffusif selon y |
---|
399 | real,dimension(nx,ny), intent(in) :: Advx ! terme advectif selon x |
---|
400 | real,dimension(nx,ny), intent(in) :: Advy ! terme advectif selon x |
---|
401 | real,dimension(nx,ny), intent(in) :: vieuxH ! terme advectif selon x |
---|
402 | |
---|
403 | real,dimension(nx,ny), intent(in) :: H_presc ! H value if prescribed |
---|
404 | integer,dimension(nx,ny), intent(in) :: i_Hpresc ! 1 if H is prescribed on this node, else 0 |
---|
405 | |
---|
406 | |
---|
407 | ! tableaux de travail. resolution M H = Frelax |
---|
408 | real,dimension(nx,ny) :: crelax ! diagnonale de M |
---|
409 | real,dimension(nx,ny) :: arelax ! sous diagonale selon x |
---|
410 | real,dimension(nx,ny) :: brelax ! sur diagonale selon x |
---|
411 | real,dimension(nx,ny) :: drelax ! sous diagonale selon y |
---|
412 | real,dimension(nx,ny) :: erelax ! sur diagonale selon y |
---|
413 | real,dimension(nx,ny) :: frelax ! vecteur |
---|
414 | real,dimension(nx,ny) :: c_west ! sur demi mailles Ux |
---|
415 | real,dimension(nx,ny) :: c_east ! sur demi mailles Ux |
---|
416 | real,dimension(nx,ny) :: c_north ! sur demi mailles Uy |
---|
417 | real,dimension(nx,ny) :: c_south !sur demi mailles Ux |
---|
418 | |
---|
419 | real,dimension(nx,ny) :: bdx ! pente socle |
---|
420 | real,dimension(nx,ny) :: bdy ! pente socle |
---|
421 | |
---|
422 | real,dimension(nx,ny) :: hdx ! pente epaisseur |
---|
423 | real,dimension(nx,ny) :: hdy ! pente epaisseur |
---|
424 | |
---|
425 | real :: frdx,frdy ! pour calcul frelax : termes diffusion |
---|
426 | real :: fraxw,fraxe,frays,frayn ! termes advection |
---|
427 | |
---|
428 | real,dimension(nx,ny) :: deltah ! dans calcul relax |
---|
429 | real :: delh ! dans calcul relax |
---|
430 | real :: testh ! dans calcul relax |
---|
431 | |
---|
432 | logical :: stopp |
---|
433 | integer :: ntour |
---|
434 | real :: reste |
---|
435 | |
---|
436 | |
---|
437 | integer :: it1,it2,jt1,jt2 ! pour des tests d'asym??trie |
---|
438 | |
---|
439 | ! attention H et bm sont pass??s par le module module3D_phy |
---|
440 | |
---|
441 | if (itracebug.eq.1) call tracebug(' Entree dans routine resolution_diffusion') |
---|
442 | |
---|
443 | ! calcul de bdx et hdx |
---|
444 | hdx(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0.,dim=1)) |
---|
445 | hdy(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0.,dim=2)) |
---|
446 | bdx(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0.,dim=1)) |
---|
447 | bdy(:,:)=dx1*(B(:,:)-eoshift(B(:,:),shift=-1,boundary=0.,dim=2)) |
---|
448 | |
---|
449 | |
---|
450 | ! initialisations (qui feront aussi les conditions aux limites) |
---|
451 | Arelax(:,:)=0. |
---|
452 | Brelax(:,:)=0. |
---|
453 | Drelax(:,:)=0. |
---|
454 | Erelax(:,:)=0. |
---|
455 | Crelax(:,:)=1. |
---|
456 | Frelax(:,:)=0. |
---|
457 | DeltaH(:,:)=0. |
---|
458 | |
---|
459 | |
---|
460 | H(1,:)=0. ! bord gauche |
---|
461 | H(nx,:)=0. ! bord droit |
---|
462 | H(:,1)=0. ! bord bas |
---|
463 | H(:,ny)=0. ! bord haut |
---|
464 | |
---|
465 | |
---|
466 | ! schema spatial |
---|
467 | |
---|
468 | if (upwind.eq.1.) then !schema amont |
---|
469 | |
---|
470 | where (Advx(:,:).ge.0.) |
---|
471 | c_west(:,:)=1. |
---|
472 | c_east(:,:)=0. |
---|
473 | elsewhere |
---|
474 | c_west(:,:)=0. |
---|
475 | c_east(:,:)=1. |
---|
476 | end where |
---|
477 | |
---|
478 | where (Advy(:,:).ge.0.) |
---|
479 | c_south(:,:)=1. |
---|
480 | c_north(:,:)=0. |
---|
481 | elsewhere |
---|
482 | c_south(:,:)=0. |
---|
483 | c_north(:,:)=1. |
---|
484 | end where |
---|
485 | |
---|
486 | else if (upwind.lt.1.) then ! schema centre |
---|
487 | c_west(:,:)=0.5 |
---|
488 | c_east(:,:)=0.5 |
---|
489 | c_south(:,:)=0.5 |
---|
490 | c_north(:,:)=0.5 |
---|
491 | end if |
---|
492 | |
---|
493 | ! attribution des elements des diagonales |
---|
494 | do j=2,ny-1 |
---|
495 | do i=2,nx-1 |
---|
496 | |
---|
497 | ! sous diagonale en x |
---|
498 | arelax(i,j)=-omega*Dtdx2*Dfx(i,j) & ! partie diffusive en x |
---|
499 | -mu_adv*dtdx*c_west(i,j)*Advx(i,j) ! partie advective en x |
---|
500 | |
---|
501 | ! sur diagonale en x |
---|
502 | brelax(i,j)=-omega*Dtdx2*Dfx(i+1,j) & ! partie diffusive |
---|
503 | +mu_adv*dtdx*c_east(i+1,j)*Advx(i+1,j) ! partie advective |
---|
504 | |
---|
505 | ! sous diagonale en y |
---|
506 | drelax(i,j)=-omega*Dtdx2*Dfy(i,j) & ! partie diffusive en y |
---|
507 | -mu_adv*dtdx*c_south(i,j)*Advy(i,j) ! partie advective en y |
---|
508 | |
---|
509 | ! sur diagonale en y |
---|
510 | erelax(i,j)=-omega*Dtdx2*Dfy(i,j+1) & ! partie diffusive |
---|
511 | +mu_adv*dtdx*c_north(i,j+1)*Advy(i,j+1) ! partie advective |
---|
512 | |
---|
513 | |
---|
514 | |
---|
515 | ! diagonale |
---|
516 | crelax(i,j)=omega*Dtdx2*((Dfx(i+1,j)+Dfx(i,j)) & ! partie diffusive en x |
---|
517 | +(Dfy(i,j+1)+Dfy(i,j))) ! partie diffusive en y |
---|
518 | crelax(i,j)=crelax(i,j)+mu_adv*dtdx* & |
---|
519 | ( (c_west(i+1,j)*Advx(i+1,j)-c_east(i,j)*Advx(i,j)) & !partie advective en x |
---|
520 | +(c_south(i,j+1)*Advy(i,j+1)-c_north(i,j)*Advy(i,j))) !partie advective en y |
---|
521 | crelax(i,j)=1.+crelax(i,j) ! partie temporelle |
---|
522 | |
---|
523 | |
---|
524 | ! terme du vecteur |
---|
525 | |
---|
526 | frdx= -Dfx(i,j) * (Bdx(i,j) +(1.-omega)*Hdx(i,j)) & ! partie diffusive en x |
---|
527 | +Dfx(i+1,j)*(Bdx(i+1,j)+(1.-omega)*Hdx(i+1,j)) |
---|
528 | |
---|
529 | frdy= -Dfy(i,j) * (Bdy(i,j) +(1.-omega)*Hdy(i,j)) & ! partie diffusive en y |
---|
530 | +Dfy(i,j+1)*(Bdy(i,j+1)+(1.-omega)*Hdy(i,j+1)) |
---|
531 | |
---|
532 | fraxw= -c_west(i,j)* Advx(i,j) * vieuxH(i-1,j) & ! partie advective en x |
---|
533 | +c_west(i+1,j)*Advx(i+1,j)*vieuxH(i,j) ! venant de l'west |
---|
534 | |
---|
535 | fraxe= -c_east(i,j) * Advx(i,j) * vieuxH(i,j) & ! partie advective en x |
---|
536 | +c_east(i+1,j)*Advx(i+1,j)*vieuxH(i+1,j) ! venant de l'est |
---|
537 | |
---|
538 | frays= -c_south(i,j) * Advy(i,j) * vieuxH(i,j-1) & ! partie advective en y |
---|
539 | +c_south(i,j+1)*Advy(i,j+1)*vieuxH(i,j) ! venant du sud |
---|
540 | |
---|
541 | frayn= -c_north(i,j) * Advy(i,j) * vieuxH(i,j) & ! partie advective en y |
---|
542 | +c_north(i,j+1)*Advy(i,j+1)*vieuxH(i,j+1) ! venant du nord |
---|
543 | |
---|
544 | |
---|
545 | |
---|
546 | |
---|
547 | frelax(i,j)=vieuxH(i,j)+Dt*(bm(i,j)-bmelt(i,j))+dtdx*(frdx+frdy) & |
---|
548 | + (1.-mu_adv)*dtdx*((fraxw+fraxe)+(frays+frayn)) |
---|
549 | |
---|
550 | |
---|
551 | end do |
---|
552 | end do |
---|
553 | |
---|
554 | debug_3D(:,:,44)=bm(:,:)-bmelt(:,:) |
---|
555 | |
---|
556 | where (i_hpresc(:,:).eq.1) ! thickness prescribed |
---|
557 | frelax(:,:) = H_presc(:,:) |
---|
558 | arelax(:,:) = 0. |
---|
559 | brelax(:,:) = 0. |
---|
560 | crelax(:,:) = 1. |
---|
561 | drelax(:,:) = 0. |
---|
562 | erelax(:,:) = 0. |
---|
563 | end where |
---|
564 | |
---|
565 | |
---|
566 | |
---|
567 | |
---|
568 | stopp = .false. |
---|
569 | ntour=0 |
---|
570 | |
---|
571 | |
---|
572 | |
---|
573 | relax_loop: do while(.not.stopp) |
---|
574 | ntour=ntour+1 |
---|
575 | |
---|
576 | |
---|
577 | do j=2,ny-1 |
---|
578 | do i=2,nx-1 |
---|
579 | |
---|
580 | |
---|
581 | reste = (((arelax(i,j)*H(i-1,j) +drelax(i,j)*H(i,j-1)) & |
---|
582 | + (brelax(i,j)*H(i+1,j) + erelax(i,j)*H(i,j+1))) & |
---|
583 | + crelax(i,j)*H(i,j))- frelax(i,j) |
---|
584 | |
---|
585 | if (ntour.eq.1) debug_3D(i,j,49)=(((arelax(i,j)*H(i-1,j) +drelax(i,j)*H(i,j-1)) & |
---|
586 | + (brelax(i,j)*H(i+1,j) + erelax(i,j)*H(i,j+1))) & |
---|
587 | + crelax(i,j)*H(i,j)) |
---|
588 | |
---|
589 | |
---|
590 | deltaH(i,j) = reste/crelax(i,j) |
---|
591 | |
---|
592 | end do |
---|
593 | end do |
---|
594 | |
---|
595 | debug_3D(:,:,50)=arelax(:,:) |
---|
596 | debug_3D(:,:,51)=brelax(:,:) |
---|
597 | debug_3D(:,:,52)=crelax(:,:) |
---|
598 | debug_3D(:,:,53)=drelax(:,:) |
---|
599 | debug_3D(:,:,54)=erelax(:,:) |
---|
600 | debug_3D(:,:,55)=frelax(:,:) |
---|
601 | |
---|
602 | debug_3D(:,:,49)=-deltaH(:,:) |
---|
603 | |
---|
604 | |
---|
605 | H(:,:)=H(:,:)-deltaH(:,:) |
---|
606 | |
---|
607 | !!$if (ntour.eq.1) then |
---|
608 | !!$write(166,*) 'deltaH' |
---|
609 | !!$ |
---|
610 | !!$do j=43,45 |
---|
611 | !!$write(166,*) 'j=',j |
---|
612 | !!$ write(166,*) deltaH(118:122,j) |
---|
613 | !!$end do |
---|
614 | !!$ |
---|
615 | !!$end if |
---|
616 | |
---|
617 | ! critere d'arret: |
---|
618 | ! ---------------- |
---|
619 | |
---|
620 | delh=0 |
---|
621 | |
---|
622 | |
---|
623 | do j=2,ny-1 |
---|
624 | do i=2,nx-1 |
---|
625 | delh=delh+deltaH(i,j)**2 |
---|
626 | end do |
---|
627 | end do |
---|
628 | |
---|
629 | !!$if (time.gt.22102.) then |
---|
630 | !!$ write(166,*) delh,maxval(deltaH(:,:)),maxloc(deltaH(:,:)) |
---|
631 | !!$end if |
---|
632 | |
---|
633 | if (delh.gt.0.) then |
---|
634 | testh=sqrt(delh)/((nx-2)*(ny-2)) |
---|
635 | else |
---|
636 | testh=0. |
---|
637 | endif |
---|
638 | stopp = (testh.lt.1.e-4).or.(ntour.gt.100) |
---|
639 | |
---|
640 | !write(166,*) 'sortie relax',time, ntour,testh,maxval(deltaH(:,:)) |
---|
641 | end do relax_loop |
---|
642 | |
---|
643 | !write(166,*) time,testh,ntour |
---|
644 | |
---|
645 | !!$! conditions aux limites |
---|
646 | !!$ |
---|
647 | |
---|
648 | ! tests d'asym??trie |
---|
649 | !!$it1=2 |
---|
650 | !!$it2=80 |
---|
651 | !!$jt1=41 |
---|
652 | !!$jt2=41 |
---|
653 | !!$if (H(it1,jt1).ne.H(it2,jt2)) then |
---|
654 | !!$ write(6,*)'asymetrie H time',time,'dans resol' |
---|
655 | !!$ write(6,*) 'bm',bm(it1,jt1),bm(it2,jt2) |
---|
656 | !!$ write(6,*) 'vieuxH',vieuxH(it1,jt1),vieuxH(it2,jt2) |
---|
657 | !!$ write(6,*) 'H',H(it1,jt1),H(it2,jt2) |
---|
658 | !!$ write(6,*) 'dfx',dfx(it1,jt1),dfx(it1+1,jt1),dfx(it2,jt2),dfx(it2+1,jt2) |
---|
659 | !!$ write(6,*) 'diffmx',diffmx(it1,jt1),diffmx(it1+1,jt1),diffmx(it2,jt2),diffmx(it2+1,jt2) |
---|
660 | !!$ write(6,*) 'uxbar',uxbar(it1,jt1),uxbar(it1+1,jt1),uxbar(it2,jt2),uxbar(it2+1,jt2) |
---|
661 | !!$ write(6,*) 'dfy',dfy(it1,jt1),dfy(it2,jt2),dfy(it1,jt1+1),dfy(it2,jt2+1) |
---|
662 | !!$ write(6,*) 'sdx',sdx(it1,jt1),sdx(it2,jt2),sdx(it1+1,jt1),sdx(it2-1,jt2) |
---|
663 | !!$ write(6,*) 'sdy',sdy(it1,jt1),sdy(it2,jt2),sdy(it1,jt1+1),sdy(it2,jt2+1) |
---|
664 | !!$ write(6,*) 'sdxmy',sdxmy(it1,jt1),sdxmy(it2,jt2),sdxmy(it1+1,jt1),sdxmy(it2+1,jt2) |
---|
665 | !!$ write(6,*) 'sdymx',sdymx(it1,jt1),sdymx(it2,jt2),sdymx(it1,jt1+1),sdymx(it2,jt2+1) |
---|
666 | !!$ write(6,*)'ddbx',ddbx(it1,jt1),ddby(it2,jt2) |
---|
667 | !!$ write(6,*)'ddx(i,j,1)',ddx(it1,it2,1),ddy(it2,it1,1) |
---|
668 | !!$ write(6,*) 'arelax',arelax(it1,jt1),arelax(it2,jt2),arelax(it1+1,jt1),arelax(it2-1,jt2) |
---|
669 | !!$ write(6,*) 'brelax',brelax(it1,jt1),brelax(it2,jt2),brelax(it1+1,jt1),brelax(it2-1,jt2) |
---|
670 | !!$ write(6,*) 'drelax',drelax(it1,jt1),drelax(it2,jt2) |
---|
671 | !!$ write(6,*) 'erelax',erelax(it1,jt1),erelax(it2,jt2) |
---|
672 | !!$ write(6,*) 'crelax',crelax(it1,jt1),crelax(it2,jt2) |
---|
673 | !!$ write(6,*) 'frelax',frelax(it1,jt1),frelax(it2,jt2) |
---|
674 | !!$ stop |
---|
675 | !!$endif |
---|
676 | |
---|
677 | |
---|
678 | |
---|
679 | ! write(166,*) '---------------------------------------------------------------' |
---|
680 | |
---|
681 | |
---|
682 | |
---|
683 | return |
---|
684 | |
---|
685 | end subroutine resol_adv_diff_2D_Hpresc |
---|
686 | |
---|
687 | end module reso_adv_diff_2D |
---|
688 | |
---|
689 | |
---|
690 | |
---|
691 | |
---|