!> \file bilan_flux_output_mod.f90
!! Clacul du bilan des flux : sorties initMIP
!<
      

!> MODULE: bilan_flux
!! Calcule le bilan des flux : sortie initMIP
!! \author C. Dumas
!! \date 10/10/2017
!! @note Calcul des flux output initMIP
!>

module bilan_flux_mod


use module3D_phy
implicit none

real :: dtoutflux                       ! pas de temps des sorties flux
real :: dtsum                           ! somme des pas de temps (pour verification)
!real :: alpha_flot                      ! ro/row 

! variables flux :
real,dimension(nx,ny) :: bm_flux        ! surface smb on ice mask
real,dimension(nx,ny) :: bmeltgr_flux   ! bmelt on ice grounded mask
real,dimension(nx,ny) :: bmeltfl_flux   ! bmelt on floating ice mask
double precision,dimension(nx,ny) :: dhdt_flux      ! dhdt on ice mask
real,dimension(nx,ny) :: calving_flux   ! calving on ice mask
real,dimension(nx,ny) :: grline_flux    ! grounding line flux
real,dimension(nx,ny) :: calv_abl_bord_flux ! calving + abl_bord sur points en contact avec la mer

double precision,dimension(nx,ny) :: vieuxH_flux    ! H at previous time step
!real,dimension(nx,ny) :: archimtab      ! point pose si > 0
!integer,dimension(nx,ny) :: gr_line        ! point grounding line (pose sous niveau de la mer avec point flottant a cote)

real,parameter :: ice_density=910.    !< densite de la glace pour conversion en masse

contains

subroutine init_bilan_flux
 ! 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.

  dtoutflux=1. ! initMIP sorties flux tous les 1 an
  dtsum=0.
  bm_flux(:,:)=0.
  bmeltgr_flux(:,:)=0.
  bmeltfl_flux(:,:)=0.
  dhdt_flux(:,:)=0.
  calving_flux(:,:)=0.
  grline_flux(:,:)=0.
  calv_abl_bord_flux(:,:)=0.
  vieuxH_flux(:,:) = H(:,:)
  
end subroutine init_bilan_flux

	


subroutine bilan_flux_output

dtsum = dtsum + dt


! points englacés
where (ice(2:nx-1,2:ny-1).eq.1)
   bm_flux(2:nx-1,2:ny-1) = bm_flux(2:nx-1,2:ny-1) + bm(2:nx-1,2:ny-1)*dt          ! somme Bm sur dt
   where (flot(2:nx-1,2:ny-1))
     bmeltfl_flux(2:nx-1,2:ny-1) = bmeltfl_flux(2:nx-1,2:ny-1) + bmelt(2:nx-1,2:ny-1)*dt ! + ablbord(2:nx-1,2:ny-1)  ! somme bmelt flot sur dt
   elsewhere
     bmeltgr_flux(2:nx-1,2:ny-1) = bmeltgr_flux(2:nx-1,2:ny-1) + bmelt(2:nx-1,2:ny-1)*dt ! + ablbord(2:nx-1,2:ny-1)  ! somme bmelt grounded sur dt
   endwhere
   
   dhdt_flux(2:nx-1,2:ny-1) = dhdt_flux(2:nx-1,2:ny-1) + (H(2:nx-1,2:ny-1)-vieuxH_flux(2:nx-1,2:ny-1))       ! somme dhdt  sur dt    
endwhere

! points non englacés et flottant
where (ice(2:nx-1,2:ny-1).eq.0.and.flot(2:nx-1,2:ny-1))
     calv_abl_bord_flux(2:nx-1,2:ny-1) = calv_abl_bord_flux(2:nx-1,2:ny-1) - ablbord(2:nx-1,2:ny-1)
endwhere

calving_flux(2:nx-1,2:ny-1) = calving_flux(2:nx-1,2:ny-1) + calv(2:nx-1,2:ny-1)

where (gr_line(:,:).eq.1)
!~   grline_flux(:,:)= (((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_flux(:,:)
   grline_flux(:,:)= - 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 /(dx*dy) + grline_flux(:,:)
endwhere


! calcul de la moyenne des flux sur dtoutflux
!if ((mod(abs(time),dtoutflux).lt.dtmin).and.(abs(dtsum-dtoutflux).lt.dtmin)) then
if (mod(abs(time),dtoutflux).lt.dtmin) then
  bm_flux(:,:)=(bm_flux(:,:)/dtsum)*ice_density/secyear
  bmeltgr_flux(:,:)=(bmeltgr_flux(:,:)/dtsum)*ice_density/secyear
  bmeltfl_flux(:,:)=(bmeltfl_flux(:,:)/dtsum)*ice_density/secyear
  dhdt_flux(:,:)=(dhdt_flux(:,:)/dtsum)/secyear
  calving_flux(:,:)=(calving_flux(:,:)/dtsum)*ice_density/secyear
  grline_flux(:,:)= (grline_flux(:,:)/dtsum)*ice_density/secyear 
  calv_abl_bord_flux(:,:)=(calv_abl_bord_flux(:,:)/dtsum)*ice_density/secyear

! sorties dans debug_3d:
  debug_3d(:,:,106) = bm_flux(:,:)         ! acabf
  debug_3d(:,:,107) = -bmeltgr_flux(:,:)   ! libmassbfgr
  debug_3d(:,:,108) = dhdt_flux(:,:)       ! dlithkdt
  debug_3d(:,:,109) = calving_flux(:,:)    ! licalvf
  debug_3d(:,:,110) = grline_flux(:,:)     ! ligroundf
  debug_3d(:,:,121) = -bmeltfl_flux(:,:)   ! libmassbffl
  debug_3d(:,:,123) = calv_abl_bord_flux(:,:) + calving_flux(:,:) ! lifmassbf
! remise a 0
  dtsum = 0.
  bm_flux(:,:)      = 0.
  bmeltgr_flux(:,:) = 0.
  dhdt_flux(:,:)    = 0.
  calving_flux(:,:) = 0.
  grline_flux(:,:)  = 0.
  bmeltfl_flux(:,:) = 0.
  calv_abl_bord_flux(:,:) = 0.
  
!~ else if (abs(dtsum-dtoutflux).lt.dtmin) then
!~   print*,'pb sorties flux dtoutflux : time=', time
!~   print*,'dtsum=',dtsum,' dtoutflux', dtoutflux
!~ else if (mod(abs(time),dtoutflux).lt.dtmin) then
!~   print*,'pb sorties flux dtsum : time=', time
!~   print*,'dtsum=',dtsum,' dtoutflux', dtoutflux
endif

! reinitialisation des variables
vieuxH_flux(:,:)=H(:,:)
ablbord(:,:)=0.

end subroutine bilan_flux_output

end module bilan_flux_mod