Index: trunk/SOURCES/3D-physique-gen_mod.f90
===================================================================
--- trunk/SOURCES/3D-physique-gen_mod.f90	(revision 236)
+++ trunk/SOURCES/3D-physique-gen_mod.f90	(revision 237)
@@ -141,5 +141,5 @@
   real ::  SECYEAR                      !< for relation an/seconds
   !real ::  S22                          ! ct for PDD calculation
-  real ::  sealevel                     !< niveau des mers 
+  real ::  sealevel                     !< niveau des mers, afq: now used as eustatic sea level
   ! real ::  SIGMA                        ! variabilite Tday
   real ::  SURF                         !< 
@@ -235,5 +235,4 @@
   integer,dimension(nx,ny) :: Mk_init      !< initial mask (with islands, outcrops, ...
   integer,dimension(nx,ny) :: MK           !< masks (ice sheet, max, above water, below water, 1)
-  integer,dimension(nx,ny) :: M_sealev     !< masks for sea level calculation (present ice sheet)
   integer,dimension(nx,ny) :: MNEG         !< masks (ice sheet, max, above water, below water, 1)
   integer,dimension(nx,ny) :: IBASE     !< type de base (froide, temperee)
@@ -245,6 +244,4 @@
   real,dimension(nx,ny) :: ACQUA        !< Surface des surfaces en eau
   real,dimension(nx,ny) :: B            !< Altitude de la base de la glace  'o'
-  real,dimension(nx,ny) :: B_sealev     !< Altitude de la base de la glace pour le
-  !< calcul du nivx de la mer (/present)
   real,dimension(nx,ny) :: BDOT         !< derivee de B / t
   real,dimension(nx,ny) :: B1           !< 
@@ -310,5 +307,4 @@
 
 
-  real,dimension(nx,ny) :: H_sealev     !< ice thickness for sea level calculation        ! (present ice sheet)
   real,dimension(nx,ny) :: HDOT         !< ice thickness derivee / t
   real,dimension(nx,ny) :: HDOTWATER
@@ -334,6 +330,5 @@
   real,dimension(nx,ny) :: SW           !< for bedrock isostasy
   real,dimension(nx,ny) :: S            !< altitude of ice sheet surface
-  real,dimension(nx,ny) :: S_sealev     !< altitude of ice sheet surface for sea
-  !< level calculation (present ice sheet)
+  real,dimension(nx,ny) :: sealevel_2d  !< local sea surface elevation
   real,dimension(nx,ny) :: SLOPE        !< slope 'o'
   real,dimension(nx,ny) :: SDX          !< slope derivee / x '>'
@@ -344,5 +339,5 @@
   real,dimension(nx,ny) :: SLOPE2my     !< = Sdy**2 + Sdxmy**2 '^'
   real,dimension(nx,ny) :: S0           !< altitude actuelle de la surface
-  real,dimension(nx,ny) :: slv          !< niveau de flottaison (sealevel et lakes)
+!afq -- replaced by sealevel_2D:  real,dimension(nx,ny) :: slv          !< niveau de flottaison (sealevel et lakes)
   real,dimension(nx,ny) :: TJULY        !< Ground air temperature July
   real,dimension(nx,ny) :: TANN         !< Ground air temperature annual
Index: trunk/SOURCES/Ant15_CISM_files/output_anta_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Ant15_CISM_files/output_anta_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/Ant15_CISM_files/output_anta_mod-0.4.f90	(revision 237)
@@ -83,8 +83,8 @@
 
 !        calcul de la hauteur au dessus de la flottaison
-         if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+         if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
                volf=volf+h(i,j)
          else
-            volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
+            volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j))
          endif
 
Index: trunk/SOURCES/Ant16_files/climat_InitMIP_years_perturb_mod.f90
===================================================================
--- trunk/SOURCES/Ant16_files/climat_InitMIP_years_perturb_mod.f90	(revision 236)
+++ trunk/SOURCES/Ant16_files/climat_InitMIP_years_perturb_mod.f90	(revision 237)
@@ -16,5 +16,5 @@
 
 
-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,coefbmshelf
+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,sealevel_2d,coefbmshelf
 use netcdf
 use io_netcdf_grisli
@@ -285,4 +285,6 @@
 100 continue
 
+  sealevel_2d(:,:) = sealevel
+
   Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor
   Ts(:,:)    = Tann(:,:)
Index: trunk/SOURCES/Ant40_files/lect-anteis_mod.f90
===================================================================
--- trunk/SOURCES/Ant40_files/lect-anteis_mod.f90	(revision 236)
+++ trunk/SOURCES/Ant40_files/lect-anteis_mod.f90	(revision 237)
@@ -228,5 +228,5 @@
     do J=1,NY
        do I=1,NX
-          if ((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL).LT.0.) then
+          if ((BSOC(I,J)+H(I,J)*RO/ROW -sealevel_2d(i,j)).LT.0.) then
              FLOT(I,J)=.TRUE.
           else
@@ -260,20 +260,4 @@
     ymax=ycc(nx,ny)/1000.
 
-!!$    !     lecture du fichier de reference pour le calcul du niveau des mers : etat actuel
-!!$    open(88,file=TRIM(DIRNAMEINP)//'grzone56-k000.pl')
-!!$
-!!$    read(88,*)
-!!$    read(88,*)
-!!$    read(88,*)
-!!$    do j=1,ny
-!!$       do i=1,nx
-!!$          read(88,*) S_sealev(i,j),H_sealev(i,j),B_sealev(i,j),M_sealev(i,j)
-!!$       enddo
-!!$    enddo
-!!$    close(88)
-!!$
-!!$    write(num_rep_42,*) 'fichier reference pour le niveau des mers : ', & 
-!!$         TRIM(DIRNAMEINP)//'grzone56-k000.pl'
-
     ! lecture du flux geothermique de Shapiro
     open(88,file=ghf_fich)
Index: trunk/SOURCES/Ant40_files/output_anta40_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Ant40_files/output_anta40_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/Ant40_files/output_anta40_mod-0.4.f90	(revision 237)
@@ -191,8 +191,8 @@
 
 !         calcul de la hauteur au dessus de la flottaison
-            if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+            if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
               volf=volf+h(i,j)*corrsurf(i,j)                          ! volume au-dessus de la flottaison
             else
-              volf=volf+(h(i,j)-row/ro*(sealevel-b(i,j)))*corrsurf(i,j) ! volume au-dessus de la flottaison
+              volf=volf+(h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)))*corrsurf(i,j) ! volume au-dessus de la flottaison
             endif
 
Index: trunk/SOURCES/Ant45_CISM_files/output_anta_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Ant45_CISM_files/output_anta_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/Ant45_CISM_files/output_anta_mod-0.4.f90	(revision 237)
@@ -84,8 +84,8 @@
 
 !        calcul de la hauteur au dessus de la flottaison
-         if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+         if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
                volf=volf+h(i,j)
          else
-            volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
+            volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j))
          endif
 
Index: trunk/SOURCES/Antarctique_general_files/output_anta_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Antarctique_general_files/output_anta_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/Antarctique_general_files/output_anta_mod-0.4.f90	(revision 237)
@@ -192,8 +192,8 @@
              
              !        calcul de la hauteur au dessus de la flottaison
-             if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+             if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
                 volf=volf+h(i,j)
              else
-                volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
+                volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j))
              endif
              
@@ -289,5 +289,5 @@
                 scal_np         =  scal_np + 1
                 scal_vol        =  scal_vol + H(i,j)
-                scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel-b(i,j))
+                scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel_2d(i,j)-b(i,j))
                 scal_meanhdot   =  scal_meanhdot + Hdot(i,j)
                 scal_sigma_hdot =  scal_sigma_hdot + Hdot(i,j)**2
