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

Last change on this file since 467 was 467, checked in by aquiquet, 5 months ago

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