Index: /trunk/SOURCES/Temperature-routines/advec_icetemp.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/advec_icetemp.f90	(revision 72)
+++ /trunk/SOURCES/Temperature-routines/advec_icetemp.f90	(revision 73)
@@ -14,17 +14,43 @@
 !>
 
-Subroutine  Advec_icetemp(Ii,Jj,Prodq_m,Cro_m,Advecx_m,Advecy_m,Advec_m,Ct_m)
+Subroutine  Advec_icetemp(Nz,Nzm,Nn,Prodq_m,Cro_m,Advecx_m,Advecy_m,Advec_m,Ct_m,Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n, &
+	Uxij,Uxi_plus_un,Uyij,Uyj_plus_un,Tij,Ti_moins_un,Ti_plus_un,Tj_moins_un,Tj_plus_un,Hij,Tsij,Aa,Bb,Cc,Rr, &
+	Dx11,Dou,DTT,Dee,Uzrij,Dttdx)
 
-  use icetemp_declar 
+!  use icetemp_declar 
 
   Implicit None
     !< Arguments
-  Integer,                     Intent(In)  ::  Ii,Jj          !< Indice Selon X Et Y
-  Real, Dimension(Nz),         Intent(In)   :: Prodq_m        !< Tableau 1d Vert. heat production
-  Real, Dimension(Nz),         Intent(In)   :: Cro_m          !< Tableau 1d Vert. heat capcity
-  Real, Dimension(Nz),         Intent(Inout):: Advecx_m       !< Advection selon x
-  Real, Dimension(Nz),         Intent(Inout):: Advecy_m       !< Advection selon y
-  Real, Dimension(Nz),         Intent(Inout):: Advec_m        !< Advection total
-  Real, Dimension(Nz),         Intent(In)   :: Ct_m           !< Tableau 1d Vert.  thermal cond.
+  Integer,                     intent(in)  ::  Nz,Nzm,Nn      !< Taille des tableaux
+  Real, Dimension(Nz),         intent(in)   :: Prodq_m        !< Tableau 1d Vert. heat production
+  Real, Dimension(Nz),         intent(in)   :: Cro_m          !< Tableau 1d Vert. heat capcity
+  Real, Dimension(Nz),         intent(inout):: Advecx_m       !< Advection selon x
+  Real, Dimension(Nz),         intent(inout):: Advecy_m       !< Advection selon y
+  Real, Dimension(Nz),         intent(inout):: Advec_m        !< Advection total
+  Real, Dimension(Nz),         intent(in)   :: Ct_m           !< Tableau 1d Vert.  thermal cond.
+  
+  Integer, intent(in) :: Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n
+  real, dimension(Nz), Intent(in) :: Uxij, Uxi_plus_un
+  real, dimension(Nz), Intent(in) :: Uyij, Uyj_plus_un
+  real, dimension(Nz+Nzm), Intent(inout) :: Tij
+  real, dimension(Nz+Nzm), Intent(in) :: Ti_moins_un,Ti_plus_un,Tj_moins_un,Tj_plus_un
+  real, intent(in) :: Hij
+  real,intent(in) :: Tsij
+  Real,Dimension(Nn),intent(inout) :: Aa    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn),intent(inout) :: Bb    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn),intent(inout) :: Cc    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn),intent(inout) :: Rr    !< Work Arrays For Tridiag  !Dim Nn
+  Real, intent(in) :: Dx11
+  Real, intent(inout) :: Dou
+  real, intent(in) :: DTT
+  real, intent(in) :: Dee
+  real, dimension(Nz) :: Uzrij
+  real :: Dttdx
+  
+
+  
+  ! variables locales
+  Integer :: K
+  Real :: Dzz,Dah,Ct_haut,Ct_bas
 
  
@@ -35,9 +61,9 @@
      !         Avection Selon X  (advecx_m en general > 0)
      !         ----------------
-     Advecx_m(K) = Iadvec_w(Ii,Jj) * Ux(Ii,Jj,K) *          & ! ux west if upwind
-          (T(Ii,Jj,K) - T(Ii-1,Jj,K))  & ! west T gradient
+     Advecx_m(K) = Iadvec_w * Uxij(K) *          & ! ux west if upwind
+          (Tij(K) - Ti_moins_un(K))  & ! west T gradient
             
-          + Iadvec_e(Ii+1,Jj)* Ux(Ii+1,Jj,K)*      & ! ux east if upwind
-          (T(Ii+1,Jj,K)-T(Ii,Jj,K))     ! east T gradient
+          + Iadvec_e * Uxi_plus_un(K)*      & ! ux east if upwind
+          (Ti_plus_un(K)-Tij(K))     ! east T gradient
 
      Advecx_m(K) =  Advecx_m(K) * Dx11                                        !Dx11=1/Dx
@@ -46,9 +72,9 @@
      !         Avection Selon Y
      !         ----------------
-     Advecy_m(K) =   Iadvec_s(Ii,Jj) * Uy(Ii,Jj,K) *        & ! uy sud if upwind
-          (T(Ii,Jj,K) - T(Ii,Jj-1,K))& ! south T gradient
+     Advecy_m(K) =   Iadvec_s * Uyij(K) *        & ! uy sud if upwind
+          (Tij(K) - Tj_moins_un(K))& ! south T gradient
           
