Context Navigation

source: trunk/SOURCES/Recul_force_grounding_line/toy_retreat_mod.f90 @ 22

Last change on this file since 22 was 22, checked in by roche, 8 years ago
Petites adaptations diverses du code pour compilation en gfortran. Ajout d un Makefile flexible a option pour choisir ifort ou gfortran.
File size: 52.4 KB

Line
1	! Construit les cartes qui servent dans le recul de grounding line
2	! et teste les methodes de recul.
3	! la diminution de H est lineaire avec le recul
4
5	module toy_retreat_mod
6
7	use io_netcdf_grisli
8	use declar_toy_retreat
9
10	implicit none
11
12
13	contains
14
15	!----------------------------------------------------------------------------------
16	! subroutine time_step_recu : ce qui se fait a chaque pas de temps
17	! subroutine init_proto_recul : lit les donnees et initialise
18	! subroutine init_masques : initialise les masques
19	! subroutine pos_gr_line_init : Calcule les positions initiales de Grounding line
20	! subroutine update_retreat_rates : en fonction du temps
21	! subroutine calc_delHdt_maj : calcul le dHdt et le new H sur chaque maille
22	! subroutine calc_H_float : calcule la hauteur de flottaison
23	! subroutine update_flottants : points qui deviennent flottant
24	! subroutine update_a_traiter : points qui deviennent a traiter
25	! subroutine prescribe_Hp : renvoie Hp et I_Hp vers GRISLI
26	!----------------------------------------------------------------------------------
27
28	!----------------------------------------------------------------------------------
29	! time_step_recul : ce qui se fait a chaque pas de temps
30	!----------------------------------------------------------------------------------
31
32	subroutine time_step_recul
33
34	implicit none
35
36	call update_retreat_rates ! update les vitesses de retrait (autorise que si time > time_dep)
37	call calc_eps_max ! pour le sanity check
38	call calc_delHdt_maj ! calcule delHdt et nouveau H pour chaque point
39	call calc_H_float ! calcule la valeur de H_flot
40	call Schoof_flux ! calcul coef_schoof, bf_x_schoof, bf_y_schoof
41	call update_flottants ! update les nouveaux points flottants et xgl, delHdx associes
42	! call update_ramollo ! ramolli eventuellement les ice shelves apres le temps de demarrage
43	call prescribe_Hp ! renvoie Hp et I_Hp vers GRISLI
44	call update_a_traiter ! update le masque des points a traiter
45
46	! quelques sorties
47	debug_3D(:,:,94) = time_float(:,:)
48	debug_3D(:,:,95) = ux_schoof(:,:)
49	debug_3D(:,:,96) = uy_schoof(:,:)
50
51	end subroutine time_step_recul
52
53
54
55	!----------------------------------------------------------------------------------
56	! init_proto_recul : lit les donnees et initialise
57	!----------------------------------------------------------------------------------
58
59	subroutine init_proto_recul
60
61	use declar_toy_retreat
62	implicit none
63
64	namelist/retreat_ice2sea/region_file,file_exp
65
66
67	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
68	428 format(A)
69	read(num_param,retreat_ice2sea)
70
71	write(num_rep_42,428) '!___________________________________________________________'
72	write(num_rep_42,428) '! read parameters of retreat'
73	write(num_rep_42,retreat_ice2sea)
74	write(num_rep_42,428) '! file_exp est le fichier qui contient toutes les experiences'
75	write(num_rep_42,428) '! le numero de l experience correspond au numero du job'
76	write(num_rep_42,428) '!___________________________________________________________'
77
78	call lect_retreat_files
79	! call lect_retreat_nolinfiles ! pour pouvoir lire aussi l'exposant de la loi de frottement
80
81
82	call init_masques ! initialise mk_traiter, retreat0_x, retreat0_y, Deltatot_H0,
83	call pos_gr_line_init ! calcule x_gl_mx, y_gl_my, delHdx_mx, delHdy_my a l'initialisation
84	! ensuite cette position est updated dans la boucle
85
86	call calc_eps_max
87	delHdt_sanity(:,:) = 5000.
88
89	call init_Schoof_flux
90
91
92	end subroutine init_proto_recul
93
94
95
96
97
98	!-------------------------------------------------------------------------
99	!< init_masques : initialise les masques
100	!-------------------------------------------------------------------------
101	subroutine init_masques
102
103	!--------------------------------------------------------------------------
104	! Pour chaque noeud pose dont le socle est sous le niveau des mers
105	!
106	! definition des noeuds a traiter
107	!------------------------------------
108	!
109	! mk_traiter = 2 : pose loin de gl
110	! mk_traiter = 1 : en cours de recul
111	! mk_traiter = -1: flottant implique dans le recul
112	! mk_traiter = -2: flottant loin de GL
113	! mk_traiter = 0 : pas a traiter
114
115	use declar_toy_retreat
116	implicit none
117	integer :: som_voisins !< pour des tests
118
119
120	! definition du masque a traiter
121
122	where ((Mk_gr(:,:).eq.1).and.(Bsoc(:,:).lt.0.)) ! points à traiter
123	mk_traiter(:,:) = 1
124	H_float(:,:) = Bsoc(:,:) / Coef_bflot
125	elsewhere ((Mk_gr(:,:).eq.1).and.(Bsoc(:,:).ge.0.)) ! pose et stable
126	mk_traiter(:,:) = 0
127	H_float(:,:) = -10000.
128	elsewhere (Mk_gr(:,:).eq.0)
129	mk_traiter(:,:) = -1
130	H_float (:,:) = Bsoc(:,:) / Coef_bflot
131	end where
132
133	! calcul de la hauteur de flottaison pour le socle test
134
135	where ((Mk_gr(:,:).eq.1).and.(B_test(:,:).lt.0.)) ! points à traiter
136	H_float_test(:,:) = B_test(:,:) / Coef_bflot
137	elsewhere ((Mk_gr(:,:).eq.1).and.(B_test(:,:).ge.0.)) ! pose et stable
138	H_float_test(:,:) = -10000.
139	elsewhere (Mk_gr(:,:).eq.0)
140	H_float_test(:,:) = B_test(:,:) / Coef_bflot
141	end where
142
143
144
145	do j=2,ny-1
146	do i =2,nx-1
147
148	if (mk_traiter(i,j).eq.1) then ! pose et a traiter
149	if (any(mk_gr(i-1:i+1,j-1:j+1).eq.0)) then
150	mk_traiter(i,j) = 1 ! tout de suite
151	else
152	mk_traiter(i,j) = 2 ! apres propagation
153	end if
154	end if
155
156	if (mk_traiter(i,j).eq.-1) then
157	if (any(mk_gr(i-1:i+1,j-1:j+1).eq.1)) then
158	mk_traiter(i,j) = -1 ! proche de la grounding line
159	else
160	mk_traiter(i,j) = -2 ! loin de la grounding line
161	end if
162	end if
163
164	end do
165	end do
166
167
168	call init_retreat0_alpha
169
170	call init_time_depart
171
172	call update_retreat_rates ! retreat rate depend du temps
173
174	ramollo(:,:) = 0. ! les ice shelves ont la viscosite standard
175
176	mk_traiter0(:,:) = mk_traiter(:,:)
177
178	end subroutine init_masques
179
180
181
182	!---------------------------------------------------------------------------
183	!< pos_gr_line : Calcule les positions initiales de Grounding line
184	!---------------------------------------------------------------------------
185	subroutine pos_gr_line_init
186
187	! Calcul des positions de Grounding line
188	!__________________________________________
189	!
190	! les points sont sur les mailles mineures.
191	! la coordonnee est positive quand le point grounded est à l'Est ou au Nord
192	!
193	! ~--------------x------------o
194	! i-1 i x_gl_mx >0
195	! float grounded
196	! <- West -> East
197	!
198	! si deux points flottants: x_gl_mx = -10 dx
199	! si deux points poses : x_gl_mx = 10 dx
200
201	use declar_toy_retreat
202	implicit none
203
204	! variables locales
205	real :: H_0 ! epaisseur en amont gl
206	real :: H_1 ! epaisseur en aval gl
207	real :: B_0 ! socle en amont gl
208	real :: B_1 ! socle en aval gl
209	real :: dyy ! variable de travail longueur de la maille en diagonale
210
211	dyy = dx(2.*0.5) !longueur de la diagonale
212
213	do j=2,ny-1
214	do i=2,nx-1
215
216	! recherche position sous maille selon x -------------------------------------------
217
218	x_gl: if ((mk_gr(i,j).eq.1).and.(mk_gr(i-1,j).eq.0)) then ! grounded a l'Est
219
220	H_1 = H(i-1,j)
221	H_0 = H(i,j)
222	B_1 = Bsoc(i-1,j)
223	B_0 = Bsoc(i,j)
224
225	call calc_xgl(dx,coef_Bflot,H_0,B_0,H_1,B_1, x_gl_mx(i,j),delHdx_mx(i,j))
226
227	else if ((mk_gr(i,j).eq.0).and.(mk_gr(i-1,j).eq.1)) then ! grounded a l'West
228
229	H_1 = H(i,j)
230	H_0 = H(i-1,j)
231	B_1 = Bsoc(i,j)
232	B_0 = Bsoc(i-1,j)
233
234	call calc_xgl(dx,coef_Bflot,H_0,B_0,H_1,B_1, x_gl_mx(i,j),delHdx_mx(i,j))
235
236	x_gl_mx(i,j) = - x_gl_mx(i,j) ! on change le signe de x_gl
237
238
239	else if ((mk_gr(i,j).eq.0).and.(mk_gr(i-1,j).eq.0)) then ! tout le monde flottant
240
241	x_gl_mx(i,j) = -10.* dx
242	delHdx_mx(i,j) = 0.
243
244
245	else ! tout le monde pose
246
247	x_gl_mx(i,j) = 10.*dx
248	delHdx_mx(i,j) = 0.
249
250	end if x_gl
251
252	! recherche position sous maille selon y -------------------------------------------
253
254	y_gl: if ((mk_gr(i,j).eq.1).and.(mk_gr(i,j-1).eq.0)) then ! grounded au Nord
255
256	H_1 = H(i,j-1)
257	H_0 = H(i,j)
258	B_1 = Bsoc(i,j-1)
259	B_0 = Bsoc(i,j)
260
261	call calc_xgl(dx,coef_Bflot,H_0,B_0,H_1,B_1, y_gl_my(i,j),delHdy_my(i,j))
262
263
264	else if ((mk_gr(i,j).eq.0).and.(mk_gr(i,j-1).eq.1)) then ! grounded au Sud
265
266	H_1 = H(i,j)
267	H_0 = H(i,j-1)
268	B_1 = Bsoc(i,j)
269	B_0 = Bsoc(i,j-1)
270
271	call calc_xgl(dx,coef_Bflot,H_0,B_0,H_1,B_1, y_gl_my(i,j),delHdy_my(i,j))
272
273	y_gl_my(i,j) = - y_gl_my(i,j) ! on change le signe de y_gl
274
275
276	else if ((mk_gr(i,j).eq.0).and.(mk_gr(i,j-1).eq.0)) then ! tout le monde flottant
277
278	y_gl_my(i,j) = -10.*dx
279	delHdy_my(i,j) = 0.
280
281
282	else ! tout le monde pose
283
284	y_gl_my(i,j) = 10.*dx
285	delHdy_my(i,j) = 0.
286
287	end if y_gl
288
289	! recherche position sous maille selon diagonale SW-NE -------------------------------------------
290
291	SW_gl: if ((mk_gr(i,j).eq.1).and.(mk_gr(i-1,j-1).eq.0)) then ! grounded au nord-est
292
293	H_1 = H(i-1,j-1)
294	H_0 = H(i,j)
295	B_1 = Bsoc(i-1,j-1)
296	B_0 = Bsoc(i,j)
297
298	call calc_xgl(dyy,coef_Bflot,H_0,B_0,H_1,B_1, SW_gl_m(i,j),delHdx_SW(i,j))
299
300
301	else if ((mk_gr(i,j).eq.0).and.(mk_gr(i-1,j-1).eq.1)) then ! grounded au Sud-west
302
303	H_1 = H(i,j)
304	H_0 = H(i-1,j-1)
305	B_1 = Bsoc(i,j)
306	B_0 = Bsoc(i-1,j-1)
307
308	call calc_xgl(dyy,coef_Bflot,H_0,B_0,H_1,B_1, SW_gl_m(i,j),delHdx_SW(i,j))
309
310	SW_gl_m(i,j) = - SW_gl_m(i,j) ! on change le signe de SW_gl
311
312
313	else if ((mk_gr(i,j).eq.0).and.(mk_gr(i-1,j-1).eq.0)) then ! tout le monde flottant
314
315	SW_gl_m(i,j) = -10.*dx
316	delHdx_SW(i,j) = 0.
317
318
319	else ! tout le monde pose
320
321	SW_gl_m(i,j) = 10.*dx
322	delHdx_SW(i,j) = 0.
323
324
325	end if SW_gl
326
327	! recherche position sous maille selon diagonale SE-NW -------------------------------------------
328
329
330	SE_gl: if ((mk_gr(i,j).eq.1).and.(mk_gr(i+1,j-1).eq.0)) then ! grounded au nord-west
331
332	H_1 = H(i+1,j-1)
333	H_0 = H(i,j)
334	B_1 = Bsoc(i+1,j-1)
335	B_0 = Bsoc(i,j)
336
337	call calc_xgl(dyy,coef_Bflot,H_0,B_0,H_1,B_1,SE_gl_m(i,j),delHdx_SE(i,j))
338
339
340	else if ((mk_gr(i,j).eq.0).and.(mk_gr(i+1,j-1).eq.1)) then ! grounded au Sud-Est
341
342	H_1 = H(i,j)
343	H_0 = H(i+1,j-1)
344	B_1 = Bsoc(i,j)
345	B_0 = Bsoc(i+1,j-1)
346
347	call calc_xgl(dyy,coef_Bflot,H_0,B_0,H_1,B_1,SE_gl_m(i,j),delHdx_SE(i,j))
348
349
350	SE_gl_m(i,j) = - SE_gl_m(i,j) ! on change le signe de SE_gl
351
352
353	else if ((mk_gr(i,j).eq.0).and.(mk_gr(i+1,j-1).eq.0)) then ! tout le monde flottant
354
355	SE_gl_m(i,j) = -10.*dx
356	delHdx_SE(i,j) = 0.
357
358
359	else ! tout le monde pose
360
361	SE_gl_m(i,j) = 10.*dx
362	delHdx_SE(i,j) = 0.
363
364
365
366	end if SE_gl
367
368
369	end do
370	end do
371
372	return
373	end subroutine pos_gr_line_init
374
375
376
377	!-------------------------------------------------------------------------
378	!< update_retreat_rates : en fonction du temps
379	!-------------------------------------------------------------------------
380
381	subroutine update_retreat_rates
382
383	use declar_toy_retreat
384	implicit none
385
386	call init_retreat0_alpha ! reinitialise retreat0 pour tenir compte des variations beta
387
388	where (time_dep_mx(:,:).lt.time)
389	retreat_x(:,:) = retreat0_x(:,:)
390	elsewhere
391	retreat_x(:,:) = 0.
392	end where
393
394	where (time_dep_my(:,:).lt.time)
395	retreat_y(:,:) = retreat0_y(:,:)
396	elsewhere
397	retreat_y(:,:) = 0.
398	end where
399
400	where (time_dep_SW(:,:).le.time)
401	retreat_SW(:,:) = retreat0_SW(:,:)
402	elsewhere
403	retreat_SW(:,:) = 0.
404	end where
405
406	where (time_dep_SE(:,:).le.time)
407	retreat_SE(:,:) = retreat0_SE(:,:)
408	elsewhere
409	retreat_SE(:,:) = 0.
410	end where
411
412
413
414	return
415	end subroutine update_retreat_rates
416	!------------------------------------------------------------------------
417	!< pour ramollir eventuellement les ice shelves
418	!------------------------------------------------------------------------
419
420	subroutine update_ramollo
421
422	use declar_toy_retreat
423	implicit none
424
425	where ((time_dep(:,:).lt.time).and.(flot(:,:)))
426	ramollo(:,:) = 1.
427	elsewhere
428	ramollo(:,:) = 0.
429	end where
430
431
432	return
433	end subroutine update_ramollo
434
435
436	!-------------------------------------------------------------------------
437	!< calc_delHdt_maj : calcul le dHdt et le new H sur chaque maille
438	!-------------------------------------------------------------------------
439
440	subroutine calc_delHdt_maj
441
442	use declar_toy_retreat
443	implicit none
444
445	real,dimension(2) :: dH_x_voisin ! quand on balaye les voisins
446	real,dimension(2) :: dH_y_voisin ! quand on balaye les voisins
447	real,dimension(2) :: dH_SW_voisin ! quand on balaye les voisins
448	real,dimension(2) :: dH_SE_voisin ! quand on balaye les voisins
449	real :: max_x ! valeur max de dH_x_voisin
450	real :: max_y ! valeur max de dH_y_voisin
451	real :: max_SW ! valeur max de dH_SW_voisin
452	real :: max_SE ! valeur max de dH_SE_voisin
453	integer :: som_voisins !< pour des tests
454
455	! conditions grounding line du point de vue du noeud majeur i
456
457	! ~--------------x--------------o---------------x-----------------~
458	! i-1 i i+1
459	! 0 < x_gl(i) < 10 dx -10 dx < x_gl(i+1) < 0
460
461
462	do j=2,ny-1
463	do i=2,nx-1
464
465	point_a_traiter: if (mk_traiter(i,j).eq.1) then ! Point à traiter
466
467	som_voisins = 0 ! balaye les voisins en croix
468	dH_x_voisin(:) = 0.
469	dH_y_voisin(:) = 0.
470
471	dH_SW_voisin(:) = 0. ! et en diagonale
472	dH_SE_voisin(:) = 0. !
473
474
475	if ((mk_traiter(i-1,j).eq.-1).and. & ! voisin west---------
476	(x_gl_mx(i,j).lt.10.*dx).and.(x_gl_mx(i,j).ge.0.)) then
477
478	dH_x_voisin(1) = delHdx_mx(i,j)*retreat_x(i,j)
479
480	som_voisins = som_voisins + 1
481
482	end if
483
484
485	if ((mk_traiter(i+1,j).eq.-1).and. & ! voisin East---------
486	(x_gl_mx(i+1,j).gt.-10.*dx).and.(x_gl_mx(i+1,j).le.0.)) then
487
488
489	dH_x_voisin(2) = delHdx_mx(i+1,j)*retreat_x(i+1,j)
490	som_voisins = som_voisins + 1
491
492	end if
493
494
495	if ((mk_traiter(i,j-1).eq.-1).and. & ! voisin Sud---------
496	(y_gl_my(i,j).lt.10.*dx).and.(y_gl_my(i,j).ge.0.)) then
497
498	dH_y_voisin(1) = delHdy_my(i,j)*retreat_y(i,j)
499	som_voisins = som_voisins + 1
500
501	end if
502
503	if ((mk_traiter(i,j+1).eq.-1).and. & ! voisin Nord---------
504	(y_gl_my(i,j+1).gt.-10.*dx).and.(y_gl_my(i,j+1).le.0.)) then
505
506	dH_y_voisin(2) = delHdy_my(i,j+1)*retreat_y(i,j+1)
507	som_voisins = som_voisins + 1
508
509	endif
510
511
512	! en diagonale
513
514	if ((mk_traiter(i-1,j-1).eq.-1).and. & ! voisin SW---------
515	(SW_gl_m(i,j).lt.10.*dx).and.(SW_gl_m(i,j).ge.0.)) then
516
517	dH_SW_voisin(1) = delHdx_SW(i,j)*retreat_SW(i,j)
518
519	som_voisins = som_voisins + 1
520
521	end if
522
523
524	if ((mk_traiter(i+1,j+1).eq.-1).and. & ! voisin NE---------
525	(SW_gl_m(i+1,j+1).gt.-10.*dx).and.(SW_gl_m(i+1,j+1).le.0.)) then
526
527	dH_SW_voisin(2) = delHdx_SW(i+1,j+1)*retreat_SW(i+1,j+1)
528	som_voisins = som_voisins + 1
529
530	end if
531
532
533	if ((mk_traiter(i+1,j-1).le.-1).and. & ! voisin SE---------
534	(SE_gl_m(i,j).lt.10.*dx).and.(SE_gl_m(i,j).ge.0.)) then
535
536	dH_SE_voisin(1) = delHdx_SE(i,j)*retreat_SE(i,j)
537	som_voisins = som_voisins + 1
538
539	end if
540
541	if ((mk_traiter(i-1,j+1).le.-1).and. & ! voisin NW---------
542	(SE_gl_m(i-1,j+1).gt.-10.*dx).and.(SE_gl_m(i-1,j+1).le.0.)) then
543
544	dH_SE_voisin(2) = delHdx_SE(i-1,j+1)*retreat_SE(i-1,j+1)
545	som_voisins = som_voisins + 1
546
547	endif
548
549	if (som_voisins.eq.0) then
550	write(6,'("probleme pour i,j=",11(i3,x))'), i,j,mk_traiter(i,j)
551	write(6,*)'j-1', mk_gr(i-1:i+1,j-1)
552	write(6,*)'j ', mk_gr(i-1:i+1,j)
553	write(6,*)'j+1', mk_gr(i-1:i+1,j+1)
554	end if
555
556	! seuls les voisins au dessus peuvent venir
557	! pour les differentes version VERIFER LA LECTURE de B_TEST et B_VOISIN
558
559	! utilise un autre socle que celui du modele pour faire le test
560	! call teste_socle_Btest(i,j,dH_x_voisin,dH_y_voisin,dH_SW_voisin,dH_SE_voisin)
561	call teste_Btest_schoof(i,j,dH_x_voisin,dH_y_voisin,dH_SW_voisin,dH_SE_voisin)
562
563	! utilise les socles min et max pour faire les tests
564	! call teste_socle_Bminmax(i,j,dH_x_voisin,dH_y_voisin,dH_SW_voisin,dH_SE_voisin)
565
566	! utilise le socle du model pour le test
567	! call teste_socle_Bsoc(i,j,dH_x_voisin,dH_y_voisin,dH_SW_voisin,dH_SE_voisin) ! only for tests
568
569	! calcul de la variation d'epaisseur pour ce point
570
571	max_x = maxval(dH_x_voisin)
572	max_y = maxval(dH_y_voisin)
573	max_SW = maxval(dH_SW_voisin)
574	max_SE = maxval(dH_SE_voisin)
575	delHdt(i,j) = max(max_x,max_y,max_SE,max_SW)
576
577	! call sanity_check
578	call sanity_check(i,j)
579	! call no_sanity_check(i,j)
580	delHdt(i,j) = min(delHdt(i,j),delHdt_sanity(i,j))
581
582	! pour tester un recul maxi : le sanity check est aussi utilis comme minimum
583	! delHdt(i,j) = max(delHdt(i,j),delHdt_sanity(i,j)*0.01)
584	! delHdt(i,j) = max(delHdt(i,j),0.)
585
586
587
588	H(i,j) = H(i,j)-delHdt(i,j)*dt ! attention peut etre sous-flottaison
589
590
591	! sanity check
592
593
594
595
596	end if point_a_traiter
597	end do
598	end do
599	end subroutine calc_delHdt_maj
600
601
602
603	!-------------------------------------------------------------------------
604	!< calc_H_float : calcule la hauteur de flottaison
605	!-------------------------------------------------------------------------
606
607	subroutine calc_H_float
608
609	use declar_toy_retreat
610	implicit none
611
612	H_float(:,:) = 0.
613
614	where ((Mk_gr(:,:).gt.0).and.(Bsoc(:,:).lt.0.))
615	H_float(:,:) = Bsoc(:,:) / Coef_bflot
616
617	elsewhere ((Mk_gr(:,:).gt.0).and.(Bsoc(:,:).ge.0.))
618	H_float(:,:) = -10000.
619
620	elsewhere (Mk_gr(:,:).eq.0)
621	H_float(:,:) = Bsoc(:,:) / Coef_bflot
622
623	end where
624
625	end subroutine calc_H_float
626
627
628
629	!-------------------------------------------------------------------------
630	!< update_flottants : points qui deviennent flottant
631	!-------------------------------------------------------------------------
632
633	subroutine update_flottants
634
635	use declar_toy_retreat
636
637	implicit none
638
639	integer :: som_voisins !< pour des tests
640	real :: Hf !< Hf = 0 si socle > 0 sinon H_float
641	real :: dyy ! variable de travail longueur de la maille en diagonale
642
643	dyy = dx(2.*0.5) !longueur de la diagonale
644
645	! update des masques : points qui deviennent flottant
646	do j=2,ny-1
647	do i=2,nx-1
648
649
650
651	new_float: if ((mk_traiter(i,j).eq.1).and.(H(i,j).le.H_float(i,j))) then ! le noeud est passe flottant
652	mk_traiter(i,j) = -1
653	time_float(i,j) = time
654	H(i,j) = H_float(i,j)-20.
655	mk_gr(i,j) = 0
656
657
658	! update vers les noeuds voisins
659
660	if (mk_gr(i-1,j).eq.1) then ! noeud west est grounded
661	x_gl_mx(i,j) = -(dx-1.) ! position gl proche de i,j
662	Hf = max(0.,H_float(i-1,j))
663	delHdx_mx(i,j) = ( H(i-1,j) - Hf )/dx ! taux d'amincissement de (i-1,j)
664	else
665	x_gl_mx(i,j) = -10.*dx
666	delHdx_mx(i,j) = 0.
667	end if
668
669	if (mk_gr(i+1,j).eq.1) then ! noeud Est est grounded
670	x_gl_mx(i+1,j) = dx-1. ! position gl proche de i,j
671	Hf = max(0.,H_float(i+1,j))
672	delHdx_mx(i+1,j) = ( H(i+1,j) - Hf )/dx ! taux d'amincissement de (i+1,j)
673	else
674	x_gl_mx(i+1,j) = -10.*dx
675	delHdx_mx(i+1,j) = 0.
676	end if
677
678	if (mk_gr(i,j-1).eq.1) then ! noeud Sud est grounded
679	y_gl_my(i,j) = -(dx-1.) ! position gl proche de i,j
680	Hf = max(0.,H_float(i,j-1))
681	delHdy_my(i,j) = ( H(i,j-1) - Hf )/dx ! taux d'amincissement de (i-1,j)
682	else
683	y_gl_my(i,j) = -10.*dx
684	delHdy_my(i,j) = 0.
685	end if
686
687	if (mk_gr(i,j+1).eq.1) then ! noeud Nord est grounded
688	y_gl_my(i,j+1) = dx-1. ! position gl proche de i,j
689	Hf = max(0.,H_float(i,j+1))
690	delHdy_my(i,j+1) = ( H(i,j+1) - Hf )/dx ! taux d'amincissement de (i+1,j)
691	else
692	y_gl_my(i,j+1) = -10.*dx
693	delHdy_my(i,j+1) = 0.
694	end if
695
696	! en diagonale
697
698	if (mk_gr(i-1,j-1).eq.1) then ! noeud SW grounded majeur (i-1,j-1)
699	SW_gl_m(i,j) = -(dyy-1.) ! position gl proche de i,j <0
700	Hf = max(0.,H_float(i-1,j-1)) ! noeud mineur i,j
701	delHdx_SW(i,j) = ( H(i-1,j-1) - Hf )/dyy ! taux d'amincissement de (i-1,j-1)
702	else
703	SW_gl_m(i,j) = -10.*dx
704	delHdx_SW(i,j) = 0.
705	end if
706
707	if (mk_gr(i+1,j+1).eq.1) then ! noeud NE grounded majeur (i+1,j+1)
708	SW_gl_m(i+1,j+1) = dyy-1. ! position gl proche de i,j
709	Hf = max(0.,H_float(i+1,j+1)) ! noeud mineur (i+1,j+1)
710	delHdx_SW(i+1,j+1) = ( H(i+1,j+1) - Hf )/dyy ! taux d'amincissement de (i+1,j+1)
711	else
712	SW_gl_m(i+1,j+1) = -10.*dx
713	delHdx_SW(i+1,j+1) = 0.
714	end if
715
716
717	if (mk_gr(i+1,j-1).eq.1) then ! noeud SE grounded majeur (i+1,j-1)
718	SE_gl_m(i,j) = -(dyy-1.) ! position gl proche de i,j <0
719	Hf = max(0.,H_float(i+1,j-1)) ! noeud mineur i,j
720	delHdx_SE(i,j) = ( H(i+1,j-1) - Hf )/dyy ! taux d'amincissement de (i+1,j-1)
721
722	else
723	SE_gl_m(i,j) = -10.*dx
724	delHdx_SE(i,j) = 0.
725	end if
726
727	if (mk_gr(i-1,j+1).eq.1) then ! noeud NW grounded majeur (i-1,j+1)
728	SE_gl_m(i-1,j+1) = dyy-1. ! position gl proche de i,j
729	Hf = max(0.,H_float(i-1,j+1)) ! noeud mineur (i-1,j+1)
730	delHdx_SE(i-1,j+1) = ( H(i-1,j+1) - Hf )/dyy ! taux d'amincissement de (i-1,j+1)
731	else
732	SE_gl_m(i-1,j+1) = -10.*dx
733	delHdx_SE(i-1,j+1) = 0.
734	end if
735
736
737	end if new_float
738	end do
739	end do
740
741	end subroutine update_flottants
742
743
744
745	!-------------------------------------------------------------------------
746	!< update_a_traiter : points qui deviennent a traiter
747	!-------------------------------------------------------------------------
748
749	subroutine update_a_traiter
750
751	use declar_toy_retreat
752	implicit none
753
754	integer :: som_voisins !< pour des tests
755
756	! les nouveaux points a traiter
757	do j=2,ny-1
758	do i =2,nx-1
759
760	! formulation en croix
761	! som_voisins = mk_gr(i-1,j)+mk_gr(i+1,j)+mk_gr(i,j-1)+mk_gr(i,j+1)
762
763	! if ((mk_traiter(i,j).eq.2).and.(som_voisins.gt.0).and.(som_voisins.lt.4)) then
764	! mk_traiter(i,j) = 1
765	! end if
766
767	! Formulation avec diagonales
768	if ((mk_traiter(i,j).eq.2).and.(any(mk_gr(i-1:i+1,j-1:j+1).eq.0))) then
769	mk_traiter(i,j) = 1
770	end if
771
772	end do
773	end do
774	end subroutine update_a_traiter
775
776	!-------------------------------------------------------------------------
777	!< prescribe_Hp : renvoie Hp et I_Hp vers GRISLI
778	!-------------------------------------------------------------------------
779
780	subroutine prescribe_Hp
781
782	use declar_toy_retreat
783	implicit none
784
785	! remet les masques Hp et i_Hp aux valeurs fixes
786
787	i_HP(:,:) = i_Hp0(:,:)
788	Hp(:,:) = Hp0(:,:)
789
790	i=208
791	j=176
792
793	where (mk_gr(:,:).eq.0) ! points flottants
794	i_Hp(:,:) = 1
795	Hp (:,:) = H(:,:) ! les points flottants ne sont pas modifies
796	end where
797
798	where(mk_traiter(:,:).eq.1) ! points imposes Epaisseur en cours de decroissance
799	i_Hp(:,:) = 1
800	Hp (:,:) = H(:,:)
801	end where
802
803	end subroutine prescribe_Hp
804
805	!-------------------------------------------------------------------------
806	!< calc_eps_max : calcule le epsilon max (sanity check)
807	!-------------------------------------------------------------------------
808
809	subroutine calc_eps_max
810
811	use declar_toy_retreat, only:
812	implicit none
813	real :: gamma !< coefficient de flottaison
814
815	gamma = rog(1.+coef_Bflot)
816
817	epsmax(:,:) = Abar(:,:)H(:,:)3(gamma/4.)**3
818
819	end subroutine calc_eps_max
820
821	!-------------------------------------------------------------------------
822	!< calc_xgl : calcule la position sous maille de la grounding line
823	!< en supposant que l'epaisseur varie lineairement
824	!-------------------------------------------------------------------------
825
826	subroutine calc_xgl(dy,alpha,H_0,B_0,H_1,B_1,xpos,delHdx)
827
828	implicit none
829
830	! dummy
831	real,intent(in) :: dy !< longueur de la maille
832	real,intent(in) :: alpha !< coefficient de flottaison = coef_Bflot
833	real,intent(in) :: H_0 !< epaisseur au point pose
834	real,intent(in) :: B_0 !< altitude socle au point pose
835	real,intent(in) :: H_1 !< epaisseur au point flottant
836	real,intent(in) :: B_1 !< altitude socle au point flottant
837	real,intent(out) :: xpos !< position de la ligne (en distance depuis le point pose
838	real,intent(out) :: delHdx !< variation d'epaisseur au noeud pose en fonction d'une variation de xpos
839
840
841	real :: Cgl ! variable de travail
842
843	Cgl = (alpha * (H_1-H_0) - (B_1-B_0))
844
845	if (abs(Cgl).gt.1.e-5) then
846	xpos = dy * (B_0 - alpha * H_0) / Cgl
847	else
848	xpos = dy-1. ! verifier
849	end if
850
851	! la variation d'epaisseur est proportionnelle au recul
852
853	if (xpos.gt.1.) then
854	delHdx = (H_0 - B_0 /alpha) / xpos
855	else
856	delHdx = (H_0 - B_0 /alpha)
857	end if
858
859	end subroutine calc_xgl
860
861	!----------------------------------------------------------------------
862	! lect_retreat_files : lit les fichiers specifiques au retrait
863	!----------------------------------------------------------------------
864	subroutine lect_retreat_files
865
866	use netcdf
867	use io_netcdf_grisli
868	use declar_toy_retreat
869
870	implicit none
871
872	integer :: num ! pour la lecture
873	integer :: num_job ! numero du job
874	integer :: lenrun ! longueur du run name
875	integer :: lenb ! longueur du nom de fichier
876	integer :: pos ! position d'un sous_string
877
878	real*8, dimension(:,:), pointer :: tab !< tableau 2d real ecrit dans le fichier
879	character(len=4) :: char_num ! pour la lecture
880	character(len=80) :: file_B_test !< file name of the bedrock map to test topo instability
881	character(len=80) :: file_voisin !< file name of the bedrock map min or max to test topo instability
882
883
884	! recupere le numero du job dans le runname
885	lenrun = len(trim(runname))
886	char_num = runname(lenrun-3:lenrun)
887
888	read(char_num,*) num_job
889	write(6,*) num_job
890
891	! Nom du fichier d'experience
892	file_exp = trim(dirnameinp)//trim(file_exp)
893
894	! Read the experiment list and keep the line corresponding to the job number
895
896	open(99,file=file_exp)
897	!read(99,*) ! lit la ligne de titre
898
899	do i=1,2000 ! each line is a parameter set
900
901	read(99,*,end=200) num, alpha_lim_1,alpha_lim_2,retreat_1,retreat_2,(time_region(j),j=1,nb_regions),file_B_test
902
903	if (num.eq.num_job) exit ! stop at the job number
904
905	end do
906
907	close(99)
908
909	! read the bedrock map on which we test topo instability
910
911	file_B_test = trim(dirnameinp)//'Socles_Tony/'//file_B_test
912
913	call Read_Ncdf_var('z',trim(file_B_test),tab) ! lit la variable 'z'
914	B_test(:,:) = tab(:,:)
915
916	! lit le fichier region (en general le même que celui des sorties)
917
918	region_file = trim(dirnameinp)//region_file
919	call Read_Ncdf_var('z',region_file,tab) ! lit la variable 'z'
920	map_region(:,:) = nint(tab(:,:))
921
922
923
924	write(6,*) 'runname ',runname,' num',num
925	write(6,*) 'alpha_lim_1,alpha_lim_2,retreat1,retreat2', alpha_lim_1,alpha_lim_2,retreat_1,retreat_2
926	write(6,*) 'time',time_region(:)
927	write(6,*) 'B_test file :', file_B_test
928	!write(6,*) 'B_voisin file :', file_voisin
929	write(6,*) '----------------------------------------------------------'
930
931	return
932	200 write(6,*) ' This job number',num,' runname',runname,' is not in the file',file_exp
933
934	close(99)
935	end subroutine lect_retreat_files
936
937
938	!----------------------------------------------------------------------
939	! init_retreat_rates : calcule les taux de retrait en fonction de divers criteres
940	! ici avec une limite une variation lieaire en alpha
941	!----------------------------------------------------------------------
942	subroutine init_retreat0_alpha
943
944	use declar_toy_retreat
945	implicit none
946	real :: A_alpha !pour calcul retreat = A_alpha * alpha + B_alpha
947	real :: B_alpha !
948
949	A_alpha = (retreat_2 - retreat_1)/(alpha_lim_2 - alpha_lim_1)
950	B_alpha = retreat_1 - A_alpha * alpha_lim_1
951
952
953	! critere sur le alpha (log10 de beta)
954	!------------------------
955
956	travail_centre(:,:) = log10(max(beta_centre(:,:),1.e-5))
957
958	call moyennes_demi_mailles
959
960	where (travail_mx(:,:).gt.alpha_lim_2) !---- en x
961	retreat0_x(:,:) = retreat_2*1000.
962	elsewhere (travail_mx(:,:).lt.alpha_lim_1)
963	retreat0_x(:,:) = retreat_1*1000.
964	elsewhere
965	retreat0_x(:,:) = A_alpha * travail_mx(:,:) + B_alpha
966	end where
967
968	where (travail_my(:,:).gt.alpha_lim_2) !---- en y
969	retreat0_y(:,:) = retreat_2*1000.
970	elsewhere (travail_my(:,:).lt.alpha_lim_1)
971	retreat0_y(:,:) = retreat_1*1000.
972	elsewhere
973	retreat0_y(:,:) = A_alpha * travail_my(:,:) + B_alpha
974	end where
975
976	where (travail_SW(:,:).gt.alpha_lim_2) !---- en SW
977	retreat0_SW(:,:) = retreat_2*1000.
978	elsewhere (travail_SW(:,:).lt.alpha_lim_1)
979	retreat0_SW(:,:) = retreat_1*1000.
980	elsewhere
981	retreat0_SW(:,:) = A_alpha * travail_SW(:,:) + B_alpha
982	end where
983
984	where (travail_SE(:,:).gt.alpha_lim_2) !---- en SE
985	retreat0_SE(:,:) = retreat_2*1000.
986	elsewhere (travail_SE(:,:).lt.alpha_lim_1)
987	retreat0_SE(:,:) = retreat_1*1000.
988	elsewhere
989	retreat0_SE(:,:) = A_alpha * travail_SE(:,:) + B_alpha
990	end where
991
992
993	! enleve les regions qui ne sont pas a traiter
994
995	do j=2,ny-1
996	do i =2,nx-1
997	if (mk_traiter(i,j)*mk_traiter(i-1,j).eq.0) then !---- en x
998	retreat0_x(i,j) = 0.
999	end if
1000	if (mk_traiter(i,j)*mk_traiter(i,j-1).eq.0) then !---- en y
1001	retreat0_y(i,j) = 0.
1002	end if
1003	if (mk_traiter(i,j)*mk_traiter(i-1,j-1).eq.0) then !---- en SW
1004	retreat0_SW(i,j) = 0.
1005	end if
1006	if (mk_traiter(i,j)*mk_traiter(i+1,j-1).eq.0) then !---- en SE
1007	retreat0_SE(i,j) = 0.
1008	end if
1009	end do
1010	end do
1011
1012
1013	end subroutine init_retreat0_alpha
1014
1015
1016	!----------------------------------------------------------------------
1017	! init_retreat_rates : calcule les taux de retrait en fonction de divers criteres
1018	! ici avec une limite abrupte en beta
1019	!----------------------------------------------------------------------
1020	subroutine init_retreat0_beta
1021
1022	use declar_toy_retreat
1023	implicit none
1024
1025	! critere sur le beta
1026	!------------------------
1027
1028	travail_centre(:,:) = beta_centre(:,:)
1029
1030	call moyennes_demi_mailles
1031
1032	where (travail_mx(:,:).gt.beta_lim_ret) !---- en x
1033	retreat0_x(:,:) = retreat_2*1000.
1034	elsewhere
1035	retreat0_x(:,:) = retreat_1*1000.
1036	end where
1037
1038	where (travail_my(:,:).gt.beta_lim_ret) !---- en y
1039	retreat0_y(:,:) = retreat_2*1000.
1040	elsewhere
1041	retreat0_y(:,:) = retreat_1*1000.
1042	end where
1043
1044	where (travail_SW(:,:).gt.beta_lim_ret) !---- en SW
1045	retreat0_SW(:,:) = retreat_2*1000.
1046	elsewhere
1047	retreat0_SW(:,:) = retreat_1*1000.
1048	end where
1049
1050	where (travail_SE(:,:).gt.beta_lim_ret) !---- en SE
1051	retreat0_SE(:,:) = retreat_2*1000.
1052	elsewhere
1053	retreat0_SE(:,:) = retreat_1*1000.
1054	end where
1055
1056	! enleve les regions qui ne sont pas a traiter
1057
1058	do j=2,ny-1
1059	do i =2,nx-1
1060	if (mk_traiter(i,j)*mk_traiter(i-1,j).eq.0) then !---- en x
1061	retreat0_x(i,j) = 0.
1062	end if
1063	if (mk_traiter(i,j)*mk_traiter(i,j-1).eq.0) then !---- en y
1064	retreat0_y(i,j) = 0.
1065	end if
1066	if (mk_traiter(i,j)*mk_traiter(i-1,j-1).eq.0) then !---- en SW
1067	retreat0_SW(i,j) = 0.
1068	end if
1069	if (mk_traiter(i,j)*mk_traiter(i+1,j-1).eq.0) then !---- en SE
1070	retreat0_SE(i,j) = 0.
1071	end if
1072	end do
1073	end do
1074
1075
1076	end subroutine init_retreat0_beta
1077
1078
1079	!----------------------------------------------------------------------
1080	! moyennes_demi_mailles : calcule les moyennes d'un tableau ! travail
1081
1082	!----------------------------------------------------------------------
1083	subroutine moyennes_demi_mailles
1084	use declar_toy_retreat
1085	implicit none
1086	do j=2,ny-1
1087	do i =2,nx-1
1088	travail_mx(i,j) = 0.5* ( travail_centre(i,j) + travail_centre(i-1,j) )
1089	travail_my(i,j) = 0.5* ( travail_centre(i,j) + travail_centre(i,j-1) )
1090	travail_SW(i,j) = 0.5* ( travail_centre(i,j) +travail_centre(i-1,j-1))
1091	travail_SE(i,j) = 0.5* ( travail_centre(i,j) +travail_centre(i+1,j-1))
1092	end do
1093	end do
1094
1095	end subroutine moyennes_demi_mailles
1096
1097	!----------------------------------------------------------------------
1098	! max_demi_mailles : calcule les max d'un tableau ! travail
1099
1100	!----------------------------------------------------------------------
1101	subroutine max_demi_mailles
1102	use declar_toy_retreat
1103	implicit none
1104	do j=2,ny-1
1105	do i =2,nx-1
1106	travail_mx(i,j) = max ( travail_centre(i,j) , travail_centre(i-1,j) )
1107	travail_my(i,j) = max ( travail_centre(i,j) , travail_centre(i,j-1) )
1108	travail_SW(i,j) = max ( travail_centre(i,j) , travail_centre(i-1,j-1))
1109	travail_SE(i,j) = max ( travail_centre(i,j) , travail_centre(i+1,j-1))
1110	end do
1111	end do
1112
1113	end subroutine max_demi_mailles
1114
1115	!----------------------------------------------------------------------
1116	! init_time_depart : donne les temps de depart en fonction des regions
1117	!----------------------------------------------------------------------
1118	! init_time_depart : donne les temps de depart en fonction des regions
1119	!----------------------------------------------------------------------
1120
1121	subroutine init_time_depart
1122
1123
1124	use declar_toy_retreat
1125	implicit none
1126	real :: time_inacc = 10000. ! temps donne aux noeuds inaccessibles
1127	real :: t_already_float = -10. ! temps donne aux noeuds qui flottent deja
1128	real :: time_run_begin = 0. ! debut du run
1129	real :: time_acc = 1000. ! time for the regions below sea level (to see in graphic output)
1130
1131	! region par region
1132
1133	time_region(0) = t_already_float-time_run_begin
1134
1135	do j=2,ny-1
1136	do i=2,nx-1
1137	time_dep(i,j) = (time_region(map_region(i,j)) - time_run_begin)
1138	end do
1139	end do
1140
1141	where (mk_traiter(:,:).eq.0)
1142	time_dep(:,:) = time_inacc
1143	end where
1144
1145	! calcul sur les demi-mailles
1146	travail_centre(:,:) = time_dep(:,:)
1147
1148	! call moyennes_demi_mailles
1149	call max_demi_mailles
1150
1151	time_dep_mx(:,:) = travail_mx(:,:)
1152	time_dep_my(:,:) = travail_my(:,:)
1153	time_dep_SW(:,:) = travail_SW(:,:)
1154	time_dep_SE(:,:) = travail_SE(:,:)
1155
1156
1157	! le temps est tres eleve pour les zones inaccessibles
1158	do j=2,ny-1
1159	do i =2,nx-1
1160	if (mk_traiter(i,j)*mk_traiter(i-1,j).eq.0) then !---- en x
1161	time_dep_mx(i,j) = time_inacc
1162	end if
1163	if (mk_traiter(i,j)*mk_traiter(i,j-1).eq.0) then !---- en y
1164	time_dep_my(i,j) = time_inacc
1165	end if
1166	if (mk_traiter(i,j)*mk_traiter(i-1,j-1).eq.0) then !---- en SW
1167	time_dep_SW(i,j) = time_inacc
1168	end if
1169	if (mk_traiter(i,j)*mk_traiter(i+1,j-1).eq.0) then !---- en SE
1170	time_dep_SE(i,j) = time_inacc
1171	end if
1172	end do
1173	end do
1174
1175	! pour bien voir les régions dans time_float
1176	where(mk_traiter(:,:).eq.2)
1177	time_float(:,:) = time_acc
1178	elsewhere (mk_traiter(:,:).eq.0)
1179	time_float(:,:) = time_inacc
1180	elsewhere (mk_traiter(:,:).lt.0)
1181	time_float(:,:) = t_already_float
1182	end where
1183
1184
1185	end subroutine init_time_depart
1186
1187
1188	!----------------------------------------------------------------------
1189	!< teste_socle_Btest : pour un point i,j : teste si les voisins sont plus haut
1190	!----------------------------------------------------------------------
1191
1192	subroutine teste_socle_Btest(ii,jj,A_x,A_y,A_SW,A_SE)
1193
1194	use declar_toy_retreat
1195	implicit none
1196
1197	integer,intent(in) :: ii,jj ! indice du point considere
1198	real,dimension(2),intent(inout) :: A_x ! quand on balaye les voisins
1199	real,dimension(2),intent(inout) :: A_y ! quand on balaye les voisins
1200	real,dimension(2),intent(inout) :: A_SW ! quand on balaye les voisins
1201	real,dimension(2),intent(inout) :: A_SE ! quand on balaye les voisins
1202
1203
1204
1205	! Btest est plus haut = Bsoc + Bsoc_sigma
1206
1207
1208	! Mais pas suffisamment
1209
1210	if (B_test(ii-1,jj ).le.B_test(ii,jj)) A_x(1) = 0. !---- en x
1211	if (B_test(ii+1,jj ).le.B_test(ii,jj)) A_x(2) = 0.
1212	if (B_test(ii ,jj-1).le.B_test(ii,jj)) A_y(1) = 0. !---- en y
1213	if (B_test(ii ,jj+1).le.B_test(ii,jj)) A_y(2) = 0.
1214
1215	if (B_test(ii-1,jj-1).le.B_test(ii,jj)) A_SW(1) = 0. !---- en SW
1216	if (B_test(ii+1,jj+1).le.B_test(ii,jj)) A_SW(2) = 0.
1217	if (B_test(ii+1,jj-1).le.B_test(ii,jj)) A_SE(1) = 0. !---- en SE
1218	if (B_test(ii-1,jj+1).le.B_test(ii,jj)) A_SE(2) = 0.
1219
1220	end subroutine teste_socle_Btest
1221	!----------------------------------------------------------------------
1222	!< teste_socle_Btest : pour un point i,j : teste si les voisins sont plus haut
1223	!----------------------------------------------------------------------
1224
1225	subroutine teste_Btest_schoof(ii,jj,A_x,A_y,A_SW,A_SE)
1226
1227	use declar_toy_retreat
1228	implicit none
1229
1230	integer,intent(in) :: ii,jj ! indice du point considere
1231	real,dimension(2),intent(inout) :: A_x ! quand on balaye les voisins
1232	real,dimension(2),intent(inout) :: A_y ! quand on balaye les voisins
1233	real,dimension(2),intent(inout) :: A_SW ! quand on balaye les voisins
1234	real,dimension(2),intent(inout) :: A_SE ! quand on balaye les voisins
1235	real :: U2_schoof ! variable de calcul
1236
1237
1238	! test associe sur la pente du socle et le flux de schoof
1239
1240
1241	! le fait de faire le test de schoof sur les vitesses sous-entends
1242	! qu'on fait le flux avec l'epaisseur H_float_test(i,j) ce qui est
1243	! la condiition la plus dure
1244
1245	! si la pente remonte vers le noeud et le flux de Schoof < flux dynamique. Arret.
1246
1247	if ((B_test(ii-1,jj ).le.B_test(ii,jj)) .and. & !---- en x
1248	(Ux_schoof(ii,jj).lt.abs(Uxbar(ii,jj)))) & ! test schoof
1249
1250	A_x(1) = 0. ! noeud mineur_x ii,jj
1251
1252
1253	if ((B_test(ii+1,jj ).le.B_test(ii,jj)) .and. &
1254	(Ux_schoof(ii,jj).lt.abs(Uxbar(ii+1,jj)))) & ! test schoof
1255
1256	A_x(2) = 0. ! noeud mineur_x ii+1
1257
1258
1259	if ((B_test(ii ,jj-1).le.B_test(ii,jj)).and. & !---- en y
1260	(Uy_schoof(ii,jj).lt.abs(Uybar(ii,jj)) )) &
1261
1262	A_y(1) = 0. ! noeud mineur_y ii,jj
1263
1264	if ((B_test(ii ,jj+1).le.B_test(ii,jj)).and. &
1265	(Uy_schoof(ii,jj).lt.abs(Uybar(ii,jj+1)) )) & ! test schoof
1266
1267	A_y(2) = 0. ! noeud mineur_y ii,jj+1
1268
1269	U2_schoof= Ux_schoof(ii,jj)2+Uy_schoof(ii,jj)2
1270
1271	if ((B_test(ii-1,jj-1).le.B_test(ii,jj)) .and. & !---- en SW
1272	(U2_schoof.lt.(Uxbar(ii,jj)2+Uybar(ii,jj)2))) & ! test Schoof
1273
1274
1275	A_SW(1) = 0. ! noeuds mineur_SW (ii,jj)
1276
1277
1278	if ((B_test(ii+1,jj+1).le.B_test(ii,jj)).and. & !---- en NE
1279	(U2_schoof.lt.(Uxbar(ii+1,jj)2+Uybar(ii,jj+1)2))) & ! test Schoof
1280
1281	A_SW(2) = 0. ! noeuds mineurs_SW ii+1,j et ii,jj+1
1282
1283
1284
1285	if ((B_test(ii+1,jj-1).le.B_test(ii,jj)).and. & !---- en SE
1286	(U2_schoof.lt.(Uxbar(ii+1,jj)2+Uybar(ii,jj)2))) & ! test Schoof
1287
1288
1289	A_SE(1) = 0. ! noeuds mineurs_SE ii+1,j et ii,jj
1290
1291
1292
1293	if ((B_test(ii-1,jj+1).le.B_test(ii,jj)).and. & !--- en NW
1294	(U2_schoof.lt.(Uxbar(ii,jj)2+Uybar(ii,jj+1)2))) & ! test Schoof
1295
1296	A_SE(2) = 0. ! noeuds mineurs_SE ii,j et ii,jj+1
1297
1298
1299	end subroutine teste_Btest_schoof
1300
1301	!----------------------------------------------------------------------
1302
1303	!----------------------------------------------------------------------
1304	!< teste_socle_Btest : pour un point i,j : teste si les voisins sont plus haut
1305	!----------------------------------------------------------------------
1306
1307	subroutine teste_socle_Bminmax(ii,jj,A_x,A_y,A_SW,A_SE)
1308
1309	use declar_toy_retreat
1310	implicit none
1311
1312	integer,intent(in) :: ii,jj ! indice du point considere
1313	real,dimension(2),intent(inout) :: A_x ! quand on balaye les voisins
1314	real,dimension(2),intent(inout) :: A_y ! quand on balaye les voisins
1315	real,dimension(2),intent(inout) :: A_SW ! quand on balaye les voisins
1316	real,dimension(2),intent(inout) :: A_SE ! quand on balaye les voisins
1317
1318	! B_voisin est le socle test du point voisin (par ex. le max de la grille 15)
1319	! B_test est le socle test du point considere (par ex. le min de la grille 15)
1320
1321	! Btest est plus haut = Bsoc + Bsoc_sigma
1322
1323
1324	! Mais pas suffisamment
1325
1326	if (B_voisin(ii-1,jj ).le.B_test(ii,jj)) A_x(1) = 0. !---- en x
1327	if (B_voisin(ii+1,jj ).le.B_test(ii,jj)) A_x(2) = 0.
1328	if (B_voisin(ii ,jj-1).le.B_test(ii,jj)) A_y(1) = 0. !---- en y
1329	if (B_voisin(ii ,jj+1).le.B_test(ii,jj)) A_y(2) = 0.
1330
1331	if (B_voisin(ii-1,jj-1).le.B_test(ii,jj)) A_SW(1) = 0. !---- en SW
1332	if (B_voisin(ii+1,jj+1).le.B_test(ii,jj)) A_SW(2) = 0.
1333	if (B_voisin(ii+1,jj-1).le.B_test(ii,jj)) A_SE(1) = 0. !---- en SE
1334	if (B_voisin(ii-1,jj+1).le.B_test(ii,jj)) A_SE(2) = 0.
1335
1336	end subroutine teste_socle_Bminmax
1337
1338
1339	!----------------------------------------------------------------------
1340	! idem previous subroutine but the test is done on Bsoc
1341	! only for tests
1342	!< teste_socle_Btest : pour un point i,j : teste si les voisins sont plus haut
1343	!----------------------------------------------------------------------
1344
1345	subroutine teste_socle_Bsoc(ii,jj,A_x,A_y,A_SW,A_SE)
1346
1347	use declar_toy_retreat
1348	implicit none
1349
1350	integer,intent(in) :: ii,jj ! indice du point considere
1351	real,dimension(2),intent(inout) :: A_x ! quand on balaye les voisins
1352	real,dimension(2),intent(inout) :: A_y ! quand on balaye les voisins
1353	real,dimension(2),intent(inout) :: A_SW ! quand on balaye les voisins
1354	real,dimension(2),intent(inout) :: A_SE ! quand on balaye les voisins
1355
1356
1357
1358	! Btest est plus haut = Bsoc + Bsoc_sigma
1359
1360
1361	! Mais pas suffisamment
1362
1363	if (B_test(ii-1,jj ).le.Bsoc(ii,jj)) A_x(1) = 0. !---- en x
1364	if (B_test(ii+1,jj ).le.Bsoc(ii,jj)) A_x(2) = 0.
1365	if (B_test(ii ,jj-1).le.Bsoc(ii,jj)) A_y(1) = 0. !---- en y
1366	if (B_test(ii ,jj+1).le.Bsoc(ii,jj)) A_y(2) = 0.
1367
1368	if (B_test(ii-1,jj-1).le.Bsoc(ii,jj)) A_SW(1) = 0. !---- en SW
1369	if (B_test(ii+1,jj+1).le.Bsoc(ii,jj)) A_SW(2) = 0.
1370	if (B_test(ii+1,jj-1).le.Bsoc(ii,jj)) A_SE(1) = 0. !---- en SE
1371	if (B_test(ii-1,jj+1).le.Bsoc(ii,jj)) A_SE(2) = 0.
1372
1373	end subroutine teste_socle_Bsoc
1374
1375
1376
1377
1378	!----------------------------------------------------------------------
1379	!< sanity_check : pour un point i,j limite le dhdt
1380	!----------------------------------------------------------------------
1381
1382	subroutine sanity_check(ii,jj)
1383
1384	use declar_toy_retreat
1385	implicit none
1386
1387	integer,intent(in) :: ii,jj ! indice du point considere
1388
1389	real :: dH_x ! variable de travail
1390	real :: dH_y ! variable de travail
1391
1392	DelHdt_sanity(ii,jj) = 0.
1393
1394	! en x -----------------------------------------------------------------------------------
1395	if ((Uxbar(ii,jj).ge.0.).and.(Uxbar(ii+1,jj).ge.0.)) then ! ecoulement --->
1396
1397	dH_x = Uxbar(ii,jj)*(H(ii,jj)-H(ii-1,jj)) / dx
1398	dH_x = Dh_x + epsmax(ii,jj) * H(ii,jj)
1399
1400	else if ((Uxbar(ii,jj).le.0.).and.(Uxbar(ii+1,jj).le.0.)) then ! ecoulement <---
1401
1402	dH_x = Uxbar(ii+1,jj)*(H(ii+1,jj)-H(ii,jj)) / dx
1403	dH_x = Dh_x + epsmax(ii,jj) * H(ii,jj)
1404
1405	else if ((Uxbar(ii,jj).le.0.).and.(Uxbar(ii+1,jj).ge.0.)) then ! ecoulement <------>
1406
1407	dH_x = Dh_x + epsmax(ii,jj) * H(ii,jj)
1408	else ! ecoulement ---><-----
1409	dH_x = Uxbar(ii,jj)(H(ii,jj)-H(ii-1,jj)) / dx + Uxbar(ii+1,jj)(H(ii+1,jj)-H(ii,jj)) / dx
1410	! mais ici pas de terme en epsmax car confine
1411	endif
1412
1413
1414	! en y -----------------------------------------------------------------------------------
1415	if ((Uybar(ii,jj).ge.0.).and.(Uybar(ii,jj+1).ge.0.)) then ! ecoulement ^
1416
1417	dH_y = Uybar(ii,jj)*(H(ii,jj)-H(ii,jj-1)) / dx
1418	dH_y = Dh_y + epsmax(ii,jj) * H(ii,jj)
1419
1420	else if ((Uybar(ii,jj).le.0.).and.(Uybar(ii,jj+1).le.0.)) then ! ecoulement v
1421
1422	dH_y = Uybar(ii,jj+1)*(H(ii,jj+1)-H(ii,jj)) / dx
1423	dH_y = Dh_y + epsmax(ii,jj) * H(ii,jj)
1424
1425	else if ((Uybar(ii,jj).le.0.).and.(Uybar(ii,jj+1).ge.0.)) then ! ecoulement divergent
1426
1427	dH_y = Dh_y + epsmax(ii,jj) * H(ii,jj)
1428	else ! ecoulement convergent
1429	dH_y = Uybar(ii,jj)(H(ii,jj)-H(ii,jj-1)) / dx + Uybar(ii,jj+1)(H(ii,jj+1)-H(ii,jj)) / dx
1430	! mais ici pas de terme en epsmax car confine
1431	endif
1432
1433	DelHdt_sanity(ii,jj) = max(dH_x + dH_y,0.)
1434
1435	end subroutine sanity_check
1436
1437	!-----------------------------------------------------------------------------------
1438	! no_sanity_check pour tester l'influence du sanity check
1439	!-----------------------------------------------------------------------------------
1440
1441	subroutine no_sanity_check(ii,jj)
1442	use declar_toy_retreat
1443	implicit none
1444
1445	integer,intent(in) :: ii,jj ! indice du point considere
1446
1447	real :: dH_x ! variable de travail
1448	real :: dH_y ! variable de travail
1449
1450	DelHdt_sanity(ii,jj) = 5000.
1451	end subroutine no_sanity_check
1452
1453	!--------------------------------------------------------------------
1454	! Schoof_flux : calcule le coefficient et le flux de Schoof
1455	!--------------------------------------------------------------------
1456
1457	subroutine init_Schoof_flux
1458
1459	use declar_toy_retreat
1460	implicit none
1461
1462	nschoof = 3 ! loi de Glen
1463	mschoof = 1./alpha_drag ! exposant loi de friction = alpha_drag m_scoof = 1/alpha_drag
1464	m1_schoof = 1/(mschoof+1)
1465	exp_schoof = (nschoof+3+mschoof)* m1_schoof ! exposant de l'epaisseur
1466
1467
1468	cst_schoof = ( (rog)(nschoof+1) (1+coef_Bflot)nschoof / 4nschoof ) **(m1_schoof) ! formulation no lin
1469
1470	end subroutine Init_Schoof_flux
1471
1472
1473	!--------------------------------------------------------------------
1474	! Schoof_flux : calcule le coefficient et le flux de Schoof
1475	!--------------------------------------------------------------------
1476
1477	subroutine Schoof_flux
1478
1479	use declar_toy_retreat
1480	implicit none
1481
1482	where (H_float_test(:,:).gt.0.) ! epaisseur de flottaison avec le socle test
1483
1484	! coef_schoof(:,:) = (Abar(:,:)/max(1.,beta_centre(:,:)))**m1_schoof
1485	! coef_schoof(:,:) = min(1000.,beta_centre(:,:))
1486
1487	! coef_schoof(:,:) = beta_centre(:,:) ! seulement vrai pour dragging lineaire
1488
1489	coef_schoof(:,:) = coef_drag(:,:) ! coef_drag calcule dans init_dragging
1490
1491
1492	coef_schoof(:,:) = (Abar(:,:)/max(1.,coef_schoof(:,:)))**m1_schoof
1493
1494	coef_schoof(:,:) = coef_schoof(:,:) * cst_schoof * H_float_test(:,:)**exp_schoof
1495	tauf_schoof(:,:) = 0.5 * ro * g * H_float_test(:,:) * (1+coef_Bflot)
1496
1497	elsewhere
1498	coef_schoof(:,:) = 0.
1499	tauf_schoof(:,:) = 0.
1500	end where
1501	debug_3D(:,:,105) = coef_drag(:,:)
1502	! calcul du terme de confinement. A revoir
1503
1504	!!$do j=2,ny-1
1505	!!$ do i=2,nx-1
1506	!!$ if (H_float_test(i,j).gt.0.) then
1507	!!$
1508	!!$ ! 2 lignes en dessous, H et pas H_float_test parce que c'est pour tirer la viscosite moyenne
1509	!!$ if (H(i,j).gt.0.) then
1510	!!$ bf_x_schoof(i,j) = 2.(Uxbar(i+1,j)-Uxbar(i,j)) /dx pvi(i,j) /H(i,j) ! tau'xx,
1511	!!$ bf_y_schoof(i,j) = 2.(Uybar(i,j+1)-Uybar(i,j)) /dx pvi(i,j) /H(i,j) ! tau'yy
1512	!!$ end if
1513	!!$
1514	!!$ if (tauf_schoof(i,j).gt.0.) then
1515	!!$ bf_x_schoof(i,j) = abs(bf_x_schoof(i,j))/tauf_schoof(i,j) ! rapport tauxx/tauf
1516	!!$ bf_y_schoof(i,j) = abs(bf_y_schoof(i,j))/tauf_schoof(i,j) ! rapport tauyy/tauf
1517	!!$ else
1518	!!$ bf_x_schoof(i,j) = 1.
1519	!!$ bf_y_schoof(i,j) = 1.
1520	!!$ end if
1521	!!$
1522	!!$ if (bf_x_schoof(i,j).gt.0) then ! a la puissance nschoof*m1_schoof
1523	!!$ bf_x_schoof(i,j) = bf_x_schoof(i,j) *(nschoofm1_schoof)
1524	!!$ else
1525	!!$ bf_x_schoof(i,j) = 0.
1526	!!$ end if
1527	!!$
1528	!!$ if (bf_y_schoof(i,j).gt.0) then ! a la puissance nschoof*m1_schoof
1529	!!$ bf_y_schoof(i,j) = bf_y_schoof(i,j) *(nschoofm1_schoof)
1530	!!$ else
1531	!!$ bf_y_schoof(i,j) = 0.
1532	!!$ end if
1533	!!$ else ! pas de flottaison
1534	!!$ bf_x_schoof(i,j) = 0.
1535	!!$ bf_y_schoof(i,j) = 0.
1536	!!$ end if
1537	!!$
1538	!!$ end do
1539	!!$end do
1540
1541	bf_x_schoof(:,:) = 1.
1542	bf_y_schoof(:,:) = 1.
1543
1544	where ((H_float_test(:,:).gt.0.).and.(H(:,:).gt.2.))
1545	Ux_schoof(:,:) = bf_x_schoof(:,:)*coef_schoof(:,:)/H_float_test(:,:)
1546	Uy_schoof(:,:) = bf_y_schoof(:,:)*coef_schoof(:,:)/H_float_test(:,:)
1547	elsewhere
1548	Ux_schoof(:,:) = 0.
1549	Uy_schoof(:,:) = 0.
1550	end where
1551
1552	end subroutine Schoof_flux
1553
1554
1555
1556	end module toy_retreat_mod
1557	!----------------------------------------------------------------------------------
1558
1559
1560
1561
1562

Note: See TracBrowser for help on using the repository browser.

Download in other formats: