Index: /trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc
===================================================================
--- /trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc	(revision 75)
+++ /trunk/SOURCES/Fichiers-parametres/Makefile.tof-lsce3130.inc	(revision 76)
@@ -40,5 +40,5 @@
 	IFORT= ifort
 
-	ARITHMi = -O2 -fp-model precise
+	ARITHMi = -O2 -fp-model precise  #-openmp
                                                                       # (normalement reproductible)
     ifeq ($(debug), 1)
Index: /trunk/SOURCES/Temperature-routines/prop_th_icetemp.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/prop_th_icetemp.f90	(revision 75)
+++ /trunk/SOURCES/Temperature-routines/prop_th_icetemp.f90	(revision 76)
@@ -15,4 +15,5 @@
 
 Subroutine Thermal_prop_icetemp
+!$ USE OMP_LIB
   Use Icetemp_declar 
 
@@ -21,4 +22,6 @@
   If (Itracebug.Eq.1) Write(Num_tracebug,*)' Entree Dans Routine Thermal_prop_icetemp'
 
+!$OMP PARALLEL
+!$OMP DO COLLAPSE(2)
   Do K=1,Nz  
      Do J=1,Ny
@@ -29,6 +32,6 @@
            !  Attention Pour La Conductivite C'Est La Formule De Yin-Chao Yen
 
-           Ct(I,J,K)=6.62E7
-           Cp(I,J,K)=2009. 
+!           Ct(I,J,K)=6.62E7
+!           Cp(I,J,K)=2009. 
            Tpmp(I,J,K)=-0.00087*(K-1)*1./(Nz-1)*H(I,J)!De=1/(Nz_m-1)
            Cp(I,J,K)=(2115.3+7.79293*T(I,J,K)) 
@@ -38,5 +41,6 @@
      End Do
   End Do
-
+!$OMP END DO
+!$OMP END PARALLEL
 End Subroutine Thermal_prop_icetemp
 
Index: /trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90
===================================================================
--- /trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90	(revision 75)
+++ /trunk/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90	(revision 76)
@@ -80,4 +80,6 @@
 subroutine icethick3
 
+!$ USE OMP_LIB
+
 implicit none
 real,dimension(nx,ny) :: Dminx,Dminy
@@ -92,4 +94,6 @@
 ! Copie de H sur vieuxH        
 ! --------------------
+!$OMP PARALLEL
+!$OMP WORKSHARE
 where (mk0(:,:).eq.0)
    vieuxH(:,:)=0.
@@ -99,4 +103,5 @@
    vieuxH(:,:)=H(:,:)
 end where
+!$OMP END WORKSHARE
 
 ! mk0 est le masque à ne jamais dépasser
@@ -110,24 +115,34 @@
 jmaxdia=0
 
-! attention limitation ! 
+! attention limitation !
+!$OMP WORKSHARE 
 uxbar(:,:)=max(uxbar(:,:),-V_limit)
 uxbar(:,:)=min(uxbar(:,:),V_limit)
 uybar(:,:)=max(uybar(:,:),-V_limit)
 uybar(:,:)=min(uybar(:,:),V_limit)
+!$OMP END WORKSHARE
 
 ! le pas de temps est choisi pour rendre la matrice advective diagonale dominante
-do i=2,nx-1
-   do j=2,ny-1
+!$OMP END PARALLEL
+!$OMP PARALLEL
+!$OMP DO PRIVATE(aux)
+!do i=2,nx-1
+!   do j=2,ny-1
+do j=2,ny-1
+  do i=2,nx-1
       aux = (abs(uxbar(i,j)) + abs(uxbar(i+1,j))) &
            + (abs(uybar(i,j)) + abs(uybar(i,j+1)))
       aux=aux*dx1*0.5
+      !$OMP CRITICAL
       if (aux.gt.maxdia) then
          maxdia=aux
-         imaxdia=i
-         jmaxdia=j
+!         imaxdia=i
+!         jmaxdia=j
       endif
+      !$OMP END CRITICAL
    end do
 end do
-
+!$OMP END DO
+!$OMP END PARALLEL
 
 if (maxdia.ge.(testdiag/dtmax)) then
@@ -141,10 +156,11 @@
 ! avec le pas de temps long (dtt)
 