-          + Iadvec_n(Ii,Jj+1) * Uy(Ii,Jj+1,K) *   & ! uy nord is upwind
-          (T(Ii,Jj+1,K)-T(Ii,Jj,K))    ! north T gradient
+          + Iadvec_n * Uyj_plus_un(K) *   & ! uy nord is upwind
+          (Tj_plus_un(K)-Tij(K))    ! north T gradient
 
      Advecy_m(K) =  Advecy_m(K) * Dx11
@@ -62,12 +88,12 @@
 
   !      -----------------------------Cas General (H>10m)
-  thick_ice:  If (H(Ii,Jj).Gt.10.) Then
+  thick_ice:  If (Hij.Gt.10.) Then
 
      !        Variables De Calcul Dans La Glace
      ! dou = dtt/dz^2
-     Dou=Dtt/Dee/Dee/ H(Ii,Jj) / H(Ii,Jj)  
+     Dou=Dtt/Dee/Dee/ Hij / Hij  
 
      ! Dah = dtt/dz
-     Dah=Dtt /Dee /H(Ii,Jj)
+     Dah=Dtt /Dee /Hij
 
      ! thermal conductivity at mid point just below the surface                   
@@ -75,5 +101,5 @@
 
      ! surface temperature : cannot go above 0 celsius
-     T(Ii,Jj,1)=Min(0.,Ts(Ii,Jj))
+     Tij(1)=Min(0.,Tsij)
 
      ! surface boundary condition
@@ -81,5 +107,5 @@
      Bb(1)=1.
      Cc(1)=0.
-     Rr(1)=T(Ii,Jj,1)
+     Rr(1)=Tij(1)
 
      Do K=2,Nz-1
@@ -91,9 +117,9 @@
 
         ! Advection Verticale Centree
-        Aa(K) = -Dzz * Ct_haut - Uzr(Ii,Jj,K) * Dah / 2.  ! lower diag
-        Cc(K) = -Dzz * Ct_bas  + Uzr(Ii,Jj,K) * Dah / 2.  ! upper diag
+        Aa(K) = -Dzz * Ct_haut - Uzrij(K) * Dah / 2.  ! lower diag
+        Cc(K) = -Dzz * Ct_bas  + Uzrij(K) * Dah / 2.  ! upper diag
         Bb(K) =  1.+ Dzz * (Ct_haut+Ct_bas)                                ! diag
         
-        Rr(K) =  T(Ii,Jj,K) + Prodq_m(K) * Dtt   ! vector
+        Rr(K) =  Tij(K) + Prodq_m(K) * Dtt   ! vector
         Rr(K) =  Rr(K)- Dttdx * (Advec_m(K))
 
Index: /trunk/SOURCES/Temperature-routines/icetemp_declar_mod.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/icetemp_declar_mod.f90	(revision 72)
+++ /trunk/SOURCES/Temperature-routines/icetemp_declar_mod.f90	(revision 73)
@@ -22,6 +22,6 @@
   Integer :: Nfracq   !< Exposant Fracq
   Integer :: Iq       !< Choix Du Type De Routine Chaleur
-  Real :: Sx,Sy,Sx2,Sy2,Deh22,Tss
-  Real ::  Dou,Dah,Duu,Dzz,Dzi,Chalbed,Ct_bas,Ct_haut
+  Real :: Sx,Sy,Sx2,Sy2,Deh22 ! ,Tss
+  Real ::  Dou,Duu,Chalbed ! ,Dah,Dzz,Dzi,Ct_bas,Ct_haut
 
   Real,Parameter :: Acof1=-0.0575  !< Pour La Temperature De L'Eau De Mer
@@ -40,30 +40,30 @@
   ! _______________
 
-  Real,allocatable,Dimension(:) :: Aa    !< Work Arrays For Tridiag  !Dim Nn
-  Real,allocatable,Dimension(:) :: Bb    !< Work Arrays For Tridiag  !Dim Nn
-  Real,allocatable,Dimension(:) :: Cc    !< Work Arrays For Tridiag  !Dim Nn
-  Real,allocatable,Dimension(:) :: Rr    !< Work Arrays For Tridiag  !Dim Nn
-  Real,allocatable,Dimension(:) :: Hh    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn) :: Aa    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn) :: Bb    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn) :: Cc    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn) :: Rr    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn) :: Hh    !< Work Arrays For Tridiag  !Dim Nn
 
-  Real,allocatable,Dimension(:) :: Tdot   !< Temperature Time Derivative*Dtt  
-  Real,allocatable,Dimension(:) :: Abis,Bbis,Cbis,Rbis,Hbis                   
-  Real,allocatable,Dimension(:) :: Ee        !< Vertical Coordinate In Ice, Scaled To H Zeta  
+  Real,Dimension(Nz+Nzm) :: Tdot   !< Temperature Time Derivative*Dtt  
+!  Real,allocatable,Dimension(:) :: Abis,Bbis,Cbis,Rbis,Hbis                   
+  Real,Dimension(Nz) :: Ee        !< Vertical Coordinate In Ice, Scaled To H Zeta  
 
   ! Tableaux De Travail 2d 
   ! ___________________________
