!> \file flottab2-0.7.f90
!! Module pour determiner les endroits ou la glace
!! flotte , les iles, et la position du front de glace
!<

!> \namespace flottab_mod
!! Determine les endroits ou la glace
!! flotte , les iles, et la position du front de glace
!! \author Vincent & Cat
!! \date 10 juillet 2005
!! @note Used module
!! @note   - use module3D_phy
!! @note   - use module_choix
!! @todo "nasty Island". If the bedrock is above sealevel force grounded (mk_init)
!<
module flottab_mod
  
  !$ USE OMP_LIB
USE module3D_phy
use module_choix



 IMPLICIT NONE

 real ::  surnet !< surnet hauteur de glace au dessus de la mer
 real ::  archim !< test de flottaison
! real, parameter :: Hmin=1.001   !< Hmin pour être considere comme point ice 

 integer:: itestf

 logical,dimension(nx,ny) ::  gz1mx,gz1my
 logical,dimension(nx,ny) ::  fl1mx,fl1my


 real,dimension(nx,ny) ::  uxs1   !< uxbar a l'entree de flottab
 real,dimension(nx,ny) ::  uys1   !< uybar a l'entree de flottab


 integer pmx,pmy !pm=plus-moins -1 ou 1 pour x et y 


! Variables pour la determination des differents shelfs/stream
! (representés comme des taches ou l'on resoud l'eq elliptique)
!________________________________________________________________
 integer,parameter               :: n_ta_max=2000!< nombre de tache max
 integer,dimension(nx,ny)        :: table_out    !< pour les numeros des taches
 integer,dimension(nx,ny)        :: tablebis     !< pour les numeros des taches
 integer,dimension(0:n_ta_max)     :: compt        !< contient les equivalence entre les taches
 integer,dimension(0:n_ta_max)     :: nb_pts_tache !< indique le nombre de points par tache
 logical,dimension(0:n_ta_max)     :: iceberg1D      !< T si iceberg, F si calotte posee
 
 logical,dimension(nx,ny)        :: mask_tache_ij !< masque de travail 
                                                  !< vrai pour toute la tache de i,j
 integer,dimension(2)            :: smax_coord    !< pour le maxloc des iles

!  Variables pour determiner le point le plus haut (surf)
!  d'une ile completement stream
!_________________________________________________________

! icetrim : T si ice stream, F si calotte posee(vertical shear)
 logical,dimension(n_ta_max) :: icetrim 

 integer ::  ii,jj
 integer ::  smax_i
 integer ::  smax_j
 real    ::  smax_
 integer ::  numtache
 integer ::  nb_pt
 real    ::  petit_H=0.001 ! pour test ice sur zone flottante
contains
! -----------------------------------------------------------------------------------
!> SUBROUTINE: flottab()
!! Cette routine determine les endroits ou la glace
!! flotte , les iles, et la position du front de glace
!! @note Il y a 4 sortes de zone
!! @note    - Pose
!  @note    - Grounding zone et streams    gzmx et gzmy
!  @note    - Iles             ilemx, ilemy
!  @note    - flottant  flot sur le noeud majeur, flotmx sur le noeud mineur
!>
  subroutine flottab()
    !
    !       				Vince 5 Jan 95
    !                                       Modifie 20 Jan 95
    !                                       Modifie le 30 Novembre 98
    !    Passage f90 + determination des fronts Vincent dec 2003
    !    nettoyage et nouvelle détermination de gzmx et gzmy Cat le 10 juillet 2005
    !    Re-nettoyage par Cat en aout 2006.
    !    Le calcul de gzmx et gzmy pour les points intérieurs passe dans la subroutine dragging
    !    pour les points cotiers, toujours fait dans flottab
    !
    !                                       -----------------
    !
    !      Cette routine determine les endroits ou la glace
    !      flotte , les iles, et la position du front de glace
    !
    !      Il y a 4 sortes de zone
    !      Pose
    !      Grounding zone et streams    gzmx et gzmy
    !      Iles             ilemx, ilemy
    !      flottant  flot sur le noeud majeur, flotmx sur le noeud mineur

    !   passage dans flottab tous les pas de temps dt  ! 
    !
    ! _________________________________________________________

    !~ print*,'debut flottab',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel
    !~ print*,'debut flottab',flot(132,183),ice(132,183)