- call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug)
+call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug)
 
 
 ! calcul diffusivite - vitesse 
 !----------------------------------------------------------------- 
-
+!$OMP PARALLEL
+!$OMP WORKSHARE
 Uxdiff(:,:)=diffmx(:,:)*(-sdx(:,:))   ! vitesse qui peut s'exprimer en terme diffusif
 Uydiff(:,:)=diffmy(:,:)*(-sdy(:,:))   ! dans les where qui suivent elle est limitee a uxbar
@@ -172,5 +188,5 @@
   uydiff(:,:)=max(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot
 end where
-
+!$OMP END WORKSHARE
 
 
@@ -179,13 +195,15 @@
 
 if (adv_frac.gt.1.) then ! advection seulement
+   !$OMP WORKSHARE
    advmx(:,:)=Uxbar(:,:)
    advmy(:,:)=Uybar(:,:)
    Dmx(:,:)=0.
    Dmy(:,:)=0.
-
+   !$OMP END WORKSHARE
 else if (abs(adv_frac).lt.1.e-8) then   ! diffusion seulement
+   !$OMP WORKSHARE
    advmx(:,:)=0.
    advmy(:,:)=0.
-
+   !$OMP END WORKSHARE
    ! D reste la valeur au dessus)
 
@@ -197,5 +215,5 @@
 ! diffusion = ce qui peut être exprime en diffusion 
 !             + une partie (1-adv_frac) de l'advection
-
+   !$OMP WORKSHARE
    where ((abs(sdx(:,:)).gt.1.e-8).and.(Hmx(:,:).gt.2.))      ! cas general
        advmx(:,:)=(Uxbar(:,:)-Uxdiff(:,:))                   ! tout ce qui n'est pas diffusion
@@ -236,15 +254,12 @@
           dmy(:,:)=dminy(:,:)*H(:,:)
        end where
-
-
-
    elsewhere   ! zones trop plates ou trop fines
       Dmy(:,:)=0.
       advmy(:,:)=Uybar(:,:)
    end where
-
+   !$OMP END WORKSHARE
 
 else if (adv_frac.lt.0) then                  ! diffusif dans zones SIA, advectif ailleurs
-
+   !$OMP WORKSHARE
    zonemx(:,:)=flgzmx(:,:)
    zonemy(:,:)=flgzmy(:,:)
@@ -265,8 +280,8 @@
     advmy(:,:)=0.
    end where
-
+   !$OMP END WORKSHARE
 end if
 
-
+!$OMP WORKSHARE
 where (flot(:,:) )
    tabtest(:,:)=1.
@@ -274,5 +289,6 @@
    tabtest(:,:)=0.
 end where
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 !------------------------------------------------------------------detect
 !!$ tabtest(:,:)=dmx(:,:)
@@ -375,11 +391,14 @@
 
 
-debug_3D(:,:,45)=dmx(:,:)
-debug_3D(:,:,46)=dmy(:,:)
-debug_3D(:,:,47)=advmx(:,:)
-debug_3D(:,:,48)=advmy(:,:)
-
-
+!debug_3D(:,:,45)=dmx(:,:)
+!debug_3D(:,:,46)=dmy(:,:)
+!debug_3D(:,:,47)=advmx(:,:)
+!debug_3D(:,:,48)=advmy(:,:)
+
+!$OMP PARALLEL
+!$OMP WORKSHARE
 bilmass(:,:)=bm(:,:)-bmelt(:,:)                       ! surface and bottom mass balance
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 ! diverses precription de l'epaisseur en fonction de l'objectif
@@ -387,8 +406,5 @@
 
 
-  call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions)
-
-
-   
+call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions)
 
 
@@ -436,8 +452,8 @@
 !!$end if
 
-debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88)
-
-debug_3D(:,:,86)=Hp(:,:)
-debug_3D(:,:,85)=i_Hp(:,:)
+!debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88)
+
+!debug_3D(:,:,86)=Hp(:,:)
+!debug_3D(:,:,85)=i_Hp(:,:)
 
 !!$where (i_hp(:,:).eq.1)
@@ -463,5 +479,7 @@
 call resol_adv_diff_2D_vect (Dmx,Dmy,advmx,advmy,Hp,i_Hp,bilmass,vieuxH,H)
 
