source: branches/iLoveclim/SOURCES/calving_frange.f90 @ 187

!> \file calving_frange.f90
!! Module et routines qui calculent le calving avec la methode de Vincent
!<


!> \namespace module3D_phy
!! Module et routines qui calculent le calving avec la methode de Vincent
!! \author ...
!! \date ...
!! @note recopie direct de icethick5-ant. A revoir en particulier, accroche avec la glace posee
!! @note Used modules
!! @note   - use module3D_phy
!! @todo il faudrait que le calving soit inactif dans les zones ou l epaisseur est imposee
!<

module calving_frange

use module3D_phy
use bilan_eau_mod
implicit none

real, dimension (nx,ny) :: hmhc  ! hauteur au dessus de la coupure
real, dimension (nx,ny) :: bil_tot ! bilan surface et fond (bm-bmelt)
real, dimension (nx,ny) :: hcoup ! epaiseur de coupure au temps time
real :: hcoup_plateau  ! coupure points peu profonds
real :: hcoup_abysses  ! coupure points ocean profond
real :: prof_plateau  ! profondeur max des points peu profonds
real :: prof_abysses  ! profondeur min des points ocean profond
integer :: meth_hcoup  ! pour avoir hcoup dépendant du coefbmshelf
integer :: ifrange        !  0 pas de traitement particulier pres du bord, 1 -> franges
integer ::  iin2,jin2,iin3,jin3             ! pour la detection polynies
logical :: testmij,testpij,testimj,testipj
logical :: bilan_surf_fond                  ! vrai si bm-bmelt est positif
logical :: avalw,avale,avals,avaln,interieur

real :: mass_total
integer :: testH ! somme des flux entrants en un point ij

contains

!---------------------------------------------------------------------------------------
subroutine init_calving

namelist/calving/Hcoup_plateau,Hcoup_abysses,prof_plateau,prof_abysses,ifrange,meth_hcoup



mass_total = 0d0
calv(:,:)=0. 
calv_dtt(:,:)=0.

! formats pour les ecritures dans 42
428 format(A)

! lecture des parametres du run                      block eaubasale1
!--------------------------------------------------------------------
rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
read(num_param,calving)


write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&calving              ! nom du bloc calving méthode Vincent '
write(num_rep_42,*)
write(num_rep_42,*)   'Hcoup_plateau  = ', Hcoup_plateau
write(num_rep_42,*)   'Hcoup_abysses  = ', Hcoup_abysses
write(num_rep_42,*)   'prof_plateau  = ', prof_plateau
write(num_rep_42,*)   'prof_abysses  = ', prof_abysses
write(num_rep_42,*)   'ifrange    = ', ifrange
write(num_rep_42,*)   'meth_hcoup = ', meth_hcoup
write(num_rep_42,*)'/'
                            
write(num_rep_42,428) '! Hcoup epaisseurs de coupure pour les zones peu prodondes et profondes'
write(num_rep_42,428) '! Hcoup_plateau<Hcoup_abysses && prof_plateau<prof_abysses'
write(num_rep_42,428) '! prof profondeur delimitant les zones peu prodondes et profondes'
write(num_rep_42,428) '! ifrange=0 -> pas de traitement particulier sur les bords'
write(num_rep_42,428) '! ifrange=1 -> traitement de Vincent avec ice shelves frangeants partout'
write(num_rep_42,428) '! ifrange=2 -> ice shelves frangeants seulement si bm-bmelt positif'
write(num_rep_42,*)   '! meth_hcoup pour faire eventuellement varier Hcoup avec le climat'
write(num_rep_42,*)   '! Hcoup_plateau=',Hcoup_plateau
write(num_rep_42,*)   '! Hcoup_abysses=',Hcoup_abysses
write(num_rep_42,*)   '! prof_plateau=',prof_plateau
write(num_rep_42,*)   '! prof_abysses=',prof_abysses
write(num_rep_42,*)


! afq -- coupure depend de la profondeur:
!
!  hcoup                       prof_abysses
!   ^                               v 
!   |                                _______ hcoup_abysses
!   |                               /
!   |                              /
!   |                             /
!   |                            /
!   |      hcoup_plateau _______/
!                               ^
!                          prof_plateau     