!print*,'H(90,179) flottab 1',H(90,179),ice(90,179), flot(90,179)

    if (itracebug.eq.1)  call tracebug(' Entree dans routine flottab')

    SHELFY = .FALSE.


    ! cas particulier des runs paleo ou on impose un masque grounded

    !$OMP PARALLEL PRIVATE(archim,surnet)
    if (igrdline.eq.2) then
       !$OMP WORKSHARE
       where ( mk_init(:,:).eq.1)                 ! pose
          flot(:,:) = .False.
          H(:,:)=max(H(:,:),(10.+sealevel_2d(:,:)-Bsoc(:,:))*row/ro)  ! pour avoir archim=10
          S(:,:) = Bsoc(:,:) + H(:,:) 
       elsewhere
          flot(:,:) = .True.
       end where
       !$OMP END WORKSHARE
    end if


    ! 1-INITIALISATION
    ! ----------------
    ! initialisation des variables pour detecter les points qui se mettent
    ! a flotter entre 2 dtt

    appel_new_flot=.false.
    !$OMP DO
    do j=1,ny
       do i=1,nx
          new_flot_point(i,j)=.false.
          new_flotmx(i,j)=.false.
          new_flotmy(i,j)=.false.
       enddo
    enddo
    !$OMP END DO

    ! ICE(:,:)=(H(:,:).gt.1) ! ice=.true. si epaisseur > 1m

    !$OMP WORKSHARE
    ICE(:,:)=0
    front(:,:)=0
    frontfacex(:,:)=0
    frontfacey(:,:)=0
    isolx(:,:)=.false.
    isoly(:,:)=.false.
    cotemx(:,:)=.false.
    cotemy(:,:)=.false.
    boost=.false.
    iceberg(:,:)=.false.
    !$OMP END WORKSHARE

    ! fin de l'initialisation
    !_____________________________________________________________________

    ! 2-TESTE LES NOUVEAUX POINTS FLOTTANTS
    ! -------------------------------------

    !$OMP DO
    do j=1,ny
       do i=1,nx

          uxs1(i,j)=uxbar(i,j)
          uys1(i,j)=uybar(i,j)

          archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j)
          !      if ((i.eq.132).and.(j.eq.183)) print*,'archim=',archim


          arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.e-3)) then    ! le point flotte
             mk(i,j)=1


             ex_pose: if ((.not.FLOT(I,J)).and.(isynchro.eq.1)) then  !  il ne flottait pas avant
                FLOT(I,J)=.true.
                BOOST=.false.

                if (igrdline.eq.1) then   ! en cas de grounding line prescrite
                   flot(i,j)=.false.
                   H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10
                   new_flot_point(i,j)=.false.
                endif

             else                                       ! isynchro=0 ou il flottait déja

                if (.not.FLOT(I,J)) then               ! il ne flottait pas (isynchro=0)
                   new_flot_point(i,j)=.true.          ! signale un point qui ne flottait pas
                   ! au pas de temps precedent
                   flot(i,j)=.true.

                   if (igrdline.eq.1) then   ! en cas de grounding line prescrite
                      flot(i,j)=.false.
                      H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10
                      new_flot_point(i,j)=.false.
                   endif

                endif
             endif ex_pose


          else if ((H(i,j).ge.0.).and.(archim.GE.0.)) then    !    le point ne flotte pas et est englace
             mk(i,j)=0

             if(FLOT(I,J)) then         !  mais il flottait avant
                FLOT(I,J)=.false.
                BOOST=.false.
             endif
             !cdc correction topo pour suivre  variations sealevel
             !cdd           S(i,j)=Bsoc(i,j)+H(i,j)
             S(i,j)=Bsoc(i,j)+H(i,j)
             B(i,j)=Bsoc(i,j)

          else if  ((H(i,j).LE.0.).and.(archim.LT.0.)) then    !    terre deglace qui devient ocean 
!cdc             ice(i,j)=0
!cdc         H(i,j)=1.
!cdc  1m           H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01))
             flot(i,j)=.true.  !cdc points ocean sont flot meme sans glace
             H(i,j)=0.
             S(i,j)=sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
             B(i,j)=S(i,j)-H(i,j)
          endif arch

          !   Si la glace flotte -> PRUDENCE !!!
          !   S et B sont alors determines avec la condition de flottabilite

          if (flot(i,j)) then
             shelfy = .true.

             surnet=H(i,j)*(1.-ro/row)
             S(i,j)=surnet+sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
             B(i,j)=S(i,j)-H(i,j)
          end if

       end do
    end do
    !$OMP END DO
       
!!$ do i=1,nx
!!$    do j=1,ny
!!$       if (flot(i,j)) then
!!$          mk(i,j)=1 
!!$       else
!!$          mk(i,j)=0
!!$       endif
!!$    end do
!!$ end do



    !-----------------------------------------------------------------------
    !$OMP DO
    domain_x: do j=1,ny       
       do i=2,nx

          !    3_x    A- NOUVELLE DEFINITION DE FLOTMX, LES POINTS
          !            AYANT UN DES VOISINS FLOTTANTS SONT FLOTMX ET GZMX
          !         -------------------------------------------

          !            fl1 est l'ancienne valeur de flotmx
          gz1mx(i,j)=gzmx(i,j)         !     gz1 est l'ancienne valeur de gzmx
          fl1mx(i,j)=flotmx(i,j)       !     fl1 est l'ancienne valeur de flotmx

          flotmx(i,j)=flot(i,j).and.flot(i-1,j)

          ! test pour detecter les nouveaux flotmx entre 2 dtt :

          if (flotmx(i,j).and.(new_flot_point(i,j).or. &
               new_flot_point(i-1,j))) then
             appel_new_flot=.true.
             new_flotmx(i,j)=.true.
          endif


          !   premiere determination de gzmx
          !__________________________________________________________________________

          ! gzmx si un des deux voisins est flottant et l'autre posé
          !                                                                   i-1   i
          gzmx(i,j)=((flot(i,j).and..not.flot(i-1,j)) & !         F    P
               .or.(.not.flot(i,j).and.flot(i-1,j)))         !         P    F

          !                 A condition d'etre assez proche de la flottaison
          !                 sur le demi noeud condition archim < 100 m

          archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-(sealevel_2d(i,j)+sealevel_2d(i-1,j))*0.5+ro/row*Hmx(i,j)
          gzmx(i,j)=gzmx(i,j).and.(archim.le.100.) 
          cotemx(i,j)=gzmx(i,j)

       end do
    end do domain_x
    !$OMP END DO
    !if (itracebug.eq.1)  call tracebug('  routine flottab apres domain_x')

    !     3_y    B- NOUVELLE DEFINITION DE FLOTMY
    !         --------------------------------
    !$OMP DO
    domain_y: do j=2,ny       
       do i=1,nx

          gz1my(i,j)=gzmy(i,j)           !       gz1 est l'ancienne valeur de gzmy
          fl1my(i,j)=flotmy(i,j)         !       fl1 est l'ancienne valeur de flotmy

          flotmy(i,j)=flot(i,j).and.flot(i,j-1)

          ! test pour detecter les nouveaux flotmy entre 2 dtt :

          if (flotmy(i,j).and.(new_flot_point(i,j).or.     &
               new_flot_point(i,j-1))) then
             appel_new_flot=.true.
             new_flotmy(i,j)=.true.
          endif

          !   premiere determination de gzmy
          !__________________________________________________________________________

          ! gzmy si un des deux voisins est flottant et l'autre posé

          gzmy(i,j)=((flot(i,j).and..not.flot(i,j-1))          &
               .or.(.not.flot(i,j).and.flot(i,j-1)))                 

          !                 A condition d'etre assez proche de la flottaison
          !                 sur le demi noeud condition archim > 100 m

          archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-(sealevel_2d(i,j)+sealevel_2d(i,j-1))*0.5+ro/row*Hmy(i,j)
          gzmy(i,j)=gzmy(i,j).and.(archim.le.100.) 
          cotemy(i,j)=gzmy(i,j)

       end do
    end do domain_y
    !$OMP END DO


    !-------------------------------------------------------------------------------------
    ! attention : pour expériences Heino
    !             gzmy(i,j)=gzmy_heino(i,j)
    !             shelfy=.true.
    !             shelfy=.false.
    !____________________________________________________________________________________ 


