!> \file bilan_eau.f90
!! Clacul du bilan d'eau sur la calotte
!<
      

!> MODULE: bilan_eau
!! Calcule le bilan d'eau sur la grille GRISLI 
!! \author C. Dumas
!! \date 20/02/2017
!! @note Permet de vérifier que le bilan d'eau est ok et qu'il n'y a pas de fuite.
!>

! bilan d'eau :
! ce qui tombe : Acc
! ce qui sort : Bmelt, Abl, calving, ablbord
! ce qui est stocké : masse de glace
! le total de tous ces termes doit être constant (sur la grille, pas localement)

module bilan_eau_mod


use module3D_phy
implicit none

real :: sum_H


real,dimension(nx,ny) :: tot_water    !< bilan d'eau
real,dimension(nx,ny) :: H_beau_old, diff_H_2D, diff_H_water_bilan
real,dimension(nx,ny) :: Bm_dtt       !< mass balance on ice points accumulated during dtt
real,dimension(nx,ny) :: Bmelt_dtt    !< basal melting on ice points accumulated during dtt
real,dimension(nx,ny) :: calv_dtt     !< calving sur dtt (pour calcul bilan d'eau)
real,dimension(nx,ny) :: archimtab    ! point pose si > 0
real,dimension(nx,ny) :: grline_dtt   ! grounding line flux during dtt

real :: sum_H_old
real :: diff_H
real :: alpha_flot                      ! ro/row 

real,dimension(nx,ny) :: bm_dt,bmelt_dt


contains

subroutine init_bilan_eau
 ! initialisation des variables
        diff_H=0.
        sum_H_old = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)
        tot_water(:,:)=0.
        bm_dt(:,:)=0.
        bmelt_dt(:,:)=0.
        alpha_flot=ro/row
        archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:)
        gr_line(:,:)=0
        do j=1,ny
           do i=1,nx
              !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice
              if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice
                 if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1
                 if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) gr_line(i,j)=1
                 if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) gr_line(i,j)=1
                 if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) gr_line(i,j)=1
              endif
           enddo
        enddo

end subroutine init_bilan_eau




subroutine bilan_eau

tot_water(:,:)=0.

bm_dt(:,:)=0.
bmelt_dt(:,:)=0.

bm_dt(2:nx-1,2:ny-1)=bm(2:nx-1,2:ny-1)*dt
bmelt_dt(2:nx-1,2:ny-1)=bmelt(2:nx-1,2:ny-1)*dt
!~ where (Bm(:,:).lt.0.)
!~ 	bm_dt(:,:)=max(Bm(:,:)*dt,-H(:,:))
!~ elsewhere
!~ 	bm_dt(:,:)=bm(:,:)*dt
!~ endwhere

!~ where (Bmelt(:,:).lt.0.)
!~ 	bmelt_dt(:,:)=max(Bmelt(:,:)*dt,-H(:,:))
!~ elsewhere
!~ 	bmelt_dt(:,:)=bmelt(:,:)*dt
!~ endwhere


! on fait la somme sur la grille en excluant les bords :

sum_H = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)

diff_H = diff_H + (sum_H - sum_H_old) ! on calcul la variation de volume a chaque dt

diff_H_2D(2:nx-1,2:ny-1)=H(2:nx-1,2:ny-1)-H_beau_old(2:nx-1,2:ny-1)

where (ice(2:nx-1,2:ny-1).eq.1)
        Bm_dtt(2:nx-1,2:ny-1) = Bm_dtt(2:nx-1,2:ny-1) + Bm_dt(2:nx-1,2:ny-1) !* dt ! somme Bm sur dt
        bmelt_dtt(2:nx-1,2:ny-1) = bmelt_dtt(2:nx-1,2:ny-1) + bmelt_dt(2:nx-1,2:ny-1) ! * dt ! somme bmelt sur dt
endwhere

archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:)
gr_line(:,:)=0
do j=1,ny
  do i=1,nx
     !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice
     if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice
      if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1
      if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) gr_line(i,j)=1
      if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) gr_line(i,j)=1
      if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) gr_line(i,j)=1
    endif
  enddo
enddo

where (gr_line(:,:).eq.1)
!~   grline_dtt(:,:)= (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ &
!~                     (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5 &
!~                     *H(:,:) + grline_dtt(:,:) 
   grline_dtt(:,:)= - sqrt(                                                                        &
                           ( (uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))*dy/2. )**2 +  &
                           ( (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))*dx/2. )**2 )  &
                           * H(:,:) * dt + grline_dtt(:,:)
endwhere


if (isynchro.eq.1) then
! on raisonne en bilan annuel pour simplifier :
!~ 	where (ice(2:nx-1,2:ny-1).eq.1)
!~ 		tot_water(2:nx-1,2:ny-1) = (Bm_dtt(2:nx-1,2:ny-1) - Bmelt_dtt(2:nx-1,2:ny-1) + calv_dtt(2:nx-1,2:ny-1) - ablbord_dtt(2:nx-1,2:ny-1)) / dtt 
!~ 		elsewhere
!~ 		tot_water(2:nx-1,2:ny-1) = (calv_dtt(2:nx-1,2:ny-1) - ablbord_dtt(2:nx-1,2:ny-1)) / dtt
!~ 	endwhere
!cdc pas besoin du test sur ice ici, il a ete fait avant (et le masque ice varie a chaque dt)	
  tot_water(2:nx-1,2:ny-1) = (Bm_dtt(2:nx-1,2:ny-1) - Bmelt_dtt(2:nx-1,2:ny-1) + calv_dtt(2:nx-1,2:ny-1) - ablbord_dtt(2:nx-1,2:ny-1)) / dtt
   
! bilan d'eau sur la grille :
        water_bilan=sum(tot_water(:,:))
        diff_H = diff_H/dtt

!999 format(f0.2,1x,e15.8,1x,i10,8(1x,e15.8))
!	write(6,999),time,sum_H,count(ice(:,:)==1),diff_H,water_bilan,sum(calv_dtt(:,:))/dtt,sum(ablbord_dtt(:,:))/dtt,sum(bmelt_dtt(:,:),mask=ice(:,:)==1)/dtt,sum(bm(:,:),mask=ice(:,:)==1),sum(Bm_dtt(:,:))/dtt,sum(bmelt_dtt(:,:))/dtt
        diff_H_water_bilan(2:nx-1,2:ny-1)=tot_water(2:nx-1,2:ny-1)-diff_H_2D(2:nx-1,2:ny-1)

endif
sum_H_old=sum_H

H_beau_old(:,:)=H(:,:)

end subroutine bilan_eau

end module bilan_eau_mod