-  Real,allocatable,Dimension(:,:):: Tbmer    !< Temperature De La Mer A La Base De L'Ice Shelf  
-  Real,allocatable,Dimension(:,:) :: Chalglissx,Chalglissy !< Chaleur De Glissement             
-  Integer,allocatable,Dimension(:,:) :: Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n   
+  Real,Dimension(Nx,Ny):: Tbmer    !< Temperature De La Mer A La Base De L'Ice Shelf  
+  Real,Dimension(Nx,Ny) :: Chalglissx,Chalglissy !< Chaleur De Glissement             
+  Integer,Dimension(Nx,Ny) :: Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n   
   ! Pour L'Advection        
 
   ! Tableaux De Travail 3d
   ! __________________________
-  Real,allocatable,Dimension(:,:,:,:) :: Chalx, Chaly         !Dim Nx,Ny,Nz,N1poly:N2poly
+  Real,Dimension(Nx,Ny,Nz,size(Btt,4)) :: Chalx, Chaly         !Dim Nx,Ny,Nz,N1poly:N2poly
   ! Utilise Pour Calculer La Chaleur De Deformation Selon Xx Yy Zz Et Xy
-  Real,allocatable,Dimension(:,:,:) :: Ffx,Ffy                                       
-  Real,allocatable,Dimension(:,:,:) :: T3d_new                      
-  Real,allocatable,Dimension(:,:,:) :: Chal2_x, Chal2_y, Chal2_z, Chal2_xy
-  Real,allocatable,Dimension(:,:,:) :: Chaldef_maj   !< Chaleur De Deformation                   
-  Real,allocatable,Dimension(:,:,:) :: Advecx,Advecy,Advec
+  Real,Dimension(Nx,Ny,size(Btt,4)) :: Ffx,Ffy                                       
+  Real,Dimension(Nx,Ny,Nz+Nzm) :: T3d_new                      
+  Real,Dimension(Nx,Ny,Nz) :: Chal2_x, Chal2_y, Chal2_z, Chal2_xy
+  Real,Dimension(Nx,Ny,Nz) :: Chaldef_maj   !< Chaleur De Deformation                   
+  Real,Dimension(Nx,Ny,Nz) :: Advecx,Advecy,Advec
   Real,Dimension(Nx,Ny,Nz) :: Cp            !< Specific Heat Capacity (J/(M-3)/K)=Ro Cp  
   Real,Dimension(Nx,Ny,Nz) :: Ct            !< Thermal Conductivity (J/M/K/A)            
Index: /trunk/SOURCES/Temperature-routines/icetemp_mod.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/icetemp_mod.f90	(revision 72)
+++ /trunk/SOURCES/Temperature-routines/icetemp_mod.f90	(revision 73)
@@ -66,48 +66,21 @@
 
   subroutine icetemp
-
+    !$ use OMP_LIB 
     use Icetemp_declar
 
 
     Implicit None 
-
-    ! Allocation de la memoire pour les tableaux dans icetemp_declar
-
-    If (.Not. Allocated(Aa)) then 
-       ! tab Nn
-       Allocate (Aa(Nn),Bb(Nn),Cc(Nn),Rr(Nn),Hh(Nn))  
-       ! tab NzNzm_m
-       Allocate(Tdot(Nz+Nzm))   
-       ! tab Nzm+1
-       Allocate(Abis(Nzm+1),Bbis(Nzm+1),Cbis(Nzm+1),Rbis(Nzm+1),Hbis(Nzm+1))                  
-       ! tab nz_m
-       Allocate(Ee(Nz))        
-       !Tab Nx,Ny
-       Allocate(Tbmer(Nx,Ny),&   
-            Chalglissx(Nx,Ny),&
-            Chalglissy(Nx,Ny),&
-            Iadvec_w(Nx,Ny),&
-            Iadvec_e(Nx,Ny),&
-            Iadvec_s(Nx,Ny),&
-            Iadvec_n(Nx,Ny) )      
-       ! Tab Nx,Ny,Nz
-       Allocate(Chal2_x(Nx,Ny,Nz),&
-            Chal2_y(Nx,Ny,Nz),&
-            Chal2_z(Nx,Ny,Nz),&
-            Chal2_xy(Nx,Ny,Nz),&
-            Chaldef_maj(Nx,Ny,Nz),&
-            Advecx(Nx,Ny,Nz),&
-            Advecy(Nx,Ny,Nz),&
-            Advec(Nx,Ny,Nz))
-       ! Tab Nx,Ny,Nz,N1poly:N2poly
-       Allocate(Chalx(Nx,Ny,Nz,size(Btt,4)),&
-            Chaly(Nx,Ny,Nz,size(Btt,4)))         
-       ! Tab Nx,Ny,N1poly:N2poly
-       Allocate(Ffx( Nx,Ny,size(Btt,4)),&
-            Ffy( Nx,Ny,size(Btt,4)))                                        
-       ! Tab Nx,Ny,NzNzm
-       Allocate( T3d_new (Nx,Ny,Nz+Nzm))                  
-    end If
-
+    
+!~   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)
+    
+    
+	 !$OMP PARALLEL
+	 !$OMP WORKSHARE
     ! init
     Aa=0.    
