Index: /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90
===================================================================
--- /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90	(revision 42)
+++ /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90	(revision 42)
@@ -0,0 +1,622 @@
+!> \file climat_GrIce2sea_mod.f90
+! forcage avec BM
+!! Module ou les variations temporelles des variables climatiques
+!! sont directement imposees
+!<
+
+!> \namespace  climat_Grice2sea_mod
+!! Module ou les variations temporelles des variables climatiques
+!! sont directement imposees
+!! \author Cat
+!! \date 31 oct
+!! @note Used modules:
+!! @note   - use module3D_phy
+!<
+module climat_Grice2sea_years_perturb_mod
+
+
+use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,num_forc,dirforcage,dirnameinp,tafor,time,sealevel
+!use lect_climref_Ice2sea
+use netcdf
+use io_netcdf_grisli
+implicit none
+
+
+real :: T_lapse_rate              !< pour la temperature
+
+
+! declarations specifiques year by year
+
+real,dimension(:),allocatable        :: time_snap         !> date des snapshots  indice : nb de nb_snap)
+real,dimension(:,:,:),allocatable    :: smb_snap          !> bilan de masse des snapshots indices : nx,ny,nb_snap 
+real                                 :: ecart_snap        !> ecart temporel entre les snapshots
+real                                 :: time_depart_snaps !> temps du debut premier snapshot
+integer                              :: nb_snap           !> nombre de snapshots
+
+! declaration pour le bilan de masse
+real,dimension(nx,ny)                :: bm_0              !> bilan de masse initial equ. S_0ref de Tamsin
+real,dimension(nx,ny)                :: bm_time           !> bilan de masse interpole entre 2 snapshots (~ S_t^RCM)
+
+! calcul du gradient
+real,dimension(nx,ny)                :: bm_ref_t          !> bilan de masse de reference au temps t
+real,dimension(nx,ny,10)             :: bm_ref_10         !> tableau des bm_ref pour moyenne 10 ans
+real                                 :: time_prec         !> temps du precedent snapshot
+integer                              :: icum              !> pour le test stockage dans bm_ref_10 
+integer                              :: i_moy             !> nombre de snapshots stockes
+
+real,dimension(nx,ny)                :: grad_bm           !> gradient du bilan de masse 
+real,dimension(nx,ny)                :: S_ref             !> surface de reference
+
+real                                 :: grad_N_smb_neg    !> SMB < 0 North
+real                                 :: grad_N_smb_pos    !> SMB >= 0 North
+real                                 :: grad_S_smb_neg    !> SMB < 0 South
+real                                 :: grad_S_smb_pos    !> SMB >= 0 South
+real                                 :: coef_smb_unit     ! pour corriger l'unite
+
+real,dimension(nx,ny)                :: TA0               !> Temp annuelle sur S0
+
+
+! aurel, pour climat type perturb:
+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
+character(len=120) :: filforc    !< nom du fichier forcage
+
+
+
+! pour les lectures ncdf
+real*8, dimension(:,:,:),  pointer   :: tab3D             !< tableau 3d real pointer
+real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
+Character(len=10), dimension(3)      :: dimtestname       !> pour la sortie test netcdf
+integer                              :: ncidloc           !> pour la sortie test netcdf
+integer                              :: status            !> pour la sortie test netcdf
+
+integer                              :: massb_time        !< pour selectionner le type de calcul de smb
+! On a deux routines : celle avec un seul fichier (donnees observees) : massb_Ice2sea_fixe
+! Celle avec le bilan de masse sur plusieurs snapshots (annuels par ex.) entre lesquels on interpole.
+
+contains
+
+!--------------------------------------------------------------------------------
+!> SUBROUTINE: input_clim  
+!! Routine qui permet d'initialiser les variations temporelles des variables climatiques
+!>
+!--------------------------------------------------------------------------------
+
+subroutine input_clim  
+
+  character(len=100) :: smb_file           ! fichier smb
+  character(len=100) :: temp_annual_file   ! temperature annuelles
+  character(len=100) :: file_smb_snap      !> nom du fichier dans lequel sont les snapshots
+
+  integer            :: err                ! recuperation d'erreur
+  integer            :: i,j,k
+
+  !aurel for Tafor:
+  character(len=8) :: control      !label to check clim. forc. file (filin) is usable
+  character(len=80):: filin
+
+  namelist/clim_smb_T_gen/smb_file,coef_smb_unit,temp_annual_file
+
+
+428 format(A)
+  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
+  read(num_param,clim_smb_T_gen)
+
+  write(num_rep_42,428)'!________________________________________________________________' 
+  write(num_rep_42,428)'!  module climat_Grice2sea_years_mod lecture climat ref          '
+  write(num_rep_42,clim_smb_T_gen)
+  write(num_rep_42,428)'! smb_file          = fichier SMB (kg/m2/an)                     '
+  write(num_rep_42,428)'! coef_smb_unit     = coef passage m glace/an  (1/910 ou 1/918)  '
+  write(num_rep_42,428)'! temp_annual_file  = Temp moy annuelle  (°C)                    '
+  write(num_rep_42,428)'!________________________________________________________________'
+
+  
+! smb : surface mass balance
+  smb_file  = trim(dirnameinp)//trim(smb_file)
+
+  call Read_Ncdf_var('z',smb_file,tab)
+
+!  call lect_input(3,'smb',1,bm,smb_file,trim(dirnameinp)//trim(runname)//'.nc')
+
+  bm(:,:)  = tab(:,:) * coef_smb_unit
+
+!  where ((H(:,:).lt.1.).and.(Bsoc(:,:).gt.0.))
+!     bm(:,:) = bm(:,:) - 5.                     ! pour faire un masque a l'exterieur du Groenland actuel
+!  end where
+
+!cdc test debug Hemin15 et Greeneem15
+!  where (bm(:,:).lt.-1000) bm(:,:)=0.
+
+  acc(:,:) = 0.
+  abl(:,:) = 0.
+
+  where (bm(:,:).gt.0.)
+     acc(:,:) = bm(:,:)   ! accumulation quand positif
+  elsewhere
+     abl(:,:) = - bm(:,:) ! ablation quand negatif
+  end where
+
+
+! surface temperature  Tann
+
+  temp_annual_file = trim(dirnameinp)//trim(temp_annual_file)
+
+  call Read_Ncdf_var('z',temp_annual_file,tab)
+  Tann(:,:)  = tab(:,:)
+!  call lect_input(3,'Tann',1,Tann,temp_annual_file,trim(dirnameinp)//trim(runname)//'.nc')
+
+!cdc test debug Hemin15
+!  where (Tann(:,:).lt.-100.) Tann(:,:)=10.
+
+
+  ta0(:,:)   = Tann(:,:)
+  Tjuly(:,:) = Tann(:,:)
+
+
+
+!lecture namelist 
+
+  namelist/clim_snap/nb_snap,time_depart_snaps,ecart_snap,file_smb_snap,massb_time
+
+  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+  read(num_param,clim_snap)
+
+! formats pour les ecritures dans 42
+  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
+  read(num_param,clim_snap)
+
+  write(num_rep_42,428)'!_______________________________________________________________________' 
+  write(num_rep_42,428)'!  module climat_Grice2sea_years_mod                                    '
+  write(num_rep_42,clim_snap)
+  write(num_rep_42,428)'! nb_snap           = nombre de snapshots                               '
+  write(num_rep_42,428)'! time_depart_snaps = debut du forçage                                  '
+  write(num_rep_42,428)'! ecart_snap        = ecart entre les snapshots                         '
+  write(num_rep_42,428)'! file_smb_snap     = fichier serie temp anomalie SMB de GCM            '
+  write(num_rep_42,428)'! massb_time        = 0:fixe, 1:interp snapshots, 2:snapsh+interp vert  '
+  write(num_rep_42,428)'!_______________________________________________________________________' 
+
+  if (massb_time == 1) then ! lecture des snapshots
+! allocation dynamique de time_snap
+     if (allocated(time_snap)) then 
+        deallocate(time_snap,stat=err)
+        if (err/=0) then
+           print *,"Erreur à la desallocation de time_snap",err
+           stop 
+        end if
+     end if
+
+     allocate(time_snap(nb_snap),stat=err)
+     if (err/=0) then
+        print *,"erreur a l'allocation du tableau time_snap ",err
+        print *,"nb_snap = ",nb_snap
+        stop 
+     end if
+
+! remplissage de time_snap
+!write(6,*) 'time_snap'
+     do i=1,nb_snap
+        time_snap(i) = time_depart_snaps + ecart_snap * (i-1)
+!   write(6,*) i,time_snap(i)
+     end do
+
+! allocation dynamique de smb_snap
+     if (allocated(smb_snap)) then 
+        deallocate(smb_snap,stat=err)
+        if (err/=0) then
+           print *,"Erreur à la desallocation de smb_snap",err
+           stop 
+        end if
+     end if
+
+     allocate(smb_snap(nx,ny,nb_snap),stat=err)
+     if (err/=0) then
+        print *,"erreur a l'allocation du tableau smb_snap ",err
+        print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap
+        stop 
+     end if
+
+! lecture de smb_snap
+     file_smb_snap  = trim(dirnameinp)//trim(file_smb_snap)
+     call Read_Ncdf_var('z',file_smb_snap,tab3D)
+     smb_snap (:,:,:) = Tab3D(:,:,:) * coef_smb_unit
+
+! ce sont des anomalies : ajoute les valeurs de reference
+     do k = 1,nb_snap
+        do j = 1,ny
+           do i = 1,nx
+              smb_snap (i,j,k) = smb_snap(i,j,k) + bm(i,j)
+           end do
+        end do
+     end do
+
+! copie la valeur de reference dans bm_0
+     bm_0(:,:) = bm(:,:)
+  endif
+
+  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 (.not.allocated(tdate)) then
+     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 (.not.allocated(spert)) then
+     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 (.not.allocated(tpert)) then
+     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
+  
+
+  
+! ecriture de verification
+!file_smb_snap = 'test_sortie_smb.nc'
+!file_smb_snap  = trim(dirnameout)//trim(file_smb_snap)
+
+!ncidloc = 501
+
+!call  lect_netcdf_type
+!write(6,*) 'ncdf_type'
+
+!if (ncdf_type.eq.32) then
+!   status  = nf90_create(TRIM(file_smb_snap),NF90_WRITE,ncidloc)    ! ouverture du fichier
+!else if (ncdf_type.eq.64) then
+!   status  = nf90_create(trim(file_smb_snap),and(nf90_write,nf90_64bit_offset),ncidloc) ! r2d2
+!end if
+
+!   status  = nf90_close(ncidloc)                                 ! fermeture
+
+! call write_ncdf_dim('x',trim(file_smb_snap),nx)               ! dimensions des tableaux
+! call write_ncdf_dim('y',trim(file_smb_snap),ny)
+! call write_ncdf_dim('time',trim(file_smb_snap),nb_snap)
+
+!Tab3D(:,:,:) = smb_snap (:,:,:)
+!dimtestname(1) = 'x'
+!dimtestname(2) = 'y'
+!dimtestname(3) = 'time'
+
+!call write_ncdf_var(trim('smb'),dimtestname,trim(file_smb_snap),tab3D,'double')
+   
+end subroutine input_clim
+
+!--------------------------------------------------------------------------------
+!> SUBROUTINE: init_forclim
+!! Routine qui permet d'initialiser les variables climatiques au cours du temps
+!>
+subroutine init_forclim
+
+  implicit none
+  namelist/lapse_rates/T_lapse_rate
+  namelist/clim_pert_massb/coefT,filforc
+  
+  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+  read(num_param,lapse_rates)
+
+! formats pour les ecritures dans 42
+428 format(A)
+
+  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
+  read(num_param,lapse_rates)
+
+  write(num_rep_42,428)'!________________________________________________________________' 
+  write(num_rep_42,428)'!  module climat_Grice2sea_years_mod                             '
+  write(num_rep_42,lapse_rates)
+  write(num_rep_42,428)'!T_lapse_rate       = lapse rate temp annuelle                   '
+  write(num_rep_42,428)'!________________________________________________________________' 
+
+  
+  rewind(num_param)  
+  read(num_param,clim_pert_massb)
+  
+  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,'(A,A)') ' filforc      = ', filforc
+  write(num_rep_42,*)'/'
+  write(num_rep_42,*)
+
+
+! appelle la routine de lecture des smb annuels
+  call input_clim
+  if (massb_time == 1) then ! lecture gradients smb
+     call init_grad_smb
+  endif
+
+  return
+end subroutine init_forclim
+
+!--------------------------------------------------------------------------------
+!> SUBROUTINE: forclim
+!!  
+!!  Routine qui permet le calcul climatique au cours du temps
+!!  @note Au temps considere (time) attribue les scalaires
+!!  @note   - tafor : forcage en temperature
+!!  @note   - sealevel : forcage niveau des mers
+!!  @note   - coefbmelt : forcage fusion basale ice shelves 
+!>
+subroutine forclim               !  au temps considere (time) 
+
+  implicit none
+
+  integer i,j,ift
+  
+  select case (massb_time)
+     case(0)
+        ! smb fixe et Tann avec lapse rate + Tafor
+        
+        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
+
+        
+        Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor
+        Ts(:,:)    = Tann(:,:)
+     case(1)
+        call massb_Ice2sea_RCM
+     case default
+        print *, 'Numero de massb invalide dans climat_Grice2sea_years_mod'
+        stop
+     end select
+end subroutine forclim
+
+
+subroutine massb_Ice2sea_RCM                ! calcule le mass balance 
+
+  implicit none
+  integer             :: k_snap               ! pour calculer les indices de temps
+  real                :: time_bis             ! pour repliquer les dernieres annees
+  integer             :: k
+
+! calcul smb a partir fichier snapshots massb_Ice2sea_RCM
+! Calcule le mass balance d'apres un fichier de snapshots
+! avec la temperature parametree
+!     surface temperature et accumulation
+  Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) 
+  Ts(:,:)    = Tann(:,:)
+! calcule bm_time par interpolation entre deux snapshots
+! avant prend la valeur de reference
+! apres prend la derniere valeur
+! en general les snapshots vont de 2000 a 2200
+  if(time.lt.time_snap(1)) then              ! time avant le forcage
+     bm_time(:,:) = bm_0(:,:)
+     k_snap       = 1
+     S_ref(:,:)   = S(:,:)                   ! du coup sera la surface de reference avant le forcage
+     icum         = 0
+     i_moy        = 0
+     bm_ref_t(:,:)=  bm_0(:,:)               ! bilan de masse de reference
+     time_prec    = time
+  else if (time.ge.time_snap(nb_snap)) then  ! time apres le forcage
+     k_snap       = nb_snap
+     bm_time(:,:) =  smb_snap (:,:,k_snap)
+     if (abs(time-time_prec-1.).lt.dt) then   ! 
+        time_prec = time_prec + 1
+        icum      = 1
+     else
+        icum      = 0
+     endif
+  else                                       ! cas general
+     do k = 1 , nb_snap-1
+        if((time.ge.time_snap(k)).and.(time.lt.time_snap(k+1))) then ! entre k et k+1 
+           bm_time(:,:) = smb_snap(:,:,k)+(smb_snap(:,:,k+1)-smb_snap(:,:,k)) *   &
+                (time-time_snap(k))/(time_snap(k+1)-time_snap(k))
+! exactement sur le snapshot et avec un ecart 1 an par rapport au precedent stockage
+!           write(6,*) 'time,tests',k,time,time-time_snap(k),time-time_prec-1.
+           if ((abs(time-time_snap(k)).le.dt).and.(abs(time-time_prec-1.).lt.dt)) then    
+              k_snap    = k
+              icum      = 1
+              time_prec = time_snap(k)     ! time_prec est le temps du precedent
+                                                 ! snapshot garde
+           else
+              icum    = 0
+           endif
+           exit
+        endif
+     end do
+  endif
+  call  grad_smb              !-----------------------------> A faire
+
+  if (massb_time == 1) then ! pas d'interpolation verticale
+     bm(:,:) = bm_time(:,:)
+  else if (massb_time == 2) then ! interpolation verticale
+     ! ajuste bm en fonction du temps et du gradient
+     bm(:,:) = bm_time(:,:) + grad_bm(:,:) *(S(:,:) - S_ref(:,:))
+     write(6,897) time, time_prec, icum, i_moy
+897  format('test temps smb   ',2(f0.3,1x),2(i0,1x))
+  endif
+
+! garde les 10 dernieres annees et calcule la moyenne
+  if (icum.eq.1) then                           ! stockage dans le tableau bm_ref_10
+     do k = 9,1,-1
+        bm_ref_10(:,:,k+1) = bm_ref_10(:,:,k)   ! on decale tous les elements
+     end do
+     bm_ref_10(:,:,k+1) = bm(:,:)               ! le plus recent est en position 1
+     i_moy              = i_moy +1              ! compte combien il y en a pour la moyenne
+     i_moy              = min(i_moy,10)
+     bm_ref_t(:,:)   = 0.
+     do k = 1,i_moy
+        bm_ref_t(:,:) = bm_ref_t(:,:) + bm_ref_10(:,:,k)
+     end do
+     bm_ref_t(:,:)    = bm_ref_t(:,:)/i_moy
+     write(6,898) time, time_prec, icum, i_moy
+898  format('cumul pour gradient  ',2(f0.3,1x),2(i0,1x))
+  end if
+end subroutine massb_Ice2sea_RCM
+ 
+!------------------------------------------------------------------------------
+!> initialise le calcul du gradient vertical de smb
+subroutine init_grad_smb 
+
+  use module3D_phy
+  implicit none
+
+  character(len=120)       ::  file_grad_smb ! nom du fichier gradients de Tamsin
+  character(len=40)        ::  fin_ligne     ! fin de la ligne
+  real                     ::  grad
+
+  namelist/grad_smb/file_grad_smb
+
+! Dans lequel sont :
+! grad_N_smb_pos,grad_N_smb_neg,grad_S_smb_pos,grad_S_smb_neg,lim_lat
+
+  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+  read(num_param,grad_smb)
+
+
+! formats pour les ecritures dans 42
+428 format(A)
+
+  rewind(num_param)                     ! pour revenir au debut du fichier param_list.dat
+  read(num_param,grad_smb)
+
+  write(num_rep_42,428)'!________________________________________________________________' 
+  write(num_rep_42,428)'!  gradient smb climat_Grice2sea_years_mod                       '
+  write(num_rep_42,grad_smb)
+  write(num_rep_42,428)'!grad_smb       = fichier gradient SMB                           '
+  write(num_rep_42,428)'!________________________________________________________________' 
+
+
+  file_grad_smb=trim(dirnameinp)//'SMB-H-Feedback/gradients_18_07_2012/'//trim(file_grad_smb)
+  open(622,file=file_grad_smb)
+  do i=1,4
+     read(622,'(f9.4,A)') grad,fin_ligne
+     write(6,*) grad,fin_ligne
+
+     if (index(fin_ligne, "North").ne.0) then       ! North est dans la ligne
+        if (index(fin_ligne, "<").ne.0)  then       ! smb negatif
+           grad_N_smb_neg = grad
+        else if (index(fin_ligne, ">=").ne.0) then  ! smb positif
+           grad_N_smb_pos = grad
+        else
+           write(6,*) 'probleme lecture North fichier smb', file_grad_smb
+           STOP
+        end if
+
+     else if (index(fin_ligne, "South").ne.0) then   !South est dans la ligne
+        if (index(fin_ligne, "<").ne.0)  then        ! smb negatif
+           grad_S_smb_neg = grad
+        else if (index(fin_ligne, ">=").ne.0) then   ! smb positif
+           grad_S_smb_pos = grad
+        else
+           write(6,*) 'probleme lecture South fichier smb  ', file_grad_smb
+           STOP
+        end if
+
+     else 
+        write(6,*) 'probleme lecture ni North ni South fichier smb  ', file_grad_smb
+        write(6,*) 'fin_ligne',fin_ligne,'   index North', index(fin_ligne, "North"),'   index South', index(fin_ligne, "South")
+        STOP
+     end if
+  end do
+
+  write(6,*) 'coefficients lus '
+  write(6,*) 'grad_N_smb_pos', grad_N_smb_pos
+  write(6,*) 'grad_N_smb_neg', grad_N_smb_neg
+  write(6,*) 'grad_S_smb_pos', grad_S_smb_pos
+  write(6,*) 'grad_S_smb_neg', grad_S_smb_neg
+
+  grad_N_smb_pos = coef_smb_unit * grad_N_smb_pos
+  grad_N_smb_neg = coef_smb_unit * grad_N_smb_neg
+  grad_S_smb_pos = coef_smb_unit * grad_S_smb_pos
+  grad_S_smb_neg = coef_smb_unit * grad_S_smb_neg
+
+  write(6,*) 'coefficients *coef_smb_unit '
+  write(6,*) 'grad_N_smb_pos', grad_N_smb_pos
+  write(6,*) 'grad_N_smb_neg', grad_N_smb_neg
+  write(6,*) 'grad_S_smb_pos', grad_S_smb_pos
+  write(6,*) 'grad_S_smb_neg', grad_S_smb_neg
+
+
+end subroutine init_grad_smb
+
+
+!------------------------------------------------------------------------------
+!> Calcule le gradient vertical de smb
+
+subroutine grad_smb 
+
+  use module3D_phy
+  implicit none
+
+  do j = 1,ny
+     do i = 1,nx
+        if (Ylat(i,j).gt.77.) then                     ! region nord
+           if (bm_ref_t(i,j).lt. 0.) then              ! smb negatif
+              grad_bm (i,j) = grad_N_smb_neg
+
+           else if (bm_ref_t(i,j).ge. 0.) then         ! smb positif
+              grad_bm (i,j) = grad_N_smb_pos 
+
+           end if
+
+        else if (Ylat(i,j).le.77.) then                ! region sud
+           if (bm_ref_t(i,j).lt. 0.) then              ! smb negatif
+              grad_bm (i,j) = grad_S_smb_neg
+
+           else if (bm_ref_t(i,j).ge. 0.) then         ! smb positif
+              grad_bm (i,j) = grad_S_smb_pos 
+
+           end if
+        end if
+     end do
+  end do
+
+end subroutine grad_smb
+
+end module  climat_Grice2sea_years_perturb_mod
Index: /trunk/SOURCES/Makefile.grisli-gfortran.inc
===================================================================
--- /trunk/SOURCES/Makefile.grisli-gfortran.inc	(revision 41)
+++ /trunk/SOURCES/Makefile.grisli-gfortran.inc	(revision 42)
@@ -43,4 +43,5 @@
 mod_clim_tof = climat_forcage_mois_mod.o climat-forcage-insolation_mod_oneway.o \
 	climat-forcage-insolation_mod.o climat_GrIce2sea_years_mod.o \
