Index: /trunk/SOURCES/Makefile.grisli.inc
===================================================================
--- /trunk/SOURCES/Makefile.grisli.inc	(revision 23)
+++ /trunk/SOURCES/Makefile.grisli.inc	(revision 24)
@@ -64,10 +64,9 @@
 	massb_perturb_mois.o massb_perturb_mois_lapsecouche.o 
 
-mod_communs =   util_recovery.o \
-	geom_type.o temperature_type.o ice_flow_type.o mask_flgz_type.o \
-	interface_input.o deformation_type.o autre_pr_temp_type.o  \
+mod_communs = util_recovery.o \
+	interface_input.o \
 	iso_declar_mod-0.3.o isostasie_mod-0.3.o noisostasie_mod-0.3.o \
-        deformation_mod_2lois.o deformation_mod_2lois_isotherme.o     \
-        sliding_Bindshadler_mod.o no_sliding.o			    \
+	deformation_mod_2lois.o deformation_mod_2lois_isotherme.o     \
+	sliding_Bindshadler_mod.o no_sliding.o			    \
 	dragging_hwatermax_0.2_mod.o dragging_calc_beta_mod.o  \
 	dragging_hwat-contigu_mod.o dragging_hwat_contmaj_mod.o \
@@ -82,7 +81,9 @@
 	tridiag_mod-0.3.o out_horiz_mod.o bmelt-seuil-profondeur_mod.o 
 
+	#    geom_type.o temperature_type.o ice_flow_type.o mask_flgz_type.o \
+	#deformation_type.o autre_pr_temp_type.o  \
+	
 mod_communs_RGL =   util_recovery.o \
-        geom_type.o temperature_type.o ice_flow_type.o mask_flgz_type.o \
-        interface_input.o deformation_type.o autre_pr_temp_type.o  \
+        interface_input.o \
         iso_declar_mod-0.3.o isostasie_mod-0.3.o noisostasie_mod-0.3.o \
         deformation_mod_2lois.o deformation_mod_2lois_isotherme.o      \
@@ -100,9 +101,5 @@
         tridiag_mod-0.3.o out_horiz_mod.o bmelt-seuil-profondeur_mod.o
 
-
-
-# Modules communs Post geographiques 
-
-mod_post_geo = step_temp_type.o
+#        geom_type.o temperature_type.o ice_flow_type.o mask_flgz_type.o deformation_type.o autre_pr_temp_type.o \
 
 
@@ -133,13 +130,13 @@
 # enleve  : flow_general-0.3.o flowlaw-0.3.o  Cat 23 mars 2012
 
-routines_communes = tracebug.o ablation_bord.o allocate_types.o \
-	icetemp_declar-0.5_mod.o bmelt-grounded.o calceps2-0.2.o  \
+routines_communes = tracebug.o ablation_bord.o \
+	icetemp_declar_mod.o bmelt-grounded.o calceps2-0.2.o  \
 	diffusiv-polyn-0.6.o \
 	flottab-rescue-1.o flottab2-0.7.o   \
-	interface_advec_icetemp.o \
-	interface_prop_th_icetemp.o \
-	interface_Qprod_icetemp.o \
-	interface_temp_col.o \
-	interface_icetempmod.o \
+	advec_icetemp.o \
+	prop_th_icetemp.o \
+	Qprod_icetemp.o \
+	temp_col.o \
+	icetemp_mod.o \
 	initial-0.3.o initial-phy-2.o initial2-0.4.o \
 	lastoutput-0.2.o lect-eis.o  readinput.o lect_datfile.o write_datfile.o \
@@ -151,5 +148,7 @@
 	tab-litho-0.3.o taubed-0.3.o \
 	velocities-polyn-0.3.o \
-	vtoutput-0.2.o  courbures.o 
+	vtoutput-0.2.o  courbures.o
+
+# allocate_types.o \
 #	vtoutput-0.2.o steps_time_loop.o courbures.o 
 
@@ -471,33 +470,6 @@
 	$(FT) 3D-physique-gen_mod.f90
 
-# types Hassine
 param_phy_mod.o: param_phy_mod.f90
 	$(FT) param_phy_mod.f90
-
-geom_type.o: Types/geom_type.f90
-	$(FT) Types/geom_type.f90
-
-temperature_type.o: Types/temperature_type.f90
-	$(FT) Types/temperature_type.f90
-
-ice_flow_type.o: Types/ice_flow_type.f90
-	$(FT) Types/ice_flow_type.f90
-
-mask_flgz_type.o: Types/mask_flgz_type.f90
-	$(FT) Types/mask_flgz_type.f90
-
-deformation_type.o: Types/deformation_type.f90
-	$(FT) Types/deformation_type.f90
-
-autre_pr_temp_type.o: Types/autre_pr_temp_type.f90
-	$(FT) Types/autre_pr_temp_type.f90
-
-step_temp_type.o: Types/step_temp_type.f90
-	$(FT) Types/step_temp_type.f90
-
-allocate_types.o:  Types/allocate_types.f90
-	$(FT) Types/allocate_types.f90
-
-# fin types Hassine
 
 lect-anteis_mod.o : Ant40_files/lect-anteis_mod.f90
@@ -1320,23 +1292,23 @@
 
 
-# Routines Temperature Hassine
-
-icetemp_declar-0.5_mod.o : Temperature-routines/icetemp_declar-0.5_mod.f90
-	$(FT) Temperature-routines/icetemp_declar-0.5_mod.f90
-
-interface_advec_icetemp.o :  Temperature-routines/interface_advec_icetemp.f90
-	$(FT) Temperature-routines/interface_advec_icetemp.f90
-
-interface_icetempmod.o :  Temperature-routines/interface_icetempmod.f90
-	$(FT) Temperature-routines/interface_icetempmod.f90
-
-interface_prop_th_icetemp.o:  Temperature-routines/interface_prop_th_icetemp.f90
-	$(FT) Temperature-routines/interface_prop_th_icetemp.f90
-
-interface_Qprod_icetemp.o:  Temperature-routines/interface_Qprod_icetemp.f90
-	$(FT) Temperature-routines/interface_Qprod_icetemp.f90
-
-interface_temp_col.o :  Temperature-routines/interface_temp_col.f90
-	$(FT) Temperature-routines/interface_temp_col.f90
+# Routines Temperature
+
+icetemp_declar_mod.o : Temperature-routines/icetemp_declar_mod.f90
+	$(FT) Temperature-routines/icetemp_declar_mod.f90
+
+advec_icetemp.o :  Temperature-routines/advec_icetemp.f90
+	$(FT) Temperature-routines/advec_icetemp.f90
+
+icetemp_mod.o :  Temperature-routines/icetemp_mod.f90
+	$(FT) Temperature-routines/icetemp_mod.f90
+
+prop_th_icetemp.o:  Temperature-routines/prop_th_icetemp.f90
+	$(FT) Temperature-routines/prop_th_icetemp.f90
+
+Qprod_icetemp.o:  Temperature-routines/Qprod_icetemp.f90
+	$(FT) Temperature-routines/Qprod_icetemp.f90
+
+temp_col.o :  Temperature-routines/temp_col.f90
+	$(FT) Temperature-routines/temp_col.f90
 
 # Routines Agrif
@@ -1556,5 +1528,4 @@
 	$(mod_tracers) \
 	$(mod_ell) $(Liste_ANT15-LBq) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1568,5 +1539,4 @@
 	$(mod_ell) \
 	$(Liste_ANT15-LBq) \
-	$(mod_post_geo) \
 	$(diagnoshelf) $(Liste_Netcdf) \
 	$(routines_communes) steps_time_loop.o $(routine_elliptiques) $(NCDF_LIB)  $(MKL_LIB)
@@ -1578,5 +1548,4 @@
 	$(mod_tracers) \
 	$(mod_ell) $(Liste_ANT15-LBq) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1590,5 +1559,4 @@
 	$(mod_ell) \
 	$(Liste_ANT15-LBq) \
-	$(mod_post_geo) \
 	$(diagnoshelf) $(Liste_Netcdf) \
 	$(routines_communes) steps_time_loop.o $(routine_elliptiques) $(NCDF_LIB)  $(MKL_LIB)
@@ -1613,5 +1581,4 @@
 	$(mod_tracers) \
 	$(mod_ell) $(Liste_GrIce2sea) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1625,5 +1592,4 @@
 	$(mod_tracers) \
 	$(mod_ell) $(Liste_GrIce2sea) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1637,5 +1603,4 @@
 	$(mod_tracers) \
 	$(mod_ell) $(Liste_GrIce2sea) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1650,5 +1615,4 @@
 	$(mod_tracers) \
 	$(mod_ell) $(Liste_GrIce2sea) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1662,5 +1626,4 @@
 	$(mod_no_tracers) \
 	$(mod_ell) $(Liste_hemin40) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1676,5 +1639,4 @@
 	$(mod_no_tracers) \
 	$(mod_ell) $(Liste_hemin40) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1688,5 +1650,4 @@
 	$(mod_no_tracers) \
 	$(mod_ell) $(Liste_hemin15) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1701,5 +1662,4 @@
 	$(mod_no_tracers) \
 	$(mod_ell) $(Liste_hemin15) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1713,5 +1673,4 @@
 	$(mod_no_tracers) \
 	$(mod_ell) $(Liste_greeneem15) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
@@ -1726,5 +1685,4 @@
 	$(mod_no_tracers) \
 	$(mod_ell) $(Liste_greeneem15) \
-	$(mod_post_geo) \
 	$(diagnoshelf) \
 	$(Liste_Netcdf) \