@@ -124,4 +97,5 @@
     Chalglissx=0.
     Chalglissy=0.
+    
     !Instruction 
     Dee=1./(Nz-1) 
@@ -130,7 +104,12 @@
     Ee(1)=0.
     Ee(Nz)=1.
+    !$OMP END WORKSHARE
+    
+    !$OMP DO
     Do K=1,Nz
        If ((K.Ne.1).And.(K.Ne.Nz)) Ee(K)=(K-1.)/(Nz-1.)
     End Do
+    !$OMP END DO
+    !$OMP END PARALLEL
 
     !   Variables Dependant Du Pas De Temps Dtt
@@ -156,6 +135,8 @@
     Call Qprod_ice(Iq)
 
-
+  
     ! Dans Le Socle : Les elements de La Matrice Tridiag Sont Invariants Dans L'Espace
+    !$OMP PARALLEL
+    !$OMP DO
     Do K=Nz+1,nz+nzm-1
        Aa(K)=-Ctm
@@ -163,4 +144,5 @@
        Cc(K)=-Ctm
     End Do
+	 !$OMP END DO
 
     ! Conditions Aux Limites Pour Le Socle
@@ -174,5 +156,5 @@
     !         Advection Selon X
     !         -----------------
-
+    !$OMP WORKSHARE
     Iadvec_w(:,:)=Nint(Max(Sign(1.0,Uxbar(:,:)),0.0))
     Iadvec_e(:,:)=Nint(Abs(Min(Sign(1.0,Uxbar(:,:)),0.0)))
@@ -183,15 +165,19 @@
     Iadvec_s(:,:)=Nint(Max(Sign(1.0,Uybar(:,:)),0.0)) 
     Iadvec_n(:,:)=Nint(Abs(Min(Sign(1.0,Uybar(:,:)),0.0)))
-
+    !$OMP END WORKSHARE
     ! Boucle De Resolution Des Temperatures Colonne Par Colonne
-
+	 !$OMP SINGLE
     If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  boucle Advec_icetemp'
-
+    !$OMP END SINGLE
+
+	 !$OMP DO PRIVATE(Deh22,Duu,Aa,Bb,Cc,Rr,Hh,Dou)
     Do J=2,Ny-1
        Do I=2,Nx-1
 
  !  If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Advect_icetemp  i,j', i,j
-          call Advec_icetemp(I,J,Chaldef_maj(I,J,:),Cp(I,J,:),& 
-               Advecx(I,J,:),Advecy(I,J,:),Advec(I,J,:),Ct(I,J,:))
+          call Advec_icetemp(Nz,Nzm,Nn,Chaldef_maj(I,J,:),Cp(I,J,:),& 
+               Advecx(I,J,:),Advecy(I,J,:),Advec(I,J,:),Ct(I,J,:),Iadvec_w(I,J),Iadvec_e(I+1,J),Iadvec_s(I,J),Iadvec_n(I,J+1), &
+					Ux(i,j,:),Ux(i+1,j,:),Uy(i,j,:),Uy(i,j+1,:),T(i,j,:),T(i-1,j,:),T(i+1,j,:),T(i,j-1,:),T(i,j+1,:), &
+					H(i,j),Ts(i,j),Aa,Bb,Cc,Rr,Dx11,Dou,DTT,Dee,Uzr(i,j,:),Dttdx)
 
           ! variables de calcul dans la glace
@@ -201,6 +187,7 @@
 
 !   If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  temp_col  i,j', i,j
-          Call temp_col(I,J,T3d_new,Ct(I,J,:),&
-               flot(I,J),H(I,J),Tbmer(I,J))
+          Call temp_col(Nz,Nzm,Nn,T3d_new(i,j,:),Ct(i,j,:),&
+               flot(i,j),H(i,j),Tbmer(i,j),Ibase(i,j),T(i,j,:),Tpmp(i,j,:),Aa,Bb,Cc,Rr,Hh,Ncond,Dee,Cm,&
+               Dzm,Phid(i,j),Ghf(i,j),Ifail1,Tbdot(i,j),DTT,Dou,Ts(i,j),Bmelt(i,j))
 
           ! Test Permettant D'Ecrire Les Variables Quand Il Y A Une Erreur
@@ -218,5 +205,7 @@
        End Do
     End Do
-
+	 !$OMP END DO
+
+	 !$OMP DO
     Do J=1,Ny
        Do I=1,Nx
@@ -226,8 +215,9 @@
        End Do
     End Do
+    !$OMP END DO
 
     ! Attribution Finale De La Temperature
     ! Limitation a 3 deg De Variation Par Dtt Et Pas Plus Froid Que -70
-
+	 !$OMP DO COLLAPSE(2) PRIVATE(Tdot)
     Do K=1,nz+nzm
        Do J=1,Ny
@@ -250,6 +240,7 @@
        End Do
     End Do
-
-
+    !$OMP END DO
+
+	 !$OMP DO COLLAPSE(2)
     Do K=2,Nz
        Do J=1,Ny
@@ -262,8 +253,10 @@
        End Do
     End Do