Hcoup(:,:) = min ( max(                                         &
     (-(Bsoc0(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    &
     *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               &
     , hcoup_plateau), hcoup_abysses )

if (meth_hcoup.eq.1) then
   Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
   Hcoup(:,:)=min( max(Hcoup (:,:),Hcoup_plateau),Hcoup_abysses)
else if (meth_hcoup.eq.2) then
   Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
   Hcoup(:,:)=max(Hcoup(:,:),Hcoup_plateau)
else if (meth_hcoup.eq.3) then
   Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
   Hcoup(:,:)=min(max(Hcoup(:,:),0.),Hcoup_abysses) 
endif


end subroutine init_calving
!---------------------------------------------------------------------------------------
subroutine calving
  
  integer :: I_did_something  ! pour la boucle sur le calving

! initialisation calving
  calv(:,:)=0. 
  mass_total=0d0

! calcul du dhdt lagrangien
dhdt(:,:)=bm(:,:)-bmelt(:,:)-H(:,:)*(epsxx(:,:)+epsyy(:,:))

! calcul du bilan surface et fond : divise par 2 car utilise dans des moyennes
bil_tot(:,:)=0.5*(bm(:,:)-bmelt(:,:))

Hcoup(:,:) = min ( max(                                         &
     (-(Bsoc(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    &
     *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               &
     , hcoup_plateau), hcoup_abysses )

if (meth_hcoup.eq.1) then
   Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
   Hcoup(:,:)=min( max(Hcoup (:,:),Hcoup_plateau),Hcoup_abysses)
else if (meth_hcoup.eq.2) then
   Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
   Hcoup(:,:)=max(Hcoup(:,:),Hcoup_plateau)
else if (meth_hcoup.eq.3) then
   Hcoup(:,:)=coefbmshelf*Hcoup(:,:)
   Hcoup(:,:)=min(max(Hcoup(:,:),0.),Hcoup_abysses)
endif

! hauteur au dessus de la coupure
hmhc(:,:)=H(:,:)-Hcoup(:,:)

! coupure de l'ice shelf
!---------------------------
! on divise le test en 2 'ifext' et 'ifint' pour reduire les tests redondants  


    MULTI_CALV_LOOP : do l=1,10
      ! calcul du dhdt lagrangien
      dhdt(:,:)=bm(:,:)-bmelt(:,:)-H(:,:)*(epsxx(:,:)+epsyy(:,:))
      ! hauteur au dessus de la coupure
      hmhc(:,:)=H(:,:)-Hcoup(:,:)
      I_did_something = 0
      do j=2,ny-1
        do i=2,nx-1


          ifext:  if((flot(i,j)).and.(h(i,j).le.hcoup(i,j)).and.(h(i,j).gt.0.))  then
! ifext: pour les noeuds flottants avec h < hcoup 

!ifint:    if((front(i,j).gt.0).and.(front(i,j).lt.4)) then
!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
        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
! ifint: le point doit avoir au - 1 voisin mais ne pas etre entouré de glace
! ce qui evite la formation des polynies dans les shelfs


! on regarde dhdt minimum pour voir si 1 des voisins peut alimenter le point

!                     hmhc est l'épaisseur en plus de l'épaisseur de coupure
!                     ux/dx=1/(deltat) ou deltat est le temps nécessaire
!                     à la glace pour passer d'un noeud à l'autre
!                     la deuxieme partie du test revient à dh/dt*deltat > hmhc  

!                     Rajout vince : si le point amont est pose -> test vrai.

! rappel des differents ifrange
!-------------------------------
!Attention ifrange 1,2  ont sans doute un bug (il faudrait abs(uxbar))

! ifrange=1               Rajout vince : si un point voisin est pose -> test vrai.
!                         comme dans article fenno2007

! ifrange=2               pour les points ayant des voisins poses, garde le point 
!                         seulement si le bilan  surface fond est positif.
!                         bm(i,j)-bmelt(i,j) > 0.

! ifrange=3               pas de test sur l'epaisseur du point amont.
!                         test dhdt > -hmhc/deltat= -hmhc abs(U)/dx (correction du bug)
!                         + test sur bm(i,j)-bmelt(i,j) > 0.

! ifrange = 4             idem 3 mais avec test sur l'epaisseur du point amont

! ifrange = 0             pas de traitement specifique près du continent

! ifrange = -1            ancienne version MIS11

! ifrange = 5             semble idem 0

! ifrange = 6             ??

! ifrange = 7             calving drastique sauf si coefbmshelf < 0.5 

!cdc 
! ifrange=44              idem 4 mais avec test sur i et j couple // afq: probably useless
! ifrange=8               C. Dumas : coupure si alimentation ne permet pas H > Hcoup

if (ifrange.eq.1) then           !   Rajout vince : si un point voisin est pose -> test vrai.
                                 !   comme dans article fenno2007

   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
        .and.(dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx)))& 
        .or.(.not.flot(i-1,j)) ) !

   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) &
        .or.(.not.flot(i+1,j)) ) !

   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
        .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))&
        .or.(.not.flot(i,j-1)) ) !

   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))&
        .or.(.not.flot(i,j+1)) ) !

