Index: trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc
===================================================================
--- trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc	(revision 64)
+++ trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc	(revision 65)
@@ -39,10 +39,10 @@
 	IFORT= ifort
 
-	ARITHMi    = -O2 -fp-model precise  -heap-arrays -traceback -mcmodel=medium -diag-disable warn #  -warn -traceback   -CB  -g  # options pour une meilleure arithmetique
+	ARITHMi    = -O2 -fp-model precise  -heap-arrays -traceback -mcmodel=medium #  -warn -traceback   -CB  -g  # options pour une meilleure arithmetique
                                                                       # (normalement reproductible)
     ifeq ($(debug), 1)
 		ARITHM    = $(ARITHMi) -CB -g -traceback -warn all
 	else
-		ARITHM    = $(ARITHMi)
+		ARITHM    = $(ARITHMi) -diag-disable warn
 	endif
 		
Index: trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90
===================================================================
--- trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90	(revision 64)
+++ trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_perturb_mod.f90	(revision 65)
@@ -379,5 +379,5 @@
   implicit none
 
-  integer i,j,ift
+  integer i,ift
   
   select case (massb_time)
Index: trunk/SOURCES/Hemin40_files/output_hemin40_mod.f90
===================================================================
--- trunk/SOURCES/Hemin40_files/output_hemin40_mod.f90	(revision 64)
+++ trunk/SOURCES/Hemin40_files/output_hemin40_mod.f90	(revision 65)
@@ -157,7 +157,5 @@
      integer inp(13) ! Surface posee (nb de noeuds)
      integer inf(13) ! surface flottante
-      integer npab, npcal
-
-      REAL long, lat, rd, cx, cy, cxmin, cymin, r2
+
 !Variables pour sommer      
       real isvol(13),isvolf(13) ! volume posé et flottants
@@ -165,5 +163,5 @@
       REAL ISABLBORD(13),ABLATOT(13),TACC(13),TBM(13)
       REAL ITJJA(13)
-      REAL ABLAMEAN
+!      REAL ABLAMEAN
 !moyennes utilisées en output
       REAL HMAX_(13) , HMEAN_(13) 
@@ -178,5 +176,5 @@
 
       REAL HDOTMEAN_G
-      REAL ABLA(NX,NY)
+!      REAL ABLA(NX,NY)
       REAL DELTAVOL
 
@@ -384,14 +382,5 @@
 905 format(f9.1,1x, e11.4,1x,e11.5, 1x,i5, 1x,e10.4,1x,i5 , 6(1x,e12.5), &
          12(2(1x,e10.4,1x,i5), 9(1x,e11.4) ) ) 
-900   format(i5,1x,f9.1,1x,f5.1,1x,e10.5,1x,i4,1x,i4,1x,f6.2,1X, &
-             f6.2,1X,f6.2,1X,f6.2,1X,i5,8(1x,f6.2))
-903   format(f9.1,1x,f4.1,1x,f5.1,1x,e10.5,1x,i5,1x,i4,1x,i5,1X, &
-             f7.3,1x,f6.2,1X,i6,4(1x,e7.1),1X,f6.2,1x,f7.3)
  
-904  format(F9.1,1x,e10.5,1x,e11.5,1x,i5,6(1x,e12.6),3(1x,e10.5,        &
-         1x,i5,6(1X,e12.6)),1X,i4,1X,i5,1x,f7.3,1x,                      &
-         f6.2,1x,i6,2(1x,e7.1),2(1x,f8.4),1X,f6.2,2(1x,i5),1x,e7.1,      &
-         8(1x,e12.6),4(1x,f15.7),9(1x,e10.5,1x,i5,2(1x,e12.6),1x,f15.7), &
-         26(1x,f8.3))
     
 940   format('%%%% ',a,'   time=',f8.0,' %%%%')
Index: trunk/SOURCES/Makefile
===================================================================
--- trunk/SOURCES/Makefile	(revision 64)
+++ trunk/SOURCES/Makefile	(revision 65)
@@ -1,2 +1,3 @@
+
 # Makefile pour construire les executables grisli en fonction de la machin
 # Pour compiler : make -f Makefile-template Grice2sea
@@ -14,5 +15,5 @@
 # Choice for librairies (mkl_c = 1 = MKL,  = 0 = BLAS)
 #Â NOTA: mkl_c = 1 with ifort = 0 ignores the MKL and uses BLAS anyhow
-mkl_c ?= 0
+mkl_c ?= 1
 
 #ifeq ($(ifort),1)
Index: trunk/SOURCES/Makefile.grisli.inc
===================================================================
--- trunk/SOURCES/Makefile.grisli.inc	(revision 64)
+++ trunk/SOURCES/Makefile.grisli.inc	(revision 65)
@@ -75,5 +75,5 @@
 	calving_frange.o no_calving.o no_lakes.o \
 	out_profile_mod.o printtable_mod.o mix-SIA-L1_mod.o \
-	relaxation_mod-0.3.o relaxation_water_mod-0.4.o relaxation_water_diffusion.o \
+	relaxation_water_diffusion.o \
 	prescribe-H-i2s_mod.o  \
 	resol_adv_diff_2D-sept2009.o  \
@@ -91,5 +91,5 @@
         calving_frange.o no_calving.o no_lakes.o \
         out_profile_mod.o printtable_mod.o mix-SIA-L1_mod.o \
-        relaxation_mod-0.3.o relaxation_water_mod-0.4.o relaxation_water_diffusion.o \
+        relaxation_water_diffusion.o \
         prescribe-H-i2s_mod.o   \
         resol_adv_diff_2D-sept2009.o  \
Index: trunk/SOURCES/Netcdf-routines/sortie_netcdf_GRISLI_mod.0.2-hassine.f90
===================================================================
--- trunk/SOURCES/Netcdf-routines/sortie_netcdf_GRISLI_mod.0.2-hassine.f90	(revision 64)
+++ trunk/SOURCES/Netcdf-routines/sortie_netcdf_GRISLI_mod.0.2-hassine.f90	(revision 65)
@@ -157,5 +157,5 @@
     integer :: num_file=22
     integer :: i1,i2,i3
-    integer :: i,j,k,l     
+    integer :: i,j,k     
     character(len=20) :: name1
     character(len=20) :: name2
@@ -383,5 +383,5 @@
     integer   :: ifin
     integer   :: npr
-    integer  :: i,j,k !< indices de travail
+    integer  :: i !< indices de travail
 
     if (itracebug.eq.1)  call tracebug(' Entree dans routine testsort_time_ncdf')
Index: trunk/SOURCES/New-remplimat/diagno-L2_mod.f90
===================================================================
--- trunk/SOURCES/New-remplimat/diagno-L2_mod.f90	(revision 64)
+++ trunk/SOURCES/New-remplimat/diagno-L2_mod.f90	(revision 65)
@@ -86,5 +86,4 @@
 subroutine diagnoshelf !      Resolution numerique des equations diagnostiques
 
-  integer :: itest
 
   if (itracebug.eq.1)  call tracebug(' Entree dans diagnoshelf')
Index: trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90
===================================================================
--- trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90	(revision 64)
+++ trunk/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90	(revision 65)
@@ -87,9 +87,9 @@
 !--------------------
 
-real    :: epsilon
 real    :: dx2=dx*dx     ! variable de travail
 real    :: beta          ! pour le frottement
-real    :: moteur        ! pour le frottement
 real    :: scal          ! pour le conditionnement (diagonale=1)
+
+character(len=80) :: fileout ! pour sortie fichier verification
 
 
@@ -233,11 +233,11 @@
 
 if (iter_beta.eq.1) then
-   call calc_beta(Vcol_x, Vcol_y) 
-
-   call write_datfile3(nx,ny,betamx,betamy,beta_centre,'beta-estime.dat')
-   if (itracebug.eq.1)  call tracebug(' remplimat apres calc_beta')
-   iter_beta = 0             ! pour faire un calcul de vitesse le coup d'apres
-   corr_def = .true.
-   goto 888
+	call calc_beta(Vcol_x, Vcol_y) 
+	fileout='beta-estime.dat'
+	call write_datfile3(nx,ny,betamx,betamy,beta_centre,fileout)
+	if (itracebug.eq.1)  call tracebug(' remplimat apres calc_beta')
+	iter_beta = 0             ! pour faire un calcul de vitesse le coup d'apres
+	corr_def = .true.
+	goto 888
 
 else if ((iter_beta.eq.0).and.(corr_def)) then