!!$
!!$
!!$! 4- Condition sur les bords
!!$
!!$      
!!$       do i=2,nx
!!$         flotmx(i,1) = (flot(i,1).or.flot(i-1,1))
!!$         flotmy(i,1) = .false. 
!!$       end do
!!$       
!!$       do j=2,ny
!!$          flotmy(1,j) = (flot(1,j).or.flot(1,j-1))
!!$          flotmx(1,j) = .false. 
!!$       end do
!!$
!!$       flotmx(1,1) = .false.
!!$       flotmy(1,1) = .false.
!!$


    !      4- determination des iles
    !      -------------------------
    !$OMP WORKSHARE
    ilemx(:,:)=.false.
    ilemy(:,:)=.false.
    !$OMP END WORKSHARE

! afq -- 17/01/19: on supprime les iles.
!    !       selon x
!    !$OMP DO
!    ilesx:  do j=2,ny-1
!       do i=3,nx-2
!          !                       F   G   F   
!          !                           x
!          ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
!          if ((flot(i-1,j).and..not.flot(i,j).and.flot(i+1,j)).and. &
!               (sdx(i,j).LT.1.E-02)) then
!             ilemx(i,j)=.true.
!             ilemx(i+1,j)=.true.

!             !                      F  G   G   F
!             !                         x
!             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
!          else if ((flot(i-1,j).and..not.flot(i,j)                  &
!               .and..not.flot(i+1,j)).and.flot(i+2,j).and.           &
!               (sdx(i,j).LT.1.E-02.and.sdx(i+1,j).LT.1.E-02)) then
!             ilemx(i,j)=.true.
!             ilemx(i+1,j)=.true.
!             ilemx(i+2,j)=.true.

!             !                      F   G   G   F
!             !                              x
!             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
!          else if ((flot(i-2,j).and..not.flot(i-1,j)                &
!               .and..not.flot(i,j)).and.flot(i+1,j).and.             &
!               (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02)) then
!             ilemx(i-1,j)=.true.
!             ilemx(i,j)=.true.
!             ilemx(i+1,j)=.true.

!             !                      F   G   G   G    F
!             !                              x
!             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
!          else if ((i.lt.nx-2)                                      &
!               .and.(flot(i-2,j).and..not.flot(i-1,j)                &
!               .and..not.flot(i,j)).and..not.flot(i+1,j)             &
!               .and.flot(i+2,j).and.                                &
!               (sdx(i,j).LT.1.E-02.and.sdx(i-1,j).LT.1.E-02         &
!               .and.sdx(i+1,j).LT.1.E-02)) then
!             ilemx(i-1,j)=.true.
!             ilemx(i,j)=.true.
!             ilemx(i+1,j)=.true.
!             ilemx(i+2,j)=.true.

!          endif

!       end do
!    end do ilesx
!    !$OMP END DO

!    !       selon y
!    !$OMP DO
!    ilesy: do j=3,ny-2
!       do i=2,nx-1
!          !                       F   G   F   
!          !                           x
!          ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
!          if ((flot(i,j-1).and..not.flot(i,j).and.flot(i,j+1)).and. &
!               (sdy(i,j).LT.1.E-02)) then
!             ilemy(i,j)=.true.
!             ilemy(i,j+1)=.true.

!             !                      F  G   G   F
!             !                         x
!             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
!          else if ((flot(i,j-1).and..not.flot(i,j)                  &
!               .and..not.flot(i,j+1)).and.flot(i,j+2).and.           &
!               (sdy(i,j).LT.1.E-02.and.sdy(i,j+1).LT.1.E-02)) then
!             ilemy(i,j)=.true.
!             ilemy(i,j+1)=.true.
!             ilemy(i,j+2)=.true.

!             !                      F   G   G   F
!             !                              x
!             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
!          else if ((flot(i,j-2).and..not.flot(i,j-1)                 &
!               .and..not.flot(i,j)).and.flot(i,j+1).and.              &
!               (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02)) then
!             ilemy(i,j-1)=.true.
!             ilemy(i,j)=.true.
!             ilemy(i,j+1)=.true.

!             !                      F   G   G   G    F
!             !                              x
!             ! modif tof 26/08/02 limite sur la pente (si diff S > 400 m)
!          else if ((j.lt.ny-2)                                      &
!               .and.(flot(i,j-2).and..not.flot(i,j-1)                &
!               .and..not.flot(i,j)).and..not.flot(i,j+1)             &
!               .and.flot(i,j+2).and.                                &
!               (sdy(i,j).LT.1.E-02.and.sdy(i,j-1).LT.1.E-02         &
!               .and.sdy(i,j+1).LT.1.E-02)) then 
!             ilemy(i,j-1)=.true.
!             ilemy(i,j)=.true.
!             ilemy(i,j+1)=.true.
!             ilemy(i,j+2)=.true.
!          endif
!       end do
!    end do ilesy
!    !$OMP END DO
!    ! fin des iles

    !$OMP END PARALLEL


    ! 5- calcule les noeuds qui sont streams a l'interieur et donne le betamx et betamy
    !----------------------------------------------------------------------------------
    call mstream_dragging          
    
    if (itracebug.eq.1)  call tracebug(' routine flottab apres call dragging')