else if (ifrange.eq.2) then       ! pour les points ayant des voisins posés, seulement si le bilan 
                                  ! surface fond est positif. bm(i,j)-bmelt(i,j) > 0.

   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)

   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
        .and.  (dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx))) & 
        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !

   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx))) &
        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !

   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
        .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))&
        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !

   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))&
        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !

else if (ifrange.eq.3) then       ! nouvelle formulation Cat mars 08

! pas de test sur l'epaisseur du point amont.
! test dhdt > -hmhc/deltat= -hmhc abs(U)/dx

! pour les points ayant des voisins posés, seulement si le bilan 
! surface fond est positif. bm(i,j)-bmelt(i,j) > 0.


   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)

   testmij=( (uxbar(i,j).ge.0.) &  ! voisin (i-1,j) amont et > hcoup
        .and.(dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j))/dx))) & 
        .or.((.not.flot(i-1,j)).and.bilan_surf_fond ) !

   testpij=( (uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx))) &
        .or.((.not.flot(i+1,j)).and.bilan_surf_fond )  !

   testimj=(( uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx))) &
        .or.((.not.flot(i,j-1)).and.bilan_surf_fond )  !

   testipj=((uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx))) &
        .or.((.not.flot(i,j+1)).and.bilan_surf_fond )  !

else if (ifrange.eq.4) then       ! idem 3 mais avec test sur l'épaisseur du point amont

   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)

   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
        .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & 
        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !

   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) &
        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !

   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))&
        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !

   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))&
        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !

else if (ifrange.eq.0) then           ! pas de traitement special pres du continent

   testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
        .and.(dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))

   testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))

   testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))

   testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))

else if (ifrange.eq.-1) then           ! ancienne version MIS11

testmij=((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
   .and.(dhdt(i-1,j).gt.(hmhc(i-1,j)*uxbar(i-1,j)/dx)))    

testpij=((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
   .and.(dhdt(i+1,j).gt.(hmhc(i+1,j)*uxbar(i+1,j)/dx)))

testimj=((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
   .and.(dhdt(i,j-1).gt.(hmhc(i,j-1)*uybar(i,j-1)/dx)))

testipj=((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
    .and.(dhdt(i,j+1).gt.(hmhc(i,j+1)*uybar(i,j+1)/dx)))


else if (ifrange.eq.5) then           ! pas de traitement special pres du continent

   testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
        .and.((bil_tot(i-1,j)+bil_tot(i,j)) &
        .gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))

   testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.((bil_tot(i+1,j)+bil_tot(i,j)) &
        .gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))

   testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
        .and.((bil_tot(i,j-1)+bil_tot(i,j)) &
        .gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))

   testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.((bil_tot(i,j+1)+bil_tot(i,j)) &
        .gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))

else if (ifrange.eq.6) then           ! pas de traitement special pres du continent

   testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
        .and.(2.*bil_tot(i,j) &
        .gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx))))

   testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.(2.*bil_tot(i,j) &
        .gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j))/dx)))

   testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
        .and.(2.*bil_tot(i,j) &
        .gt.(-hmhc(i,j-1)*abs(uybar(i,j-1))/dx)))

   testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.(2.*bil_tot(i,j) &
        .gt.(-hmhc(i,j+1)*abs(uybar(i,j+1))/dx)))