+	climat_GrIce2sea_years_perturb_mod.o \
 	climat-perturb_mod-0.4.o \
 	ablation_mod.o no_ablation_mod.o 
@@ -628,4 +629,7 @@
 climat_GrIce2sea_years_mod.o : GrIce2sea_files/climat_GrIce2sea_years_mod.f90
 	$(F_NETCDF) GrIce2sea_files/climat_GrIce2sea_years_mod.f90
+
+climat_GrIce2sea_years_perturb_mod.o : GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90
+	$(F_NETCDF) GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90
 
 climat-forcage-stat-mois_mod-0.1.o : climat-forcage-stat-mois_mod-0.1.f90
Index: /trunk/SOURCES/Makefile.grisli.inc
===================================================================
--- /trunk/SOURCES/Makefile.grisli.inc	(revision 41)
+++ /trunk/SOURCES/Makefile.grisli.inc	(revision 42)
@@ -43,4 +43,5 @@
 mod_clim_tof = climat_forcage_mois_mod.o climat-forcage-insolation_mod_oneway.o \
 	climat-forcage-insolation_mod.o climat_GrIce2sea_years_mod.o \
+	climat_GrIce2sea_years_perturb_mod.o \
 	climat-perturb_mod-0.4.o \
 	ablation_mod.o no_ablation_mod.o 
@@ -632,4 +633,7 @@
 climat_GrIce2sea_years_mod.o : GrIce2sea_files/climat_GrIce2sea_years_mod.f90
 	$(F_NETCDF) GrIce2sea_files/climat_GrIce2sea_years_mod.f90
+
+climat_GrIce2sea_years_perturb_mod.o : GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90
+	$(F_NETCDF) GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90
 
 climat-forcage-stat-mois_mod-0.1.o : climat-forcage-stat-mois_mod-0.1.f90
