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 | real,dimension(nx,ny) :: Hp !< H value if prescribed |
---|
25 | real,dimension(nx,ny) :: Hp0 !< H value if prescribed (reference value) |
---|
26 | real,dimension(nx,ny) :: Delta_H !< Delta_H value if prescribed |
---|
27 | integer,dimension(nx,ny) :: i_delta_H !< 1 if Delta_H is prescribed on this node, else 0 |
---|
28 | integer,dimension(nx,ny) :: i_Hp !< 1 if H is prescribed on this node, else 0 |
---|
29 | integer,dimension(nx,ny) :: i_Hp0 !< i_hp mask reference value does not change with time |
---|
30 | integer,dimension(nx,ny) :: MK_gl0 !< mask grounding line initial |
---|
31 | integer,dimension(nx,ny) :: MK_flot0 !< mask float initial |
---|
32 | integer :: voisin !< work variable |
---|
33 | real (kind=kind(0.d0)) :: t_break !< temps : en double precision |
---|
34 | |
---|
35 | ! for prescribed grounding line retreat : ancienne version |
---|
36 | |
---|
37 | integer,parameter :: iretreat_max = 20 !< maximum number of retreat steps |
---|
38 | integer :: iretreat |
---|
39 | real (kind=kind(0.d0)),dimension(iretreat_max) :: time_gl !< time of the retreat steps |
---|
40 | integer, dimension(nx,ny) :: tag_gl !< tag of all the retreat points |
---|
41 | |
---|
42 | !< floating -> 10000 |
---|
43 | !< grounded not gl -> 1000 |
---|
44 | |
---|
45 | integer :: floating_tag = -10000 |
---|
46 | integer :: grounded_tag = 1000 |
---|
47 | |
---|
48 | |
---|
49 | integer, dimension(nx,ny) :: upwind_gl !< 1 if gl can still retreat, else 0 |
---|
50 | real, dimension(nx,ny) :: H_gl !< flotation thickness at the time |
---|
51 | real, dimension(nx,ny) :: dH_gl_dt !< to smoothly decrease ice thickness |
---|
52 | !< for next gl points |
---|
53 | |
---|
54 | ! nouvelle version recul grounding line pour ice2sea |
---|
55 | ! ------------------------------------------------------ |
---|
56 | |
---|
57 | ! real, dimension(nx,ny) :: H_gl !< flotation thickness at the time (deja declare autre version) |
---|
58 | real, dimension(nx,ny) :: H_beggin !< thickness when the point begins to be prescribed |
---|
59 | real, dimension(nx,ny) :: time_beggin !< time when the point begins to be prescribed |
---|
60 | real, dimension(nx,ny) :: time_just_grounded !< time when the point is just grounded |
---|
61 | |
---|
62 | contains |
---|
63 | |
---|
64 | !______________________________________________________________________________________ |
---|
65 | !> initialise prescribe H |
---|
66 | !! calculate the initial grounding line position |
---|
67 | !! @return MK_flot0 and Mk_gl0 |
---|
68 | |
---|
69 | subroutine init_prescribe_H |
---|
70 | |
---|
71 | implicit none |
---|
72 | if (itracebug.eq.1) call tracebug(' Entree dans routine init_prescribe_H') |
---|
73 | |
---|
74 | ! keep the mask of floating points |
---|
75 | where (flot(:,:)) |
---|
76 | MK_flot0(:,:)= 1 |
---|
77 | elsewhere |
---|
78 | MK_flot0(:,:)= 0 |
---|
79 | end where |
---|
80 | |
---|
81 | |
---|
82 | ! determine the initial grounding line |
---|
83 | |
---|
84 | MK_gl0(:,:)=0 |
---|
85 | HP0(:,:)=H0(:,:) |
---|
86 | |
---|
87 | |
---|
88 | do j=2,ny-1 |
---|
89 | do i=2,nx-1 |
---|
90 | voisin = MK_flot0(i-1,j) + MK_flot0(i+1,j) + MK_flot0(i,j-1) + MK_flot0(i,j+1) |
---|
91 | if ((MK_flot0(i,j).eq.0).and.(voisin.gt.0)) then ! grounded and at least one neighbour floating |
---|
92 | MK_gl0(i,j) = 1 ! tagged mask grounding line |
---|
93 | |
---|
94 | end if |
---|
95 | end do |
---|
96 | end do |
---|
97 | |
---|
98 | |
---|
99 | call prescribe_grid_boundary_0 |
---|
100 | |
---|
101 | !call init_prescribe_retreat |
---|
102 | time_gl(1)= 1000000 ! pour un run sans recul |
---|
103 | t_break = 1000000. ! pour un run avec ice shelves prescrits |
---|
104 | |
---|
105 | |
---|
106 | ! Dans le cas du Groenland ice2sea, mk_init = 3 -> noeuds qui ne bougent pas |
---|
107 | |
---|
108 | if (geoplace(1:4).eq.'GI2S') then |
---|
109 | |
---|
110 | do j = 1,ny |
---|
111 | do i = 1,nx |
---|
112 | if (mk_init(i,j).eq.3) then |
---|
113 | i_hp0 (i,j) = 1 |
---|
114 | Hp0 (i,j) = H0(i,j) |
---|
115 | endif |
---|
116 | end do |
---|
117 | end do |
---|
118 | end if |
---|
119 | |
---|
120 | |
---|
121 | |
---|
122 | end subroutine init_prescribe_H |
---|
123 | !______________________________________________________________________________________ |
---|
124 | !> function prescribe_present_H_gl |
---|
125 | !! calculate the initial grounding line position |
---|
126 | !! @return i_hp(:,:) and hp(:,:) |
---|
127 | |
---|
128 | subroutine prescribe_present_H_gl |
---|
129 | |
---|
130 | implicit none |
---|
131 | |
---|
132 | if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_present_H_gl') |
---|
133 | |
---|
134 | where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line |
---|
135 | i_hp(:,:) = 1 ! thickness prescribed to present value |
---|
136 | hp(:,:) = Hp0(:,:) |
---|
137 | end where |
---|
138 | |
---|
139 | |
---|
140 | end subroutine prescribe_present_H_gl |
---|
141 | |
---|
142 | !______________________________________________________________________________________ |
---|
143 | !> function prescribe_fixed_points |
---|
144 | !! fixes thickness for some points |
---|
145 | !! @return i_hp(:,:) and hp(:,:) |
---|
146 | |
---|
147 | subroutine prescribe_fixed_points |
---|
148 | |
---|
149 | where (i_hp0(:,:).eq.1) ! les points i_hp0 le sont pour tout le run |
---|
150 | i_hp(:,:) = i_hp0(:,:) |
---|
151 | Hp(:,:) = Hp0(:,:) |
---|
152 | end where |
---|
153 | end subroutine prescribe_fixed_points |
---|
154 | |
---|
155 | !______________________________________________________________________________________ |
---|
156 | !> function prescribe_paleo_gl_shelf |
---|
157 | !! calculate the grounding line thickness |
---|
158 | !! no ice shelves |
---|
159 | !! @return i_hp(:,:) and hp(:,:) |
---|
160 | |
---|
161 | subroutine prescribe_paleo_gl_shelf |
---|
162 | |
---|
163 | implicit none |
---|
164 | |
---|
165 | if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_present_H_gl') |
---|
166 | |
---|
167 | |
---|
168 | ! noeuds qui doivent être imposés |
---|
169 | |
---|
170 | where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line |
---|
171 | i_hp(:,:) = 1 ! thickness prescribed |
---|
172 | end where |
---|
173 | |
---|
174 | ! valeur imposee |
---|
175 | where (MK_flot0(:,:).eq. 1) ! paleo shelf epaisseur a 1 |
---|
176 | hp(:,:) = 1. |
---|
177 | end where |
---|
178 | |
---|
179 | |
---|
180 | where (MK_gl0(:,:) .eq. 1) |
---|
181 | hp(:,:) = (sealevel-Bsoc(:,:))*row/ro +1. |
---|
182 | end where |
---|
183 | |
---|
184 | |
---|
185 | end subroutine prescribe_paleo_gl_shelf |
---|
186 | |
---|
187 | |
---|
188 | !______________________________________________________________________________________ |
---|
189 | !> function prescribe_grid_boundary |
---|
190 | !! prescribe the grid boundary thickness to 0 |
---|
191 | !! @return i_hp(:,:) and hp(:,:) |
---|
192 | |
---|
193 | subroutine prescribe_grid_boundary_0 |
---|
194 | |
---|
195 | implicit none |
---|
196 | if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_grid_boundary_0') |
---|
197 | |
---|
198 | ! prescribe thickness on the grid boundaries |
---|
199 | i_hp(:,:) = 0 |
---|
200 | |
---|
201 | i_hp(1,:) = 1 |
---|
202 | i_hp(nx,:)= 1 |
---|
203 | hp(1,:) = 0 |
---|
204 | hp(nx,:) = 0 |
---|
205 | |
---|
206 | i_hp(:,1) = 1 |
---|
207 | i_hp(:,ny)= 1 |
---|
208 | hp(:,1) = 0 |
---|
209 | hp(:,ny) = 0 |
---|
210 | |
---|
211 | ! valeurs de reference |
---|
212 | i_hp0(1,:) = 1 |
---|
213 | i_hp0(nx,:) = 1 |
---|
214 | hp0(1,:) = 0 |
---|
215 | hp0(nx,:) = 0 |
---|
216 | |
---|
217 | i_hp0(:,1) = 1 |
---|
218 | i_hp0(:,ny) = 1 |
---|
219 | hp0(:,1) = 0 |
---|
220 | hp0(:,ny) = 0 |
---|
221 | |
---|
222 | |
---|
223 | |
---|
224 | |
---|
225 | end subroutine prescribe_grid_boundary_0 |
---|
226 | |
---|
227 | !______________________________________________________________________________________ |
---|
228 | !______________________________________________________________________________________ |
---|
229 | !> lecture of a prescribe field of delta H |
---|
230 | !! |
---|
231 | !! @return Delta_H, i_delta_H |
---|
232 | |
---|
233 | subroutine lect_prescribed_deltaH |
---|
234 | |
---|
235 | implicit none |
---|
236 | character(len=100) :: delta_H_file ! variation d'epaisseur imposee |
---|
237 | |
---|
238 | if (itracebug.eq.1) call tracebug(' Entree dans routine lect_prescribed_deltaH') |
---|
239 | if (itracebug.eq.1) call tracebug(' Attention ne dvrait pas passer par la') |
---|
240 | ! lecture |
---|
241 | delta_H_file = 'sensibxg-25km.dat' |
---|
242 | delta_H_file = trim(dirnameinp)//trim(delta_H_file) |
---|
243 | |
---|
244 | |
---|
245 | |
---|
246 | |
---|
247 | call lect_input(3,'Delta_H',1,Delta_H,delta_H_file,trim(dirnameinp)//trim(runname)//'.nc') |
---|
248 | !call lect_datfile(nx,ny,Delta_H,1,delta_H_file) |
---|
249 | |
---|
250 | do j=1,ny |
---|
251 | do i=1,nx |
---|
252 | if (Delta_H(i,j).ne.0.) then |
---|
253 | i_delta_H(i,j)=1 |
---|
254 | i_Hp(i,j)=1 ! i_delta_H might be not completely included in i_Hp |
---|
255 | Hp0(i,j) = H(i,j)-Delta_H(i,j) ! warning Delta_H is substracted |
---|
256 | |
---|
257 | if (MK_flot0(i,j).eq. 1) then ! si flotte, pas de changement d'epaisseur |
---|
258 | Hp0(i,j)=H0(i,j) |
---|
259 | else if (Hp0(i,j).lt. -Bsoc(i,j)*row/ro) then ! devient flottant |
---|
260 | Hp0(i,j) = -Bsoc(i,j)*row/ro-10. ! floating thickness - 10 |
---|
261 | MK_flot0(i,j)= 1 ! est tagge masque flottant |
---|
262 | end if |
---|
263 | Hp(i,j) = Hp0(i,j) |
---|
264 | else |
---|
265 | i_delta_H(i,j)=0 |
---|
266 | end if |
---|
267 | end do |
---|
268 | end do |
---|
269 | end subroutine lect_prescribed_deltaH |
---|
270 | |
---|
271 | !______________________________________________________________________________________ |
---|
272 | !> substract Delta_H to the current thickness |
---|
273 | !! |
---|
274 | !! @return |
---|
275 | subroutine prescribe_deltaH |
---|
276 | |
---|
277 | implicit none |
---|
278 | |
---|
279 | if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_deltaH') |
---|
280 | where (i_delta_H(:,:).eq.1) |
---|
281 | Hp(:,:) = Hp0(:,:) |
---|
282 | i_Hp(:,:)=1 |
---|
283 | end where |
---|
284 | end subroutine prescribe_deltaH |
---|
285 | !______________________________________________________________________________________ |
---|
286 | !> break ice shelves |
---|
287 | !! |
---|
288 | !! |
---|
289 | subroutine break_all_ice_shelves |
---|
290 | |
---|
291 | implicit none |
---|
292 | if (itracebug.eq.1) call tracebug(' Entree dans routine break_all_ice_shelves ') |
---|
293 | |
---|
294 | |
---|
295 | debug_3D(:,:,89)=0. |
---|
296 | |
---|
297 | where (flot(:,:)) ! floating from the beginning and eventually from grounding line retreat |
---|
298 | i_hp(:,:) = 1 ! thickness prescribed to 1 m |
---|
299 | hp (:,:) = 1. |
---|
300 | hp0 (:,:) = 1. |
---|
301 | H(:,:) =1. |
---|
302 | debug_3D(:,:,89)=1. |
---|
303 | |
---|
304 | end where |
---|
305 | where (MK_flot0(:,:).eq. 1) |
---|
306 | i_hp(:,:) = 1 ! thickness prescribed to 1 m |
---|
307 | hp (:,:) = 1. |
---|
308 | hp0 (:,:) = 1. |
---|
309 | H(:,:) =1. |
---|
310 | end where |
---|
311 | |
---|
312 | |
---|
313 | end subroutine break_all_ice_shelves |
---|
314 | |
---|
315 | ! |
---|
316 | !______________________________________________________________________________________ |
---|
317 | !> break ice shelves |
---|
318 | !! |
---|
319 | !! |
---|
320 | subroutine melt_ice_shelves |
---|
321 | |
---|
322 | implicit none |
---|
323 | integer :: nbvoisins,nbdeglac,iv,jv |
---|
324 | real :: hmoy,hmin,hmax |
---|
325 | |
---|
326 | if (itracebug.eq.1) call tracebug(' Entree dans routine melt_ice_shelves') |
---|
327 | |
---|
328 | do j = 2,ny-1 |
---|
329 | do i=2, nx-1 |
---|
330 | if (flot(i,j)) then ! floating |
---|
331 | i_hp(i,j) = 1 ! thickness prescribed |
---|
332 | hmoy=0. |
---|
333 | hmin=0. |
---|
334 | hmax=0. |
---|
335 | nbvoisins=0 |
---|
336 | nbdeglac=0 |
---|
337 | do jv=-1,1 |
---|
338 | do iv=1,1 |
---|
339 | if (flot(i+iv,j+jv)) then |
---|
340 | nbvoisins=nbvoisins+1 |
---|
341 | hmoy=hmoy+h(i+iv,j+jv) |
---|
342 | hmin=min(hmin,h(i+iv,j+jv)) |
---|
343 | hmax=max(hmax,h(i+iv,j+jv)) |
---|
344 | if (h(i+jv,j+jv).lt.10.) nbdeglac=nbdeglac+1 |
---|
345 | end if |
---|
346 | end do |
---|
347 | end do |
---|
348 | hmoy=hmoy/(max(nbvoisins,1)) |
---|
349 | |
---|
350 | if(H(i,j).gt.hmoy) then ! melt only if smooth |
---|
351 | hp(i,j)=H(i,j)-1.*dt |
---|
352 | else if (hmax-hmin.lt.10.) then |
---|
353 | hp(i,j)=H(i,j)-1.*dt |
---|
354 | else if ((H(i,j).lt.50.).and.(Hmin.lt.30.)) then |
---|
355 | hp(i,j)=1. |
---|
356 | else if ((H(i,j).lt.50.).and.(nbdeglac.ge.4)) then |
---|
357 | hp(i,j)=1. |
---|
358 | endif |
---|
359 | hp(i,j) = max(hp(i,j),1.) |
---|
360 | end if |
---|
361 | end do |
---|
362 | end do |
---|
363 | |
---|
364 | end subroutine melt_ice_shelves |
---|
365 | |
---|
366 | |
---|
367 | |
---|
368 | ! |
---|
369 | !______________________________________________________________________________________ |
---|
370 | !> prescribe grounding line retreat |
---|
371 | !! |
---|
372 | !! time_gl(iretreat_max) : times of the retreat steps (in) |
---|
373 | !! tag_gl(nx,ny) : tag that gives at which time a node is prescribed |
---|
374 | !! as a grounding line node |
---|
375 | !! @return Hp |
---|
376 | !! @return i_hp |
---|
377 | |
---|
378 | subroutine prescribe_retreat |
---|
379 | |
---|
380 | implicit none |
---|
381 | integer, dimension(nx,ny) :: new_tag ! working array |
---|
382 | |
---|
383 | if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_retreat') |
---|
384 | if (itracebug.eq.1) write(num_tracebug,*) 'time_gl(1)=', time_gl(1) |
---|
385 | |
---|
386 | new_tag(:,:) = tag_gl(:,:) ! initialize new_tag |
---|
387 | |
---|
388 | start_time: if (time.lt.time_gl(1)) then ! before the perturbation |
---|
389 | |
---|
390 | if (itracebug.eq.1) call tracebug(' Avant appel prescribe_present_H_gl') |
---|
391 | |
---|
392 | else |
---|
393 | if (itracebug.eq.1) call tracebug('avant entree step_gl') |
---|
394 | |
---|
395 | iretreat_loop: do iretreat = 1,iretreat_max-1 |
---|
396 | |
---|
397 | if (itracebug.eq.1) write(num_tracebug,*) 'iretreat',iretreat |
---|
398 | |
---|
399 | find_step: if (abs(time-time_gl(iretreat)).lt.dtmin/2.) then |
---|
400 | ! exactly on a retreat step |
---|
401 | ! write(6,*) time, 'sur un step',iretreat |
---|
402 | |
---|
403 | where (tag_gl(:,:).eq.iretreat) |
---|
404 | H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness |
---|
405 | Hp(:,:) = H_gl(:,:)+5. ! +5 to be sure to be grounded |
---|
406 | ! even with small isostatic changes |
---|
407 | i_Hp(:,:) = 1 |
---|
408 | |
---|
409 | elsewhere (tag_gl(:,:).eq.iretreat+1) ! points that will be gl at next step |
---|
410 | |
---|
411 | H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness |
---|
412 | dH_gl_dt(:,:) = (H(:,:)-H_gl(:,:))-5. ! thickness to be reduced to reach gl |
---|
413 | dH_gl_dt(:,:) = dH_gl_dt(:,:) / (time_gl(iretreat+1)-time_gl(iretreat)) |
---|
414 | ! rate of decrease |
---|
415 | |
---|
416 | elsewhere ((tag_gl(:,:).eq.iretreat-1).and. & ! points that become floating |
---|
417 | (upwind_gl(:,:).eq.1)) |
---|
418 | |
---|
419 | H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness |
---|
420 | Hp(:,:) = H_gl(:,:)-10. ! to be sure to be floating |
---|
421 | i_Hp(:,:) = 1 |
---|
422 | new_tag(:,:) = floating_tag |
---|
423 | |
---|
424 | end where |
---|
425 | |
---|
426 | ! pour une raison mysterieuse la boucle suivante marche et celle en where ne marche pas |
---|
427 | ! (memory falut) |
---|
428 | do j=1,ny |
---|
429 | do i=1,nx |
---|
430 | if ((tag_gl(i,j).eq.iretreat-1).and. & ! points that stay gl |
---|
431 | (upwind_gl(i,j).eq.0)) then |
---|
432 | H_gl(i,j) = -Bsoc(i,j)*row/ro ! flotation thickness |
---|
433 | Hp(i,j) = H_gl(i,j)+5. ! to be sure to be grounded |
---|
434 | i_Hp(i,j) = 1 |
---|
435 | new_tag(i,j) = iretreat ! tagged with the new step number |
---|
436 | ! write(6,*) i,j,tag_gl(i,j), upwind_gl(i,j) |
---|
437 | endif |
---|
438 | end do |
---|
439 | end do |
---|
440 | |
---|
441 | !!$ where ((tag_gl(:,:).eq.iretreat-1).and. & ! points that stay gl |
---|
442 | !!$ (upwind_gl(:,:).eq.0)) |
---|
443 | !!$ |
---|
444 | !!$ H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness |
---|
445 | !!$ Hp(:,:) = H_gl(:,:)+5. ! to be sure to be grounded |
---|
446 | !!$ i_Hp(:,:) = 1 |
---|
447 | !!$ new_tag(:,:) = iretreat ! tagged with the new step number |
---|
448 | !!$ end where |
---|
449 | |
---|
450 | else if ((time.gt.time_gl(iretreat)+dtmin/2.).and.(time.lt.time_gl(iretreat+1))) then |
---|
451 | ! write(6,*) time, 'entre',time_gl(iretreat),time_gl(iretreat+1) |
---|
452 | |
---|
453 | where (tag_gl(:,:).eq.iretreat) |
---|
454 | H_gl(:,:) = -Bsoc(:,:)*row/ro ! flotation thickness |
---|
455 | Hp(:,:) = H_gl(:,:)+5. ! +5 to be sure to be grounded |
---|
456 | ! even with small isostatic changes |
---|
457 | i_Hp(:,:) = 1 |
---|
458 | |
---|
459 | elsewhere (tag_gl(:,:).eq.iretreat+1) ! points that will be gl at next step |
---|
460 | Hp(:,:) = H(:,:)-dH_gl_dt(:,:)*(time-time_gl(iretreat)) |
---|
461 | i_Hp(:,:) = 1 |
---|
462 | |
---|
463 | elsewhere (tag_gl(:,:).eq.floating_tag) |
---|
464 | Hp(:,:) = H(:,:) |
---|
465 | i_Hp(:,:) = 1 |
---|
466 | end where |
---|
467 | |
---|
468 | else if (time.gt.time_gl(iretreat_max)) then |
---|
469 | ! write(6,*) time, 'apres iretreat_max',iretreat,time_gl(iretreat_max) |
---|
470 | |
---|
471 | where (tag_gl(:,:).eq.iretreat_max) |
---|
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.floating_tag) |
---|
478 | Hp(:,:) = H(:,:) |
---|
479 | i_Hp(:,:) = 1 |
---|
480 | end where |
---|
481 | |
---|
482 | !faire une sortie de la boucle step_gl |
---|
483 | |
---|
484 | end if find_step |
---|
485 | end do iretreat_loop |
---|
486 | |
---|
487 | end if start_time |
---|
488 | |
---|
489 | tag_gl(:,:) = new_tag(:,:) |
---|
490 | debug_3D(:,:,84) = tag_gl(:,:) |
---|
491 | debug_3D(:,:,85) = i_Hp(:,:) |
---|
492 | debug_3D(:,:,86) = Hp(:,:) |
---|
493 | |
---|
494 | end subroutine prescribe_retreat |
---|
495 | ! |
---|
496 | !______________________________________________________________________________________ |
---|
497 | !> initialise prescribe retreat |
---|
498 | !! |
---|
499 | !! @return tag_gl |
---|
500 | !! @return upwind_gl |
---|
501 | !! @return time_gl |
---|
502 | |
---|
503 | subroutine init_prescribe_retreat |
---|
504 | |
---|
505 | implicit none |
---|
506 | real, dimension(nx,ny) :: lect_tab |
---|
507 | if (itracebug.eq.1) call tracebug(' Entree dans routine init_prescribe_retreat') |
---|
508 | |
---|
509 | ! read the prescribed gl file |
---|
510 | |
---|
511 | call lect_input(3,'lecttab',1,lect_tab,trim(dirnameinp)//'masque_gl_15km_final.dat',trim(dirnameinp)//trim(runname)//'.nc') |
---|
512 | !call lect_datfile(nx,ny,lect_tab,1,trim(dirnameinp)//'masque_gl_15km_final.dat') |
---|
513 | |
---|
514 | tag_gl(:,:)=nint(lect_tab(:,:)) |
---|
515 | |
---|
516 | ! look if neighbors could become grounding line later |
---|
517 | upwind_gl(:,:) = 0 |
---|
518 | do j=2,ny-1 |
---|
519 | do i=2,nx-1 |
---|
520 | if ((tag_gl(i,j).ne.floating_tag).and.(tag_gl(i,j).ne.grounded_tag)) then |
---|
521 | if (any(tag_gl(i-1:i+1,j-1:j+1).eq.tag_gl(i,j)+1)) then |
---|
522 | upwind_gl(i,j) = 1 |
---|
523 | end if |
---|
524 | end if |
---|
525 | end do |
---|
526 | end do |
---|
527 | |
---|
528 | ! time of the retreat steps |
---|
529 | |
---|
530 | do i=1,iretreat_max |
---|
531 | time_gl(i)=(i-1)*30.+100000. |
---|
532 | end do |
---|
533 | time_gl(iretreat)=200000. |
---|
534 | |
---|
535 | write(num_rep_42,*) '! Grounding line retreat prescribed' |
---|
536 | write(num_rep_42,*) '! time_gl begin et suivant : ',time_gl(1),time_gl(2), & |
---|
537 | ' time_glmax=',time_gl(iretreat) |
---|
538 | |
---|
539 | ! beginning of breaking or melting ice shelves |
---|
540 | |
---|
541 | t_break = 1000000.d00 |
---|
542 | |
---|
543 | write(num_rep_42,*) '! t_break=', t_break |
---|
544 | |
---|
545 | end subroutine init_prescribe_retreat |
---|
546 | !______________________________________________________________________________________ |
---|
547 | !______________________________________________________________________________________ |
---|
548 | !> function prescribe_present_H_gl |
---|
549 | !! calculate the initial grounding line position |
---|
550 | !! @return i_hp(:,:) and hp(:,:) |
---|
551 | |
---|
552 | subroutine prescribe_present_H_gl_copy |
---|
553 | |
---|
554 | implicit none |
---|
555 | if (itracebug.eq.1) call tracebug(' Entree dans routine prescribe_present_H_gl-copy') |
---|
556 | |
---|
557 | do j=1,ny |
---|
558 | do i=1,nx |
---|
559 | if ((MK_flot0(i,j).eq. 1).or.(MK_gl0(i,j) .eq. 1)) then |
---|
560 | i_hp(i,j) = 1 ! thickness prescribed to present value |
---|
561 | hp(i,j) = Hp0(i,j) |
---|
562 | end if |
---|
563 | end do |
---|
564 | end do |
---|
565 | |
---|
566 | !!$where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line |
---|
567 | !!$ i_hp(:,:) = 1 ! thickness prescribed to present value |
---|
568 | !!$! hp(:,:) = 1 !Hp0(:,:) |
---|
569 | !!$end where |
---|
570 | |
---|
571 | |
---|
572 | end subroutine prescribe_present_H_gl_copy |
---|
573 | |
---|
574 | !______________________________________________________________________________________ |
---|
575 | !> initialise prescr_ice2sea_retreat |
---|
576 | !! |
---|
577 | !! @return tag_gl |
---|
578 | !! @return upwind_gl |
---|
579 | !! @return time_gl |
---|
580 | |
---|
581 | ! nouvelle routine de recul impose de grounding line pour ice2sea : decembre 2011 |
---|
582 | ! cette routine transforme les informations distance en information temps ou un noeud sera just grounded" |
---|
583 | |
---|
584 | |
---|
585 | subroutine init_prescr_ice2sea_retreat |
---|
586 | |
---|
587 | implicit none |
---|
588 | |
---|
589 | character(len=100) :: dist_file !< fichier contenant les distances |
---|
590 | character(len=100) :: correc_gl_file !< fichier contenant les distances |
---|
591 | |
---|
592 | |
---|
593 | real, dimension(nx,ny) :: dist_gl ! distance à la grounding line ( a lire) |
---|
594 | real, dimension(nx,ny) :: correc_present_line ! nombre de sous-noeuds grounded (a lire) |
---|
595 | |
---|
596 | real :: vretrait ! vitesse de retrait |
---|
597 | real :: time_onset ! moment ou la grounding line commence a reculer |
---|
598 | |
---|
599 | |
---|
600 | namelist/retreat_ice2sea/dist_file,correc_gl_file,vretrait,time_onset |
---|
601 | |
---|
602 | |
---|
603 | if (itracebug.eq.1) call tracebug(' Entree dans routine init_prescr_ice2sea_retreat') |
---|
604 | |
---|
605 | |
---|
606 | ! read the distance file and the |
---|
607 | |
---|
608 | call lect_input(3,'lecttab',1,lect_tab,trim(dirnameinp)//'masque_gl_15km_final.dat',trim(dirnameinp)//trim(runname)//'.nc') |
---|
609 | !call lect_datfile(nx,ny,lect_tab,1,trim(dirnameinp)//'masque_gl_15km_final.dat') |
---|
610 | |
---|
611 | tag_gl(:,:)=nint(lect_tab(:,:)) |
---|
612 | |
---|
613 | ! look if neighbors could become grounding line later |
---|
614 | upwind_gl(:,:) = 0 |
---|
615 | do j=2,ny-1 |
---|
616 | do i=2,nx-1 |
---|
617 | if ((tag_gl(i,j).ne.floating_tag).and.(tag_gl(i,j).ne.grounded_tag)) then |
---|
618 | if (any(tag_gl(i-1:i+1,j-1:j+1).eq.tag_gl(i,j)+1)) then |
---|
619 | upwind_gl(i,j) = 1 |
---|
620 | end if |
---|
621 | end if |
---|
622 | end do |
---|
623 | end do |
---|
624 | |
---|
625 | ! time of the retreat steps |
---|
626 | |
---|
627 | do i=1,iretreat_max |
---|
628 | time_gl(i)=(i-1)*30.+100000. |
---|
629 | end do |
---|
630 | time_gl(iretreat)=200000. |
---|
631 | |
---|
632 | write(num_rep_42,*) '! Grounding line retreat prescribed' |
---|
633 | write(num_rep_42,*) '! time_gl begin et suivant : ',time_gl(1),time_gl(2), & |
---|
634 | ' time_glmax=',time_gl(iretreat) |
---|
635 | |
---|
636 | ! beginning of breaking or melting ice shelves |
---|
637 | |
---|
638 | t_break = 1000000.d00 |
---|
639 | |
---|
640 | write(num_rep_42,*) '! t_break=', t_break |
---|
641 | |
---|
642 | end subroutine init_prescr_ice2sea_retreat |
---|
643 | !______________________________________________________________________________________ |
---|
644 | |
---|
645 | |
---|
646 | end module prescribe_H |
---|