!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
!!$if (itestf.gt.0) then
!!$   write(6,*) 'dans flottab apres dragging asymetrie sur H  pour time=',time
!!$   stop
!!$else
!!$   write(6,*) 'dans flottab aapres dragging pas d asymetrie sur H  pour time=',time
!!$
!!$end if

    !   6- calcule les vitesses des points qui sont devenus gzm
    !$OMP PARALLEL
    !$OMP DO         
    do j=1,ny
       do i=2,nx-1
          !      si le point etait posé (non gz) et devient gzmx
          !      definir la direction de la vitesse (moyenne des points)

          if ((.not.gz1mx(i,j)).and.(.not.fl1mx(i,j)).and.gzmx(i,j).and. &
               (i.gt.2).and.(i.lt.nx))  then
             uxs1(i,j)=(uxbar(i+1,j)+uxbar(i-1,j))/2.
          endif

       end do
    end do
    !$OMP END DO

    !$OMP DO 
    do j=2,ny-1
       do i=1,nx
          !      si le point etait posé (non gz) et devient gzmy
          !      definir la direction de la vitesse (moyenne des points)
          if ((.not.gz1my(i,j)).and.(.not.fl1my(i,j)).and.gzmy(i,j).and. &
               (j.gt.2).and.(j.lt.ny))  then
             uys1(i,j)=(uybar(i,j+1)+uybar(i,j-1))/2.
          endif

       end do
    end do
    !$OMP END DO


    !  7-On determine finalement la position des noeuds stream ou shelf
    !  -------------------------------------------------------------

    if (nt.ge.2) then   ! pour ne pas faire ce calcul lors du premier passage
       !$OMP WORKSHARE
       uxbar(:,:)=uxs1(:,:)
       uybar(:,:)=uys1(:,:)
       !$OMP END WORKSHARE
    endif

    !$OMP WORKSHARE
    flgzmx(:,:)=(marine.and.(flotmx(:,:).or.gzmx(:,:).or.ilemx(:,:)))   &
         .or.(.not.marine.and.flotmx(:,:))
    where (hmx(:,:).eq.0.)
      flgzmx(:,:) = .false.
    endwhere
    flgzmy(:,:)=(marine.and.(flotmy(:,:).or.gzmy(:,:).or.ilemy(:,:)))   &
         .or.(.not.marine.and.flotmy(:,:))
    where (hmy(:,:).eq.0.)
      flgzmy(:,:) = .false.
    endwhere
    !$OMP END WORKSHARE


    ! 8- Pour la fusion basale sous les ice shelves- region proche de la grounding line 
    !---------------------------------------------------------------------------------
    !       fbm est vrai si le point est flottant mais un des voisins est pose 
    !_________________________________________________________________________
    !$OMP DO
    do j=2,ny-1
       do i=2,nx-1

          ! if (i.gt.2.AND.i.lt.nx) then   
          fbm(i,j)=flot(i,j).and.      & 
               ((.not.flot(i+1,j)).or.(.not.flot(i,j+1)) &
               .or.(.not.flot(i-1,j)).or.(.not.flot(i,j-1)))
          !     endif
       end do
    end do
    !$OMP END DO


    !  9-On determine maintenant la position du front de glace
    !  -------------------------------------------------------
    !  C'est ici que l'on determine la position et le type de front
    !       print*,'on est dans flottab pour definir les fronts'



!!$do i=3,nx-2
!!$   do j=3,ny-2
!!$      if (h(i,j).gt.1.1) ice(i,j)=1        
!!$   end do
!!$end do
!			print*, 'flolottab debug', H(71,25),flot(71,25),bm(71,25),ice(71,25),time
!print*,'H(90,179) flottab 2',H(90,179),ice(90,179),flot(90,179)
    !$OMP WORKSHARE
    where (flot(:,:))
!cdc 1m       where (H(:,:).gt.max(Hmin,Hmin+BM(:,:)-Bmelt(:,:)))
      where (H(:,:).gt.max(BM(:,:)-Bmelt(:,:)+petit_H,0.)*dt)
!cdc			 where (H(:,:).gt.0.)
          ice(:,:)=1
       elsewhere 
          ice(:,:)=0
          H(:,:)=0.
          S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:)
          B(:,:)=S(:,:) - H(:,:)
       end where
    elsewhere
       where (H(:,:).gt.0.) 
          ice(:,:)=1
       elsewhere 
          ice(:,:)=0
       end where
    end where
    !$OMP END WORKSHARE
    !$OMP END PARALLEL
!    print*,'flottab',time,H(191,81),bm(191,81),bmelt(191,81),(BM(191,81)-Bmelt(191,81)+petit_H)*dt,ice(191,81)
!print*,'H(90,179) flottab 3',H(90,179),ice(90,179),flot(90,179),Bsoc(90,179)+H(90,179)*ro/row -sealevel
 
    !call determin_front ! cette version ne conserve pas la masse !!!
    call determin_front_tof ! version simplifiee

    ! pour sorties initMIP:
    debug_3D(:,:,118) = ice(:,:)*(1-mk(:,:))
    debug_3D(:,:,119) = ice(:,:)*mk(:,:)

    
  end subroutine flottab
!--------------------------------------------------------------------  

!> SUBROUTINE: determin_tache
!! Routine pour la dtermination du numero de tache a effectuer 
!>
subroutine determin_tache

!$ USE OMP_LIB

implicit none
integer :: indice
integer :: label         ! no des taches rencontrées dans le mask
integer :: label_max     ! no temporaire maxi de tache rencontrées
!      integer :: mask_nb = 4
integer,parameter :: mask_nb = 2   ! version ou on ne compte pas les diagonales
!     integer,dimension(mask_nb) :: mask
integer,dimension(mask_nb) :: mask


! 1-initialisation
!-----------------
label_max=1      ! numero de la tache, la premiere tache est notée 1
label=1
do i=1,n_ta_max
   compt(i)=i
