Index: /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90
===================================================================
--- /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90	(revision 163)
+++ /trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90	(revision 164)
@@ -62,8 +62,11 @@
 real,dimension(:),allocatable :: tpert          !> temperature for climate forcing
 real,dimension(:),allocatable :: spert          !> sea surface perturbation
+real,dimension(:),allocatable :: bpert          !> basal melt index perturbation
 real :: coefT                    !> pour modifier l'amplitude de la perturb. T
 character(len=120) :: filforc    !> nom du fichier forcage
 integer :: pertsmb               !> boolean: do we modify the smb?
 real    :: rapsmb                !> if we modify the smb this is the equivalent of rappact
+integer :: pertbmb               !> boolean: do we modify the bmelt under ice shelves?
+real    :: coefbmb               !> if we modify the bmelt, we read a perturbation (based on dD, insolation, etc.) and apply this coefficient
 
 
@@ -129,7 +132,7 @@
   !   bm(:,:) = bm(:,:) - 5.                     ! pour faire un masque a l'exterieur du Groenland actuel
   !end where
-  where ((H(:,:).lt.1.).and.(Bsoc(:,:).gt.0.).and.(bm(:,:).gt.0.))
-     bm(:,:) = -0.1
-  end where
+  !where ((H(:,:).lt.1.).and.(Bsoc(:,:).gt.0.).and.(bm(:,:).gt.0.))
+  !   bm(:,:) = -0.1
+  !end where
 !cdc test debug Hemin15 et Greeneem15
 !  where (bm(:,:).lt.-1000) bm(:,:)=0.
@@ -278,11 +281,24 @@
   end if
   
+  if (.not.allocated(bpert)) then
+     allocate(bpert(nft),stat=err)
+     if (err/=0) then
+        print *,"erreur a l'allocation du tableau Bpert",err
+        stop 4
+     end if
+  end if
+  
   do i=1,nft
-     read(num_forc,*) tdate(i),spert(i),tpert(i)
+     if (pertbmb==1) then
+        read(num_forc,*) tdate(i),spert(i),tpert(i),bpert(i)
+     else
+        read(num_forc,*) tdate(i),spert(i),tpert(i)
+        bpert(i)=0.
+     endif
   end do
   close(num_forc)
   
   tpert(:)=tpert(:)*coefT
-  
+  bpert(:)=1.+bpert(:)*coefbmb
 
   
@@ -325,5 +341,5 @@
   implicit none
   namelist/lapse_rates/T_lapse_rate
-  namelist/clim_pert_massb/coefT,filforc,pertsmb,rapsmb
+  namelist/clim_pert_massb/coefT,filforc,pertsmb,rapsmb,pertbmb,coefbmb
   
   rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
@@ -353,4 +369,6 @@
   write(num_rep_42,*) 'pertsmb        = ', pertsmb
   write(num_rep_42,*) 'rapsmb         = ', rapsmb
+  write(num_rep_42,*) 'pertbmb        = ', pertbmb
+  write(num_rep_42,*) 'coefbmb         = ', coefbmb
   write(num_rep_42,*)'/'
   write(num_rep_42,*)
@@ -388,4 +406,5 @@
            tafor=tpert(1)
            sealevel=spert(1)
+           coefbmshelf=bpert(1)
            ift=1
 
@@ -393,4 +412,5 @@
            tafor=tpert(nft)
            sealevel=spert(nft)
+           coefbmshelf=bpert(nft)
            ift=nft
            
@@ -402,4 +422,6 @@
                  sealevel=spert(i)+(spert(i+1)-spert(i))*    &
                       (time-tdate(i))/(tdate(i+1)-tdate(i))
+                 coefbmshelf=bpert(i)+(bpert(i+1)-bpert(i))*    &
+                      (time-tdate(i))/(tdate(i+1)-tdate(i))
                  ift=i
                  goto 100
@@ -409,5 +431,4 @@
 100     continue
 
-        
         Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor
         Ts(:,:)    = Tann(:,:)
@@ -420,9 +441,14 @@
            end if
         end if
+
+        !afq -- as in Pollard and Deconto, melt is more efficient at high temperature... ->
+        if (coefbmshelf.gt.1.) then
+           coefbmshelf = 1. + ( 1. - coefbmshelf ) * 5. !afq.... tuning for the degla?...
+        end if
         
-        ! coefmshelf est un coefficient qui fait varier bmgrz et bmshelf en fonction de tafor
-        coefbmshelf=(1.+tafor/10.)       ! coefbmshelf=0 pour tafor=-10deg
-        coefbmshelf=max(coefbmshelf,0.)  
-        coefbmshelf=min(coefbmshelf,2.)
+!afq --old:        ! coefmshelf est un coefficient qui fait varier bmgrz et bmshelf en fonction de tafor
+!afq --old:        coefbmshelf=(1.+tafor/10.)       ! coefbmshelf=0 pour tafor=-10deg
+!afq --old:        coefbmshelf=max(coefbmshelf,0.)  
+!afq --old:        coefbmshelf=min(coefbmshelf,2.)
         
 !        coefbmshelf=1.  !afq for calibration Greenland