+    !$OMP END DO
 
     !     Temperature Sur Les Bords :
     !     -------------------------
   
+    !$OMP DO COLLAPSE(2)
     Do K=1,Nz
        Do J=1,Ny
@@ -281,6 +274,7 @@
        End Do
     End Do
-   
-
+    !$OMP END DO
+
+    !$OMP DO COLLAPSE(2) 
     Do K=1,Nz
        Do J=1,Ny
@@ -294,9 +288,21 @@
        End Do
     End Do
-
-    Deallocate (Aa,Bb,Cc,Rr,Hh,Tdot,Abis,Bbis,Cbis,Rbis,Hbis,Ee,Tbmer,            & 
-                Chalglissx,Chalglissy,Iadvec_w,Iadvec_e,                          & 
-                Iadvec_s,Iadvec_n,Chal2_x,Chal2_y,Chal2_z,Chal2_xy,Chaldef_maj,   &
-                Advecx,Advecy,Advec,Chalx,Chaly,Ffx,Ffy,T3d_new)        
+    !$OMP END DO
+    !$OMP END PARALLEL
+ 
+!~   ! 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
+
 
     If (Itracebug.Eq.1) Write(num_tracebug,*)' fin routine icetemp'
Index: /trunk/SOURCES/Temperature-routines/temp_col.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/temp_col.f90	(revision 72)
+++ /trunk/SOURCES/Temperature-routines/temp_col.f90	(revision 73)
@@ -14,20 +14,47 @@
 !>
 
-Subroutine  Temp_col(Ii,Jj,Newtempcol_m,Ct_m,Flot_m,H_m,Tbmer_m)
-
-  Use Icetemp_declar
+Subroutine  Temp_col(Nz,Nzm,Nn,Newtempcol,Ctij,Flotij,Hij,Tbmerij,Ibase,Tij,Tpmpij,Aa,Bb,Cc,Rr,Hh,Ncond,Dee,Cm, &
+				Dzm,Phidij,Ghfij,Ifail1,Tbdotij,DTT,Dou,Tsij,Bmeltij)
+
+!  Use Icetemp_declar
   Use Tridiagmod
   
   
   Implicit None
-  integer :: Ii,Jj
-  Real,Dimension(Nx,Ny,Nz+Nzm),Intent(Out)  :: Newtempcol_m        !< Tableau 1d Vert. Des Temperatures En Sortie
-  Real,Dimension(Nz),Intent(In)             :: Ct_m                !< Tableau 1d Vert. Conductivite Thermique
-  Logical,Intent(In)                        :: Flot_m              !< Vrai Si Flot_mtant (Test D'Archimede) 'O'
-  Real,Intent(In)                           :: H_m                 !< Ice Thickness  'O'
-  Real,Intent(In)                           :: Tbmer_m             !< Temperature De La Mer A La Base De L'Ice Shelf
-
-
-  If (H_m.Gt.10.) Then
+  integer,intent(in) :: Nz,Nzm,Nn
+  Real,Dimension(Nz+Nzm),Intent(Out)  :: Newtempcol        !< Tableau 1d Vert. Des Temperatures En Sortie
+  Real,Dimension(Nz),Intent(In)       :: Ctij                !< Tableau 1d Vert. Conductivite Thermique
+  Logical,Intent(In)                  :: Flotij              !< Vrai Si Flotijtant (Test D'Archimede) 'O'
+  Real,Intent(In)                     :: Hij                 !< Ice Thickness  'O'
+  Real,Intent(In)                     :: Tbmerij             !< Temperature De La Mer A La Base De L'Ice Shelf
+  integer,intent(inout)               :: Ibase
+  real,dimension(Nz+Nzm),intent(inout) :: Tij             !< Temperature sur la colonne
+  real,dimension(Nz),intent(in)       :: Tpmpij               !< pressure melting point temperature in ice sheet 'o'
+  Real,Dimension(Nn),intent(inout) :: Aa    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn),intent(inout) :: Bb    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn),intent(inout) :: Cc    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn),intent(inout) :: Rr    !< Work Arrays For Tridiag  !Dim Nn
+  Real,Dimension(Nn),intent(inout) :: Hh    !< Work Arrays For Tridiag  !Dim Nn
+  Integer,intent(in)               :: Ncond !< Switch : Ncond=1 Thermal Conductivity In Mantle
+  real, intent(in) :: Dee
+  real,intent(in) :: Cm
+  real,intent(in) :: Dzm
+  real,intent(inout) :: Phidij     !< flux de chaleur lie a la deformation et glissement basal
+  real,intent(in) :: Ghfij      !< flux geothermique
+  integer,intent(out) :: Ifail1 !< Diagnostique erreur dans Tridiag
+  real, intent(out) :: Tbdotij  !< variation in time of basal temperature
+  real, intent(in) :: DTT
+  Real, intent(inout) :: Dou
+  real, intent(in) :: Tsij        !< surface ice temperature  'o'
+  real, intent(out) :: Bmeltij    !< fonte basale
+  
+  
+  
+  ! variables locales
+  integer :: K
+  Real :: Dzi,Tss
+  Real,Dimension(Nzm+1) :: Abis,Bbis,Cbis,Rbis,Hbis   
+
+  If (Hij.Gt.10.) Then
 
      ! Conditions Aux Limites a La Base De La Glace