-! remise à 0 des epaisseurs negatives, on garde la difference dans ablbord 
+! remise à 0 des epaisseurs negatives, on garde la difference dans ablbord
+!$OMP PARALLEL
+!$OMP WORKSHARE 
 where (H(:,:).lt.0.)
    ablbord(:,:)=H(:,:)/dt
@@ -489,13 +507,16 @@
 endwhere
 
-
 ! eventuellement retirer apres spinup ou avoir un cas serac
 Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt
+!$OMP END WORKSHARE
 
 if (igrdline.ne.3) then
+   !$OMP WORKSHARE
    Hdot(:,:)=min(Hdot(:,:),10.)
    Hdot(:,:)=max(Hdot(:,:),-10.)
+   !$OMP END WORKSHARE
 endif
 
+!$OMP WORKSHARE
 where (i_hp(:,:).ne.1)
    H(:,:)=vieuxH(:,:)+Hdot(:,:)*dt
@@ -511,6 +532,6 @@
    Hdot(:,:)=0.
 end where
-
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 !calul de l'ablation sur les bords (pourrait n'être appelé que pour les sorties courtes)
Index: /trunk/SOURCES/diffusiv-polyn-0.6.f90
===================================================================
--- /trunk/SOURCES/diffusiv-polyn-0.6.f90	(revision 75)
+++ /trunk/SOURCES/diffusiv-polyn-0.6.f90	(revision 76)
@@ -34,5 +34,5 @@
 
 !     ===============================================================
-
+!$ USE OMP_LIB
 USE module3D_phy
 
@@ -58,7 +58,10 @@
 
 ! initialisation de la diffusion
+!$OMP PARALLEL
+!$OMP WORKSHARE
 diffmx(:,:)=0.
 diffmy(:,:)=0.
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 ! calcul de SDX, SDY et de la pente au carre en mx et en my :
@@ -79,5 +82,6 @@
 ! Calcul de la vitesse de glissement, qu'elle soit calculee par diagno ou par sliding
 !------------------------------------------------------------------------------------
-! 
+!$OMP PARALLEL
+!$OMP WORKSHARE
 where (flgzmx(:,:))
     ubx(:,:) = uxflgz(:,:)
@@ -91,6 +95,6 @@
     uby(:,:) = ddby(:,:)* (-sdy(:,:))
 endwhere
-
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 if (itracebug.eq.1) call tracebug(' diffusiv  apres calcul glissement')
 
@@ -99,32 +103,30 @@
 ! Deformation SIA : Calcul de ddy  et ddx pour chaque valeur de iglen
 !------------------------------------------------------------------
-
 do iglen=n1poly,n2poly
    glenexp=max(0.,(glen(iglen)-1.)/2.)
-
-
+  !$OMP PARALLEL 
+  !$OMP DO
    do j=1,ny
       do i=1,nx
-
          if (.not.flotmy(i,j)) then !  On calcule quand la deformation même pour les ice streams
             ddy(i,j,iglen)=((slope2my(i,j)**  &  ! slope2my calcule en debut de routine
                  glenexp)*(rog)**glen(iglen)) &
                  *hmy(i,j)**(glen(iglen)+1)
-
-
          endif
       end do
    end do
+   !$OMP END DO
+   !$OMP END PARALLEL
 end do
+
+
 
 do iglen=n1poly,n2poly
    glenexp=max(0.,(glen(iglen)-1.)/2.)
-
-
+   !$OMP PARALLEL
+   !$OMP DO
    do j=1,ny  
       do i=1,nx