enddo
!      table_in  = .false.
!$OMP PARALLEL
!$OMP WORKSHARE
table_out(:,:) = 0
iceberg1D(:)  = .true.
icetrim (:) = .true.
nb_pts_tache(:) = 0
!$OMP END WORKSHARE
!$OMP END PARALLEL
!    open(unit=100,file="tache.data",status='replace')

! 2-reperage des taches
!----------------------

do j=2,ny-1
  do i=2,nx-1

     IF (ice(i,j).ge.1) THEN ! on est sur la glace-----------------------------!
			
        if ((ice(i-1,j).ge.1).or.(ice(i,j-1).ge.1)) then   !masque de 2 cases adjacentes
           !      un des voisins est deja en glace
           mask(1) = table_out(i-1,j)
           mask(2) = table_out(i,j-1) 
           label = label_max

           !on determine la valeur de la tache minimun (>0) presente ds le masque
           do indice=1,mask_nb
              if (mask(indice).gt.0) label=min(label,mask(indice))
           enddo
!cdc       label=min(label,minval(mask(:), mask=mask > 0))

           !on fixe la valeur de la tache voisine minimun au point etudie (via label)
           table_out(i,j)=label
           !si ce noeud est posé, alors la tache n'est pas un iceberg et iceberg=.F.
           if (.not.FLOT(I,J)) then
              iceberg1D(label)=.false.
           endif

           !si ce noeud est posé, alors la tache n'est pas un ice stream et icestrim=.F.
           if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and.  &
               (.not.gzmy(i,j).and..not.gzmy(i,j+1))) then
              icetrim(label)=.false.
           endif

           ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt
           ! on lui affecte le numero de la tache fondamentale

           do indice=1,mask_nb
              if(mask(indice).gt.label) then
                 compt(mask(indice))=-label
              endif
           enddo
                                                                   ! exemple on est sur le point X : 5  X
          do indice=1,mask_nb                                      !                                   20
            if(mask(indice).gt.label) then                         ! mask(2)=20 > 5             
              compt(mask(indice))=label                            ! compt(20)=5
              if (.not.iceberg1D(mask(indice))) iceberg1D(label)=.false. ! si la tache n'etais pas un iceberg => iceberg =.false.  iceberg(5)=.false.               
              if (.not.icetrim(mask(indice))) icetrim(label)=.false.
              where (table_out(:,:).eq.mask(indice))               ! where table_out(:,:)=mask(2)=20
                table_out(:,:)=label                        ! table_out(:,:)=label=5                 
              endwhere
            endif
          enddo   

        else !aucun des voisins est une tache
           table_out(i,j)= label_max
           compt(label_max)=label_max
           !si ce noeud est posé, alors la ache n'est pas un iceberg et iceberg=.F.
           if (.not.FLOT(I,J)) then
              iceberg1D(label_max)=.false.
           endif

           !si ce noeud est posé, alors le tache n'est pas un ice stream et icestrim=.F.
           if ((.not.gzmx(i,j).and..not.gzmx(i+1,j)).and.  &
               (.not.gzmy(i,j).and..not.gzmy(i,j+1))) then
              icetrim(label)=.false.
           endif
             
           label_max  = label_max+1
           if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches icebergs=',label_max 
        endif


     else           !on est pas sur une tache----------------------------------------------
        table_out(i,j)=0    ! Pas necessaire (reecrit 0 sur 0)             
     endif          !---------------------------------------------------------------------


  enddo
enddo



! On reorganise compt en ecrivant le numero de la tache fondamentale
! i.e. du plus petit numero present sur la tache (Sans utiliser de recursivité)
! On indique aussi le nb de point que contient chaque taches (nb_pts_tache)

  !$OMP PARALLEL
  !$OMP DO
  do j=1,ny
    do i=1,nx
      if (table_out(i,j).ne.0) then
        nb_pts_tache(compt(table_out(i,j)))= nb_pts_tache(compt(table_out(i,j)))+1
      endif
    enddo
  enddo
  !$OMP END DO
  !$OMP END PARALLEL

  !On compte comme englacé uniquement les calottes dont une partie est posée      
  !$OMP PARALLEL PRIVATE(smax_,smax_coord,smax_i,smax_j,mask_tache_ij)
  !$OMP DO
  do j=3,ny-2
    do i=3,nx-2
      test1: if (.not.iceberg1D(table_out(i,j))) then ! on est pas sur un iceberg 
        if (nb_pts_tache(table_out(i,j)).ge.1) then
          ice(i,j)=1
          ! ici on est sur une tache non iceberg >= 5 points                       
          ! on teste si la tache n'est pas completement ice stream
    
          test2: if (icetrim(table_out(i,j))) then   ! on a une ile d'ice stream
    
            mask_tache_ij(:,:)=.false.
            mask_tache_ij(:,:)=(table_out(:,:).eq.table_out(i,j))       ! pour toute la tache
    
            smax_=maxval(S(:,:),MASK=mask_tache_ij(:,:))
            smax_coord(:)=maxloc(S(:,:),MASK=mask_tache_ij(:,:))
            smax_i=smax_coord(1)
            smax_j=smax_coord(2)
  
            gzmx(smax_i,smax_j)=.false. ; gzmx(smax_i+1,smax_j)=.false.
            gzmy(smax_i,smax_j)=.false. ; gzmx(smax_i,smax_j+1)=.false.
            flgzmx(smax_i,smax_j)=.false. ; flgzmx(smax_i+1,smax_j)=.false.
            flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false.
     
            if (Smax_.le.sealevel) then  ! afq -- WARNING: en fait avec le niveau marin local cest plus complique que ca!
            !if (Smax_.le.sealevel_2d(i,j)) then ! afq --WARNING: il faudrait regarder point par point sur la tache...
              write(num_tracebug,*)'Attention, une ile avec la surface sous l eau'
              write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j
            end if    
          endif test2 
        end if                                             ! endif deplace
!cdc transfere dans calving :             
      else ! on est sur un iceberg                          !   test1
        iceberg(i,j)=iceberg1D(table_out(i,j))
