!> \file ablation_month.f90
!! Pour le calcul de l'ablation mensuelle
!<

!> \namespace ablation_month
!! Calcule l ablation mensuelle
!! \author ...
!! \date ...
!! @note  Used modules:
!! @note  - declare_month
!! @note  - printtable
!<


module ablation_month

use declare_month
use printtable

implicit none

! parametres pour le bilan de masse

real,dimension(nx,ny) :: snow, melt
real                  :: frac, frac_min, frac_tot
integer               :: meth_abl_vi
! dimensionnements locaux
real ::  dtp                   !< integrating step for positive degree days (degrees)
real ::  pddct                 !< ct for PDD calculation
real ::  py                    !< ct for PDD calculation
real ::  s22                   !< ct for PDD calculation
real ::  nyear                 !< number of months in 1 year, st. dev. for temp *)

! parametres qui seront lus dans le fichier param

real ::  sigma                 !< variabilite Tday
real ::  CSI                   !< proportion of melted water that can refreeze *)
real ::  Cice                  !< melting factors for ice
real ::  Csnow                 !< melting factors for snow

!!$! aurel (greeneem) : pour comparer les deux approches de calcul de pdd ->
!!$real,dimension(nx,ny) :: pddann
!!$real,dimension(nx,ny) :: bmann


contains
!--------------------------------------------------------------------------------
!> SUBROUTINE: init_ablation
!! @note Attribution et initialisation des parametres lies a la methode de calcul
!! Used modules:
!! - module3D_phy
!>
subroutine init_ablation

namelist/ablation_month/Cice,Csnow,csi,sigma

! lecture des parametres

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

write(num_rep_42,*)'!___________________________________________________________'
write(num_rep_42,*)'!              ablation_month     PDD base Tann et Tjuly'
write(num_rep_42,ablation_month)                   ! pour ecrire les parametres lus
write(num_rep_42,*)
write(num_rep_42,*)
write(num_rep_42,*)'! Cice and Csnow, melting factors for ice and snow '
write(num_rep_42,*)'! sigma variabilite Tday'
write(num_rep_42,*)'! csi proportion of melted water that can refreeze '

! attribution des parametres lies a la methode de calcul
dtp=2.0
nyear=12
s22=0.5/sigma/sigma
py=2*pi/nyear
pddct=dtp/sigma/sqrt(2.*pi)/nyear*365.

end subroutine init_ablation

!-------------------------------------------------------
!> SUBROUTINE: Ablation 
!! @note Calcul l'ablation de chaque mois
!>
subroutine ablation

if (itracebug.eq.1)  call tracebug(' Passage dans ablation_month')

!------ Mass balance options

        meth_abl_vi=1    ! pdd Reeh
        !meth_abl_vi=2    ! Parametrization Thompson and Pollard 1997



! Calcul du PDD
!----------------------------
do i=1,nx
   do j=1,ny
        summ=0.0

! On a deja calcule Tm_surf(i,j,m) : temperature moyenne de chaque mois

        month: do mo=1,mois                    ! boucle sur les mois

            temp=0.0                             ! variable d'integration
            average_day: do k=1,50               ! pdd d'un jour par mois
               if (temp.gt.tm_surf(i,j,mo)+2.5*sigma) goto 209

               summ=summ+temp* &
                    exp(-((temp-tm_surf(i,j,mo))*(temp-tm_surf(i,j,mo)))*s22)
               
               temp=temp+dtp                   ! pdtp pas d'integration

            end do average_day

209         continue

        end do month                             !  boucle sur le mois

        pdd(i,j)=summ*pddct                      ! pdd pour toute l'annee

  end do
end do
debug_3d(:,:,58)=pdd(:,:)    

!!$
!!$!_____________________________________________________________
!!$!
!!$! aurel : je rajoute un calcul de PDD avec Tann et Tjuly pour comparer les deux approches
!!$! positive degree day
!!$do j=1,ny
!!$   do i=1,nx
!!$        summ=0.0
!!$        do nday=1,nyear   
!!$           tt(nday)=tann_par(i,j)+(tjuly_par(i,j)-tann_par(i,j))*cos(py*nday) 
!!$           temp=0.0 
!!$           do k=1,50
!!$              if (temp.gt.tt(nday)+2.5*sigma) goto 205
!!$              summ=summ+temp*exp(-(temp-tt(nday))*(temp-tt(nday))*s22) 
!!$              temp=temp+dtp
!!$           end do
!!$205        continue
!!$        end do
!!$        pddann(i,j)=summ*pddct 
!!$     end do
!!$  end do
!!$
!!$!sorties debug3d
!!$debug_3d(:,:,43)=pdd(:,:)      
!!$debug_3d(:,:,44)=pdd(:,:)-pddann(:,:)


!-----------------------

! calcul du bilan de masse :
!

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

!        Positive degrees required to melt the snow layer
!
     snow(i,j)=precip(i,j)           ! voir routine accum_month

     pds=snow(i,j)/Csnow             ! nombre de pdd pour fondre cette neige

     simax=snow(i,j)*csi             ! csi defini dans initial-phy-2.f90 (0.6)
                                      ! simax : Maximum amount of super. ice that can be formed

     pdsi=simax/cice                 ! Pos. degrees required to melt the superimposed ice



