1 | !> \file dragging_calc_beta_mod.f90 |
---|
2 | !! Module qui calcule le beta a partir de vitesses de bilan |
---|
3 | !< |
---|
4 | |
---|
5 | !> \namespace dragging_calc_beta |
---|
6 | !! Calcule le beta a partir de vitesses de bilan |
---|
7 | !! @note Il faut partir d'un cptr. |
---|
8 | !! \author Cat |
---|
9 | !! \date august 2010 |
---|
10 | !! @note Used module |
---|
11 | !! @note - use module3D_phy |
---|
12 | !< |
---|
13 | |
---|
14 | module dragging_calc_beta |
---|
15 | |
---|
16 | ! Calcule le beta a partir de vitesses de bilan |
---|
17 | ! il faut partir d'un cptr (pour avoir un bon champ de temperature) |
---|
18 | |
---|
19 | |
---|
20 | |
---|
21 | use module3d_phy |
---|
22 | use interface_input |
---|
23 | |
---|
24 | implicit none |
---|
25 | |
---|
26 | logical, dimension(nx,ny) :: gz_centre !< stream on major node |
---|
27 | real, dimension(nx,ny) :: Vcol_x !< vertically averaged velocity along x (balance) |
---|
28 | real, dimension(nx,ny) :: Vcol_y !< vertically averaged velocity along y (balance) |
---|
29 | |
---|
30 | real, dimension(nx,ny) :: Vsl_x !<sliding velocity x direction (balance) |
---|
31 | real, dimension(nx,ny) :: Vsl_y !<sliding velocity y direction (balance) |
---|
32 | |
---|
33 | real, dimension(nx,ny) :: Vcol2 !< square of vertically averaged velocity norme on major nodes (balance) |
---|
34 | ! real, dimension(nx,ny) :: beta_centre ! beta on major node (average) declare dans 3D |
---|
35 | real :: beta_limgz !< when beta gt beta_limgz -> not gzmx |
---|
36 | real :: beta_bord !< for the sides of ice streams |
---|
37 | real :: ubil_limgz !< when ubil < ubil_limgz -> not gzmx |
---|
38 | real :: coefbeta !< coefficient to ajust beta |
---|
39 | real :: ecart_quad !< somme of quadratic difference |
---|
40 | real :: umag2 !< square of veloc. at major node |
---|
41 | integer :: stagger= 0 !< iteration based on centered (stagger = 0) |
---|
42 | !< stagg nodes (stagger = 1) |
---|
43 | !< or redistributed (stagger = 2) |
---|
44 | logical :: corr_def = .true. !< for deformation correction |
---|
45 | |
---|
46 | |
---|
47 | contains |
---|
48 | !----------------------------------------------------------------------------------- |
---|
49 | subroutine init_dragging ! Cette routine fait l'initialisation du dragging. |
---|
50 | |
---|
51 | implicit none |
---|
52 | |
---|
53 | if (itracebug.eq.1) call tracebug(' Calc_beta subroutine init_dragging') |
---|
54 | mstream_mx(:,:)=1 ! autorise streams partout |
---|
55 | mstream_my(:,:)=1 |
---|
56 | |
---|
57 | beta_limgz=5.e5 ! au dela duquel on considere pas de glissement |
---|
58 | beta_bord= 1000. ! pour regler des problemes a la cote |
---|
59 | ubil_limgz=.1 |
---|
60 | |
---|
61 | ! coefficient permettant de modifier le basal drag. |
---|
62 | drag_mx(:,:)=1. |
---|
63 | drag_my(:,:)=1. |
---|
64 | |
---|
65 | call lect_vitbil |
---|
66 | |
---|
67 | iter_beta = 1 |
---|
68 | iter_loop : if (iter_beta.eq.1) then ! pour calculer initial guess de beta a partir de vitesses |
---|
69 | |
---|
70 | betamx(:,:)=0 |
---|
71 | betamy(:,:)=0 |
---|
72 | beta_centre(:,:)=0. |
---|
73 | coefbeta=5. |
---|
74 | betamax = beta_limgz |
---|
75 | |
---|
76 | call gzm_betacalc ! determine flgzmx, ... |
---|
77 | |
---|
78 | else if (iter_beta.eq.0) then ! pour calculer les vitesses dynamiques deduites de iter_beta=1 |
---|
79 | call gzm_betacalc ! determine flgzmx, ... |
---|
80 | |
---|
81 | betamx(:,:) = max(0., betamx(:,:)) ! enleve les beta negatifs |
---|
82 | betamy(:,:) = max(0., betamy(:,:)) |
---|
83 | |
---|
84 | else if (iter_beta.eq.2) then ! pour partir d'un beta uniforme, (ne marche pas) |
---|
85 | |
---|
86 | betamx(:,:)=0 ! pour les tests gzm_betaclac |
---|
87 | betamy(:,:)=0 |
---|
88 | |
---|
89 | call gzm_betacalc ! determine flgzmx, ... |
---|
90 | |
---|
91 | betamx(:,:)=1.e5 |
---|
92 | betamy(:,:)=1.e5 |
---|
93 | beta_centre(:,:)=0. |
---|
94 | coefbeta=10. |
---|
95 | |
---|
96 | else if (iter_beta.eq.3) then ! iteration Arthern en partant d'un champ lu |
---|
97 | |
---|
98 | call lect_beta_stag |
---|
99 | call average_beta_centr |
---|
100 | call distribute_stag_beta ! averaging with minimum value |
---|
101 | call gzm_betacalc ! determine flgzmx, ... |
---|
102 | |
---|
103 | coefbeta=1.e-3 |
---|
104 | |
---|
105 | end if iter_loop |
---|
106 | |
---|
107 | return |
---|
108 | |
---|
109 | end subroutine init_dragging |
---|
110 | |
---|
111 | |
---|
112 | !------------------------------------------------------------------------- |
---|
113 | subroutine dragging ! defini la localisation des streams et le frottement basal |
---|
114 | |
---|
115 | implicit none |
---|
116 | |
---|
117 | |
---|
118 | if (itracebug.eq.1) call tracebug(' Subroutine dragging calc beta') |
---|
119 | if (itracebug.eq.1) write(num_tracebug,*)' Dragging_calc_beta: iter=',iter_beta, & |
---|
120 | ' corr_def=',corr_def |
---|
121 | call gzm_betacalc |
---|
122 | |
---|
123 | |
---|
124 | if (iter_beta.eq.1) then |
---|
125 | betamx(:,:)=0. |
---|
126 | betamy(:,:)=0. |
---|
127 | |
---|
128 | else if ((iter_beta.eq.0).and.(corr_def)) then |
---|
129 | call correct_from_def |
---|
130 | call gzm_betacalc |
---|
131 | |
---|
132 | else if (iter_beta.gt.1) then ! boucle iterative |
---|
133 | |
---|
134 | call Arthern_like_iter ! update betamx and betamy |
---|
135 | iter_beta=iter_beta+1 |
---|
136 | end if |
---|
137 | |
---|
138 | |
---|
139 | |
---|
140 | if (itracebug.eq.1) call tracebug(' Dragging sortie') |
---|
141 | |
---|
142 | return |
---|
143 | end subroutine dragging |
---|
144 | |
---|
145 | !________________________________________________________________________________ |
---|
146 | |
---|
147 | subroutine lect_vitbil |
---|
148 | |
---|
149 | |
---|
150 | character(len=100) :: balance_Ux_file ! balance velocity file Ux |
---|
151 | character(len=100) :: balance_Uy_file ! balance velocity file Uy |
---|
152 | |
---|
153 | namelist/vitbil_upwind/balance_Ux_file, balance_Uy_file |
---|
154 | |
---|
155 | if (itracebug.eq.1) call tracebug(' Subroutine lect_vitbil de calc_beta') |
---|
156 | |
---|
157 | ! lecture des parametres du run |
---|
158 | !-------------------------------------------------------------------- |
---|
159 | |
---|
160 | rewind(num_param) ! pour revenir au debut du fichier param_list.dat |
---|
161 | 428 format(A) |
---|
162 | read(num_param,vitbil_upwind) |
---|
163 | |
---|
164 | write(num_rep_42,428) '!___________________________________________________________' |
---|
165 | write(num_rep_42,428) '!read balance velocities on staggered grid' |
---|
166 | write(num_rep_42,vitbil_upwind) |
---|
167 | write(num_rep_42,428) '! balance_Ux_file : nom du fichier qui contient les vit. bilan Ux' |
---|
168 | write(num_rep_42,428) '! balance_Uy_file : nom du fichier qui contient les vit. bilan Uy' |
---|
169 | write(num_rep_42,*) |
---|
170 | |
---|
171 | ! read balance velocities on staggered nodes |
---|
172 | |
---|
173 | balance_Ux_file = trim(dirnameinp)//trim(balance_Ux_file) |
---|
174 | balance_Uy_file = trim(dirnameinp)//trim(balance_Uy_file) |
---|
175 | |
---|
176 | call lect_input(1,'Vcol_x',1,Vcol_x,balance_Ux_file,"") |
---|
177 | call lect_input(1,'Vcol_y',1,Vcol_y,balance_Uy_file,"") |
---|
178 | |
---|
179 | ! limit the maximum value |
---|
180 | |
---|
181 | Vcol_x(:,:)=max(-3900.,Vcol_x(:,:)) |
---|
182 | Vcol_x(:,:)=min( 3900.,Vcol_x(:,:)) |
---|
183 | Vcol_y(:,:)=max(-3900.,Vcol_y(:,:)) |
---|
184 | Vcol_y(:,:)=min( 3900.,Vcol_y(:,:)) |
---|
185 | |
---|
186 | |
---|
187 | do j=2,ny-1 |
---|
188 | do i=2,nx-1 |
---|
189 | umag2=((Vcol_x(i,j)+Vcol_x(i+1,j))*(Vcol_x(i,j)+Vcol_x(i+1,j))) & ! ux2 |
---|
190 | +((Vcol_y(i,j)+Vcol_y(i,j+1))*(Vcol_y(i,j)+Vcol_y(i,j+1))) ! uy2 |
---|
191 | Vcol2(i,j) =umag2*0.25 |
---|
192 | end do |
---|
193 | end do |
---|
194 | |
---|
195 | |
---|
196 | debug_3D(:,:,59)=Vcol_x(:,:) |
---|
197 | debug_3D(:,:,60)=Vcol_y(:,:) |
---|
198 | debug_3D(:,:,58)=Vcol2(:,:)**0.5 |
---|
199 | |
---|
200 | end subroutine lect_vitbil |
---|
201 | !__________________________________________________________________________________________ |
---|
202 | subroutine lect_beta_stag |
---|
203 | |
---|
204 | |
---|
205 | character(len=100) :: betamx_file ! betamx file |
---|
206 | character(len=100) :: betamy_file ! betamy file |
---|
207 | |
---|
208 | namelist/beta_stag/betamx_file,betamy_file,beta_limgz |
---|
209 | |
---|
210 | if (itracebug.eq.1) call tracebug(' Subroutine lect__beta_stag de calc_beta') |
---|
211 | |
---|
212 | ! lecture des parametres du run |
---|
213 | !-------------------------------------------------------------------- |
---|
214 | |
---|
215 | rewind(num_param) ! pour revenir au debut du fichier param_list.dat |
---|
216 | 428 format(A) |
---|
217 | read(num_param,beta_stag) |
---|
218 | |
---|
219 | write(num_rep_42,428) '!___________________________________________________________' |
---|
220 | write(num_rep_42,428) '!read beta on staggered grid' |
---|
221 | write(num_rep_42,beta_stag) |
---|
222 | write(num_rep_42,428) '! betamx_file : nom du fichier qui contient les betamx' |
---|
223 | write(num_rep_42,428) '! betamy_file : nom du fichier qui contient les betamy' |
---|
224 | write(num_rep_42,428) '! above beta_limgz, gzmx is false' |
---|
225 | write(num_rep_42,*) |
---|
226 | |
---|
227 | |
---|
228 | call lect_input(1,'betamx',1,betamx,betamx_file,"") |
---|
229 | call lect_input(1,'betamy',1,betamy,betamy_file,"") |
---|
230 | |
---|
231 | ! average on major nodes |
---|
232 | |
---|
233 | beta_centre(:,:)=0. |
---|
234 | do j=2,ny-1 |
---|
235 | do i=2,nx-1 |
---|
236 | beta_centre(i,j)= ((betamx(i,j)+betamx(i+1,j)) & |
---|
237 | + (betamy(i,j)+betamy(i,j+1)))*0.25 |
---|
238 | end do |
---|
239 | end do |
---|
240 | |
---|
241 | end subroutine lect_beta_stag |
---|
242 | !____________________________________________________________________________________ |
---|
243 | subroutine average_beta_centr ! average betamx and betamy on major nodes |
---|
244 | |
---|
245 | do j=2,ny-1 |
---|
246 | do i=2,nx-1 |
---|
247 | beta_centre(i,j) = (betamx(i,j) + betamx(i+1,j)) + (betamy(i,j) + betamy(i,j+1)) |
---|
248 | beta_centre(i,j) = beta_centre(i,j) * 0.25 |
---|
249 | end do |
---|
250 | end do |
---|
251 | end subroutine average_beta_centr |
---|
252 | |
---|
253 | !----------------------------------------------------------------------- |
---|
254 | subroutine distribute_stag_beta ! redistribute on staggered grid |
---|
255 | |
---|
256 | |
---|
257 | if (stagger.eq.0) then ! MIN(betamx,average beta_centr) |
---|
258 | do j=2,ny |
---|
259 | do i=2,nx |
---|
260 | betamx(i,j)= min(betamx(i,j),(beta_centre(i-1,j)+beta_centre(i,j))*0.5) |
---|
261 | betamy(i,j)= min(betamx(i,j),(beta_centre(i,j-1)+beta_centre(i,j))*0.5) |
---|
262 | end do |
---|
263 | end do |
---|
264 | else if (stagger.eq.1) then ! usual average |
---|
265 | do j=2,ny |
---|
266 | do i=2,nx |
---|
267 | betamx(i,j)=(beta_centre(i-1,j)+beta_centre(i,j))*0.5 |
---|
268 | betamy(i,j)=(beta_centre(i,j-1)+beta_centre(i,j))*0.5 |
---|
269 | end do |
---|
270 | end do |
---|
271 | |
---|
272 | else if (stagger.eq.2) then ! MIN(betamx,beta_bord) only on ice streams V> 50 m/year |
---|
273 | |
---|
274 | do j=2,ny |
---|
275 | do i=2,nx |
---|
276 | if (abs(Vcol_x(i,j)).gt.50.) then ! fleuve de glace |
---|
277 | betamx(i,j+1)=min(betamx(i,j+1),beta_bord) |
---|
278 | betamx(i,j-1)=min(betamx(i,j-1),beta_bord) |
---|
279 | betamx(i-1,j)=min(betamx(i-1,j),beta_bord) |
---|
280 | betamx(i+1,j)=min(betamx(i+1,j),beta_bord) |
---|
281 | end if |
---|
282 | if (abs(Vcol_y(i,j)).gt.50.) then ! fleuve de glace |
---|
283 | betamy(i+1,j)=min(betamy(i+1,j),beta_bord) |
---|
284 | betamy(i-1,j)=min(betamy(i-1,j),beta_bord) |
---|
285 | betamy(i,j+1)=min(betamy(i,j+1),beta_bord) |
---|
286 | betamy(i,j-1)=min(betamy(i,j-1),beta_bord) |
---|
287 | end if |
---|
288 | end do |
---|
289 | end do |
---|
290 | end if |
---|
291 | |
---|
292 | end subroutine distribute_stag_beta |
---|
293 | |
---|
294 | !---------------------------------------------------------------------------------- |
---|
295 | subroutine correct_from_def |
---|
296 | |
---|
297 | ! update la cible vitesse pour enlever la deformation |
---|
298 | |
---|
299 | if (itracebug.eq.1) call tracebug(' Subroutine correct_from_def') |
---|
300 | |
---|
301 | !Vsl_x(:,:) = Vcol_x(:,:) - uxdef(:,:) |
---|
302 | !Vsl_y(:,:) = Vcol_y(:,:) - uydef(:,:) |
---|
303 | |
---|
304 | do j = 2, ny |
---|
305 | do i = 2, nx |
---|
306 | if (Vcol_x(i,j).gt.0.) then |
---|
307 | Vsl_x(i,j) = max(0.,Vcol_x(i,j) - uxdef(i,j)) ! si la vitesse de deformation est > vitesse de bilan : pas de glissement |
---|
308 | |
---|
309 | else if (Vcol_x(i,j).lt.0.) then |
---|
310 | Vsl_x(i,j) = min(0.,Vcol_x(i,j) - uxdef(i,j)) |
---|
311 | else |
---|
312 | Vsl_x(i,j) = 0. |
---|
313 | endif |
---|
314 | |
---|
315 | if (Vcol_y(i,j).gt.0.) then |
---|
316 | Vsl_y(i,j) = max(0.,Vcol_y(i,j) - uydef(i,j)) ! si la vitesse de deformation est > vitesse de bilan : pas de glissement |
---|
317 | |
---|
318 | else if (Vcol_y(i,j).lt.0.) then |
---|
319 | Vsl_y(i,j) = min(0.,Vcol_y(i,j) - uydef(i,j)) |
---|
320 | else |
---|
321 | Vsl_y(i,j) = 0. |
---|
322 | endif |
---|
323 | |
---|
324 | end do |
---|
325 | end do |
---|
326 | |
---|
327 | do j = 2, ny-1 |
---|
328 | do i = 2, nx-1 |
---|
329 | debug_3D(:,:,68) = (((Vsl_x(i,j)+Vsl_x(i+1,j))**2 + & |
---|
330 | (Vsl_y(i,j)+Vsl_y(i,j+1))**2)*0.25)**0.5 |
---|
331 | debug_3D(:,:,69) = Vsl_x(i,j) |
---|
332 | debug_3D(:,:,70) = Vsl_y(i,j) |
---|
333 | end do |
---|
334 | end do |
---|
335 | |
---|
336 | end subroutine correct_from_def |
---|
337 | |
---|
338 | subroutine gzm_betacalc |
---|
339 | |
---|
340 | ! calcul de gzmx qui ne sera plus modifie ensuite |
---|
341 | if (itracebug.eq.1) call tracebug(' Subroutine gzm_betacalc') |
---|
342 | |
---|
343 | gzmx(:,:)=.true. |
---|
344 | gzmy(:,:)=.true. |
---|
345 | flgzmx(:,:)=.false. |
---|
346 | flgzmy(:,:)=.false. |
---|
347 | |
---|
348 | |
---|
349 | ! points flottants : flgzmx mais pas gzmx |
---|
350 | do j=2,ny |
---|
351 | do i=2,nx |
---|
352 | if (flot(i,j).and.(flot(i-1,j))) then |
---|
353 | flgzmx(i,j)=.true. |
---|
354 | gzmx(i,j)=.false. |
---|
355 | betamx(i,j)=0. |
---|
356 | flotmx(i,j) = .true. |
---|
357 | end if |
---|
358 | if (flot(i,j).and.(flot(i,j-1))) then |
---|
359 | flgzmy(i,j)=.true. |
---|
360 | gzmy(i,j)=.false. |
---|
361 | betamy(i,j)=0. |
---|
362 | flotmy(i,j) = .true. |
---|
363 | end if |
---|
364 | end do |
---|
365 | end do |
---|
366 | |
---|
367 | where (flot(:,:)) |
---|
368 | beta_centre(:,:)= 0. |
---|
369 | end where |
---|
370 | |
---|
371 | |
---|
372 | !--------- autres criteres |
---|
373 | |
---|
374 | ! points poses |
---|
375 | |
---|
376 | ! limite le calcul elliptique aux points a faible beta |
---|
377 | ! gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz) ! Pas de calcul pour les points |
---|
378 | ! gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz) ! au fort beta |
---|
379 | |
---|
380 | ! pour être sûr d'avoir quelques points "Dirichlet" |
---|
381 | gzmx(:,:)=gzmx(:,:).and.(abs(Vcol_x(:,:)).gt.ubil_limgz) ! Pas de calcul pour les pts lents |
---|
382 | gzmy(:,:)=gzmy(:,:).and.(abs(Vcol_y(:,:)).gt.ubil_limgz) ! Pas de calcul pour les pts lents |
---|
383 | |
---|
384 | |
---|
385 | flgzmx(:,:)=flgzmx(:,:).or.gzmx(:,:) |
---|
386 | flgzmy(:,:)=flgzmy(:,:).or.gzmy(:,:) |
---|
387 | |
---|
388 | gz_centre(:,:)=.false. |
---|
389 | do j=2,ny-1 |
---|
390 | do i=2,nx-1 |
---|
391 | gz_centre(i,j)=gzmx(i,j).or.gzmx(i+1,j).or.gzmy(i,j).or.gzmy(i,j+1) |
---|
392 | end do |
---|
393 | end do |
---|
394 | |
---|
395 | where (flgzmy(:,:)) |
---|
396 | debug_3D(:,:,84)=1 |
---|
397 | elsewhere |
---|
398 | debug_3D(:,:,84)=0 |
---|
399 | endwhere |
---|
400 | |
---|
401 | end subroutine gzm_betacalc |
---|
402 | !____________________________________________ |
---|
403 | |
---|
404 | |
---|
405 | subroutine Arthern_like_iter ! ATTENTION CE BLOC A ETE DEPLACE ET PAS REVALIDE |
---|
406 | |
---|
407 | ! iteration Arthern like sur beta |
---|
408 | !________________________________________ |
---|
409 | ! On part d un beta issus de calc beta et eventuellement lisse par moyenne mobile |
---|
410 | ! |
---|
411 | ! Vcol_x tient le role de solution Dirichlet et uxbar de solution Neumann |
---|
412 | |
---|
413 | |
---|
414 | ecart_quad=0. |
---|
415 | if (iter_beta.ge.4) then |
---|
416 | |
---|
417 | stag: if (stagger.eq.1) then |
---|
418 | do j=2,ny |
---|
419 | do i=2,nx |
---|
420 | if (gzmx(i,j)) then |
---|
421 | betamx(i,j)=betamx(i,j)-coefbeta*(Vcol_x(i,j)*Vcol_x(i,j)-uxbar(i,j)*uxbar(i,j)) |
---|
422 | betamx(i,j)=max(betamx(i,j),0.) |
---|
423 | ecart_quad= ecart_quad & |
---|
424 | + abs(Vcol_x(i,j)*Vcol_x(i,j)-uxbar(i,j)*uxbar(i,j)) |
---|
425 | end if |
---|
426 | if (gzmy(i,j)) then |
---|
427 | betamy(i,j)=betamy(i,j)-coefbeta*(Vcol_y(i,j)*Vcol_y(i,j)-uybar(i,j)*uybar(i,j)) |
---|
428 | betamy(i,j)=max(betamy(i,j),0.) |
---|
429 | ecart_quad= ecart_quad & |
---|
430 | + abs(Vcol_y(i,j)*Vcol_y(i,j)-uybar(i,j)*uybar(i,j)) |
---|
431 | end if |
---|
432 | end do |
---|
433 | end do |
---|
434 | else if (stagger.eq.0) then ! le chemin se fait sur les noeuds centres |
---|
435 | do j=2,ny-1 |
---|
436 | do i=2,nx-1 |
---|
437 | if (gz_centre(i,j)) then |
---|
438 | umag2=((uxbar(i,j)+uxbar(i+1,j))*(uxbar(i,j)+uxbar(i+1,j))) & ! ux2 |
---|
439 | +((uybar(i,j)+uybar(i,j+1))*(uybar(i,j)+uybar(i,j+1))) ! uy2 |
---|
440 | umag2=umag2*0.25 |
---|
441 | |
---|
442 | beta_centre(i,j)=beta_centre(i,j)-coefbeta*(Vcol2(i,j)-umag2) |
---|
443 | beta_centre(i,j)=max(beta_centre(i,j),0.) |
---|
444 | ecart_quad= ecart_quad & |
---|
445 | + abs(Vcol2(i,j)-umag2) |
---|
446 | end if |
---|
447 | end do |
---|
448 | end do |
---|
449 | ! redistribute on staggered grid |
---|
450 | do j=2,ny |
---|
451 | do i=2,nx |
---|
452 | betamx(i,j)= (beta_centre(i-1,j)+beta_centre(i,j))*0.5 |
---|
453 | betamy(i,j)= (beta_centre(i,j-1)+beta_centre(i,j))*0.5 |
---|
454 | end do |
---|
455 | end do |
---|
456 | |
---|
457 | else if (stagger.eq.2) then ! le chemin se fait sur les noeuds centres |
---|
458 | beta_centre(:,:)=0. |
---|
459 | write(6,*) coefbeta |
---|
460 | do j=2,ny-1 |
---|
461 | do i=2,nx-1 |
---|
462 | if (gz_centre(i,j)) then |
---|
463 | umag2=((uxbar(i,j)+uxbar(i+1,j))*(uxbar(i,j)+uxbar(i+1,j))) & ! ux2 |
---|
464 | +((uybar(i,j)+uybar(i,j+1))*(uybar(i,j)+uybar(i,j+1))) ! uy2 |
---|
465 | umag2=umag2*0.25 |
---|
466 | |
---|
467 | beta_centre(i,j)=-coefbeta*(Vcol2(i,j)-umag2) |
---|
468 | beta_centre(i,j)=min(beta_centre(i,j),100.) |
---|
469 | beta_centre(i,j)=max(beta_centre(i,j),-100.) |
---|
470 | ecart_quad= ecart_quad & |
---|
471 | + abs(Vcol2(i,j)-umag2) |
---|
472 | end if |
---|
473 | end do |
---|
474 | end do |
---|
475 | debug_3D(:,:,86)=beta_centre(:,:) |
---|
476 | |
---|
477 | ! redistribute on staggered grid |
---|
478 | do j=2,ny |
---|
479 | do i=2,nx |
---|
480 | betamx(i,j)= betamx(i,j)+((beta_centre(i-1,j)+beta_centre(i,j)) & |
---|
481 | + (beta_centre(i-1,j-1)+beta_centre(i,j-1)) & |
---|
482 | + (beta_centre(i-1,j+1)+beta_centre(i,j+1)))/6. |
---|
483 | |
---|
484 | betamx(i,j)=max(betamx(i,j),0.) |
---|
485 | betamx(i,j)=min(betamx(i,j),beta_limgz) |
---|
486 | |
---|
487 | betamy(i,j)= betamy(i,j)+((beta_centre(i,j-1)+beta_centre(i,j)) & |
---|
488 | + (beta_centre(i-1,j-1)+beta_centre(i-1,j)) & |
---|
489 | + (beta_centre(i+1,j-1)+beta_centre(i+1,j)))/6. |
---|
490 | |
---|
491 | betamy(i,j)=max(betamy(i,j),0.) |
---|
492 | betamy(i,j)=min(betamy(i,j),beta_limgz) |
---|
493 | end do |
---|
494 | end do |
---|
495 | |
---|
496 | else if (stagger.eq.-1) then ! test dramatique ! |
---|
497 | write(6,*) 'stagger', stagger, iter_beta |
---|
498 | do j=2,ny |
---|
499 | do i=2,nx |
---|
500 | if (abs(uxbar(i,j))-abs(Vcol_x(i,j)).lt.-50.) then |
---|
501 | ! write(6,*) i,j, betamx(i,j),betamx(i,j+1), betamx(i,j-1),uxbar(i,j),Vcol_x(i,j) |
---|
502 | betamx(i,j)=min(betamx(i,j),100.) |
---|
503 | ! write(6,*) i,j, betamx(i,j),betamx(i,j+1), betamx(i,j-1) |
---|
504 | endif |
---|
505 | |
---|
506 | if (abs(uybar(i,j))-abs(Vcol_y(i,j)).lt.-50.) then |
---|
507 | betamy(i,j)=min(betamy(i,j),100.) |
---|
508 | endif |
---|
509 | end do |
---|
510 | end do |
---|
511 | do j=2,ny |
---|
512 | do i=2,nx |
---|
513 | if (betamx(i,j).eq.-1) then |
---|
514 | |
---|
515 | betamx(i,j)=100. |
---|
516 | betamx(i,j+1)=200 |
---|
517 | betamx(i,j-1)=200 |
---|
518 | end if |
---|
519 | if (betamy(i,j).eq.-1) then |
---|
520 | betamy(i,j)=100 |
---|
521 | betamy(i+1,j)=200 |
---|
522 | betamy(i-1,j)=200 |
---|
523 | end if |
---|
524 | end do |
---|
525 | end do |
---|
526 | endif stag |
---|
527 | else |
---|
528 | do j=2,ny |
---|
529 | do i=2,nx |
---|
530 | if (gzmx(i,j)) then |
---|
531 | ecart_quad= ecart_quad & |
---|
532 | + abs(Vcol_x(i,j)*Vcol_x(i,j)-uxbar(i,j)*uxbar(i,j)) |
---|
533 | end if |
---|
534 | if (gzmy(i,j)) then |
---|
535 | ecart_quad= ecart_quad & |
---|
536 | + abs(Vcol_y(i,j)*Vcol_y(i,j)-uybar(i,j)*uybar(i,j)) |
---|
537 | end if |
---|
538 | end do |
---|
539 | end do |
---|
540 | end if |
---|
541 | |
---|
542 | |
---|
543 | write(6,*) 'ecart_quad',ecart_quad |
---|
544 | !-------------------------- fin iteration Arthern like-------------------------- |
---|
545 | end subroutine Arthern_like_iter |
---|
546 | |
---|
547 | |
---|
548 | end module dragging_calc_beta |
---|
549 | |
---|