-
          if (.not.flotmx(i,j)) then      ! bug y->x corrige le 5 avril 2008 (
-                                    
             ddx(i,j,iglen)=((slope2mx(i,j)**  & ! slope2mx calcule en debut de routine
                  glenexp)*(rog)**glen(iglen)) &
@@ -133,9 +135,12 @@
       end do
    end do
+   !$OMP END DO
+   !$OMP END PARALLEL
 end do
 
+
 ! vitesse due a la déformation-----------------------------------------------
-
-
+!$OMP PARALLEL
+!$OMP DO
 ydef: do j=2,ny
        do i=1,nx 
@@ -164,5 +169,7 @@
     end do
    end do ydef
-
+!$OMP END DO
+
+!$OMP DO
 xdef: do j=1,ny
        do i=2,nx 
@@ -192,4 +199,5 @@
     end do
    end do xdef
+!$OMP END DO
 
 if (itracebug.eq.1) call tracebug(' diffusiv  avant limit')
@@ -199,4 +207,5 @@
 
 ! la limitation selon x
+!$OMP WORKSHARE
 where(.not.flgzmx(:,:))
    uxbar(:,:)=min(uxbar(:,:),ultot)
@@ -209,6 +218,8 @@
    uybar(:,:)=max(uybar(:,:),-ultot)
 end where
+!$OMP END WORKSHARE
 
 !    cas particulier des sommets ou il ne reste plus de neige.
+!$OMP DO
 do j=2,ny-1
    do i=2,nx-1
@@ -221,5 +232,6 @@
    end do
 end do
-
+!$OMP END DO
+!$OMP END PARALLEL
 if (itracebug.eq.1)  call tracebug(' fin diffusiv ')
 
Index: /trunk/SOURCES/dragging_neff_slope_mod.f90
===================================================================
--- /trunk/SOURCES/dragging_neff_slope_mod.f90	(revision 75)
+++ /trunk/SOURCES/dragging_neff_slope_mod.f90	(revision 76)
@@ -168,10 +168,11 @@
 !-------------------------------------------------------------------------
 subroutine dragging   ! defini la localisation des streams et le frottement basal
-
+!$ USE OMP_LIB
 
 !         les iles n'ont pas de condition neff mais ont la condition toblim
 !         (ce bloc etait avant dans flottab)
 
-
+!$OMP PARALLEL
+!$OMP DO
 do j=2,ny
    do i=2,nx
@@ -180,4 +181,5 @@
    end do
 end do
+!$OMP END DO
 
 ! le calcul des fleuves se fait sur les mailles majeures en partant de la cote
@@ -193,11 +195,13 @@
 ! Attention : ce systeme ne permet pas d'avoir des fleuves convergents
 !             et les fleuves n'ont une épaisseur que de 1 noeud.
-
+!$OMP WORKSHARE
 fleuvemx(:,:)=.false.
 fleuvemy(:,:)=.false.
 fleuve(:,:)=.false.
 cote(:,:)=.false.
+!$OMP END WORKSHARE
 
 ! detection des cotes
+!$OMP DO
 do  j=2,ny-1  
    do i=2,nx-1  
@@ -208,6 +212,8 @@
    end do
 end do
+!$OMP END DO
 
 ! calcul de neff (pression effective sur noeuds majeurs)
+!$OMP DO
 do j=1,ny-1
    do i=1,nx-1
@@ -215,7 +221,10 @@
    enddo
 enddo
+!$OMP END DO
 !aurel, for the borders:
+!$OMP WORKSHARE
 neff(:,ny)=1.e5
 neff(nx,:)=1.e5
+
 
 !!$
@@ -277,5 +286,7 @@
 ! aurel, we add the neff threshold:
 where ((neff(:,:).le.seuil_neff).and.(.not.flot(:,:)).and.(H(:,:).gt.1.)) fleuve(:,:)=.true.
-
+!$OMP END WORKSHARE
+
+!$OMP DO
 do j=1,ny-1
    do i=1,nx-1
@@ -288,9 +299,13 @@
    end do
 end do
+!$OMP END DO
 
 ! we look for the warm based points that will not be treated as stream (ub from SSA):
+!$OMP WORKSHARE
 slowssa(:,:)=.false.
 slowssamx(:,:)=.false.
 slowssamy(:,:)=.false.
+!$OMP END WORKSHARE
+!$OMP DO
 do  j=1,ny
    do i=1,nx
@@ -299,4 +314,6 @@
    end do
 end do
+!$OMP END DO
+!$OMP DO
 do j=1,ny-1
    do i=1,nx-1
@@ -309,5 +326,6 @@
    end do
 end do
-
+!$OMP END DO
+!$OMP DO
 do j=1,ny
    do i=1,nx
@@ -354,5 +372,7 @@
    end do
 end do
-
+!$OMP END DO
+
+!$OMP DO
 do j=1,ny
    do i=1,nx
@@ -399,6 +419,7 @@
    end do
 end do
-
-
+!$OMP END DO
+
+!$OMP DO
 do j=2,ny-1
    do i=2,nx-1
@@ -407,45 +428,46 @@
    end do
 end do
-
-
-where (fleuve(:,:))
-   debug_3D(:,:,1)=1
-elsewhere (flot(:,:))
-   debug_3D(:,:,1)=2
-elsewhere
-   debug_3D(:,:,1)=0
-endwhere
-
-where (cote(:,:))
-   debug_3D(:,:,2)=1
-elsewhere
-   debug_3D(:,:,2)=0
-endwhere
-
-where (fleuvemx(:,:))
-   debug_3D(:,:,3)=1
-elsewhere
-   debug_3D(:,:,3)=0
-endwhere
-
-where (flotmx(:,:))
-   debug_3D(:,:,3)=-1
-endwhere
+!$OMP END DO
+!$OMP END PARALLEL
+
+!~ where (fleuve(:,:))
+!~    debug_3D(:,:,1)=1
+!~ elsewhere (flot(:,:))
+!~    debug_3D(:,:,1)=2
+!~ elsewhere
+!~    debug_3D(:,:,1)=0
+!~ endwhere
+!~ 
+!~ where (cote(:,:))
+!~    debug_3D(:,:,2)=1
+!~ elsewhere
+!~    debug_3D(:,:,2)=0
+!~ endwhere
+!~ 
+!~ where (fleuvemx(:,:))
+!~    debug_3D(:,:,3)=1
+!~ elsewhere
+!~    debug_3D(:,:,3)=0
+!~ endwhere
+!~ 
+!~ where (flotmx(:,:))
+!~    debug_3D(:,:,3)=-1
+!~ endwhere
 
 !_____________________
 
 
-where (fleuvemy(:,:))
-   debug_3D(:,:,4)=1
-elsewhere
-   debug_3D(:,:,4)=0
-endwhere
-
-where (flotmy(:,:))
-   debug_3D(:,:,4)=-1
-end where
-
-debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
-debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
+!~ where (fleuvemy(:,:))
+!~    debug_3D(:,:,4)=1
+!~ elsewhere
+!~    debug_3D(:,:,4)=0
+!~ endwhere
+!~ 
+!~ where (flotmy(:,:))
+!~    debug_3D(:,:,4)=-1
+!~ end where
+!~ 
+!~ debug_3D(:,:,23)= abs(RO*G*HMX(:,:)*sdx(:,:)*1.e-5) 
+!~ debug_3D(:,:,24)= abs(RO*G*HMY(:,:)*sdy(:,:)*1.e-5) 
 
 return
Index: /trunk/SOURCES/eaubasale-0.5_mod.f90
===================================================================
--- /trunk/SOURCES/eaubasale-0.5_mod.f90	(revision 75)
+++ /trunk/SOURCES/eaubasale-0.5_mod.f90	(revision 76)
@@ -133,6 +133,12 @@
 subroutine eaubasale !(pwater)   version correspondant à la thèse de Vincent
   ! A terme, il faudrait en faire un module 
-
-
+!$ USE OMP_LIB
+!~   integer :: t1,t2,ir
+!~   real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme  
+!~   
+!~    ! Temps CPU de calcul initial.
+!~    call cpu_time(t_cpu_0)
+!~    ! Temps elapsed de reference.
+!~    call system_clock(count=t1, count_rate=ir)
 
 
@@ -141,4 +147,6 @@
   if (itracebug.eq.1)  call tracebug(' Entree dans routine eau basale')
 
+!$OMP PARALLEL
+!$OMP WORKSHARE
   vieuxhwater(:,:) = hwater(:,:)
   kond(:,:) = kond0*SECYEAR 
@@ -148,4 +156,6 @@
   klimit(:,:)=0
   limit_hw(:,:)=-9999.
+!$OMP END WORKSHARE
+!$OMP DO  
   do j=1,ny
      do i=1,nx
@@ -167,23 +177,25 @@
      end do
   end do
-
-
-
-  do i=1,nx
-     do j=2,ny-1
+!$OMP END DO
+
+!$OMP DO
+!  do j=2,ny-1
+  do j=1,ny
+     do i=1,nx
         bmelt_w(i,j)=bmelt(i,j)*rofresh/ro    
         hwater(i,j)=max(hwater(i,j),0.)
         hw(i,j)=min(hwater(i,j),hmax_wat) 
-
      enddo
   enddo
-
+!$OMP END DO
+!$OMP END PARALLEL
 
   ecoul:  if (ecoulement_eau) then
      !  print*,'===dans eaubasale t, dt',time,dt,'kond 1.0e-5'
      !   write(6,*) 'entree ecoulement_eau'
-
-     do J=1,NY
-        do I=1,NX
+     !$OMP PARALLEL
+     !$OMP DO
+     do j=1,ny
+        do i=1,nx
            tetar(i,j)=(xlong(i,j))*PI/180. ! pourrait etre fait une fois pour toute
            PGX(I,J)=101*sin(tetar(i,j))*1.e-2                 
@@ -200,18 +212,16 @@
         enddo
      enddo
+!$OMP END DO
 
      ! sorties debug 17 juillet 2007
-     debug_3D(:,:,5)=pot_w(:,:)
-     debug_3D(:,:,6)=pot_f(:,:)
-     debug_3D(:,:,7)=pot_w(:,:)+rofreshg*hwater(:,:)
-     debug_3D(:,:,8)=hwater(:,:)
-
+!     debug_3D(:,:,5)=pot_w(:,:)
+!     debug_3D(:,:,6)=pot_f(:,:)
+!     debug_3D(:,:,7)=pot_w(:,:)+rofreshg*hwater(:,:)
+!     debug_3D(:,:,8)=hwater(:,:)
+!$OMP DO
      do j=2,ny
         do i=2,nx
            if (H(i,j).gt.25.) then
-
               !           calcul du gradient de pression
-              ! _______________________________________
-
               if (flotmx(i,j)) then 
                  pgx(i,j)=(pot_f(i-1,j)-pot_f(i,j))/dx
@@ -225,33 +235,26 @@
                  pgy(i,j)=(pot_w(i,j-1)-pot_w(i,j))/dy
               endif
-
-
            endif
-
-
            pgx(i,j)=pgx(i,j)/rofreshg   ! pour passer pgx sans unité
            pgy(i,j)=pgy(i,j)/rofreshg      
         end do
      end do
-
-
+!$OMP END DO
+!$OMP END PARALLEL
 
      if (nt.gt.0) then
         if (dt.gt.0.)  then !!!!!!!!!!!!!!!!!!!!!!relax_water si dt>0
-
-
            !-------------------------------------------------------------------------
            ! les points de la grounding line ont une conductivité hydraulique élevée
            ! même si la base est froide.  modif catritz 18 janvier 2005
            !open(166,file='erreur_eau')
-
-           do J=2,NY-1
-              do I=2,NX-1
-
+!$OMP PARALLEL
+!$OMP DO
+           do j=2,Ny-1
+              do i=2,Nx-1
                  ! cond=0 si glace froide et pas sur la grounding line
-
-                 if ((IBASE(I,J).eq.1).and.   &
+                 if ((IBASE(i,j).eq.1).and.   &
                       (.not.flot(i,j-1)).and.(.not.flot(i,j+1)).and.  &
-                      (.not.flot(i-1,j)).and.(.not.flot(i+1,j)))  KOND(I,J)=0.! 1.0e-5
+                      (.not.flot(i-1,j)).and.(.not.flot(i+1,j)))  KOND(i,j)=0.! 1.0e-5
 
                  ! cond infinie quand epaisseur faible et glace flottante
@@ -267,25 +270,20 @@
 
                  ! conductivite effective (conductivité * taille du tuyau en m2/an) 
-
                  keff(i,j)=kond(i,j)*hw(i,j)
               end do
            end do
-
+!$OMP END DO
            ! condition sur les bords de la grille
-
+!$OMP WORKSHARE
            kond(1,:)=kondmax
            kond(nx,:)=kondmax
            kond(:,1)=kondmax
            kond(:,ny)=kondmax
-
-           ! fin de la modif 18 janvier 2005---------------------------------------
-
            vieuxhwater(:,:) = hwater(:,:)
-
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
+!!$OMP SINGLE
            call relaxation_waterdif(nxlocal,nylocal,dt,dx,vieuxhwater,limit_hw,klimit,bmelt_w,infiltr,pgx,pgy,keff,keffmax,hwater)
-
-
-
+!!$OMP END SINGLE
         else
            !         print*,'dt=',dt,'pas de relax_water'
@@ -303,7 +301,8 @@
 
      ! on le fait avant les butoirs pour justement pouvoir les estimer
-
+!$OMP PARALLEL
      if (dt.gt.0.) then
         !         print*,'dt=',dt,'pas de relax_water'
+        !$OMP DO
         do j=1,ny
            do i=1,nx
@@ -311,7 +310,8 @@
            end do
         end do
+        !$OMP END DO
      endif
 
-
+!$OMP DO PRIVATE(Zflot)
      do i=1,nx
         do j=1,ny
@@ -353,10 +353,5 @@
               !          Hwater(i,j)=hw(i,j)
               !          pwater(i,j)=rog*H(I,J)*(hw(I,J)/hmax_wat)**(3.5)
-
-
-
-
            endif
-
 
            ! bloc qui pourrait servir pour mettre l'eau encore plus sous pression
@@ -365,11 +360,7 @@
            !             pwater(i,j)=pwater(i,j)+(HWATER(i,j)-poro_till*hmax_till)/(compress_w*hmax_till)
            !            endif
-
-
-
-
         end do
      end do
-
+!$OMP END DO
 
      !   ************ valeurs de HWATER pour les coins ********
@@ -381,4 +372,5 @@
 
      ! pour les sorties de flux d'eau
+     !$OMP DO
      do j=2,ny
         do i=2,nx
@@ -390,11 +382,10 @@
               phiwx(i,j)=phiwx(i,j)*2*(keff(i,j)*keff(i-1,j))/(keff(i,j)+keff(i-1,j))
            endif
-
            ! ligne pour sortir les pgx
-
            pgx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx)
         end do
      end do
