source: branches/GRISLIv3/SOURCES/prescribe-H-i2s_mod.f90 @ 431

Last change on this file since 431 was 431, checked in by dumas, 13 months ago

Use only in module prescribe_H
File size: 10.8 KB
Line 
1!> \file prescribe-H-i2s_mod.f90
2!!
3!!  Determine mask and thickness for prescription of H
4!!  include ice2sea prescribed retreat  (Cat april 2012)
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
18module prescribe_H
19
20  ! pour les reculs de type ice2sea
21  !  use toy_retreat_mod
22
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
33  ! pour grounding line retreat, ice2sea
34  ! nouvelle version interactive pour imposer des contraintes variees.
35  ! sanity check, temps de depart en fonction des bassins, ....
36
37
38  ! nouvelle version recul grounding line pour ice2sea
39  ! ------------------------------------------------------
40
41contains
42
43  !______________________________________________________________________________________
44  !> initialise prescribe H
45  !! calculate the initial grounding line position
46  !! @return MK_flot0 and Mk_gl0
47
48  subroutine init_prescribe_H
49
50    use module3D_phy, only: flot,MK_flot0,MK_gl0,Hp0,H0,Mk_init,i_hp0
51    use geography, only: nx,ny,geoplace
52    use runparam, only: itracebug 
53
54    implicit none
55    integer    :: voisin   !< pour test sur des masques
56    integer :: i,j
57
58    if (itracebug.eq.1)  call tracebug(' Entree dans routine init_prescribe_H')
59
60    ! keep the mask of floating points
61    where (flot(:,:))
62       MK_flot0(:,:)= 1
63    elsewhere 
64       MK_flot0(:,:)= 0
65    end where
66
67    ! determine the initial grounding line
68
69    MK_gl0(:,:)=0
70    HP0(:,:)=H0(:,:)
71
72    do j=2,ny-1
73       do i=2,nx-1
74          voisin = MK_flot0(i-1,j) + MK_flot0(i+1,j) + MK_flot0(i,j-1) + MK_flot0(i,j+1)
75          if ((MK_flot0(i,j).eq.0).and.(voisin.gt.0))  then  ! grounded and at least one neighbour floating
76             MK_gl0(i,j) = 1                                ! tagged mask grounding line
77          end if
78       end do
79    end do
80
81    call prescribe_grid_boundary_0
82
83    ! Dans le cas du Groenland ice2sea, mk_init = 3 -> noeuds qui ne bougent pas
84
85    if (geoplace(1:4).eq.'GI2S') then
86
87       do j = 1,ny
88          do i = 1,nx
89             if (mk_init(i,j).eq.3) then
90                i_hp0 (i,j) = 1
91                Hp0  (i,j)  = H0(i,j)
92             endif
93          end do
94       end do
95    end if
96
97  end subroutine init_prescribe_H
98  !______________________________________________________________________________________
99  !>  function prescribe_present_H_gl
100  !! calculate the initial grounding line position
101  !! @return i_hp(:,:) and hp(:,:)
102
103  subroutine prescribe_present_H_gl
104   
105    use module3D_phy, only: MK_flot0,MK_gl0,i_hp,Hp,Hp0
106    use runparam, only: itracebug 
107   
108    implicit none
109
110    if (itracebug.eq.1)  call tracebug(' Entree dans routine  prescribe_present_H_gl')
111
112    where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
113       i_hp(:,:) = 1                                      ! thickness prescribed to present value
114       Hp(:,:)   = Hp0(:,:)
115    end where
116    if (itracebug.eq.1)  call tracebug(' fin prescribe_present_H_gl')
117
118  end subroutine prescribe_present_H_gl
119
120  !______________________________________________________________________________________
121  !>  function prescribe_present_H_gl_bmelt
122  !! calculate the initial grounding line position
123  !! only grounding line points are prescribed to present value
124  !! @return i_hp(:,:) and hp(:,:)
125
126  subroutine prescribe_present_H_gl_bmelt
127   
128    use module3D_phy, only: i_hp
129    use runparam, only: itracebug 
130   
131    implicit none
132
133    if (itracebug.eq.1)  call tracebug(' Entree dans routine  prescribe_present_H_gl_bmelt')
134
135    !    where (MK_gl0(:,:) .eq. 1) ! grounding line only !cdc pour calcule fonte basale
136    !       i_hp(:,:) = 1                                 ! thickness prescribed to present value
137    !       hp(:,:)   = Hp0(:,:)
138    !    end where
139    i_hp(:,:) = 0
140    if (itracebug.eq.1)  call tracebug(' fin prescribe_present_H_gl_bmelt')
141
142  end subroutine prescribe_present_H_gl_bmelt
143
144  !______________________________________________________________________________________
145  !>  function prescribe_fixed_points
146  !!  fixes  thickness for some points
147  !! @return i_hp(:,:) and hp(:,:)
148
149  subroutine prescribe_fixed_points
150   
151    use module3D_phy, only: i_hp0,i_hp,Hp,Hp0
152    use runparam, only: itracebug 
153   
154    if (itracebug.eq.1)  call tracebug(' Entree dans routine  prescribe_fixed_points')
155    where (i_hp0(:,:).eq.1)             ! les points i_hp0 le sont pour tout le run
156       i_hp(:,:) = i_hp0(:,:)
157       Hp(:,:)   = Hp0(:,:)
158    end where
159
160    if (itracebug.eq.1)  call tracebug(' fin prescribe_fixed_points')
161
162  end subroutine prescribe_fixed_points
163
164  !______________________________________________________________________________________
165  !>  function prescribe_paleo_gl_shelf
166  !! calculate the  grounding line thickness
167  !! no ice shelves
168  !! @return i_hp(:,:) and hp(:,:)
169
170  subroutine prescribe_paleo_gl_shelf
171   
172    use module3D_phy, only: MK_flot0,MK_gl0,i_hp,hp,sealevel_2d,Bsoc
173    use runparam, only: itracebug 
174    use param_phy_mod, only: row,ro
175   
176    implicit none
177
178    if (itracebug.eq.1)  call tracebug(' Entree dans routine prescribe_paleo_gl_shelf')
179
180
181    !   noeuds qui doivent être imposés
182
183    where ((MK_flot0(:,:).eq. 1).or.(MK_gl0(:,:) .eq. 1)) ! floating or grounding line
184       i_hp(:,:) = 1                                      ! thickness prescribed
185    end where
186
187    ! valeur imposee
188    where (MK_flot0(:,:).eq. 1)                          ! paleo shelf epaisseur a 1
189       hp(:,:)   = 1.
190    end where
191
192
193    where (MK_gl0(:,:) .eq. 1)
194       hp(:,:) = (sealevel_2d(:,:)-Bsoc(:,:))*row/ro +1.
195    end where
196
197  end subroutine prescribe_paleo_gl_shelf
198
199  !______________________________________________________________________________________
200  !>  function prescribe_grid_boundary
201  !! prescribe the grid boundary thickness to 0
202  !! @return i_hp(:,:) and hp(:,:)
203
204  subroutine prescribe_grid_boundary_0
205   
206    use module3D_phy, only: i_hp,hp,i_hp0,hp0
207    use runparam, only: itracebug 
208    use geography, only: nx,ny
209
210    implicit none
211
212    if (itracebug.eq.1)  call tracebug(' Entree dans routine prescribe_grid_boundary_0')
213
214    ! prescribe thickness on the grid boundaries
215    i_hp(:,:) = 0
216
217    i_hp(1,:) = 1
218    i_hp(nx,:)= 1
219    hp(1,:)   = 0
220    hp(nx,:)  = 0
221
222    i_hp(:,1) = 1
223    i_hp(:,ny)= 1
224    hp(:,1)   = 0
225    hp(:,ny)  = 0
226
227    ! valeurs de reference
228    i_hp0(1,:)  = 1
229    i_hp0(nx,:) = 1
230    hp0(1,:)    = 0
231    hp0(nx,:)   = 0
232
233    i_hp0(:,1)  = 1
234    i_hp0(:,ny) = 1
235    hp0(:,1)    = 0
236    hp0(:,ny)   = 0
237
238  end subroutine prescribe_grid_boundary_0
239  !______________________________________________________________________________________
240
241  !______________________________________________________________________________________
242  !> break ice shelves
243  !!
244  !!
245  subroutine break_all_ice_shelves
246   
247    use module3D_phy, only: i_hp,hp,hp0,H,debug_3D,flot,MK_flot0
248    use runparam, only: itracebug
249   
250    implicit none
251    if (itracebug.eq.1)  call tracebug(' Entree dans routine break_all_ice_shelves ')
252
253
254    debug_3D(:,:,89)=0.
255
256    where (flot(:,:))                           ! floating from the beginning and eventually from grounding line retreat
257       i_hp(:,:) = 1                                      ! thickness prescribed to 1 m
258       hp (:,:)  = 1.     
259       hp0 (:,:)  = 1.   
260       H(:,:)     =1.
261       debug_3D(:,:,89)=1.
262
263    end where
264    where (MK_flot0(:,:).eq. 1)
265       i_hp(:,:) = 1                                      ! thickness prescribed to 1 m
266       hp (:,:)  = 1. 
267       hp0 (:,:)  = 1.   
268       H(:,:)     =1.
269    end where
270
271  end subroutine break_all_ice_shelves
272
273  !
274  !______________________________________________________________________________________
275  !> break ice shelves
276  !!
277  !!
278  subroutine melt_ice_shelves
279   
280    use module3D_phy, only: i_hp,H,hp,dt,flot
281    use runparam, only: itracebug 
282    use geography, only: nx,ny
283   
284    implicit none
285   
286    integer :: nbvoisins,nbdeglac,iv,jv
287    real    :: hmoy,hmin,hmax
288    integer :: i,j
289
290    if (itracebug.eq.1)  call tracebug(' Entree dans routine melt_ice_shelves')
291
292    do j = 2,ny-1
293       do i=2, nx-1
294          if  (flot(i,j)) then                                 ! floating
295             i_hp(i,j) = 1                                      ! thickness prescribed
296             hmoy=0.
297             hmin=0.
298             hmax=0.
299             nbvoisins=0
300             nbdeglac=0
301             do jv=-1,1
302                do iv=1,1
303                   if (flot(i+iv,j+jv)) then
304                      nbvoisins=nbvoisins+1
305                      hmoy=hmoy+h(i+iv,j+jv)
306                      hmin=min(hmin,h(i+iv,j+jv))
307                      hmax=max(hmax,h(i+iv,j+jv))
308                      if (h(i+jv,j+jv).lt.10.) nbdeglac=nbdeglac+1
309                   end if
310                end do
311             end do
312             hmoy=hmoy/(max(nbvoisins,1))
313
314             if(H(i,j).gt.hmoy) then  ! melt only if smooth
315                hp(i,j)=H(i,j)-1.*dt
316             else if (hmax-hmin.lt.10.) then
317                hp(i,j)=H(i,j)-1.*dt
318             else if ((H(i,j).lt.50.).and.(Hmin.lt.30.)) then
319                hp(i,j)=1.
320             else if ((H(i,j).lt.50.).and.(nbdeglac.ge.4)) then
321                hp(i,j)=1.
322             endif
323             hp(i,j)   = max(hp(i,j),1.)
324          end if
325       end do
326    end do
327
328  end subroutine melt_ice_shelves
329
330
331  !______________________________________________________________________________________
332  !>  function prescribe_present_H_gl
333  !! calculate the initial grounding line position
334  !! @return i_hp(:,:) and hp(:,:)
335
336  subroutine prescribe_present_H_gl_copy
337   
338    use module3D_phy, only: i_hp,hp,Hp0,MK_flot0,MK_gl0
339    use runparam, only: itracebug 
340    use geography, only: nx,ny
341   
342    implicit none
343    integer :: i,j
344   
345    if (itracebug.eq.1)  call tracebug(' Entree dans routine  prescribe_present_H_gl-copy')
346
347    do j=1,ny
348       do i=1,nx
349          if ((MK_flot0(i,j).eq. 1).or.(MK_gl0(i,j) .eq. 1)) then
350             i_hp(i,j) = 1                                      ! thickness prescribed to present value
351             hp(i,j)   = Hp0(i,j)
352          end if
353       end do
354    end do
355
356  end subroutine prescribe_present_H_gl_copy
357  !______________________________________________________________________________________
358
359
360end module prescribe_H
Note: See TracBrowser for help on using the repository browser.