else if (ifrange.eq.7) then           ! drastique calving sauf si coefbmelt.lt.0.5

   testmij= ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
        .and.(coefbmshelf.lt.1.))
   testpij= ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.(coefbmshelf.lt.1.))
   testimj= ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
         .and.(coefbmshelf.lt.1.))
   testipj= ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.(coefbmshelf.lt.1.))

else if (ifrange.eq.44) then       ! idem 4 mais avec test sur i et j couple

   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)

   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.)  &  ! voisin (i-1,j) amont et > hcoup
        .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & 
        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !

   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.) & ! voisin (i+1,j) amont et > hcoup
        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) &
        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !

   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.)  &  ! voisin (i,j-1) amont et > hcoup
        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))&
        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !

   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.) & ! voisin (i,j+1) amont et > hcoup
        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))&
        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !

else if (ifrange.eq.8) then ! test C Dumas : si somme flux entrant permet H > Hcoup : vrai 
                   !   bilan_surf_fond=(bm(i,j)-bmelt(i,j).gt.0.)

   testmij=( ((hmhc(i-1,j).gt.0.).and.(uxbar(i,j).ge.0.))  &  ! voisin (i-1,j) amont et > hcoup
                        !        .and.  (dhdt(i-1,j).gt.(-hmhc(i-1,j)*abs(uxbar(i-1,j)/dx)))) & 
        .or.((.not.flot(i-1,j)).and.bilan_surf_fond )) !

   testpij=( ((hmhc(i+1,j).gt.0.).and.(uxbar(i+1,j).le.0.)) & ! voisin (i+1,j) amont et > hcoup
                        !        .and.(dhdt(i+1,j).gt.(-hmhc(i+1,j)*abs(uxbar(i+1,j)/dx)))) &
        .or.((.not.flot(i+1,j)).and.bilan_surf_fond) ) !

   testimj=( ((hmhc(i,j-1).gt.0.).and.(uybar(i,j).ge.0.))  &  ! voisin (i,j-1) amont et > hcoup
                        !        .and.(dhdt(i,j-1).gt.(-hmhc(i,j-1)*abs(uybar(i,j-1)/dx))))&
        .or.((.not.flot(i,j-1)).and.bilan_surf_fond ) ) !

   testipj=( ((hmhc(i,j+1).gt.0.).and.(uybar(i,j+1).le.0.)) & ! voisin (i,j+1) amont et > hcoup
                        !        .and.(dhdt(i,j+1).gt.(-hmhc(i,j+1)*abs(uybar(i,j+1)/dx))))&
        .or.((.not.flot(i,j+1)).and.bilan_surf_fond ) ) !

   testH=0.
   if (testmij) testH=testH + abs(uxbar(i,j))*H(i-1,j)
   if (testpij) testH=testH + abs(uxbar(i+1,j))*H(i+1,j)
   if (testimj) testH=testH + abs(uybar(i,j))*H(i,j-1)
   if (testipj) testH=testH + abs(uybar(i,j+1))*H(i,j+1)
   
   if (testH.gt.Hcoup(i,j)) then
      testmij=.true.
   else
      testmij=.false.
      testpij=.false.
      testimj=.false.
      testipj=.false.
   endif
                   
endif



! detection des polynies dans les shelfs. on regarde vers l'aval

! 1 des 3 voisins aval > hcoup du cote west (i-1)
iin2=max(1,i-2)
iin3=max(1,i-3)

avalw=((hmhc(i-1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) &
     .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i,j).le.0.)

!avalw=((hmhc(i-1,j).gt.0..and.flot(i-1,j)).or.(hmhc(iin2,j).gt.0..and.flot(iin2,j)) &
!     .or.(hmhc(iin3,j).gt.0..and.flot(iin3,j))).and.(uxbar(i,j).le.0.)   

! 1 des 3 voisins aval > hcoup du cote est (i+1)
iin2=min(i+2,nx)
iin3=min(i+3,nx)