-
+     !$OMP END DO
+     !$OMP DO
      do j=2,ny
         do i=2,nx
@@ -406,9 +397,8 @@
            endif
            pgy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx)
-
         enddo
      enddo
-
-
+     !$OMP END DO
+     !$OMP END PARALLEL
 
   else  ! partie originellement dans icetemp à mettre dans un autre module
@@ -416,4 +406,6 @@
 
      if (ISYNCHRO.eq.1) then
+        !$OMP PARALLEL
+        !$OMP DO
         do j=1,ny
            do i=1,nx
@@ -429,6 +421,23 @@
            end do
         end do
+        !$OMP END DO
+        !$OMP END PARALLEL
      end if
+     
   endif ecoul
+
+!~   ! Temps elapsed final
+!~   call system_clock(count=t2, count_rate=ir)
+!~   temps=real(t2 - t1,kind=4)/real(ir,kind=4)
+!~   ! Temps CPU de calcul final
+!~   call cpu_time(t_cpu_1)
+!~   t_cpu = t_cpu_1 - t_cpu_0
+
+!~   ! Impression du resultat.
+!~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              &
+!~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, &
+!~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, &
+!~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', &
+!~            nx,ny,temps,t_cpu,norme
 
   return
Index: /trunk/SOURCES/relaxation_water_diffusion.f90
===================================================================
--- /trunk/SOURCES/relaxation_water_diffusion.f90	(revision 75)
+++ /trunk/SOURCES/relaxation_water_diffusion.f90	(revision 76)
@@ -18,5 +18,8 @@
 
 subroutine relaxation_waterdif(NXX,NYY,DT,DX,vieuxHWATER,limit_hw,klimit,BMELT,INFILTR,PGMX,PGMY,KOND,KONDMAX,HWATER)