Index: /trunk/SOURCES/Temperature-routines/Qprod_icetemp.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/Qprod_icetemp.f90	(revision 24)
+++ /trunk/SOURCES/Temperature-routines/Qprod_icetemp.f90	(revision 24)
@@ -0,0 +1,907 @@
+!> \file Qprod_icetemp.f90
+!! Subroutines for heat production 
+!<
+
+!> \namespace Qprod_icetemp
+!! Subroutines for heat production 
+!<
+
+
+!> SUBROUTINE: Qprod 
+!! Subroutine for heat pruduction
+!! Used modules:
+!!    - use Icetemp_declar
+!>
+
+subroutine Qprod_ice(Iq1)
+
+  use Icetemp_declar
+
+  implicit none
+
+  !<Arguments
+
+  integer,intent(In)          :: Iq1           !< Choice Of Method For Heat Production
+
+  !< Local variables needed by the various case of Iq1
+
+  real, dimension(Nx,Ny)      :: Pente2_maj    !< Pente Au Carre Sur Les Noeuds Majeurs
+  real, dimension(Nx,Ny)      :: Pente_maj     !< Pente Au Carre Sur Les Noeuds Majeurs
+  real, dimension(Nx,Ny)      :: Vit_maj       !< Ubar Moyennee Sur Les Noeuds Majeurs
+  real, dimension(Nx,Ny)      :: Uslid_maj     !< Uslid Moyennee Sur Les Noeuds Majeurs
+  real, dimension(Nx,Ny)      :: Vit_stag2     !< Ubar Moyennee Sur Les Noeuds Stag
+  real, dimension(Nx,Ny,Nz)   :: Vit_stag3     !< Ubar Moyennee Sur Les Noeuds Stag
+  real, dimension(Nx,Ny)      :: Hmxy          !< Epaisseur Moyenne Sur Les Noeuds Stag
+  real, dimension(Nx,Ny)      :: Qslid         !< Chaleur Produite Par Glissement Sur Les Noeuds Stag
+  real, dimension(Nx,Ny)      :: Qdef2         !< Chaleur Produite Par Deformation Sur Les Noeuds Stag
+  real, dimension(Nx,Ny,Nz)   :: Qdef3         !< Chaleur Produite Par Deformation Sur Les Noeuds Stag
+  real, dimension(Nx,Ny)      :: Pente_stag    !< Pente Au Carre Sur Les Noeuds Stag
+  real, dimension(Nx,Ny)      :: Uslid_stag    !< Uslid Moyennee Sur Les Noeuds Stag
+  real, dimension(Nx,Ny)      :: Tob_stag      !< Contrainte Basale
+  real, dimension(Nx,Ny)      :: Tox           !< Contraintes Sur Maille Mx
+  real, dimension(Nx,Ny)      :: Toy           !< Contraintes Sur Maille Mx
+
+
+  select case (Iq1)
+
+  case (1)                ! iq1 = 1      Q_prod_demi, Ancienne Methode Sur Les Demi-Mailles
+     !------------------------------------------------------------------------------------------
+
+     ! Calcul Avec L'Ancienne Methode Sur Les Demi-Mailles
+     ! Nom Du Sub Routine Avant Et Var Utilisee  
+     ! Q_prod_demi(T_m%T,T_m%Tpmp,Deform_m%Tobmx,Deform_m%Tobmy,Geom_m%H,Geom_m%Hmx,Geom_m%Hmy,mask_flot_m%Flot,mask_flot_m%Flotmx,mask_flot_m%Flotmy,mask_flot_m%Flgzmx ,mask_flot_m%Flgzmy,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Ddx,Ice_flow_m%Ddy,&
+     !      Ice_flow_m%Ddbx,Ice_flow_m%Ddby,Deform_m%Epsxx,Deform_m%Epsyy,Deform_m%Epsxy,mask_flot_m%Gzmx,mask_flot_m%Gzmy,Deform_m%Btt,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid,Deform_m%Glen,Deform_m%Visc)
+
+     do K=2,Nz
+        do J=2,Ny-1
+           do I=2,Nx-1
+
+              ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
+              ! Pour Les Ice-Streams Et Ice-Shelves
+              ! Les Divers Eps Sont Calcules Dans Strain-Rate Pour L'Ensemble De La Grille
+
+              if ((flot(I,J).or.flgzmx(I,J).or.flgzmx(I+1,J)).or.  &
+                   (flgzmy(I,J).or.flgzmy(I,J+1))) then
+
+                 Chal2_x(I,J,K)=2.*Visc(I,J,K)*Epsxx(I,J)**2
+                 Chal2_y(I,J,K)=2.*Visc(I,J,K)*Epsyy(I,J)**2
+                 Chal2_z(I,J,K)=2.*Visc(I,J,K)*(-Epsxx(I,J)-Epsyy(I,J))**2
+
+                 !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
+                 Chal2_xy(I,J,K)=(Epsxy(I,J)+Epsxy(I+1,J)+Epsxy(I+1,J+1)+Epsxy(I,J+1))
+                 Chal2_xy(I,J,K)=2.*Visc(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
+
+              else   ! Glace Posee
+                 Chal2_x(I,J,K)=0.
+                 Chal2_y(I,J,K)=0.
+                 Chal2_z(I,J,K)=0.
+                 Chal2_xy(I,J,K)=0.
+              endif
+           end do
+        end do
+     end do
+
+     !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
+     do L=1,size(Btt,4) !N1poly,N2poly
+        do J=2,Ny
+           do I=2,Nx
+              !          Ffx A 3 Dimensions !
+              Ffx(I,J,L)=ddx(I,J,L)*Sdx(I,J)*Sdx(I,J)
+              Ffy(I,J,L)=ddy(I,J,L)*Sdy(I,J)*Sdy(I,J)
+           end do
+        end do
+     end do
+
+     do L=1,size(Btt,4)!N1poly,N2poly
+        do K=2,Nz
+           do J=2,Ny
+              do I=2,Nx
+                 if ((.not.Flotmx(I,J)).and.(.not.Gzmx(I,J))) then
+                    Chalx(I,J,K,L)=(Btt(I-1,J,K,L)+Btt(I,J,K,L))*Ffx(I,J,L) !&
+                    !                      *Ro*G*Ee(K)**(Glen(L)+1)/Cp(I,J,K)
+
+                 else if (Gzmx(I,J)) then ! Ice Streams
+                    Chalx(I,J,K,L)=0.
+                 else                     ! Ice Shelves  
+                    Chalx(I,J,K,L)=0.
+                 endif
+                 if ((.not.Flotmy(I,J)).and.(.not.Gzmy(I,J))) then
+                    Chaly(I,J,K,L)=(Btt(I,J-1,K,L)+Btt(I,J,K,L))*Ffy(I,J,L) !&
+                    !                      *Ro*G*Ee(K)**(Glen(L)+1)/Cp(I,J,K)
+                 else if (Gzmy(I,J)) then  ! Ice Streams
+                    Chaly(I,J,K,L)=0.
+
+                 else                      ! Ice Shelves
+                    Chaly(I,J,K,L)=0.
+                 endif
+
+              end do
+           end do
+        end do
+     end do
+
+     !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
+     !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
+     !         Les Productions X Et Y
+     !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy) 
+
+     do K=2,Nz
+        do J=1,Ny-1
+           do I=1,Nx-1
+
+              ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
+              Chaldef_maj(I,J,K)= 0.
+
+              ! Chalk_2 Pour Ice Shelves Et Ice Streams
+              Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
+                   Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
+
+              ! Chalk_1 Pour La Partie PosÃ©e
+              do L=1,size(Btt,4)!N1poly,N2poly
+
+                 ! On Somme La Chaleur Due Aux Diverses Lois De DÃ©formation Et Celles En X Et Y
+                 Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))+ &
+                      (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
+
+                 Chalk_1=Ro*G*Chalk_1/4.*(Ee(K)**(Glen(L)+1))/Cp(I,J,K)
+                 Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
+              enddo
+
+              ! Pour Shelves Et Streams,  On Ajoute Chalk_2
+              Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
+           end do
+        end do
+     end do
+     ! Chaleur Produite A La Base Par Le Glissement
+     do J=2,Ny
+        do I=2,Nx
+
+           if (Gzmx(I,J)) then
+              Chalglissx(I,J)= abs(Uxbar(I,J)*Tobmx(I,J))
+           else
+              Chalglissx(I,J)=ddbx(I,J)*Sdx(I,J)**2*Ro*G*Hmx(I,J)
+           endif
+
+           if (Gzmy(I,J)) then
+              Chalglissy(I,J)= abs(Uybar(I,J)*Tobmy(I,J))
+           else
+              Chalglissy(I,J)=ddby(I,J)*sdy(I,J)**2*Ro*G*Hmy(I,J)
+           endif
+
+        end do
+     end do
+
+     !  Boundary Condition Ice-Rock Interface 
+     K=Nz
+
+     do J=2,Ny-1
+        do I=2,Nx-1
+
+           ! L'Ancien Chalbed Correspond A Chaldef_maj(I,J,Nz)  : Pas La Peine De Recalculer
+           !----------------------------------------------------
+           !      Chalbed=0.
+           !      Do L=1,size(Btt,4)!N1poly,N2poly
+           !         Chalbed_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L)) &
+           !              + (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
+           !         Chalbed_1=Ro*G*Chalbed_1/4.*H(I,J)
+           !         Chalbed= Chalbed_1 +Chalbed
+           !      Enddo
+
+           !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
+           !        Dans La Derniere 1/2 Maille Et Par Le Glissement
+           !        Attention Phid Est >0 Et Ghf Est <0
+           if (.not.flot(I,J)) then
+
+              !           Phid Avec Fonte Sous Les Streams
+              !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
+              !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
+              !                        Chalbed*Fracq
+
+
+              ! Moyenne Phid Sur 4 Points : Attention 0.5 Pour Moyenne Gauche - Droite 
+              ! Mais Les Chaleurs X,Y S'Ajoutent
+
+              !              Phid(I,J)=0.5*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
+              !                        Chalglissy(I,J)+Chalglissy(I,J+1)) 
+
+              ! Moyenne Phid Sur 12 Points. On Moyenne Chaque Chaleur Et On Somme X Et Y
+
+              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
+                   Chalglissy(I,J)+Chalglissy(I,J+1)) 
+              Phid(I,J)=Phid(I,J)+0.125*( &
+                   ((Chalglissx(I,J-1)+Chalglissx(I+1,J-1))+ &
+                   (Chalglissx(I,J+1)+Chalglissx(I+1,J+1))) + &
+                   ((Chalglissy(I-1,J)+Chalglissy(I-1,J+1)) +  &
+                   (Chalglissy(I+1,J)+Chalglissy(I+1,J+1))))
+
+              ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
+              Phid(I,J)=Phid(I,J)*exp((T(I,J,Nz)-Tpmp(I,J,Nz))*Ecart_phid)
+
+              ! Pour Sorties Heino
+              Chalgliss_maj(I,J)=Phid(I,J)
+
+              Phid(I,J)=Phid(I,J)+ Chaldef_maj(I,J,Nz)*Fracq*H(I,J)*Cp(I,J,Nz)
+
+           else
+              Phid(I,J)=0.
+           endif
+        end do
+     end do
+
+
+  case (2)             ! Q_prod_centre : Calcul Avec La Somme Des Carres
+     !-------------------------------------------------------------------------------------
+     ! Calcul Avec La Somme Des Carres
+     ! Nom Du Sub Routine Avant Et Var Utilisee  
+     ! Q_prod_centre(T_m%T,T_m%Tpmp,Deform_m%Tobmx,Deform_m%Tobmy,Geom_m%H,Geom_m%Hmx,Geom_m%Hmy,mask_flot_m%Flot,mask_flot_m%Flotmx,mask_flot_m%Flotmy,mask_flot_m%Flgzmx,mask_flot_m%Flgzmy,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Ddx,Ice_flow_m%Ddy,&
+     !     Ice_flow_m%Ddbx,Ice_flow_m%Ddby,Deform_m%Epsxx,Deform_m%Epsyy,Deform_m%Epsxy,mask_flot_m%Gzmx,mask_flot_m%Gzmy,Deform_m%Btt,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid,Deform_m%Glen,Deform_m%Visc)
+
+     do K=2,Nz  
+        do J=2,Ny-1
+           do I=2,Nx-1
+
+
+              ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
+              ! Pour Les Ice-Streams Et Ice-Shelves
+              ! Les Divers Eps Sont CalculÃ©s Dans Strain-Rate Pour L'Ensemble De La Grille
+
+              if ((flot(I,J).or.Flgzmx(I,J).or.Flgzmx(I+1,J)).or.  &
+                   (Flgzmy(I,J).or.Flgzmy(I,J+1))) then
+
+                 Chal2_x(I,J,K)=2.*Visc(I,J,K)*Epsxx(I,J)**2
+                 Chal2_y(I,J,K)=2.*Visc(I,J,K)*Epsyy(I,J)**2
+                 Chal2_z(I,J,K)=2.*Visc(I,J,K)*(-Epsxx(I,J)-Epsyy(I,J))**2
+
+                 !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
+                 Chal2_xy(I,J,K)=(Epsxy(I,J)+Epsxy(I+1,J)+Epsxy(I+1,J+1)+Epsxy(I,J+1))
+                 Chal2_xy(I,J,K)=2.*Visc(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
+
+              else   ! Glace Posee
+                 Chal2_x(I,J,K)=0.
+                 Chal2_y(I,J,K)=0.
+                 Chal2_z(I,J,K)=0.
+                 Chal2_xy(I,J,K)=0.
+              endif
+           end do
+        end do
+     end do
+
+     !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
+     do L=1,size(Btt,4)!N1poly,N2poly
+        do J=2,Ny
+           do I=2,Nx
+
+              !          Ffx A 3 Dimensions !
+              Ffx(I,J,L)=ddx(I,J,L)*Sdx(I,J)*Sdx(I,J)
+              Ffy(I,J,L)=ddy(I,J,L)*Sdy(I,J)*Sdy(I,J)
+           end do
+        end do
+     end do
+
+
+     do L=1,size(Btt,4) !N1poly,N2poly
+        do K=2,Nz
+           do J=2,Ny
+              do I=2,Nx
+                 if ((.not.Flotmx(I,J)).and.(.not.Gzmx(I,J))) then
+                    Chalx(I,J,K,L)=(Btt(I-1,J,K,L)+Btt(I,J,K,L))*Ffx(I,J,L) !&
+                    !                      *Ro*G*Ee(K)**(Glen(L)+1)/Cp(I,J,K)
+
+                 else if (Gzmx(I,J)) then ! Ice Streams
+                    Chalx(I,J,K,L)=0.
+
+                 else                     ! Ice Shelves  
+                    Chalx(I,J,K,L)=0.
+
+                 endif
+
+                 if ((.not.Flotmy(I,J)).and.(.not.Gzmy(I,J))) then
+                    Chaly(I,J,K,L)=(Btt(I,J-1,K,L)+Btt(I,J,K,L))*Ffy(I,J,L) !&
+                    !                      *Ro*G*Ee(K)**(Glen(L)+1)/Cp(I,J,K)
+
+                 else if (Gzmy(I,J)) then  ! Ice Streams
+                    Chaly(I,J,K,L)=0.
+
+                 else                      ! Ice Shelves
+                    Chaly(I,J,K,L)=0.
+                 endif
+
+              end do
+           end do
+        end do
+     end do
+
+     !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
+     !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
+     !         Les Productions X Et Y
+     !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy) 
+
+     do K=2,Nz
+        do J=1,Ny-1
+           do I=1,Nx-1
+
+              ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
+              Chaldef_maj(I,J,K)= 0.
+
+              ! Chalk_2 Pour Ice Shelves Et Ice Streams
+              Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
+                   Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
+
+              ! Chalk_1 Pour La Partie PosÃ©e
+              do L=1,size(Btt,4)!N1poly,N2poly
+
+                 ! On Somme La Chaleur Due Aux Diverses Lois De DÃ©formation Et Celles En X Et Y
+                 ! La Somme Se Fait Par (Cx2+Cy2)** 0.5. Le 4 Vient Des Moyennes Deform_m%Btt + Des Moyennes Gauche-Droite
+
+                 Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))**2+ &
+                      (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))**2
+                 Chalk_1=Chalk_1**0.5
+
+                 Chalk_1=Ro*G*Chalk_1/4.*(Ee(K)**(Glen(L)+1))/Cp(I,J,K)
+                 Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
+              enddo
+
+              ! Pour Shelves Et Streams,  On Ajoute Chalk_2
+              Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
+           end do
+        end do
+     end do
+
+     ! Chaleur Produite A La Base Par Le Glissement
+
+     do J=2,Ny
+        do I=2,Nx
+
+           if (Gzmx(I,J)) then
+              Chalglissx(I,J)= abs(Uxbar(I,J)*Tobmx(I,J))
+           else
+              Chalglissx(I,J)=ddbx(I,J)*Sdx(I,J)**2*Ro*G*Hmx(I,J)
+           endif
+
+           if (Gzmy(I,J)) then
+              Chalglissy(I,J)= abs(Uybar(I,J)*Tobmy(I,J))
+           else
+              Chalglissy(I,J)=ddby(I,J)*Hmy(I,J)**2*Ro*G*Hmy(I,J)
+           endif
+
+        end do
+     end do
+
+     !  Boundary Condition Ice-Rock Interface 
+
+     K=Nz
+
+     do J=2,Ny-1
+        do I=2,Nx-1
+
+           !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
+           !        Dans La Derniere 1/2 Maille Et Par Le Glissement
+           !        Attention Phid Est >0 Et Ghf Est <0
+           if (.not.flot(I,J)) then
+
+              !           Phid Avec Fonte Sous Les Streams
+              !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
+              !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
+              !                        Chalbed*Fracq
+
+
+              ! Moyenne Phid Sur 4 Points Formulation (A2+B2)**0.5
+              !
+              Chalgliss_maj(I,J)=(Chalglissx(I,J)+Chalglissx(I+1,J))**2+  &
+                   (Chalglissy(I,J)+Chalglissy(I,J+1))**2
+
+              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)**0.5
+
+              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.5 ! Les Moyennes Droite Gauche
+
+              ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
+              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T(I,J,Nz)-Tpmp(I,J,Nz))*Ecart_phid)
+
+
+              ! Flux Total A Rajouter Ã La Base
+              Phid(I,J)=Chalgliss_maj(I,J)+Chaldef_maj(I,J,Nz)*Fracq*H(I,J)*Cp(I,J,Nz)
+
+
+           else
+              Phid(I,J)=0.
+           endif
+        end do
+     end do
+
+
+  case (3)     ! Q_prod_pente :  Pour Essayer D'Avoir La Chaleur En Alpha4
+     !-------------------------------------------------------------------------------------
+     ! Calcul Pour Essayer D'Avoir La Chaleur En Alpha4
+     ! Nom Du Sub Routine Avant Et Var Utilisee  
+     ! Q_prod_pente(T_m%T,T_m%Tpmp,Deform_m%Tobmx,Deform_m%Tobmy,Geom_m%H,Geom_m%Hmx,Geom_m%Hmy,mask_flot_m%Flot,mask_flot_m%Flotmx,mask_flot_m%Flotmy,mask_flot_m%Flgzmx,mask_flot_m%Flgzmy,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Ddx,Ice_flow_m%Ddy,&
+     !     Ice_flow_m%Ddbx,Ice_flow_m%Ddby,Deform_m%Epsxx,Deform_m%Epsyy,Deform_m%Epsxy,mask_flot_m%Gzmx,mask_flot_m%Gzmy,Deform_m%Btt,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid,Deform_m%Glen,Deform_m%Visc)
+
+
+     do J=2,Ny
+        do I=2,Nx
+           Pente2_maj(I,J)=(Sdx(I,J)+Sdx(I+1,J))**2 + &
+                (Sdy(I,J)+Sdy(I,J+1))**2
+           Pente2_maj(I,J)=Pente2_maj(I,J)*0.25    !Pour La Moyenne Sur Sdx Et Sdy
+        end do
+     end do
+
+     do K=2,Nz
+        do J=2,Ny-1
+           do I=2,Nx-1
+              ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
+              ! Pour Les Ice-Streams Et Ice-Shelves
+              ! Les Divers Eps Sont CalculÃ©s Dans Strain-Rate Pour L'Ensemble De La Grille
+              if ((flot(I,J).or.Flgzmx(I,J).or.Flgzmx(I+1,J)).or.  &
+                   (Flgzmy(I,J).or.Flgzmy(I,J+1))) then
+
+                 Chal2_x(I,J,K)=2.*Visc(I,J,K)*Epsxx(I,J)**2
+                 Chal2_y(I,J,K)=2.*Visc(I,J,K)*Epsyy(I,J)**2
+                 Chal2_z(I,J,K)=2.*Visc(I,J,K)*(-Epsxx(I,J)-Epsyy(I,J))**2
+
+                 !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
+                 Chal2_xy(I,J,K)=(Epsxy(I,J)+Epsxy(I+1,J)+Epsxy(I+1,J+1)+Epsxy(I,J+1))
+                 Chal2_xy(I,J,K)=2.*Visc(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
+
+              else   ! Glace Posee
+                 Chal2_x(I,J,K)=0.
+                 Chal2_y(I,J,K)=0.
+                 Chal2_z(I,J,K)=0.
+                 Chal2_xy(I,J,K)=0.
+              endif
+           end do
+        end do
+     end do
+
+     !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
+     !  Do L=1,size(Deform_m%Btt,4)!N1poly,N2poly
+     !     Do J=2,Ny_m
+     !        Do I=2,Nx_m
+     !
+     !          Ffx A 3 Dimensions !  
+     !           Ffx(I,J,L)=Ddx(I,J,L)*Sdx(I,J)*Sdx(I,J)
+     !           Ffy(I,J,L)=Ddy(I,J,L)*Sdy(I,J)*Sdy(I,J)
+     !        End Do
+     !     End Do
+     !  End Do
+
+     do L=1,size(Btt,4)!N1poly,N2poly
+        do K=2,Nz
+           do J=2,Ny
+              do I=2,Nx
+                 if ((.not.Flotmx(I,J)).and.(.not.Gzmx(I,J))) then
+                    Chalx(I,J,K,L)=(Btt(I-1,J,K,L)+Btt(I,J,K,L))*ddx(I,J,L)
+
+                 else if (Gzmx(I,J)) then ! Ice Streams
+                    Chalx(I,J,K,L)=0.
+
+                 else                     ! Ice Shelves  
+                    Chalx(I,J,K,L)=0.
+
+                 endif
+
+                 if ((.not.Flotmy(I,J)).and.(.not.Gzmy(I,J))) then
+                    Chaly(I,J,K,L)=(Btt(I,J-1,K,L)+Btt(I,J,K,L))*ddy(I,J,L)
+
+                 else if (Gzmy(I,J)) then  ! Ice Streams
+                    Chaly(I,J,K,L)=0.
+
+                 else                      ! Ice Shelves
+                    Chaly(I,J,K,L)=0.
+                 endif
+
+              end do
+           end do
+        end do
+     end do
+
+     !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
+     !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
+     !         Les Productions X Et Y
+     !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy) 
+
+     do K=2,Nz
+        do J=1,Ny-1
+           do I=1,Nx-1
+
+              ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
+              Chaldef_maj(I,J,K)= 0.
+
+              ! Chalk_2 Pour Ice Shelves Et Ice Streams
+              Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
+                   Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
+
+              ! Chalk_1 Pour La Partie PosÃ©e
+              do L=1,size(Btt,4)!N1poly,N2poly
+
+                 ! On Somme La Chaleur Due Aux Diverses Lois De DÃ©formation Et Celles En X Et Y
+                 ! La Somme Se Fait Par (Cx2+Cy2)** 0.5. Le 4 Vient Des Moyennes Deform_m%Btt + Des Moyennes Gauche-Droite
+
+                 ! On Fait La Moyenne Des Termes Ddx*Btt (*0.25 Pour Cette Moyenne, Le 0.5 Est Pour Les Btt)
+                 Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))+ &
+                      (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
+
+                 Chalk_1=Chalk_1*0.25*0.5
+
+                 ! On Multiplie Par La Pente Moyenne Au Carre Sur Le Noeud Majeur
+                 Chalk_1=Chalk_1*Pente2_maj(I,J)
+
+                 Chalk_1=Ro*G*Chalk_1*(Ee(K)**(Glen(L)+1))/Cp(I,J,K)  ! Attention Plus De /4
+                 Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
+              enddo
+
+              ! Pour Shelves Et Streams,  On Ajoute Chalk_2
+              Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
+           end do
+        end do
+     end do
+
+     ! Chaleur Produite A La Base Par Le Glissement
+
+     do J=2,Ny
+        do I=2,Nx
+
+           if (Gzmx(I,J)) then
+              Chalglissx(I,J)= abs(Uxbar(I,J)*Tobmx(I,J))
+           else
+              Chalglissx(I,J)=ddbx(I,J)*Sdx(I,J)**2*Ro*G*Hmx(I,J)
+           endif
+
+           if (Gzmy(I,J)) then
+              Chalglissy(I,J)= abs(Uybar(I,J)*Tobmy(I,J))
+           else
+              Chalglissy(I,J)=ddby(I,J)*Sdy(I,J)**2*Ro*G*Hmy(I,J)
+           endif
+
+        end do
+     end do
+
+
+     !  Boundary Condition Ice-Rock Interface 
+
+     K=Nz
+
+     do J=2,Ny-1
+        do I=2,Nx-1
+
+           !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
+           !        Dans La Derniere 1/2 Maille Et Par Le Glissement
+           !        Attention Phid Est >0 Et Ghf Est <0
+           if (.not.flot(I,J)) then
+
+              !           Phid Avec Fonte Sous Les Streams
+              !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
+              !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
+              !                        Chalbed*Fracq
+
+
+              ! Moyenne Phid Sur 4 Points Formulation (A2+B2)**0.5  ! Je Garde Pour L'Instant A Revoir
+              !
+              Chalgliss_maj(I,J)=(Chalglissx(I,J)+Chalglissx(I+1,J))**2+  &
+                   (Chalglissy(I,J)+Chalglissy(I,J+1))**2
+
+              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)**0.5
+
+              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.5 ! Les Moyennes Droite Gauche
+
+              ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
+              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T(I,J,Nz)-Tpmp(I,J,Nz))*Ecart_phid)
+
+
+              ! Flux Total A Rajouter a La Base
+
+              Phid(I,J)=Chalgliss_maj(I,J)+Chaldef_maj(I,J,Nz)*Fracq*H(I,J)*Cp(I,J,Nz)
+
+
+           else
+              Phid(I,J)=0.
+           endif
+        end do
+     end do
+
+
+  case (4)   ! Q_u_taub : Routine Qui concentre La Chaleur  au Socle
+     !-------------------------------------------------------------------------------------
+     ! Routine Qui Prend La Chaleur Remise Au Socle
+     ! Nom Du Sub Routine Avant Et Var Utilisee  
+     ! Q_u_taub(T_m%T,T_m%Tpmp,Geom_m%H,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid)
+
+
+     do J=2,Ny
+        do I=2,Nx
+
+
+           Pente_maj(I,J)=(Sdx(I,J)+Sdx(I+1,J))**2 + &          ! Pente
+                (Sdy(I,J)+Sdy(I,J+1))**2
+
+           Pente_maj(I,J)=Pente_maj(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
+           Pente_maj(I,J)=Pente_maj(I,J)**0.5
+
+           Vit_maj(I,J)=(Uxbar(I,J)+Uxbar(I+1,J))**2 + &        ! Vitesse De Bilan
+                (Uybar(I,J)+Uybar(I,J+1))**2
+           Vit_maj(I,J)=Vit_maj(I,J)*0.25
+           Vit_maj(I,J)=Vit_maj(I,J)**0.5
+
+           Uslid_maj(I,J)=(Ux(I,J,Nz)+Ux(I+1,J,Nz))**2 + &        ! Vitesse De Bilan
+                (Uy(I,J,Nz)+Uy(I,J+1,Nz))**2
+           Uslid_maj(I,J)=Uslid_maj(I,J)*0.25
+           Uslid_maj(I,J)=Uslid_maj(I,J)**0.5
+
+        end do
+     end do
+
+
+
+     Chaldef_maj(:,:,:)=0.
+     Chalgliss_maj(:,:)=Ro*G*H(:,:)*Pente_maj(:,:)*Uslid_maj(:,:)
+     Chaldef_maj(:,:,Nz)=Ro*G*H(:,:)*Pente_maj(:,:)*Vit_maj(:,:)-Chalgliss_maj(:,:)
+
+     ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
+
+     Chalgliss_maj(:,:)=Chalgliss_maj(:,:)*exp((T(:,:,Nz)-Tpmp(:,:,Nz))*Ecart_phid)
+
+     ! Flux Total A Rajouter a Base
+     Phid(:,:)=Chaldef_maj(:,:,Nz)+Chalgliss_maj(:,:)
+
+
+  case (5)    ! Q_u_taub_stag : concentre La Chaleur Au Socle mais Sur Les Mailles Staggered
+     !----------------------------------------------------------------------------------
+     ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
+     ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
+     ! Nom Du Sub Routine Avant Et Var UtilisÃ©  
+     ! Q_u_taub_stag(T_m%T,T_m%Tpmp,Geom_m%H,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid)
+
+
+     do J=2,Ny
+        do I=2,Nx
+
+           Pente_stag(I,J)=(Sdx(I,J)+Sdx(I,J-1))**2 + &          ! Pente
+                (Sdy(I,J)+Sdy(I-1,J))**2
+
+           Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
+           Pente_stag(I,J)=Pente_stag(I,J)**0.5
+
+           Vit_stag2(I,J)=(Uxbar(I,J)+Uxbar(I,J-1))**2 + &        ! Vitesse De Bilan
+                (Uybar(I,J)+Uybar(I-1,J))**2
+           Vit_stag2(I,J)=Vit_stag2(I,J)*0.25
+           Vit_stag2(I,J)=Vit_stag2(I,J)**0.5
+
+           Uslid_stag(I,J)=(Ux(I,J,Nz)+Ux(I,J-1,Nz))**2 + &        ! Vitesse De Bilan
+                (Uy(I,J,Nz)+Uy(I-1,J,Nz))**2
+           Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
+           Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
+
+           Hmxy(I,J)=((H(I,J)+H(I-1,J-1))+(H(I,J-1)+H(I-1,J)))*0.25
+
+        end do
+     end do
+
+
+     Qslid(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Uslid_stag(:,:)
+     Qdef2(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Vit_stag2(:,:)-Qslid(:,:)
+     Chaldef_maj(:,:,:)=0.
+
+
+     ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
+
+     do J=2,Ny-1
+        do I=2,Nx-1
+           Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
+           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
+
+           Chaldef_maj(I,J,Nz)=(Qdef2(I,J)+Qdef2(I+1,J+1))+(Qdef2(I,J+1)+Qdef2(I+1,J))
+           Chaldef_maj(I,J,Nz)=Chaldef_maj(I,J,Nz)*0.25
+
+           ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
+           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T(I,J,Nz)-Tpmp(I,J,Nz))*Ecart_phid)
+
+        end do
+     end do
+
+
+
+
+     !
+     !Write(126,*)'Time=',Time
+     !J=41
+     !Do I=20,30
+     !   Write(126,'(I3,6(E14.4,1x))') I,Chalgliss_maj(I,J),Chaldef_maj(I,J,Nz_m),Qdef2(I,J),Qdef2(I+1,J+1) &
+     !       , Qdef2(I,J+1),Qdef2(I+1,J)
+     !End Do
+
+     ! Flux Total A Rajouter Ã La Base
+     Phid(:,:)= Chaldef_maj(:,:,Nz)+Chalgliss_maj(:,:)
+
+  case (6)    ! Q_sia_stag  : concentre La Chaleur  Au Socle Mais Sur Les Mailles Staggered
+     ! Au Milieux Des Mailles Vitesses
+
+     !----------------------------------------------------------------------------------
+     ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
+     ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
+     ! Nom Du Sub Routine Avant Et Var UtilisÃ©  
+     ! Q_sia_stag(T_m%T,T_m%Tpmp,Geom_m%H,Geom_m%Slop_x,Geom_m%Slop_y,therm_var_m%Phid)
+
+
+     Pente_stag = 0.
+     Vit_stag3   = 0.
+     Uslid_stag = 0.
+     Hmxy  = 0.
+     Qslid = 0.
+     Qdef3 = 0.
+
+     do J=2,Ny
+        do I=2,Nx
+
+           !     Variables 2d
+           Pente_stag(I,J)=(Sdx(I,J)+Sdx(I,J-1))**2 + &          ! Pente
+                (Sdy(I,J)+Sdy(I-1,J))**2
+
+           Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
+           Pente_stag(I,J)=Pente_stag(I,J)**0.5
+
+           Uslid_stag(I,J)=(Ux(I,J,Nz)+Ux(I,J-1,Nz))**2 + &        ! Vitesse De Bilan
+                (Uy(I,J,Nz)+Uy(I-1,J,Nz))**2
+           Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
+           Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
+
+           Hmxy(I,J)=((H(I,J)+H(I-1,J-1))+(H(I,J-1)+H(I-1,J)))*0.25
+
+           !     En Vertical : Calcul De La Deformation Par Differentiation Des Vitesses
+
+           do K=1,Nz
+
+              Vit_stag3(I,J,K)=(Ux(I,J,K)+Ux(I,J-1,K))**2 + &        ! Magnitude Vitesse 
+                   (Uy(I,J,K)+Uy(I-1,J,K))**2            ! A Tous Niveaux K
+
+              Vit_stag3(I,J,K)=Vit_stag3(I,J,K)*0.25
+              Vit_stag3(I,J,K)=Vit_stag3(I,J,K)**0.5
+           end do
+
+           Qdef3(I,J,1)=0.
+           do K=2,Nz-1
+              Qdef3(I,J,K)=Vit_stag3(I,J,K-1)-Vit_stag3(I,J,K+1)             ! Difference Des Vitesses
+              Qdef3(I,J,K)=Ro*G*Pente_stag(I,J)*Ee(K)*Qdef3(I,J,K)/2./Dee    ! Gamma Tau
+           end do
+
+
+           Qdef3(I,J,Nz)=Vit_stag3(I,J,Nz-1)-Vit_stag3(I,J,Nz)         ! Pour Le Fond Differentiation
+           Qdef3(I,J,Nz)=Ro*G*Pente_stag(I,J)*Qdef3(I,J,K)/Dee             ! Sur Une Seule Maille
+        end do
+     end do
+
+
+     Qslid(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Uslid_stag(:,:)
+
+     ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
+
+     do J=2,Ny-1
+        do I=2,Nx-1
+           Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
+           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
+
+           ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
+           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T(I,J,Nz)-Tpmp(I,J,Nz))*Ecart_phid)
+
+           ! Chaleur Due A La Defromation
+           do K=1,Nz
+              Chaldef_maj(I,J,K)=(Qdef3(I,J,K)+Qdef3(I+1,J+1,K))+(Qdef3(I,J+1,K)+Qdef3(I+1,J,K))
+              Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)*0.25
+              Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)/Cp(I,J,K)
+           end do
+
+
+        end do
+     end do
+
+     !Write(126,*)'Time=',Time
+     !J=41
+     !Do I=20,30
+     !   Write(126,'(I3,6(E14.4,1x))') I,Chalgliss_maj(I,J),Chaldef_maj(I,J,Nz_m),Qdef2(I,J),Qdef2(I+1,J+1) &
+     !       , Qdef2(I,J+1),Qdef2(I+1,J)
+     !End Do
+
+     ! Flux Total A Rajouter Ã La Base
+
+     Phid(:,:)=Chalgliss_maj(:,:)+Chaldef_maj(:,:,Nz)*Fracq*H(:,:)*Cp(:,:,Nz)
+
+  case (7)   !  Q_all_stag idem case 6 mais Le Calcul De Taub Tient Compte Des Noeuds Grzmx
+     !----------------------------------------------------------------------------------
+     ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
+     ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
+     ! Le Calcul De Taub Tient Compte Des Noeuds Grzmx     
+     ! Nom Du Sub Routine Avant Et Var UtilisÃ©  
+     ! Q_all_stag(T_m%T,T_m%Tpmp,Deform_m%Tobmx,Deform_m%Tobmy,Geom_m%H,Geom_m%Hmx,Geom_m%Hmy,mask_flot_m%Flgzmx,mask_flot_m%Flgzmy,Geom_m%Slop_x,Geom_m%Slop_y,therm_var_m%Phid)
+
+     where (Flgzmx(:,:))
+        Tox(:,:)=Tobmx(:,:)
+     elsewhere
+        Tox(:,:)=Sdx(:,:)*Hmx(:,:)*Ro*G
+     end where
+
+     where (Flgzmy(:,:))
+        Toy(:,:)=Tobmy(:,:)
+     elsewhere
+        Toy(:,:)=Sdy(:,:)*Hmy(:,:)*Ro*G
+     end where
+
+     do J=2,Ny
+        do I=2,Nx
+
+           !     Variables 2d
+           Pente_stag(I,J)=(Sdx(I,J)+Sdx(I,J-1))**2 + &    ! Pente
+                (Sdy(I,J)+Sdy(I-1,J))**2
+
+           Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
+           Pente_stag(I,J)=Pente_stag(I,J)**0.5
+
+
+           Tob_stag(I,J)=(Tox(I,J)+Tox(I,J-1))**2 + &          ! Pente
+                (Toy(I,J)+Toy(I-1,J))**2
+
+           Tob_stag(I,J)=Tob_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
+           Tob_stag(I,J)=Tob_stag(I,J)**0.5
+
+           Uslid_stag(I,J)=(Ux(I,J,Nz)+Ux(I,J-1,Nz))**2 + &        ! Vitesse De Bilan
+                (Uy(I,J,Nz)+Uy(I-1,J,Nz))**2
+           Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
+           Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
+
+           Hmxy(I,J)=((H(I,J)+H(I-1,J-1))+(H(I,J-1)+H(I-1,J)))*0.25
+
+           !     En Vertical : Calcul De La Deformation Par Differentiation Des Vitesses
+
+           do K=1,Nz
+
+              Vit_stag3(I,J,K)=(Ux(I,J,K)+Ux(I,J-1,K))**2 + &        ! Magnitude Vitesse 
+                   (Uy(I,J,K)+Uy(I-1,J,K))**2            ! A Tous Niveaux K
+
+              Vit_stag3(I,J,K)=Vit_stag3(I,J,K)*0.25
+              Vit_stag3(I,J,K)=Vit_stag3(I,J,K)**0.5
+           end do
+
+           Qdef3(I,J,1)=0.
+           do K=2,Nz-1
+              Qdef3(I,J,K)=Vit_stag3(I,J,K-1)-Vit_stag3(I,J,K+1)            ! Difference Des Vitesses
+              Qdef3(I,J,K)=Ro*G*Pente_stag(I,J)*Ee(K)*Qdef3(I,J,K)/2./Dee   ! Gamma Tau
+           end do
+
+
+           Qdef3(I,J,Nz)=Vit_stag3(I,J,Nz-1)-Vit_stag3(I,J,Nz)        ! Pour Le Fond Differentiation
+           Qdef3(I,J,Nz)=Ro*G*Pente_stag(I,J)*Qdef3(I,J,K)/Dee            ! Sur Une Seule Maille
+        end do
+     end do
+
+     Qslid(:,:)=Tob_stag*Uslid_stag(:,:)
+
+     ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
+
+     do J=2,Ny-1
+        do I=2,Nx-1
+           Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
+           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
+
+           ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
+           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T(I,J,Nz)-Tpmp(I,J,Nz))*Ecart_phid)
+
+           ! Chaleur Due A La Defromation
+           do K=1,Nz         
+              Chaldef_maj(I,J,K)=(Qdef3(I,J,K)+Qdef3(I+1,J+1,K))+(Qdef3(I,J+1,K)+Qdef3(I+1,J,K))
+              Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)*0.25
+              Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)/Cp(I,J,K)
+           end do
+
+        end do
+     end do
+
+     !
+     !Write(126,*)'Time=',Time
+     !J=41
+     !Do I=20,30
+     !   Write(126,'(I3,6(E14.4,1x))') I,Chalgliss_maj(I,J),Chaldef_maj(I,J,Nz_m),Qdef2(I,J),Qdef2(I+1,J+1) &
+     !       , Qdef2(I,J+1),Qdef2(I+1,J)
+     !End Do
+
+     ! Flux Total A Rajouter Ã La Base
+     Phid(:,:)=Chalgliss_maj(:,:)+Chaldef_maj(:,:,Nz)*Fracq*H(:,:)*Cp(:,:,Nz)
+
+  end select
+end subroutine Qprod_ice
Index: /trunk/SOURCES/Temperature-routines/advec_icetemp.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/advec_icetemp.f90	(revision 24)
+++ /trunk/SOURCES/Temperature-routines/advec_icetemp.f90	(revision 24)
@@ -0,0 +1,104 @@
+!> \file advec_icetemp.f90
+!! Subroutines for temperature advection in the heat transfer equation
+!<
+
+!> \namespace advec_icetemp
+!! Subroutines to calculate the temperature advection 
+!<
+
+!> SUBROUTINE: Advec_icetemp
+!! Implementation of Routines to calculate advection in the heat transfer equation
+!! for a column
+!! Used modules:
+!!    - use Icetemp_declar
+!>
+
+Subroutine  Advec_icetemp(Ii,Jj,Prodq_m,Cro_m,Advecx_m,Advecy_m,Advec_m,Ct_m)
+
+  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.
+
+ 
+  Do K=2,Nz-1
+
+
+     ! On Remplit Deux Tableaux Pour eviter Des Problemes D'Asymetrie
+     !         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
+            
+          + 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
+
+     Advecx_m(K) =  Advecx_m(K) * Dx11                                        !Dx11=1/Dx
+
+
+     !         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
+          
+          + 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
+
+     Advecy_m(K) =  Advecy_m(K) * Dx11
+
+     Advec_m(K)=Advecx_m(K)+Advecy_m(K)  ! Advection Totale
+
+  End Do
+
+  !     Fin Boucle K, Ecriture Matrice Tridiag Dans La Glace
+
+
+  !      -----------------------------Cas General (H>10m)
+  thick_ice:  If (H(Ii,Jj).Gt.10.) Then
+
+     !        Variables De Calcul Dans La Glace
+     ! dou = dtt/dz^2
+     Dou=Dtt/Dee/Dee/ H(Ii,Jj) / H(Ii,Jj)  
+
+     ! Dah = dtt/dz
+     Dah=Dtt /Dee /H(Ii,Jj)
+
+     ! thermal conductivity at mid point just below the surface                   
+     Ct_bas=2.*(Ct_m(1)*Ct_m(2))/(Ct_m(1)+Ct_m(2))           
+
+     ! surface temperature : cannot go above 0 celsius
+     T(Ii,Jj,1)=Min(0.,Ts(Ii,Jj))
+
+     ! surface boundary condition
+     Aa(1)=0.
+     Bb(1)=1.
+     Cc(1)=0.
+     Rr(1)=T(Ii,Jj,1)
+
+     Do K=2,Nz-1
+        Dzz=Dou/Cro_m(K)
+
+        ! Conductivite Au Milieu Des Mailles
+        Ct_haut=Ct_bas
+        Ct_bas=2.*(Ct_m(K)*Ct_m(K+1))/(Ct_m(K)+Ct_m(K+1))
+
+        ! 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
+        Bb(K) =  1.+ Dzz * (Ct_haut+Ct_bas)                                ! diag
+        
+        Rr(K) =  T(Ii,Jj,K) + Prodq_m(K) * Dtt   ! vector
+        Rr(K) =  Rr(K)- Dttdx * (Advec_m(K))
+
+     End DO
+  End If thick_ice
+
+End Subroutine Advec_icetemp
+
Index: unk/SOURCES/Temperature-routines/icetemp_declar-0.5_mod.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/icetemp_declar-0.5_mod.f90	(revision 23)
+++ 	(revision )
@@ -1,90 +1,0 @@
-!> \file icetemp_declar-0.5_mod.f90
-!! Declar local variables for the temperature calculation in ice 
-!<
-
-
-!> \namespace icetemp_declar 
-!! module to declar local variables for the temperature calculation in ice 
-!! @note used modules:
-!! @note - Module_choix
-!! @note - Param_phy_mod
-!<
-module icetemp_declar     
-
-  Use Module_choix
-  Use Param_phy_mod
-
-  Implicit None 
-
-  !> \Var
-  integer ::  Nx_m,Ny_m,Nz_m,Nzm_m,NzNzm_m 
-  Integer :: Iec,Jec,Ngr
-  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,Parameter :: Acof1=-0.0575  !< Pour La Temperature De L'Eau De Mer
-  Real,Parameter :: Bcof1=0.0901   !< Pour La Temperature De L'Eau De Mer
-  Real,Parameter :: Ccof1=7.61e-4  !< Pour La Temperature De L'Eau De Mer
-  Real,Parameter :: S0mer=34.75    !< Pour La Temperature De L'Eau De Mer
-
-  ! a la Base De L'Ice Shelf
-  Integer :: Ifail1             !< Permet De Detecter Les Erreurs 
-  Real :: Chalk_1               !< Utilise Pour Le Calcul De Chalk : Glace PosÃ©e
-  Real :: Chalk_2               !< Utilise Pour Le Calcul De Chalk : Ice Streams Et Ice Shelves
-  Real :: Chalbed_1             !< Utilise Pour Le Calcul De Chalbed
-  Real :: Coefadv               !< Pour Limiter Le Flux De Chaleur Horiz.
-  Real :: Ecart_phid            !< Pour Limiter L'Impact De La Chaleur De Glissement Si Base Froide
-  ! Ancienne Valeur 0.5.
-
-  ! Tableaux 1d
-  ! _______________
-
-  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,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  
-
-  ! Tableaux De Travail 2d 
-  ! ___________________________
-  Real,Allocatable,Dimension(:,:):: Tbmer    !< Temperature De La Mer A La Base De L'Ice Shelf  
-  Real,Allocatable,Dimension(:,:) :: Alpha   !< Pente Locale Sur Noeud Majeur                    
-  Real,Allocatable,Dimension(:,:) :: Ubar    !< Vitesse Sur Noeud Majeur                       
-  Real,Allocatable,Dimension(:,:) :: Chalglissx,Chalglissy !< Chaleur De Glissement             
-  Integer,Allocatable,Dimension(:,:) :: 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
-  ! 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(:,:,:) :: Cp            !< Specific Heat Capacity (J/(M-3)/K)=Ro Cp  
-  Real,Allocatable,Dimension(:,:,:) :: Ct            !< Thermal Conductivity (J/M/K/A)            
-  Real,Allocatable,Dimension(:,:,:) :: Chaldef_maj   !< Chaleur De Deformation                   
-  Real,Allocatable,Dimension(:,:,:) :: Advecx,Advecy,Advec
-
-  !Variables De Travail
-  Real ::  Fracq                   !< For Temp Calculation = (1.-(1.-De/2.)**5)/5.
-  Integer,Parameter :: Ncond=0     !< Switch : Ncond=1 Thermal Conductivity In Mantle
-  Real :: Dee                      !< Vertical Step In Ice And Mantle
-  Real,Parameter :: Dzm=600        !< Grid Step In Mantle
-  Real :: Dttdx                    !< For Temp Calculation = Dtt/Dx
-  Real :: Dx11                     !< 1/Dx
-
-  !Prop Thermique 
-  Real,Parameter ::  Cm=1.04e8              !< Thermal Conductivity Of Mantle J/M/K/A
-  Real,Parameter ::  Cpm=1000.0             !< Specific Heat Capacity Of Mantle J/Kg/K
-  Real:: Ctm                                !< For Temp Calculation
-
-
-end module icetemp_declar
-
Index: /trunk/SOURCES/Temperature-routines/icetemp_declar_mod.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/icetemp_declar_mod.f90	(revision 24)
+++ /trunk/SOURCES/Temperature-routines/icetemp_declar_mod.f90	(revision 24)
@@ -0,0 +1,91 @@
+!> \file icetemp_declar-0.5_mod.f90
+!! Declar local variables for the temperature calculation in ice 
+!<
+
+
+!> \namespace icetemp_declar 
+!! module to declar local variables for the temperature calculation in ice 
+!! @note used modules:
+!! @note - Module_choix
+!! @note - Param_phy_mod
+!<
+module icetemp_declar     
+
+  Use Module_choix
+  Use Param_phy_mod
+
+  Implicit None 
+
+  !> \Var
+!  integer ::  Nx_m,Ny_m,Nz_m,Nzm_m,NzNzm_m 
+  Integer :: Iec,Jec,Ngr
+  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,Parameter :: Acof1=-0.0575  !< Pour La Temperature De L'Eau De Mer
+  Real,Parameter :: Bcof1=0.0901   !< Pour La Temperature De L'Eau De Mer
+  Real,Parameter :: Ccof1=7.61e-4  !< Pour La Temperature De L'Eau De Mer
+  Real,Parameter :: S0mer=34.75    !< Pour La Temperature De L'Eau De Mer
+
+  ! a la Base De L'Ice Shelf
+  Integer :: Ifail1             !< Permet De Detecter Les Erreurs 
+  Real :: Chalk_1               !< Utilise Pour Le Calcul De Chalk : Glace PosÃ©e
+  Real :: Chalk_2               !< Utilise Pour Le Calcul De Chalk : Ice Streams Et Ice Shelves
+  Real :: Chalbed_1             !< Utilise Pour Le Calcul De Chalbed
+  Real :: Coefadv               !< Pour Limiter Le Flux De Chaleur Horiz.
+  Real :: Ecart_phid            !< Pour Limiter L'Impact De La Chaleur De Glissement Si Base Froide
+  ! Ancienne Valeur 0.5.
+
+  ! Tableaux 1d
+  ! _______________
+
+  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,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  
+
+  ! Tableaux De Travail 2d 
+  ! ___________________________
+  Real,allocatable,Dimension(:,:):: Tbmer    !< Temperature De La Mer A La Base De L'Ice Shelf  
+  Real,allocatable,Dimension(:,:) :: Alpha   !< Pente Locale Sur Noeud Majeur                    
+  Real,allocatable,Dimension(:,:) :: Ubar    !< Vitesse Sur Noeud Majeur                       
+  Real,allocatable,Dimension(:,:) :: Chalglissx,Chalglissy !< Chaleur De Glissement             
+  Integer,allocatable,Dimension(:,:) :: 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
+  ! 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,Nz) :: Cp            !< Specific Heat Capacity (J/(M-3)/K)=Ro Cp  
+  Real,Dimension(Nx,Ny,Nz) :: Ct            !< Thermal Conductivity (J/M/K/A)            
+
+
+  !Variables De Travail
+  Real ::  Fracq                   !< For Temp Calculation = (1.-(1.-De/2.)**5)/5.
+  Integer,Parameter :: Ncond=0     !< Switch : Ncond=1 Thermal Conductivity In Mantle
+  Real :: Dee                      !< Vertical Step In Ice And Mantle
+  Real,Parameter :: Dzm=600        !< Grid Step In Mantle
+  Real :: Dttdx                    !< For Temp Calculation = Dtt/Dx
+  Real :: Dx11                     !< 1/Dx
+
+  !Prop Thermique 
+  Real,Parameter ::  Cm=1.04e8              !< Thermal Conductivity Of Mantle J/M/K/A
+  Real,Parameter ::  Cpm=1000.0             !< Specific Heat Capacity Of Mantle J/Kg/K
+  Real:: Ctm                                !< For Temp Calculation
+
+
+end module icetemp_declar
+
Index: /trunk/SOURCES/Temperature-routines/icetemp_mod.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/icetemp_mod.f90	(revision 24)
+++ /trunk/SOURCES/Temperature-routines/icetemp_mod.f90	(revision 24)
@@ -0,0 +1,303 @@
+!> \file icetemp_mod.f90
+!! Interfaces of functions and subroutines for temperature calculation 
+!<
+
+!> \namespace icetempmod
+!! Interfaces of functions and subroutines for the temperature calculation 
+!<
+
+module icetempmod
+
+contains
+
+  !> SUBROUTINE: init_icetemp 
+  !! Initialise the temperature calculation
+  !! Used modules:
+  !! - Icetemp_declar
+  !>
+
+  Subroutine init_icetemp(Num_rep_42_m)
+    Use Icetemp_declar, only : Nfracq,Iq,Ecart_phid
+
+    Implicit None
+
+    !< Arguments
+    Integer,                       Intent(In)::Num_rep_42_m !< Id Of Reponse
+    ! Character(Len=80) :: Filin
+    !________________________________________________________________
+    ! La Temperature Du Point De Congelation De L'Eau De Mer
+    ! Est Tiree De Jenkins (1991), Tb Est En Degres C
+    ! _______________________________________________________________
+
+    !  Nouvelle Formulation De Fracq
+    !  Nfracq=15  ! Ancienne Formulation
+
+    Nfracq=1
+
+    Write(Num_rep_42_m,*) 'Calcul Des Temperatures'
+    Write(Num_rep_42_m,*) 'Nfracq=',Nfracq
+
+    Iq=1
+    If (Iq.Eq.1) Write(Num_rep_42_m,*) 'Chaleur Demi Maille, Iq=',Iq
+    If (Iq.Eq.2) Write(Num_rep_42_m,*) 'Chaleur Centree, Iq=',Iq
+    If (Iq.Eq.3) Write(Num_rep_42_m,*) 'Chaleur Pente, Iq=',Iq
+    If (Iq.Eq.4) Write(Num_rep_42_m,*) 'Chaleur U_taub, Iq=',Iq
+    If (Iq.Eq.5) Write(Num_rep_42_m,*) 'Chaleur U_taub_stag, Iq=',Iq
+    If (Iq.Eq.6) Write(Num_rep_42_m,*) 'Chaleur Q_sia_stag,Iq=',Iq
+    If (Iq.Eq.7) Write(Num_rep_42_m,*) 'Chaleur Q_all_stag,Iq=',Iq
+    Write(Num_rep_42_m,*)
+
+    Ecart_phid=0.5
+    Write(Num_rep_42_m,*) ' Si Base Froide, La Chaleur De Glissement Est En Exp(Delta T * Ecart_phid)'
+    Write(Num_rep_42_m,*) 'Ecart_phid=',Ecart_phid
+    Write(Num_rep_42_m,*) '---------------------------------------------------------------------------'
+
+  End Subroutine init_icetemp
+  
+
+  !-----------------------------------------------------------------------------------------------------------------------
+
+  !> SUBROUTINE: icetemp
+  !! Calculate the temperature
+  !!
+  !! Used modules:
+  !!  -  use Icetemp_declar 
+  !>
+
+  subroutine icetemp
+
+    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),&   
+            Alpha(Nx,Ny),&                   
+            Ubar(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
+
+    ! init
+    Aa=0.    
+    Bb=0.              
+    Cc=1.             
+    Rr=0.              
+    Hh=0.              
+    Tdot=0  
+    !Instruction 
+    Dee=1./(Nz-1) 
+    Fracq=(1.-(1.-Dee/2.)**Nfracq)/Nfracq
+    Fracq=(1.-(1.-Dee/2.)**5.)/5.
+    Ee(1)=0.
+    Ee(Nz)=1.
+    Do K=1,Nz
+       If ((K.Ne.1).And.(K.Ne.Nz)) Ee(K)=(K-1.)/(Nz-1.)
+    End Do
+
+    !   Variables Dependant Du Pas De Temps Dtt
+    Ctm=Dtt*Cm/Rom/Cpm/Dzm/Dzm
+    Dttdx=Dtt/Dx
+    Dx11=1./Dx
+
+    If (Itracebug.Eq.1) Write(num_tracebug,*)' Entree Dans Routine Ice_temp Time=',Time
+
+    !________________________________________________________________
+    ! La Temperature Du Point De Congelation De L'Eau De Mer
+    ! Est Tiree De Jenkins (1991), Tb Est En Degres C
+    ! _______________________________________________________________
+    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Routine Temp_mer'
+    Call Temp_mer
+
+    ! Appel Aux Proprietes Thermiques De La Glace : Conductivite, Capacite Claorifique -Tpmp
+
+    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Routine Thermal_property'
+    Call Thermal_prop_icetemp
+
+    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Routine Qprod'
+    Call Qprod_ice(Iq)
+
+
+    ! Dans Le Socle : Les elements de La Matrice Tridiag Sont Invariants Dans L'Espace
+    Do K=Nz+1,nz+nzm-1
+       Aa(K)=-Ctm
+       Bb(K)=1.+2.*Ctm
+       Cc(K)=-Ctm
+    End Do
+
+    ! Conditions Aux Limites Pour Le Socle
+    Aa(nz+nzm)=-1.
+    Bb(nz+nzm)=1.
+    Cc(nz+nzm)=0.
+
+    ! Calcul D'Advection : On Traite Les Tableaux En Une Seule Ligne De Commandes
+    ! les tableaux iadvec valent 1 quand la vitesse consideree est upwind et 0 sinon
+
+    !         Advection Selon X
+    !         -----------------
+
+    Iadvec_w(:,:)=Nint(Max(Sign(1.0,Uxbar(:,:)),0.0))
+    Iadvec_e(:,:)=Nint(Abs(Min(Sign(1.0,Uxbar(:,:)),0.0)))
+
+    !         Advection Selon Y
+    !         -----------------
+
+    Iadvec_s(:,:)=Nint(Max(Sign(1.0,Uybar(:,:)),0.0)) 
+    Iadvec_n(:,:)=Nint(Abs(Min(Sign(1.0,Uybar(:,:)),0.0)))
+
+    ! Boucle De Resolution Des Temperatures Colonne Par Colonne
+
+    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  boucle Advec_icetemp'
+
+    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,:))
+
+          ! variables de calcul dans la glace
+          Deh22=2.*Dee*H(I,J)**2
+          Duu=Dtt/2./Dee
+          Rr(nz+nzm)=-Dzm*Ghf(I,J)/Cm  ! depend de I,J Pour Un Flux Geothermique Variable
+
+!   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))
+
+          ! Test Permettant D'Ecrire Les Variables Quand Il Y A Une Erreur
+          ! Ifail1 = 1 : Erreur Dans Tridiag
+          If (Ifail1.Eq.1) Then
+             Write(6,*)'Erreur Dans Tridiag'
+             Write(6,*)'I = ',I,' J = ',J
+             Write(6,*) 'A, B, C, R, U'
+             Do K=1,nz+nzm
+                Write(6,*)'K = ',K
+                !Write(6,*) Aa(K),Bb(K),Cc(K),Rr(K),Hh(K)
+             End Do
+             Stop
+          End If
+       End Do
+    End Do
+
+    Do J=1,Ny
+       Do I=1,Nx
+          H1(I,J)=H(I,J)
+          B1(I,J)=B(I,J)
+          Tpmp(I,J,1)=0.
+       End Do
+    End Do
+
+    ! Attribution Finale De La Temperature
+    ! Limitation a 3 deg De Variation Par Dtt Et Pas Plus Froid Que -70
+
+    Do K=1,nz+nzm
+       Do J=1,Ny
+          Do I=1,Nx
+             If (H(I,J).Gt.10.) Then
+
+                Tdot(K)=T3d_new(I,J,K)-T(I,J,K)
+
+                If ((K.Gt.1).And.(Tdot(K).Lt.-3.)) Tdot(K)=-3.
+                If ((K.Gt.1).And.(Tdot(K).Gt.3.)) Tdot(K)=3.
+
+                T(I,J,K)=T(I,J,K)+Tdot(K)
+
+                If (T(I,J,K).Lt.-70.) Then
+                   T(I,J,K)=-70.
+                Endif
+             Endif
+
+          End Do
+       End Do
+    End Do
+
+
+    Do K=2,Nz
+       Do J=1,Ny
+          Do I=1,Nx
+             If (T(I,J,K).Gt.Tpmp(I,J,K)) Then
+                T(I,J,K)=Tpmp(I,J,K)
+                Ibase(I,J)=2
+             Endif
+          End Do
+       End Do
+    End Do
+
+    !     Temperature Sur Les Bords :
+    !     -------------------------
+  
+    Do K=1,Nz
+       Do J=1,Ny
+          !        T(Nx,J,K)= T(Nx-1,J,K)
+          !        T(1,J,K)= T(2,J,K) 
+          T(Nx,J,K)= Ts(Nx,J)
+          T(1,J,K)= Ts(1,J)
+       End Do
+       !
+       Do I=1,Nx
+          !        T(I,1,K)=T(I,2,K) 
+          !        T(I,Ny,K)=T(I,Ny-1,K) 
+          T(I,1,K)=Ts(I,1)
+          T(I,Ny,K)=Ts(I,Ny)
+       End Do
+    End Do
+   
+
+    Do K=1,Nz
+       Do J=1,Ny
+          Do I=1,Nx
+             If (H(I,J).le.1.) then 
+                T(I,J,K)=Ts(I,J)
+                if ((k.ge.nz).and.(flot(i,j))) T(I,J,K) =Tbmer(I,J)
+             end If
+                
+          End Do
+       End Do
+    End Do
+
+    Deallocate (Aa,Bb,Cc,Rr,Hh,Tdot,Abis,Bbis,Cbis,Rbis,Hbis,Ee,Tbmer,            & 
+                Alpha,Ubar,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)        
+
+    If (Itracebug.Eq.1) Write(num_tracebug,*)' fin routine icetemp'
+    return
+  End Subroutine icetemp
+
+  !-----------------------------------------------------------------------------------------------------------------------
+
+end module icetempmod
Index: unk/SOURCES/Temperature-routines/interface_Qprod_icetemp.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/interface_Qprod_icetemp.f90	(revision 23)
+++ 	(revision )
@@ -1,947 +1,0 @@
-!> \file interface_Qprod_icetemp.f90
-!!Interfaces of functions and subroutines for heat production 
-!<
-
-!> \namespace interface_Qprod_icetemp
-!!Interfaces of functions and subroutines for heat production 
-!<
-
-module interface_Qprod_icetemp
-  !> \interface Qprod_ice
-  !!
-  !! Interface of functions and subroutines for heat production 
-  !!\author ...
-  !!\date 2010
-  !<
-  Interface  Qprod_ice
-     module procedure Qprod 
-  End Interface
-
-
-contains
-
-
-  !> SUBROUTINE: Qprod 
-  !! Subroutine for heat pruduction
-  !! Used modules:
-  !!    - use Icetemp_declar
-  !!    - use geom_type
-  !!    - use temperature_type
-  !!    - use Ice_flow_type
-  !!    - use Mask_flgz_type
-  !!    - use Deformation_type
-  !!    - use Autre_pr_temp_type
-  !>
-
-
-
-  subroutine Qprod (Iq1,geom_m,T_m,Ice_flow_m,mask_flot_m,deform_m,therm_var_m,debug_3D_m)
-
-    use Icetemp_declar
-    use geom_type
-    use temperature_type
-    use Ice_flow_type
-    use Mask_flgz_type
-    use Deformation_type
-    use Autre_pr_temp_type
-
-    implicit none
-
-    !<Arguments
-
-    integer,              intent(In)    :: Iq1              !< Choice Of Method For Heat Production
-    type (Geom),          intent(in)    :: geom_m           !< variable geometrie 
-    type (Temperature),   intent(in)    :: T_m              !< variable temperature 
-    type (Ice_flow),      intent(in)    :: Ice_flow_m       !< ice_flow variables  
-    type (Mask_flgz),    intent(in)    :: mask_flot_m      !< mask variable 
-    type (Deformation),   intent(in)    :: deform_m         !< deformation variable
-    type (Autre_pr_temp), intent(inout) :: therm_var_m      !< thermal variable
-    real,dimension(:,:,:),intent(inout) :: debug_3D_m       !< Debug array
-
-    !< Local variables needed by the various case of Iq1
-
-    real, dimension(Nx_m,Ny_m)      :: Pente2_maj    !< Pente Au Carre Sur Les Noeuds Majeurs
-    real, dimension(Nx_m,Ny_m)      :: Pente_maj     !< Pente Au Carre Sur Les Noeuds Majeurs
-    real, dimension(Nx_m,Ny_m)      :: Vit_maj       !< Ubar Moyennee Sur Les Noeuds Majeurs
-    real, dimension(Nx_m,Ny_m)      :: Uslid_maj     !< Uslid Moyennee Sur Les Noeuds Majeurs
-    real, dimension(Nx_m,Ny_m)      :: Vit_stag2     !< Ubar Moyennee Sur Les Noeuds Stag
-    real, dimension(Nx_m,Ny_m,Nz_m) :: Vit_stag3     !< Ubar Moyennee Sur Les Noeuds Stag
-    real, dimension(Nx_m,Ny_m)      :: Hmxy          !< Epaisseur Moyenne Sur Les Noeuds Stag
-    real, dimension(Nx_m,Ny_m)      :: Qslid         !< Chaleur Produite Par Glissement Sur Les Noeuds Stag
-    real, dimension(Nx_m,Ny_m)      :: Qdef2         !< Chaleur Produite Par Deformation Sur Les Noeuds Stag
-    real, dimension(Nx_m,Ny_m,Nz_m) :: Qdef3         !< Chaleur Produite Par Deformation Sur Les Noeuds Stag
-    real, dimension(Nx_m,Ny_m)      :: Pente_stag    !< Pente Au Carre Sur Les Noeuds Stag
-    real, dimension(Nx_m,Ny_m)      :: Uslid_stag    !< Uslid Moyennee Sur Les Noeuds Stag
-    real, dimension(Nx_m,Ny_m)      :: Tob_stag      !< Contrainte Basale
-    real, dimension(Nx_m,Ny_m)      :: Tox           !< Contraintes Sur Maille Mx
-    real, dimension(Nx_m,Ny_m)      :: Toy           !< Contraintes Sur Maille Mx
-
-
-
-    select case (Iq1)
-
-    case (1)                ! iq1 = 1      Q_prod_demi, Ancienne Methode Sur Les Demi-Mailles
-       !------------------------------------------------------------------------------------------
-
-       ! Calcul Avec L'Ancienne Methode Sur Les Demi-Mailles
-       ! Nom Du Sub Routine Avant Et Var Utilisee  
-       ! Q_prod_demi(T_m%T,T_m%Tpmp,Deform_m%Tobmx,Deform_m%Tobmy,Geom_m%H,Geom_m%Hmx,Geom_m%Hmy,mask_flot_m%Flot,mask_flot_m%Flotmx,mask_flot_m%Flotmy,mask_flot_m%Flgzmx ,mask_flot_m%Flgzmy,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Ddx,Ice_flow_m%Ddy,&
-       !      Ice_flow_m%Ddbx,Ice_flow_m%Ddby,Deform_m%Epsxx,Deform_m%Epsyy,Deform_m%Epsxy,mask_flot_m%Gzmx,mask_flot_m%Gzmy,Deform_m%Btt,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid,Deform_m%Glen,Deform_m%Visc)
-
-       do K=2,Nz_m
-          do J=2,Ny_m-1
-             do I=2,Nx_m-1
-
-                ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
-                ! Pour Les Ice-Streams Et Ice-Shelves
-                ! Les Divers Eps Sont Calcules Dans Strain-Rate Pour L'Ensemble De La Grille
-
-                if ((mask_flot_m%Flottant(I,J).or.mask_flot_m%Flgz_mx(I,J).or.mask_flot_m%Flgz_mx(I+1,J)).or.  &
-                     (mask_flot_m%Flgz_my(I,J).or.mask_flot_m%Flgz_my(I,J+1))) then
-
-                   Chal2_x(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_xx(I,J)**2
-                   Chal2_y(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_yy(I,J)**2
-                   Chal2_z(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(-Deform_m%Veloc_Deform_xx(I,J)-Deform_m%Veloc_Deform_yy(I,J))**2
-
-                   !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
-                   Chal2_xy(I,J,K)=(Deform_m%Veloc_Deform_xy(I,J)+Deform_m%Veloc_Deform_xy(I+1,J)+Deform_m%Veloc_Deform_xy(I+1,J+1)+Deform_m%Veloc_Deform_xy(I,J+1))
-                   Chal2_xy(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
-
-                else   ! Glace Posee
-                   Chal2_x(I,J,K)=0.
-                   Chal2_y(I,J,K)=0.
-                   Chal2_z(I,J,K)=0.
-                   Chal2_xy(I,J,K)=0.
-                endif
-             end do
-          end do
-       end do
-
-       !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
-       do L=1,size(Deform_m%Btt_Deform,4) !N1poly,N2poly
-          do J=2,Ny_m
-             do I=2,Nx_m
-                !          Ffx A 3 Dimensions !
-                Ffx(I,J,L)=Ice_flow_m%Difus_x(I,J,L)*Geom_m%Slop_x(I,J)*Geom_m%Slop_x(I,J)
-                Ffy(I,J,L)=Ice_flow_m%Difus_y(I,J,L)*Geom_m%Slop_y(I,J)*Geom_m%Slop_y(I,J)
-             end do
-          end do
-       end do
-
-       do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
-          do K=2,Nz_m
-             do J=2,Ny_m
-                do I=2,Nx_m
-                   if ((.not.mask_flot_m%Flot_mx(I,J)).and.(.not.mask_flot_m%Gz_mx(I,J))) then
-                      Chalx(I,J,K,L)=(Deform_m%Btt_Deform(I-1,J,K,L)+Deform_m%Btt_Deform(I,J,K,L))*Ffx(I,J,L) !&
-                      !                      *Ro*G*Ee(K)**(Glen(L)+1)/Cp(I,J,K)
-
-                   else if (mask_flot_m%Gz_mx(I,J)) then ! Ice Streams
-                      Chalx(I,J,K,L)=0.
-                   else                     ! Ice Shelves  
-                      Chalx(I,J,K,L)=0.
-                   endif
-                   if ((.not.mask_flot_m%Flot_my(I,J)).and.(.not.mask_flot_m%Gz_my(I,J))) then
-                      Chaly(I,J,K,L)=(Deform_m%Btt_Deform(I,J-1,K,L)+Deform_m%Btt_Deform(I,J,K,L))*Ffy(I,J,L) !&
-                      !                      *Ro*G*Ee(K)**(Glen(L)+1)/Cp(I,J,K)
-                   else if (mask_flot_m%Gz_my(I,J)) then  ! Ice Streams
-                      Chaly(I,J,K,L)=0.
-
-                   else                      ! Ice Shelves
-                      Chaly(I,J,K,L)=0.
-                   endif
-
-                end do
-             end do
-          end do
-       end do
-
-       !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
-       !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
-       !         Les Productions X Et Y
-       !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy) 
-
-       do K=2,Nz_m
-          do J=1,Ny_m-1
-             do I=1,Nx_m-1
-
-                ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
-                Chaldef_maj(I,J,K)= 0.
-
-                ! Chalk_2 Pour Ice Shelves Et Ice Streams
-                Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
-                     Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
-
-                ! Chalk_1 Pour La Partie PosÃ©e
-                do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
-
-                   ! On Somme La Chaleur Due Aux Diverses Lois De DÃ©formation Et Celles En X Et Y
-                   Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))+ &
-                        (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
-
-                   Chalk_1=Ro*G*Chalk_1/4.*(Ee(K)**(Deform_m%Glen_exp(L)+1))/Cp(I,J,K)
-                   Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
-                enddo
-
-                ! Pour Shelves Et Streams,  On Ajoute Chalk_2
-                Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
-             end do
-          end do
-       end do
-       ! Chaleur Produite A La Base Par Le Glissement
-       do J=2,Ny_m
-          do I=2,Nx_m
-
-             if (mask_flot_m%Gz_mx(I,J)) then
-                Chalglissx(I,J)= abs(Ice_flow_m%U_Column_x(I,J)*Deform_m%Cisail_Basal_mx(I,J))
-             else
-                Chalglissx(I,J)=Ice_flow_m%Difus_bx(I,J)*Geom_m%Slop_x(I,J)**2*Ro*G*Geom_m%Thick_mx(I,J)
-             endif
-
-             if (mask_flot_m%Gz_my(I,J)) then
-                Chalglissy(I,J)= abs(Ice_flow_m%U_Column_y(I,J)*Deform_m%Cisail_Basal_my(I,J))
-             else
-                Chalglissy(I,J)=Ice_flow_m%Difus_by(I,J)*Geom_m%Slop_y(I,J)**2*Ro*G*Geom_m%Thick_my(I,J)
-             endif
-
-          end do
-       end do
-
-       !  Boundary Condition Ice-Rock Interface 
-       K=Nz_m
-
-       do J=2,Ny_m-1
-          do I=2,Nx_m-1
-
-             ! L'Ancien Chalbed Correspond A Chaldef_maj(I,J,Nz_m)  : Pas La Peine De Recalculer
-             !----------------------------------------------------
-             !      Chalbed=0.
-             !      Do L=1,size(Deform_m%Btt,4)!N1poly,N2poly
-             !         Chalbed_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L)) &
-             !              + (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
-             !         Chalbed_1=Ro*G*Chalbed_1/4.*H(I,J)
-             !         Chalbed= Chalbed_1 +Chalbed
-             !      Enddo
-
-             !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
-             !        Dans La Derniere 1/2 Maille Et Par Le Glissement
-             !        Attention Phid Est >0 Et Ghf Est <0
-             if (.not.mask_flot_m%Flottant(I,J)) then
-
-                !           Phid Avec Fonte Sous Les Streams
-                !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
-                !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
-                !                        Chalbed*Fracq
-
-
-                ! Moyenne Phid Sur 4 Points : Attention 0.5 Pour Moyenne Gauche - Droite 
-                ! Mais Les Chaleurs X,Y S'Ajoutent
-
-                !              Phid(I,J)=0.5*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
-                !                        Chalglissy(I,J)+Chalglissy(I,J+1)) 
-
-                ! Moyenne Phid Sur 12 Points. On Moyenne Chaque Chaleur Et On Somme X Et Y
-
-                therm_var_m%Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
-                     Chalglissy(I,J)+Chalglissy(I,J+1)) 
-                therm_var_m%Phid(I,J)=therm_var_m%Phid(I,J)+0.125*( &
-                     ((Chalglissx(I,J-1)+Chalglissx(I+1,J-1))+ &
-                     (Chalglissx(I,J+1)+Chalglissx(I+1,J+1))) + &
-                     ((Chalglissy(I-1,J)+Chalglissy(I-1,J+1)) +  &
-                     (Chalglissy(I+1,J)+Chalglissy(I+1,J+1))))
-
-                ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
-                therm_var_m%Phid(I,J)=therm_var_m%Phid(I,J)*exp((T_m%Temperature(I,J,Nz_m)-T_m%Temp_melting(I,J,Nz_m))*Ecart_phid)
-
-                ! Pour Sorties Heino
-                Chalgliss_maj(I,J)=therm_var_m%Phid(I,J)
-
-                therm_var_m%Phid(I,J)=therm_var_m%Phid(I,J)+ Chaldef_maj(I,J,Nz_m)*Fracq*Geom_m%Thickness(I,J)*Cp(I,J,Nz_m)
-
-             else
-                therm_var_m%Phid(I,J)=0.
-             endif
-          end do
-       end do
-
-
-    case (2)             ! Q_prod_centre : Calcul Avec La Somme Des Carres
-       !-------------------------------------------------------------------------------------
-       ! Calcul Avec La Somme Des Carres
-       ! Nom Du Sub Routine Avant Et Var Utilisee  
-       ! Q_prod_centre(T_m%T,T_m%Tpmp,Deform_m%Tobmx,Deform_m%Tobmy,Geom_m%H,Geom_m%Hmx,Geom_m%Hmy,mask_flot_m%Flot,mask_flot_m%Flotmx,mask_flot_m%Flotmy,mask_flot_m%Flgzmx,mask_flot_m%Flgzmy,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Ddx,Ice_flow_m%Ddy,&
-       !     Ice_flow_m%Ddbx,Ice_flow_m%Ddby,Deform_m%Epsxx,Deform_m%Epsyy,Deform_m%Epsxy,mask_flot_m%Gzmx,mask_flot_m%Gzmy,Deform_m%Btt,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid,Deform_m%Glen,Deform_m%Visc)
-
-       do K=2,Nz_m  
-          do J=2,Ny_m-1
-             do I=2,Nx_m-1
-
-
-                ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
-                ! Pour Les Ice-Streams Et Ice-Shelves
-                ! Les Divers Eps Sont CalculÃ©s Dans Strain-Rate Pour L'Ensemble De La Grille
-
-                if ((mask_flot_m%Flottant(I,J).or.mask_flot_m%Flgz_mx(I,J).or.mask_flot_m%Flgz_mx(I+1,J)).or.  &
-                     (mask_flot_m%Flgz_my(I,J).or.mask_flot_m%Flgz_my(I,J+1))) then
-
-                   Chal2_x(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_xx(I,J)**2
-                   Chal2_y(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_yy(I,J)**2
-                   Chal2_z(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(-Deform_m%Veloc_Deform_xx(I,J)-Deform_m%Veloc_Deform_yy(I,J))**2
-
-                   !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
-                   Chal2_xy(I,J,K)=(Deform_m%Veloc_Deform_xy(I,J)+Deform_m%Veloc_Deform_xy(I+1,J)+Deform_m%Veloc_Deform_xy(I+1,J+1)+Deform_m%Veloc_Deform_xy(I,J+1))
-                   Chal2_xy(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
-
-                else   ! Glace Posee
-                   Chal2_x(I,J,K)=0.
-                   Chal2_y(I,J,K)=0.
-                   Chal2_z(I,J,K)=0.
-                   Chal2_xy(I,J,K)=0.
-                endif
-             end do
-          end do
-       end do
-
-       !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
-       do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
-          do J=2,Ny_m
-             do I=2,Nx_m
-
-                !          Ffx A 3 Dimensions !
-                Ffx(I,J,L)=Ice_flow_m%Difus_x(I,J,L)*Geom_m%Slop_x(I,J)*Geom_m%Slop_x(I,J)
-                Ffy(I,J,L)=Ice_flow_m%Difus_y(I,J,L)*Geom_m%Slop_y(I,J)*Geom_m%Slop_y(I,J)
-             end do
-          end do
-       end do
-
-
-       do L=1,size(Deform_m%Btt_Deform,4) !N1poly,N2poly
-          do K=2,Nz_m
-             do J=2,Ny_m
-                do I=2,Nx_m
-                   if ((.not.mask_flot_m%Flot_mx(I,J)).and.(.not.mask_flot_m%Gz_mx(I,J))) then
-                      Chalx(I,J,K,L)=(Deform_m%Btt_Deform(I-1,J,K,L)+Deform_m%Btt_Deform(I,J,K,L))*Ffx(I,J,L) !&
-                      !                      *Ro*G*Ee(K)**(Glen(L)+1)/Cp(I,J,K)
-
-                   else if (mask_flot_m%Gz_mx(I,J)) then ! Ice Streams
-                      Chalx(I,J,K,L)=0.
-
-                   else                     ! Ice Shelves  
-                      Chalx(I,J,K,L)=0.
-
-                   endif
-
-                   if ((.not.mask_flot_m%Flot_my(I,J)).and.(.not.mask_flot_m%Gz_my(I,J))) then
-                      Chaly(I,J,K,L)=(Deform_m%Btt_Deform(I,J-1,K,L)+Deform_m%Btt_Deform(I,J,K,L))*Ffy(I,J,L) !&
-                      !                      *Ro*G*Ee(K)**(Glen(L)+1)/Cp(I,J,K)
-
-                   else if (mask_flot_m%Gz_my(I,J)) then  ! Ice Streams
-                      Chaly(I,J,K,L)=0.
-
-                   else                      ! Ice Shelves
-                      Chaly(I,J,K,L)=0.
-                   endif
-
-                end do
-             end do
-          end do
-       end do
-
-       !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
-       !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
-       !         Les Productions X Et Y
-       !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy) 
-
-       do K=2,Nz_m
-          do J=1,Ny_m-1
-             do I=1,Nx_m-1
-
-                ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
-                Chaldef_maj(I,J,K)= 0.
-
-                ! Chalk_2 Pour Ice Shelves Et Ice Streams
-                Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
-                     Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
-
-                ! Chalk_1 Pour La Partie PosÃ©e
-                do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
-
-                   ! On Somme La Chaleur Due Aux Diverses Lois De DÃ©formation Et Celles En X Et Y
-                   ! La Somme Se Fait Par (Cx2+Cy2)** 0.5. Le 4 Vient Des Moyennes Deform_m%Btt + Des Moyennes Gauche-Droite
-
-                   Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))**2+ &
-                        (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))**2
-                   Chalk_1=Chalk_1**0.5
-
-                   Chalk_1=Ro*G*Chalk_1/4.*(Ee(K)**(Deform_m%Glen_exp(L)+1))/Cp(I,J,K)
-                   Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
-                enddo
-
-                ! Pour Shelves Et Streams,  On Ajoute Chalk_2
-                Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
-             end do
-          end do
-       end do
-
-       ! Chaleur Produite A La Base Par Le Glissement
-
-       do J=2,Ny_m
-          do I=2,Nx_m
-
-             if (mask_flot_m%Gz_mx(I,J)) then
-                Chalglissx(I,J)= abs(Ice_flow_m%U_Column_x(I,J)*Deform_m%Cisail_Basal_mx(I,J))
-             else
-                Chalglissx(I,J)=Ice_flow_m%Difus_bx(I,J)*Geom_m%Slop_x(I,J)**2*Ro*G*Geom_m%Thick_mx(I,J)
-             endif
-
-             if (mask_flot_m%Gz_my(I,J)) then
-                Chalglissy(I,J)= abs(Ice_flow_m%U_Column_y(I,J)*Deform_m%Cisail_Basal_my(I,J))
-             else
-                Chalglissy(I,J)=Ice_flow_m%Difus_by(I,J)*Geom_m%Slop_y(I,J)**2*Ro*G*Geom_m%Thick_my(I,J)
-             endif
-
-          end do
-       end do
-
-       !  Boundary Condition Ice-Rock Interface 
-
-       K=Nz_m
-
-       do J=2,Ny_m-1
-          do I=2,Nx_m-1
-
-             !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
-             !        Dans La Derniere 1/2 Maille Et Par Le Glissement
-             !        Attention Phid Est >0 Et Ghf Est <0
-             if (.not.mask_flot_m%Flottant(I,J)) then
-
-                !           Phid Avec Fonte Sous Les Streams
-                !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
-                !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
-                !                        Chalbed*Fracq
-
-
-                ! Moyenne Phid Sur 4 Points Formulation (A2+B2)**0.5
-                !
-                Chalgliss_maj(I,J)=(Chalglissx(I,J)+Chalglissx(I+1,J))**2+  &
-                     (Chalglissy(I,J)+Chalglissy(I,J+1))**2
-
-                Chalgliss_maj(I,J)=Chalgliss_maj(I,J)**0.5
-
-                Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.5 ! Les Moyennes Droite Gauche
-
-                ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
-                Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T_m%Temperature(I,J,Nz_m)-T_m%Temp_melting(I,J,Nz_m))*Ecart_phid)
-
-
-                ! Flux Total A Rajouter Ã La Base
-                therm_var_m%Phid(I,J)=Chalgliss_maj(I,J)+Chaldef_maj(I,J,Nz_m)*Fracq*Geom_m%Thickness(I,J)*Cp(I,J,Nz_m)
-
-
-             else
-                therm_var_m%Phid(I,J)=0.
-             endif
-          end do
-       end do
-
-
-    case (3)     ! Q_prod_pente :  Pour Essayer D'Avoir La Chaleur En Alpha4
-       !-------------------------------------------------------------------------------------
-       ! Calcul Pour Essayer D'Avoir La Chaleur En Alpha4
-       ! Nom Du Sub Routine Avant Et Var Utilisee  
-       ! Q_prod_pente(T_m%T,T_m%Tpmp,Deform_m%Tobmx,Deform_m%Tobmy,Geom_m%H,Geom_m%Hmx,Geom_m%Hmy,mask_flot_m%Flot,mask_flot_m%Flotmx,mask_flot_m%Flotmy,mask_flot_m%Flgzmx,mask_flot_m%Flgzmy,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Ddx,Ice_flow_m%Ddy,&
-       !     Ice_flow_m%Ddbx,Ice_flow_m%Ddby,Deform_m%Epsxx,Deform_m%Epsyy,Deform_m%Epsxy,mask_flot_m%Gzmx,mask_flot_m%Gzmy,Deform_m%Btt,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid,Deform_m%Glen,Deform_m%Visc)
-
-
-       do J=2,Ny_m
-          do I=2,Nx_m
-             Pente2_maj(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I+1,J))**2 + &
-                  (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I,J+1))**2
-             Pente2_maj(I,J)=Pente2_maj(I,J)*0.25    !Pour La Moyenne Sur Sdx Et Sdy
-          end do
-       end do
-
-       do K=2,Nz_m  
-          do J=2,Ny_m-1
-             do I=2,Nx_m-1
-                ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
-                ! Pour Les Ice-Streams Et Ice-Shelves
-                ! Les Divers Eps Sont CalculÃ©s Dans Strain-Rate Pour L'Ensemble De La Grille
-                if ((mask_flot_m%Flottant(I,J).or.mask_flot_m%Flgz_mx(I,J).or.mask_flot_m%Flgz_mx(I+1,J)).or.  &
-                     (mask_flot_m%Flgz_my(I,J).or.mask_flot_m%Flgz_my(I,J+1))) then
-
-                   Chal2_x(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_xx(I,J)**2
-                   Chal2_y(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_yy(I,J)**2
-                   Chal2_z(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(-Deform_m%Veloc_Deform_xx(I,J)-Deform_m%Veloc_Deform_yy(I,J))**2
-
-                   !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
-                   Chal2_xy(I,J,K)=(Deform_m%Veloc_Deform_xy(I,J)+Deform_m%Veloc_Deform_xy(I+1,J)+Deform_m%Veloc_Deform_xy(I+1,J+1)+Deform_m%Veloc_Deform_xy(I,J+1))
-                   Chal2_xy(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
-
-                else   ! Glace Posee
-                   Chal2_x(I,J,K)=0.
-                   Chal2_y(I,J,K)=0.
-                   Chal2_z(I,J,K)=0.
-                   Chal2_xy(I,J,K)=0.
-                endif
-             end do
-          end do
-       end do
-
-       !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
-       !  Do L=1,size(Deform_m%Btt,4)!N1poly,N2poly
-       !     Do J=2,Ny_m
-       !        Do I=2,Nx_m
-       !
-       !          Ffx A 3 Dimensions !  
-       !           Ffx(I,J,L)=Ddx(I,J,L)*Sdx(I,J)*Sdx(I,J)
-       !           Ffy(I,J,L)=Ddy(I,J,L)*Sdy(I,J)*Sdy(I,J)
-       !        End Do
-       !     End Do
-       !  End Do
-
-       do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
-          do K=2,Nz_m
-             do J=2,Ny_m
-                do I=2,Nx_m
-                   if ((.not.mask_flot_m%Flot_mx(I,J)).and.(.not.mask_flot_m%Gz_mx(I,J))) then
-                      Chalx(I,J,K,L)=(Deform_m%Btt_Deform(I-1,J,K,L)+Deform_m%Btt_Deform(I,J,K,L))*Ice_flow_m%Difus_x(I,J,L)
-
-                   else if (mask_flot_m%Gz_mx(I,J)) then ! Ice Streams
-                      Chalx(I,J,K,L)=0.
-
-                   else                     ! Ice Shelves  
-                      Chalx(I,J,K,L)=0.
-
-                   endif
-
-                   if ((.not.mask_flot_m%Flot_my(I,J)).and.(.not.mask_flot_m%Gz_my(I,J))) then
-                      Chaly(I,J,K,L)=(Deform_m%Btt_Deform(I,J-1,K,L)+Deform_m%Btt_Deform(I,J,K,L))*Ice_flow_m%Difus_y(I,J,L)
-
-                   else if (mask_flot_m%Gz_my(I,J)) then  ! Ice Streams
-                      Chaly(I,J,K,L)=0.
-
-                   else                      ! Ice Shelves
-                      Chaly(I,J,K,L)=0.
-                   endif
-
-                end do
-             end do
-          end do
-       end do
-
-       !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
-       !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
-       !         Les Productions X Et Y
-       !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy) 
-
-       do K=2,Nz_m
-          do J=1,Ny_m-1
-             do I=1,Nx_m-1
-
-                ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
-                Chaldef_maj(I,J,K)= 0.
-
-                ! Chalk_2 Pour Ice Shelves Et Ice Streams
-                Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
-                     Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
-
-                ! Chalk_1 Pour La Partie PosÃ©e
-                do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
-
-                   ! On Somme La Chaleur Due Aux Diverses Lois De DÃ©formation Et Celles En X Et Y
-                   ! La Somme Se Fait Par (Cx2+Cy2)** 0.5. Le 4 Vient Des Moyennes Deform_m%Btt + Des Moyennes Gauche-Droite
-
-                   ! On Fait La Moyenne Des Termes Ddx*Btt (*0.25 Pour Cette Moyenne, Le 0.5 Est Pour Les Btt)
-                   Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))+ &
-                        (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
-
-                   Chalk_1=Chalk_1*0.25*0.5
-
-                   ! On Multiplie Par La Pente Moyenne Au Carre Sur Le Noeud Majeur
-                   Chalk_1=Chalk_1*Pente2_maj(I,J)
-
-                   Chalk_1=Ro*G*Chalk_1*(Ee(K)**(Deform_m%Glen_exp(L)+1))/Cp(I,J,K)  ! Attention Plus De /4
-                   Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
-                enddo
-
-                ! Pour Shelves Et Streams,  On Ajoute Chalk_2
-                Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
-             end do
-          end do
-       end do
-
-       ! Chaleur Produite A La Base Par Le Glissement
-
-       do J=2,Ny_m
-          do I=2,Nx_m
-
-             if (mask_flot_m%Gz_mx(I,J)) then
-                Chalglissx(I,J)= abs(Ice_flow_m%U_Column_x(I,J)*Deform_m%Cisail_Basal_mx(I,J))
-             else
-                Chalglissx(I,J)=Ice_flow_m%Difus_bx(I,J)*Geom_m%Slop_x(I,J)**2*Ro*G*Geom_m%Thick_mx(I,J)
-             endif
-
-             if (mask_flot_m%Gz_my(I,J)) then
-                Chalglissy(I,J)= abs(Ice_flow_m%U_Column_y(I,J)*Deform_m%Cisail_Basal_my(I,J))
-             else
-                Chalglissy(I,J)=Ice_flow_m%Difus_by(I,J)*Geom_m%Slop_y(I,J)**2*Ro*G*Geom_m%Thick_my(I,J)
-             endif
-
-          end do
-       end do
-
-
-       !  Boundary Condition Ice-Rock Interface 
-
-       K=Nz_m
-
-       do J=2,Ny_m-1
-          do I=2,Nx_m-1
-
-             !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
-             !        Dans La Derniere 1/2 Maille Et Par Le Glissement
-             !        Attention Phid Est >0 Et Ghf Est <0
-             if (.not.mask_flot_m%Flottant(I,J)) then
-
-                !           Phid Avec Fonte Sous Les Streams
-                !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
-                !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
-                !                        Chalbed*Fracq
-
-
-                ! Moyenne Phid Sur 4 Points Formulation (A2+B2)**0.5  ! Je Garde Pour L'Instant A Revoir
-                !
-                Chalgliss_maj(I,J)=(Chalglissx(I,J)+Chalglissx(I+1,J))**2+  &
-                     (Chalglissy(I,J)+Chalglissy(I,J+1))**2
-
-                Chalgliss_maj(I,J)=Chalgliss_maj(I,J)**0.5
-
-                Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.5 ! Les Moyennes Droite Gauche
-
-                ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
-                Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T_m%Temperature(I,J,Nz_m)-T_m%Temp_melting(I,J,Nz_m))*Ecart_phid)
-
-
-                ! Flux Total A Rajouter a La Base
-
-                therm_var_m%Phid(I,J)=Chalgliss_maj(I,J)+Chaldef_maj(I,J,Nz_m)*Fracq*Geom_m%Thickness(I,J)*Cp(I,J,Nz_m)
-
-
-             else
-                therm_var_m%Phid(I,J)=0.
-             endif
-          end do
-       end do
-
-
-    case (4)   ! Q_u_taub : Routine Qui concentre La Chaleur  au Socle
-       !-------------------------------------------------------------------------------------
-       ! Routine Qui Prend La Chaleur Remise Au Socle
-       ! Nom Du Sub Routine Avant Et Var Utilisee  
-       ! Q_u_taub(T_m%T,T_m%Tpmp,Geom_m%H,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid)
-
-
-       do J=2,Ny_m
-          do I=2,Nx_m
-
-
-             Pente_maj(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I+1,J))**2 + &          ! Pente
-                  (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I,J+1))**2
-
-             Pente_maj(I,J)=Pente_maj(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
-             Pente_maj(I,J)=Pente_maj(I,J)**0.5
-
-             Vit_maj(I,J)=(Ice_flow_m%U_Column_x(I,J)+Ice_flow_m%U_Column_x(I+1,J))**2 + &        ! Vitesse De Bilan
-                  (Ice_flow_m%U_Column_y(I,J)+Ice_flow_m%U_Column_y(I,J+1))**2
-             Vit_maj(I,J)=Vit_maj(I,J)*0.25
-             Vit_maj(I,J)=Vit_maj(I,J)**0.5
-
-             Uslid_maj(I,J)=(Ice_flow_m%Veloc_x(I,J,Nz_m)+Ice_flow_m%Veloc_x(I+1,J,Nz_m))**2 + &        ! Vitesse De Bilan
-                  (Ice_flow_m%Veloc_y(I,J,Nz_m)+Ice_flow_m%Veloc_y(I,J+1,Nz_m))**2
-             Uslid_maj(I,J)=Uslid_maj(I,J)*0.25
-             Uslid_maj(I,J)=Uslid_maj(I,J)**0.5
-
-          end do
-       end do
-
-
-
-       Chaldef_maj(:,:,:)=0.
-       Chalgliss_maj(:,:)=Ro*G*Geom_m%Thickness(:,:)*Pente_maj(:,:)*Uslid_maj(:,:)
-       Chaldef_maj(:,:,Nz_m)=Ro*G*Geom_m%Thickness(:,:)*Pente_maj(:,:)*Vit_maj(:,:)-Chalgliss_maj(:,:)
-
-       ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
-
-       Chalgliss_maj(:,:)=Chalgliss_maj(:,:)*exp((T_m%Temperature(:,:,Nz_m)-T_m%Temp_melting(:,:,Nz_m))*Ecart_phid)
-
-       ! Flux Total A Rajouter a Base
-       therm_var_m%Phid(:,:)=Chaldef_maj(:,:,Nz_m)+Chalgliss_maj(:,:)
-
-
-    case (5)    ! Q_u_taub_stag : concentre La Chaleur Au Socle mais Sur Les Mailles Staggered
-       !----------------------------------------------------------------------------------
-       ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
-       ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
-       ! Nom Du Sub Routine Avant Et Var UtilisÃ©  
-       ! Q_u_taub_stag(T_m%T,T_m%Tpmp,Geom_m%H,Geom_m%Slop_x,Geom_m%Slop_y,Ice_flow_m%Uxbar,Ice_flow_m%Uybar,therm_var_m%Phid)
-
-
-       do J=2,Ny_m
-          do I=2,Nx_m
-
-             Pente_stag(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I,J-1))**2 + &          ! Pente
-                  (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I-1,J))**2
-
-             Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
-             Pente_stag(I,J)=Pente_stag(I,J)**0.5
-
-             Vit_stag2(I,J)=(Ice_flow_m%U_Column_x(I,J)+Ice_flow_m%U_Column_x(I,J-1))**2 + &        ! Vitesse De Bilan
-                  (Ice_flow_m%U_Column_y(I,J)+Ice_flow_m%U_Column_y(I-1,J))**2
-             Vit_stag2(I,J)=Vit_stag2(I,J)*0.25
-             Vit_stag2(I,J)=Vit_stag2(I,J)**0.5
-
-             Uslid_stag(I,J)=(Ice_flow_m%Veloc_x(I,J,Nz_m)+Ice_flow_m%Veloc_x(I,J-1,Nz_m))**2 + &        ! Vitesse De Bilan
-                  (Ice_flow_m%Veloc_y(I,J,Nz_m)+Ice_flow_m%Veloc_y(I-1,J,Nz_m))**2
-             Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
-             Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
-
-             Hmxy(I,J)=((Geom_m%Thickness(I,J)+Geom_m%Thickness(I-1,J-1))+(Geom_m%Thickness(I,J-1)+Geom_m%Thickness(I-1,J)))*0.25
-
-          end do
-       end do
-
-
-       Qslid(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Uslid_stag(:,:)
-       Qdef2(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Vit_stag2(:,:)-Qslid(:,:)
-       Chaldef_maj(:,:,:)=0.
-
-
-       ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
-
-       do J=2,Ny_m-1
-          do I=2,Nx_m-1
-             Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
-             Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
-
-             Chaldef_maj(I,J,Nz_m)=(Qdef2(I,J)+Qdef2(I+1,J+1))+(Qdef2(I,J+1)+Qdef2(I+1,J))
-             Chaldef_maj(I,J,Nz_m)=Chaldef_maj(I,J,Nz_m)*0.25
-
-             ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
-             Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T_m%Temperature(I,J,Nz_m)-T_m%Temp_melting(I,J,Nz_m))*Ecart_phid)
-
-          end do
-       end do
-
-
-
-
-       !
-       !Write(126,*)'Time=',Time
-       !J=41
-       !Do I=20,30
-       !   Write(126,'(I3,6(E14.4,1x))') I,Chalgliss_maj(I,J),Chaldef_maj(I,J,Nz_m),Qdef2(I,J),Qdef2(I+1,J+1) &
-       !       , Qdef2(I,J+1),Qdef2(I+1,J)
-       !End Do
-
-       ! Flux Total A Rajouter Ã La Base
-       therm_var_m%Phid(:,:)= Chaldef_maj(:,:,Nz_m)+Chalgliss_maj(:,:)
-
-    case (6)    ! Q_sia_stag  : concentre La Chaleur  Au Socle Mais Sur Les Mailles Staggered
-       ! Au Milieux Des Mailles Vitesses
-
-       !----------------------------------------------------------------------------------
-       ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
-       ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
-       ! Nom Du Sub Routine Avant Et Var UtilisÃ©  
-       ! Q_sia_stag(T_m%T,T_m%Tpmp,Geom_m%H,Geom_m%Slop_x,Geom_m%Slop_y,therm_var_m%Phid)
-
-
-       Pente_stag = 0.
-       Vit_stag3   = 0.
-       Uslid_stag = 0.
-       Hmxy  = 0.
-       Qslid = 0.
-       Qdef3 = 0.
-
-       do J=2,Ny_m
-          do I=2,Nx_m
-
-             !     Variables 2d
-             Pente_stag(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I,J-1))**2 + &          ! Pente
-                  (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I-1,J))**2
-
-             Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
-             Pente_stag(I,J)=Pente_stag(I,J)**0.5
-
-             Uslid_stag(I,J)=(Ice_flow_m%Veloc_x(I,J,Nz_m)+Ice_flow_m%Veloc_x(I,J-1,Nz_m))**2 + &        ! Vitesse De Bilan
-                  (Ice_flow_m%Veloc_y(I,J,Nz_m)+Ice_flow_m%Veloc_y(I-1,J,Nz_m))**2
-             Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
-             Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
-
-             Hmxy(I,J)=((Geom_m%Thickness(I,J)+Geom_m%Thickness(I-1,J-1))+(Geom_m%Thickness(I,J-1)+Geom_m%Thickness(I-1,J)))*0.25
-
-             !     En Vertical : Calcul De La Deformation Par Differentiation Des Vitesses
-
-             do K=1,Nz_m
-
-                Vit_stag3(I,J,K)=(Ice_flow_m%Veloc_x(I,J,K)+Ice_flow_m%Veloc_x(I,J-1,K))**2 + &        ! Magnitude Vitesse 
-                     (Ice_flow_m%Veloc_y(I,J,K)+Ice_flow_m%Veloc_y(I-1,J,K))**2            ! A Tous Niveaux K
-
-                Vit_stag3(I,J,K)=Vit_stag3(I,J,K)*0.25
-                Vit_stag3(I,J,K)=Vit_stag3(I,J,K)**0.5
-             end do
-
-             Qdef3(I,J,1)=0.
-             do K=2,Nz_m-1
-                Qdef3(I,J,K)=Vit_stag3(I,J,K-1)-Vit_stag3(I,J,K+1)             ! Difference Des Vitesses
-                Qdef3(I,J,K)=Ro*G*Pente_stag(I,J)*Ee(K)*Qdef3(I,J,K)/2./Dee    ! Gamma Tau
-             end do
-
-
-             Qdef3(I,J,Nz_m)=Vit_stag3(I,J,Nz_m-1)-Vit_stag3(I,J,Nz_m)         ! Pour Le Fond Differentiation
-             Qdef3(I,J,Nz_m)=Ro*G*Pente_stag(I,J)*Qdef3(I,J,K)/Dee             ! Sur Une Seule Maille
-          end do
-       end do
-
-
-       Qslid(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Uslid_stag(:,:)
-
-       ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
-
-       do J=2,Ny_m-1
-          do I=2,Nx_m-1
-             Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
-             Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
-
-             ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
-             Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T_m%Temperature(I,J,Nz_m)-T_m%Temp_melting(I,J,Nz_m))*Ecart_phid)
-
-             ! Chaleur Due A La Defromation
-             do K=1,Nz_m
-                Chaldef_maj(I,J,K)=(Qdef3(I,J,K)+Qdef3(I+1,J+1,K))+(Qdef3(I,J+1,K)+Qdef3(I+1,J,K))
-                Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)*0.25
-                Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)/Cp(I,J,K)
-             end do
-
-
-          end do
-       end do
-
-       !Write(126,*)'Time=',Time
-       !J=41
-       !Do I=20,30
-       !   Write(126,'(I3,6(E14.4,1x))') I,Chalgliss_maj(I,J),Chaldef_maj(I,J,Nz_m),Qdef2(I,J),Qdef2(I+1,J+1) &
-       !       , Qdef2(I,J+1),Qdef2(I+1,J)
-       !End Do
-
-       ! Flux Total A Rajouter Ã La Base
-
-       therm_var_m%Phid(:,:)=Chalgliss_maj(:,:)+Chaldef_maj(:,:,Nz_m)*Fracq*Geom_m%Thickness(:,:)*Cp(:,:,Nz_m)
-
-    case (7)   !  Q_all_stag idem case 6 mais Le Calcul De Taub Tient Compte Des Noeuds Grzmx
-       !----------------------------------------------------------------------------------
-       ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
-       ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
-       ! Le Calcul De Taub Tient Compte Des Noeuds Grzmx     
-       ! Nom Du Sub Routine Avant Et Var UtilisÃ©  
-       ! Q_all_stag(T_m%T,T_m%Tpmp,Deform_m%Tobmx,Deform_m%Tobmy,Geom_m%H,Geom_m%Hmx,Geom_m%Hmy,mask_flot_m%Flgzmx,mask_flot_m%Flgzmy,Geom_m%Slop_x,Geom_m%Slop_y,therm_var_m%Phid)
-
-       where (mask_flot_m%Flgz_mx(:,:))
-          Tox(:,:)=Deform_m%Cisail_Basal_mx(:,:)
-       elsewhere
-          Tox(:,:)=Geom_m%Slop_x(:,:)*Geom_m%Thick_mx(:,:)*Ro*G
-       end where
-
-       where (mask_flot_m%Flgz_my(:,:))
-          Toy(:,:)=Deform_m%Cisail_Basal_my(:,:)
-       elsewhere
-          Toy(:,:)=Geom_m%Slop_y(:,:)*Geom_m%Thick_my(:,:)*Ro*G
-       end where
-
-       do J=2,Ny_m
-          do I=2,Nx_m
-
-             !     Variables 2d
-             Pente_stag(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I,J-1))**2 + &    ! Pente
-                  (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I-1,J))**2
-
-             Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
-             Pente_stag(I,J)=Pente_stag(I,J)**0.5
-
-
-             Tob_stag(I,J)=(Tox(I,J)+Tox(I,J-1))**2 + &          ! Pente
-                  (Toy(I,J)+Toy(I-1,J))**2
-
-             Tob_stag(I,J)=Tob_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
-             Tob_stag(I,J)=Tob_stag(I,J)**0.5
-
-             Uslid_stag(I,J)=(Ice_flow_m%Veloc_x(I,J,Nz_m)+Ice_flow_m%Veloc_x(I,J-1,Nz_m))**2 + &        ! Vitesse De Bilan
-                  (Ice_flow_m%Veloc_y(I,J,Nz_m)+Ice_flow_m%Veloc_y(I-1,J,Nz_m))**2
-             Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
-             Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
-
-             Hmxy(I,J)=((Geom_m%Thickness(I,J)+Geom_m%Thickness(I-1,J-1))+(Geom_m%Thickness(I,J-1)+Geom_m%Thickness(I-1,J)))*0.25
-
-             !     En Vertical : Calcul De La Deformation Par Differentiation Des Vitesses
-
-             do K=1,Nz_m
-
-                Vit_stag3(I,J,K)=(Ice_flow_m%Veloc_x(I,J,K)+Ice_flow_m%Veloc_x(I,J-1,K))**2 + &        ! Magnitude Vitesse 
-                     (Ice_flow_m%Veloc_y(I,J,K)+Ice_flow_m%Veloc_y(I-1,J,K))**2            ! A Tous Niveaux K
-
-                Vit_stag3(I,J,K)=Vit_stag3(I,J,K)*0.25
-                Vit_stag3(I,J,K)=Vit_stag3(I,J,K)**0.5
-             end do
-
-             Qdef3(I,J,1)=0.
-             do K=2,Nz_m-1
-                Qdef3(I,J,K)=Vit_stag3(I,J,K-1)-Vit_stag3(I,J,K+1)            ! Difference Des Vitesses
-                Qdef3(I,J,K)=Ro*G*Pente_stag(I,J)*Ee(K)*Qdef3(I,J,K)/2./Dee   ! Gamma Tau
-             end do
-
-
-             Qdef3(I,J,Nz_m)=Vit_stag3(I,J,Nz_m-1)-Vit_stag3(I,J,Nz_m)        ! Pour Le Fond Differentiation
-             Qdef3(I,J,Nz_m)=Ro*G*Pente_stag(I,J)*Qdef3(I,J,K)/Dee            ! Sur Une Seule Maille
-          end do
-       end do
-
-       Qslid(:,:)=Tob_stag*Uslid_stag(:,:)
-
-       ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
-
-       do J=2,Ny_m-1
-          do I=2,Nx_m-1
-             Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
-             Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
-
-             ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
-             Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*exp((T_m%Temperature(I,J,Nz_m)-T_m%Temp_melting(I,J,Nz_m))*Ecart_phid)
-
-             ! Chaleur Due A La Defromation
-             do K=1,Nz_m         
-                Chaldef_maj(I,J,K)=(Qdef3(I,J,K)+Qdef3(I+1,J+1,K))+(Qdef3(I,J+1,K)+Qdef3(I+1,J,K))
-                Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)*0.25
-                Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)/Cp(I,J,K)
-             end do
-
-          end do
-       end do
-
-       !
-       !Write(126,*)'Time=',Time
-       !J=41
-       !Do I=20,30
-       !   Write(126,'(I3,6(E14.4,1x))') I,Chalgliss_maj(I,J),Chaldef_maj(I,J,Nz_m),Qdef2(I,J),Qdef2(I+1,J+1) &
-       !       , Qdef2(I,J+1),Qdef2(I+1,J)
-       !End Do
-
-       ! Flux Total A Rajouter Ã La Base
-       therm_var_m%Phid(:,:)=Chalgliss_maj(:,:)+Chaldef_maj(:,:,Nz_m)*Fracq*Geom_m%Thickness(:,:)*Cp(:,:,Nz_m)
-
-    end select
-  end subroutine Qprod
-
-  !---------------------------------------------------------------------------------------------------------
-
-end module Interface_qprod_icetemp
Index: unk/SOURCES/Temperature-routines/interface_advec_icetemp.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/interface_advec_icetemp.f90	(revision 23)
+++ 	(revision )
@@ -1,142 +1,0 @@
-!> \file interface_advec_icetemp.f90
-!! Interfaces of functions and subroutines for temperature advection in the heat transfer equation
-!<
-
-!> \namespace interface_advec_icetemp
-!! Interfaces of functions and subroutines to calculate the temperature advection 
-!<
-module interface_advec_icetemp
-
-  !> \interface Advec_icetemp
-  !! Interfaces of functions and subroutines to calculate the temperature advection 
-  !! \author 
-  !! \date 2010
-  !<
-  Interface  Advec_icetemp
-     module procedure Advec_1d_col 
-  End Interface
-
-contains
-
-  !> SUBROUTINE: Advec_icetemp_3d
-  !! Implementation of Routines to calculate advection in the heat transfer equation
-  !! for a column
-  !! Used modules:
-  !!    - use Icetemp_declar
-  !!    - use geom_type
-  !!    - use temperature_type
-  !!    - use Ice_flow_type
-  !!    - use Autre_pr_temp_type
-  !!    - use step_temp_type
-  !>
-
-  Subroutine  Advec_1d_col (Ii,Jj,Prodq_m,Cro_m,Advecx_m,Advecy_m,Advec_m,Ct_m,   &
-       geom_m,T_m,Ice_flow_m,Autre_pr_temp_m,step_temp_m,debug_3D_m)
-
-
-    use icetemp_declar 
-    use geom_type
-    use temperature_type
-    use Ice_flow_type
-    use Autre_pr_temp_type
-    use step_temp_type
-
-    Implicit None
-    !< Arguments
-    Integer,                    Intent(In)  ::  Ii,Jj          !< Indice Selon X Et Y
-    Real, Dimension(:),         Intent(In)   :: Prodq_m        !< Tableau 1d Vert. heat production
-    Real, Dimension(:),         Intent(In)   :: Cro_m          !< Tableau 1d Vert. heat capcity
-    Real, Dimension(:),         Intent(Inout):: Advecx_m       !< Advection selon x
-    Real, Dimension(:),         Intent(Inout):: Advecy_m       !< Advection selon y
-    Real, Dimension(:),         Intent(Inout):: Advec_m        !< Advection total
-    Real, Dimension(:),         Intent(In)   :: Ct_m           !< Tableau 1d Vert.  thermal cond.
-
-    type (Geom),               intent(in)    :: geom_m         !< Variable geometrie 
-    type (Ice_flow),           intent(in)    :: Ice_flow_m     !< Variable ice_flow  
-    type (Temperature),        intent(inout) :: T_m            !< Variable temperature 
-    type (Autre_pr_temp),      intent(inout) :: Autre_pr_temp_m  
-    type (step_temp),          intent(in)    :: step_temp_m    
-    
-    real,dimension(:,:,:),intent(inout) :: debug_3D_m          !< debug array
-
- 
-    Do K=2,Nz_m-1
-
-
-       ! On Remplit Deux Tableaux Pour eviter Des Problemes D'Asymetrie
-       !         Avection Selon X  (advecx_m en general > 0)
-       !         ----------------
-       Advecx_m(K) = Iadvec_w(Ii,Jj) * Ice_flow_m%Veloc_x(Ii,Jj,K) *          & ! ux west if upwind
-            (T_m%Temperature(Ii,Jj,K) - T_m%Temperature(Ii-1,Jj,K))  & ! west T gradient
-            
-            + Iadvec_e(Ii+1,Jj)* Ice_flow_m%Veloc_x(Ii+1,Jj,K)*      & ! ux east if upwind
-            (T_m%Temperature(Ii+1,Jj,K)-T_m%Temperature(Ii,Jj,K))     ! east T gradient
-
-       Advecx_m(K) =  Advecx_m(K) * Dx11                                        !Dx11=1/Dx
-
-
-       !         Avection Selon Y
-       !         ----------------
-       Advecy_m(K) =   Iadvec_s(Ii,Jj) * Ice_flow_m%Veloc_y(Ii,Jj,K) *        & ! uy sud if upwind
-            (T_m%Temperature(Ii,Jj,K) - T_m%Temperature(Ii,Jj-1,K))& ! south T gradient
-            
-            + Iadvec_n(Ii,Jj+1) * Ice_flow_m%Veloc_y(Ii,Jj+1,K) *   & ! uy nord is upwind
-            (T_m%Temperature(Ii,Jj+1,K)-T_m%Temperature(Ii,Jj,K))    ! north T gradient
-
-       Advecy_m(K) =  Advecy_m(K) * Dx11
-
-       Advec_m(K)=Advecx_m(K)+Advecy_m(K)  ! Advection Totale
-
-    End Do
-
-    !     Fin Boucle K, Ecriture Matrice Tridiag Dans La Glace
-
-
-    !      -----------------------------Cas General (H>10m)
-    thick_ice:  If (Geom_m%Thickness(Ii,Jj).Gt.10.) Then
-
-       !        Variables De Calcul Dans La Glace
-       ! dou = dtt/dz^2
-       Dou=step_temp_m%Dtt/Dee/Dee/ Geom_m%Thickness(Ii,Jj) / Geom_m%Thickness(Ii,Jj)  
-
-       ! Dah = dtt/dz
-       Dah=step_temp_m%Dtt /Dee /Geom_m%Thickness(Ii,Jj)
-
-       ! thermal conductivity at mid point just below the surface                   
-       Ct_bas=2.*(Ct_m(1)*Ct_m(2))/(Ct_m(1)+Ct_m(2))           
-
-       ! surface temperature : cannot go above 0 celsius
-       T_m%Temperature(Ii,Jj,1)=Min(0.,Autre_pr_temp_m%Ts(Ii,Jj))
-
-       ! surface boundary condition
-       Aa(1)=0.
-       Bb(1)=1.
-       Cc(1)=0.
-       Rr(1)=T_m%Temperature(Ii,Jj,1)
-
-       Do K=2,Nz_m-1
-          Dzz=Dou/Cro_m(K)
-
-          ! Conductivite Au Milieu Des Mailles
-          Ct_haut=Ct_bas
-          Ct_bas=2.*(Ct_m(K)*Ct_m(K+1))/(Ct_m(K)+Ct_m(K+1))
-
-          ! Advection Verticale Centree
-          Aa(K) = -Dzz * Ct_haut - Ice_flow_m%Veloc_z_r(Ii,Jj,K) * Dah / 2.  ! lower diag
-          Cc(K) = -Dzz * Ct_bas  + Ice_flow_m%Veloc_z_r(Ii,Jj,K) * Dah / 2.  ! upper diag
-          Bb(K) =  1.+ Dzz * (Ct_haut+Ct_bas)                                ! diag
-
-          Rr(K) =  T_m%Temperature(Ii,Jj,K) + Prodq_m(K) * step_temp_m%Dtt   ! vector
-          Rr(K) =  Rr(K)- Dttdx * (Advec_m(K))
-
-       End DO
-    End If thick_ice
-
-
-  End Subroutine Advec_1d_col
-
-  !-------------------------------------------------------------------------
-
-  
-end module interface_advec_icetemp
-
Index: unk/SOURCES/Temperature-routines/interface_icetempmod.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/interface_icetempmod.f90	(revision 23)
+++ 	(revision )
@@ -1,371 +1,0 @@
-!> \file interface_icetempmod.f90
-!! Interfaces of functions and subroutines for temperature calculation 
-!<
-
-!> \namespace interface_icetempmod
-!! Interfaces of functions and subroutines for the temperature calculation 
-!<
-
-module interface_icetempmod
-
-  !> \interface icetemp_mod
-  !! Interfaces of functions and subroutines for the temperature calculation 
-  !!\author ...
-  !!\date 2010
-  !<
-
-  Interface init_icetemp
-     module procedure init_temp
-  End Interface
-  
-  Interface icetemp
-     module procedure ice_temp
-  End Interface
-
-contains
-
-  !> SUBROUTINE: Init_temp 
-  !! Initialise the temperature calculation
-  !! Used modules:
-  !! - Icetemp_declar
-  !>
-
-  Subroutine Init_temp(Num_rep_42_m)
-    Use Icetemp_declar
-
-    Implicit None
-
-    !< Arguments
-    Integer,                       Intent(In)::Num_rep_42_m !< Id Of Reponse
-    ! Character(Len=80) :: Filin
-    !________________________________________________________________
-    ! La Temperature Du Point De Congelation De L'Eau De Mer
-    ! Est Tiree De Jenkins (1991), Tb Est En Degres C
-    ! _______________________________________________________________
-
-    !  Nouvelle Formulation De Fracq
-    !  Nfracq=15  ! Ancienne Formulation
-
-    Nfracq=1
-
-    Write(Num_rep_42_m,*) 'Calcul Des Temperatures'
-    Write(Num_rep_42_m,*) 'Nfracq=',Nfracq
-
-    Iq=1
-    If (Iq.Eq.1) Write(Num_rep_42_m,*) 'Chaleur Demi Maille, Iq=',Iq
-    If (Iq.Eq.2) Write(Num_rep_42_m,*) 'Chaleur Centree, Iq=',Iq
-    If (Iq.Eq.3) Write(Num_rep_42_m,*) 'Chaleur Pente, Iq=',Iq
-    If (Iq.Eq.4) Write(Num_rep_42_m,*) 'Chaleur U_taub, Iq=',Iq
-    If (Iq.Eq.5) Write(Num_rep_42_m,*) 'Chaleur U_taub_stag, Iq=',Iq
-    If (Iq.Eq.6) Write(Num_rep_42_m,*) 'Chaleur Q_sia_stag,Iq=',Iq
-    If (Iq.Eq.7) Write(Num_rep_42_m,*) 'Chaleur Q_all_stag,Iq=',Iq
-    Write(Num_rep_42,*)
-
-    Ecart_phid=0.5
-    Write(Num_rep_42_m,*) ' Si Base Froide, La Chaleur De Glissement Est En Exp(Delta T * Ecart_phid)'
-    Write(Num_rep_42_m,*) 'Ecart_phid=',Ecart_phid
-    Write(Num_rep_42_m,*) '---------------------------------------------------------------------------'
-
-  End Subroutine Init_temp
-  
-
-  !-----------------------------------------------------------------------------------------------------------------------
-
-  !> SUBROUTINE: Ice_temp()
-  !! Calculate the temperature
-  !!
-  !! Used modules:
-  !!  -  use Icetemp_declar 
-  !!  -  use Interface_prop_th_icetemp
-  !!  -  use Interface_qprod_icetemp
-  !!  -  use Interface_advec_icetemp
-  !!  -  use Interface_temp_col
-  !!  -  use geom_type
-  !!  -  use temperature_type
-  !!  -  use ice_flow_type
-  !!  -  use mask_flgz_type
-  !!  -  use deformation_type
-  !!  -  use autre_pr_temp_type
-  !!  -  use step_temp_type 
-  !>
-
-  Subroutine ice_temp(geom_m,temperature_m,Ice_flow_m,Mask_flgz_m,&
-       deformation_m,Autre_pr_temp_m,step_temp_m,num_tracebug_m,debug_3D_m)
-
-    use Icetemp_declar
-    use Interface_prop_th_icetemp
-    use Interface_qprod_icetemp
-    use Interface_advec_icetemp
-    use Interface_temp_col
-
-    use geom_type
-    use temperature_type
-    use ice_flow_type
-    use mask_flgz_type
-    use deformation_type
-    use autre_pr_temp_type
-    use step_temp_type 
-
-    Implicit None 
-
-    !< Arguments
-
-    type (Geom),          intent(in)    :: geom_m             
-    type (Temperature),   intent(inout) :: temperature_m 
-    type (Ice_flow),      intent(in)    :: Ice_flow_m        
-    type (Mask_flgz),    intent(in)     :: Mask_flgz_m   
-    type (Deformation),   intent(in)    :: deformation_m
-    type (Autre_pr_temp), intent(inout) :: Autre_pr_temp_m
-    type (step_temp),     intent(in)    :: step_temp_m
-    integer,              intent(in)    :: num_tracebug_m
-    real,dimension(:,:,:),intent(inout) :: debug_3D_m
-
-    Nx_m=size(Temperature_m%Temp_melting,1)
-    Ny_m=size(Temperature_m%Temp_melting,2)
-    Nz_m=size(Temperature_m%Temp_melting,3)
-    NzNzm_m=size(Temperature_m%Temperature,3)
-    Nzm_m= NzNzm_m-Nz_m
-
-    ! Allocation de la memoire pour les tableaux dans icetemp_declar
-
-    If (.Not. Allocated(Ct)) then 
-       Allocate(Ct(Nx_m,Ny_m,Nz_m))
-    end if
-
-
-    If (.Not. Allocated(Aa)) then 
-       ! tab Nn
-       Allocate (Aa(Nn),Bb(Nn),Cc(Nn),Rr(Nn),Hh(Nn))  
-       ! tab Nznzm_m
-       Allocate(Tdot(Nznzm_m))   
-       ! tab Nzm+1
-       Allocate(Abis(Nzm_m+1),Bbis(Nzm_m+1),Cbis(Nzm_m+1),Rbis(Nzm_m+1),Hbis(Nzm_m+1))                  
-       ! tab nz_m
-       Allocate(Ee(Nz_m))        
-       !Tab Nx_m,Ny_m
-       Allocate(Tbmer(Nx_m,Ny_m),&   
-            Alpha(Nx_m,Ny_m),&                   
-            Ubar(Nx_m,Ny_m),&                    
-            Chalglissx(Nx_m,Ny_m),&
-            Chalglissy(Nx_m,Ny_m),&
-            Iadvec_w(Nx_m,Ny_m),&
-            Iadvec_e(Nx_m,Ny_m),&
-            Iadvec_s(Nx_m,Ny_m),&
-            Iadvec_n(Nx_m,Ny_m) )      
-       ! Tab Nx_m,Ny_m,Nz_m
-       Allocate(Chal2_x(Nx_m,Ny_m,Nz_m),&
-            Chal2_y(Nx_m,Ny_m,Nz_m),&
-            Chal2_z(Nx_m,Ny_m,Nz_m),&
-            Chal2_xy(Nx_m,Ny_m,Nz_m),&
-            Chaldef_maj(Nx_m,Ny_m,Nz_m),&
-            Advecx(Nx_m,Ny_m,Nz_m),&
-            Advecy(Nx_m,Ny_m,Nz_m),&
-            Advec(Nx_m,Ny_m,Nz_m))
-       ! Tab Nx_m,Ny_m,Nz_m,N1poly:N2poly
-       Allocate(Chalx(Nx_m,Ny_m,Nz_m,size(Deformation_m%Btt_Deform,4)),&
-            Chaly(Nx_m,Ny_m,Nz_m,size(Deformation_m%Btt_Deform,4)))         
-       ! Tab Nx_m,Ny_m,N1poly:N2poly
-       Allocate(Ffx( Nx_m,Ny_m,size(Deformation_m%Btt_Deform,4)),&
-            Ffy( Nx_m,Ny_m,size(Deformation_m%Btt_Deform,4)))                                        
-       ! Tab Nx_m,Ny_m,Nznzm_m
-       Allocate( T3d_new (Nx_m,Ny_m,Nznzm_m))                  
-    end If
-
-    ! init
-    Aa=0.    
-    Bb=0.              
-    Cc=1.             
-    Rr=0.              
-    Hh=0.              
-    Tdot=0  
-    !Instruction 
-    Dee=1./(Nz_m-1) 
-    Fracq=(1.-(1.-Dee/2.)**Nfracq)/Nfracq
-    Fracq=(1.-(1.-Dee/2.)**5.)/5.
-    Ee(1)=0.
-    Ee(Nz_m)=1.
-    Do K=1,Nz_m
-       If ((K.Ne.1).And.(K.Ne.Nz_m)) Ee(K)=(K-1.)/(Nz_m-1.)
-    End Do
-
-    !   Variables Dependant Du Pas De Temps Dtt
-    Ctm=step_temp_m%Dtt*Cm/Rom/Cpm/Dzm/Dzm
-    Dttdx=step_temp_m%Dtt/step_temp_m%Dx
-    Dx11=1./step_temp_m%Dx
-
-    If (Itracebug.Eq.1) Write(Num_tracebug_m,*)' Entree Dans Routine Ice_temp Time=',Time
-
-    !________________________________________________________________
-    ! La Temperature Du Point De Congelation De L'Eau De Mer
-    ! Est Tiree De Jenkins (1991), Tb Est En Degres C
-    ! _______________________________________________________________
-    If (Itracebug.Eq.1) Write(Num_tracebug_m,*)' avant appel  Routine Temp_mer'
-    Call Temp_mer(Geom_m%Thickness,Mask_flgz_m%Flottant,debug_3D_m)
-
-    ! Appel Aux Proprietes Thermiques De La Glace : Conductivite, Capacite Claorifique -Temperature_m%Tpmp
-
-    If (Itracebug.Eq.1) Write(Num_tracebug_m,*)' avant appel  Routine Thermal_property'
-    Call Thermal_property (Temperature_m%Temperature,Temperature_m%Temp_melting,&
-         Geom_m%Thickness,debug_3D_m)
-
-    If (Itracebug.Eq.1) Write(Num_tracebug_m,*)' avant appel  Routine Qprod'
-    Call Qprod_ice (Iq,geom_m,temperature_m,Ice_flow_m,&
-         Mask_flgz_m,deformation_m,Autre_pr_temp_m,debug_3D_m)
-
-
-    ! Dans Le Socle : Les elements de La Matrice Tridiag Sont Invariants Dans L'Espace
-    Do K=Nz_m+1,Nznzm_m-1
-       Aa(K)=-Ctm
-       Bb(K)=1.+2.*Ctm
-       Cc(K)=-Ctm
-    End Do
-
-    ! Conditions Aux Limites Pour Le Socle
-    Aa(Nznzm_m)=-1.
-    Bb(Nznzm_m)=1.
-    Cc(Nznzm_m)=0.
-
-    ! Calcul D'Advection : On Traite Les Tableaux En Une Seule Ligne De Commandes
-    ! les tableaux iadvec valent 1 quand la vitesse consideree est upwind et 0 sinon
-
-    !         Advection Selon X
-    !         -----------------
-
-    Iadvec_w(:,:)=Nint(Max(Sign(1.0,Ice_flow_m%U_Column_x(:,:)),0.0))
-    Iadvec_e(:,:)=Nint(Abs(Min(Sign(1.0,Ice_flow_m%U_Column_x(:,:)),0.0)))
-
-    !         Advection Selon Y
-    !         -----------------
-
-    Iadvec_s(:,:)=Nint(Max(Sign(1.0,Ice_flow_m%U_Column_y(:,:)),0.0)) 
-    Iadvec_n(:,:)=Nint(Abs(Min(Sign(1.0,Ice_flow_m%U_Column_y(:,:)),0.0)))
-
-    ! Boucle De Resolution Des Temperatures Colonne Par Colonne
-
-    If (Itracebug.Eq.1) Write(Num_tracebug_m,*)' avant appel  boucle Advec_icetemp'
-
-    Do J=2,Ny_m-1
-       Do I=2,Nx_m-1
-
- !  If (Itracebug.Eq.1) Write(Num_tracebug_m,*)' 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,:),&
-               geom_m,temperature_m,Ice_flow_m,Autre_pr_temp_m,& 
-               step_temp_m,debug_3D_m)
-
-          ! variables de calcul dans la glace
-          Deh22=2.*Dee*Geom_m%Thickness(I,J)**2
-          Duu=step_temp_m%Dtt/2./Dee
-          Rr(Nznzm_m)=-Dzm*Autre_pr_temp_m%Ghf(I,J)/Cm  ! depend de I,J Pour Un Flux Geothermique Variable
-
-!   If (Itracebug.Eq.1) Write(Num_tracebug_m,*)' avant appel  temp_col  i,j', i,j
-          Call temp_col(I,J,T3d_new,Ct(I,J,:),&
-               Mask_flgz_m%Flottant(I,J),Geom_m%Thickness(I,J),Tbmer(I,J),&
-               Temperature_m,Autre_pr_temp_m,step_temp_m,debug_3D_m)
-
-          ! Test Permettant D'Ecrire Les Variables Quand Il Y A Une Erreur
-          ! Ifail1 = 1 : Erreur Dans Tridiag
-          If (Ifail1.Eq.1) Then
-             Write(6,*)'Erreur Dans Tridiag'
-             Write(6,*)'I = ',I,' J = ',J
-             Write(6,*) 'A, B, C, R, U'
-             Do K=1,Nznzm_m
-                Write(6,*)'K = ',K
-                !Write(6,*) Aa(K),Bb(K),Cc(K),Rr(K),Hh(K)
-             End Do
-             Stop
-          End If
-       End Do
-    End Do
-
-    Do J=1,Ny_m
-       Do I=1,Nx_m
-          H1(I,J)=Geom_m%Thickness(I,J)
-          B1(I,J)=Geom_m%Base(I,J)
-          Temperature_m%Temp_melting(I,J,1)=0.
-       End Do
-    End Do
-
-    ! Attribution Finale De La Temperature
-    ! Limitation a 3 deg De Variation Par Dtt Et Pas Plus Froid Que -70
-
-    Do K=1,Nznzm_m
-       Do J=1,Ny_m
-          Do I=1,Nx_m
-             If (Geom_m%Thickness(I,J).Gt.10.) Then
-
-                Tdot(K)=T3d_new(I,J,K)-Temperature_m%Temperature(I,J,K)
-
-                If ((K.Gt.1).And.(Tdot(K).Lt.-3.)) Tdot(K)=-3.
-                If ((K.Gt.1).And.(Tdot(K).Gt.3.)) Tdot(K)=3.
-
-                Temperature_m%Temperature(I,J,K)=Temperature_m%Temperature(I,J,K)+Tdot(K)
-
-                If (Temperature_m%Temperature(I,J,K).Lt.-70.) Then
-                   Temperature_m%Temperature(I,J,K)=-70.
-                Endif
-             Endif
-
-          End Do
-       End Do
-    End Do
-
-
-    Do K=2,Nz_m
-       Do J=1,Ny_m
-          Do I=1,Nx_m
-             If (Temperature_m%Temperature(I,J,K).Gt.Temperature_m%Temp_melting(I,J,K)) Then
-                Temperature_m%Temperature(I,J,K)=Temperature_m%Temp_melting(I,J,K)
-                Autre_pr_temp_m%Ibase(I,J)=2
-             Endif
-          End Do
-       End Do
-    End Do
-
-    !     Temperature Sur Les Bords :
-    !     -------------------------
-  
-    Do K=1,Nz_m
-       Do J=1,Ny_m
-          !        T(Nx_m,J,K)= T(Nx_m-1,J,K)
-          !        T(1,J,K)= T(2,J,K) 
-          Temperature_m%Temperature(Nx_m,J,K)= Autre_pr_temp_m%Ts(Nx_m,J)
-          Temperature_m%Temperature(1,J,K)= Autre_pr_temp_m%Ts(1,J)
-       End Do
-       !
-       Do I=1,Nx_m
-          !        T(I,1,K)=T(I,2,K) 
-          !        T(I,Ny_m,K)=T(I,Ny_m-1,K) 
-          Temperature_m%Temperature(I,1,K)=Autre_pr_temp_m%Ts(I,1)
-          Temperature_m%Temperature(I,Ny_m,K)=Autre_pr_temp_m%Ts(I,Ny_m)
-       End Do
-    End Do
-   
-
-    Do K=1,Nz_m
-       Do J=1,Ny_m
-          Do I=1,Nx_m
-             If (Geom_m%Thickness(I,J).le.1.) then 
-                Temperature_m%Temperature(I,J,K)=Autre_pr_temp_m%Ts(I,J)
-                if ((k.ge.nz).and.(Mask_flgz_m%flottant(i,j))) Temperature_m%Temperature(I,J,K) =Tbmer(I,J)
-             end If
-                
-          End Do
-       End Do
-    End Do
-
-    Deallocate (Aa,Bb,Cc,Rr,Hh,Tdot,Abis,Bbis,Cbis,Rbis,Hbis,Ee,Tbmer,            & 
-                Alpha,Ubar,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)        
-
- If (Itracebug.Eq.1) Write(Num_tracebug_m,*)' fin routine icetemp'
-    return
-  End Subroutine Ice_temp
-
-  !-----------------------------------------------------------------------------------------------------------------------
-
-
-
-
-end module interface_icetempmod
Index: unk/SOURCES/Temperature-routines/interface_prop_th_icetemp.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/interface_prop_th_icetemp.f90	(revision 23)
+++ 	(revision )
@@ -1,106 +1,0 @@
-!> \file interface_prop_th_icetemp.f90
-!! Interfaces of functions and subroutines for thermal propretys calculation
-!<
-
-!> \namespace interface_prop_th_icetemp
-!! Interfaces of functions and subroutines to calculate thermal prop 
-!<
-
-module interface_prop_th_icetemp
-
-  !> \interface Thermal_prop
-  !! Interfaces of functions and subroutines to calculate thermal prop 
-  !!\author ...
-  !!\date 2010
-  !!
-  !<
-
-  Interface  Thermal_property 
-     module procedure  Thermal_prop_icetemp
-  End Interface
-
-  Interface Temp_mer
-     module procedure Temp_mer_icetemp
-  End Interface
-
-contains
-
-  !> SUBROUTINE: Thermal_prop_icetemp
-  !! 
-  !! Calculate the thermal prop of ice
-  !!
-  !! Used modules:
-  !! - Icetemp_declar
-  !>
-
-  Subroutine Thermal_prop_icetemp(T_m,Tpmp_m,H_m,debug_3D_m)
-    Use Icetemp_declar 
-
-    Implicit None
-    !< Arguments
-    Real,Dimension(:,:,:),         Intent(In) ::     T_m         !< Temperature In Ice And Mantle 'O'
-    Real,Dimension(:,:,:),         Intent(Inout) ::  Tpmp_m      !< Pressure Melting Point Temperature In Ice Sheet 'O'
-    Real,Dimension(:,:),           Intent(In)::      H_m         !< Ice Thickness  'O'
-    real,dimension(:,:,:),         Intent(Inout) ::  debug_3D_m  !< Debug array
-    Nx_m=size(Tpmp_m,1)
-    Ny_m=size(Tpmp_m,2)
-    Nz_m=size(Tpmp_m,3)
-
-    If (.Not. Allocated(Ct)) then 
-       Allocate(Cp(Nx_m,Ny_m,Nz_m),&
-            Ct(Nx_m,Ny_m,Nz_m))
-    end if
-
-    If (Itracebug.Eq.1) Write(Num_tracebug,*)' Entree Dans Routine Thermal_prop_icetemp'
-
-    Do K=1,Nz_m  
-       Do J=1,Ny_m
-          Do I=1,Nx_m
-             !  Temperature De Fusion:Tpmp (En Deg C)
-             !  Capacite Calorifique : Cp En J/Kg/K Puis Multiplie Par Ro En J/M3/K
-             !  Conductivite : Ct  En J/M/K/An
-             !  Attention Pour La Conductivite C'Est La Formule De Yin-Chao Yen
-
-             Ct(I,J,K)=6.62E7
-             Cp(I,J,K)=2009. 
-             Tpmp_m(I,J,K)=-0.00087*(K-1)*1./(Nz_m-1)*H_m(I,J)!De=1/(Nz_m-1)
-             Cp(I,J,K)=(2115.3+7.79293*T_m(I,J,K)) 
-             Cp(I,J,K)=Cp(I,J,K)*Ro              ! En J/M3/K
-             Ct(I,J,K)=3.1014e8*Exp(-0.0057*(T_m(I,J,K)+273.15)) !En J/M/K/An
-
-          End Do
-       End Do
-    End Do
-
-  End Subroutine Thermal_prop_icetemp
-  !-----------------------------------------------------------------------------------------------------------
-
-  !> SUBROUTINE: Temp_mer
-  !! Calculate temperature of the sea
-  !! Used modules:
-  !! - Icetemp_declar
-  !>
-  Subroutine Temp_mer_icetemp(H_m,Flot_m,debug_3D_m)
-    Use Icetemp_declar
-    Implicit None 
-    !< Arguments
-    Real,Dimension(:,:),   Intent(In)::H_m                          !< Ice Thickness  'O'
-    Logical,Dimension(:,:),Intent(In)::Flot_m                       !< Vrai Si Flottant (Test D'Archimede) 'O'
-    Real,dimension(:,:,:), Intent(Inout) :: debug_3D_m              !< Debug array
-
-    !Instructions
-    If (Itracebug.Eq.1) Write(Num_tracebug,*)' Entree Dans Routine Temp_mer_icetemp '
-
-    Do J=1,Ny_m
-       Do I=1,Nx_m
-          If (Flot_m(I,J)) &
-               Tbmer(I,J) = Acof1*S0mer + Bcof1 + Ccof1*H_m(I,J)*Ro/Row
-
-       End Do
-    End Do
-  End Subroutine Temp_mer_icetemp
-
-  !-----------------------------------------------------------------------------------------------------------
-
-
-end module interface_prop_th_icetemp
Index: unk/SOURCES/Temperature-routines/interface_temp_col.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/interface_temp_col.f90	(revision 23)
+++ 	(revision )
@@ -1,228 +1,0 @@
-!> \file interface_temp_col.f90
-!! Interfaces of functions and subroutines for vertical temperature calculation 
-!<
-
-!> \namespace interface_temp_col
-!!  Interface of functions and subroutines for vertical temperature calculation 
-!<
-module interface_temp_col
-  !> \interface temp_col
-  !! Interface of functions and subroutines for vertical temperature calculation 
-  !!\author ...
-  !!\date 2010
-  !<
-  Interface temp_col
-     module procedure Temp_column
-  End Interface
-
-
-contains
-
-  !> SUBROUTINE: Temp_column 
-  !! Subroutine for vertical temperature calculation
-  !! Used modules:
-  !! - use Icetemp_declar
-  !! - use Tridiagmod
-  !! - use temperature_type
-  !! - use autre_pr_temp_type
-  !! - use step_temp_type
-  !>
-
-  Subroutine  Temp_column(Ii,Jj,Newtempcol_m,Ct_m,Flot_m,H_m,Tbmer_m,&
-       Temperature_m,Autre_pr_temp_m,step_temp_m,debug_3D_m)
-
-    Use Icetemp_declar
-    Use Tridiagmod
-    use temperature_type
-    use autre_pr_temp_type
-    use step_temp_type
-
-
-    Implicit None
-    integer :: Ii,Jj
-    Real,Dimension(:,:,:),     Intent(Out)   :: Newtempcol_m        !< Tableau 1d Vert. Des Temperatures En Sortie
-    Real,Dimension(:),         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
-    type (Temperature),        Intent(Inout) :: temperature_m       !< Variable temperature
-    type (Autre_pr_temp),      Intent(Inout) :: Autre_pr_temp_m     !< Variable 
-    type (step_temp),          Intent(In)    :: step_temp_m         !< Variable step pour temperature
-    real,dimension(:,:,:),     Intent(Inout) :: debug_3D_m
-
-    If (H_m.Gt.10.) Then
-
-       ! Conditions Aux Limites a La Base De La Glace
-       !________________________________________________
-
-       ! Base Froide 
-       If ( ((Autre_pr_temp_m%Ibase(Ii,Jj).Eq.1).Or.(Autre_pr_temp_m%Ibase(Ii,Jj).Eq.4) &
-            .Or. ((Autre_pr_temp_m%Ibase(Ii,Jj).Eq.5).And.(Temperature_m%Temperature(Ii,Jj,Nz_m).Lt.Temperature_m%Temp_melting(Ii,Jj,Nz_m)))) &
-            .And..Not.Flot_m) Then
-
-          Autre_pr_temp_m%Ibase(Ii,Jj)=1
-          Bb(Nz_m)=1.
-
-          If (Ncond.Eq.1) Then ! Avec Socle
-             Dzi=H_m*Dee*Cm/Ct_m(Nz_m)
-             Aa(Nz_m)=-Dzm/(Dzm+Dzi)
-             Cc(Nz_m)=-Dzi/(Dzm+Dzi)
-             Rr(Nz_m)=H_m*Dee*Autre_pr_temp_m%Phid(Ii,Jj)*Dzm/Ct_m(Nz_m)/(Dzm+Dzi)
-          Else                 ! Sans Socle
-             Aa(Nz_m)=-1.
-             Cc(Nz_m)=0.
-             Rr(Nz_m)=-(Autre_pr_temp_m%Ghf(Ii,Jj)-Autre_pr_temp_m%Phid(Ii,Jj))/Ct_m(Nz_m)*H_m*Dee
-          Endif
-
-       Else
-
-          ! Base Temperee Ou Shelf
-          ! ----------------------               
-          If (Autre_pr_temp_m%Ibase(Ii,Jj).Eq.5)  Autre_pr_temp_m%Ibase(Ii,Jj)=2
-          Autre_pr_temp_m%Ibase(Ii,Jj)=Max(Autre_pr_temp_m%Ibase(Ii,Jj),2)
-          Aa(Nz_m)=0.
-          Bb(Nz_m)=1.
-          Cc(Nz_m)=0.
-
-          If (.Not.Flot_m) Then
-             Rr(Nz_m)=Temperature_m%Temp_melting(Ii,Jj,Nz_m)
-          Else
-             Rr(Nz_m)=Tbmer_m
-          End If
-
-       Endif
-
-
-       If (Ncond.Eq.1) Then ! Avec Socle
-          Do K=Nz_m+1,NzNzm_m-1
-             Rr(K)=Temperature_m%Temperature(Ii,Jj,K)
-          End Do
-          Call Tridiag (Aa,Bb,Cc,Rr,Hh,NzNzm_m,Ifail1)
-
-       Else                 ! Sans Socle
-
-          Call Tridiag (Aa,Bb,Cc,Rr,Hh,Nz_m,Ifail1)
-       Endif
-
-
-       ! Temperature Dans Le Socle Si Ncond=0, Lineaire Avec Le Gradient Ghf
-       If (Ncond.Eq.0) Then
-          Do K=Nz_m+1,NzNzm_m
-             Hh(K)=Hh(Nz_m)-Dzm*(K-Nz_m)*Autre_pr_temp_m%Ghf(Ii,Jj)/Cm
-          End Do
-       Endif
-
-       Do K=1,NzNzm_m
-          Newtempcol_m(Ii,Jj,K)=Hh(K)
-       End Do
-       Autre_pr_temp_m%Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz_m)-Temperature_m%Temperature(Ii,Jj,Nz_m))/step_temp_m%Dtt
-
-       !      ------------------------------------------------- Glace Tres Fine (H<10m Ou H=0)
-    Else If (H_m.Le.10.) Then
-
-       !       Pour Eviter Des Problemes Lorsque H Est Inferieur A 10 M
-       !       Profil Lineaire, Pente=Flux Geothermique
-       !       Le Cas Sans Glace Est Traite De La Meme Facon  (Dou=0)
-
-       ! ........................................ Avec Calcul Dans Le Socle
-       If (Ncond.Eq.1.And..Not.Flot_m) Then   
-          If (H_m.Gt.0.) Then
-             ! Gradient Dans Le Socle
-             Dou=(Temperature_m%Temperature(Ii,Jj,Nz_m+1)-Temperature_m%Temperature(Ii,Jj,Nz_m))/Dzm*Cm     
-             Dou=Dou/Ct_m(Nz_m)*Dee*H_m
-
-             Tss=Min(0.,Autre_pr_temp_m%Ts(Ii,Jj))
-             Do K=1,Nz_m
-                Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
-             End Do
-
-          Else
-             ! Pas De Glace
-             Tss=Autre_pr_temp_m%Ts(Ii,Jj)
-             Do K=1,Nz_m
-                Newtempcol_m(Ii,Jj,K)=Tss
-             End Do
-          End If
-
-          ! Calcul Dans Le Socle Meme S'Il N'Y A Pas De Glace
-          Rr(Nz_m)=Newtempcol_m(Ii,Jj,Nz_m)
-          Aa(Nz_m)=0
-          Bb(Nz_m)=1
-          Cc(Nz_m)=0
-
-          Do K=Nz_m+1,NzNzm_m-1
-             Rr(K)=Temperature_m%Temperature(Ii,Jj,K)
-          End Do
-
-          ! Creation De Nouveaux Tableaux Juste Pour Tridiag      
-          Do K=1,Nzm_m+1
-             Abis(K)=Aa(Nz_m-1+K)
-             Bbis(K)=Bb(Nz_m-1+K)
-             Cbis(K)=Cc(Nz_m-1+K)
-             Rbis(K)=Rr(Nz_m-1+K)
-          End Do
-
-          Call Tridiag (Abis,Bbis,Cbis,Rbis,Hbis,Nzm_m+1,Ifail1)
-
-          Do K=1,Nzm_m+1
-             Hh(Nz_m-1+K)=Hbis(K)
-          End Do
-
-          ! ........................................ Sans Calcul Dans Le Socle
-       Else
-          ! Calotte Posee                    
-          If ((H_m.Gt.0.).And..Not.Flot_m) Then
-             Dou=-Autre_pr_temp_m%Ghf(Ii,Jj)/Ct_m(Nz_m)*Dee*H_m
-             Tss=Min(0.,Autre_pr_temp_m%Ts(Ii,Jj))
-             Do K=1,Nz_m
-                Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
-             End Do
-
-             !  Shelf
-          Else If ((H_m.Gt.0.).And.Flot_m) Then
-             Tss=Min(0.,Autre_pr_temp_m%Ts(Ii,Jj))
-             Dou=(Tbmer_m-Tss)*Dee
-             Do K=1,Nz_m
-                Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.)
-             End Do
-
-             ! Pas De Glace
-          Else
-             Tss=Autre_pr_temp_m%Ts(Ii,Jj)      
-             Do K=1,Nz_m
-                Newtempcol_m(Ii,Jj,K)=Tss
-             End Do
-          End If
-
-          ! Temperature Dans Le Socle, Lineaire Avec Le Gradient Ghf
-          If (.Not.Flot_m) Then
-             Do K=Nz_m+1,NzNzm_m
-                Hh(K)=Newtempcol_m(Ii,Jj,Nz_m)-Dzm*(K-Nz_m)*Autre_pr_temp_m%Ghf(Ii,Jj)/Cm
-             End Do
-          Else
-             Do K=Nz_m+1,NzNzm_m
-                Hh(K)=Tbmer_m-Dzm*(K-Nz_m)*Autre_pr_temp_m%Ghf(Ii,Jj)/Cm
-             End Do
-          End If
-
-       Endif
-       ! Fin Du Test Socle
-
-
-       Do K=Nz_m+1,NzNzm_m
-          Newtempcol_m(Ii,Jj,K)=Hh(K)
-       End Do
-
-       Autre_pr_temp_m%Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz_m)-Temperature_m%Temperature(Ii,Jj,Nz_m))/step_temp_m%Dtt
-       Autre_pr_temp_m%Bmelt(Ii,Jj)=0.
-       Autre_pr_temp_m%Ibase(Ii,Jj)=5
-       Autre_pr_temp_m%Phid(Ii,Jj)=0.
-
-
-    End If
-
-  End Subroutine Temp_column
-
-  !--------------------------------------------------------------------------------------------------------------------
-
-end module interface_temp_col
Index: /trunk/SOURCES/Temperature-routines/prop_th_icetemp.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/prop_th_icetemp.f90	(revision 24)
+++ /trunk/SOURCES/Temperature-routines/prop_th_icetemp.f90	(revision 24)
@@ -0,0 +1,62 @@
+!> \file prop_th_icetemp.f90
+!! Subroutines for thermal propreties calculation
+!<
+
+!> \namespace prop_th_icetemp
+!! Subroutines to calculate thermal prop 
+
+!> SUBROUTINE: Thermal_prop_icetemp
+!! 
+!! Calculate the thermal prop of ice
+!!
+!! Used modules:
+!! - Icetemp_declar
+!>
+
+Subroutine Thermal_prop_icetemp
+  Use Icetemp_declar 
+
+  Implicit None
+
+  If (Itracebug.Eq.1) Write(Num_tracebug,*)' Entree Dans Routine Thermal_prop_icetemp'
+
+  Do K=1,Nz  
+     Do J=1,Ny
+        Do I=1,Nx
+           !  Temperature De Fusion:Tpmp (En Deg C)
+           !  Capacite Calorifique : Cp En J/Kg/K Puis Multiplie Par Ro En J/M3/K
+           !  Conductivite : Ct  En J/M/K/An
+           !  Attention Pour La Conductivite C'Est La Formule De Yin-Chao Yen
+
+           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)) 
+           Cp(I,J,K)=Cp(I,J,K)*Ro              ! En J/M3/K
+           Ct(I,J,K)=3.1014e8*Exp(-0.0057*(T(I,J,K)+273.15)) !En J/M/K/An
+        End Do
+     End Do
+  End Do
+
+End Subroutine Thermal_prop_icetemp
+
+!-----------------------------------------------------------------------------------------------------------
+
+!> SUBROUTINE: Temp_mer
+!! Calculate temperature of the sea
+!! Used modules:
+!! - Icetemp_declar
+!>
+Subroutine Temp_mer
+  Use Icetemp_declar
+  Implicit None 
+
+  If (Itracebug.Eq.1) Write(Num_tracebug,*)' Entree Dans Routine Temp_mer_icetemp '
+
+  Do J=1,Ny
+     Do I=1,Nx
+        If (Flot(I,J)) &
+             Tbmer(I,J) = Acof1*S0mer + Bcof1 + Ccof1*H(I,J)*Ro/Row
+     End Do
+  End Do
+End Subroutine Temp_mer
Index: /trunk/SOURCES/Temperature-routines/temp_col.f90
===================================================================
--- /trunk/SOURCES/Temperature-routines/temp_col.f90	(revision 24)
+++ /trunk/SOURCES/Temperature-routines/temp_col.f90	(revision 24)
@@ -0,0 +1,202 @@
+!> \file temp_col.f90
+!! Subroutines for vertical temperature calculation 
+!<
+
+!> \namespace temp_col
+!!  Subroutines for vertical temperature calculation 
+!<
+
+!> SUBROUTINE: Temp_column 
+!! Subroutine for vertical temperature calculation
+!! Used modules:
+!! - use Icetemp_declar
+!! - use Tridiagmod
+!>
+
+Subroutine  Temp_col(Ii,Jj,Newtempcol_m,Ct_m,Flot_m,H_m,Tbmer_m)
+
+  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
+
+     ! Conditions Aux Limites a La Base De La Glace
+     !________________________________________________
+
+     ! 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
+        Bb(Nz)=1.
+
+        If (Ncond.Eq.1) Then ! Avec Socle
+           Dzi=H_m*Dee*Cm/Ct_m(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)
+        Else                 ! Sans Socle
+           Aa(Nz)=-1.
+           Cc(Nz)=0.
+           Rr(Nz)=-(Ghf(Ii,Jj)-Phid(Ii,Jj))/Ct_m(Nz)*H_m*Dee
+        Endif
+
+     Else
+
+        ! Base Temperee Ou Shelf
+        ! ----------------------               
+        If (Ibase(Ii,Jj).Eq.5)  Ibase(Ii,Jj)=2
+        Ibase(Ii,Jj)=Max(Ibase(Ii,Jj),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
+        End If
+
+     Endif
+
+
+     If (Ncond.Eq.1) Then ! Avec Socle
+        Do K=Nz+1,nz+nzm-1
+           Rr(K)=T(Ii,Jj,K)
+        End Do
+        Call Tridiag (Aa,Bb,Cc,Rr,Hh,nz+nzm,Ifail1)
+
+     Else                 ! Sans Socle
+
+        Call Tridiag (Aa,Bb,Cc,Rr,Hh,Nz,Ifail1)
+     Endif
+
+
+     ! Temperature Dans Le Socle Si Ncond=0, Lineaire Avec Le Gradient Ghf
+     If (Ncond.Eq.0) Then
+        Do K=Nz+1,nz+nzm
+           Hh(K)=Hh(Nz)-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
+        End Do
+     Endif
+
+     Do K=1,nz+nzm
+        Newtempcol_m(Ii,Jj,K)=Hh(K)
+     End Do
+     Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz)-T(Ii,Jj,Nz))/Dtt
+
+     !      ------------------------------------------------- Glace Tres Fine (H<10m Ou H=0)
+  Else If (H_m.Le.10.) Then
+
+     !       Pour Eviter Des Problemes Lorsque H Est Inferieur A 10 M
+     !       Profil Lineaire, Pente=Flux Geothermique
+     !       Le Cas Sans Glace Est Traite De La Meme Facon  (Dou=0)
+
+     ! ........................................ Avec Calcul Dans Le Socle
+     If (Ncond.Eq.1.And..Not.Flot_m) Then   
+        If (H_m.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.)
+           End Do
+
+        Else
+           ! Pas De Glace
+           Tss=Ts(Ii,Jj)
+           Do K=1,Nz
+              Newtempcol_m(Ii,Jj,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)
+        Aa(Nz)=0
+        Bb(Nz)=1
+        Cc(Nz)=0
+
+        Do K=Nz+1,nz+nzm-1
+           Rr(K)=T(Ii,Jj,K)
+        End Do
+
+        ! Creation De Nouveaux Tableaux Juste Pour Tridiag      
+        Do K=1,Nzm+1
+           Abis(K)=Aa(Nz-1+K)
+           Bbis(K)=Bb(Nz-1+K)
+           Cbis(K)=Cc(Nz-1+K)
+           Rbis(K)=Rr(Nz-1+K)
+        End Do
+
+        Call Tridiag (Abis,Bbis,Cbis,Rbis,Hbis,Nzm+1,Ifail1)
+
+        Do K=1,Nzm+1
+           Hh(Nz-1+K)=Hbis(K)
+        End Do
+
+        ! ........................................ Sans Calcul Dans Le Socle
+     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.)
+           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.)
+           End Do
+
+           ! Pas De Glace
+        Else
+           Tss=Ts(Ii,Jj)      
+           Do K=1,Nz
+              Newtempcol_m(Ii,Jj,K)=Tss
+           End Do
+        End If
+
+        ! Temperature Dans Le Socle, Lineaire Avec Le Gradient Ghf
+        If (.Not.Flot_m) Then
+           Do K=Nz+1,nz+nzm
+              Hh(K)=Newtempcol_m(Ii,Jj,Nz)-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
+           End Do
+        Else
+           Do K=Nz+1,nz+nzm
+              Hh(K)=Tbmer_m-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm
+           End Do
+        End If
+
+     Endif
+     ! Fin Du Test Socle
+
+
+     Do K=Nz+1,nz+nzm
+        Newtempcol_m(Ii,Jj,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.
+
+
+  End If
+
+End Subroutine Temp_col
Index: /trunk/SOURCES/initial-0.3.f90
===================================================================
--- /trunk/SOURCES/initial-0.3.f90	(revision 23)
+++ /trunk/SOURCES/initial-0.3.f90	(revision 24)
@@ -12,6 +12,5 @@
 !! @note    - use module3D_phy
 !! @note    - use module_choix
-!! @note    - use interface_icetempmod
-!! @note    - use interface_prop_th_icetemp
+!! @note    - use icetempmod
 !! @note    - use diagno_mod 
 !! @note    - use resolmeca_SIA_L1
@@ -25,6 +24,5 @@
   use module3D_phy
   use module_choix
-  use interface_icetempmod
-  use interface_prop_th_icetemp
+  use icetempmod
   use diagno_mod 
   use resolmeca_SIA_L1
@@ -88,5 +86,5 @@
   if (geoplace(1:5).ne.'mism3') then
      call init_icetemp(Num_rep_42)        ! initialise le calcul de temperature
-     call thermal_prop_icetemp(T,TPMP,H,debug_3D)  ! appel aux proprietes thermiques
+     call thermal_prop_icetemp            ! appel aux proprietes thermiques
      call init_bmelt                      ! Appel aux tableaux d'initialisation des ice shelves
 
@@ -137,5 +135,4 @@
   call init_calving
 
-  call allocate_types  ! cette routine alloue tous les types derive 
 !
 !-------------------------------------------------------------------------------------
Index: /trunk/SOURCES/main3D-0.4-40km.f90
===================================================================
--- /trunk/SOURCES/main3D-0.4-40km.f90	(revision 23)
+++ /trunk/SOURCES/main3D-0.4-40km.f90	(revision 24)
@@ -55,5 +55,5 @@
 !! @note use module_choix 
 !! @note use flottab_mod
-!! @note use interface_icetempmod
+!! @note use icetempmod
 !! @note use sorties_ncdf_grisli
 !! @note use diagno_mod 
@@ -85,5 +85,5 @@
 !! @note use module_choix 
 !! @note use flottab_mod
-!! @note use interface_icetempmod
+!! @note use icetempmod
 !! @note use sorties_ncdf_grisli
 !! @note use diagno_mod 
@@ -107,5 +107,5 @@
   !  de declaration des packages
   use flottab_mod
-  use interface_icetempmod
+  use icetempmod
   use sorties_ncdf_grisli
   use diagno_mod 
@@ -147,5 +147,5 @@
   !  de declaration des packages
   use flottab_mod
-  use interface_icetempmod
+  use icetempmod
   use sorties_ncdf_grisli
   use util_recovery
Index: /trunk/SOURCES/steps_time_loop.f90
===================================================================
--- /trunk/SOURCES/steps_time_loop.f90	(revision 23)
+++ /trunk/SOURCES/steps_time_loop.f90	(revision 24)
@@ -16,11 +16,4 @@
   use module3d_phy
 
-  use geom_type
-  use temperature_type
-  use ice_flow_type
-  use mask_flgz_type
-  use deformation_type
-  use autre_pr_temp_type
-
   use module_choix ! module de choix du type de run
   !  module_choix donne acces a tous les modules
@@ -29,5 +22,5 @@
   use sorties_ncdf_grisli
   use flottab_mod
-  use interface_icetempmod
+  use icetempmod
   use diagno_mod 
   !  use track_debug 
@@ -199,5 +192,5 @@
   !  module_choix donne acces a tous les modules
   !  de declaration des packages
-  use interface_icetempmod
+  use icetempmod
   use sorties_ncdf_grisli
   use flottab_mod
@@ -205,12 +198,4 @@
   use resolmeca_SIA_L1
 !  use track_debug 
-
-  use geom_type
-  use temperature_type
-  use ice_flow_type
-  use mask_flgz_type
-  use deformation_type
-  use autre_pr_temp_type
-  use step_temp_type
 
 
@@ -303,29 +288,6 @@
 
      calc_temp:     if ((nt.gt.2).and.(geoplace(1:5).ne.'mism3')) then
-        ! Hassine ca va etre modifier bien sur apres: l appel sera beaucoup moin long que ca 
-        call Update_geom(Geom_g)               
-        call Update_Temperature(Temperature_g)
-        call Update_Ice_flow(Ice_flow_g)
-        call Update_Mask_flgz(Mask_flgz_g)
-        call Update_Deformation(Deformation_g)
-        call Update_autre_pr_temp(Autre_pr_temp_g)
-        call Update_step_temp(step_temp_g)
-
         if (itracebug.eq.1)  call tracebug('avant appel icetemp')
-
-        call icetemp(geom_g,temperature_g, &
-             Ice_flow_g,Mask_flgz_g,&
-             deformation_g,Autre_pr_temp_g,&
-             step_temp_g,num_tracebug,debug_3D)
-
-        ! update arrays T, Tpmp
-        T= Temperature_g%Temperature 
-        Tpmp= Temperature_g%Temp_melting
-        Tbdot= Autre_pr_temp_g%Tbdot
-        Bmelt= Autre_pr_temp_g%BMELT
-        Ibase=Autre_pr_temp_g%IBASE
-        Phid=Autre_pr_temp_g%Phid
-
-        !      write(6,*) time,'call icetemp'
+        call icetemp
         ! subroutines pour le calcul de la fusion basale
         call bmeltshelf
