!> \file ablation_ann_july_mod.f90
!! Module pour l'initialisation et le calcul de l'ablation annuelle 
!<

!> \namespace ablation_ann
!! Initialisation et calcul de l'ablation annuelle
!! \author ...
!! \date ...
!! @note Used module
!! @note use module3D_phy
!! @todo essayer de  mettre certaines autres variables en local.summ, tt,pdd,pds,simax,pdsi, .....
!<
     
module ablation_ann             ! positive degree day  base annuelle

                                ! remarque : on pourrait mettre certaines autres variables en local.
                                ! summ, tt,pdd,pds,simax,pdsi, .....

use module3d_phy
implicit none

! 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 

contains
!-------------------------------------------------------------------------------

!> SUBROUTINE: init_ablation
!! @note Attribution des parametres lies a la methode de calcul 
!>
subroutine init_ablation

namelist/ablation_ann/Cice,Csnow,csi,sigma

! lecture des parametres

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

write(num_rep_42,*)'!___________________________________________________________' 
write(num_rep_42,*)'!             module ablation_ann     PDD base Tann et Tjuly'
write(num_rep_42,ablation_ann)                   ! 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 de l'ablation annuelle july 
!>
subroutine ablation

! positive degree day
do j=1,ny
   do i=1,nx
        summ=0.0
        do nday=1,nyear   
           tt(nday)=tann(i,j)+(tjuly(i,j)-tann(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
        pdd(i,j)=summ*pddct 
     end do
  end do

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 (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)

!      mass balance
        if (pdd(i,j).le.csi*pds)  bm(i,j)=acc(i,j) 

        if ((csi*pds.lt.pdd(i,j)).and.(pdd(i,j).le.pds))       &
             bm(i,j)=acc(i,j)+simax-pdd(i,j)*csnow 

        if ((pds.lt.pdd(i,j)).and.(pdd(i,j).le.pds+pdsi))      &
             bm(i,j)=simax-(pdd(i,j)-pds)*cice

        if (pds+pdsi.le.pdd(i,j))   bm(i,j)=(pds+pdsi-pdd(i,j))*cice

     end do
  end do

debug_3D(:,:,31)=bm(:,:)-bmelt(:,:)

end subroutine ablation

end module ablation_ann