source: trunk/SOURCES/Recul_force_grounding_line/proto_recul_mod.f90 @ 4

Last change on this file since 4 was 4, checked in by dumas, 10 years ago

initial import GRISLI trunk
File size: 36.7 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
5module proto_recul_mod
6
7use io_netcdf_grisli
8use declar_proto_recul
9implicit none
10
11
12contains   
13
14!----------------------------------------------------------------------------------
15! subroutine  time_step_recu       : ce qui se fait a chaque pas de temps
16! subroutine  init_proto_recul     : lit les donnees et initialise
17! subroutine  init_masques         : initialise les masques
18! subroutine  pos_gr_line_init     : Calcule les positions initiales de Grounding line
19! subroutine  update_retreat_rates : en fonction du temps
20! subroutine  calc_delHdt_maj      : calcul le dHdt et le new H sur chaque maille
21! subroutine  calc_H_float         : calcule la hauteur de flottaison
22! subroutine  update_flottants     :  points qui deviennent flottant
23! subroutine  update_a_traiter     :  points qui deviennent a traiter
24! subroutine  prescribe_Hp         :  renvoie Hp et I_Hp vers GRISLI
25!----------------------------------------------------------------------------------
26
27!----------------------------------------------------------------------------------
28! time_step_recul : ce qui se fait a chaque pas de temps
29!----------------------------------------------------------------------------------
30
31subroutine time_step_recul
32
33implicit none
34
35  call update_retreat_rates      !  update les vitesses de retrait (autorise que si time > time_dep)
36  Btest(:,:) = Bsoc(:,:) + Bsoc_sigma(:,:)
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 update_flottants          !  update les nouveaux points flottants et xgl, delHdx associes
41  call prescribe_Hp              !  renvoie Hp et I_Hp vers GRISLI
42  call update_a_traiter          !  update le masque des points a traiter
43
44end subroutine time_step_recul
45
46
47
48!----------------------------------------------------------------------------------
49! init_proto_recul : lit les donnees et initialise
50!----------------------------------------------------------------------------------
51
52subroutine init_proto_recul
53
54use declar_proto_recul
55implicit none
56
57
58 namelist/retreat/number_data_B_file,region_file,file_exp,bed_sigma_1,bed_sigma_2
59
60
61    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
62428 format(A)
63    read(num_param,retreat)
64
65    write(num_rep_42,428) '!___________________________________________________________' 
66    write(num_rep_42,428) '!  read parameters of retreat'
67    write(num_rep_42,retreat)                   
68    write(num_rep_42,428) '!___________________________________________________________' 
69
70  call lect_retreat_files
71
72
73  call init_masques                     ! initialise mk_traiter, retreat0_x, retreat0_y, Deltatot_H0,
74  call pos_gr_line_init                 ! calcule x_gl_mx,  y_gl_my, delHdx_mx, delHdy_my a l'initialisation
75                                        ! ensuite cette position est updated dans la boucle
76
77  call calc_eps_max
78  delHdt_sanity(:,:) = 5000.
79
80
81end subroutine init_proto_recul
82
83
84
85
86
87!-------------------------------------------------------------------------
88!<  init_masques : initialise les masques
89!-------------------------------------------------------------------------
90subroutine init_masques
91
92!--------------------------------------------------------------------------
93! Pour chaque noeud pose dont le socle est sous le niveau des mers
94!
95! definition des noeuds a traiter
96!------------------------------------
97!
98! mk_traiter = 2 : pose loin de gl
99! mk_traiter = 1 : en cours de recul
100! mk_traiter = -1: flottant implique dans le recul
101! mk_traiter = -2: flottant loin de GL
102! mk_traiter = 0 : pas a traiter
103
104  use declar_proto_recul
105  implicit none
106  integer                           :: som_voisins  !< pour des tests
107
108
109! definition du masque a traiter
110
111  where ((Mk_gr(:,:).eq.1).and.(Bsoc(:,:).lt.0.))      ! points à traiter
112     mk_traiter(:,:) = 1
113     H_float(:,:)    = Bsoc(:,:) / Coef_bflot
114  elsewhere ((Mk_gr(:,:).eq.1).and.(Bsoc(:,:).ge.0.))  ! pose et stable
115     mk_traiter(:,:) = 0
116     H_float(:,:)     = -10000. 
117  elsewhere (Mk_gr(:,:).eq.0)                         
118     mk_traiter(:,:) = -1
119     H_float   (:,:) = Bsoc(:,:) / Coef_bflot
120  end where
121
122
123  do j=2,ny-1
124     do i =2,nx-1 
125
126        if (mk_traiter(i,j).eq.1) then              ! pose et a traiter
127           if (any(mk_gr(i-1:i+1,j-1:j+1).eq.0)) then
128              mk_traiter(i,j) = 1                   ! tout de suite
129           else
130              mk_traiter(i,j) = 2                   ! apres propagation
131           end if
132        end if
133
134        if (mk_traiter(i,j).eq.-1) then
135           if (any(mk_gr(i-1:i+1,j-1:j+1).eq.1)) then
136              mk_traiter(i,j) = -1                  ! proche de la grounding line
137           else
138              mk_traiter(i,j) = -2                  ! loin de la grounding line
139           end if
140        end if
141   
142     end do
143  end do
144
145
146call init_retreat0
147
148call init_time_depart
149
150call update_retreat_rates                  ! retreat rate depend du temps
151
152call bedrock_sigma
153
154
155
156
157
158mk_traiter0(:,:) = mk_traiter(:,:)
159
160end subroutine init_masques
161
162
163
164!---------------------------------------------------------------------------
165!< pos_gr_line : Calcule les positions initiales de Grounding line
166!---------------------------------------------------------------------------
167subroutine pos_gr_line_init
168
169! Calcul des positions de Grounding line
170!__________________________________________
171!
172! les points sont sur les mailles mineures.
173! la coordonnee est positive quand le point grounded est à l'Est ou au Nord
174!
175!         ~--------------x------------o           
176!         i-1                         i                 x_gl_mx >0
177!        float                      grounded
178!  <-  West                               -> East     
179!
180! si deux points flottants:  x_gl_mx = -10 dx
181! si deux points poses :      x_gl_mx = 10 dx
182
183  use declar_proto_recul
184  implicit none
185
186  ! variables locales
187  real      :: H_0         ! epaisseur en amont gl
188  real      :: H_1         ! epaisseur en aval gl
189  real      :: B_0         ! socle en amont gl
190  real      :: B_1         ! socle en aval gl
191  real      :: dyy         ! variable de travail  longueur de la maille en diagonale       
192
193  dyy = dx*(2.**0.5)       !longueur de la diagonale
194
195  do j=2,ny-1
196     do i=2,nx-1   
197     
198!     recherche position sous maille  selon x -------------------------------------------
199
200        x_gl: if ((mk_gr(i,j).eq.1).and.(mk_gr(i-1,j).eq.0)) then      ! grounded a l'Est
201
202           H_1 = H(i-1,j)
203           H_0 = H(i,j)
204           B_1 = Bsoc(i-1,j)
205           B_0 = Bsoc(i,j)
206
207           call calc_xgl(dx,coef_Bflot,H_0,B_0,H_1,B_1, x_gl_mx(i,j),delHdx_mx(i,j)) 
208
209        else if ((mk_gr(i,j).eq.0).and.(mk_gr(i-1,j).eq.1)) then       ! grounded a l'West
210
211           H_1 = H(i,j)
212           H_0 = H(i-1,j)
213           B_1 = Bsoc(i,j)
214           B_0 = Bsoc(i-1,j)
215
216           call calc_xgl(dx,coef_Bflot,H_0,B_0,H_1,B_1, x_gl_mx(i,j),delHdx_mx(i,j)) 
217
218           x_gl_mx(i,j)   = - x_gl_mx(i,j)                            ! on change le signe de x_gl
219         
220
221        else if ((mk_gr(i,j).eq.0).and.(mk_gr(i-1,j).eq.0)) then      ! tout le monde flottant
222
223           x_gl_mx(i,j)   = -10.* dx
224           delHdx_mx(i,j) =  0.
225
226
227        else                                                          ! tout le monde pose
228
229           x_gl_mx(i,j)   = 10.*dx
230           delHdx_mx(i,j) = 0.
231
232        end if x_gl
233
234!     recherche position sous maille  selon y -------------------------------------------
235
236        y_gl: if ((mk_gr(i,j).eq.1).and.(mk_gr(i,j-1).eq.0)) then      ! grounded au Nord
237
238           H_1 = H(i,j-1)
239           H_0 = H(i,j)
240           B_1 = Bsoc(i,j-1)
241           B_0 = Bsoc(i,j)
242
243           call calc_xgl(dx,coef_Bflot,H_0,B_0,H_1,B_1, y_gl_my(i,j),delHdy_my(i,j)) 
244
245
246        else if ((mk_gr(i,j).eq.0).and.(mk_gr(i,j-1).eq.1)) then       ! grounded au Sud
247
248           H_1 = H(i,j)
249           H_0 = H(i,j-1)
250           B_1 = Bsoc(i,j)
251           B_0 = Bsoc(i,j-1)
252
253           call calc_xgl(dx,coef_Bflot,H_0,B_0,H_1,B_1, y_gl_my(i,j),delHdy_my(i,j)) 
254
255           y_gl_my(i,j)   = - y_gl_my(i,j)                            ! on change le signe de y_gl
256         
257
258        else if ((mk_gr(i,j).eq.0).and.(mk_gr(i,j-1).eq.0)) then      ! tout le monde flottant
259
260           y_gl_my(i,j)   = -10.*dx
261           delHdy_my(i,j) = 0.
262
263
264        else                                                          ! tout le monde pose
265
266           y_gl_my(i,j)   = 10.*dx
267           delHdy_my(i,j) = 0.
268
269        end if y_gl
270
271!     recherche position sous maille  selon diagonale SW-NE -------------------------------------------
272
273        SW_gl: if ((mk_gr(i,j).eq.1).and.(mk_gr(i-1,j-1).eq.0)) then      ! grounded au nord-est
274
275           H_1 = H(i-1,j-1)
276           H_0 = H(i,j)
277           B_1 = Bsoc(i-1,j-1)
278           B_0 = Bsoc(i,j)
279
280           call calc_xgl(dyy,coef_Bflot,H_0,B_0,H_1,B_1, SW_gl_m(i,j),delHdx_SW(i,j)) 
281
282
283        else if ((mk_gr(i,j).eq.0).and.(mk_gr(i-1,j-1).eq.1)) then       ! grounded au Sud-west
284
285           H_1 = H(i,j)
286           H_0 = H(i-1,j-1)
287           B_1 = Bsoc(i,j)
288           B_0 = Bsoc(i-1,j-1)
289
290           call calc_xgl(dyy,coef_Bflot,H_0,B_0,H_1,B_1, SW_gl_m(i,j),delHdx_SW(i,j)) 
291
292           SW_gl_m(i,j) = - SW_gl_m(i,j)                                 ! on change le signe de SW_gl
293         
294
295        else if ((mk_gr(i,j).eq.0).and.(mk_gr(i-1,j-1).eq.0)) then      ! tout le monde flottant
296
297           SW_gl_m(i,j)   = -10.*dx
298           delHdx_SW(i,j) = 0.
299
300
301        else                                                          ! tout le monde pose
302
303           SW_gl_m(i,j)   = 10.*dx
304           delHdx_SW(i,j) = 0.
305
306
307        end if SW_gl
308
309!     recherche position sous maille  selon diagonale SE-NW -------------------------------------------
310
311
312        SE_gl: if ((mk_gr(i,j).eq.1).and.(mk_gr(i+1,j-1).eq.0)) then      ! grounded au nord-west
313
314           H_1 = H(i+1,j-1)
315           H_0 = H(i,j)
316           B_1 = Bsoc(i+1,j-1)
317           B_0 = Bsoc(i,j)
318
319           call calc_xgl(dyy,coef_Bflot,H_0,B_0,H_1,B_1,SE_gl_m(i,j),delHdx_SE(i,j)) 
320
321
322        else if ((mk_gr(i,j).eq.0).and.(mk_gr(i+1,j-1).eq.1)) then       ! grounded au Sud-Est
323
324           H_1 = H(i,j)
325           H_0 = H(i+1,j-1)
326           B_1 = Bsoc(i,j)
327           B_0 = Bsoc(i+1,j-1)
328
329           call calc_xgl(dyy,coef_Bflot,H_0,B_0,H_1,B_1,SE_gl_m(i,j),delHdx_SE(i,j)) 
330         
331
332           SE_gl_m(i,j) = - SE_gl_m(i,j)                                 ! on change le signe de SE_gl
333         
334
335        else if ((mk_gr(i,j).eq.0).and.(mk_gr(i+1,j-1).eq.0)) then       ! tout le monde flottant
336
337           SE_gl_m(i,j)   = -10.*dx
338           delHdx_SE(i,j) = 0.
339
340
341        else                                                              ! tout le monde pose
342
343           SE_gl_m(i,j)   = 10.*dx
344           delHdx_SE(i,j) = 0.
345
346
347
348        end if SE_gl
349
350     
351     end do
352  end do
353
354  return
355end subroutine pos_gr_line_init
356
357
358
359!-------------------------------------------------------------------------
360!<  update_retreat_rates  : en fonction du temps
361!-------------------------------------------------------------------------
362
363subroutine update_retreat_rates
364
365  use declar_proto_recul
366  implicit none
367
368
369  where (time_dep_mx(:,:).lt.time)                                     
370     retreat_x(:,:) =  retreat0_x(:,:)
371  elsewhere
372      retreat_x(:,:) = 0.
373   end where
374
375   where (time_dep_my(:,:).lt.time)                                     
376      retreat_y(:,:) =  retreat0_y(:,:)
377   elsewhere
378      retreat_y(:,:) = 0.
379   end where
380
381  where (time_dep_SW(:,:).le.time)                                     
382     retreat_SW(:,:) =  retreat0_SW(:,:)
383  elsewhere
384     retreat_SW(:,:) = 0.
385  end where
386
387  where (time_dep_SE(:,:).le.time)                                     
388     retreat_SE(:,:) =  retreat0_SE(:,:)
389  elsewhere
390     retreat_SE(:,:) = 0.
391  end where
392
393
394
395
396   return
397 end subroutine update_retreat_rates
398
399
400
401!-------------------------------------------------------------------------
402!<  calc_delHdt_maj  : calcul le dHdt et le new H sur chaque maille
403!-------------------------------------------------------------------------
404
405subroutine calc_delHdt_maj
406
407  use declar_proto_recul
408  implicit none
409
410  real,dimension(2)                  :: dH_x_voisin   ! quand on balaye les voisins
411  real,dimension(2)                  :: dH_y_voisin   ! quand on balaye les voisins
412  real,dimension(2)                  :: dH_SW_voisin  ! quand on balaye les voisins
413  real,dimension(2)                  :: dH_SE_voisin  ! quand on balaye les voisins
414  real                               :: max_x         ! valeur max de  dH_x_voisin
415  real                               :: max_y         ! valeur max de  dH_y_voisin
416  real                               :: max_SW        ! valeur max de  dH_SW_voisin
417  real                               :: max_SE        ! valeur max de  dH_SE_voisin
418  integer                            :: som_voisins   !< pour des tests
419
420! conditions grounding line du point de vue du noeud majeur i
421
422!  ~--------------x--------------o---------------x-----------------~
423!  i-1                           i                                i+1
424!         0 < x_gl(i) < 10 dx          -10 dx < x_gl(i+1) < 0
425
426 
427  do j=2,ny-1
428     do i=2,nx-1
429       
430        point_a_traiter:   if (mk_traiter(i,j).eq.1) then                ! Point à traiter
431
432           som_voisins     = 0                                           ! balaye les voisins en croix
433           dH_x_voisin(:)  = 0.                                                         
434           dH_y_voisin(:)  = 0.                                         
435
436           dH_SW_voisin(:) = 0.                                          ! et en diagonale
437           dH_SE_voisin(:) = 0.                                          !
438
439
440           if ((mk_traiter(i-1,j).eq.-1).and.                      &     ! voisin west---------
441                (x_gl_mx(i,j).lt.10.*dx).and.(x_gl_mx(i,j).ge.0.)) then
442             
443              dH_x_voisin(1) = delHdx_mx(i,j)*retreat_x(i,j) 
444
445              som_voisins  =  som_voisins + 1
446             
447           end if
448
449
450           if ((mk_traiter(i+1,j).eq.-1).and.                      &     ! voisin East---------
451                (x_gl_mx(i+1,j).gt.-10.*dx).and.(x_gl_mx(i+1,j).le.0.)) then
452
453               
454              dH_x_voisin(2) = delHdx_mx(i+1,j)*retreat_x(i+1,j)
455              som_voisins  =  som_voisins + 1
456
457           end if
458
459
460           if ((mk_traiter(i,j-1).eq.-1).and.                      &     ! voisin Sud---------
461                (y_gl_my(i,j).lt.10.*dx).and.(y_gl_my(i,j).ge.0.)) then
462
463              dH_y_voisin(1) = delHdy_my(i,j)*retreat_y(i,j)
464              som_voisins  =  som_voisins + 1
465
466           end if
467
468           if ((mk_traiter(i,j+1).eq.-1).and.                      &     ! voisin Nord---------
469                (y_gl_my(i,j+1).gt.-10.*dx).and.(y_gl_my(i,j+1).le.0.)) then
470
471              dH_y_voisin(2) = delHdy_my(i,j+1)*retreat_y(i,j+1)
472              som_voisins  =  som_voisins + 1
473
474           endif
475
476
477! en diagonale
478
479           if ((mk_traiter(i-1,j-1).eq.-1).and.                     &     ! voisin SW---------
480                (SW_gl_m(i,j).lt.10.*dx).and.(SW_gl_m(i,j).ge.0.)) then
481             
482              dH_SW_voisin(1) = delHdx_SW(i,j)*retreat_SW(i,j) 
483
484              som_voisins  =  som_voisins + 1
485             
486           end if
487
488
489           if ((mk_traiter(i+1,j+1).eq.-1).and.                     &     ! voisin NE---------
490                (SW_gl_m(i+1,j+1).gt.-10.*dx).and.(SW_gl_m(i+1,j+1).le.0.)) then
491
492              dH_SW_voisin(2) = delHdx_SW(i+1,j+1)*retreat_SW(i+1,j+1)
493              som_voisins  =  som_voisins + 1
494
495           end if
496
497
498           if ((mk_traiter(i+1,j-1).le.-1).and.                     &     ! voisin SE---------
499                (SE_gl_m(i,j).lt.10.*dx).and.(SE_gl_m(i,j).ge.0.)) then
500
501              dH_SE_voisin(1) = delHdx_SE(i,j)*retreat_SE(i,j)
502              som_voisins  =  som_voisins + 1
503
504           end if
505
506           if ((mk_traiter(i-1,j+1).le.-1).and.                     &     ! voisin NW---------
507                (SE_gl_m(i-1,j+1).gt.-10.*dx).and.(SE_gl_m(i-1,j+1).le.0.)) then
508
509              dH_SE_voisin(2) = delHdx_SE(i-1,j+1)*retreat_SE(i-1,j+1)
510              som_voisins  =  som_voisins + 1
511
512           endif
513
514           if (som_voisins.eq.0) then
515              write(6,'("probleme pour i,j=",11(i3,x))'), i,j,mk_traiter(i,j)
516              write(6,*)'j-1',     mk_gr(i-1:i+1,j-1)
517              write(6,*)'j  ',     mk_gr(i-1:i+1,j)
518              write(6,*)'j+1',     mk_gr(i-1:i+1,j+1)
519           end if
520
521!     seuls les voisins au dessus peuvent venir
522           call teste_socle_voisins(i,j,dH_x_voisin,dH_y_voisin,dH_SW_voisin,dH_SE_voisin)
523
524
525!     calcul de la variation d'epaisseur pour ce point         
526
527           max_x  = maxval(dH_x_voisin)
528           max_y  = maxval(dH_y_voisin)
529           max_SW = maxval(dH_SW_voisin)
530           max_SE = maxval(dH_SE_voisin)
531           delHdt(i,j) = max(max_x,max_y,max_SE,max_SW) 
532
533! call sanity_check
534           call sanity_check(i,j)
535           delHdt(i,j) = min(delHdt(i,j),delHdt_sanity(i,j))
536
537           H(i,j) = H(i,j)-delHdt(i,j)*dt                            ! attention peut etre sous-flottaison         
538
539
540!      sanity check
541
542
543
544
545        end if point_a_traiter
546     end do
547  end do
548end subroutine calc_delHdt_maj
549
550
551
552!-------------------------------------------------------------------------
553!<  calc_H_float : calcule la hauteur de flottaison
554!-------------------------------------------------------------------------
555
556subroutine calc_H_float
557
558  use declar_proto_recul
559  implicit none
560
561  H_float(:,:) = 0.
562
563  where ((Mk_gr(:,:).gt.0).and.(Bsoc(:,:).lt.0.))     
564     H_float(:,:)    = Bsoc(:,:) / Coef_bflot
565
566  elsewhere ((Mk_gr(:,:).gt.0).and.(Bsoc(:,:).ge.0.))
567     H_float(:,:) = -10000.
568
569  elsewhere (Mk_gr(:,:).eq.0)
570     H_float(:,:)    = Bsoc(:,:) / Coef_bflot
571
572  end where
573
574end subroutine calc_H_float
575
576
577
578!-------------------------------------------------------------------------
579!<  update_flottants :  points qui deviennent flottant
580!-------------------------------------------------------------------------
581
582subroutine update_flottants
583
584  use declar_proto_recul
585
586  implicit none
587
588  integer                           :: som_voisins  !< pour des tests
589  real                              :: Hf           !< Hf = 0 si socle > 0 sinon H_float
590  real                              :: dyy         ! variable de travail  longueur de la maille en diagonale       
591
592  dyy = dx*(2.**0.5)       !longueur de la diagonale
593
594! update des masques : points qui deviennent flottant
595  do j=2,ny-1
596     do i=2,nx-1
597
598
599
600        new_float: if ((mk_traiter(i,j).eq.1).and.(H(i,j).le.H_float(i,j))) then    ! le noeud est passe flottant
601           mk_traiter(i,j) = -1
602           time_float(i,j) = time
603           H(i,j)          = H_float(i,j)-20.
604           mk_gr(i,j)      = 0
605
606
607           ! update vers les noeuds voisins
608
609           if (mk_gr(i-1,j).eq.1)  then                              ! noeud west est grounded
610              x_gl_mx(i,j)     = -(dx-1.)                            ! position gl proche de i,j
611              Hf               = max(0.,H_float(i-1,j))
612              delHdx_mx(i,j)   = ( H(i-1,j) - Hf )/dx                ! taux d'amincissement de (i-1,j)
613           else
614              x_gl_mx(i,j)     = -10.*dx
615              delHdx_mx(i,j)   =   0.
616           end if
617
618           if (mk_gr(i+1,j).eq.1)  then                             ! noeud Est est grounded
619              x_gl_mx(i+1,j)   = dx-1.                              ! position gl proche de i,j
620              Hf               = max(0.,H_float(i+1,j))
621              delHdx_mx(i+1,j) = ( H(i+1,j) - Hf )/dx               ! taux d'amincissement de (i+1,j)
622           else
623              x_gl_mx(i+1,j)   = -10.*dx
624              delHdx_mx(i+1,j) =   0.
625           end if
626
627           if (mk_gr(i,j-1).eq.1)  then                              ! noeud Sud est grounded
628              y_gl_my(i,j)     = -(dx-1.)                            ! position gl proche de i,j
629              Hf               = max(0.,H_float(i,j-1))
630              delHdy_my(i,j)   = ( H(i,j-1) - Hf )/dx                ! taux d'amincissement de (i-1,j)
631           else
632              y_gl_my(i,j)     = -10.*dx
633              delHdy_my(i,j)   =   0.
634           end if
635
636           if (mk_gr(i,j+1).eq.1)  then                              ! noeud Nord est grounded
637              y_gl_my(i,j+1)   = dx-1.                               ! position gl proche de i,j
638              Hf               = max(0.,H_float(i,j+1))
639              delHdy_my(i,j+1) = ( H(i,j+1) - Hf )/dx                ! taux d'amincissement de (i+1,j)
640           else
641              y_gl_my(i,j+1)   = -10.*dx
642              delHdy_my(i,j+1) =   0.
643           end if
644
645! en diagonale
646 
647           if (mk_gr(i-1,j-1).eq.1)  then                            ! noeud SW grounded majeur (i-1,j-1)
648              SW_gl_m(i,j)       = -(dyy-1.)                          ! position gl proche de i,j <0
649              Hf                 = max(0.,H_float(i-1,j-1))          ! noeud mineur i,j
650              delHdx_SW(i,j)     = ( H(i-1,j-1) - Hf )/dyy           ! taux d'amincissement de (i-1,j-1)
651           else
652              SW_gl_m(i,j)       = -10.*dx
653              delHdx_SW(i,j)     =   0.
654           end if
655
656           if (mk_gr(i+1,j+1).eq.1)  then                            ! noeud NE grounded majeur (i+1,j+1)
657              SW_gl_m(i+1,j+1)   = dyy-1.                             ! position gl proche de i,j
658              Hf                 = max(0.,H_float(i+1,j+1))          ! noeud mineur (i+1,j+1)
659              delHdx_SW(i+1,j+1) = ( H(i+1,j+1) - Hf )/dyy           ! taux d'amincissement de (i+1,j+1)
660           else
661              SW_gl_m(i+1,j+1)   = -10.*dx
662              delHdx_SW(i+1,j+1) =   0.
663           end if
664
665
666           if (mk_gr(i+1,j-1).eq.1)  then                            ! noeud SE grounded majeur (i+1,j-1)
667              SE_gl_m(i,j)       = -(dyy-1.)                         ! position gl proche de i,j <0
668              Hf                 = max(0.,H_float(i+1,j-1))          ! noeud mineur i,j
669              delHdx_SE(i,j)     = ( H(i+1,j-1) - Hf )/dyy           ! taux d'amincissement de (i+1,j-1)
670
671           else
672              SE_gl_m(i,j)       = -10.*dx
673              delHdx_SE(i,j)     =   0.
674           end if
675
676           if (mk_gr(i-1,j+1).eq.1)  then                            ! noeud NW grounded majeur (i-1,j+1)
677              SE_gl_m(i-1,j+1)   = dyy-1.                            ! position gl proche de i,j
678              Hf                 = max(0.,H_float(i-1,j+1))          ! noeud mineur (i-1,j+1)
679              delHdx_SE(i-1,j+1) = ( H(i-1,j+1) - Hf )/dyy           ! taux d'amincissement de (i-1,j+1)
680           else
681              SE_gl_m(i-1,j+1)   = -10.*dx
682              delHdx_SE(i-1,j+1) =   0.
683           end if
684
685           
686        end if new_float
687     end do
688  end do
689
690end subroutine update_flottants
691
692
693
694!-------------------------------------------------------------------------
695!<  update_a_traiter :  points qui deviennent a traiter
696!-------------------------------------------------------------------------
697
698subroutine update_a_traiter
699
700  use declar_proto_recul
701  implicit none
702
703  integer                           :: som_voisins  !< pour des tests
704
705  ! les nouveaux points a traiter
706  do j=2,ny-1
707     do i =2,nx-1
708
709!       formulation en croix
710!        som_voisins = mk_gr(i-1,j)+mk_gr(i+1,j)+mk_gr(i,j-1)+mk_gr(i,j+1)     
711
712!        if ((mk_traiter(i,j).eq.2).and.(som_voisins.gt.0).and.(som_voisins.lt.4)) then
713!           mk_traiter(i,j) = 1
714!        end if
715
716!       Formulation avec diagonales
717       if ((mk_traiter(i,j).eq.2).and.(any(mk_gr(i-1:i+1,j-1:j+1).eq.0))) then
718          mk_traiter(i,j) = 1
719       end if
720
721     end do
722  end do
723end subroutine update_a_traiter
724
725!-------------------------------------------------------------------------
726!<  prescribe_Hp :  renvoie Hp et I_Hp vers GRISLI
727!-------------------------------------------------------------------------
728
729subroutine prescribe_Hp
730
731  use declar_proto_recul
732  implicit none
733
734! remet les masques Hp et i_Hp aux valeurs fixes
735
736      i_HP(:,:) = i_Hp0(:,:)
737      Hp(:,:)   = Hp0(:,:)
738
739i=208
740j=176
741
742  where (mk_gr(:,:).eq.0)       ! points flottants
743     i_Hp(:,:) = 1
744     Hp (:,:)  = H(:,:)         ! les points flottants ne sont pas modifies
745  end where
746
747  where(mk_traiter(:,:).eq.1)   ! points imposes Epaisseur en cours de decroissance
748     i_Hp(:,:) = 1
749     Hp (:,:)  = H(:,:)       
750  end where
751
752end subroutine prescribe_Hp
753
754!-------------------------------------------------------------------------
755!<  calc_eps_max :  calcule le epsilon max (sanity check)
756!-------------------------------------------------------------------------
757
758subroutine calc_eps_max 
759
760  use declar_proto_recul
761  implicit none
762  real              :: gamma   !< coefficient de flottaison
763 
764  gamma = ro*g*(1.+coef_Bflot)
765
766  epsmax(:,:) = Abar(:,:)*H(:,:)**3*(gamma/4.)**3
767
768end subroutine calc_eps_max
769
770!-------------------------------------------------------------------------
771!<  calc_xgl :  calcule la position sous maille de la grounding line
772!<              en supposant que l'epaisseur varie lineairement
773!-------------------------------------------------------------------------
774
775subroutine calc_xgl(dy,alpha,H_0,B_0,H_1,B_1,xpos,delHdx) 
776
777  implicit none
778
779! dummy
780  real,intent(in)        ::    dy        !< longueur de la maille
781  real,intent(in)        ::    alpha     !< coefficient de flottaison = coef_Bflot
782  real,intent(in)        ::    H_0       !< epaisseur au point pose
783  real,intent(in)        ::    B_0       !< altitude socle au point pose
784  real,intent(in)        ::    H_1       !< epaisseur au point flottant
785  real,intent(in)        ::    B_1       !< altitude socle au point flottant
786  real,intent(out)       ::    xpos      !< position de la ligne (en distance depuis le point pose 
787  real,intent(out)       ::    delHdx    !< variation d'epaisseur au noeud pose en fonction d'une variation de xpos     
788 
789
790  real                   ::    Cgl       ! variable de travail 
791
792  Cgl   = (alpha * (H_1-H_0) - (B_1-B_0)) 
793
794  if (abs(Cgl).gt.1.e-5) then   
795     xpos    = dy * (B_0 - alpha * H_0) / Cgl 
796  else
797     xpos    = dy-1.   ! verifier
798  end if
799
800! la variation d'epaisseur est proportionnelle au recul
801
802  if (xpos.gt.1.) then
803     delHdx = (H_0 - B_0 /alpha) / xpos 
804  else
805     delHdx = (H_0 - B_0 /alpha)
806  end if
807
808end subroutine calc_xgl
809
810!----------------------------------------------------------------------
811!  lect_retreat_files : lit les fichiers specifiques au retrait
812!----------------------------------------------------------------------
813subroutine  lect_retreat_files
814
815use netcdf
816use io_netcdf_grisli
817use declar_proto_recul
818
819implicit none
820
821integer    :: num               ! pour la lecture
822integer    :: num_job           ! numero du job
823real*8, dimension(:,:),   pointer  :: tab              !< tableau 2d real ecrit dans le fichier
824
825character(len=2) :: char_num    ! pour la lecture
826
827! lit le numero du job
828
829char_num = runname(7:8)
830read(char_num,*) num_job
831write(6,*) num_job
832
833! noms des fichiers d'entree
834
835number_data_B_file = trim(dirnameinp)//trim(number_data_B_file)
836region_file        = trim(dirnameinp)//trim(region_file)
837file_exp           = trim(dirnameinp)//trim(file_exp)
838
839
840! lit le fichier sigma_socle
841call Read_Ncdf_var('z',number_data_B_file,tab)      ! lit la variable 'z'
842Bsoc_sigma(:,:)  = tab(:,:)
843
844! lit le fichier region
845call Read_Ncdf_var('z',region_file,tab)             ! lit la variable 'z'
846map_region(:,:)  = nint(tab(:,:))
847
848
849! lit le fichier d'experience pour ne garder que la ligne concernee
850
851open(99,file=file_exp)
852read(99,*)                                               ! lit la ligne de titre
853
854do i=1,200
855   read(99,*,end=200) num, beta_lim_ret,retreat_1,retreat_2,(time_region(j),j=1,nb_regions)
856   if (num.eq.num_job) exit
857end do
858
859write(6,*) 'runname ',runname,'   num',num
860write(6,*) 'beta_lim,retrat1,retreat2', beta_lim_ret,retreat_1,retreat_2
861write(6,*) 'time',time_region(:)
862write(6,*) '----------------------------------------------------------'
863
864return
865200 write(6,*) " ce numero de job n est pas dans le fichier"
866
867
868
869end subroutine lect_retreat_files
870
871
872
873!----------------------------------------------------------------------
874!  init_retreat_rates : calcule les taux de retrait en fonction de divers criteres
875!----------------------------------------------------------------------
876subroutine  init_retreat0
877
878use declar_proto_recul
879implicit none
880
881! critere sur le beta
882!------------------------
883
884travail_centre(:,:) =  beta_centre(:,:)
885
886call moyennes_demi_mailles 
887
888where (travail_mx(:,:).gt.beta_lim_ret)       !---- en x
889   retreat0_x(:,:) = retreat_2*1000.
890elsewhere
891   retreat0_x(:,:) = retreat_1*1000.
892end where
893
894where (travail_my(:,:).gt.beta_lim_ret)       !---- en y
895   retreat0_y(:,:) = retreat_2*1000.
896elsewhere
897   retreat0_y(:,:) = retreat_1*1000.
898end where
899
900where (travail_SW(:,:).gt.beta_lim_ret)      !---- en SW
901   retreat0_SW(:,:) = retreat_2*1000.
902elsewhere
903   retreat0_SW(:,:) = retreat_1*1000.
904end where
905
906where (travail_SE(:,:).gt.beta_lim_ret)      !---- en SE
907   retreat0_SE(:,:) = retreat_2*1000.
908elsewhere
909   retreat0_SE(:,:) = retreat_1*1000.
910end where
911
912! enleve les regions qui ne sont pas a traiter
913
914  do j=2,ny-1
915     do i =2,nx-1
916        if (mk_traiter(i,j)*mk_traiter(i-1,j).eq.0) then     !---- en x
917           retreat0_x(i,j)  = 0.
918        end if
919        if (mk_traiter(i,j)*mk_traiter(i,j-1).eq.0) then     !---- en y
920           retreat0_y(i,j)  = 0.
921        end if
922        if (mk_traiter(i,j)*mk_traiter(i-1,j-1).eq.0) then   !---- en SW
923           retreat0_SW(i,j) = 0.
924        end if
925        if (mk_traiter(i,j)*mk_traiter(i+1,j-1).eq.0) then   !---- en SE
926           retreat0_SE(i,j) = 0.
927        end if
928     end do
929  end do
930
931
932end subroutine init_retreat0
933
934
935!----------------------------------------------------------------------
936!  moyennes_demi_mailles : calcule les moyennes d'un tableau ! travail
937
938!----------------------------------------------------------------------
939subroutine moyennes_demi_mailles 
940use declar_proto_recul
941implicit none
942  do j=2,ny-1
943     do i =2,nx-1
944        travail_mx(i,j) = 0.5* ( travail_centre(i,j)  + travail_centre(i-1,j) )
945        travail_my(i,j) = 0.5* ( travail_centre(i,j)  + travail_centre(i,j-1) )
946        travail_SW(i,j) = 0.5* ( travail_centre(i,j)  +travail_centre(i-1,j-1))
947        travail_SE(i,j) = 0.5* ( travail_centre(i,j)  +travail_centre(i+1,j-1))
948     end do
949  end do
950
951end subroutine moyennes_demi_mailles
952
953!----------------------------------------------------------------------
954!  init_time_depart : donne les temps de depart en fonction des regions
955!----------------------------------------------------------------------
956
957subroutine init_time_depart
958
959
960use declar_proto_recul
961implicit none
962real               :: time_inacc      = 10000.   ! temps donne aux noeuds inaccessibles
963real               :: t_already_float = -10.     ! temps donne aux noeuds qui flottent deja
964real               :: time_run_begin  = 2000.    ! debut du run             
965
966
967! region par region
968
969  time_region(0) =  t_already_float-time_run_begin
970
971  do j=2,ny-1
972     do i=2,nx-1
973        time_dep(i,j) = (time_region(map_region(i,j)) - time_run_begin)
974     end do
975  end do
976
977  where (mk_traiter(:,:).eq.0)
978     time_dep(:,:) = time_inacc
979  end where
980
981! calcul sur les demi-mailles
982  travail_centre(:,:) = time_dep(:,:)
983
984  call moyennes_demi_mailles
985
986  time_dep_mx(:,:) = travail_mx(:,:)
987  time_dep_my(:,:) = travail_my(:,:)
988  time_dep_SW(:,:) = travail_SW(:,:)
989  time_dep_SE(:,:) = travail_SE(:,:)
990
991
992! le temps est tres eleve pour les zones inaccessibles
993  do j=2,ny-1
994     do i =2,nx-1
995        if (mk_traiter(i,j)*mk_traiter(i-1,j).eq.0) then     !---- en x
996           time_dep_mx(i,j)  = time_inacc
997        end if
998        if (mk_traiter(i,j)*mk_traiter(i,j-1).eq.0) then     !---- en y
999           time_dep_my(i,j)  = time_inacc
1000        end if
1001        if (mk_traiter(i,j)*mk_traiter(i-1,j-1).eq.0) then   !---- en SW
1002           time_dep_SW(i,j) = time_inacc
1003        end if
1004        if (mk_traiter(i,j)*mk_traiter(i+1,j-1).eq.0) then   !---- en SE
1005           time_dep_SE(i,j) = time_inacc
1006        end if
1007     end do
1008  end do
1009
1010! pour bien voir les régions dans time_float
1011  where(mk_traiter(:,:).eq.2)
1012     time_float(:,:) = tend + 100.
1013  elsewhere (mk_traiter(:,:).eq.0)
1014     time_float(:,:) = time_inacc
1015  elsewhere (mk_traiter(:,:).lt.0)
1016     time_float(:,:) = t_already_float
1017  end where
1018
1019
1020end subroutine init_time_depart
1021
1022!----------------------------------------------------------------------
1023!<  bedrock_sigma : calcul la tolerance sur le socle
1024!                   en fonction de bed_sigma_1 et bed_sigma_2!                 
1025!----------------------------------------------------------------------
1026
1027subroutine  bedrock_sigma
1028
1029
1030use declar_proto_recul
1031implicit none
1032
1033real        :: prec_min = 20.          ! precision min de la mesure
1034
1035
1036where (Bsoc_sigma (:,:).lt.1.)         ! pas de donnee socle
1037   travail_centre(:,:) = bed_sigma_2
1038elsewhere (Bsoc_sigma (:,:).ge.1.)     !
1039   travail_centre(:,:) = max(prec_min,bed_sigma_1/sqrt(Bsoc_sigma(:,:)))
1040end where
1041Bsoc_sigma (:,:) = travail_centre(:,:)
1042
1043! zones particulieres
1044
1045where(mk_traiter(:,:).eq.0) 
1046   Bsoc_sigma(:,:) = -10000            ! zones interdites
1047elsewhere (mk_traiter(:,:).lt.0)
1048   Bsoc_sigma(:,:) = 0.
1049end where
1050   
1051end subroutine bedrock_sigma
1052
1053
1054!----------------------------------------------------------------------
1055!< teste_socle_voisins : pour un point i,j : teste si les voisins sont plus haut                           
1056!----------------------------------------------------------------------
1057
1058subroutine teste_socle_voisins(ii,jj,A_x,A_y,A_SW,A_SE)
1059
1060use declar_proto_recul
1061implicit none
1062
1063integer,intent(in)                 :: ii,jj         ! indice du point considere
1064real,dimension(2),intent(inout)    :: A_x           ! quand on balaye les voisins
1065real,dimension(2),intent(inout)    :: A_y           ! quand on balaye les voisins
1066real,dimension(2),intent(inout)    :: A_SW          ! quand on balaye les voisins
1067real,dimension(2),intent(inout)    :: A_SE          ! quand on balaye les voisins
1068
1069
1070
1071! Btest est plus haut = Bsoc + Bsoc_sigma
1072
1073
1074! Mais pas suffisamment
1075
1076if (Btest(ii-1,jj  ).le.Bsoc(ii,jj))  A_x(1) = 0.  !---- en x
1077if (Btest(ii+1,jj  ).le.Bsoc(ii,jj))  A_x(2) = 0.
1078if (Btest(ii  ,jj-1).le.Bsoc(ii,jj))  A_y(1) = 0.  !---- en y
1079if (Btest(ii  ,jj+1).le.Bsoc(ii,jj))  A_y(2) = 0.
1080
1081if (Btest(ii-1,jj-1).le.Bsoc(ii,jj))  A_SW(1) = 0. !---- en SW
1082if (Btest(ii+1,jj+1).le.Bsoc(ii,jj))  A_SW(2) = 0.
1083if (Btest(ii+1,jj-1).le.Bsoc(ii,jj))  A_SE(1) = 0. !---- en SE
1084if (Btest(ii-1,jj+1).le.Bsoc(ii,jj))  A_SE(2) = 0.
1085
1086end subroutine teste_socle_voisins
1087
1088
1089!----------------------------------------------------------------------
1090!< sanity_check : pour un point i,j  limite le dhdt                         
1091!----------------------------------------------------------------------
1092
1093subroutine sanity_check(ii,jj)
1094
1095use declar_proto_recul
1096implicit none
1097
1098integer,intent(in)                 :: ii,jj         ! indice du point considere
1099
1100real                               :: dH_x          ! variable de travail
1101real                               :: dH_y          ! variable de travail
1102
1103DelHdt_sanity(ii,jj) = 0.
1104
1105! en x -----------------------------------------------------------------------------------
1106if ((Uxbar(ii,jj).ge.0.).and.(Uxbar(ii+1,jj).ge.0.)) then             ! ecoulement --->
1107
1108   dH_x   =  Uxbar(ii,jj)*(H(ii,jj)-H(ii-1,jj)) / dx
1109   dH_x   =  Dh_x +  epsmax(ii,jj) * H(ii,jj)
1110
1111else if ((Uxbar(ii,jj).le.0.).and.(Uxbar(ii+1,jj).le.0.)) then        ! ecoulement <---
1112
1113   dH_x  =  Uxbar(ii+1,jj)*(H(ii+1,jj)-H(ii,jj)) / dx
1114   dH_x  =  Dh_x +  epsmax(ii,jj) * H(ii,jj)
1115
1116else if ((Uxbar(ii,jj).le.0.).and.(Uxbar(ii+1,jj).ge.0.)) then        ! ecoulement  <------>
1117
1118   dH_x   =  Dh_x +  epsmax(ii,jj) * H(ii,jj)
1119else                                                                  ! ecoulement ---><-----
1120   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
1121!  mais ici pas de terme en epsmax car confine
1122endif
1123
1124
1125! en y -----------------------------------------------------------------------------------
1126  if ((Uybar(ii,jj).ge.0.).and.(Uybar(ii,jj+1).ge.0.)) then           ! ecoulement  ^
1127                                                       
1128   dH_y   =  Uybar(ii,jj)*(H(ii,jj)-H(ii,jj-1)) / dx
1129   dH_y   =  Dh_y +  epsmax(ii,jj) * H(ii,jj)
1130
1131else if ((Uybar(ii,jj).le.0.).and.(Uybar(ii,jj+1).le.0.)) then        ! ecoulement  v
1132
1133   dH_y  =  Uybar(ii,jj+1)*(H(ii,jj+1)-H(ii,jj)) / dx
1134   dH_y  =  Dh_y +  epsmax(ii,jj) * H(ii,jj)
1135
1136else if ((Uybar(ii,jj).le.0.).and.(Uybar(ii,jj+1).ge.0.)) then        ! ecoulement divergent
1137
1138   dH_y   =  Dh_y +  epsmax(ii,jj) * H(ii,jj)
1139else                                                                  ! ecoulement convergent
1140   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
1141!  mais ici pas de terme en epsmax car confine
1142endif
1143
1144DelHdt_sanity(ii,jj) = max(dH_x + dH_y,0.)
1145
1146end subroutine sanity_check
1147
1148end module proto_recul_mod
1149!----------------------------------------------------------------------------------
1150
1151
1152
1153
1154
Note: See TracBrowser for help on using the repository browser.