!-------- Ablation et Bilan de masse Reeh (1991)
!-----------------------------------------------------------------

     if (meth_abl_vi==1) then


        ! aurel 11/09 on rajoute l'impact du regel sur la temperature de surface.
        if (pdd(i,j).le.csi*pds) then
           sif=pdd(i,j)*csnow
        else 
           sif=simax 
        endif
        !     impact of refreesing on surface temperature 
        ts(i,j)=(tann(i,j)+26.6*sif)
        ts(i,j)=min(0.0,ts(i,j))
        tshelf(i,j)=ts(i,j)


           ! if pdd*csnow < thickness of melted snow that refreze: no ablation
           if (pdd(i,j).le.csi*pds) &
              bm(i,j)=acc(i,j)                                  ! toute la neige fondue regele

           ! if melted snow that refreze < pdd*csnow < snow
           if ((csi*pds.lt.pdd(i,j)).and.(pdd(i,j).le.pds))  &
              bm(i,j)=acc(i,j)+simax-pdd(i,j)*csnow             ! simax regel, le reste s'en va

           ! if snow < pdd*csnow <snow (pds) + surimposed ice (pdsi)
           if ((pds.lt.pdd(i,j)).and.(pdd(i,j).le.pds+pdsi)) &
              bm(i,j)=simax-(pdd(i,j)-pds)*cice                 ! attaque la glace de regele

           ! melting of old ice bm<0 with the remaining pdd
           if (pds+pdsi.le.pdd(i,j)) &
              bm(i,j)=(pds+pdsi-pdd(i,j))*cice                  ! attaque la vieille glace


           abl(i,j) = acc(i,j) - bm(i,j)                        ! Ablation
           frac=0.6


!-------- Energy surface mass balance (Thompson and Pollard 1997)
!----------------------------------------------------------------

    else if (meth_abl_vi==2) then               ! abla vince 2006

! aurel 11/09 on rajoute l'impact du regel sur la temperature de surface.
        if (pdd(i,j).le.csi*pds) then
           sif=pdd(i,j)*csnow
        else 
           sif=simax 
        endif
!     impact of refreesing on surface temperature 
        ts(i,j)=(tann(i,j)+26.6*sif)
        ts(i,j)=min(0.0,ts(i,j))
        tshelf(i,j)=ts(i,j)


        melt(i,j)=pdd(i,j)*Csnow
        melt(i,j)=min(snow(i,j),melt(i,j))

!----- refreezing fraction of precipiations (liquid and solid)

        if (snow(i,j)==0.) then
           frac=1
        else
           frac=1-((melt(i,j)/snow(i,j))-.7)/.3
           frac_min=min(1.,frac)
           frac_tot=max(0.,frac_min)

        endif

!------ Ablation

        if (frac_tot.lt.1.) then
           abl(i,j)= frac_tot*melt(i,j)

           if (pdd(i,j).gt.pds ) then   ! enough pdd to totaly melt the snow

              melt(i,j) = melt(i,j)+(pdd(i,j)-pds)*cice
              abl(i,j)  = abl(i,j)+(pdd(i,j)-pds)*cice

           endif

        else
           abl(i,j)=melt(i,j)
        endif

!------ Accumulation and mass balance

    acc(i,j)=snow(i,j)

       if (frac_tot.lt.1. .and. frac_tot.gt.0.) then
          acc(i,j)= acc(i,j) + (precip(i,j)-snow(i,j))
          !acc(i,j)= acc(i,j) + (precip(i,j)*48-precip(i,j))*(1-frac_tot)  ! voir la routine accum_month
       endif

    bm(i,j) = acc(i,j) - abl(i,j)

    endif

  end do
end do

!!$!_____________________________________________________________
!!$!
!!$! aurel : je rajoute un calcul de bm avec Tann et Tjuly pour comparer les deux approches
!!$do j=1,ny
!!$   do i=1,nx
!!$      pds=acc(i,j)/csnow        !  positive degrees required to melt the snow laye
!!$      simax=acc(i,j)*csi        !  maximum amount of super. ice that can be formed
!!$      pdsi=simax/cice           !  pos. degrees required to melt the superimposed ice
!!$      if (pddann(i,j).le.csi*pds) then
!!$         sif=pddann(i,j)*csnow
!!$      else 
!!$         sif=simax 
!!$      endif
!!$!      mass balance
!!$        if (pddann(i,j).le.csi*pds)  bmann(i,j)=acc(i,j) 
!!$        if ((csi*pds.lt.pddann(i,j)).and.(pddann(i,j).le.pds))       &
!!$             bmann(i,j)=acc(i,j)+simax-pddann(i,j)*csnow 
!!$        if ((pds.lt.pddann(i,j)).and.(pddann(i,j).le.pds+pdsi))      &
!!$             bmann(i,j)=simax-(pddann(i,j)-pds)*cice
!!$        if (pds+pdsi.le.pddann(i,j))   bmann(i,j)=(pds+pdsi-pddann(i,j))*cice
!!$     end do
!!$  end do
!!$! sortie debug
!!$debug_3d(:,:,45)=bm(:,:)-bmann(:,:)


!!$!aurel : pour experience de calibration : je force avec les donnees de bilan de masse MAR
!!$do j=1,ny
!!$   do i=1,nx
!!$      bm(i,j)=sum(bm_mar(i,j,:))/12
!!$   end do
!!$end do


end subroutine ablation

end module ablation_month