!~        ice(i,j)=0
!~        h(i,j)=0. !1. afq, we should put everything in calving!
!~        surnet=H(i,j)*(1.-ro/row)
!~        S(i,j)=surnet+sealevel
!~        B(i,j)=S(i,j)-H(i,j)             
      endif test1    
    end do
  end do    
  !$OMP END DO
  !$OMP END PARALLEL

debug_3D(:,:,124)=real(table_out(:,:))

end subroutine determin_tache
!----------------------------------------------------------------------
!> SUBROUTINE: determin_front
!!Routine pour la determination du front
!>
subroutine determin_front
!!$ USE OMP_LIB
integer :: i_moins1,i_plus1,i_plus2
integer :: j_moins1,j_plus1,j_plus2
      
      !$OMP PARALLEL
      !$OMP DO
      do j=3,ny-2
        do i=3,nx-2

 surice:if  (ice(i,j).eq.0) then
         do pmx=-1,1,2
          do pmy=-1,1,2

 diagice :  if (ice(i+pmx,j+pmy).eq.1) then

              if ((ice(i+pmx,j)+ice(i,j+pmy).eq.2)) then    ! test (i) pour eviter les langues
                                                            ! de glaces diagonales en coin(26dec04)
                   if ((ice(i+2*pmx,j).eq.1.and.ice(i+2*pmx,j+pmy).eq.0).or.&
                       (ice(i,j+2*pmy).eq.1.and.ice(i+pmx,j+2*pmy).eq.0)) then
                           ice(i,j)=1
                           h(i,j)=max(1.,h(i,j))
                   endif         

! test (i) pour eviter les langues de glaces diagonales :
!                        mouvement du cheval aux echecs  

               if ((ice(i+2*pmx,j+pmy)+ice(i+pmx,j+2*pmy).eq.1)) then
                   if (ice(i+2*pmx,j+pmy).eq.1.and. &
                      (ice(i+2*pmx,j+2*pmy)+ice(i,j+2*pmy)).ge.1) then
                           ice(i,j)=1
                           h(i,j)=max(1.,h(i,j))
                   endif        
                   if (ice(i+pmx,j+2*pmy).eq.1.and. &
                      (ice(i+2*pmx,j+2*pmy)+ice(i+2*pmx,j)).ge.1) then
                           ice(i,j)=1
                           h(i,j)=max(1.,h(i,j))
                   endif        

! test (ii) pour eviter les langues de glaces diagonales :
!            le point glace ice(i+pmx,j+pmy) a : 
!                          - ses 4 voisins frontaux en glace
!                          - mais 2 voisins vides diagonalement opposes   
            
               elseif ((ice(i+2*pmx,j+pmy)+ice(i+pmx,j+2*pmy).eq.2)  &
                                .and.ice(i+2*pmx,j+2*pmy).eq.0) then

! test (iii) pour faire les tests (i) et (ii)
                    ice(i,j)=1
                    h(i,j)=max(1.,h(i,j))
                    ice(i+2*pmx,j+2*pmy)=1
                    h(i+2*pmx,j+2*pmy)=max(1.,h(i+2*pmx,j+2*pmy))
               endif
            endif
           endif diagice   
          enddo
         enddo
         endif surice       
       end do
      end do
      !$OMP END DO
      !$OMP ENd PARALLEL

!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
!!$if (itestf.gt.0) then
!!$   write(6,*) 'dans front avant remplissage baies asymetrie sur H  pour time=',time
!!$   stop
!!$else
!!$   write(6,*) 'dans front avant remplissage baies pas d asymetrie sur H  pour time=',time
!!$
!!$end if


!     print*,'dans remplissage baies',time
       
baies: do k=1,2
!$OMP PARALLEL
!$OMP DO PRIVATE(i_moins1,j_moins1,i_plus1,j_plus1,i_plus2,j_plus2)
do j=1,ny
   do i=1,nx

surice_xy: if  (ice(i,j).eq.0) then
   i_moins1=max(i-1,1)
   j_moins1=max(j-1,1)
   i_plus1=min(i+1,nx)
   j_plus1=min(j+1,ny)
   i_plus2=min(i+2,nx)
   j_plus2=min(j+2,ny)
  
! test (iii) pour trouver les baies de largeur 1 ou 2 cases 
! et combler les trous si ce sont des baies
! si ce ne sont pas des baies, ne pas combler et creer des langues de glaces artificielles             
! baies horizontales  
 
         if (ice(i_moins1,j).eq.1.and.(ice(i_plus1,j).eq.1.or.ice(i_plus2,j).eq.1)) then 
            if (ice(i,j_moins1).eq.1.or.ice(i,j_plus1).eq.1)  then    ! ice(i,j)=1
               ice(i,j)=1
               H(i,j)=max(1.,H(i,j))
            endif
         endif


         if (ice(i,j_moins1).eq.1.and.(ice(i,j_plus1).eq.1.or.ice(i,j_plus2).eq.1)) then 
            if (ice(i_moins1,j).eq.1.or.ice(i_plus1,j).eq.1) then !ice(i,j)=1
               ice(i,j)=1
               H(i,j)=max(1.,H(i,j))
            endif
         endif

      endif surice_xy
   end do
end do
!$OMP END DO
!$OMP END PARALLEL
end do baies

!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
!!$if (itestf.gt.0) then
!!$   write(6,*) 'dans front apres remplissage baies asymetrie sur H  pour time=',time
!!$   stop
!!$else
!!$   write(6,*) 'dans front apres remplissage baies pas d asymetrie sur H  pour time=',time
!!$
!!$end if

!$OMP PARALLEL
!$OMP DO
do j=2,ny-1
   do i=2,nx-1

      if (ice(i,j).eq.1) then     !         test si ice=1

! if ice, on determine front...
! ainsi, front=0 sur les zones = 0

         front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
         !front= le nb de faces en contact avec un voisin englacé
      endif
   end do
end do
!$OMP END DO

! traitement des bords. On considere que l'exterieur n'a pas de glace
! attention ce n'est vrai que sur la grande grille