Index: trunk/SOURCES/Draggings_modules/dragging_prescr_beta_buoyency_mod.f90
===================================================================
--- trunk/SOURCES/Draggings_modules/dragging_prescr_beta_buoyency_mod.f90	(revision 236)
+++ trunk/SOURCES/Draggings_modules/dragging_prescr_beta_buoyency_mod.f90	(revision 237)
@@ -500,9 +500,9 @@
 
 
-    where ( sealevel-B1(:,:).le.0.)             ! socle au dessus du niveau des mers
+    where ( sealevel_2d(:,:)-B1(:,:).le.0.)             ! socle au dessus du niveau des mers
        H_buoy(:,:) = H1(:,:)
 
     elsewhere
-       H_buoy(:,:) = H1(:,:)-row/ro*(sealevel-B1(:,:))
+       H_buoy(:,:) = H1(:,:)-row/ro*(sealevel_2d(:,:)-B1(:,:))
        
     end where
Index: trunk/SOURCES/Draggings_modules/dragging_prescr_beta_nolin_mod.f90
===================================================================
--- trunk/SOURCES/Draggings_modules/dragging_prescr_beta_nolin_mod.f90	(revision 236)
+++ trunk/SOURCES/Draggings_modules/dragging_prescr_beta_nolin_mod.f90	(revision 237)
@@ -606,9 +606,9 @@
 
 
-    where ( sealevel-B1(:,:).le.0.)             ! socle au dessus du niveau des mers
+    where ( sealevel_2d(:,:)-B1(:,:).le.0.)             ! socle au dessus du niveau des mers
        H_buoy(:,:) = H1(:,:)
 
     elsewhere
-       H_buoy(:,:) = H1(:,:)-row/ro*(sealevel-B1(:,:))
+       H_buoy(:,:) = H1(:,:)-row/ro*(sealevel_2d(:,:)-B1(:,:))
        
     end where
Index: trunk/SOURCES/Eurasie40_files/bmelt-eurasie-depth-lake_mod.f90
===================================================================
--- trunk/SOURCES/Eurasie40_files/bmelt-eurasie-depth-lake_mod.f90	(revision 236)
+++ trunk/SOURCES/Eurasie40_files/bmelt-eurasie-depth-lake_mod.f90	(revision 237)
@@ -85,5 +85,5 @@
          do i=1,nx
 
-            if ((sealevel-BSOC(i,j)).lt.450) then
+            if ((sealevel_2d(i,j)-BSOC(i,j)).lt.450) then
                     bmshelf(i,j)=2.00 !3.00! 2.55! 0.45!   0.65  
                     bmgrz(i,j)=2.00 !3.00! 2.55! 0.45 !  0.65
@@ -127,5 +127,5 @@
  
 
-                   if (slv(i,j).gt.0) then
+                   if (sealevel_2d(i,j).gt.0) then
                          bmshelf(i,j)=0.05!5.0    
                          bmgrz(i,j)=0.05!5.0
@@ -135,5 +135,5 @@
         
                      
-             if (slv(i,j).gt.sealevel) then
+             if (sealevel_2d(i,j).gt.sealevel) then
                      bmelt(i,j)=bmshelf(i,j)
                      if (fbm(i,j))  bmelt(i,j)=bmgrz(i,j)
@@ -166,5 +166,5 @@
 !   en fonction du nombre de points flottants
 
-           if (slv(i,j).gt.sealevel) then
+           if (sealevel_2d(i,j).gt.sealevel) then
                  bmelt(i,j)= ngr/4.*bmgrz(i,j)+(1.-ngr/4.)*bmelt(i,j)         
            else
Index: trunk/SOURCES/Eurasie40_files/climat-forcage-eurasie_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Eurasie40_files/climat-forcage-eurasie_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/Eurasie40_files/climat-forcage-eurasie_mod-0.4.f90	(revision 237)
@@ -493,5 +493,7 @@
    nom_table='lvflo'
    call printtable_r(levelflot,nom_table)
-    
+   
+   sealevel_2d(:,:) = sealevel
+
 !   coefbmshelf(:,:)=(1.+delTatime(:,:)/7.)  
 !coefbmshelf=1. ! modif pour test de sensibilite (tof 20 avril 01)
Index: trunk/SOURCES/Eurasie40_files/lakes-prescribed_mod-0.1.f90
===================================================================
--- trunk/SOURCES/Eurasie40_files/lakes-prescribed_mod-0.1.f90	(revision 236)
+++ trunk/SOURCES/Eurasie40_files/lakes-prescribed_mod-0.1.f90	(revision 237)
@@ -54,7 +54,7 @@
 subroutine lake_to_slv
 
-! inclut les lacs dans le tableau slv des niveaux de flottaison
+! inclut les lacs dans le tableau sealevel_2d des niveaux de flottaison
 
-slv(:,:)=sealevel
+sealevel_2d(:,:)=sealevel
 
 ! VOIR COMMENT ACTIVER-DESACTIVER LES LACS EN FONCTION DU TEMPS
@@ -64,5 +64,5 @@
       do i_lake=1,nb_lakes
          if (lake(i,j).eq.i_lake) then 
-            slv(i,j)=max(level_lakes(i_lake),Bsoc(i,j))
+            sealevel_2d(i,j)=max(level_lakes(i_lake),Bsoc(i,j))
          end if
       end do
Index: trunk/SOURCES/Eurasie40_files/output_eurasie40_mod-0.1.f90
===================================================================
--- trunk/SOURCES/Eurasie40_files/output_eurasie40_mod-0.1.f90	(revision 236)
+++ trunk/SOURCES/Eurasie40_files/output_eurasie40_mod-0.1.f90	(revision 237)
@@ -85,8 +85,8 @@
 
 !        calcul de la hauteur au dessus de la flottaison
-         if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+         if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
                volf=volf+h(i,j)
          else
-            volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
+            volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j))
          endif
 
Index: trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90
===================================================================
--- trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90	(revision 236)
+++ trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90	(revision 237)
@@ -16,5 +16,5 @@
 
 
-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,coefbmshelf
+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,sealevel_2d,coefbmshelf
 !use lect_climref_Ice2sea
 use netcdf
@@ -382,4 +382,6 @@
         endif
 100     continue
+
+        sealevel_2d(:,:) = sealevel
 
         Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor
Index: trunk/SOURCES/GrIce2sea_files/lect_GrIce2sea_gen_nc.f90
===================================================================
--- trunk/SOURCES/GrIce2sea_files/lect_GrIce2sea_gen_nc.f90	(revision 236)
+++ trunk/SOURCES/GrIce2sea_files/lect_GrIce2sea_gen_nc.f90	(revision 237)
@@ -207,11 +207,4 @@
 
 
-    !     lecture du fichier de reference pour le calcul du niveau des mers : etat actuel
-    !     pas de version correcte pour l'instant -> valeurs initiales
-    S_sealev(:,:) = S0(:,:)
-    H_sealev(:,:) = H0(:,:)
-    B_sealev(:,:) = Bsoc(:,:)
-    M_sealev(:,:) = Mk0(:,:)
-
     if (itracebug.eq.1)  call tracebug(' apres lect entete')
   
Index: trunk/SOURCES/GrIce2sea_files/output_Grice2sea_mod.f90
===================================================================
--- trunk/SOURCES/GrIce2sea_files/output_Grice2sea_mod.f90	(revision 236)
+++ trunk/SOURCES/GrIce2sea_files/output_Grice2sea_mod.f90	(revision 237)
@@ -196,8 +196,8 @@
              
              !        calcul de la hauteur au dessus de la flottaison
-             if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+             if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
                 volf=volf+h(i,j)
              else
-                volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
+                volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j))
              endif
              
@@ -298,5 +298,5 @@
                 scal_np         =  scal_np + 1
                 scal_vol        =  scal_vol + H(i,j)
-                scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel-b(i,j))
+                scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel_2d(i,j)-b(i,j))
                 scal_meanhdot   =  scal_meanhdot + Hdot(i,j)
                 scal_sigma_hdot =  scal_sigma_hdot + Hdot(i,j)**2