-    implicit none
+
+  !$ USE OMP_LIB
+
+  implicit none
 
 
@@ -57,11 +60,13 @@
 
 ! write(166,*)' entree  relaxation waterdif'
+!$OMP PARALLEL
+!$OMP WORKSHARE
     HWATER(:,:)= vieuxHWATER(:,:)
-
+!$OMP END WORKSHARE
     !   calcul de kmx et kmx a partir de KOND
     !   conductivite hyrdraulique sur les noeuds mineurs
     !   moyenne harmonique
     !   ----------------------------------------
-
+!$OMP DO
     do j=2,nyy
        do i=2,nxx
@@ -75,5 +80,7 @@
        end do
     end do
-
+!$OMP END DO
+
+!$OMP DO
     do j=2,nyy
        do i=2,nxx
@@ -86,5 +93,5 @@
        enddo
     enddo
-
+!$OMP END DO
 
     !   attribution des coefficients  arelax ....
@@ -100,24 +107,22 @@
     dtwdx2=dt/dx/dx
 
-arelax(:,:)=0.
-brelax(:,:)=0.
-crelax(:,:)=1.
-drelax(:,:)=0.
-erelax(:,:)=0.
-frelax(:,:)=limit_hw(:,:)
-
-
-
+!$OMP WORKSHARE
+  arelax(:,:)=0.
+  brelax(:,:)=0.
+  crelax(:,:)=1.
+  drelax(:,:)=0.
+  erelax(:,:)=0.
+  frelax(:,:)=limit_hw(:,:)
+!$OMP END WORKSHARE
+
+!$OMP DO
     do J=2,NYY-1
        do I=2,NXX-1
 
           if (klimit(i,j).eq.0) then
