module  climat_perturb_mois_lapsecouche_mod


use module3d_phy
implicit none

integer nft             ! nombre de lignes a lire dans le fichier forcage climatique
real,dimension(:),allocatable :: tdate          ! time for climate forcing
real,dimension(:),allocatable :: tpert          ! temperature for climate forcing
real,dimension(:),allocatable :: spert          ! sea surface perturbation

real :: coefT                   ! pour modifier l'amplitude de la perturb. T
real :: rappact                 ! pour le calcul du rapport 'accumulation
integer :: retroac              ! 1-> full retroactions accum
real :: rapbmshelf              ! pour calcul coefbmshelf
real :: mincoefbmelt            ! butoirs mini 
real :: maxcoefbmelt            ! butoirs maxi de coefbmelt
character(len=80) :: filforc    ! nom du fichier forcage

! Pour l'instant tafor est global meme si inutile si on utilise
! un forcage  (variation spatiales)

contains

!--------------------------------------------------------------------------------
subroutine input_clim !routine qui permet d'initialiser les variables climatiques

implicit none
character(len=8) :: control      !label to check clim. forc. file (filin) is usable
character(len=80):: filin


! Lecture du forcage
!-----------------------
! Le fichier de forcage est lu dans le fichier entree parametres

filin=trim(dirforcage)//trim(filforc)

open(num_forc,file=filin,status='old')

read(num_forc,*) control,nft

! Determination of file size (line nb), allocation of perturbation array

if (control.ne.'nb_lines') then
   write(6,*) filin,'indiquer le nb de ligne en debut de fichier:'
   write(6,*) 'le nb de lignes et le label de control nb_lines'
   stop
endif

! Dimensionnement des tableaux tdate, ....
#if defined key_agrif
   if (.not.associated(tdate)) then
#else
   if (.not.allocated(tdate)) then
#endif

      allocate(tdate(nft),stat=err)
      if (err/=0) then
         print *,"erreur a l'allocation du tableau Tdate",err
         stop 4
      end if
   end if

#if defined key_agrif
  if (.not.associated(spert)) then
#else
  if (.not.allocated(spert)) then
#endif
     allocate(spert(nft),stat=err)
     if (err/=0) then
        print *,"erreur a l'allocation du tableau Spert",err
        stop 4
     end if
  end if

#if defined key_agrif
  if (.not.associated(tpert)) then
#else
  if (.not.allocated(tpert)) then
#endif
     allocate(tpert(nft),stat=err)
     if (err/=0) then
        print *,"erreur a l'allocation du tableau Tpert",err
        stop 4
     end if
  end if

do i=1,nft
   read(num_forc,*) tdate(i),spert(i),tpert(i)
end do
close(num_forc)

tpert(:)=tpert(:)*coefT

end subroutine input_clim

!--------------------------------------------------------------------------------
subroutine init_forclim


namelist/clim_pert/coefT,rappact,retroac,rapbmshelf,mincoefbmelt,maxcoefbmelt,filforc

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

! formats pour les ecritures dans 42
428 format(A)

write(num_rep_42,428)'!___________________________________________________________' 
write(num_rep_42,428) '&clim_pert                      ! module climat_perturb_mod '
write(num_rep_42,*)
write(num_rep_42,*) 'coefT        = ', coefT 
write(num_rep_42,*) 'rappact      = ', rappact 
write(num_rep_42,*) 'rretroac     = ', retroac 
write(num_rep_42,*) 'rapbmshelf   = ', rapbmshelf  
write(num_rep_42,*) 'mincoefbmelt = ', mincoefbmelt
write(num_rep_42,*) 'maxcoefbmelt = ', maxcoefbmelt
write(num_rep_42,'(A,A)') ' filforc      = ', filforc
write(num_rep_42,*)'/'                            
write(num_rep_42,*)

return
end subroutine init_forclim

!--------------------------------------------------------------------------------
subroutine forclim               !  au temps considere (time) attribue les scalaires
                                 !  tafor : forcage en temperature
                                 !  sealevel : forcage niveau des mers
                                 !  coefbmelt : forcage fusion basale ice shelves 

use module3d_phy
implicit none

!       time en dehors des limites du fichier forcage
if(time.lt.tdate(1)) then
   tafor=tpert(1)
   sealevel=spert(1)
   ift=1

else if (time.ge.tdate(nft)) then
   tafor=tpert(nft)
   sealevel=spert(nft)
   ift=nft

else
   do i=1,nft-1 
      if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then ! entre i et i+1 : cas general
         tafor=tpert(i)+(tpert(i+1)-tpert(i))*       &
              (time-tdate(i))/(tdate(i+1)-tdate(i))
         sealevel=spert(i)+(spert(i+1)-spert(i))*    &
              (time-tdate(i))/(tdate(i+1)-tdate(i))
         ift=i
         goto 100
      endif
   end do
endif
100   continue

        
!  coefmshelf est un coefficient qui fait varier bmgrz et bmshelf
!  en fonction de tafor 

!   coefbmshelf=(1.+tafor/10.)           ! coefbmshelf=0 pour tafor=-10deg

if ((tafor.le.0).and.(tafor.gt.-5.)) then
   coefbmshelf=(1.+tafor/rapbmshelf)     ! coefbmshelf=0 pour tafor=-7 standard precedent

else if (tafor.le.-5) then               ! lineaire en 0 a -10 et la valeur a -5
   coefbmshelf=(1.-5./rapbmshelf)/5.*(tafor+10)

else
   coefbmshelf=(1.+5.*tafor/rapbmshelf)  ! 5 fois plus efficace vers le chaud
endif

coefbmshelf=max(coefbmshelf,mincoefbmelt)  
coefbmshelf=min(coefbmshelf,maxcoefbmelt)

#if defined key_agrif
   if ( .not. agrif_root() ) then
      geoplace = agrif_parent(geoplace)
   endif
#endif



   call massb_perturb_mois_lapsecouche

end subroutine forclim

end module  climat_perturb_mois_lapsecouche_mod