!> \file Qprod_icetemp_3D_mod.f90
!! Implementation of functions and subroutines for heat production 
!<

!
!> SUBROUTINE: Qprod
!! 
!! Subroutine for heat production
!!
!! 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)

  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 global
  type (Temperature),   intent(in)    :: T_m              !  variable temperature global
  type (Ice_flow),      intent(in)    :: Ice_flow_m       !  ice_flow variables global 
  type (Mask_flgz),    intent(in)    :: mask_flot_m      !  
  type (Deformation),   intent(in)    :: deform_m
  type (Autre_pr_temp), intent(inout) :: therm_var_m


  !< 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 Posée
                 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*E(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*E(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.*(E(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*E(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*E(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.*(E(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*(E(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)*E(K)*Qdef3(I,J,K)/2./De    ! 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)/De               ! 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)*E(K)*Qdef3(I,J,K)/2./De    ! 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)/De               ! 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


!///////////////////////////////////////////////////////////////////////
!///////////////////////////////////////////////////////////////////////