Index: trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mar_mod.f90
===================================================================
--- trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mar_mod.f90	(revision 236)
+++ trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mar_mod.f90	(revision 237)
@@ -55,5 +55,5 @@
    do i=1,nx
       read (20,*) S0CLIM(i,j)
-      S0CLIM(i,j)=max(S0CLIM(i,j),slv(i,j))   ! pour etre au niveau des mers
+      S0CLIM(i,j)=max(S0CLIM(i,j),sealevel_2d(i,j))   ! pour etre au niveau des mers
    end do
 end do
@@ -146,5 +146,5 @@
   do j=1,ny
      do i=1,nx
-        alti=max(slv(i,j),S(I,J))
+        alti=max(sealevel_2d(i,j),S(I,J))
         Tann(i,j)=TA0(i,j)+lapse_ann*(alti-S0CLIM(i,j))
         Tjuly(i,j)=TJ0(i,j)+lapse_july*(alti-S0CLIM(i,j))
Index: trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_lapsecouche_mod.f90
===================================================================
--- trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_lapsecouche_mod.f90	(revision 236)
+++ trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_lapsecouche_mod.f90	(revision 237)
@@ -114,5 +114,5 @@
             Tm_0(i,j,mo)=Tm(i,j,mo)      ! temp sur la calotte du GCM a l'instant initial
             
-            zel=max(slv(i,j),S0CLIM(i,j))
+            zel=max(sealevel_2d(i,j),S0CLIM(i,j))
 ! on calcule une temperature au sol a partir de la temperature sur la calotte de ref du GCM
 
@@ -131,5 +131,5 @@
             end if
 
-            !zel=max(slv(i,j),S(i,j))    
+            !zel=max(sealevel_2d(i,j),S(i,j))    
             zel=max(0.,S(i,j))    
 ! puis on calcule la temperature sur la calotte S
Index: trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_mod.f90
===================================================================
--- trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_mod.f90	(revision 236)
+++ trunk/SOURCES/Greeneem_files/lect-clim-act-greeneem_mois_mod.f90	(revision 237)
@@ -74,5 +74,5 @@
       do i=1,nx
          read (20,*) S0CLIM(i,j)
-         S0CLIM(i,j)=max(S0CLIM(i,j),slv(i,j))   ! pour etre au niveau des mers
+         S0CLIM(i,j)=max(S0CLIM(i,j),sealevel_2d(i,j))   ! pour etre au niveau des mers
       end do
    end do
@@ -120,5 +120,5 @@
       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),S0CLIM(I,J)) !slv=sealev pour nolake
+         ZS(i,j)=max(sealevel_2d(i,j),S0CLIM(I,J)) !sealevel_2d=sealev pour nolake
 
          do mo=1,mois
@@ -261,5 +261,5 @@
 do j=1,ny
    do i=1,nx
-      ZS(i,j)=max(slv(i,j),S(i,j))
+      ZS(i,j)=max(sealevel_2d(i,j),S(i,j))
       Tann(i,j)=d_ann+gamma_ann*S(i,j)+c_ann*ylat(i,j)+kappa_ann*(360-xlong(i,j))
       Tjuly(i,j)=d_july+gamma_july*S(i,j)+c_july*ylat(i,j)+kappa_july*(360-xlong(i,j))
Index: trunk/SOURCES/Greeneem_files/lect-greeneem_mod.f90
===================================================================
--- trunk/SOURCES/Greeneem_files/lect-greeneem_mod.f90	(revision 236)
+++ trunk/SOURCES/Greeneem_files/lect-greeneem_mod.f90	(revision 237)
@@ -86,5 +86,5 @@
        do i=1,nx
 
-          archimed = Bsoc(i,j)+H(i,j)*ro/row -sealevel
+          archimed = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j)
 
           if ((archimed.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte
Index: trunk/SOURCES/Greeneem_files/output_greeneem_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Greeneem_files/output_greeneem_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/Greeneem_files/output_greeneem_mod-0.4.f90	(revision 237)
@@ -141,8 +141,8 @@
      do j=1,ny
         do i=1,nx
-           if (mask_cal(i,j,k).and..not.flot(i,j).and.H(i,j).GT.1..and.((sealevel-B(i,j)).LE.0.)) then
+           if (mask_cal(i,j,k).and..not.flot(i,j).and.H(i,j).GT.1..and.((sealevel_2d(i,j)-B(i,j)).LE.0.)) then
               isvolf(kk) = isvolf(kk) + H(i,j)
            else
-              isvolf(kk) = isvolf(kk) + H(i,j)-row/ro*(sealevel-B(i,j))
+              isvolf(kk) = isvolf(kk) + H(i,j)-row/ro*(sealevel_2d(i,j)-B(i,j))
            endif
         enddo
Index: trunk/SOURCES/Greenmint40_files/lect-greenmint_mod.f90
===================================================================
--- trunk/SOURCES/Greenmint40_files/lect-greenmint_mod.f90	(revision 236)
+++ trunk/SOURCES/Greenmint40_files/lect-greenmint_mod.f90	(revision 237)
@@ -112,5 +112,5 @@
    do i=1,nx
 
-      archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
+      archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j)
 
       if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte
Index: trunk/SOURCES/Greenmint40_files/output_green_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Greenmint40_files/output_green_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/Greenmint40_files/output_green_mod-0.4.f90	(revision 237)
@@ -89,8 +89,8 @@
 !        calcul de la hauteur au dessus de la flottaison
 
-         if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+         if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
                volf=volf+h(i,j)
          else
-            volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
+            volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j))
          endif
 
Index: trunk/SOURCES/Heino_files/flottab2-0.5-heino.f90
===================================================================
--- trunk/SOURCES/Heino_files/flottab2-0.5-heino.f90	(revision 236)
+++ trunk/SOURCES/Heino_files/flottab2-0.5-heino.f90	(revision 237)
@@ -135,5 +135,5 @@
       uys1(i,j)=uybar(i,j)
 
-      archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
+      archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j)
  
 arch: if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte
@@ -146,5 +146,5 @@
              if (igrdline.eq.1) then   ! en cas de grounding line prescrite
                 flot(i,j)=.false.
-                H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10
+                H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10
                 new_flot_point(i,j)=.false.
              endif
@@ -162,5 +162,5 @@
                 if (igrdline.eq.1) then   ! en cas de grounding line prescrite
                     flot(i,j)=.false.
-                    H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10
+                    H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10
                    new_flot_point(i,j)=.false.
                endif
@@ -184,5 +184,5 @@
 
                 surnet=H(i,j)*(1.-ro/row)
-                S(i,j)=surnet+sealevel
+                S(i,j)=surnet+sealevel_2d(i,j)
                 B(i,j)=S(i,j)-H(i,j)
              end if
@@ -234,5 +234,5 @@
 !                 sur le demi noeud condition archim > 100 m
 
-             archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel+ro/row*Hmx(i,j)
+             archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel_2d(i,j)+ro/row*Hmx(i,j)
              gzmx(i,j)=gzmx(i,j).and.(archim.le.100.) 
 
@@ -301,5 +301,5 @@
 !                 sur le demi noeud condition archim > 100 m
 
-             archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel+ro/row*Hmy(i,j)
+             archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel_2d(i,j)+ro/row*Hmy(i,j)
              gzmy(i,j)=gzmy(i,j).and.(archim.le.100.) 
 
@@ -554,5 +554,5 @@
 !!!!!!!!!                      on teste si la tache n'est pas completement ice stream
       if (icetrim(table_out(i,j))) THEN ! on a une ile d'ice stream   !!!! TEST2
-              smax_i=0 ; smax_j=0 ; smax_=sealevel
+              smax_i=0 ; smax_j=0 ; smax_=sealevel !afq -- warning: what about local sealevel?
               DO II=3,NX-2
                 DO JJ=3,NY-2