Index: trunk/SOURCES/Old-sources/relaxation_mod-0.3.f90
===================================================================
--- trunk/SOURCES/Old-sources/relaxation_mod-0.3.f90	(revision 65)
+++ trunk/SOURCES/Old-sources/relaxation_mod-0.3.f90	(revision 65)
@@ -0,0 +1,243 @@
+!> \file relaxation_mod-0.3.f90
+!! Module pour la resolution de l'equation de relaxation.
+!<
+
+!> \namespace RELAXATION_MOD
+!! Module pour la resolution de l'equation de relaxation.
+!! \author Christophe DUMAS
+!! \date fevrier 2000
+!! @note Cette routine est appellee par ICETHICK3
+!! @note Modele couple "ice sheet - ice shelves"
+!<
+
+MODULE RELAXATION_MOD
+
+CONTAINS
+!> SUBROUTINE: relaxation
+!! Routine pour la methode de relaxation
+!! \param NXX     [in]         defini la taille des tableaux
+!! \param NYY     [in]         defini la taille des tableaux 
+!! \param Mk      [in]         masks (ice sheet, max, above water, below water, 1)
+!! \param MARINE
+!! \param FLOT
+!! \param DT      [in]         pas de temps court 
+!! \param DX      [in]         pas en x
+!! \param vieuxH  [in]         H au pas de temps precedent
+!! \param BM      [in]         mass balance   'o'
+!! \param BMELT   [in]         basal melting  'o'
+!! \param UXBAR   [in]         vertically integrated velocity '>'
+!! \param UYBAR   [in]         vertically integrated velocity '^'
+!! \param H       [out]        ice thickness  'o'
+!! \param HDOT    [out]        ice thickness derivee / t  'o'
+!>
+  subroutine RELAXATION(NXX,NYY,MK,MARINE,FLOT,DT,DX,vieuxH,BM,BMELT,UXBAR,UYBAR,H,HDOT)
+
+    IMPLICIT NONE
+
+
+! declaration des variables en entree
+    INTEGER, intent(in) :: NXX, NYY !< defini la taille des tableaux
+    INTEGER,dimension(NXX,NYY), intent(in) :: MK     !< masks (ice sheet, max, above water, below water, 1)
+    LOGICAL, intent(in) :: MARINE
+    LOGICAL,dimension(NXX,NYY), intent(in) :: FLOT
+    REAL, intent(in) ::  DT                          !< pas de temps court
+    REAL, intent(in) ::  DX                          !< pas en x
+    REAL,dimension(NXX,NYY), intent(in) :: vieuxH    !< H au pas de temps precedent
+    REAL,dimension(NXX,NYY), intent(in) :: BM        !< mass balance   'o'
+    REAL,dimension(NXX,NYY), intent(in) :: BMELT     !< basal melting  'o'
+    REAL,dimension(NXX,NYY), intent(in) :: UXBAR     !< vertically integrated velocity '>'
+    REAL,dimension(NXX,NYY), intent(in) :: UYBAR     !< vertically integrated velocity '^'
+
+
+! declaration des variables en sortie
+    REAL,dimension(NXX,NYY), intent(out) :: H      ! ice thickness  'o'
+    REAL,dimension(NXX,NYY), intent(out) :: HDOT     ! ice thickness derivee / t  'o'
+
+
+! declaration des variables locales
+    INTEGER :: I,J
+    REAL  :: TESTH 
+    REAL,dimension(NXX,NYY) :: ARELAX,BRELAX,CRELAX,DRELAX,ERELAX,FRELAX
+    REAL,dimension(NXX,NYY) :: DELTAH
+    REAL :: RESTE,DELH,VH
+    INTEGER :: ntour
+    integer :: mbord
+    REAL  :: DTSDX 
+    LOGICAL :: STOPP
+
+
+!   print*,'dans relaxation mod'
+    H(:,:)= vieuxH(:,:)
+
+!   attribution des coefficients  arelax ....
+!   ----------------------------------------
+
+    dtsdx=dt/DX
+  
+    do J=2,NYY-1
+       do I=2,NXX-1
+
+
+             FRELAX(I,J)= VIEUXH(I,J)+(BM(I,J)-BMELT(I,J))*DT
+!          selon x
+!                                     ecoulement vers la droite
+             if ((UXBAR(I,J).ge.0.).and.(UXBAR(I+1,J).ge.0.)) then
+                ARELAX(I,J)=-UXBAR(I,J)*dtsdx
+                BRELAX(I,J)=0.
+                CRELAX(I,J)=UXBAR(I+1,J)*dtsdx
+              
+!                                     ecoulement vers la gauche
+             else if ((UXBAR(I,J).le.0.).and.(UXBAR(I+1,J).le.0.)) then
+                ARELAX(I,J)=0.
+                BRELAX(I,J)=UXBAR(I+1,J)*dtsdx
+                CRELAX(I,J)=-UXBAR(I,J)*dtsdx
+              
+!                                     ecoulement divergent
+             else if ((UXBAR(I,J).le.0.).and.(UXBAR(I+1,J).ge.0.)) then
+                ARELAX(I,J)=0.
+                BRELAX(I,J)=0.
+                CRELAX(I,J)=(UXBAR(I+1,J)-UXBAR(I,J))*dtsdx
+              
+!                                     ecoulement convergent
+             else if ((UXBAR(I,J).ge.0.).and.(UXBAR(I+1,J).le.0.)) then
+                ARELAX(I,J)=-UXBAR(I,J)*dtsdx
+                BRELAX(I,J)=UXBAR(I+1,J)*dtsdx
+                CRELAX(I,J)=0.
+             endif
+
+!          selon y 
+!                                 ecoulement vers le haut (sens axe y)
+             if ((UYBAR(I,J).ge.0.).and.(UYBAR(I,J+1).ge.0.)) then
+                DRELAX(I,J)=-UYBAR(I,J)*dtsdx
+                ERELAX(I,J)=0.
+                CRELAX(I,J)=1.+(CRELAX(I,J)+(UYBAR(I,J+1)*dtsdx))
+                
+!                                     ecoulement vers le bas
+             else if ((UYBAR(I,J).le.0.).and.(UYBAR(I,J+1).le.0.)) then
+                DRELAX(I,J)=0.
+                ERELAX(I,J)=UYBAR(I,J+1)*dtsdx
+                CRELAX(I,J)=1.+(CRELAX(I,J)-(UYBAR(I,J)*dtsdx))
+           
+!                                     ecoulement divergent
+             else if ((UYBAR(I,J).le.0.).and.(UYBAR(I,J+1).ge.0.)) then
+                DRELAX(I,J)=0.
+                ERELAX(I,J)=0.
+                CRELAX(I,J)=1.+(CRELAX(I,J)+(UYBAR(I,J+1)-UYBAR(I,J))*dtsdx)
+           
+!                                     ecoulement convergent
+             else if ((UYBAR(I,J).ge.0.).and.(UYBAR(I,J+1).le.0.)) then
+                DRELAX(I,J)=-UYBAR(I,J)*dtsdx
+                ERELAX(I,J)=UYBAR(I,J+1)*dtsdx
+                CRELAX(I,J)=1.+CRELAX(I,J)
+           
+             endif
+!          endif
+       end do
+    end do 
+
+! Boucle de relaxation :
+! ----------------------
+
+    STOPP = .false.
+    dtsdx=dt/dx
+    ntour=0
+
+    Do 362 WHILE(.NOT.STOPP)
+       ntour=ntour+1
+     
+       do j=2,NYY-1
+          do i=2,NXX-1
+            
+             RESTE = (((ARELAX(I,J)*H(I-1,J) + BRELAX(I,J)*H(I+1,J)) &
+                       + (DRELAX(I,J)*H(I,J-1) + ERELAX(I,J)*H(I,J+1))) &
+                       + CRELAX(I,J)*H(I,J))- FRELAX(I,J)
+     
+             DELTAH(I,J) = RESTE/CRELAX(I,J)
+!           H(I,J) = H(I,J) - DELTAH(I,J)
+          end do
+       end do
+
+       do j=2,NYY-1
+          do i=2,NXX-1
+             H(I,J) = H(I,J) - DELTAH(I,J)
+          end do
+       end do
+         
+! critere d'arret:
+! ----------------         
+
+       Delh=0
+       Vh=0
+
+       DO j=2,NYY-1
+          DO i=2,NXX-1
+!   write(6,*) I,J,delh,deltah(i,j)
+             Delh=Delh+deltah(i,j)**2
+!            Vh=Vh+h(i,j)**2.
+          END DO
+       END DO
+
+!      testh=SQRT(Delh/Vh)
+       if (delh.gt.0.) then
+          testh=SQRT(Delh)/((NXX-2)*(NYY-2))
+       else
+          testh=0.
+       endif
+       STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000)
+     
+     
+362  Continue
+
+! Conditions aux limites:
+! -----------------------
+
+!   ************ valeurs de H pour les coins ********
+
+       DO i=2,NXX-1
+          H(i,1) = 2.*H(i,2) - H(i,3)
+          H(i,NYY) = 2.*H(i,NYY-1) - H(i,NYY-2)
+       END DO
+       
+       DO J=2,NYY-1
+          H(1,J) = 2.*H(2,j) - H(3,j)
+          H(NXX,J) = 2.*H(NXX-1,j) - H(NXX-2,j)
+       END DO
+
+     
+       h(1,1)=(h(1,2)+h(2,1))/2.
+       Hdot(1,1)=0.
+       h(1,NYY)=(h(1,NYY-1)+h(2,NYY))/2.
+       hdot(1,NYY)=0.
+       h(NXX,1)=(h(NXX,2)+h(NXX-1,1))/2.
+       hdot(nxx,1)=0.
+       h(NXX,NYY)=(h(NXX,NYY-1)+h(NXX-1,NYY))/2.
+       hdot(nxx,nyy)=0.
+! =========================================================================
+
+
+!     open(87,file='test-H')
+!     write(87,*)'time=',time
+
+       do J=2,NYY-1
+          do I=2,NXX-1
+             HDOT(I,J)=(H(i,j)-vieuxH(i,j))/DT 
+
+
+
+!         if (abs(HDOT(I,J)).gt.hdotmax) then
+!  hdotmax=hdot(i,j)
+!  idotmax=i
+! jdotmax=j
+! end if
+
+!         write(87,*)i,j,h(i,j),hdot(i,j)
+          end do
+       end do
+
+
+!         close(87) 
+
+
+  end subroutine relaxation
+
+end module relaxation_mod
Index: trunk/SOURCES/Old-sources/relaxation_water_mod-0.4.f90
===================================================================
--- trunk/SOURCES/Old-sources/relaxation_water_mod-0.4.f90	(revision 65)
+++ trunk/SOURCES/Old-sources/relaxation_water_mod-0.4.f90	(revision 65)
@@ -0,0 +1,233 @@
+!> \file relaxation_water_mod-0.4.f90
+!! Module pour le calcule de l'ecoulement de l'eau sous la glace  
+!! selon une equation de Darcy avec un shema amont.     
+!<
+
+!> \namespace RELAXATION_WATER_MOD
+!! Module pour le calcule de l'ecoulement de l'eau sous la glace  
+!! selon une equation de Darcy avec un shema amont.    
+!! \author Vincent Peyaud 
+!! \date septembre 2004
+!! @note  reprise de relaxation_mod (pour le calcul de l'epaisseur)            
+!!     et adaptation pour le calcul de l'ecoulement de l'eau sous la glace  
+!!     selon une equation de Darcy avec un shema amont.                     
+!<
+
+
+! ************************************************************************ !
+!                                                                          !                                      
+!                       Vincent Peyaud septembre 2004                      !
+!     reprise de relaxation_mod (pour le calcul de l'epaisseur)            !
+!     et adaptation pour le calcul de l'ecoulement de l'eau sous la glace  !
+!     selon une equation de Darcy avec un shema amont.                     !
+!                                                                          !
+! ************************************************************************ !
+
+MODULE RELAXATION_WATER_MOD
+
+CONTAINS
+
+
+  subroutine RELAXATION_WAT(NXX,NYY,DT,DX,vieuxHWATER,BMELT,INFILTR,PGMX,PGMY,KOND,KONDMAX,HWATER)
+    implicit none
+
+
+    ! declaration des variables en entree
+    !------------------------------------------------
+    INTEGER, intent(in) :: NXX, NYY       !< defini la taille des tableaux
+    REAL,    intent(in) ::  DT            !< pas de temps court
+    REAL,    intent(in) ::  DX            !< pas en x
+    REAL,    intent(in) :: INFILTR        !< basal infiltration (lose of water) 
+    REAL,    intent(in) :: KONDMAX        !< maximum hydaulic conductivity (outside ice sheet)
+
+    REAL,dimension(NXX,NYY), intent(in) :: vieuxHWATER    !< H au pas de temps precedent 'o'
+    REAL,dimension(NXX,NYY), intent(in) :: BMELT     !< basal melting  'o'
+    REAL,dimension(NXX,NYY), intent(in) :: PGMX     !< hydaulic potential gratient '>' 
+    REAL,dimension(NXX,NYY), intent(in) :: PGMY     !< hydaulic potential gratient '^'
+    REAL,dimension(NXX,NYY), intent(in) :: KOND     !< hydaulic conductivity 'o'
+
+    ! declaration des variables en sortie
+    !-------------------------------------
+    REAL,dimension(NXX,NYY), intent(out):: HWATER      !< basal water thickness  'o' (pressure equivalent)
+
+
+    ! declaration des variables locales
+    !----------------------------------
+    INTEGER :: I,J
+    REAL  :: TESTH 
+    REAL,dimension(NXX,NYY) :: ARELAX,BRELAX,CRELAX,DRELAX,ERELAX,FRELAX
+    REAL,dimension(NXX,NYY) :: DELTAH
+    REAL :: RESTE,DELH,VH
+    INTEGER :: ntour
+    INTEGER :: mbord
+    REAL  :: DTSRGDX 
+    LOGICAL :: STOPP
+    REAL,dimension(NXX,NYY) :: KMX, KMY
+    !    REAL  :: RHO,RHOW,RHOG  !,SECYEAR!! ice density, water density, density*acceleration
+
+    HWATER(:,:)= vieuxHWATER(:,:)
+
+    !   calcul de kmx et kmx a partir de KOND
+    !   conductivite hyrdraulique sur les noeuds mineurs
+    !   moyenne harmonique
+    !   ----------------------------------------
+
+    do j=1,nyy
+       do i=2,nxx
+          if  (kond(i,j)==0.or.kond(i-1,j)==0) then
+             kmx(i,j)=0. ! to avoid division by o
+          else
+             kmx(i,j)=2*(kond(i,j)*kond(i-1,j))/(kond(i,j)+kond(i-1,j))
+          endif
+
+       end do
+    end do
+
+    do j=2,nyy
+       do i=2,nxx
+          if  (kond(i,j)==0.or.kond(i,j-1)==0) then
+             kmy(i,j)=0. ! to avoid division by o
+          else
+             kmy(i,j)=2*(kond(i,j)*kond(i,j-1))/(kond(i,j)+kond(i,j-1))
+          endif
+
+       enddo
+    enddo
+
+
+    !   attribution des coefficients  arelax ....
+    !   ----------------------------------------
+    !    SECYEAR=365.*24.*3600.
+    !   rho=910.
+    !    rhow=1000.
+    !    rhog=rhow*9.81
+    !    dtsrgdx=dt/(rhog*DX) a mon avis c'est rhow qu'il fallait utiliser. Maintenant cette 
+    !    division est faite dans eaubasale 
+
+    dtsrgdx=dt/DX
+
+
+    do J=2,NYY-1
+       do I=2,NXX-1
+
+          FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT !
+
+          !          selon x
+          !                                     ecoulement vers la droite (sens axe x)
+          if ((PGMX(I,J).ge.0.).and.(PGMX(I+1,J).ge.0.)) then
+             ARELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
+             BRELAX(I,J)=0.
+             CRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
+
+             !                                     ecoulement vers la gauche
+          else if ((PGMX(I,J).le.0.).and.(PGMX(I+1,J).le.0.)) then
+             ARELAX(I,J)=0.
+             BRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
+             CRELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
+
+             !                                     ecoulement divergent
+          else if ((PGMX(I,J).le.0.).and.(PGMX(I+1,J).ge.0.)) then
+             ARELAX(I,J)=0.
+             BRELAX(I,J)=0.
+             CRELAX(I,J)=(KMX(I+1,J)*PGMX(I+1,J)-KMX(I,J)*PGMX(I,J))*dtsrgdx
+
+             !                                     ecoulement convergent
+          else if ((PGMX(I,J).ge.0.).and.(PGMX(I+1,J).le.0.)) then
+             ARELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
+             BRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
+             CRELAX(I,J)=0.
+          endif
+  
+
+
+          !          selon y 
+          !                                 ecoulement vers le haut (sens axe y)
+          if ((PGMY(I,J).ge.0.).and.(PGMY(I,J+1).ge.0.)) then
+             DRELAX(I,J)=-KMY(I,J)*PGMY(I,J)*dtsrgdx
+             ERELAX(I,J)=0.
+             CRELAX(I,J)=CRELAX(I,J)+1. +(KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx)
+
+             !                                     ecoulement vers le bas
+          else if ((PGMY(I,J).le.0.).and.(PGMY(I,J+1).le.0.)) then
+             DRELAX(I,J)=0.
+             ERELAX(I,J)=KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx
+             CRELAX(I,J)=CRELAX(I,J)+1. -(KMY(I,J)*PGMY(I,J)*dtsrgdx)
+
+             !                                     ecoulement divergent
+          else if ((PGMY(I,J).le.0.).and.(PGMY(I,J+1).ge.0.)) then
+             DRELAX(I,J)=0.
+             ERELAX(I,J)=0.
+             CRELAX(I,J)=CRELAX(I,J)+1.+((KMY(I,J+1)*PGMY(I,J+1)-KMY(I,J)*PGMY(I,J))*dtsrgdx)
+
+             !                                     ecoulement convergent
+          else if ((PGMY(I,J).ge.0.).and.(PGMY(I,J+1).le.0.)) then
+             DRELAX(I,J)=-KMY(I,J)*PGMY(I,J)*dtsrgdx
+             ERELAX(I,J)=KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx
+             CRELAX(I,J)=CRELAX(I,J)+1.
+          endif
+
+       end do
+    end do
+
+    ! Boucle de relaxation :
+    ! ----------------------
+    STOPP = .false.
+
+    ntour=0
+
+    Do 362 WHILE(.NOT.STOPP)
+       ntour=ntour+1
+
+       do j=2,NYY-1
+          do i=2,NXX-1
+
+             RESTE = (((ARELAX(I,J)*HWATER(I-1,J) + BRELAX(I,J)*HWATER(I+1,J)) &
+                  + (DRELAX(I,J)*HWATER(I,J-1) + ERELAX(I,J)*HWATER(I,J+1))) &
+                  + CRELAX(I,J)*HWATER(I,J))- FRELAX(I,J)
+
+             DELTAH(I,J) = RESTE/CRELAX(I,J)             
+
+          end do
+       end do
+
+      
+
+       ! il faut faire le calcul suivant dans une autre boucle car RESTE est fonction
+       ! de hwater sur les points voisins.
+       do j=2,NYY-1
+          do i=2,NXX-1
+             HWATER(I,J) = HWATER(I,J) - DELTAH(I,J)
+          end do
+       end do
+
+       ! critere d'arret:
+
+      
+
+       Delh=0
+       Vh=0
+
+       DO j=2,NYY-1
+          DO i=2,NXX-1
+
+             !   write(6,*) I,J,delh,deltah(i,j)
+             Delh=Delh+deltah(i,j)**2
+             !            Vh=Vh+h(i,j)**2.
+          END DO
+       END DO
+
+!       write(6,*) delh,maxval(deltah),minval(deltah)
+       !      testh=SQRT(Delh/Vh)
+       if (delh.gt.0.) then
+          testh=SQRT(Delh)/((NXX-2)*(NYY-2))
+       else
+          testh=0.
+       endif
+       STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000)
+
+
+362    Continue
+
+     end subroutine relaxation_wat
+
+end module relaxation_water_mod
Index: trunk/SOURCES/Recul_force_grounding_line/toy_retreat_mod.f90
===================================================================
--- trunk/SOURCES/Recul_force_grounding_line/toy_retreat_mod.f90	(revision 64)
+++ trunk/SOURCES/Recul_force_grounding_line/toy_retreat_mod.f90	(revision 65)
@@ -115,5 +115,5 @@
   use declar_toy_retreat
   implicit none
-  integer                           :: som_voisins  !< pour des tests
+!  integer                           :: som_voisins  !< pour des tests
 
 
@@ -637,5 +637,5 @@
   implicit none
 
-  integer                           :: som_voisins  !< pour des tests
+!  integer                           :: som_voisins  !< pour des tests
   real                              :: Hf           !< Hf = 0 si socle > 0 sinon H_float
   real                              :: dyy         ! variable de travail  longueur de la maille en diagonale       
@@ -752,5 +752,5 @@
   implicit none
 
-  integer                           :: som_voisins  !< pour des tests
+!  integer                           :: som_voisins  !< pour des tests
 
   ! les nouveaux points a traiter
@@ -873,11 +873,11 @@
 integer    :: num_job           ! numero du job
 integer    :: lenrun            ! longueur du run name
-integer    :: lenb              ! longueur du nom de fichier 
-integer    :: pos               ! position d'un sous_string
+!integer    :: lenb              ! longueur du nom de fichier 
+!integer    :: pos               ! position d'un sous_string
 
 real*8, dimension(:,:),   pointer  :: tab         !< tableau 2d real ecrit dans le fichier
 character(len=4)                   :: char_num    ! pour la lecture
 character(len=80)                  :: file_B_test !< file name of the bedrock map to test topo instability
-character(len=80)                  :: file_voisin !< file name of the bedrock map min or max to test topo instability
+!character(len=80)                  :: file_voisin !< file name of the bedrock map min or max to test topo instability
 
 
@@ -1445,6 +1445,6 @@
 integer,intent(in)                 :: ii,jj         ! indice du point considere
 
-real                               :: dH_x          ! variable de travail
-real                               :: dH_y          ! variable de travail
+!real                               :: dH_x          ! variable de travail
+!real                               :: dH_y          ! variable de travail
 
 DelHdt_sanity(ii,jj) = 5000.
Index: trunk/SOURCES/climat-forcage-insolation_mod.f90
===================================================================
--- trunk/SOURCES/climat-forcage-insolation_mod.f90	(revision 64)
+++ trunk/SOURCES/climat-forcage-insolation_mod.f90	(revision 65)
@@ -149,14 +149,11 @@
 
   implicit none
-  integer :: mo,ti,tj  
-  character(len=8) :: control	! label to check clim. forc. file (filin) is usable
-  integer :: l 			! In snapshot files:the first column is the mask, read but not used
-  real :: T_surf_ref            ! variable de travail calcul temp a l'instant t a la surface S
+  
   integer :: intr
   integer :: igtr
   integer :: ictr
-  integer :: k,d1,d2
+  integer :: d1,d2
   integer :: i,j
-  character(len=100) :: file_ncdf      !< fichier netcdf 
+!  character(len=100) :: file_ncdf      !< fichier netcdf 
   real*8, dimension(:,:,:), pointer :: data_3D => null() ! donnees lues dans le netcdf
   real*8, dimension(:,:),pointer :: data_2D => null()    ! donnees lues dans le netcdf
@@ -450,5 +447,5 @@
   write(num_rep_42,'(A,2(f7.1,","))')   'palier_ice = ', palier_ice(:,:)
   write(num_rep_42,'(A,A)')   'surf_ice = ', surf_ice
-  write(num_rep_42,'(A,2(f3.1,","))')   'palier_CO2 = ', palier_CO2(:)
+  write(num_rep_42,'(A,2(f6.1,","))')   'palier_CO2 = ', palier_CO2(:)
   write(num_rep_42,'(A,A)')   'orb_file = ', orb_file
   write(num_rep_42,'(A,A)')   'co2_file = ', co2_file
@@ -494,5 +491,4 @@
 implicit none
 
-real COEFT,COEFP       !
 !integer l              ! dumm index for loops on snapshots files  l=ITR,NTR-1
 !cdc integer itr            ! index of the current snapshot file (change with time)
@@ -517,5 +513,5 @@
 !real,dimension(nx,ny) :: coeff1,coeff2
 real,dimension(nx,ny) :: coeff2
-integer :: k1,k2,m1,m2
+integer :: m1,m2
 integer,dimension(nx,ny) :: Mask,Mask2
 integer :: update_H
@@ -523,5 +519,5 @@
 
 !Definition des variables pour le calcul de la fonction de CO2
-double precision :: coefCO2
+!double precision :: coefCO2
 
 !*****************************************************************************
Index: trunk/SOURCES/climat-forcage-insolation_mod_oneway.f90
===================================================================
--- trunk/SOURCES/climat-forcage-insolation_mod_oneway.f90	(revision 64)
+++ trunk/SOURCES/climat-forcage-insolation_mod_oneway.f90	(revision 65)
@@ -85,13 +85,9 @@
 
   implicit none
-  integer :: mo,ti,tj  
-  character(len=8) :: control	! label to check clim. forc. file (filin) is usable
-  integer :: l 			! In snapshot files:the first column is the mask, read but not used
-  real :: T_surf_ref            ! variable de travail calcul temp a l'instant t a la surface S
+
   integer :: intr
   integer :: igtr
   integer :: ictr
-  integer :: i,j
-  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
+!  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
   real*8, dimension(:,:,:), pointer :: data_3D => null() ! donnees lues dans le netcdf
   real*8, dimension(:,:),pointer :: data_2D => null()    ! donnees lues dans le netcdf
@@ -231,11 +227,8 @@
 implicit none
 
-real COEFT,COEFP       !
 !integer l              ! dumm index for loops on snapshots files  l=ITR,NTR-1
 !cdc integer itr            ! index of the current snapshot file (change with time)
 
 integer mo
-integer :: ictr
-integer :: igtr
 integer :: i,j
 
Index: trunk/SOURCES/climat-perturb_mod-0.4.f90
===================================================================
--- trunk/SOURCES/climat-perturb_mod-0.4.f90	(revision 64)
+++ trunk/SOURCES/climat-perturb_mod-0.4.f90	(revision 65)
@@ -197,5 +197,5 @@
 !  use module3d_phy
   implicit none
-  integer :: i,j,ift
+  integer :: i,ift
 
   !       time en dehors des limites du fichier forcage
Index: trunk/SOURCES/climat_forcage_mois_mod.f90
===================================================================
--- trunk/SOURCES/climat_forcage_mois_mod.f90	(revision 64)
+++ trunk/SOURCES/climat_forcage_mois_mod.f90	(revision 65)
@@ -64,9 +64,8 @@
 
   implicit none
-  integer :: i,j,mo,ti,tj  
-  character(len=8) :: control	! label to check clim. forc. file (filin) is usable
-  integer :: l 			! In snapshot files:the first column is the mask, read but not used
-  real :: T_surf_ref            ! variable de travail calcul temp a l'instant t a la surface S
-  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
+!  integer :: i,j,mo,ti,tj  
+!  integer :: l 			! In snapshot files:the first column is the mask, read but not used
+!  real :: T_surf_ref            ! variable de travail calcul temp a l'instant t a la surface S
+!  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
   real*8, dimension(:,:,:), pointer :: data_3D => null() ! donnees lues dans le netcdf
 
Index: trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90
===================================================================
--- trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90	(revision 64)
+++ trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90	(revision 65)
@@ -83,6 +83,4 @@
 real,dimension(nx,ny) :: Dminx,Dminy
 real,dimension(nx,ny) :: Uxdiff,Uydiff           ! vitesse due a la diffusion
-integer :: it1,it2,it3
-
 real aux                                         ! pour le calcul du critere
 real maxdia                                      ! sur le pas de temps
Index: trunk/SOURCES/diffusiv-polyn-0.6.f90
===================================================================
--- trunk/SOURCES/diffusiv-polyn-0.6.f90	(revision 64)
+++ trunk/SOURCES/diffusiv-polyn-0.6.f90	(revision 65)
@@ -44,6 +44,4 @@
 REAL :: INV_4DX ! inverse de dx pour eviter les divisions =1/(4*dx)
 REAL :: INV_4DY ! inverse de dy pour eviter les divisions =1/(4*dy)
-INTEGER :: ii
-integer :: ll
 
 if (itracebug.eq.1)  call tracebug(' Entree dans routine diffusiv')
Index: trunk/SOURCES/initial-0.3.f90
===================================================================
--- trunk/SOURCES/initial-0.3.f90	(revision 64)
+++ trunk/SOURCES/initial-0.3.f90	(revision 65)
@@ -30,5 +30,4 @@
   use sorties_ncdf_grisli
   use util_recovery
-  character(len=80) :: filinit
 
   !-------------------------------------------------------------------------------------
Index: trunk/SOURCES/litho-0.4.f90
===================================================================
--- trunk/SOURCES/litho-0.4.f90	(revision 64)
+++ trunk/SOURCES/litho-0.4.f90	(revision 65)
@@ -45,5 +45,5 @@
       implicit none
 
-      INTEGER :: IP,JP,LPX,LPY,II,SOM1,SOM2
+      INTEGER :: II,SOM1,SOM2
       REAL, dimension(:,:), allocatable :: WLOC
 
Index: trunk/SOURCES/main3D-0.4-40km.f90
===================================================================
--- trunk/SOURCES/main3D-0.4-40km.f90	(revision 64)
+++ trunk/SOURCES/main3D-0.4-40km.f90	(revision 65)
@@ -116,7 +116,4 @@
 
 
-
-  integer :: m
-
   ! good luck
 
Index: trunk/SOURCES/neffect-0.4.f90
===================================================================
--- trunk/SOURCES/neffect-0.4.f90	(revision 64)
+++ trunk/SOURCES/neffect-0.4.f90	(revision 65)
@@ -24,10 +24,5 @@
       implicit none
 
-!!!!  real pwater,N1,N2,XFRAC,SLOPEM,uxmy,uymx
-
-      real XFRAC,SLOPEM,uxmy,uymx
-!     real,dimension(nx,ny) :: pwater !!! global var
-      real pentem,toslid,topente,neffoto,fonction_N
-      real,dimension(nx,ny) :: table_interm
+!      real,dimension(nx,ny) :: table_interm
       real :: Nefmin=1.e5                  ! Pression effective minimum (~ 10 m de glace)
 
Index: trunk/SOURCES/next-time.f90
===================================================================
--- trunk/SOURCES/next-time.f90	(revision 64)
+++ trunk/SOURCES/next-time.f90	(revision 65)
@@ -21,5 +21,4 @@
   real  (kind=kind(0.d0)) :: time
   real (kind=kind(0.d0)) :: time_loc
-  real (kind=kind(0.d0)) :: test_time
   real (kind=kind(0.d0)) :: time2_loc
   real (kind=kind(0.d0)) :: timesup
Index: trunk/SOURCES/no_spinup_mod.f90
===================================================================
--- trunk/SOURCES/no_spinup_mod.f90	(revision 64)
+++ trunk/SOURCES/no_spinup_mod.f90	(revision 65)
@@ -14,4 +14,5 @@
 module no_spinup
 
+use module3d_phy
 
 contains
@@ -22,5 +23,4 @@
 subroutine init_spinup
 
-use module3d_phy
 implicit none
 
Index: trunk/SOURCES/out_cptr_mod.f90
===================================================================
--- trunk/SOURCES/out_cptr_mod.f90	(revision 64)
+++ trunk/SOURCES/out_cptr_mod.f90	(revision 65)
@@ -387,6 +387,6 @@
   character(len=20),dimension(4)    :: dimnames3d      !< pour 3d troisieme dim est nz
   character(len=20),dimension(4)    :: dimnames3dT     !< pour 3d troisieme dim est nz+nzm
-  integer                           :: nrecs=1         !< compteur pour les enregistrements temps des variables  
-  integer                           :: idef=0          !< pour savoir si la variable a ete definie ou non
+!  integer                           :: nrecs=1         !< compteur pour les enregistrements temps des variables  
+!  integer                           :: idef=0          !< pour savoir si la variable a ete definie ou non
   integer                           :: isortie         !< pour le choix de type de fichier en sortie 
 
Index: trunk/SOURCES/out_horiz_mod.f90
===================================================================
--- trunk/SOURCES/out_horiz_mod.f90	(revision 64)
+++ trunk/SOURCES/out_horiz_mod.f90	(revision 65)
@@ -76,5 +76,5 @@
     integer :: i2
     integer :: i3
-    integer  :: i,j,k !< indices de travail
+    integer  :: k !< indices de travail
 
     integer :: num_dat = 21
@@ -194,5 +194,5 @@
 integer   :: ifin
 integer   :: npr
-integer  :: i,j,k ! indices de travail
+integer  :: i ! indices de travail
 
 
@@ -281,5 +281,4 @@
 real :: maxcol
 real :: coef
- integer  :: i,j,k ! indices de travail
 
 npos=npos+1
Index: trunk/SOURCES/out_profile_mod.f90
===================================================================
--- trunk/SOURCES/out_profile_mod.f90	(revision 64)
+++ trunk/SOURCES/out_profile_mod.f90	(revision 65)
@@ -146,5 +146,5 @@
 integer :: numtime
 !integer :: nmax
-integer :: n,l
+integer :: n
 integer, dimension(:), allocatable :: nbr_ligne_prof
 integer :: nbr_ligne_total
Index: trunk/SOURCES/relaxation_mod-0.3.f90
===================================================================
--- trunk/SOURCES/relaxation_mod-0.3.f90	(revision 64)
+++ 	(revision )
@@ -1,243 +1,0 @@
-!> \file relaxation_mod-0.3.f90
-!! Module pour la resolution de l'equation de relaxation.
-!<
-
-!> \namespace RELAXATION_MOD
-!! Module pour la resolution de l'equation de relaxation.
-!! \author Christophe DUMAS
-!! \date fevrier 2000
-!! @note Cette routine est appellee par ICETHICK3
-!! @note Modele couple "ice sheet - ice shelves"
-!<
-
-MODULE RELAXATION_MOD
-
-CONTAINS
-!> SUBROUTINE: relaxation
-!! Routine pour la methode de relaxation
-!! \param NXX     [in]         defini la taille des tableaux
-!! \param NYY     [in]         defini la taille des tableaux 
-!! \param Mk      [in]         masks (ice sheet, max, above water, below water, 1)
-!! \param MARINE
-!! \param FLOT
-!! \param DT      [in]         pas de temps court 
-!! \param DX      [in]         pas en x
-!! \param vieuxH  [in]         H au pas de temps precedent
-!! \param BM      [in]         mass balance   'o'
-!! \param BMELT   [in]         basal melting  'o'
-!! \param UXBAR   [in]         vertically integrated velocity '>'
-!! \param UYBAR   [in]         vertically integrated velocity '^'
-!! \param H       [out]        ice thickness  'o'
-!! \param HDOT    [out]        ice thickness derivee / t  'o'
-!>
-  subroutine RELAXATION(NXX,NYY,MK,MARINE,FLOT,DT,DX,vieuxH,BM,BMELT,UXBAR,UYBAR,H,HDOT)
-
-    IMPLICIT NONE
-
-
-! declaration des variables en entree
-    INTEGER, intent(in) :: NXX, NYY !< defini la taille des tableaux
-    INTEGER,dimension(NXX,NYY), intent(in) :: MK     !< masks (ice sheet, max, above water, below water, 1)
-    LOGICAL, intent(in) :: MARINE
-    LOGICAL,dimension(NXX,NYY), intent(in) :: FLOT
-    REAL, intent(in) ::  DT                          !< pas de temps court
-    REAL, intent(in) ::  DX                          !< pas en x
-    REAL,dimension(NXX,NYY), intent(in) :: vieuxH    !< H au pas de temps precedent
-    REAL,dimension(NXX,NYY), intent(in) :: BM        !< mass balance   'o'
-    REAL,dimension(NXX,NYY), intent(in) :: BMELT     !< basal melting  'o'
-    REAL,dimension(NXX,NYY), intent(in) :: UXBAR     !< vertically integrated velocity '>'
-    REAL,dimension(NXX,NYY), intent(in) :: UYBAR     !< vertically integrated velocity '^'
-
-
-! declaration des variables en sortie
-    REAL,dimension(NXX,NYY), intent(out) :: H      ! ice thickness  'o'
-    REAL,dimension(NXX,NYY), intent(out) :: HDOT     ! ice thickness derivee / t  'o'
-
-
-! declaration des variables locales
-    INTEGER :: I,J
-    REAL  :: TESTH 
-    REAL,dimension(NXX,NYY) :: ARELAX,BRELAX,CRELAX,DRELAX,ERELAX,FRELAX
-    REAL,dimension(NXX,NYY) :: DELTAH
-    REAL :: RESTE,DELH,VH
-    INTEGER :: ntour
-    integer :: mbord
-    REAL  :: DTSDX 
-    LOGICAL :: STOPP
-
-
-!   print*,'dans relaxation mod'
-    H(:,:)= vieuxH(:,:)
-
-!   attribution des coefficients  arelax ....
-!   ----------------------------------------
-
-    dtsdx=dt/DX
-  
-    do J=2,NYY-1
-       do I=2,NXX-1
-
-
-             FRELAX(I,J)= VIEUXH(I,J)+(BM(I,J)-BMELT(I,J))*DT
-!          selon x
-!                                     ecoulement vers la droite
-             if ((UXBAR(I,J).ge.0.).and.(UXBAR(I+1,J).ge.0.)) then
-                ARELAX(I,J)=-UXBAR(I,J)*dtsdx
-                BRELAX(I,J)=0.
-                CRELAX(I,J)=UXBAR(I+1,J)*dtsdx
-              
-!                                     ecoulement vers la gauche
-             else if ((UXBAR(I,J).le.0.).and.(UXBAR(I+1,J).le.0.)) then
-                ARELAX(I,J)=0.
-                BRELAX(I,J)=UXBAR(I+1,J)*dtsdx
-                CRELAX(I,J)=-UXBAR(I,J)*dtsdx
-              
-!                                     ecoulement divergent
-             else if ((UXBAR(I,J).le.0.).and.(UXBAR(I+1,J).ge.0.)) then
-                ARELAX(I,J)=0.
-                BRELAX(I,J)=0.
-                CRELAX(I,J)=(UXBAR(I+1,J)-UXBAR(I,J))*dtsdx
-              
-!                                     ecoulement convergent
-             else if ((UXBAR(I,J).ge.0.).and.(UXBAR(I+1,J).le.0.)) then
-                ARELAX(I,J)=-UXBAR(I,J)*dtsdx
-                BRELAX(I,J)=UXBAR(I+1,J)*dtsdx
-                CRELAX(I,J)=0.
-             endif
-
-!          selon y 
-!                                 ecoulement vers le haut (sens axe y)
-             if ((UYBAR(I,J).ge.0.).and.(UYBAR(I,J+1).ge.0.)) then
-                DRELAX(I,J)=-UYBAR(I,J)*dtsdx
-                ERELAX(I,J)=0.
-                CRELAX(I,J)=1.+(CRELAX(I,J)+(UYBAR(I,J+1)*dtsdx))
-                
-!                                     ecoulement vers le bas
-             else if ((UYBAR(I,J).le.0.).and.(UYBAR(I,J+1).le.0.)) then
-                DRELAX(I,J)=0.
-                ERELAX(I,J)=UYBAR(I,J+1)*dtsdx
-                CRELAX(I,J)=1.+(CRELAX(I,J)-(UYBAR(I,J)*dtsdx))
-           
-!                                     ecoulement divergent
-             else if ((UYBAR(I,J).le.0.).and.(UYBAR(I,J+1).ge.0.)) then
-                DRELAX(I,J)=0.
-                ERELAX(I,J)=0.
-                CRELAX(I,J)=1.+(CRELAX(I,J)+(UYBAR(I,J+1)-UYBAR(I,J))*dtsdx)
-           
-!                                     ecoulement convergent
-             else if ((UYBAR(I,J).ge.0.).and.(UYBAR(I,J+1).le.0.)) then
-                DRELAX(I,J)=-UYBAR(I,J)*dtsdx
-                ERELAX(I,J)=UYBAR(I,J+1)*dtsdx
-                CRELAX(I,J)=1.+CRELAX(I,J)
-           
-             endif
-!          endif
-       end do
-    end do 
-
-! Boucle de relaxation :
-! ----------------------
-
-    STOPP = .false.
-    dtsdx=dt/dx
-    ntour=0
-
-    Do 362 WHILE(.NOT.STOPP)
-       ntour=ntour+1
-     
-       do j=2,NYY-1
-          do i=2,NXX-1
-            
-             RESTE = (((ARELAX(I,J)*H(I-1,J) + BRELAX(I,J)*H(I+1,J)) &
-                       + (DRELAX(I,J)*H(I,J-1) + ERELAX(I,J)*H(I,J+1))) &
-                       + CRELAX(I,J)*H(I,J))- FRELAX(I,J)
-     
-             DELTAH(I,J) = RESTE/CRELAX(I,J)
-!           H(I,J) = H(I,J) - DELTAH(I,J)
-          end do
-       end do
-
-       do j=2,NYY-1
-          do i=2,NXX-1
-             H(I,J) = H(I,J) - DELTAH(I,J)
-          end do
-       end do
-         
-! critere d'arret:
-! ----------------         
-
-       Delh=0
-       Vh=0
-
-       DO j=2,NYY-1
-          DO i=2,NXX-1
-!   write(6,*) I,J,delh,deltah(i,j)
-             Delh=Delh+deltah(i,j)**2
-!            Vh=Vh+h(i,j)**2.
-          END DO
-       END DO
-
-!      testh=SQRT(Delh/Vh)
-       if (delh.gt.0.) then
-          testh=SQRT(Delh)/((NXX-2)*(NYY-2))
-       else
-          testh=0.
-       endif
-       STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000)
-     
-     
-362  Continue
-
-! Conditions aux limites:
-! -----------------------
-
-!   ************ valeurs de H pour les coins ********
-
-       DO i=2,NXX-1
-          H(i,1) = 2.*H(i,2) - H(i,3)
-          H(i,NYY) = 2.*H(i,NYY-1) - H(i,NYY-2)
-       END DO
-       
-       DO J=2,NYY-1
-          H(1,J) = 2.*H(2,j) - H(3,j)
-          H(NXX,J) = 2.*H(NXX-1,j) - H(NXX-2,j)
-       END DO
-
-     
-       h(1,1)=(h(1,2)+h(2,1))/2.
-       Hdot(1,1)=0.
-       h(1,NYY)=(h(1,NYY-1)+h(2,NYY))/2.
-       hdot(1,NYY)=0.
-       h(NXX,1)=(h(NXX,2)+h(NXX-1,1))/2.
-       hdot(nxx,1)=0.
-       h(NXX,NYY)=(h(NXX,NYY-1)+h(NXX-1,NYY))/2.
-       hdot(nxx,nyy)=0.
-! =========================================================================
-
-
-!     open(87,file='test-H')
-!     write(87,*)'time=',time
-
-       do J=2,NYY-1
-          do I=2,NXX-1
-             HDOT(I,J)=(H(i,j)-vieuxH(i,j))/DT 
-
-
-
-!         if (abs(HDOT(I,J)).gt.hdotmax) then
-!  hdotmax=hdot(i,j)
-!  idotmax=i
-! jdotmax=j
-! end if
-
-!         write(87,*)i,j,h(i,j),hdot(i,j)
-          end do
-       end do
-
-
-!         close(87) 
-
-
-  end subroutine relaxation
-
-end module relaxation_mod
Index: trunk/SOURCES/relaxation_water_diffusion.f90
===================================================================
--- trunk/SOURCES/relaxation_water_diffusion.f90	(revision 64)
+++ trunk/SOURCES/relaxation_water_diffusion.f90	(revision 65)
@@ -10,4 +10,6 @@
 !<
 module relaxation_waterdif_mod