@@ -35,20 +62,20 @@
 
      ! Base Froide 
-     If ( ((Ibase(Ii,Jj).Eq.1).Or.(Ibase(Ii,Jj).Eq.4) &
-          .Or. ((Ibase(Ii,Jj).Eq.5).And.(T(Ii,Jj,Nz).Lt.Tpmp(Ii,Jj,Nz)))) &
-          .And..Not.Flot_m) Then
-
-        Ibase(Ii,Jj)=1
+     If ( ((Ibase.Eq.1).Or.(Ibase.Eq.4) &
+          .Or. ((Ibase.Eq.5).And.(Tij(Nz).Lt.Tpmpij(Nz)))) &
+          .And..Not.Flotij) Then
+
+        Ibase=1
         Bb(Nz)=1.
 
         If (Ncond.Eq.1) Then ! Avec Socle
-           Dzi=H_m*Dee*Cm/Ct_m(Nz)
+           Dzi=Hij*Dee*Cm/Ctij(Nz)
            Aa(Nz)=-Dzm/(Dzm+Dzi)
            Cc(Nz)=-Dzi/(Dzm+Dzi)
-           Rr(Nz)=H_m*Dee*Phid(Ii,Jj)*Dzm/Ct_m(Nz)/(Dzm+Dzi)
+           Rr(Nz)=Hij*Dee*Phidij*Dzm/Ctij(Nz)/(Dzm+Dzi)
         Else                 ! Sans Socle
            Aa(Nz)=-1.
            Cc(Nz)=0.
-           Rr(Nz)=-(Ghf(Ii,Jj)-Phid(Ii,Jj))/Ct_m(Nz)*H_m*Dee
+           Rr(Nz)=-(Ghfij-Phidij)/Ctij(Nz)*Hij*Dee
         Endif
 
@@ -57,14 +84,14 @@
         ! Base Temperee Ou Shelf
         ! ----------------------               
-        If (Ibase(Ii,Jj).Eq.5)  Ibase(Ii,Jj)=2
-        Ibase(Ii,Jj)=Max(Ibase(Ii,Jj),2)
+        If (Ibase.Eq.5)  Ibase=2
+        Ibase=Max(Ibase,2)
         Aa(Nz)=0.
         Bb(Nz)=1.
         Cc(Nz)=0.
 
-        If (.Not.Flot_m) Then
-           Rr(Nz)=Tpmp(Ii,Jj,Nz)
-        Else
-           Rr(Nz)=Tbmer_m
+        If (.Not.Flotij) Then
+           Rr(Nz)=Tpmpij(Nz)
+        Else
+           Rr(Nz)=Tbmerij
         End If
 
@@ -74,5 +101,5 @@
      If (Ncond.Eq.1) Then ! Avec Socle
         Do K=Nz+1,nz+nzm-1
-           Rr(K)=T(Ii,Jj,K)
+           Rr(K)=Tij(K)
         End Do
         Call Tridiag (Aa,Bb,Cc,Rr,Hh,nz+nzm,Ifail1)
@@ -87,15 +114,15 @@
      If (Ncond.Eq.0) Then
         Do K=Nz+1,nz+nzm
-           Hh(K)=Hh(Nz)-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
+           Hh(K)=Hh(Nz)-Dzm*(K-Nz)*Ghfij/Cm
         End Do
      Endif
 
      Do K=1,nz+nzm
-        Newtempcol_m(Ii,Jj,K)=Hh(K)
+        Newtempcol(K)=Hh(K)
      End Do
-     Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz)-T(Ii,Jj,Nz))/Dtt
+     Tbdotij=(Newtempcol(Nz)-Tij(Nz))/Dtt
 
      !      ------------------------------------------------- Glace Tres Fine (H<10m Ou H=0)
-  Else If (H_m.Le.10.) Then
+  Else If (Hij.Le.10.) Then
 
      !       Pour Eviter Des Problemes Lorsque H Est Inferieur A 10 M
@@ -104,25 +131,25 @@
 
      ! ........................................ Avec Calcul Dans Le Socle
-     If (Ncond.Eq.1.And..Not.Flot_m) Then   
-        If (H_m.Gt.0.) Then
+     If (Ncond.Eq.1.And..Not.Flotij) Then   
+        If (Hij.Gt.0.) Then
            ! Gradient Dans Le Socle
-           Dou=(T(Ii,Jj,Nz+1)-T(Ii,Jj,Nz))/Dzm*Cm     
-           Dou=Dou/Ct_m(Nz)*Dee*H_m
-
-           Tss=Min(0.,Ts(Ii,Jj))
-           Do K=1,Nz
-              Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
+           Dou=(Tij(Nz+1)-Tij(Nz))/Dzm*Cm     
+           Dou=Dou/Ctij(Nz)*Dee*Hij
+
+           Tss=Min(0.,Tsij)
+           Do K=1,Nz
+              Newtempcol(K)=Tss+Dou*(K-1.)
            End Do
 
         Else
            ! Pas De Glace