Index: trunk/SOURCES/Hemin40_files/bmelt-hemin40-depth_mod.f90
===================================================================
--- trunk/SOURCES/Hemin40_files/bmelt-hemin40-depth_mod.f90	(revision 236)
+++ trunk/SOURCES/Hemin40_files/bmelt-hemin40-depth_mod.f90	(revision 237)
@@ -85,5 +85,5 @@
          do i=1,nx
 
-            if ((sealevel-bsoc(i,j)).lt.600) then
+            if ((sealevel_2d(i,j)-bsoc(i,j)).lt.600) then
                     bmshelf(i,j)=0.2 ! 0.45!   0.65  
                     bmgrz(i,j)=0.2 ! 0.45 !  0.65
@@ -129,5 +129,5 @@
 shelf:    if (flot(i,j)) then    ! partie flottante
 
-            if ((sealevel-bsoc(i,j)).lt.600) then
+            if ((sealevel_2d(i,j)-bsoc(i,j)).lt.600) then
                     if ((j.lt.75).and.(i.lt.145) ) then  !Atlantique 
                     bmelt(i,j)=(.6+.4*coefbmshelf)*bmshelf(i,j)
Index: trunk/SOURCES/Hudson_files/bmelt_hudson_mod.f90
===================================================================
--- trunk/SOURCES/Hudson_files/bmelt_hudson_mod.f90	(revision 236)
+++ trunk/SOURCES/Hudson_files/bmelt_hudson_mod.f90	(revision 237)
@@ -52,5 +52,5 @@
 !             bmgrz(i,j)= 10.
 
-          bmshelf(i,j)=(sealevel-BSOC(i,j))/200  ! jalv: eq lineaire en fonction de la profondeur de la mer
+          bmshelf(i,j)=(sealevel_2d(i,j)-BSOC(i,j))/200  ! jalv: eq lineaire en fonction de la profondeur de la mer
           bmgrz(i,j)=bmshelf(i,j)
 !          endif
Index: trunk/SOURCES/Hudson_files/climat-hudson_mod.f90
===================================================================
--- trunk/SOURCES/Hudson_files/climat-hudson_mod.f90	(revision 236)
+++ trunk/SOURCES/Hudson_files/climat-hudson_mod.f90	(revision 237)
@@ -102,11 +102,11 @@
 !jalv: test pour que le bilan de masse soit 0 sur l'ocean plus profond de 600 m:
 !      c'est le bon critere??:
-where (sealevel-BSOC(:,:).gt.600) acc(:,:)=0.
-where (sealevel-BSOC(:,:).gt.600) abl(:,:)=0.
+where (sealevel_2d(:,:)-BSOC(:,:).gt.600) acc(:,:)=0.
+where (sealevel_2d(:,:)-BSOC(:,:).gt.600) abl(:,:)=0.
 
-where ((sealevel-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2500.)) abl(:,:)=-1.
-where ((sealevel-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2600.)) abl(:,:)=-5.
-where ((sealevel-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2650.)) abl(:,:)=-10.
-where ((sealevel-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2700.)) abl(:,:)=-20.
+where ((sealevel_2d(:,:)-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2500.)) abl(:,:)=-1.
+where ((sealevel_2d(:,:)-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2600.)) abl(:,:)=-5.
+where ((sealevel_2d(:,:)-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2650.)) abl(:,:)=-10.
+where ((sealevel_2d(:,:)-BSOC(:,:).gt.100).and.(dist_heino(:,:).gt.2700.)) abl(:,:)=-20.
 
 bm(:,:)=acc(:,:)+abl(:,:)
@@ -135,5 +135,5 @@
 !      c'est le bon critere??:
 !where (BSOC(:,:).lt.B(:,:)) Tann(:,:)=0.
-where (sealevel-BSOC(:,:).gt.600) Tann(:,:)=0. 
+where (sealevel_2d(:,:)-BSOC(:,:).gt.600) Tann(:,:)=0. 
 
 
Index: trunk/SOURCES/Hudson_files/eaubasale-0.5_hudson_mod.f90
===================================================================
--- trunk/SOURCES/Hudson_files/eaubasale-0.5_hudson_mod.f90	(revision 236)
+++ trunk/SOURCES/Hudson_files/eaubasale-0.5_hudson_mod.f90	(revision 237)
@@ -139,5 +139,5 @@
         if (flot(i,j))then  ! points flottants
            klimit(i,j)=1
-           limit_hw(i,j)=(sealevel-Bsoc(i,j))*rowg/rofreshg
+           limit_hw(i,j)=(sealevel_2d(i,j)-Bsoc(i,j))*rowg/rofreshg
 
         else if (IBASE(I,J).eq.1) then ! base froide
@@ -182,5 +182,5 @@
            pot_w(I,J)=pot_w(I,J)+rog*H(I,J)
 