+
+
 
 
@@ -48,5 +50,4 @@
     REAL :: RESTE,DELH,VH
     INTEGER :: ntour
-    INTEGER :: mbord
     REAL  :: DTSRGDX,dtwdx2
     LOGICAL :: STOPP
Index: trunk/SOURCES/relaxation_water_mod-0.4.f90
===================================================================
--- trunk/SOURCES/relaxation_water_mod-0.4.f90	(revision 64)
+++ 	(revision )
@@ -1,233 +1,0 @@
-!> \file relaxation_water_mod-0.4.f90
-!! Module pour le calcule de l'ecoulement de l'eau sous la glace  
-!! selon une equation de Darcy avec un shema amont.     
-!<
-
-!> \namespace RELAXATION_WATER_MOD
-!! Module pour le calcule de l'ecoulement de l'eau sous la glace  
-!! selon une equation de Darcy avec un shema amont.    
-!! \author Vincent Peyaud 
-!! \date septembre 2004
-!! @note  reprise de relaxation_mod (pour le calcul de l'epaisseur)            
-!!     et adaptation pour le calcul de l'ecoulement de l'eau sous la glace  
-!!     selon une equation de Darcy avec un shema amont.                     
-!<
-
-
-! ************************************************************************ !
-!                                                                          !                                      
-!                       Vincent Peyaud septembre 2004                      !
-!     reprise de relaxation_mod (pour le calcul de l'epaisseur)            !
-!     et adaptation pour le calcul de l'ecoulement de l'eau sous la glace  !
-!     selon une equation de Darcy avec un shema amont.                     !
-!                                                                          !
-! ************************************************************************ !
-
-MODULE RELAXATION_WATER_MOD
-
-CONTAINS
-
-
-  subroutine RELAXATION_WAT(NXX,NYY,DT,DX,vieuxHWATER,BMELT,INFILTR,PGMX,PGMY,KOND,KONDMAX,HWATER)
-    implicit none
-
-
-    ! declaration des variables en entree
-    !------------------------------------------------
-    INTEGER, intent(in) :: NXX, NYY       !< defini la taille des tableaux
-    REAL,    intent(in) ::  DT            !< pas de temps court
-    REAL,    intent(in) ::  DX            !< pas en x
-    REAL,    intent(in) :: INFILTR        !< basal infiltration (lose of water) 
-    REAL,    intent(in) :: KONDMAX        !< maximum hydaulic conductivity (outside ice sheet)
-
-    REAL,dimension(NXX,NYY), intent(in) :: vieuxHWATER    !< H au pas de temps precedent 'o'
-    REAL,dimension(NXX,NYY), intent(in) :: BMELT     !< basal melting  'o'
-    REAL,dimension(NXX,NYY), intent(in) :: PGMX     !< hydaulic potential gratient '>' 
-    REAL,dimension(NXX,NYY), intent(in) :: PGMY     !< hydaulic potential gratient '^'
-    REAL,dimension(NXX,NYY), intent(in) :: KOND     !< hydaulic conductivity 'o'
-
-    ! declaration des variables en sortie
-    !-------------------------------------
-    REAL,dimension(NXX,NYY), intent(out):: HWATER      !< basal water thickness  'o' (pressure equivalent)
-
-
-    ! declaration des variables locales
-    !----------------------------------
-    INTEGER :: I,J
-    REAL  :: TESTH 
-    REAL,dimension(NXX,NYY) :: ARELAX,BRELAX,CRELAX,DRELAX,ERELAX,FRELAX
-    REAL,dimension(NXX,NYY) :: DELTAH
-    REAL :: RESTE,DELH,VH
-    INTEGER :: ntour
-    INTEGER :: mbord
-    REAL  :: DTSRGDX 
-    LOGICAL :: STOPP
-    REAL,dimension(NXX,NYY) :: KMX, KMY
-    !    REAL  :: RHO,RHOW,RHOG  !,SECYEAR!! ice density, water density, density*acceleration
-
-    HWATER(:,:)= vieuxHWATER(:,:)
-
-    !   calcul de kmx et kmx a partir de KOND
-    !   conductivite hyrdraulique sur les noeuds mineurs
-    !   moyenne harmonique
-    !   ----------------------------------------
-
-    do j=1,nyy
-       do i=2,nxx
-          if  (kond(i,j)==0.or.kond(i-1,j)==0) then
-             kmx(i,j)=0. ! to avoid division by o
-          else
-             kmx(i,j)=2*(kond(i,j)*kond(i-1,j))/(kond(i,j)+kond(i-1,j))
-          endif
-
-       end do
-    end do
-
-    do j=2,nyy
-       do i=2,nxx
-          if  (kond(i,j)==0.or.kond(i,j-1)==0) then
-             kmy(i,j)=0. ! to avoid division by o
-          else
-             kmy(i,j)=2*(kond(i,j)*kond(i,j-1))/(kond(i,j)+kond(i,j-1))
-          endif
-
-       enddo
-    enddo
-
-
-    !   attribution des coefficients  arelax ....
-    !   ----------------------------------------
-    !    SECYEAR=365.*24.*3600.
-    !   rho=910.
-    !    rhow=1000.
-    !    rhog=rhow*9.81
-    !    dtsrgdx=dt/(rhog*DX) a mon avis c'est rhow qu'il fallait utiliser. Maintenant cette 
-    !    division est faite dans eaubasale 
-
-    dtsrgdx=dt/DX
-
-
-    do J=2,NYY-1
-       do I=2,NXX-1
-
-          FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT !
-
-          !          selon x
-          !                                     ecoulement vers la droite (sens axe x)
-          if ((PGMX(I,J).ge.0.).and.(PGMX(I+1,J).ge.0.)) then
-             ARELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
-             BRELAX(I,J)=0.
-             CRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
-
-             !                                     ecoulement vers la gauche
-          else if ((PGMX(I,J).le.0.).and.(PGMX(I+1,J).le.0.)) then
-             ARELAX(I,J)=0.
-             BRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
-             CRELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
-
-             !                                     ecoulement divergent
-          else if ((PGMX(I,J).le.0.).and.(PGMX(I+1,J).ge.0.)) then
-             ARELAX(I,J)=0.
-             BRELAX(I,J)=0.
-             CRELAX(I,J)=(KMX(I+1,J)*PGMX(I+1,J)-KMX(I,J)*PGMX(I,J))*dtsrgdx
-
-             !                                     ecoulement convergent
-          else if ((PGMX(I,J).ge.0.).and.(PGMX(I+1,J).le.0.)) then
-             ARELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
-             BRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
-             CRELAX(I,J)=0.
-          endif
-  
-
-
-          !          selon y 
-          !                                 ecoulement vers le haut (sens axe y)
-          if ((PGMY(I,J).ge.0.).and.(PGMY(I,J+1).ge.0.)) then
-             DRELAX(I,J)=-KMY(I,J)*PGMY(I,J)*dtsrgdx
-             ERELAX(I,J)=0.
-             CRELAX(I,J)=CRELAX(I,J)+1. +(KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx)
-
-             !                                     ecoulement vers le bas
-          else if ((PGMY(I,J).le.0.).and.(PGMY(I,J+1).le.0.)) then
-             DRELAX(I,J)=0.
-             ERELAX(I,J)=KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx
-             CRELAX(I,J)=CRELAX(I,J)+1. -(KMY(I,J)*PGMY(I,J)*dtsrgdx)
-
-             !                                     ecoulement divergent
-          else if ((PGMY(I,J).le.0.).and.(PGMY(I,J+1).ge.0.)) then
-             DRELAX(I,J)=0.
-             ERELAX(I,J)=0.
-             CRELAX(I,J)=CRELAX(I,J)+1.+((KMY(I,J+1)*PGMY(I,J+1)-KMY(I,J)*PGMY(I,J))*dtsrgdx)
-
-             !                                     ecoulement convergent
-          else if ((PGMY(I,J).ge.0.).and.(PGMY(I,J+1).le.0.)) then
-             DRELAX(I,J)=-KMY(I,J)*PGMY(I,J)*dtsrgdx
-             ERELAX(I,J)=KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx
-             CRELAX(I,J)=CRELAX(I,J)+1.
-          endif
-
-       end do
-    end do
-
-    ! Boucle de relaxation :
-    ! ----------------------
-    STOPP = .false.
-
-    ntour=0
-
-    Do 362 WHILE(.NOT.STOPP)
-       ntour=ntour+1
-
-       do j=2,NYY-1
-          do i=2,NXX-1
-
-             RESTE = (((ARELAX(I,J)*HWATER(I-1,J) + BRELAX(I,J)*HWATER(I+1,J)) &
-                  + (DRELAX(I,J)*HWATER(I,J-1) + ERELAX(I,J)*HWATER(I,J+1))) &
-                  + CRELAX(I,J)*HWATER(I,J))- FRELAX(I,J)
-
-             DELTAH(I,J) = RESTE/CRELAX(I,J)             
-
-          end do
-       end do
-
-      
-
-       ! il faut faire le calcul suivant dans une autre boucle car RESTE est fonction
-       ! de hwater sur les points voisins.
-       do j=2,NYY-1
-          do i=2,NXX-1
-             HWATER(I,J) = HWATER(I,J) - DELTAH(I,J)
-          end do
-       end do
-
-       ! critere d'arret:
-
-      
-
-       Delh=0
-       Vh=0
-
-       DO j=2,NYY-1
-          DO i=2,NXX-1
-
-             !   write(6,*) I,J,delh,deltah(i,j)
-             Delh=Delh+deltah(i,j)**2
-             !            Vh=Vh+h(i,j)**2.
-          END DO
-       END DO
-
-!       write(6,*) delh,maxval(deltah),minval(deltah)
-       !      testh=SQRT(Delh/Vh)
-       if (delh.gt.0.) then
-          testh=SQRT(Delh)/((NXX-2)*(NYY-2))
-       else
-          testh=0.
-       endif
-       STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000)
-
-
-362    Continue
-
-     end subroutine relaxation_wat
-
-end module relaxation_water_mod
Index: trunk/SOURCES/resol_adv_diff_2D-sept2009.f90
===================================================================
--- trunk/SOURCES/resol_adv_diff_2D-sept2009.f90	(revision 64)
+++ trunk/SOURCES/resol_adv_diff_2D-sept2009.f90	(revision 65)
@@ -146,7 +146,4 @@
 
 
-integer :: it1,it2,jt1,jt2           ! pour des tests d'asymétrie
-
-
 if (itracebug.eq.1)  call tracebug(' Entree dans routine resolution_diffusion')
 
Index: trunk/SOURCES/snaptime.f90
===================================================================
--- trunk/SOURCES/snaptime.f90	(revision 64)
+++ trunk/SOURCES/snaptime.f90	(revision 65)
@@ -24,5 +24,5 @@
 double precision :: tunit
 integer :: ipuissance
-integer :: itime,m
+integer :: itime
 character(len=1) :: signe,unit,mm 
 character(len=1),dimension (4) :: unittab
Index: trunk/SOURCES/spinup_mod.f90
===================================================================
--- trunk/SOURCES/spinup_mod.f90	(revision 64)
+++ trunk/SOURCES/spinup_mod.f90	(revision 65)
@@ -516,5 +516,4 @@
 
     implicit none
-    real,dimension(nx,ny) :: corr_ghf        !< correction du gflux geothermique
     real,dimension(nx,ny) :: coefdef_maj     !< coefficient deformation
     real :: increment_ghf
Index: trunk/SOURCES/write_datfile.f90
===================================================================
--- trunk/SOURCES/write_datfile.f90	(revision 64)
+++ trunk/SOURCES/write_datfile.f90	(revision 65)
@@ -22,6 +22,4 @@
 
 integer  :: i,j                ! working integers
-integer  :: lx,ly              ! nxx, nyy read in the file
-
 
 open(22,file=trim(filename))
@@ -55,6 +53,4 @@
 
 integer  :: i,j                ! working integers
-integer  :: lx,ly              ! nxx, nyy read in the file
-
 
 open(22,file=trim(filename))