-
 ! calcul du vecteur
-
           FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT
           frelax(i,j)=frelax(i,j)+(kmx(i,j)*pgmx(i,j)-kmx(i+1,j)*pgmx(i+1,j))*dtsrgdx
           frelax(i,j)=frelax(i,j)+(kmy(i,j)*pgmy(i,j)-kmy(i,j+1)*pgmy(i,j+1))*dtsrgdx
-
 ! calcul des diagonales      
           arelax(i,j)=-kmx(i,j)*dtwdx2        ! arelax : diagonale i-1,j 
@@ -131,15 +136,12 @@
           crelax(i,j)=1.+((kmx(i,j)+kmx(i+1,j))+(kmy(i,j)+kmy(i,j+1)))*dtwdx2 
                                               !crelax : diagonale i,j 
-
           else if (klimit(i,j).eq.1) then
              hwater(i,j)=limit_hw(i,j)
 !             write(6,*) i,j,hwater(i,j),crelax(i,j),frelax(i,j),arelax(i,j)
-
           endif
-
        end do
     end do
-
-
+!$OMP END DO
+!$OMP END PARALLEL
 
     ! Boucle de relaxation :
@@ -152,5 +154,6 @@
        ntour=ntour+1
 !       write(6,*) 'boucle de relaxation numero',ntour
