subroutine massb_perturb_mois_lapsecouche              ! calcule le mass balance en mode perturbation
                                         
  use climat_perturb_mois_lapsecouche_mod
  use declare_month_lapsecouche

  implicit none


  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 utilise un lapserate par tranche d'altitude
            if(ZS(i,j).gt.3000) then
               T_surf_ref=Tsol(i,j,mo)+(lapserate(mo,4)*(ZS(i,j)-3000)+lapserate(mo,3)*1000+lapserate(mo,2)*1000+lapserate(mo,1)*1000)
            end if
            if((ZS(i,j).le.3000).and.(zel.gt.2000)) then
               T_surf_ref=Tsol(i,j,mo)+(lapserate(mo,3)*(ZS(i,j)-2000)+lapserate(mo,2)*1000+lapserate(mo,1)*1000)
            end if
            if ((ZS(i,j).le.2000).and.(zel.gt.1000)) then
               T_surf_ref=Tsol(i,j,mo)+(lapserate(mo,2)*(ZS(i,j)-1000)+lapserate(mo,1)*1000)
            end if
            if (ZS(i,j).le.1000) then
               T_surf_ref=Tsol(i,j,mo)+lapserate(mo,1)*ZS(i,j)
            end if
            
            T_surf_ref=T_surf_ref+Tafor  ! on applique la perturbation homogene

            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)))       
            else if (retroac.eq.0) then
               Pm_surf(i,j,mo)=Pm_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
         Tjuly(i,j)= Tm_surf(i,j,7)
      end do
   end do
   
   
   debug_3D(:,:,59)=Tjuly(:,:)



call accum_month_lapsecouche()    ! 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_lapsecouche