!$OMP DO PRIVATE(i)
do j=2,ny-1
   i=1
   front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
   i=nx
   front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
end do
!$OMP END DO 

!$OMP DO PRIVATE(j)
do i=2,nx-1
   j=1 
   front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))
   j=ny
   front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
end do
!$OMP END DO 

! traitement des coins

front(1,1)=ice(2,1)+ice(2,1)
front(1,ny)=ice(2,ny)+ice(1,ny-1)
front(nx,1)=ice(nx,2)+ice(nx-1,1)
front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)

!!$call detect_assym(nx,ny,0,41,1,0,1,0,H,itestf)
!!$if (itestf.gt.0) then
!!$   write(6,*) 'dans front apres front asymetrie sur H  pour time=',time
!!$   stop
!!$else
!!$   write(6,*) 'dans front apres front pas d asymetrie sur H  pour time=',time
!!$
!!$end if

!   on ne compte pas les taches de glace de 2 cases (horizontales ou verticales)
!   en fait, si ces deux cases sont flottantes, il faut enlever les icebergs 
!   de n'importe quelle taille).
!   si ces deux taches sont posées (ou une des deux), il n'y a pas assez de conditions aux limites

!$OMP DO 
do j=1,ny
   do i=1,nx-1
      if (front(i,j).eq.1) then
         if (front(i+1,j).eq.1) then
            ice(i,j)=0
            ice(i+1,j)=0
            front(i,j)=0
            front(i+1,j)=0
         endif
      endif
   end do
end do
!$OMP END DO

!$OMP DO 
do j=1,ny-1
   do i=1,nx
      if (front(i,j).eq.1) then
         if (front(i,j+1).eq.1) then
            ice(i,j)=0
            ice(i,j+1)=0
            front(i,j)=0
            front(i,j+1)=0
         endif
      end if
   end do
end do
!$OMP END DO

!isolx signifie pas de voisins en x
!isoly signifie pas de voisins en y
!remarque : 
!si isolx/y=.true. alors frontfacex/y=0 (a la fois +1 & -1 or +1-1=0)

! calcul de frontfacex et isolx
!$OMP DO
do j=1,ny
   do i=2,nx-1

      if (front(i,j).ge.1.and.front(i,j).le.3) then    !front(entre 1 et 3)
         
         if ((ice(i-1,j)+ice(i+1,j)).lt.2) then        ! il y a un front // a x
          
            if ((ice(i-1,j)+ice(i+1,j)).eq.0) then
               isolx(i,j)=.true.
            elseif (ice(i-1,j).eq.0) then
               frontfacex(i,j)=-1                      ! front  i-1 |i  i+1
            else
               frontfacex(i,j)=+1                      ! front  i-1  i| i+1
            endif
         endif
      end if !fin du test il y a un front
         
   end do
end do
!$OMP END DO

! calcul de frontfacey et isoly
!$OMP DO
do j=2,ny-1
   do i=1,nx

      if (front(i,j).ge.1.and.front(i,j).le.3) then   !front(entre 1 et 3)
         
         if ((ice(i,j-1)+ice(i,j+1)).lt.2) then       ! il y a un front // a y
             
            if ((ice(i,j-1)+ice(i,j+1)).eq.0) then
               isoly(i,j)=.true.                      !front   j-1 |j| j+1
            elseif (ice(i,j-1).eq.0) then
               frontfacey(i,j)=-1                     !front   j-1 |j j+1
            else
               frontfacey(i,j)=+1                     !front   j-1  j| j+1
            endif
         endif
     end if !fin du test il y a un front
         
   end do
end do
!$OMP END DO


! traitement des bords. On considere que l'exterieur n'a pas de glace
! attention ce n'est vrai que sur la grande grille

!$OMP DO PRIVATE(i)
do j=2,ny-1
   i=1
   if (front(i,j).ge.1)  then 
      if (ice(i+1,j).eq.0) then 
         isolx(i,j)=.true.
      else
         frontfacex(i,j)=-1
      endif
   end if
   i=nx
   if (front(i,j).ge.1)  then 
      if (ice(i-1,j).eq.0) then 
         isolx(i,j)=.true.
      else
         frontfacex(i,j)=1
      endif
   end if
end do
!$OMP END DO

!$OMP DO PRIVATE(j)
do i=2,nx-1
   j=1 
   if (front(i,j).ge.1)  then 
      if (ice(i,j+1).eq.0) then 
         isoly(i,j)=.true.
      else
         frontfacey(i,j)=-1
      endif
   end if
   j=ny
   if (front(i,j).ge.1)  then 
      if (ice(i,j-1).eq.0) then 
         isoly(i,j)=.true.
      else
         frontfacey(i,j)=1
      endif
   end if
end do
!$OMP END DO
!$OMP END PARALLEL

return
end subroutine determin_front
!------------------------------------------------------------------------------

subroutine determin_front_tof

integer,dimension(nx,ny) :: nofront ! tableau de travail (points au dela du front)
!$OMP PARALLEL
!$OMP DO
do j=2,ny-1
   do i=2,nx-1

      if (ice(i,j).eq.1) then     !         test si ice=1

! if ice, on determine front...
! ainsi, front=0 sur les zones = 0

         front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
         !front= le nb de faces en contact avec un voisin englacé
      endif
   end do
end do
!$OMP END DO

! traitement des bords. On considere que l'exterieur n'a pas de glace
! attention ce n'est vrai que sur la grande grille

!$OMP DO PRIVATE(i)
do j=2,ny-1
   i=1
   front(i,j)=(ice(i+1,j)+ice(i,j+1)+ice(i,j-1))
   i=nx
   front(i,j)=(ice(i-1,j)+ice(i,j+1)+ice(i,j-1))
end do
!$OMP END DO 

!$OMP DO PRIVATE(j)
do i=2,nx-1
   j=1 
   front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j+1))
   j=ny
   front(i,j)=(ice(i-1,j)+ice(i+1,j)+ice(i,j-1))
end do
!$OMP END DO 
!$OMP BARRIER
! traitement des coins

front(1,1)=ice(2,1)+ice(2,1)
front(1,ny)=ice(2,ny)+ice(1,ny-1)
front(nx,1)=ice(nx,2)+ice(nx-1,1)
front(nx,ny)=ice(nx,ny-1)+ice(nx-1,ny)

! Les points à plus d'un point de grille du bord sont front=-1
!$OMP WORKSHARE
nofront(:,:)=0
!$OMP END WORKSHARE
!$OMP BARRIER

!$OMP DO
do j=2,ny-1
	do i=2,nx-1
		if ((ice(i,j).eq.0).and.(front(i-1,j)+front(i+1,j)+front(i,j-1)+front(i,j+1)).eq.0) then
			nofront(i,j)=-1
		endif
	enddo
enddo
!$OMP END DO
!$OMP BARRIER

!$OMP WORKSHARE
where (nofront(:,:).eq.-1)
	front(:,:)=-1
endwhere
!$OMP END WORKSHARE
!$OMP END PARALLEL		

end subroutine determin_front_tof


!> SUBROUTINE: determin_marais
!! afq -- Routine pour l'identification des "marais"
!! Un marais est une tache "shelf" entouré de points grounded
!! Copie sauvage de determin_tache, adapte au probleme du marais
!>
subroutine determin_marais

  !$ USE OMP_LIB

  implicit none
  
  integer :: indice
  integer :: label         ! no des taches rencontrées dans le mask
  integer :: label_max     ! no temporaire maxi de tache rencontrées
  integer,parameter :: mask_nb = 2   ! version ou on ne compte pas les diagonales
  integer,dimension(mask_nb) :: mask ! numero de tache des points adjacents

  integer,dimension(nx,ny)          :: table_out_marais    !< numeros de tache d'un point ij
  integer,dimension(0:n_ta_max)     :: compt_marais        !< contient les equivalence entre les taches
  integer,dimension(0:n_ta_max)     :: nb_pts_marais       !< indique le nombre de points par tache
  logical,dimension(0:n_ta_max)     :: marais              !< T si flottants entoure de poses, F sinon


! 1-initialisation
!-----------------
  label_max=1      ! numero de la tache, la premiere tache est notée 1
  label=1
  do i=1,n_ta_max
    compt_marais(i)=i
  enddo
  !$OMP PARALLEL
  !$OMP WORKSHARE
  table_out_marais(:,:) = 0
  marais(:)  = .true.
  nb_pts_marais(:) = 0
  !$OMP END WORKSHARE
  !$OMP END PARALLEL

! 2-reperage des taches
!----------------------

  do j=2,ny-1
    do i=2,nx-1
      if ((ice(i,j).ge.1).and.flot(i,j)) then ! on est sur la glace qui flotte-------------------!

        if (((ice(i-1,j).ge.1).and.flot(i-1,j)).or.((ice(i,j-1).ge.1).and.flot(i,j-1))) then   !masque de 2 cases adjacentes
          !  un des voisins est deja en glace
          mask(1) = table_out_marais(i-1,j)
          mask(2) = table_out_marais(i,j-1) 
          label = label_max

          ! on determine la valeur de la tache minimun (>0) presente ds le masque
          do indice=1,mask_nb
            if (mask(indice).gt.0) label=min(label,mask(indice))
          enddo

          ! on fixe la valeur de la tache voisine minimun au point etudie (via label)
          table_out_marais(i,j)=label
           
          !si un des voisins n'est pas glace alors la tache n'est pas un marais
          if ( (ice(i+1,j).eq.0) .or. (ice(i,j+1).eq.0) .or. (ice(i-1,j).eq.0) .or. (ice(i,j-1).eq.0) ) then
            marais(label)=.false.
          endif

          ! si 2 taches differentes sont dans le masque, il faut les identifier dans compt_marais
          ! on lui affecte le numero de la tache fondamentale

                                                                   ! exemple on est sur le point X : 5  X
          do indice=1,mask_nb                                      !                                   20
            if(mask(indice).gt.label) then                         ! mask(2)=20 > 5             
              compt_marais(mask(indice))=label                     ! compt_marais(20)=5
              if (.not.marais(mask(indice))) marais(label)=.false. ! si la tache n'etais pas un marais => marais =.false.  marais(-(-5))=.false.               
              where (table_out_marais(:,:).eq.mask(indice))        ! where table_out_marais(:,:)=mask(2)=20
                table_out_marais(:,:)=label                        ! table_out_marais(:,:)=label=5                 
              endwhere                                                
            endif
          enddo
           
        else !aucun des voisins est une tache
          table_out_marais(i,j)= label_max
          compt_marais(label_max)=label_max
           
          ! si un des voisins n'est pas glace alors la tache n'est pas un marais
          if ( (ice(i+1,j).eq.0) .or. (ice(i,j+1).eq.0) .or. (ice(i-1,j).eq.0) .or. (ice(i,j-1).eq.0) ) then
            marais(label_max)=.false.
          endif
          label_max  = label_max+1
          if (label_max.gt.n_ta_max) print*,'ATTENTION trop de taches=',label_max 
        endif
      else          ! on est pas sur une tache--------------------------------------------
        table_out_marais(i,j)=0    ! Pas necessaire (reecrit 0 sur 0)             
      endif         !---------------------------------------------------------------------
    enddo
  enddo

  marais(0)=.false.

  !$OMP PARALLEL
  !$OMP DO
  do j=1,ny
    do i=1,nx
      if (table_out_marais(i,j).ne.0) then
        nb_pts_marais(compt_marais(table_out_marais(i,j)))= nb_pts_marais(compt_marais(table_out_marais(i,j)))+1   
      endif
      flot_marais(i,j) = marais(table_out_marais(i,j))
    enddo
  enddo
  !$OMP END DO
  !$OMP END PARALLEL

  debug_3D(:,:,125)=real(table_out_marais(:,:))

end subroutine determin_marais

end module flottab_mod