-
+       !$OMP PARALLEL
+       !$OMP DO PRIVATE(reste)
        do j=2,NYY-1
           do i=2,NXX-1
@@ -161,39 +164,37 @@
 
              DELTAH(I,J) = RESTE/CRELAX(I,J)             
-
           end do
        end do
-
-deltah(:,:)=min(deltah(:,:),10.)
-deltah(:,:)=max(deltah(:,:),-10.)
-     
-
-
+       !$OMP END DO
+
+       !$OMP WORKSHARE
+       deltah(:,:)=min(deltah(:,:),10.)
+       deltah(:,:)=max(deltah(:,:),-10.)
+       !$OMP END WORKSHARE
        ! il faut faire le calcul suivant dans une autre boucle car RESTE est fonction
        ! de hwater sur les points voisins.
+       !$OMP DO
        do j=2,NYY-1
           do i=2,NXX-1
              HWATER(I,J) = HWATER(I,J) - DELTAH(I,J)
-
           end do
        end do
-
+       !$OMP END DO
 
        ! critere d'arret:
-
-      
-
        Delh=0
        Vh=0
-
+       
+       !$OMP DO REDUCTION(+:Delh)
        DO j=2,NYY-1
           DO i=2,NXX-1
              
              !   write(166,*) I,J,delh,deltah(i,j)
-                         Delh=Delh+deltah(i,j)**2
+             Delh=Delh+deltah(i,j)**2
              !            Vh=Vh+h(i,j)**2.
           END DO
        END DO
-
+       !$OMP END DO
+       !$OMP END PARALLEL
 !       write(6,*) delh,maxval(deltah),minval(deltah)
        !      testh=SQRT(Delh/Vh)