-           Tss=Ts(Ii,Jj)
-           Do K=1,Nz
-              Newtempcol_m(Ii,Jj,K)=Tss
+           Tss=Tsij
+           Do K=1,Nz
+              Newtempcol(K)=Tss
            End Do
         End If
 
         ! Calcul Dans Le Socle Meme S'Il N'Y A Pas De Glace
-        Rr(Nz)=Newtempcol_m(Ii,Jj,Nz)
+        Rr(Nz)=Newtempcol(Nz)
         Aa(Nz)=0
         Bb(Nz)=1
@@ -130,5 +157,5 @@
 
         Do K=Nz+1,nz+nzm-1
-           Rr(K)=T(Ii,Jj,K)
+           Rr(K)=Tij(K)
         End Do
 
@@ -150,35 +177,35 @@
      Else
         ! Calotte Posee                    
-        If ((H_m.Gt.0.).And..Not.Flot_m) Then
-           Dou=-Ghf(Ii,Jj)/Ct_m(Nz)*Dee*H_m
-           Tss=Min(0.,Ts(Ii,Jj))
-           Do K=1,Nz
-              Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
+        If ((Hij.Gt.0.).And..Not.Flotij) Then
+           Dou=-Ghfij/Ctij(Nz)*Dee*Hij
+           Tss=Min(0.,Tsij)
+           Do K=1,Nz
+              Newtempcol(K)=Tss+Dou*(K-1.)
            End Do
 
            !  Shelf
-        Else If ((H_m.Gt.0.).And.Flot_m) Then
-           Tss=Min(0.,Ts(Ii,Jj))
-           Dou=(Tbmer_m-Tss)*Dee
-           Do K=1,Nz
-              Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
+        Else If ((Hij.Gt.0.).And.Flotij) Then
+           Tss=Min(0.,Tsij)
+           Dou=(Tbmerij-Tss)*Dee
+           Do K=1,Nz
+              Newtempcol(K)=Tss+Dou*(K-1.)
            End Do
 
            ! Pas De Glace
         Else
-           Tss=Ts(Ii,Jj)      
-           Do K=1,Nz
-              Newtempcol_m(Ii,Jj,K)=Tss
+           Tss=Tsij      
+           Do K=1,Nz
+              Newtempcol(K)=Tss
            End Do
         End If
 
         ! Temperature Dans Le Socle, Lineaire Avec Le Gradient Ghf
-        If (.Not.Flot_m) Then
+        If (.Not.Flotij) Then
            Do K=Nz+1,nz+nzm
-              Hh(K)=Newtempcol_m(Ii,Jj,Nz)-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
+              Hh(K)=Newtempcol(Nz)-Dzm*(K-Nz)*Ghfij/Cm
            End Do
         Else
            Do K=Nz+1,nz+nzm
-              Hh(K)=Tbmer_m-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
+              Hh(K)=Tbmerij-Dzm*(K-Nz)*Ghfij/Cm
            End Do
         End If
@@ -189,11 +216,11 @@
 
      Do K=Nz+1,nz+nzm
-        Newtempcol_m(Ii,Jj,K)=Hh(K)
+        Newtempcol(K)=Hh(K)
      End Do
 
-     Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz)-T(Ii,Jj,Nz))/Dtt
-     Bmelt(Ii,Jj)=0.
-     Ibase(Ii,Jj)=5
-     Phid(Ii,Jj)=0.
+     Tbdotij=(Newtempcol(Nz)-Tij(Nz))/Dtt
+     Bmeltij=0.
+     Ibase=5
+     Phidij=0.
 
 
Index: /trunk/SOURCES/deformation_mod_2lois.f90
===================================================================
--- /trunk/SOURCES/deformation_mod_2lois.f90	(revision 72)
+++ /trunk/SOURCES/deformation_mod_2lois.f90	(revision 73)
@@ -468,19 +468,4 @@
 !debug_3d(:,:,90+iiglen-1) = s2a(:,:,1,iiglen)
 
-
-!~   ! 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
-
 end subroutine flowlaw
 
Index: /trunk/SOURCES/flottab2-0.7.f90
===================================================================
--- /trunk/SOURCES/flottab2-0.7.f90	(revision 72)
+++ /trunk/SOURCES/flottab2-0.7.f90	(revision 73)
@@ -105,12 +105,4 @@
 
 
-!~   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)
-
 
 if (itracebug.eq.1)  call tracebug(' Entree dans routine flottab')
@@ -725,18 +717,4 @@
 !fin de routine flottab2
 !print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j)
-
-!~   ! 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
 
 end subroutine flottab
Index: /trunk/SOURCES/resol_adv_diff_2D-sept2009.f90
===================================================================
--- /trunk/SOURCES/resol_adv_diff_2D-sept2009.f90	(revision 72)
+++ /trunk/SOURCES/resol_adv_diff_2D-sept2009.f90	(revision 73)
@@ -101,4 +101,5 @@
 subroutine resol_adv_diff_2D_vect(Dfx,Dfy,advx,advy,H_presc,i_Hpresc,bil,vieuxH,newH)
 
+!$ USE OMP_LIB
 implicit none
 
