Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

prescribe-H_mod.f90 @ 25

Last change on this file since 25 was 4, checked in by dumas, 10 years ago
initial import GRISLI trunk
File size: 26.5 KB

Line
1	!> \file precribe-H_mod.f90
2	!!
3	!! Determine mask and thickness for prescription of H
4	!!
5	!<
6
7	!> \namespace prescribe_H
8	!!
9	!! Determine mask and thickness for prescription of H
10	!!
11	!! \author CatRitz
12	!! \date june 2009
13	!!
14	!! @note use module3D_phy
15	!! @note use param_phy_mod
16	!<
17
18	module prescribe_H
19
20	use module3D_phy
21	use param_phy_mod
22	use interface_input
23	implicit none
24
25	! ces dimensionnements sont maintenant dans 3D
26
27	! real,dimension(nx,ny) :: Hp !< H value if prescribed
28	! real,dimension(nx,ny) :: Hp0 !< H value if prescribed (reference value)
29	! real,dimension(nx,ny) :: Delta_H !< Delta_H value if prescribed
30	! integer,dimension(nx,ny) :: i_delta_H !< 1 if Delta_H is prescribed on this node, else 0
31	! integer,dimension(nx,ny) :: i_Hp !< 1 if H is prescribed on this node, else 0
32	! integer,dimension(nx,ny) :: i_Hp0 !< i_hp mask reference value does not change with time
33	! integer,dimension(nx,ny) :: MK_gl0 !< mask grounding line initial
34	! integer,dimension(nx,ny) :: MK_flot0 !< mask float initial
35
36	! pour grounding line retrea, ice2sea
37	integer,dimension(nx,ny) :: i_status !< 0 pas prescrit, 1 prescrit retrait, 2 prescrit shelf
38	real :: time_onset ! moment ou la grounding line commence a reculer
39
40	! for prescribed grounding line retreat : ancienne version
41	integer :: voisin !< work variable
42	real (kind=kind(0.d0)) :: t_break !< temps : en double precision
43
44	integer,parameter :: iretreat_max = 20 !< maximum number of retreat steps
45	integer :: iretreat
46	real (kind=kind(0.d0)),dimension(iretreat_max) :: time_gl !< time of the retreat steps
47	integer, dimension(nx,ny) :: tag_gl !< tag of all the retreat points
48
49	!< floating -> 10000
50	!< grounded not gl -> 1000
51
52	integer :: floating_tag = -10000
53	integer :: grounded_tag = 1000
54
55
56	integer, dimension(nx,ny) :: upwind_gl !< 1 if gl can still retreat, else 0
57	real, dimension(nx,ny) :: H_gl !< flotation thickness at the time
58	real, dimension(nx,ny) :: dH_gl_dt !< to smoothly decrease ice thickness
59	!< for next gl points
60
61	! nouvelle version recul grounding line pour ice2sea
62	! ------------------------------------------------------
63
64	! real, dimension(nx,ny) :: H_gl !< flotation thickness at the time (deja declare autre version)
65	real, dimension(nx,ny) :: H_beggin !< thickness when the point begins to be prescribed
66	real, dimension(nx,ny) :: time_beggin !< time when the point begins to be prescribed
67	real, dimension(nx,ny) :: time_just_gr !< time when the point is just grounded
68
69	contains
70
71	!______________________________________________________________________________________
72	!> initialise prescribe H
73	!! calculate the initial grounding line position
74	!! @return MK_flot0 and Mk_gl0
75
76	subroutine init_prescribe_H
77
78	implicit none
79	if (itracebug.eq.1) call tracebug(' Entree dans routine init_prescribe_H')
80
81	! keep the mask of floating points
82	where (flot(:,:))
83	MK_flot0(:,:)= 1
84	elsewhere
85	MK_flot0(:,:)= 0
86	end where
87
88
89	! determine the initial grounding line
90
91	MK_gl0(:,:)=0
92	HP0(:,:)=H0(:,:)
93
94
95	do j=2,ny-1
96	do i=2,nx-1
97	voisin = MK_flot0(i-1,j) + MK_flot0(i+1,j) + MK_flot0(i,j-1) + MK_flot0(i,j+1)
98	if ((MK_flot0(i,j).eq.0).and.(voisin.gt.0)) then ! grounded and at least one neighbour floating
99	MK_gl0(i,j) = 1 ! tagged mask grounding line
100
101	end if
102	end do
103	end do
104
105
106	call prescribe_grid_boundary_0
107
108	!call init_prescribe_retreat
109	time_gl(1)= 1000000 ! pour un run sans recul
110	t_break = 1000000. ! pour un run avec ice shelves prescrits
111
112
113	! Dans le cas du Groenland ice2sea, mk_init = 3 -> noeuds qui ne bougent pas
114
115	if (geoplace(1:4).eq.'GI2S') then
116
117	do j = 1,ny
118	do i = 1,nx
119	if (mk_init(i,j).eq.3) then
120	i_hp0 (i,j) = 1
121	Hp0 (i,j) = H0(i,j)
122	endif
123	end do
124	end do
125	end if
126
127
128
129	if (igrdline.eq.3) then
130	call init_prescr_ice2sea_retreat
131	end if
132
133
134	end subroutine init_prescribe_H
135	!______________________________________________________________________________________
136	!> function prescribe_present_H_gl
137	!! calculate the initial grounding line position
138	!! @return i_hp(:,:) and hp(:,:)
139
140	subroutine prescribe_present_H_gl
141
142	implicit none
143
144	if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_present_H_gl')
145
146	where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
147	i_hp(:,:) = 1 ! thickness prescribed to present value
148	hp(:,:) = Hp0(:,:)
149	end where
150	if (itracebug.eq.1) call tracebug(' fin prescribe_present_H_gl ')
151
152
153	end subroutine prescribe_present_H_gl
154
155	!______________________________________________________________________________________
156	!> function prescribe_fixed_points
157	!! fixes thickness for some points
158	!! @return i_hp(:,:) and hp(:,:)
159
160	subroutine prescribe_fixed_points
161
162	if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_fixed_points')
163	where (i_hp0(:,:).eq.1) ! les points i_hp0 le sont pour tout le run
164	i_hp(:,:) = i_hp0(:,:)
165	Hp(:,:) = Hp0(:,:)
166	end where
167
168	if (itracebug.eq.1) call tracebug(' fin prescribe_fixed_points')
169
170
171	end subroutine prescribe_fixed_points
172
173	!______________________________________________________________________________________
174	!> function prescribe_paleo_gl_shelf
175	!! calculate the grounding line thickness
176	!! no ice shelves
177	!! @return i_hp(:,:) and hp(:,:)
178
179	subroutine prescribe_paleo_gl_shelf
180
181	implicit none
182
183	if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_paleo_gl_shelf')
184
185
186	! noeuds qui doivent être imposés
187
188	where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
189	i_hp(:,:) = 1 ! thickness prescribed
190	end where
191
192	! valeur imposee
193	where (MK_flot0(:,:).eq. 1) ! paleo shelf epaisseur a 1
194	hp(:,:) = 1.
195	end where
196
197
198	where (MK_gl0(:,:) .eq. 1)
199	hp(:,:) = (sealevel-Bsoc(:,:))*row/ro +1.
200	end where
201
202
203	end subroutine prescribe_paleo_gl_shelf
204
205
206	!______________________________________________________________________________________
207	!> function prescribe_grid_boundary
208	!! prescribe the grid boundary thickness to 0
209	!! @return i_hp(:,:) and hp(:,:)
210
211	subroutine prescribe_grid_boundary_0
212
213	implicit none
214	if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_grid_boundary_0')
215
216	! prescribe thickness on the grid boundaries
217	i_hp(:,:) = 0
218
219	i_hp(1,:) = 1
220	i_hp(nx,:)= 1
221	hp(1,:) = 0
222	hp(nx,:) = 0
223
224	i_hp(:,1) = 1
225	i_hp(:,ny)= 1
226	hp(:,1) = 0
227	hp(:,ny) = 0
228
229	! valeurs de reference
230	i_hp0(1,:) = 1
231	i_hp0(nx,:) = 1
232	hp0(1,:) = 0
233	hp0(nx,:) = 0
234
235	i_hp0(:,1) = 1
236	i_hp0(:,ny) = 1
237	hp0(:,1) = 0
238	hp0(:,ny) = 0
239
240
241
242
243	end subroutine prescribe_grid_boundary_0
244
245	!______________________________________________________________________________________
246	!______________________________________________________________________________________
247	!> lecture of a prescribe field of delta H
248	!!
249	!! @return Delta_H, i_delta_H
250
251	subroutine lect_prescribed_deltaH
252
253	implicit none
254	character(len=100) :: delta_H_file ! variation d'epaisseur imposee
255
256	if (itracebug.eq.1) call tracebug(' Entree dans routine lect_prescribed_deltaH')
257	if (itracebug.eq.1) call tracebug(' Attention ne dvrait pas passer par la')
258	! lecture
259	delta_H_file = 'sensibxg-25km.dat'
260	delta_H_file = trim(dirnameinp)//trim(delta_H_file)
261
262
263
264
265	call lect_input(3,'Delta_H',1,Delta_H,delta_H_file,trim(dirnameinp)//trim(runname)//'.nc')
266	!call lect_datfile(nx,ny,Delta_H,1,delta_H_file)
267
268	do j=1,ny
269	do i=1,nx
270	if (Delta_H(i,j).ne.0.) then
271	i_delta_H(i,j)=1
272	i_Hp(i,j)=1 ! i_delta_H might be not completely included in i_Hp
273	Hp0(i,j) = H(i,j)-Delta_H(i,j) ! warning Delta_H is substracted
274
275	if (MK_flot0(i,j).eq. 1) then ! si flotte, pas de changement d'epaisseur
276	Hp0(i,j)=H0(i,j)
277	else if (Hp0(i,j).lt. -Bsoc(i,j)*row/ro) then ! devient flottant
278	Hp0(i,j) = -Bsoc(i,j)*row/ro-10. ! floating thickness - 10
279	MK_flot0(i,j)= 1 ! est tagge masque flottant
280	end if
281	Hp(i,j) = Hp0(i,j)
282	else
283	i_delta_H(i,j)=0
284	end if
285	end do
286	end do
287	end subroutine lect_prescribed_deltaH
288
289	!______________________________________________________________________________________
290	!> substract Delta_H to the current thickness
291	!!
292	!! @return
293	subroutine prescribe_deltaH
294
295	implicit none
296
297	if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_deltaH')
298	where (i_delta_H(:,:).eq.1)
299	Hp(:,:) = Hp0(:,:)
300	i_Hp(:,:)=1
301	end where
302	end subroutine prescribe_deltaH
303	!______________________________________________________________________________________
304	!> break ice shelves
305	!!
306	!!
307	subroutine break_all_ice_shelves
308
309	implicit none
310	if (itracebug.eq.1) call tracebug(' Entree dans routine break_all_ice_shelves ')
311
312
313	debug_3D(:,:,89)=0.
314
315	where (flot(:,:)) ! floating from the beginning and eventually from grounding line retreat
316	i_hp(:,:) = 1 ! thickness prescribed to 1 m
317	hp (:,:) = 1.
318	hp0 (:,:) = 1.
319	H(:,:) =1.
320	debug_3D(:,:,89)=1.
321
322	end where
323	where (MK_flot0(:,:).eq. 1)
324	i_hp(:,:) = 1 ! thickness prescribed to 1 m
325	hp (:,:) = 1.
326	hp0 (:,:) = 1.
327	H(:,:) =1.
328	end where
329
330
331	end subroutine break_all_ice_shelves
332
333	!
334	!______________________________________________________________________________________
335	!> break ice shelves
336	!!
337	!!
338	subroutine melt_ice_shelves
339
340	implicit none
341	integer :: nbvoisins,nbdeglac,iv,jv
342	real :: hmoy,hmin,hmax
343
344	if (itracebug.eq.1) call tracebug(' Entree dans routine melt_ice_shelves')
345
346	do j = 2,ny-1
347	do i=2, nx-1
348	if (flot(i,j)) then ! floating
349	i_hp(i,j) = 1 ! thickness prescribed
350	hmoy=0.
351	hmin=0.
352	hmax=0.
353	nbvoisins=0
354	nbdeglac=0
355	do jv=-1,1
356	do iv=1,1
357	if (flot(i+iv,j+jv)) then
358	nbvoisins=nbvoisins+1
359	hmoy=hmoy+h(i+iv,j+jv)
360	hmin=min(hmin,h(i+iv,j+jv))
361	hmax=max(hmax,h(i+iv,j+jv))
362	if (h(i+jv,j+jv).lt.10.) nbdeglac=nbdeglac+1
363	end if
364	end do
365	end do
366	hmoy=hmoy/(max(nbvoisins,1))
367
368	if(H(i,j).gt.hmoy) then ! melt only if smooth
369	hp(i,j)=H(i,j)-1.*dt
370	else if (hmax-hmin.lt.10.) then
371	hp(i,j)=H(i,j)-1.*dt
372	else if ((H(i,j).lt.50.).and.(Hmin.lt.30.)) then
373	hp(i,j)=1.
374	else if ((H(i,j).lt.50.).and.(nbdeglac.ge.4)) then
375	hp(i,j)=1.
376	endif
377	hp(i,j) = max(hp(i,j),1.)
378	end if
379	end do
380	end do
381
382	end subroutine melt_ice_shelves
383
384
385
386	!
387	!______________________________________________________________________________________
388	!> prescribe grounding line retreat
389	!!
390	!! time_gl(iretreat_max) : times of the retreat steps (in)
391	!! tag_gl(nx,ny) : tag that gives at which time a node is prescribed
392	!! as a grounding line node
393	!! @return Hp
394	!! @return i_hp
395
396	subroutine prescribe_retreat
397
398	implicit none
399	integer, dimension(nx,ny) :: new_tag ! working array
400
401	if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_retreat')
402	if (itracebug.eq.1) write(num_tracebug,*) 'time_gl(1)=', time_gl(1)
403
404	new_tag(:,:) = tag_gl(:,:) ! initialize new_tag
405
406	start_time: if (time.lt.time_gl(1)) then ! before the perturbation
407
408	if (itracebug.eq.1) call tracebug(' Avant appel prescribe_present_H_gl')
409
410	else
411	if (itracebug.eq.1) call tracebug('avant entree step_gl')
412
413	iretreat_loop: do iretreat = 1,iretreat_max-1
414
415	if (itracebug.eq.1) write(num_tracebug,*) 'iretreat',iretreat
416
417	find_step: if (abs(time-time_gl(iretreat)).lt.dtmin/2.) then
418	! exactly on a retreat step
419	! write(6,*) time, 'sur un step',iretreat
420
421	where (tag_gl(:,:).eq.iretreat)
422	H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness
423	Hp(:,:) = H_gl(:,:)+5. ! +5 to be sure to be grounded
424	! even with small isostatic changes
425	i_Hp(:,:) = 1
426
427	elsewhere (tag_gl(:,:).eq.iretreat+1) ! points that will be gl at next step
428
429	H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness
430	dH_gl_dt(:,:) = (H(:,:)-H_gl(:,:))-5. ! thickness to be reduced to reach gl
431	dH_gl_dt(:,:) = dH_gl_dt(:,:) / (time_gl(iretreat+1)-time_gl(iretreat))
432	! rate of decrease
433
434	elsewhere ((tag_gl(:,:).eq.iretreat-1).and. & ! points that become floating
435	(upwind_gl(:,:).eq.1))
436
437	H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness
438	Hp(:,:) = H_gl(:,:)-10. ! to be sure to be floating
439	i_Hp(:,:) = 1
440	new_tag(:,:) = floating_tag
441
442	end where
443
444	! pour une raison mysterieuse la boucle suivante marche et celle en where ne marche pas
445	! (memory falut)
446	do j=1,ny
447	do i=1,nx
448	if ((tag_gl(i,j).eq.iretreat-1).and. & ! points that stay gl
449	(upwind_gl(i,j).eq.0)) then
450	H_gl(i,j) = -Bsoc(i,j)*row/ro ! flotation thickness
451	Hp(i,j) = H_gl(i,j)+5. ! to be sure to be grounded
452	i_Hp(i,j) = 1
453	new_tag(i,j) = iretreat ! tagged with the new step number
454	! write(6,*) i,j,tag_gl(i,j), upwind_gl(i,j)
455	endif
456	end do
457	end do
458
459	!!$ where ((tag_gl(:,:).eq.iretreat-1).and. & ! points that stay gl
460	!!$ (upwind_gl(:,:).eq.0))
461	!!$
462	!!$ H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness
463	!!$ Hp(:,:) = H_gl(:,:)+5. ! to be sure to be grounded
464	!!$ i_Hp(:,:) = 1
465	!!$ new_tag(:,:) = iretreat ! tagged with the new step number
466	!!$ end where
467
468	else if ((time.gt.time_gl(iretreat)+dtmin/2.).and.(time.lt.time_gl(iretreat+1))) then
469	! write(6,*) time, 'entre',time_gl(iretreat),time_gl(iretreat+1)
470
471	where (tag_gl(:,:).eq.iretreat)
472	H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness
473	Hp(:,:) = H_gl(:,:)+5. ! +5 to be sure to be grounded
474	! even with small isostatic changes
475	i_Hp(:,:) = 1
476
477	elsewhere (tag_gl(:,:).eq.iretreat+1) ! points that will be gl at next step
478	Hp(:,:) = H(:,:)-dH_gl_dt(:,:)*(time-time_gl(iretreat))
479	i_Hp(:,:) = 1
480
481	elsewhere (tag_gl(:,:).eq.floating_tag)
482	Hp(:,:) = H(:,:)
483	i_Hp(:,:) = 1
484	end where
485
486	else if (time.gt.time_gl(iretreat_max)) then
487	! write(6,*) time, 'apres iretreat_max',iretreat,time_gl(iretreat_max)
488
489	where (tag_gl(:,:).eq.iretreat_max)
490	H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness
491	Hp(:,:) = H_gl(:,:)+5. ! +5 to be sure to be grounded
492	! even with small isostatic changes
493	i_Hp(:,:) = 1
494
495	elsewhere (tag_gl(:,:).eq.floating_tag)
496	Hp(:,:) = H(:,:)
497	i_Hp(:,:) = 1
498	end where
499
500	!faire une sortie de la boucle step_gl
501
502	end if find_step
503	end do iretreat_loop
504
505	end if start_time
506
507	tag_gl(:,:) = new_tag(:,:)
508	debug_3D(:,:,84) = tag_gl(:,:)
509	debug_3D(:,:,85) = i_Hp(:,:)
510	debug_3D(:,:,86) = Hp(:,:)
511
512	end subroutine prescribe_retreat
513	!
514	!______________________________________________________________________________________
515	!> initialise prescribe retreat
516	!!
517	!! @return tag_gl
518	!! @return upwind_gl
519	!! @return time_gl
520
521	subroutine init_prescribe_retreat
522
523	implicit none
524	real, dimension(nx,ny) :: lect_tab
525	if (itracebug.eq.1) call tracebug(' Entree dans routine init_prescribe_retreat')
526
527	! read the prescribed gl file
528
529	call lect_input(3,'lecttab',1,lect_tab,trim(dirnameinp)//'masque_gl_15km_final.dat',trim(dirnameinp)//trim(runname)//'.nc')
530	!call lect_datfile(nx,ny,lect_tab,1,trim(dirnameinp)//'masque_gl_15km_final.dat')
531
532	tag_gl(:,:)=nint(lect_tab(:,:))
533
534	! look if neighbors could become grounding line later
535	upwind_gl(:,:) = 0
536	do j=2,ny-1
537	do i=2,nx-1
538	if ((tag_gl(i,j).ne.floating_tag).and.(tag_gl(i,j).ne.grounded_tag)) then
539	if (any(tag_gl(i-1:i+1,j-1:j+1).eq.tag_gl(i,j)+1)) then
540	upwind_gl(i,j) = 1
541	end if
542	end if
543	end do
544	end do
545
546	! time of the retreat steps
547
548	do i=1,iretreat_max
549	time_gl(i)=(i-1)*30.+100000.
550	end do
551	time_gl(iretreat)=200000.
552
553	write(num_rep_42,*) '! Grounding line retreat prescribed'
554	write(num_rep_42,*) '! time_gl begin et suivant : ',time_gl(1),time_gl(2), &
555	' time_glmax=',time_gl(iretreat)
556
557	! beginning of breaking or melting ice shelves
558
559	t_break = 1000000.d00
560
561	write(num_rep_42,*) '! t_break=', t_break
562
563	end subroutine init_prescribe_retreat
564	!______________________________________________________________________________________
565	!______________________________________________________________________________________
566	!> function prescribe_present_H_gl
567	!! calculate the initial grounding line position
568	!! @return i_hp(:,:) and hp(:,:)
569
570	subroutine prescribe_present_H_gl_copy
571
572	implicit none
573	if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_present_H_gl-copy')
574
575	do j=1,ny
576	do i=1,nx
577	if ((MK_flot0(i,j).eq. 1).or.(MK_gl0(i,j) .eq. 1)) then
578	i_hp(i,j) = 1 ! thickness prescribed to present value
579	hp(i,j) = Hp0(i,j)
580	end if
581	end do
582	end do
583
584	!!$where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
585	!!$ i_hp(:,:) = 1 ! thickness prescribed to present value
586	!!$! hp(:,:) = 1 !Hp0(:,:)
587	!!$end where
588
589
590	end subroutine prescribe_present_H_gl_copy
591
592	!______________________________________________________________________________________
593	!> initialise prescr_ice2sea_retreat
594	!!
595	!! @return tag_gl
596	!! @return upwind_gl
597	!! @return time_gl
598
599	! nouvelle routine de recul impose de grounding line pour ice2sea : decembre 2011
600	! cette routine transforme les informations distance en information temps ou un noeud sera just grounded"
601
602
603	subroutine init_prescr_ice2sea_retreat
604
605	implicit none
606
607	character(len=100) :: dist_file !< nom fichier contenant les distances
608	character(len=100) :: file_ncdf !< nom bidon pour lecture netcdf
609
610
611	real, dimension(nx,ny) :: dist_gl ! distance à la grounding line ( a lire)
612	real :: vretrait ! vitesse de retrait
613	real :: delta_t
614
615	real :: time_max ! pour calcul time_max
616	integer :: i1
617	integer :: j1
618
619
620	namelist/retreat_ice2sea/dist_file,vretrait,time_onset
621
622
623	if (itracebug.eq.1) call tracebug(' Entree dans routine init_prescr_ice2sea_retreat')
624
625	! lecture des noms de fichiers et parametres
626	!--------------------------------------------
627	428 format(A)
628
629	rewind(num_param) ! pour revenir au debut du fichier param_list.dat
630	read(num_param,retreat_ice2sea)
631
632	write(num_rep_42,428)'!___________________________________________________________'
633	write(num_rep_42,428)'! prescr_ice2sea_retreat '
634	write(num_rep_42,retreat_ice2sea)
635	write(num_rep_42,428)'!___________________________________________________________'
636
637
638	! lecture de la carte de distance dites
639	!-------------------------------------
640
641	! distance file : distances en km
642
643	dist_file = trim(dirnameinp)//trim(dist_file)
644	call lect_input(1,'dist_gl',1,dist_gl,dist_file,file_ncdf)
645
646	! les zones interdites sont taggees -1 -> les transformer en +100000
647	where ((dist_gl(:,:).eq.-1.).or.(Bsoc0(:,:).gt.0.))
648	dist_gl(:,:) = 100000.
649	end where
650
651
652	! calcul de l'epaisseur H_gl, correction des distances pour mieux prendre en compte le premier recul
653
654	where(Mk_flot0(:,:).eq.0) ! points poses y compris GL
655
656	h_gl(:,:) = (sealevel-Bsoc(:,:))*row/ro ! epaisseur flottaison
657	dist_gl(:,:) = dist_gl(:,:)+2.+2.5 ! 2 Pour avoir au moins 0 sur la GL
658	! 2.5 = demi maille 5 km
659
660	time_just_gr(:,:) = dist_gl(:,:)/Vretrait + time_onset ! aussi valable pour les points GL
661	H_beggin(:,:) = H0(:,:) ! pour initialiser, strictement vrai pour la GL
662
663	elsewhere ! points flottants toujours prescrit
664	time_just_gr(:,:) = time_onset-100000.
665	time_beggin(:,:) = time_onset-110000
666	i_status(:,:) = 2
667	H_beggin(:,:) = H0(:,:)
668
669	end where
670
671
672
673	! calcul specifique a la grounding line
674
675	where(Mk_gl0(:,:).eq.1)
676	time_beggin(:,:) = time_onset
677	H_beggin(:,:) = H0(:,:)
678	end where
679
680	! recherche de Time_beggin : pour les points grounded mais pas gl
681
682	delta_t = dx/1000./Vretrait -0.1 ! temps pour traverser une maille
683
684	do j=2,ny-1
685	do i=2,nx-1
686	if((Mk_flot0(i,j).eq.0).and.(Mk_gl0(i,j).eq.0)) then ! grounded mais pas GL
687	time_max = time_onset
688	do j1 = j-1,j+1
689	do i1 = i-1,i+1
690	if (Time_just_gr(i1,j1).lt.Time_just_gr(i,j)-delta_t) then ! 0.1 par securite
691	time_max = max(time_max,Time_just_gr(i1,j1))
692	endif
693	end do
694	end do
695
696	if (Time_just_gr(i,j)-time_max.lt.2.*delta_t) then
697	time_beggin(i,j) = time_max
698	else
699	time_beggin(i,j) = Time_just_gr(i,j)-2.*delta_t
700	end if
701
702
703	end if
704	end do
705	end do
706
707	! impression pour verifier ( a enlever apres)
708	!do j=1,ny
709	! do i=1,nx
710	! write(666,'(i3,1x,i3,3(f8.0,1x))') i,j,Time_just_gr(i,j),time_beggin(i,j),Time_just_gr(i,j)-time_beggin(i,j)
711	! end do
712	!end do
713
714
715	end subroutine init_prescr_ice2sea_retreat
716	!______________________________________________________________________________________
717
718	subroutine prescr_ice2sea_retreat
719
720	implicit none
721	real :: t_local ! pour les tests car time est en double precision
722
723	if (itracebug.eq.1) call tracebug(' Entree dans routine prescr_ice2sea_retreat')
724
725	!temps dans la boucle
726	t_local = time
727
728	! avant le demarrage la grounding line est fixee a sa valeur actuelle
729
730	if (t_local.lt.time_onset) then
731
732	where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
733	i_hp(:,:) = 1 ! thickness prescribed to present value
734	hp(:,:) = Hp0(:,:)
735	end where
736
737	else ! calcul du H prescrit avec recul
738
739	do j=2,ny-1
740	do i= 2,nx-1
741
742	! point pas encore impose
743
744	time_test: if (t_local.lt.time_beggin(i,j)) then
745	i_Hp(i,j) = 0
746	i_status(i,j) = 0
747
748	! point en periode de retrait
749
750	else if ( (t_local.ge.time_beggin(i,j)).and.(t_local.le.time_just_gr(i,j)) ) then
751
752	retrait: if (i_status(i,j).eq.0) then ! on vient de rentrer dans la periode retrait
753	i_status(i,j) = 1
754	H_beggin(i,j) = H(i,j)
755	i_Hp(i,j) = 0 ! on n'impose pas H le premier coup
756
757	else if (i_status(i,j) .eq. 1) then
758	i_Hp(i,j) = 1 ! on impose ce point
759	h_gl(i,j) = (sealevel-Bsoc(i,j))*row/ro ! epaisseur flottaison
760
761
762	! decroissance lineaire entre time_beggin et time_just_gr
763
764	Hp(i,j) = (H_beggin(i,j) - H_gl(i,j)) * ( (time_just_gr(i,j) - t_local) &
765	/ (time_just_gr(i,j) - time_beggin(i,j))) + H_gl(i,j)
766	end if retrait
767
768	! point en periode flottantte
769
770	else if (t_local.gt.time_just_gr(i,j)) then
771
772	h_gl(i,j) = (sealevel-Bsoc(i,j))*row/ro ! epaisseur flottaison
773	i_status(i,j) = 2
774	i_Hp(i,j) = 1 ! on impose ce point
775	Hp(i,j) = min(H_gl(i,j)-10., H0(i,j)) ! minimum flottaison ou epaisseur initiale
776
777	end if time_test
778	end do
779	end do
780	end if
781	end subroutine prescr_ice2sea_retreat
782
783
784	end module prescribe_H

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

Context Navigation

source: trunk/SOURCES/prescribe-H_mod.f90 @ 25

Download in other formats: