!> \file bmelt-seuil-profondeur_mod.f90
!! Prametrise la fusion basale (ice shelves)
!<

!> \namespace bmelt_seuil_prof
!! Prametrise la fusion basale (ice shelves)
!! \author ...
!! \date ...
!! @note Pour l'actuel : 2 valeurs  pour 2 domaines de profondeur
!! + une valeur pour bmgrz
!! A choisir dans le module_choix
!! @note Used module
!! @note   - use module3D_phy
!<

module  bmelt_seuil_prof         

  ! prametrise la fusion basale (ice shelves)
  ! Pour l'actuel : 2 valeurs  pour 2 domaines de profondeur
  ! + une valeur pour bmgrz
  ! A choisir dans le module_choix

  use module3d_phy


  implicit none

  real :: bm_grz                        !< valeur prescrite a la grounding zone
  real,dimension(nx,ny) ::  bmgrz       !< tabelau fusion basale a la grounding zone

  real :: bmshelf_plateau               !< valeur prescrite sur le plateau cont.
  real :: bmshelf_abysses               !< valeur prescrite au dessus de l'ocean profond
  real :: depth_talus                   !< profondeur de transition 
  real,dimension(nx,ny) ::  bmshelf     !< tableau fusion basale sous shelf


contains
  !-------------------------------------------------------------------------------
  !> SUBROUTINE: init_bmelt
  !! Cette routine fait l'initialisation pour la fusion basale.
  !! @note Elle est appelée par inputfile-vec-0.5.f90
  !!
  !>
  subroutine init_bmelt



    ! Cette routine fait l'initialisation pour la fusion basale.
    ! Elle est appelée par inputfile-vec-0.5.f90


    namelist/bmelt_seuil/bm_grz,bmshelf_plateau,bmshelf_abysses,depth_talus

    if (itracebug.eq.1)  call tracebug('entree dans init_bmelt de bmelt_seuil_prof')

    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
    read(num_param,bmelt_seuil)
!    write(num_rep_42,bmelt_seuil)

    !    ecriture dans le fichier parametres

    write(num_rep_42,'(A)')'!___________________________________________________________' 
    write(num_rep_42,'(A)') '&bmelt_seuil                     ! module  bmelt_seuil_prof'
    write(num_rep_42,*)   'bm_grz           =',bm_grz
    write(num_rep_42,*)   'bmshelf_plateau  =',bmshelf_plateau
    write(num_rep_42,*)   'bmshelf_abysses  =',bmshelf_abysses
    write(num_rep_42,*)   'depth_talus      =',depth_talus
    write(num_rep_42,*)'/'                      
    write(num_rep_42,'(A)') '! Pour l actuel : bm_grz a la grounding line'
    write(num_rep_42,'(A)') '!                 bmshelf_plateau sur le plateau continental'
    write(num_rep_42,'(A)') '!                 bmshelf_abysses pour les grandes profondeurs'
    write(num_rep_42,'(A)') '!       depth_talus, negative, separation entre les 2 domaines'
    write(num_rep_42,*)
    write(num_rep_42,'(A)')'!___________________________________________________________' 



    bmgrz(:,:) = bm_grz

    where (Bsoc0(:,:).lt.depth_talus)
        bmshelf(:,:)=bmshelf_abysses
    elsewhere
        bmshelf(:,:)=bmshelf_plateau
    end where

    debug_3D(:,:,34)=bmshelf(:,:)
    return
  end subroutine init_bmelt

  !________________________________________________________________________________

  !> SUBROUTINE: bmeltshelf
  !! Cette routine calcule la fusion basale proprement dite pour les points flottants
  !! @note coefbmshelf a ete calcule par forclim et permet de faire varier en fonction du climat
  !>

  subroutine bmeltshelf


    ! cette routine calcule la fusion basale proprement dite pour les points flottants
    ! coefbmshelf a ete calcule par forclim et permet de faire varier en fonction du climat

    integer :: ngr           ! nombre de voisins flottants
    real    :: coef_talus    ! pour ne pas changer la fusion au dessus de l'ocean profond

    if (itracebug.eq.1)  call tracebug('entree dans bmeltshelf de bmelt_seuil_prof')

    where (Bsoc0(:,:).lt.depth_talus)
        bmshelf(:,:)=bmshelf_abysses
    elsewhere
        bmshelf(:,:)=bmshelf_plateau
    end where

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

        talus_nochange: if (Bsoc0(i,j).lt.depth_talus) then 
            coef_talus = 1         ! pas de changement au dessus de l'ocecan profond
        else
            coef_talus = coefbmshelf
        endif talus_nochange


        shelf: if (flot(i,j)) then    ! partie flottante

            bmelt(i,j)=coef_talus*bmshelf(i,j)

            ! fbm est vrai si le point est flottant mais un des voisins est pose
            ! calcule dans flottab

            if (fbm(i,j)) then 
                bmelt(i,j)=coef_talus*bmgrz(i,j)
            endif



            ! ATTENTION le bloc suivant sert a determiner la fusion basale d'equilibre
            ! pour les shelves stationnaires

            if ((igrdline.eq.1).and.(ibmelt_inv.eq.1)) then
                !     corrbmelt(i,j)=hdot(i,j)+bmelt(i,j)   ! le bmelt d'equilibre
                !     bmelt(i,j)=corrbmelt(i,j)
                corrbmelt(i,j)=corrbmelt(i,j)+hdot(i,j)*0.85
                bmelt(i,j)=bmelt(i,j)+corrbmelt(i,j)     
            endif


        else                   ! point posé, on compte le nombre de voisins flottants
            ngr=0
            if ((i.ne.1).and.(i.ne.nx).and.(j.ne.1).and.(j.ne.ny)) then 
                if (flot(i+1,j)) ngr=ngr+1
                if (flot(i-1,j)) ngr=ngr+1
                if (flot(i,j+1)) ngr=ngr+1
                if (flot(i,j-1)) ngr=ngr+1
            end if


            !   la fusion des points limites est une combinaison entre valeur posée et valeur flottante
            !   en fonction du nombre de points flottants

            bmelt(i,j)= ngr/4.*bmgrz(i,j)*coef_talus+(1.-ngr/4.)*bmelt(i,j)

            if ((igrdline.eq.1).and.(ibmelt_inv.eq.1)) then
                corrbmelt(i,j)=corrbmelt(i,j)+hdot(i,j)  !*0.85
                bmelt(i,j)=bmelt(i,j)+corrbmelt(i,j)
            endif
        endif shelf
      end do
    end do

if ((igrdline.eq.1).and.(ibmelt_inv.eq.0)) then
    where (i_Hp(:,:).eq.1)
        bmelt(:,:)=0.    ! hdot donne alors le -bmelt
    endwhere
endif

return
end subroutine  bmeltshelf



end module  bmelt_seuil_prof