avale=((hmhc(i+1,j).gt.0.).or.(hmhc(iin2,j).gt.0.) &    
     .or.(hmhc(iin3,j).gt.0.)).and.(uxbar(i+1,j).ge.0)                          

!avale=((hmhc(i+1,j).gt.0..and.flot(i+1,j)).or.(hmhc(iin2,j).gt.0..and.flot(iin2,j)) &    
!     .or.(hmhc(iin3,j).gt.0..and.flot(iin3,j))).and.(uxbar(i+1,j).ge.0)
                
! 1 des 3 voisins aval > hcoup du cote sud (j-1)
jin2=max(1,j-2)
jin3=max(1,j-3)

avals=((hmhc(i,j-1).gt.0.).or.(hmhc(i,jin2).gt.0.) &    
     .or.(hmhc(i,jin3).gt.0.)).and.(uybar(i,j).le.0.)  
                          
!avals=((hmhc(i,j-1).gt.0..and.flot(i,j-1)).or.(hmhc(i,jin2).gt.0..and.flot(i,jin2)) &    
!     .or.(hmhc(i,jin3).gt.0..and.flot(i,jin3))).and.(uybar(i,j).le.0.)
                
! 1 des 3 voisins aval > hcoup du cote nord (j-1)
jin2=min(j+2,ny)
jin3=min(j+3,ny)

avaln=((hmhc(i,j+1).gt.0.).or.(hmhc(i,jin2).gt.0.) &    
     .or.(hmhc(i,jin2).gt.0.)) .and.(uybar(i,j+1).ge.0.)                         

!avaln=((hmhc(i,j+1).gt.0..and.flot(i,j+1)).or.(hmhc(i,jin2).gt.0..and.flot(i,jin2)) &    
!     .or.(hmhc(i,jin2).gt.0..and.flot(i,jin2))) .and.(uybar(i,j+1).ge.0.)

interieur=(avalw.or.avale).and.(avals.or.avaln)



 if ((.not.(testmij.or.testpij.or.testimj.or.testipj))  &  ! pas suffisament alimente
                .and.(.not.interieur)) then                ! et pas interieur

!cdc pour ne pas faire de calving si H > 350m :   if ((.not.(testmij.or.testpij.or.testimj.or.testipj)).and.(H(i,j).LT.350)) then
! mab: if none of the neighbours doesn't provide ice for the considered point,
! mab: then the ice is cut at this point and given to calv(i,j)
                !if ((.not.(testmij.or.testpij.or.testimj.or.testipj))) then
! mab: new version so that grid points with h=1m (all points were calving takes place are afterwards
! mab: set to a height of 1m) are not taken into account,  provides a problem as soon as you start 
! mab: coupling GRISLI to the iceberg modell
                   !if(h(i,j).gt.1) then
                     calv(i,j)=-h(i,j) 
!cdc                     H(i,j)=1. 
!cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01))
                     H(i,j)=0.
                     S(i,j)=H(i,j)*(1.-ro/row) + sealevel
                     B(i,j)=S(i,j) - H(i,j)
          ! ATTENTION ne pas mettre ice=0 sinon degradation bilan d'eau (bm et bmelt non comptabilises dans ce cas)

                      I_did_something = I_did_something + 1
                   endif

                  end if ifint
              end if ifext
         end do
      end do
      if (I_did_something.eq.0) then
          !        print*,'stop MULTI_CALV_LOOP l=',l
          EXIT MULTI_CALV_LOOP
      else
          !        print*,'calving continue l',l,I_did_something
      endif
    enddo MULTI_CALV_LOOP
         
      
! on met en calving les points detectes iceberg :
       where (iceberg(:,:).and.(H(:,:).gt.0.))
             calv(:,:)=-h(:,:)
             ice(:,:)=0
             H(:,:)=0.
             S(:,:)=H(:,:)*(1.-ro/row) + sealevel
             B(:,:)=S(:,:) - H(:,:)
       endwhere   
      
       calv_dtt(:,:) = calv_dtt(:,:) + calv(:,:) ! somme du calving sur dtt  
       ! calv_dtt est remis a 0 dans steps_time_loop (tous les dtt)
       
    end subroutine calving
!------------------------------------------------------------------------------------------
  end module calving_frange