-           pot_f(I,J)=rowg*(sealevel-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf
+           pot_f(I,J)=rowg*(sealevel_2d(i,j)-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf
         enddo
      enddo
@@ -308,5 +308,5 @@
 
               if (flot(i,j)) then  ! if flot > hwater=hwater in ocean 
-                 hwater(i,j)= sealevel - bsoc(i,j)
+                 hwater(i,j)= sealevel_2d(i,j) - bsoc(i,j)
                  !             hwater(i,j)= max(0.,hwater(i,j))
                  if (hwater(i,j).lt.0.) hwater(i,j)=0.
@@ -320,5 +320,5 @@
 
               ! Attention le bloc suivant est pour le run 20 
-              !           Zflot=row/ro*(sealevel-Bsoc(i,j))-10.
+              !           Zflot=row/ro*(sealevel_2d(i,j)-Bsoc(i,j))-10.
               Zflot=H(i,j)*rog/rofreshg
 
Index: trunk/SOURCES/MISMIP3D_files/lect-mismip3d_mod.f90
===================================================================
--- trunk/SOURCES/MISMIP3D_files/lect-mismip3d_mod.f90	(revision 236)
+++ trunk/SOURCES/MISMIP3D_files/lect-mismip3d_mod.f90	(revision 237)
@@ -111,5 +111,5 @@
   do j=1,ny
      do i=1,nx
-        if ((bsoc(i,j)+h(i,j)*ro/row -sealevel).lt.0.) then
+        if ((bsoc(i,j)+h(i,j)*ro/row -sealevel_2d(i,j)).lt.0.) then
            flot(i,j)=.true.
         else
Index: trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90
===================================================================
--- trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90	(revision 236)
+++ trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90	(revision 237)
@@ -483,5 +483,5 @@
 
          case(-2)                             ! bord Est
-            call vel_U_East(pvi(i,j),pvi(i-1,j),rog,H(i-1,j),rowg,slv(i-1,j),B(i-1,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
+            call vel_U_East(pvi(i,j),pvi(i-1,j),rog,H(i-1,j),rowg,sealevel_2d(i-1,j),B(i-1,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
        
          case(-3)                             ! bord nord
@@ -489,5 +489,5 @@
 
          case(-4)                             ! bord West
-            call vel_U_West(pvi(i,j),pvi(i-1,j),rog,H(i,j),rowg,slv(i,j),B(i,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
+            call vel_U_West(pvi(i,j),pvi(i-1,j),rog,H(i,j),rowg,sealevel_2d(i,j),B(i,j),dx,Tu(i,j,:,:),Tv(i,j,:,:),opposx(i,j))
 
          case(-12)                            ! coin SE
@@ -527,5 +527,5 @@
             
          case(-1)                             ! bord sud
-            call vel_V_South(pvi(i,j),pvi(i,j-1),rog,H(i,j),rowg,slv(i,j),B(i,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
+            call vel_V_South(pvi(i,j),pvi(i,j-1),rog,H(i,j),rowg,sealevel_2d(i,j),B(i,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
             
          case(-2)                             ! bord Est
@@ -533,5 +533,5 @@
             
          case(-3)                             ! bord nord
-            call vel_V_North(pvi(i,j),pvi(i,j-1),rog,H(i-1,j),rowg,slv(i-1,j),B(i-1,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
+            call vel_V_North(pvi(i,j),pvi(i,j-1),rog,H(i-1,j),rowg,sealevel_2d(i-1,j),B(i-1,j),dx,Sv(i,j,:,:),Su(i,j,:,:),opposy(i,j))
             
          case(-4)                             ! bord West
Index: trunk/SOURCES/Snowball_files/bmelt-snowball-depth_mod.f90
===================================================================
--- trunk/SOURCES/Snowball_files/bmelt-snowball-depth_mod.f90	(revision 236)
+++ trunk/SOURCES/Snowball_files/bmelt-snowball-depth_mod.f90	(revision 237)
@@ -66,5 +66,5 @@
          do i=1,nx
 
-            if ((sealevel-bsoc(i,j)).lt.600) then
+            if ((sealevel_2d(i,j)-bsoc(i,j)).lt.600) then
                     bmshelf(i,j)=0.2 ! 0.45!   0.65  
                     bmgrz(i,j)=0.2 ! 0.45 !  0.65
@@ -106,5 +106,5 @@
 shelf:    if (flot(i,j)) then    ! partie flottante
 
-            if ((sealevel-bsoc(i,j)).lt.600) then
+            if ((sealevel_2d(i,j)-bsoc(i,j)).lt.600) then
                     if ((j.lt.75).and.(i.lt.145) ) then  !Atlantique 
                     bmelt(i,j)=(.6+.4*coefbmshelf)*bmshelf(i,j)
Index: trunk/SOURCES/Snowball_files/output_snowball_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Snowball_files/output_snowball_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/Snowball_files/output_snowball_mod-0.4.f90	(revision 237)
@@ -83,8 +83,8 @@
                   vol=vol+h(i,j) 
 !        calcul de la hauteur au dessus de la flottaison
-                  if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
+                  if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
                      volf=volf+h(i,j)
                   else
-                     volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
+                     volf=volf+h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j))
                   endif
 
Index: trunk/SOURCES/bilan_eau_mod.f90
===================================================================
--- trunk/SOURCES/bilan_eau_mod.f90	(revision 236)
+++ trunk/SOURCES/bilan_eau_mod.f90	(revision 237)
@@ -45,15 +45,15 @@
 subroutine init_bilan_eau
  ! initialisation des variables
-	diff_H=0.
-	sum_H_old = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)
-	tot_water(:,:)=0.
-	bm_dt(:,:)=0.
-	bmelt_dt(:,:)=0.
+        diff_H=0.
+        sum_H_old = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)
+        tot_water(:,:)=0.
+        bm_dt(:,:)=0.
+        bmelt_dt(:,:)=0.
         alpha_flot=ro/row
-        archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel
+        archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:)
         gr_line(:,:)=0
         do j=1,ny
            do i=1,nx
-              !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded with ice
+              !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice
               if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice
                  if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1
@@ -67,5 +67,5 @@
 end subroutine init_bilan_eau
 
-	
+
 
 
@@ -101,13 +101,13 @@
 
 where (ice(2:nx-1,2:ny-1).eq.1)
-	Bm_dtt(2:nx-1,2:ny-1) = Bm_dtt(2:nx-1,2:ny-1) + Bm_dt(2:nx-1,2:ny-1) !* dt ! somme Bm sur dt
-	bmelt_dtt(2:nx-1,2:ny-1) = bmelt_dtt(2:nx-1,2:ny-1) + bmelt_dt(2:nx-1,2:ny-1) ! * dt ! somme bmelt sur dt
+        Bm_dtt(2:nx-1,2:ny-1) = Bm_dtt(2:nx-1,2:ny-1) + Bm_dt(2:nx-1,2:ny-1) !* dt ! somme Bm sur dt
+        bmelt_dtt(2:nx-1,2:ny-1) = bmelt_dtt(2:nx-1,2:ny-1) + bmelt_dt(2:nx-1,2:ny-1) ! * dt ! somme bmelt sur dt
 endwhere
 
-archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel
+archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:)
 gr_line(:,:)=0
 do j=1,ny
   do i=1,nx
-     !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded with ice
+     !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice
      if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice
       if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1
@@ -141,10 +141,10 @@
    
 ! bilan d'eau sur la grille :
-	water_bilan=sum(tot_water(:,:))
-	diff_H = diff_H/dtt
+        water_bilan=sum(tot_water(:,:))
+        diff_H = diff_H/dtt
 
 !999 format(f0.2,1x,e15.8,1x,i10,8(1x,e15.8))
 !	write(6,999),time,sum_H,count(ice(:,:)==1),diff_H,water_bilan,sum(calv_dtt(:,:))/dtt,sum(ablbord_dtt(:,:))/dtt,sum(bmelt_dtt(:,:),mask=ice(:,:)==1)/dtt,sum(bm(:,:),mask=ice(:,:)==1),sum(Bm_dtt(:,:))/dtt,sum(bmelt_dtt(:,:))/dtt
-	diff_H_water_bilan(2:nx-1,2:ny-1)=tot_water(2:nx-1,2:ny-1)-diff_H_2D(2:nx-1,2:ny-1)
+        diff_H_water_bilan(2:nx-1,2:ny-1)=tot_water(2:nx-1,2:ny-1)-diff_H_2D(2:nx-1,2:ny-1)
 
 endif
Index: trunk/SOURCES/calving_frange.f90
===================================================================
--- trunk/SOURCES/calving_frange.f90	(revision 236)
+++ trunk/SOURCES/calving_frange.f90	(revision 237)
@@ -89,5 +89,5 @@
 
     Hcoup(:,:) = min ( max(                                         &
-       (-(Bsoc0(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    &
+       (-(Bsoc0(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    &
        *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               &
        , hcoup_plateau), hcoup_abysses )
@@ -121,5 +121,5 @@
 
     Hcoup(:,:) = min ( max(                                         &
-       (-(Bsoc(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    &
+       (-(Bsoc(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    &
        *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               &
        , hcoup_plateau), hcoup_abysses )
@@ -409,7 +409,7 @@
                       calv(i,j)=-h(i,j) 
                       !cdc                     H(i,j)=1. 
-                      !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01))
+                      !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel_2d(i,j) - Bsoc(i,j))*row/ro-0.01))
                       H(i,j)=0.
-                      S(i,j)=H(i,j)*(1.-ro/row) + sealevel
+                      S(i,j)=H(i,j)*(1.-ro/row) + sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
                       B(i,j)=S(i,j) - H(i,j)
                       ! ATTENTION ne pas mettre ice=0 sinon degradation bilan d'eau (bm et bmelt non comptabilises dans ce cas)
@@ -439,5 +439,5 @@
         ice(:,:)=0
         H(:,:)=0.
-        S(:,:)=H(:,:)*(1.-ro/row) + sealevel
+        S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
         B(:,:)=S(:,:) - H(:,:)
     endwhere
Index: trunk/SOURCES/calving_frange_abuk.f90
===================================================================
--- trunk/SOURCES/calving_frange_abuk.f90	(revision 236)
+++ trunk/SOURCES/calving_frange_abuk.f90	(revision 237)
@@ -96,5 +96,5 @@
 
     Hcoup(:,:) = min ( max(                                         &
-       (-(Bsoc0(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    &
+       (-(Bsoc0(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    &
        *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               &
        , hcoup_plateau), hcoup_abysses )
@@ -128,5 +128,5 @@
 
     Hcoup(:,:) = min ( max(                                         &
-       (-(Bsoc(:,:)-sealevel) - prof_plateau)/(prof_abysses-prof_plateau)    &
+       (-(Bsoc(:,:)-sealevel_2d(:,:)) - prof_plateau)/(prof_abysses-prof_plateau)    &
        *(hcoup_abysses-hcoup_plateau)+hcoup_plateau               &
        , hcoup_plateau), hcoup_abysses )
@@ -416,7 +416,7 @@
                       calv(i,j)=-h(i,j) 
                       !cdc                     H(i,j)=1. 
-                      !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel - Bsoc(i,j))*row/ro-0.01))
+                      !cdc 1m                     H(i,j)=min(1.,max(0.,(sealevel_2d(i,j) - Bsoc(i,j))*row/ro-0.01))
                       H(i,j)=0.
-                      S(i,j)=H(i,j)*(1.-ro/row) + sealevel
+                      S(i,j)=H(i,j)*(1.-ro/row) + sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
                       B(i,j)=S(i,j) - H(i,j)
                       ! ATTENTION ne pas mettre ice=0 sinon degradation bilan d'eau (bm et bmelt non comptabilises dans ce cas)
@@ -446,5 +446,5 @@
 !~         ice(:,:)=0
 !~         H(:,:)=0.
-!~         S(:,:)=H(:,:)*(1.-ro/row) + sealevel
+!~         S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:)
 !~         B(:,:)=S(:,:) - H(:,:)
 !~     endwhere
@@ -455,5 +455,5 @@
         ice(:,:)=0
         H(:,:)=0.
-        S(:,:)=H(:,:)*(1.-ro/row) + sealevel
+        S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
         B(:,:)=S(:,:) - H(:,:)
     endwhere    
Index: trunk/SOURCES/climat-forcage-insolation_mod.f90
===================================================================
--- trunk/SOURCES/climat-forcage-insolation_mod.f90	(revision 236)
+++ trunk/SOURCES/climat-forcage-insolation_mod.f90	(revision 237)
@@ -12,5 +12,5 @@
 ! C. Dumas 07/2015
 
-use module3D_phy,only:nx,ny,S,H,flot,slv,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp,time,pi,dx,xlong,ylat,bsoc		
+use module3D_phy,only:nx,ny,S,H,flot,Tann,Tjuly,Tmois,acc,coefbmshelf,sealevel_2d,ro,num_param,num_rep_42,dirnameinp,time,pi,dx,xlong,ylat,bsoc		
 use netcdf
 use io_netcdf_grisli
@@ -637,5 +637,5 @@
 do j=1,ny
    do i=1,nx
-      Zs(i,j)=max(slv(i,j),S(i,j))
+      Zs(i,j)=max(sealevel_2d(i,j),S(i,j))
       do ictr=CO2SnapMin,CO2SnapMax,1
          do igtr=1,gtr
Index: trunk/SOURCES/climat-forcage-insolation_mod_oneway.f90
===================================================================
--- trunk/SOURCES/climat-forcage-insolation_mod_oneway.f90	(revision 236)
+++ trunk/SOURCES/climat-forcage-insolation_mod_oneway.f90	(revision 237)
@@ -12,5 +12,5 @@
 ! C. Dumas 06/2015
 
-USE module3D_phy,only:nx,ny,S,slv,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp,time
+USE module3D_phy,only:nx,ny,S,sealevel_2d,Tann,Tjuly,Tmois,acc,coefbmshelf,ro,num_param,num_rep_42,dirnameinp,time
 !use interface_input
 use netcdf
@@ -245,5 +245,5 @@
 do j=1,ny
    do i=1,nx
-      Zs(i,j)=max(slv(i,j),S(i,j))
+      Zs(i,j)=max(sealevel_2d(i,j),S(i,j))
       !Il faut mettre S0(i,j) si pas d'interpolation sinon Ssnap(i,j,ictr,igtr)
       !if (Zs(i,j).ge.Ssnap(i,j,ictr,igtr)) then
Index: trunk/SOURCES/climat-perturb_mod-0.4.f90
===================================================================
--- trunk/SOURCES/climat-perturb_mod-0.4.f90	(revision 236)
+++ trunk/SOURCES/climat-perturb_mod-0.4.f90	(revision 237)
@@ -13,5 +13,5 @@
 
 
-use module3d_phy,only:nx,ny,S,S0,Tann,Tjuly,precip,acc,Ylat,num_forc,num_param,num_rep_42,dirforcage,dirnameinp,tafor,time,sealevel,coefbmshelf
+use module3d_phy,only:nx,ny,S,S0,Tann,Tjuly,precip,acc,Ylat,num_forc,num_param,num_rep_42,dirforcage,dirnameinp,tafor,time,sealevel,sealevel_2d,coefbmshelf
 use netcdf
 use io_netcdf_grisli
@@ -307,4 +307,5 @@
 100 continue
 
+  sealevel_2d(:,:)=sealevel
 
   !  coefmshelf est un coefficient qui fait varier bmgrz et bmshelf
Index: trunk/SOURCES/climat_forcage_mod.f90
===================================================================
--- trunk/SOURCES/climat_forcage_mod.f90	(revision 236)
+++ trunk/SOURCES/climat_forcage_mod.f90	(revision 237)
@@ -6,5 +6,5 @@
 ! possibilite d'utiliser autant de snapshots que l'on veut (ntr)
 ! fichiers de forcage, Snapshot et valeurs de lapse rate definis dans fichier param
-  use module3d_phy,only:nx,ny,sealevel,slv,S,S0,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,coefbmshelf,PI,ro,time,dirforcage,dirnameinp
+  use module3d_phy,only:nx,ny,sealevel,sealevel_2d,S,S0,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,coefbmshelf,PI,ro,time,dirforcage,dirnameinp
   use netcdf
   use io_netcdf_grisli
@@ -279,5 +279,5 @@
 100 continue
 
-  slv(:,:)=sealevel
+  sealevel_2d(:,:)=sealevel
 
   if(timeloc.le.ttr(1)) then   ! time en dehors des limites du fichier forcage
@@ -333,5 +333,5 @@
 !     calcul de Tjuly et et Tann
   !$OMP WORKSHARE
-  Zs(:,:)=max(slv(:,:),S(:,:))
+  Zs(:,:)=max(sealevel_2d(:,:),S(:,:))
   deltaZs(:,:)= -lapserate*(Zs(:,:)-S0(:,:)) ! difference de temperature entre Zs et S0
   !$OMP END WORKSHARE
Index: trunk/SOURCES/eaubasale-0.5_mod.f90
===================================================================
--- trunk/SOURCES/eaubasale-0.5_mod.f90	(revision 236)
+++ trunk/SOURCES/eaubasale-0.5_mod.f90	(revision 237)
@@ -158,5 +158,5 @@
         if (flot(i,j))then  ! points flottants
            klimit(i,j)=1
-           limit_hw(i,j)=(sealevel-Bsoc(i,j))*rowg/rofreshg
+           limit_hw(i,j)=(sealevel_2d(i,j)-Bsoc(i,j))*rowg/rofreshg
 
         else if (IBASE(I,J).eq.1) then ! base froide
@@ -200,5 +200,5 @@
            pot_w(I,J)=pot_w(I,J)+rog*H(I,J)
 
-           pot_f(I,J)=rowg*(sealevel-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf
+           pot_f(I,J)=rowg*(sealevel_2d(i,j)-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf
         enddo
      enddo
@@ -314,5 +314,5 @@
 
               if (flot(i,j)) then  ! if flot > hwater=hwater in ocean 
-                 hwater(i,j)= sealevel - bsoc(i,j)
+                 hwater(i,j)= sealevel_2d(i,j) - bsoc(i,j)
                  !             hwater(i,j)= max(0.,hwater(i,j))
                  if (hwater(i,j).lt.0.) hwater(i,j)=0.
Index: trunk/SOURCES/flottab2-0.7.f90
===================================================================
--- trunk/SOURCES/flottab2-0.7.f90	(revision 236)
+++ trunk/SOURCES/flottab2-0.7.f90	(revision 237)
@@ -122,5 +122,5 @@
        where ( mk_init(:,:).eq.1)                 ! pose
           flot(:,:) = .False.
-          H(:,:)=max(H(:,:),(10.+sealevel-Bsoc(:,:))*row/ro)  ! pour avoir archim=10
+          H(:,:)=max(H(:,:),(10.+sealevel_2d(:,:)-Bsoc(:,:))*row/ro)  ! pour avoir archim=10
           S(:,:) = Bsoc(:,:) + H(:,:) 
        elsewhere
@@ -175,5 +175,5 @@
           uys1(i,j)=uybar(i,j)
 
-          archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
+          archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel_2d(i,j)
           !      if ((i.eq.132).and.(j.eq.183)) print*,'archim=',archim
 
@@ -189,5 +189,5 @@
                 if (igrdline.eq.1) then   ! en cas de grounding line prescrite
                    flot(i,j)=.false.
-                   H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10
+                   H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10
                    new_flot_point(i,j)=.false.
                 endif
@@ -202,5 +202,5 @@
                    if (igrdline.eq.1) then   ! en cas de grounding line prescrite
                       flot(i,j)=.false.
-                      H(i,j)=(10.+sealevel-Bsoc(i,j))*row/ro  ! pour avoir archim=10
+                      H(i,j)=(10.+sealevel_2d(i,j)-Bsoc(i,j))*row/ro  ! pour avoir archim=10
                       new_flot_point(i,j)=.false.
                    endif
@@ -228,5 +228,5 @@
              flot(i,j)=.true.  !cdc points ocean sont flot meme sans glace
              H(i,j)=0.
-             S(i,j)=sealevel
+             S(i,j)=sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
              B(i,j)=S(i,j)-H(i,j)
           endif arch
@@ -239,5 +239,5 @@
 
              surnet=H(i,j)*(1.-ro/row)
-             S(i,j)=surnet+sealevel
+             S(i,j)=surnet+sealevel_2d(i,j) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
              B(i,j)=S(i,j)-H(i,j)
           end if
@@ -246,7 +246,4 @@
     end do
     !$OMP END DO
-    !~ print*,'flottab 2',S(132,183),H(132,183),BSOC(132,183),B(132,183),sealevel
-    !~ print*,'flottab 2',flot(132,183),ice(132,183)
-    !~ if (S(132,183).LT.sealevel) print*,'BUGGGGGGGGGGGGG !!!!!!!!!!!!!!!!'
        
 !!$ do i=1,nx
@@ -262,5 +259,4 @@
 
 
-
     !-----------------------------------------------------------------------
     !$OMP DO
@@ -298,5 +294,5 @@
           !                 sur le demi noeud condition archim < 100 m
 
-          archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-sealevel+ro/row*Hmx(i,j)
+          archim=(Bsoc(i,j)+Bsoc(i-1,j))*0.5-(sealevel_2d(i,j)+sealevel_2d(i-1,j))*0.5+ro/row*Hmx(i,j)
           gzmx(i,j)=gzmx(i,j).and.(archim.le.100.) 
           cotemx(i,j)=gzmx(i,j)
@@ -337,5 +333,5 @@
           !                 sur le demi noeud condition archim > 100 m
 
-          archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-sealevel+ro/row*Hmy(i,j)
+          archim=(Bsoc(i,j)+Bsoc(i,j-1))*0.5-(sealevel_2d(i,j)+sealevel_2d(i,j-1))*0.5+ro/row*Hmy(i,j)
           gzmy(i,j)=gzmy(i,j).and.(archim.le.100.) 
           cotemy(i,j)=gzmy(i,j)
@@ -605,5 +601,5 @@
           ice(:,:)=0
           H(:,:)=0.
-          S(:,:)=H(:,:)*(1.-ro/row) + sealevel
+          S(:,:)=H(:,:)*(1.-ro/row) + sealevel_2d(:,:)
           B(:,:)=S(:,:) - H(:,:)
        end where
@@ -790,5 +786,6 @@
             flgzmy(smax_i,smax_j)=.false. ; flgzmx(smax_i,smax_j+1)=.false.
      
-            if (Smax_.le.sealevel) then
+            if (Smax_.le.sealevel) then  ! afq -- WARNING: en fait avec le niveau marin local cest plus complique que ca!
+            !if (Smax_.le.sealevel_2d(i,j)) then ! afq --WARNING: il faudrait regarder point par point sur la tache...
               write(num_tracebug,*)'Attention, une ile avec la surface sous l eau'
               write(num_tracebug,*)'time=',time,'   coord:',smax_i,smax_j
Index: trunk/SOURCES/furst_schoof_mod.f90
===================================================================
--- trunk/SOURCES/furst_schoof_mod.f90	(revision 236)
+++ trunk/SOURCES/furst_schoof_mod.f90	(revision 237)
@@ -85,5 +85,5 @@
   gr_line_schoof(:,:)=0
 
-  archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel
+  archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot -sealevel_2d(:,:)
 
   imx_diag_in(:,:) = imx_diag(:,:)
@@ -93,5 +93,5 @@
      do i=3,nx-3
         !archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel
-        if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then    ! grounded and with ice
+        if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then    ! grounded and with ice
            !           print*,'schoof gr line ij',i,j
            ! selon x (indice i)
@@ -100,5 +100,5 @@
            !if (flot(i-1,j).and.Uxbar(i,j).LT.0.) then    ! grounding line between i-1,i    CAS WEST (ecoulement vers grille West)
 
-              call  calc_xgl4schoof(-dx,alpha_flot,H(i,j),Bsoc(i,j),H(i-1,j),Bsoc(i-1,j),Sealevel,xpos,Hgl)
+              call  calc_xgl4schoof(-dx,alpha_flot,H(i,j),Bsoc(i,j),sealevel_2d(i,j),H(i-1,j),Bsoc(i-1,j),sealevel_2d(i-1,j),xpos,Hgl)
               xpos_tab(i,j)=xpos
               Hglx_tab(i,j)=Hgl
@@ -194,5 +194,5 @@
            !else if (flot(i+1,j).and.Uxbar(i+1,j).GT.0.) then    ! grounding line between i,i+1    CAS EST (ecoulement vers grille Est)
 
-              call  calc_xgl4schoof(dx,alpha_flot,H(i,j),Bsoc(i,j),H(i+1,j),Bsoc(i+1,j),Sealevel,xpos,Hgl)
+              call  calc_xgl4schoof(dx,alpha_flot,H(i,j),Bsoc(i,j),sealevel_2d(i,j),H(i+1,j),Bsoc(i+1,j),sealevel_2d(i+1,j),xpos,Hgl)
               xpos_tab(i,j)=xpos
               Hglx_tab(i,j)=Hgl
@@ -285,5 +285,5 @@
 
 
-              call  calc_xgl4schoof(-dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j-1),Bsoc(i,j-1),Sealevel,ypos,Hgl)
+              call  calc_xgl4schoof(-dy,alpha_flot,H(i,j),Bsoc(i,j),sealevel_2d(i,j),H(i,j-1),Bsoc(i,j-1),sealevel_2d(i,j-1),ypos,Hgl)
               ypos_tab(i,j)=ypos
               Hgly_tab(i,j)=Hgl
@@ -373,5 +373,5 @@
            !else if (flot(i,j+1).and.Uybar(i,j+1).GT.0.) then    ! grounding line between j,j+1    CAS NORD (ecoulement vers grille Nord)
 
-              call  calc_xgl4schoof(dy,alpha_flot,H(i,j),Bsoc(i,j),H(i,j+1),Bsoc(i,j+1),Sealevel,ypos,Hgl)
+              call  calc_xgl4schoof(dy,alpha_flot,H(i,j),Bsoc(i,j),sealevel_2d(i,j),H(i,j+1),Bsoc(i,j+1),sealevel_2d(i,j+1),ypos,Hgl)
               ypos_tab(i,j)=ypos
               Hgly_tab(i,j)=Hgl
@@ -507,5 +507,5 @@
 !-------------------------------------------------------------------------
 ! calcule la position sous maille de la ligne dechouage
- subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,H_1,B_1,SL,xpos,Hgl)
+ subroutine calc_xgl4schoof(dy,alpha,H_0,B_0,SL_0,H_1,B_1,SL_1,xpos,Hgl)
                 
 
@@ -516,9 +516,10 @@
   real,intent(in)        ::    H_0       !< epaisseur au point pose
   real,intent(in)        ::    B_0       !< altitude socle au point pose
+  real,intent(in)        ::    SL_0      !< sea level au point pose
   real,intent(in)        ::    H_1       !< epaisseur au point flottant
   real,intent(in)        ::    B_1       !< altitude socle au point flottant
-  real,intent(in)        ::    SL        !< Sea level
+  real,intent(in)        ::    SL_1      !< sea level au point flottant
   real,intent(out)       ::    xpos      !< position de la ligne (en distance depuis le point pose)
-  real,intent(out)       ::    Hgl       !< Epaisseur de glace a la grounding line
+  real,intent(out)       ::    Hgl       !< epaisseur de glace a la grounding line
    
   
@@ -541,8 +542,8 @@
 !~ 	endif
 		
-  Cgl   = (B_1-B_0) + alpha * (H_1-H_0)
+  Cgl   = (B_1-B_0) + alpha * (H_1-H_0) - (SL_1-SL_0)
 
   if (abs(Cgl).gt.1.e-5) then   
-     xpos    = (SL - (B_0 + alpha * H_0)) / Cgl 
+     xpos    = (SL_0 - (B_0 + alpha * H_0)) / Cgl 
   else
      xpos    = 0.5  ! verifier
@@ -558,6 +559,6 @@
      print*,'calc_xgl4schoof : xpos value error i,j=',i,j
      print*, 'xpos,Cgl', xpos,Cgl
-     print*,'B_0, alpha, H_0, SL', B_0, alpha, H_0, SL
-     print*,'archim=',B_1+H_1*alpha - SL
+     print*,'B_0, alpha, H_0, SL_0', B_0, alpha, H_0, SL_0
+     print*,'archim=',B_1+H_1*alpha - SL_1
      !stop
      xpos = min(max(0.,xpos),1.)
Index: trunk/SOURCES/initial-phy-2-job.f90
===================================================================
--- trunk/SOURCES/initial-phy-2-job.f90	(revision 236)
+++ trunk/SOURCES/initial-phy-2-job.f90	(revision 237)
@@ -192,4 +192,5 @@
   TJFOR=0.0
   SEALEVEL=0.0
+  sealevel_2d(:,:)=0.
 
   SECYEAR=365.*24.*3600.
Index: trunk/SOURCES/initial-phy-2.f90
===================================================================
--- trunk/SOURCES/initial-phy-2.f90	(revision 236)
+++ trunk/SOURCES/initial-phy-2.f90	(revision 237)
@@ -185,4 +185,5 @@
   TJFOR=0.0
   SEALEVEL=0.0
+  sealevel_2d(:,:)=0.
 
   SECYEAR=365.*24.*3600.
Index: trunk/SOURCES/initial2-0.4.f90
===================================================================
--- trunk/SOURCES/initial2-0.4.f90	(revision 236)
+++ trunk/SOURCES/initial2-0.4.f90	(revision 237)
@@ -71,5 +71,5 @@
   do J=1,NY
      do I=1,NX
-        if ((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL).LT.0.) then
+        if ((BSOC(I,J)+H(I,J)*RO/ROW -SEALEVEL_2D(I,J)).LT.0.) then
            FLOT(I,J)=.TRUE.
         else
Index: trunk/SOURCES/isostasie_mod-0.3.f90
===================================================================
--- trunk/SOURCES/isostasie_mod-0.3.f90	(revision 236)
+++ trunk/SOURCES/isostasie_mod-0.3.f90	(revision 237)
@@ -103,5 +103,5 @@
     do i=1,nx
        do j=1,ny
-          if (ro*h0(i,j).ge.-row*(BSOC0(i,j)-sealevel)) then
+          if (ro*h0(i,j).ge.-row*(BSOC0(i,j)-sealevel_2d(i,j))) then
 
              !        glace ou terre
@@ -110,5 +110,5 @@
           else
              !        ocean
-             charge(i,j)=-rowg*(Bsoc0(i,j)-sealevel)
+             charge(i,j)=-rowg*(Bsoc0(i,j)-sealevel_2d(i,j))
           endif
        end do
Index: trunk/SOURCES/no_lakes.f90
===================================================================
--- trunk/SOURCES/no_lakes.f90	(revision 236)
+++ trunk/SOURCES/no_lakes.f90	(revision 237)
@@ -31,5 +31,5 @@
 subroutine lake_to_slv
 
-slv(:,:)=sealevel
+sealevel_2d(:,:)=sealevel
 
 end subroutine lake_to_slv
Index: trunk/SOURCES/prescribe-H-i2s_mod.f90
===================================================================
--- trunk/SOURCES/prescribe-H-i2s_mod.f90	(revision 236)
+++ trunk/SOURCES/prescribe-H-i2s_mod.f90	(revision 237)
@@ -198,5 +198,5 @@
 
     where (MK_gl0(:,:) .eq. 1)
-       hp(:,:) = (sealevel-Bsoc(:,:))*row/ro +1.
+       hp(:,:) = (sealevel_2d(:,:)-Bsoc(:,:))*row/ro +1.
     end where
 
Index: trunk/SOURCES/steps_time_loop.f90
===================================================================
--- trunk/SOURCES/steps_time_loop.f90	(revision 236)
+++ trunk/SOURCES/steps_time_loop.f90	(revision 237)
@@ -68,5 +68,5 @@
      where(.not.flot(:,:))
         S(:,:)=Bsoc(:,:)+H(:,:)
-        S(:,:)=max(S(:,:),sealevel)
+        S(:,:)=max(S(:,:),sealevel_2d(:,:)) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
         B(:,:)=Bsoc(:,:)
      end where
Index: trunk/SOURCES/steps_time_loop_avec_iterbeta.f90
===================================================================
--- trunk/SOURCES/steps_time_loop_avec_iterbeta.f90	(revision 236)
+++ trunk/SOURCES/steps_time_loop_avec_iterbeta.f90	(revision 237)
@@ -69,5 +69,5 @@
      where(.not.flot(:,:))
         S(:,:)=Bsoc(:,:)+H(:,:)
-        S(:,:)=max(S(:,:),sealevel)
+        S(:,:)=max(S(:,:),sealevel_2d(:,:)) !afq -- WARNING: est-ce qu'on veut vraiment mettre S a la valeur locale du niveau marin?
         B(:,:)=Bsoc(:,:)
      end where
Index: trunk/SOURCES/taubed-0.3.f90
===================================================================
--- trunk/SOURCES/taubed-0.3.f90	(revision 236)
+++ trunk/SOURCES/taubed-0.3.f90	(revision 237)
@@ -60,10 +60,10 @@
      do J=1,NY 
         do I=1,NX
-           if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL)) then
+           if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then
               !           glace ou terre 
               CHARGE(I,J)=ROG*H(I,J)
            else
               !           ocean
-              CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL)
+              CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL_2D(i,j))
            endif
         end do
@@ -102,10 +102,10 @@
      do J=1,NY
         do I=1,NX
-           if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL)) then
+           if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then
               !           glace ou terre 
               W1(I,J)=RO/ROM*H(I,J)
            else
               !           ocean
-              W1(I,J)=-ROW/ROM*(BSOC(I,J)-SEALEVEL)
+              W1(I,J)=-ROW/ROM*(BSOC(I,J)-SEALEVEL_2D(i,j))
            endif
         end do