@@ -148,5 +149,6 @@
 if (itracebug.eq.1)  call tracebug(' Entree dans routine resolution_diffusion')
 
-
+!$OMP PARALLEL
+!$OMP WORKSHARE
 ! calcul de bdx et hdx
 hdx(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0.,dim=1))
@@ -164,10 +166,10 @@
 Frelax(:,:)=0.
 DeltaH(:,:)=0.
-
+!$OMP END WORKSHARE
 
 ! schema spatial
 
 if (upwind.eq.1.) then                 !schema amont
-
+!$OMP WORKSHARE
    where (Advx(:,:).ge.0.)
       c_west(:,:)=1.
@@ -185,13 +187,17 @@
       c_north(:,:)=1.
    end where
-
+!$OMP END WORKSHARE
 else if (upwind.lt.1.) then             ! schema centre
+!$OMP WORKSHARE
       c_west(:,:)=0.5
       c_east(:,:)=0.5
       c_south(:,:)=0.5
       c_north(:,:)=0.5
+!$OMP END WORKSHARE      
 end if
 
+
 ! attribution des elements des diagonales
+!$OMP DO PRIVATE(frdx,frdy,fraxw,fraxe,frays,frayn)
 do j=2,ny-1
   do i=2,nx-1
@@ -254,6 +260,7 @@
      end do
   end do
-
-
+!$OMP END DO
+
+!$OMP WORKSHARE
 where (i_hpresc(:,:) .eq.1)          ! thickness prescribed
    frelax(:,:) = H_presc(:,:)
@@ -264,14 +271,14 @@
    erelax(:,:) = 0.
 end where
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 stopp = .false.
 ntour=0
 
-
-
 relax_loop: do while(.not.stopp)
 ntour=ntour+1
-
-
+!$OMP PARALLEL
+!$OMP DO PRIVATE(reste)
 do j=2,ny-1
    do i=2,nx-1
@@ -281,7 +288,7 @@
            + crelax(i,j)*newH(i,j))- frelax(i,j)
 
-if (ntour.eq.1)  debug_3D(i,j,49)=(((arelax(i,j)*newH(i-1,j) +drelax(i,j)*newH(i,j-1)) &
-           + (brelax(i,j)*newH(i+1,j) + erelax(i,j)*newH(i,j+1))) &
-           + crelax(i,j)*newH(i,j))
+!if (ntour.eq.1)  debug_3D(i,j,49)=(((arelax(i,j)*newH(i-1,j) +drelax(i,j)*newH(i,j-1)) &
+!           + (brelax(i,j)*newH(i+1,j) + erelax(i,j)*newH(i,j+1))) &
+!           + crelax(i,j)*newH(i,j))
 
 
@@ -290,17 +297,18 @@
    end do
 end do
-
-debug_3D(:,:,50)=arelax(:,:)
-debug_3D(:,:,51)=brelax(:,:)
-debug_3D(:,:,52)=crelax(:,:)
-debug_3D(:,:,53)=drelax(:,:)
-debug_3D(:,:,54)=erelax(:,:)
-debug_3D(:,:,55)=frelax(:,:)
-
-
-
+!$OMP END DO
+
+!debug_3D(:,:,50)=arelax(:,:)
+!debug_3D(:,:,51)=brelax(:,:)
+!debug_3D(:,:,52)=crelax(:,:)
+!debug_3D(:,:,53)=drelax(:,:)
+!debug_3D(:,:,54)=erelax(:,:)
+!debug_3D(:,:,55)=frelax(:,:)
+
+
+!$OMP WORKSHARE
 newH(:,:)=newH(:,:)-deltaH(:,:)
-
-
+!$OMP END WORKSHARE
+!$OMP END PARALLEL
 
 
@@ -310,5 +318,6 @@
 delh=0
 
-
+!$OMP PARALLEL
+!$OMP DO REDUCTION(+:delh)
 do j=2,ny-1
    do i=2,nx-1
@@ -316,4 +325,6 @@
    end do
 end do
+!$OMP END DO
+!$OMP END PARALLEL
 
 if (delh.gt.0.) then
@@ -330,16 +341,18 @@
 
 ! thickness at the upwind node
-do j = 1, ny
-   do i = 2, nx
-      debug_3D(i,j,92) = c_west(i,j) * newH(i-1,j) + c_east(i,j) * newH(i,j)
-   end do
-end do
-do j = 2, ny
-   do i = 1, nx
-      debug_3D(i,j,93) = c_south(i,j) * newH(i,j-1) + c_north(i,j) * newH(i,j)
-   end do
-end do
-
- if (itracebug.eq.1)  call tracebug(' fin routine resolution_diffusion')
+!do j = 1, ny
+!   do i = 2, nx
+!      debug_3D(i,j,92) = c_west(i,j) * newH(i-1,j) + c_east(i,j) * newH(i,j)
+!   end do
+!end do
+!do j = 2, ny
+!   do i = 1, nx
+!      debug_3D(i,j,93) = c_south(i,j) * newH(i,j-1) + c_north(i,j) * newH(i,j)
+!   end do
+!end do
+
+if (itracebug.eq.1)  call tracebug(' fin routine resolution_diffusion')
+
+ 
 return
 
