source: trunk/SOURCES/calving_frange_ISMIP_glaciers.f90 @ 334

Last change on this file since 334 was 294, checked in by aquiquet, 4 years ago

Calving version used for ISMIP6-Greenland
File size: 25.4 KB
Line 
1!> \file calving_frange_glaciers.f90
2!! Module et routines qui calculent le calving avec la methode de Vincent
3! Pour experience Ice Shelf glaciers retreat ISMIP6 Greenland
4!<
5
6
7!> \namespace module3D_phy
8!! Module et routines qui calculent le calving avec la methode de Vincent
9!! \author ...
10!! \date ...
11!! @note recopie direct de icethick5-ant. A revoir en particulier, accroche avec la glace posee
12!! @note Used modules
13!! @note   - use module3D_phy
14!! @todo il faudrait que le calving soit inactif dans les zones ou l epaisseur est imposee
15!<
16
17module calving_frange_glaciers
18
19  use module3D_phy
20  use netcdf
21  use io_netcdf_grisli
22  use bilan_eau_mod
23  implicit none
24
25  real, dimension (nx,ny) :: hmhc  ! hauteur au dessus de la coupure
26  real, dimension (nx,ny) :: bil_tot ! bilan surface et fond (bm-bmelt)
27  real, dimension (nx,ny) :: hcoup ! epaiseur de coupure au temps time
28  real :: hcoup_plateau  ! coupure points peu profonds
29  real :: hcoup_abysses  ! coupure points ocean profond
30  real :: prof_plateau  ! profondeur max des points peu profonds
31  real :: prof_abysses  ! profondeur min des points ocean profond
32  integer :: meth_hcoup  ! pour avoir hcoup dépendant du coefbmshelf
33  integer :: ifrange        !  0 pas de traitement particulier pres du bord, 1 -> franges
34  integer ::  iin2,jin2,iin3,jin3             ! pour la detection polynies
35  logical :: testmij,testpij,testimj,testipj
36  logical :: bilan_surf_fond                  ! vrai si bm-bmelt est positif
37  logical :: avalw,avale,avals,avaln,interieur
38  real,dimension(:),allocatable        :: time_snap     !> date des snapshots  indice : nb de nb_snap
39  integer                              :: nb_snap       !> nombre de snapshots
40  real,dimension(:,:,:),allocatable :: mask_glaciers !> mask zones glaciers retreat : indices nx, ny, nb_snap
41  real, dimension (nx,ny) :: href_glaciers !> epaisseur de la topo de reference
42 
43contains
44
45  !---------------------------------------------------------------------------------------
46  subroutine init_calving
47
48    character(len=100) :: file_IS_glaciers      !> nom du fichier avec les glaciers retreat
49    character(len=100) :: file_IS_toporef      !> nom du fichier avec les glaciers retreat
50    real               :: time_depart_snaps     !> temps du debut premier snapshot
51    ! pour les lectures ncdf
52    real*8, dimension(:), pointer        :: tab1d => null()         !< tableau 1d real pointer
53    real*8, dimension(:,:), pointer      :: tab2d => null()         !< tableau 2d real pointer
54    real*8, dimension(:,:,:), pointer    :: tab3d => null()         !< tableau 3d real pointer
55
56    namelist/calving_glaciers/Hcoup_plateau,Hcoup_abysses,prof_plateau,prof_abysses,ifrange,meth_hcoup,file_IS_glaciers,file_IS_toporef,nb_snap,time_depart_snaps
57
58    calv(:,:)=0. 
59    calv_dtt(:,:)=0.
60
61    ! formats pour les ecritures dans 42
62428 format(A)
63
64    ! lecture des parametres du run                      block eaubasale1
65    !--------------------------------------------------------------------
66    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
67    read(num_param,calving_glaciers)
68
69
70    write(num_rep_42,428)'!___________________________________________________________' 
71    write(num_rep_42,428) '&calving_glaciers'
72    write(num_rep_42,calving_glaciers)
73    write(num_rep_42,428) '! Hcoup epaisseurs de coupure pour les zones peu prodondes et profondes'
74    write(num_rep_42,428) '! Hcoup_plateau<Hcoup_abysses && prof_plateau<prof_abysses'
75    write(num_rep_42,428) '! prof profondeur delimitant les zones peu prodondes et profondes'
76    write(num_rep_42,428) '! ifrange=0 -> pas de traitement particulier sur les bords'
77    write(num_rep_42,428) '! ifrange=1 -> traitement de Vincent avec ice shelves frangeants partout'
78    write(num_rep_42,428) '! ifrange=2 -> ice shelves frangeants seulement si bm-bmelt positif'
79    write(num_rep_42,428) '! meth_hcoup pour faire eventuellement varier Hcoup avec le climat'
80    write(num_rep_42,428) '! file_IS_glaciers = fichier glaciers retreat'
81    write(num_rep_42,428) '! file_IS_toporef  = fichier avec epaisseur de reference'
82    write(num_rep_42,428) '! nb_snap          = nombre de snapshots'
83    write(num_rep_42,428) '! time_depart_snaps = time depart fichier glaciers'
84
85    ! afq -- coupure depend de la profondeur:
86    !
87    !  hcoup                       prof_abysses
88    !   ^                               v
89    !   |                                _______ hcoup_abysses
90    !   |                               /
91    !   |                              /
92    !   |                             /
93    !   |                            /
94    !   |      hcoup_plateau _______/
95    !                               ^
96    !                          prof_plateau
97   
98    ! lecture fichier glaciers
99    file_IS_glaciers  = trim(dirnameinp)//trim(file_IS_glaciers)
100    file_IS_toporef   = trim(dirnameinp)//trim(file_IS_toporef)
101   
102    ! allocation dynamique de time_snap
103    if (allocated(time_snap)) then
104      deallocate(time_snap,stat=err)
105      if (err/=0) then
106        print *,"Erreur à la desallocation de time_snap",err
107        stop
108      end if
109    end if
110
111    allocate(time_snap(nb_snap),stat=err)
112    if (err/=0) then
113      print *,"erreur a l'allocation du tableau time_snap ",err
114      print *,"nb_snap = ",nb_snap
115      stop
116    end if
117! allocation dynamique de smb_snap
118    if (allocated(mask_glaciers)) then
119      deallocate(mask_glaciers,stat=err)
120      if (err/=0) then
121        print *,"Erreur à la desallocation de mask_glaciers",err
122        stop
123      end if
124    end if
125
126    allocate(mask_glaciers(nx,ny,nb_snap),stat=err)
127    if (err/=0) then
128      print *,"erreur a l'allocation du tableau mask_glaciers",err
129      print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap
130      stop
131    end if
132
133! lecture de mask_glaciers
134    call Read_Ncdf_var('sftgif',file_IS_glaciers,tab3d)
135    mask_glaciers(:,:,:) = tab3d(:,:,:)
136
137! lecture de time_snap
138    call Read_Ncdf_var('time',file_IS_glaciers,tab1d)
139    time_snap(:) = tab1d(:)
140    time_snap(:) = time_snap(:) + time_depart_snaps
141
142! lecture epaisseur de reference
143    call Read_Ncdf_var('H',file_IS_toporef,tab2d)
144    href_glaciers(:,:) = tab2d(:,:)
145    href_glaciers(:,:) = -100.
146
147    Hcoup(:,:) = min ( max(                                         &
148       (-(Bsoc0(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    &
149       *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               &
150       , hcoup_plateau), hcoup_abysses )
151
152    if (meth_hcoup.eq.1) then
153        Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
154        Hcoup(:,:)=min( max(Hcoup (:,:),Hcoup_plateau),Hcoup_abysses)
155    else if (meth_hcoup.eq.2) then
156        Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
157        Hcoup(:,:)=max(Hcoup(:,:),Hcoup_plateau)
158    else if (meth_hcoup.eq.3) then
159        Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
160        Hcoup(:,:)=min(max(Hcoup(:,:),0.),Hcoup_abysses)   
161    endif
162
163  end subroutine init_calving
164  !---------------------------------------------------------------------------------------
165  subroutine calving
166
167    integer :: I_did_something  ! pour la boucle sur le calving
168    integer :: k
169    real,dimension(nx,ny) :: mask_glaciers_time ! mask_glaciers au pas de temps courant
170
171    ! selection du mask_glaciers fonction de time
172    ! /!\ the outputs are written after a dt... so, to be consistent we have to add a dt
173    if(time.lt.(time_snap(1)+dt)) then              ! time avant le forcage
174      mask_glaciers_time(:,:) = 0 ! pas de mask
175    else if (time.ge.(time_snap(nb_snap)+dt)) then  ! time apres le forcage
176      mask_glaciers_time(:,:) = mask_glaciers(:,:,nb_snap) ! mask du dernier pas de temps
177    else ! cas general
178      do k = 1 , nb_snap-1
179        if((time.ge.(time_snap(k)+dt)).and.(time.lt.(time_snap(k+1)+dt))) then ! entre k et k+1
180!          print*,'time et k',time,k, time_snap(k),time_snap(k+1)
181          mask_glaciers_time(:,:)=mask_glaciers(:,:,k)          ! mask de time k
182          exit
183        endif
184     end do
185  endif
186
187    ! initialisation calving
188    calv(:,:)=0. 
189
190    ! calcul du dhdt lagrangien
191    dhdt(:,:)=bm(:,:)-bmelt(:,:)-H(:,:)*(epsxx(:,:)+epsyy(:,:))
192
193    ! calcul du bilan surface et fond : divise par 2 car utilise dans des moyennes
194    bil_tot(:,:)=0.5*(bm(:,:)-bmelt(:,:))
195
196    Hcoup(:,:) = min ( max(                                         &
197       (-(Bsoc(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    &
198       *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               &
199       , hcoup_plateau), hcoup_abysses )
200
201    if (meth_hcoup.eq.1) then
202        Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
203        Hcoup(:,:)=min( max(Hcoup (:,:),Hcoup_plateau),Hcoup_abysses)
204    else if (meth_hcoup.eq.2) then
205        Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
206        Hcoup(:,:)=max(Hcoup(:,:),Hcoup_plateau)
207    else if (meth_hcoup.eq.3) then
208        Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
209        Hcoup(:,:)=min(max(Hcoup(:,:),0.),Hcoup_abysses)
210    endif
211
212    ! hauteur au dessus de la coupure
213    hmhc(:,:)=H(:,:)-Hcoup(:,:)
214
215    ! coupure de l'ice shelf
216    !---------------------------
217    ! on divise le test en 2 'ifext' et 'ifint' pour reduire les tests redondants 
218
219
220
221    MULTI_CALV_LOOP : do l=1,10
222      ! calcul du dhdt lagrangien
223      dhdt(:,:)=bm(:,:)-bmelt(:,:)-H(:,:)*(epsxx(:,:)+epsyy(:,:))
224      ! hauteur au dessus de la coupure
225      hmhc(:,:)=H(:,:)-Hcoup(:,:)
226      I_did_something = 0
227      do j=2,ny-1
228        do i=2,nx-1
229
230
231          ifext:  if((flot(i,j)).and.(h(i,j).le.hcoup(i,j)).and.(h(i,j).gt.0.))  then
232              ! ifext: pour les noeuds flottants englaces avec h < hcoup
233
234              !ifint:    if((front(i,j).gt.0).and.(front(i,j).lt.4)) then
235              !cdc pb avec front    ifint:    if((front(i,j).lt.4).or.((front(i-1,j)+front(i+1,j)+front(i,j-1)+front(i,j+1)).lt.16)) then
236              ifint:    if((H(i-1,j).lt.2.).or.(H(i+1,j).lt.2.).or.(H(i,j-1).lt.2.).or.(H(i,j+1).lt.2)) then ! si on est au bord avec test sur H
237                  ! ifint: le point doit avoir au - 1 voisin mais ne pas etre entouré de glace
238                  ! ce qui evite la formation des polynies dans les shelfs
239
240
241                  ! on regarde dhdt minimum pour voir si 1 des voisins peut alimenter le point
242
243                  !                     hmhc est l'épaisseur en plus de l'épaisseur de coupure
244                  !                     ux/dx=1/(deltat) ou deltat est le temps nécessaire
245                  !                     à la glace pour passer d'un noeud à l'autre
246                  !                     la deuxieme partie du test revient à dh/dt*deltat > hmhc 
247
248                  !                     Rajout vince : si le point amont est pose -> test vrai.
249
250                  ! rappel des differents ifrange
251                  !-------------------------------
252                  !Attention ifrange 1,2  ont sans doute un bug (il faudrait abs(uxbar))
253
254                  ! ifrange=1               Rajout vince : si un point voisin est pose -> test vrai.
255                  !                         comme dans article fenno2007
256
257                  ! ifrange=2               pour les points ayant des voisins poses, garde le point
258                  !                         seulement si le bilan  surface fond est positif.
259                  !                         bm(i,j)-bmelt(i,j) > 0.
260
261                  ! ifrange=3               pas de test sur l'epaisseur du point amont.
262                  !                         test dhdt > -hmhc/deltat= -hmhc abs(U)/dx (correction du bug)
263                  !                         + test sur bm(i,j)-bmelt(i,j) > 0.
264
265                  ! ifrange = 4             idem 3 mais avec test sur l'epaisseur du point amont
266
267                  ! ifrange = 0             pas de traitement specifique près du continent
268
269                  ! ifrange = -1            ancienne version MIS11
270
271                  ! ifrange = 5             semble idem 0
272
273                  ! ifrange = 6             ??
274
275                  ! ifrange = 7             calving drastique sauf si coefbmshelf < 0.5
276
277
278                  if (ifrange.eq.1) then           !   Rajout vince : si un point voisin est pose -> test vrai.
279                      !   comme dans article fenno2007
280
281                      testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
282                         .and.(dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx)))& 
283                         .or.(.not.flot(i-1,j)) ) !
284
285                      testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
286                         .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) &
287                         .or.(.not.flot(i+1,j)) ) !
288
289                      testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
290                         .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))&
291                         .or.(.not.flot(i,j-1)) ) !
292
293                      testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
294                         .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))&
295                         .or.(.not.flot(i,j+1)) ) !
296
297                  else if (ifrange.eq.2) then       ! pour les points ayant des voisins posés, seulement si le bilan
298                      ! surface fond est positif. bm(i,j)-bmelt(i,j) > 0.
299
300                      bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)
301
302                      testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
303                         .and.  (dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx))) & 
304                         .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !
305
306                      testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
307                         .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) &
308                         .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !
309
310                      testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
311                         .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))&
312                         .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !
313
314                      testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
315                         .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))&
316                         .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !
317
318                  else if (ifrange.eq.3) then       ! nouvelle formulation Cat mars 08
319
320                      ! pas de test sur l'epaisseur du point amont.
321                      ! test dhdt > -hmhc/deltat= -hmhc abs(U)/dx
322
323                      ! pour les points ayant des voisins posés, seulement si le bilan
324                      ! surface fond est positif. bm(i,j)-bmelt(i,j) > 0.
325
326
327                      bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)
328
329                      testmij=( (uxbar(i,j).ge.0.) &  ! voisin (i-1,j) amont et > hcoup
330                         .and.(dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j))/dx))) & 
331                         .or.((.not.flot(i-1,j)).and.bilan_surf_fond ) !
332
333                      testpij=( (uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
334                         .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx))) &
335                         .or.((.not.flot(i+1,j)).and.bilan_surf_fond )  !
336
337                      testimj=(( uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
338                         .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx))) &
339                         .or.((.not.flot(i,j-1)).and.bilan_surf_fond )  !
340
341                      testipj=((uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
342                         .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx))) &
343                         .or.((.not.flot(i,j+1)).and.bilan_surf_fond )  !
344
345                  else if (ifrange.eq.4) then       ! idem 3 mais avec test sur l'épaisseur du point amont
346
347                      bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)
348
349                      testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
350                         .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & 
351                         .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !
352
353                      testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
354                         .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) &
355                         .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !
356
357                      testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
358                         .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))&
359                         .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !
360
361                      testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
362                         .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))&
363                         .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !
364
365                  else if (ifrange.eq.0) then           ! pas de traitement special pres du continent
366
367                      testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
368                         .and.(dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))
369
370                      testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
371                         .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))
372
373                      testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
374                         .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))
375
376                      testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
377                         .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))
378
379                  else if (ifrange.eq.-1) then           ! ancienne version MIS11
380
381                      testmij=((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
382                         .and.(dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx)))   
383
384                      testpij=((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
385                         .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx)))
386
387                      testimj=((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
388                         .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))
389
390                      testipj=((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
391                         .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))
392
393
394                  else if (ifrange.eq.5) then           ! pas de traitement special pres du continent
395
396                      testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
397                         .and.((bil_tot(i-1,j)+bil_tot(i,j)) &
398                         .gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))
399
400                      testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
401                         .and.((bil_tot(i+1,j)+bil_tot(i,j)) &
402                         .gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))
403
404                      testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
405                         .and.((bil_tot(i,j-1)+bil_tot(i,j)) &
406                         .gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))
407
408                      testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
409                         .and.((bil_tot(i,j+1)+bil_tot(i,j)) &
410                         .gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))
411
412                  else if (ifrange.eq.6) then           ! pas de traitement special pres du continent
413
414                      testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
415                         .and.(2.*bil_tot(i,j) &
416                         .gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))
417
418                      testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
419                         .and.(2.*bil_tot(i,j) &
420                         .gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))
421
422                      testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
423                         .and.(2.*bil_tot(i,j) &
424                         .gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))
425
426                      testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
427                         .and.(2.*bil_tot(i,j) &
428                         .gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))
429
430                  else if (ifrange.eq.7) then           ! drastique calving sauf si coefbmelt.lt.0.5
431
432                      testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
433                         .and.(coefbmshelf.lt.1.))
434                      testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
435                         .and.(coefbmshelf.lt.1.))
436                      testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
437                         .and.(coefbmshelf.lt.1.))
438                      testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
439                         .and.(coefbmshelf.lt.1.))
440
441
442                  endif
443
444
445
446                  ! detection des polynies dans les shelfs. on regarde vers l'aval
447
448                  ! 1 des 3 voisins aval > hcoup du cote west (i-1)
449                  iin2=max(1,i-2)
450                  iin3=max(1,i-3)
451
452                  avalw=((hmhc(i-1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) &
453                     .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i,j).le.0.)                         
454
455
456                  ! 1 des 3 voisins aval > hcoup du cote est (i+1)
457                  iin2=min(i+2,nx)
458                  iin3=min(i+3,nx)
459
460                  avale=((hmhc(i+1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) &   
461                     .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i+1,j).ge.0)                         
462
463                  ! 1 des 3 voisins aval > hcoup du cote sud (j-1)
464                  jin2=max(1,j-2)
465                  jin3=max(1,j-3)
466
467                  avals=((hmhc(i,j-1).gt.0.).or.(hmhc(i,jin2).gt.0.) &   
468                     .or.(hmhc(i,jin3).gt.0.)).and.(uybar(i,j).le.0.) 
469
470                  ! 1 des 3 voisins aval > hcoup du cote nord (j-1)
471                  jin2=min(j+2,ny)
472                  jin3=min(j+3,ny)
473
474                  avaln=((hmhc(i,j+1).gt.0.).or.(hmhc(i,jin2).gt.0.) &   
475                     .or.(hmhc(i,jin2).gt.0.)) .and.(uybar(i,j+1).ge.0.)                         
476
477                  interieur=(avalw.or.avale).and.(avals.or.avaln)
478
479
480
481                  if ((.not.(testmij.or.testpij.or.testimj.or.testipj))  &  ! pas suffisament alimente
482                     .and.(.not.interieur)) then                ! et pas interieur
483                      calv(i,j)=-h(i,j) 
484                      !cdc                     H(i,j)=1.
485                      !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel_2d(i,j) - Bsoc(i,j))*row/ro-0.01))
486                      H(i,j)=0.
487                      S(i,j)=H(i,j)*(1.-ro/row) + sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
488                      B(i,j)=S(i,j) - H(i,j)
489                      ! ATTENTION ne pas mettre ice=0 sinon degradation bilan d'eau (bm et bmelt non comptabilises dans ce cas)
490
491                      I_did_something = I_did_something + 1
492                      !   if (l.ge.2) then
493                      !     print*,'calving l ij',l,i,j
494                      !   endif
495                  endif
496
497              end if ifint
498          end if ifext
499        end do
500      end do
501      if (I_did_something.eq.0) then
502          !        print*,'stop MULTI_CALV_LOOP l=',l
503          EXIT MULTI_CALV_LOOP
504      else
505          !        print*,'calving continue l',l,I_did_something
506      endif
507    enddo MULTI_CALV_LOOP
508
509
510    ! on met en calving les points detectes iceberg et les points mask_glaciers :
511    where ((iceberg(:,:).and.(H(:,:).gt.0.)))
512        calv(:,:)=-h(:,:)
513        ice(:,:)=0
514        H(:,:)=0.
515        S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
516        B(:,:)=S(:,:) - H(:,:)
517    endwhere
518    !where ( mask_glaciers_time(:,:)*href_glaciers(:,:) .lt. H(:,:) )
519        !calv(:,:) = mask_glaciers_time(:,:)*href_glaciers(:,:) - H(:,:)
520        !H(:,:) = mask_glaciers_time(:,:)*href_glaciers(:,:)
521    !endwhere
522   
523    where ( (href_glaciers(:,:).lt.-10.) .and. (mask_glaciers_time(:,:).lt.1.) )
524       href_glaciers(:,:) = H(:,:)
525    endwhere
526    where ( (H(:,:) .gt. mask_glaciers_time(:,:)*href_glaciers(:,:)) .and. (mask_glaciers_time(:,:) .lt. 1.) .and. (mask_glaciers(:,:,1).gt.0.) )
527       calv(:,:) = -(H(:,:) - mask_glaciers_time(:,:)*href_glaciers(:,:))
528       H(:,:)    = mask_glaciers_time(:,:)*href_glaciers(:,:)
529    endwhere
530
531    where ( (mask_glaciers_time(:,:) .le. 0. ) .and. ( mask_glaciers(:,:,1) .gt. 0. ) )
532        ice(:,:) = 0
533    endwhere
534
535    calv_dtt(:,:) = calv_dtt(:,:) + calv(:,:) ! somme du calving sur dtt 
536    ! calv_dtt est remis a 0 dans steps_time_loop (tous les dtt)
537  end subroutine calving
538  !------------------------------------------------------------------------------------------
539end module calving_frange_glaciers
540
Note: See TracBrowser for help on using the repository browser.