!> \file massb_perturb_mois.f90
!!Module pour le calcul du bilan de masse en mode perturbation
!<

!> SUBROUTINE: massb_perturb_mois()
!! \author ... 
!! \date ...
!! @note Cette routine permet de calculer le mass balance en mode perturbation
!! @note Used modules:
!! @note    - use climat_perturb_mois_mod
!! @note    - use declare_month
!<
subroutine massb_perturb_mois              ! calcule le mass balance en mode perturbation
                                         
  use climat_perturb_mois_mod
  use declare_month

  implicit none


! variable de travail local
!  real,dimension(nx,ny,mois+1) :: Tm_time
!  real,dimension(nx,ny,mois+1) :: Pm_time ! Pm_time n'est pas utilise


  real zel     ! variable de travail
  real T_surf_ref



   do j=1,ny
      do i=1,nx
! Zs est l'altitude de la surface qu'elle soit mer, glace ou lac
         ZS(i,j)=max(slv(i,j),S(I,J))

         do mo=1,mois

! on corrige la temperature par rapport a l'alti de ref
            T_surf_ref = Tm_0(i,j,mo) + lapserate(i,j,mo) * (ZS(i,j)-S0CLIM(i,j))

! on ajoute le terme de climat-perturb
            T_surf_ref = T_surf_ref+Tafor

! on corrige egalement les precipes si on considere les retroactions
            if(retroac.eq.1) then      ! full retroaction acc.
               Pm_surf(i,j,mo)= Pm_surf(i,j,mo)*exp(rappact*(T_surf_ref-Tm_surf(i,j,mo)))
            end if
            Tm_surf(i,j,mo)=T_surf_ref

         end do
      end do
   end do

! Pour les modulos
   Tm_surf(:,:,mois+1)=Tm_surf(:,:,1)
   Pm_surf(:,:,mois+1)=Pm_surf(:,:,1)


   do j=1,ny
      do i=1,nx
         Tann(i,j)=(sum(Tm_surf(i,j,:))-Tm_surf(i,j,1))/12   ! janvier compte en double
         !Tann(i,j)=sum(Tm_surf(i,j,:))/12
         Tjuly(i,j)= Tm_surf(i,j,7)
      end do
   end do

   debug_3D(:,:,59)=Tjuly(:,:)
! aurel : certainement pas a mettre la mais pour le moment...
   do j=1,ny
      do i=1,nx
         debug_3D(i,j,60)=sqrt(ux(i,j,1)**2+uy(i,j,1)**2)
      end do
   end do


call accum_month()    ! on calcule l'accumulation deduite de ce champ (temp, precipes)



!     ablation et bilan de masse sont maintenant appeles par forclim



!--------------------------------------------------------------------------------------
! etape 3 : calcul tann et tjuly // sortie de debug3D

!!$do j=1,ny
!!$   do i=1,nx
!!$      Tann(i,j)=sum(Tm_surf(i,j,:))/12
!!$      Tjuly(i,j)=(Tm_surf(i,j,6)+Tm_surf(i,j,7)+Tm_surf(i,j,8))/3
!!$      zel=max(S(i,j),20.*(ylat(i,j)-65.))
!!$      Tann_par(i,j)=49.13-0.007992*Zel-0.7576*ylat(i,j)
!!$      Tjuly_par(i,j)= 30.78-0.006277*S(i,j)-0.3262*ylat(i,j)
!!$   end do
!!$end do

!!$debug_3D(:,:,29)=Tann(:,:)-Ta0(:,:)
!!$debug_3D(:,:,30)=Acc(:,:)-precip(:,:)
!!$debug_3D(:,:,31)=Acc(:,:)-bmelt(:,:)
!!$debug_3D(:,:,46)=Tann(:,:)-Tann_par(:,:)
!!$debug_3D(:,:,47)=Tjuly(:,:)-Tjuly_par(:,:)
!!$debug_3D(:,:,48)=Tm_surf(:,:,7)-Tjuly_par(:,:)

return
end subroutine massb_perturb_mois