source: trunk/SOURCES/Temperature-routines/Old/Qprod_icetemp_3D_mod.f90 @ 154

Last change on this file since 154 was 4, checked in by dumas, 10 years ago

initial import GRISLI trunk
File size: 37.8 KB
Line 
1!> \file Qprod_icetemp_3D_mod.f90
2!! Implementation of functions and subroutines for heat production
3!<
4
5!
6!> SUBROUTINE: Qprod
7!!
8!! Subroutine for heat production
9!!
10!! Used modules:
11!!    - use Icetemp_declar
12!!    - use geom_type
13!!    - use temperature_type
14!!    - use Ice_flow_type
15!!    - use Mask_flgz_type
16!!    - use Deformation_type
17!!    - use Autre_pr_temp_type
18!>
19
20
21subroutine Qprod (Iq1,geom_m,T_m,Ice_flow_m,mask_flot_m,deform_m,therm_var_m)
22
23  use Icetemp_declar
24  use geom_type
25  use temperature_type
26  use Ice_flow_type
27  use Mask_flgz_type
28  use Deformation_type
29  use Autre_pr_temp_type
30 
31  implicit none
32
33  !<Arguments
34
35  integer,              intent(In)    :: Iq1              !< Choice Of Method For Heat Production
36  type (Geom),          intent(in)    :: geom_m           !  variable geometrie global
37  type (Temperature),   intent(in)    :: T_m              !  variable temperature global
38  type (Ice_flow),      intent(in)    :: Ice_flow_m       !  ice_flow variables global
39  type (Mask_flgz),    intent(in)    :: mask_flot_m      ! 
40  type (Deformation),   intent(in)    :: deform_m
41  type (Autre_pr_temp), intent(inout) :: therm_var_m
42
43
44  !< Local variables needed by the various case of Iq1
45
46  real, dimension(Nx_m,Ny_m)      :: Pente2_maj    !< Pente Au Carre Sur Les Noeuds Majeurs
47  real, dimension(Nx_m,Ny_m)      :: Pente_maj     !< Pente Au Carre Sur Les Noeuds Majeurs
48  real, dimension(Nx_m,Ny_m)      :: Vit_maj       !< Ubar Moyennee Sur Les Noeuds Majeurs
49  real, dimension(Nx_m,Ny_m)      :: Uslid_maj     !< Uslid Moyennee Sur Les Noeuds Majeurs
50  real, dimension(Nx_m,Ny_m)      :: Vit_stag2     !< Ubar Moyennee Sur Les Noeuds Stag
51  real, dimension(Nx_m,Ny_m,Nz_m) :: Vit_stag3     !< Ubar Moyennee Sur Les Noeuds Stag
52  real, dimension(Nx_m,Ny_m)      :: Hmxy          !< Epaisseur Moyenne Sur Les Noeuds Stag
53  real, dimension(Nx_m,Ny_m)      :: Qslid         !< Chaleur Produite Par Glissement Sur Les Noeuds Stag
54  real, dimension(Nx_m,Ny_m)      :: Qdef2         !< Chaleur Produite Par Deformation Sur Les Noeuds Stag
55  real, dimension(Nx_m,Ny_m,Nz_m) :: Qdef3         !< Chaleur Produite Par Deformation Sur Les Noeuds Stag
56  real, dimension(Nx_m,Ny_m)      :: Pente_stag    !< Pente Au Carre Sur Les Noeuds Stag
57  real, dimension(Nx_m,Ny_m)      :: Uslid_stag    !< Uslid Moyennee Sur Les Noeuds Stag
58  real, dimension(Nx_m,Ny_m)      :: Tob_stag      !< Contrainte Basale
59  real, dimension(Nx_m,Ny_m)      :: Tox           !< Contraintes Sur Maille Mx
60  real, dimension(Nx_m,Ny_m)      :: Toy           !< Contraintes Sur Maille Mx
61
62  select case (Iq1)
63
64  case (1)                ! iq1 = 1      Q_prod_demi, Ancienne Methode Sur Les Demi-Mailles
65!------------------------------------------------------------------------------------------
66
67     ! Calcul Avec L'Ancienne Methode Sur Les Demi-Mailles
68     ! Nom Du Sub Routine Avant Et Var Utilisee 
69     ! 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,&
70     !      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)
71
72     do K=2,Nz_m
73        do J=2,Ny_m-1
74           do I=2,Nx_m-1
75
76              ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
77              ! Pour Les Ice-Streams Et Ice-Shelves
78              ! Les Divers Eps Sont Calcules Dans Strain-Rate Pour L'Ensemble De La Grille
79
80              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.  &
81                   (mask_flot_m%Flgz_my(I,J).or.mask_flot_m%Flgz_my(I,J+1))) then
82
83                 Chal2_x(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_xx(I,J)**2
84                 Chal2_y(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_yy(I,J)**2
85                 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
86
87                 !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
88                 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))
89                 Chal2_xy(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
90
91              else   ! Glace Posée
92                 Chal2_x(I,J,K)=0.
93                 Chal2_y(I,J,K)=0.
94                 Chal2_z(I,J,K)=0.
95                 Chal2_xy(I,J,K)=0.
96              endif
97           end do
98        end do
99     end do
100
101     !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
102     do L=1,size(Deform_m%Btt_Deform,4) !N1poly,N2poly
103        do J=2,Ny_m
104           do I=2,Nx_m
105              !          Ffx A 3 Dimensions !
106              Ffx(I,J,L)=Ice_flow_m%Difus_x(I,J,L)*Geom_m%Slop_x(I,J)*Geom_m%Slop_x(I,J)
107              Ffy(I,J,L)=Ice_flow_m%Difus_y(I,J,L)*Geom_m%Slop_y(I,J)*Geom_m%Slop_y(I,J)
108           end do
109        end do
110     end do
111
112     do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
113        do K=2,Nz_m
114           do J=2,Ny_m
115              do I=2,Nx_m
116                 if ((.not.mask_flot_m%Flot_mx(I,J)).and.(.not.mask_flot_m%Gz_mx(I,J))) then
117                    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) !&
118                    !                      *Ro*G*E(K)**(Glen(L)+1)/Cp(I,J,K)
119
120                 else if (mask_flot_m%Gz_mx(I,J)) then ! Ice Streams
121                    Chalx(I,J,K,L)=0.
122                 else                     ! Ice Shelves 
123                    Chalx(I,J,K,L)=0.
124                 endif
125                 if ((.not.mask_flot_m%Flot_my(I,J)).and.(.not.mask_flot_m%Gz_my(I,J))) then
126                    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) !&
127                    !                      *Ro*G*E(K)**(Glen(L)+1)/Cp(I,J,K)
128                 else if (mask_flot_m%Gz_my(I,J)) then  ! Ice Streams
129                    Chaly(I,J,K,L)=0.
130
131                 else                      ! Ice Shelves
132                    Chaly(I,J,K,L)=0.
133                 endif
134
135              end do
136           end do
137        end do
138     end do
139
140     !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
141     !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
142     !         Les Productions X Et Y
143     !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy)
144
145     do K=2,Nz_m
146        do J=1,Ny_m-1
147           do I=1,Nx_m-1
148
149              ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
150              Chaldef_maj(I,J,K)= 0.
151
152              ! Chalk_2 Pour Ice Shelves Et Ice Streams
153              Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
154                   Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
155
156              ! Chalk_1 Pour La Partie Posée
157              do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
158
159                 ! On Somme La Chaleur Due Aux Diverses Lois De Déformation Et Celles En X Et Y
160                 Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))+ &
161                      (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
162
163                 Chalk_1=Ro*G*Chalk_1/4.*(E(K)**(Deform_m%Glen_exp(L)+1))/Cp(I,J,K)
164                 Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
165              enddo
166
167              ! Pour Shelves Et Streams,  On Ajoute Chalk_2
168              Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
169           end do
170        end do
171     end do
172     ! Chaleur Produite A La Base Par Le Glissement
173     do J=2,Ny_m
174        do I=2,Nx_m
175
176           if (mask_flot_m%Gz_mx(I,J)) then
177              Chalglissx(I,J)= abs(Ice_flow_m%U_Column_x(I,J)*Deform_m%Cisail_Basal_mx(I,J))
178           else
179              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)
180           endif
181
182           if (mask_flot_m%Gz_my(I,J)) then
183              Chalglissy(I,J)= abs(Ice_flow_m%U_Column_y(I,J)*Deform_m%Cisail_Basal_my(I,J))
184           else
185              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)
186           endif
187
188        end do
189     end do
190
191     !  Boundary Condition Ice-Rock Interface
192     K=Nz_m
193
194     do J=2,Ny_m-1
195        do I=2,Nx_m-1
196
197           ! L'Ancien Chalbed Correspond A Chaldef_maj(I,J,Nz_m)  : Pas La Peine De Recalculer
198           !----------------------------------------------------
199           !      Chalbed=0.
200           !      Do L=1,size(Deform_m%Btt,4)!N1poly,N2poly
201           !         Chalbed_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L)) &
202           !              + (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
203           !         Chalbed_1=Ro*G*Chalbed_1/4.*H(I,J)
204           !         Chalbed= Chalbed_1 +Chalbed
205           !      Enddo
206
207           !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
208           !        Dans La Derniere 1/2 Maille Et Par Le Glissement
209           !        Attention Phid Est >0 Et Ghf Est <0
210           if (.not.mask_flot_m%Flottant(I,J)) then
211
212              !           Phid Avec Fonte Sous Les Streams
213              !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
214              !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
215              !                        Chalbed*Fracq
216
217
218              ! Moyenne Phid Sur 4 Points : Attention 0.5 Pour Moyenne Gauche - Droite
219              ! Mais Les Chaleurs X,Y S'Ajoutent
220
221              !              Phid(I,J)=0.5*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
222              !                        Chalglissy(I,J)+Chalglissy(I,J+1))
223
224              ! Moyenne Phid Sur 12 Points. On Moyenne Chaque Chaleur Et On Somme X Et Y
225
226              therm_var_m%Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
227                   Chalglissy(I,J)+Chalglissy(I,J+1)) 
228              therm_var_m%Phid(I,J)=therm_var_m%Phid(I,J)+0.125*( &
229                   ((Chalglissx(I,J-1)+Chalglissx(I+1,J-1))+ &
230                   (Chalglissx(I,J+1)+Chalglissx(I+1,J+1))) + &
231                   ((Chalglissy(I-1,J)+Chalglissy(I-1,J+1)) +  &
232                   (Chalglissy(I+1,J)+Chalglissy(I+1,J+1))))
233
234              ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
235              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)
236
237              ! Pour Sorties Heino
238              Chalgliss_maj(I,J)=therm_var_m%Phid(I,J)
239
240              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)
241
242           else
243              therm_var_m%Phid(I,J)=0.
244           endif
245        end do
246     end do
247
248
249  case (2)             ! Q_prod_centre : Calcul Avec La Somme Des Carres
250     !-------------------------------------------------------------------------------------
251     ! Calcul Avec La Somme Des Carres
252     ! Nom Du Sub Routine Avant Et Var Utilisee 
253     ! 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,&
254     !     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)
255
256     do K=2,Nz_m 
257        do J=2,Ny_m-1
258           do I=2,Nx_m-1
259
260
261              ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
262              ! Pour Les Ice-Streams Et Ice-Shelves
263              ! Les Divers Eps Sont Calculés Dans Strain-Rate Pour L'Ensemble De La Grille
264
265              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.  &
266                   (mask_flot_m%Flgz_my(I,J).or.mask_flot_m%Flgz_my(I,J+1))) then
267
268                 Chal2_x(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_xx(I,J)**2
269                 Chal2_y(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_yy(I,J)**2
270                 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
271
272                 !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
273                 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))
274                 Chal2_xy(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
275
276              else   ! Glace Posee
277                 Chal2_x(I,J,K)=0.
278                 Chal2_y(I,J,K)=0.
279                 Chal2_z(I,J,K)=0.
280                 Chal2_xy(I,J,K)=0.
281              endif
282           end do
283        end do
284     end do
285
286     !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
287     do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
288        do J=2,Ny_m
289           do I=2,Nx_m
290
291              !          Ffx A 3 Dimensions !
292              Ffx(I,J,L)=Ice_flow_m%Difus_x(I,J,L)*Geom_m%Slop_x(I,J)*Geom_m%Slop_x(I,J)
293              Ffy(I,J,L)=Ice_flow_m%Difus_y(I,J,L)*Geom_m%Slop_y(I,J)*Geom_m%Slop_y(I,J)
294           end do
295        end do
296     end do
297
298
299     do L=1,size(Deform_m%Btt_Deform,4) !N1poly,N2poly
300        do K=2,Nz_m
301           do J=2,Ny_m
302              do I=2,Nx_m
303                 if ((.not.mask_flot_m%Flot_mx(I,J)).and.(.not.mask_flot_m%Gz_mx(I,J))) then
304                    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) !&
305                    !                      *Ro*G*E(K)**(Glen(L)+1)/Cp(I,J,K)
306
307                 else if (mask_flot_m%Gz_mx(I,J)) then ! Ice Streams
308                    Chalx(I,J,K,L)=0.
309
310                 else                     ! Ice Shelves 
311                    Chalx(I,J,K,L)=0.
312
313                 endif
314
315                 if ((.not.mask_flot_m%Flot_my(I,J)).and.(.not.mask_flot_m%Gz_my(I,J))) then
316                    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) !&
317                    !                      *Ro*G*E(K)**(Glen(L)+1)/Cp(I,J,K)
318
319                 else if (mask_flot_m%Gz_my(I,J)) then  ! Ice Streams
320                    Chaly(I,J,K,L)=0.
321
322                 else                      ! Ice Shelves
323                    Chaly(I,J,K,L)=0.
324                 endif
325
326              end do
327           end do
328        end do
329     end do
330
331     !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
332     !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
333     !         Les Productions X Et Y
334     !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy)
335
336     do K=2,Nz_m
337        do J=1,Ny_m-1
338           do I=1,Nx_m-1
339
340              ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
341              Chaldef_maj(I,J,K)= 0.
342
343              ! Chalk_2 Pour Ice Shelves Et Ice Streams
344              Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
345                   Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
346
347              ! Chalk_1 Pour La Partie Posée
348              do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
349
350                 ! On Somme La Chaleur Due Aux Diverses Lois De Déformation Et Celles En X Et Y
351                 ! La Somme Se Fait Par (Cx2+Cy2)** 0.5. Le 4 Vient Des Moyennes Deform_m%Btt + Des Moyennes Gauche-Droite
352
353                 Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))**2+ &
354                      (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))**2
355                 Chalk_1=Chalk_1**0.5
356
357                 Chalk_1=Ro*G*Chalk_1/4.*(E(K)**(Deform_m%Glen_exp(L)+1))/Cp(I,J,K)
358                 Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
359              enddo
360
361              ! Pour Shelves Et Streams,  On Ajoute Chalk_2
362              Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
363           end do
364        end do
365     end do
366
367     ! Chaleur Produite A La Base Par Le Glissement
368
369     do J=2,Ny_m
370        do I=2,Nx_m
371
372           if (mask_flot_m%Gz_mx(I,J)) then
373              Chalglissx(I,J)= abs(Ice_flow_m%U_Column_x(I,J)*Deform_m%Cisail_Basal_mx(I,J))
374           else
375              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)
376           endif
377
378           if (mask_flot_m%Gz_my(I,J)) then
379              Chalglissy(I,J)= abs(Ice_flow_m%U_Column_y(I,J)*Deform_m%Cisail_Basal_my(I,J))
380           else
381              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)
382           endif
383
384        end do
385     end do
386
387     !  Boundary Condition Ice-Rock Interface
388
389     K=Nz_m
390
391     do J=2,Ny_m-1
392        do I=2,Nx_m-1
393
394           !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
395           !        Dans La Derniere 1/2 Maille Et Par Le Glissement
396           !        Attention Phid Est >0 Et Ghf Est <0
397           if (.not.mask_flot_m%Flottant(I,J)) then
398
399              !           Phid Avec Fonte Sous Les Streams
400              !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
401              !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
402              !                        Chalbed*Fracq
403
404
405              ! Moyenne Phid Sur 4 Points Formulation (A2+B2)**0.5
406              !
407              Chalgliss_maj(I,J)=(Chalglissx(I,J)+Chalglissx(I+1,J))**2+  &
408                   (Chalglissy(I,J)+Chalglissy(I,J+1))**2
409
410              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)**0.5
411
412              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.5 ! Les Moyennes Droite Gauche
413
414              ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
415              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)
416
417
418              ! Flux Total A Rajouter À La Base
419              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)
420
421
422           else
423              therm_var_m%Phid(I,J)=0.
424           endif
425        end do
426     end do
427
428
429  case (3)     ! Q_prod_pente :  Pour Essayer D'Avoir La Chaleur En Alpha4
430     !-------------------------------------------------------------------------------------
431     ! Calcul Pour Essayer D'Avoir La Chaleur En Alpha4
432     ! Nom Du Sub Routine Avant Et Var Utilisee 
433     ! 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,&
434     !     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)
435
436
437     do J=2,Ny_m
438        do I=2,Nx_m
439           Pente2_maj(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I+1,J))**2 + &
440                (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I,J+1))**2
441           Pente2_maj(I,J)=Pente2_maj(I,J)*0.25    !Pour La Moyenne Sur Sdx Et Sdy
442        end do
443     end do
444
445     do K=2,Nz_m 
446        do J=2,Ny_m-1
447           do I=2,Nx_m-1
448              ! Calcul De La Chaleur De Deformation Selon Xx Yy Zz Et Xy
449              ! Pour Les Ice-Streams Et Ice-Shelves
450              ! Les Divers Eps Sont Calculés Dans Strain-Rate Pour L'Ensemble De La Grille
451              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.  &
452                   (mask_flot_m%Flgz_my(I,J).or.mask_flot_m%Flgz_my(I,J+1))) then
453
454                 Chal2_x(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_xx(I,J)**2
455                 Chal2_y(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*Deform_m%Veloc_Deform_yy(I,J)**2
456                 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
457
458                 !  Epsxy Est Calcule Sur Les Noeuds Mineur 1/2,1/2, Faire La Moyenne
459                 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))
460                 Chal2_xy(I,J,K)=2.*Deform_m%Viscosite(I,J,K)*(Chal2_xy(I,J,K)*0.25)**2
461
462              else   ! Glace Posee
463                 Chal2_x(I,J,K)=0.
464                 Chal2_y(I,J,K)=0.
465                 Chal2_z(I,J,K)=0.
466                 Chal2_xy(I,J,K)=0.
467              endif
468           end do
469        end do
470     end do
471
472     !  Partie Sia   Calcul De La Chaleur Produite Sur Chaque Demi Maille
473     !  Do L=1,size(Deform_m%Btt,4)!N1poly,N2poly
474     !     Do J=2,Ny_m
475     !        Do I=2,Nx_m
476     !
477     !          Ffx A 3 Dimensions ! 
478     !           Ffx(I,J,L)=Ddx(I,J,L)*Sdx(I,J)*Sdx(I,J)
479     !           Ffy(I,J,L)=Ddy(I,J,L)*Sdy(I,J)*Sdy(I,J)
480     !        End Do
481     !     End Do
482     !  End Do
483
484     do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
485        do K=2,Nz_m
486           do J=2,Ny_m
487              do I=2,Nx_m
488                 if ((.not.mask_flot_m%Flot_mx(I,J)).and.(.not.mask_flot_m%Gz_mx(I,J))) then
489                    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)
490
491                 else if (mask_flot_m%Gz_mx(I,J)) then ! Ice Streams
492                    Chalx(I,J,K,L)=0.
493
494                 else                     ! Ice Shelves 
495                    Chalx(I,J,K,L)=0.
496
497                 endif
498
499                 if ((.not.mask_flot_m%Flot_my(I,J)).and.(.not.mask_flot_m%Gz_my(I,J))) then
500                    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)
501
502                 else if (mask_flot_m%Gz_my(I,J)) then  ! Ice Streams
503                    Chaly(I,J,K,L)=0.
504
505                 else                      ! Ice Shelves
506                    Chaly(I,J,K,L)=0.
507                 endif
508
509              end do
510           end do
511        end do
512     end do
513
514     !         Nouvelle Formulation De Chaldef_maj(I,J,K), Le 4 Vient Des Moyennes
515     !         Deform_m%Btt Et Gauche Et Droite (Ou Haut Et Bas) Mais Il Faut Sommer
516     !         Les Productions X Et Y
517     !         Ancienne Formulation Chal=(Ro*G*H(I,J))**4*(Sx2+Sy2)*(Sx*Sx+Sy*Sy)
518
519     do K=2,Nz_m
520        do J=1,Ny_m-1
521           do I=1,Nx_m-1
522
523              ! Modif Christophe Mars 2000 : Chalx Et Chaly Sont A 4 Dim
524              Chaldef_maj(I,J,K)= 0.
525
526              ! Chalk_2 Pour Ice Shelves Et Ice Streams
527              Chalk_2=(Chal2_x(I,J,K)+Chal2_y(I,J,K) + &
528                   Chal2_z(I,J,K)+Chal2_xy(I,J,K))/Cp(I,J,K)
529
530              ! Chalk_1 Pour La Partie Posée
531              do L=1,size(Deform_m%Btt_Deform,4)!N1poly,N2poly
532
533                 ! On Somme La Chaleur Due Aux Diverses Lois De Déformation Et Celles En X Et Y
534                 ! La Somme Se Fait Par (Cx2+Cy2)** 0.5. Le 4 Vient Des Moyennes Deform_m%Btt + Des Moyennes Gauche-Droite
535
536                 ! On Fait La Moyenne Des Termes Ddx*Btt (*0.25 Pour Cette Moyenne, Le 0.5 Est Pour Les Btt)
537                 Chalk_1=(Chalx(I,J,K,L)+Chalx(I+1,J,K,L))+ &
538                      (Chaly(I,J,K,L)+Chaly(I,J+1,K,L))
539
540                 Chalk_1=Chalk_1*0.25*0.5
541
542                 ! On Multiplie Par La Pente Moyenne Au Carre Sur Le Noeud Majeur
543                 Chalk_1=Chalk_1*Pente2_maj(I,J)
544
545                 Chalk_1=Ro*G*Chalk_1*(E(K)**(Deform_m%Glen_exp(L)+1))/Cp(I,J,K)  ! Attention Plus De /4
546                 Chaldef_maj(I,J,K)= Chalk_1 + Chaldef_maj(I,J,K) 
547              enddo
548
549              ! Pour Shelves Et Streams,  On Ajoute Chalk_2
550              Chaldef_maj(I,J,K) = Chaldef_maj(I,J,K) + Chalk_2
551           end do
552        end do
553     end do
554
555     ! Chaleur Produite A La Base Par Le Glissement
556
557     do J=2,Ny_m
558        do I=2,Nx_m
559
560           if (mask_flot_m%Gz_mx(I,J)) then
561              Chalglissx(I,J)= abs(Ice_flow_m%U_Column_x(I,J)*Deform_m%Cisail_Basal_mx(I,J))
562           else
563              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)
564           endif
565
566           if (mask_flot_m%Gz_my(I,J)) then
567              Chalglissy(I,J)= abs(Ice_flow_m%U_Column_y(I,J)*Deform_m%Cisail_Basal_my(I,J))
568           else
569              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)
570           endif
571
572        end do
573     end do
574
575
576     !  Boundary Condition Ice-Rock Interface
577
578     K=Nz_m
579
580     do J=2,Ny_m-1
581        do I=2,Nx_m-1
582
583           !        Rajouter Un Flux De Chaleur Pour La Production Par Deformation
584           !        Dans La Derniere 1/2 Maille Et Par Le Glissement
585           !        Attention Phid Est >0 Et Ghf Est <0
586           if (.not.mask_flot_m%Flottant(I,J)) then
587
588              !           Phid Avec Fonte Sous Les Streams
589              !              Phid(I,J)=0.25*(Chalglissx(I,J)+Chalglissx(I+1,J)+ &
590              !                        Chalglissy(I,J)+Chalglissy(I,J+1)) + &
591              !                        Chalbed*Fracq
592
593
594              ! Moyenne Phid Sur 4 Points Formulation (A2+B2)**0.5  ! Je Garde Pour L'Instant A Revoir
595              !
596              Chalgliss_maj(I,J)=(Chalglissx(I,J)+Chalglissx(I+1,J))**2+  &
597                   (Chalglissy(I,J)+Chalglissy(I,J+1))**2
598
599              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)**0.5
600
601              Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.5 ! Les Moyennes Droite Gauche
602
603              ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
604              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)
605
606
607              ! Flux Total A Rajouter a La Base
608
609              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)
610
611
612           else
613              therm_var_m%Phid(I,J)=0.
614           endif
615        end do
616     end do
617
618
619  case (4)   ! Q_u_taub : Routine Qui concentre La Chaleur  au Socle
620     !-------------------------------------------------------------------------------------
621     ! Routine Qui Prend La Chaleur Remise Au Socle
622     ! Nom Du Sub Routine Avant Et Var Utilisee 
623     ! 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)
624
625
626     do J=2,Ny_m
627        do I=2,Nx_m
628
629
630           Pente_maj(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I+1,J))**2 + &          ! Pente
631                (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I,J+1))**2
632
633           Pente_maj(I,J)=Pente_maj(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
634           Pente_maj(I,J)=Pente_maj(I,J)**0.5
635
636           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
637                (Ice_flow_m%U_Column_y(I,J)+Ice_flow_m%U_Column_y(I,J+1))**2
638           Vit_maj(I,J)=Vit_maj(I,J)*0.25
639           Vit_maj(I,J)=Vit_maj(I,J)**0.5
640
641           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
642                (Ice_flow_m%Veloc_y(I,J,Nz_m)+Ice_flow_m%Veloc_y(I,J+1,Nz_m))**2
643           Uslid_maj(I,J)=Uslid_maj(I,J)*0.25
644           Uslid_maj(I,J)=Uslid_maj(I,J)**0.5
645
646        end do
647     end do
648
649
650
651     Chaldef_maj(:,:,:)=0.
652     Chalgliss_maj(:,:)=Ro*G*Geom_m%Thickness(:,:)*Pente_maj(:,:)*Uslid_maj(:,:)
653     Chaldef_maj(:,:,Nz_m)=Ro*G*Geom_m%Thickness(:,:)*Pente_maj(:,:)*Vit_maj(:,:)-Chalgliss_maj(:,:)
654
655     ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
656
657     Chalgliss_maj(:,:)=Chalgliss_maj(:,:)*exp((T_m%Temperature(:,:,Nz_m)-T_m%Temp_melting(:,:,Nz_m))*Ecart_phid)
658
659     ! Flux Total A Rajouter a Base
660     therm_var_m%Phid(:,:)=Chaldef_maj(:,:,Nz_m)+Chalgliss_maj(:,:)
661
662
663  case (5)    ! Q_u_taub_stag : concentre La Chaleur Au Socle mais Sur Les Mailles Staggered
664     !----------------------------------------------------------------------------------
665     ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
666     ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
667     ! Nom Du Sub Routine Avant Et Var Utilisé 
668     ! 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)
669
670
671     do J=2,Ny_m
672        do I=2,Nx_m
673
674           Pente_stag(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I,J-1))**2 + &          ! Pente
675                (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I-1,J))**2
676
677           Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
678           Pente_stag(I,J)=Pente_stag(I,J)**0.5
679
680           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
681                (Ice_flow_m%U_Column_y(I,J)+Ice_flow_m%U_Column_y(I-1,J))**2
682           Vit_stag2(I,J)=Vit_stag2(I,J)*0.25
683           Vit_stag2(I,J)=Vit_stag2(I,J)**0.5
684
685           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
686                (Ice_flow_m%Veloc_y(I,J,Nz_m)+Ice_flow_m%Veloc_y(I-1,J,Nz_m))**2
687           Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
688           Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
689
690           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
691
692        end do
693     end do
694
695
696     Qslid(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Uslid_stag(:,:)
697     Qdef2(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Vit_stag2(:,:)-Qslid(:,:)
698     Chaldef_maj(:,:,:)=0.
699
700
701     ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
702
703     do J=2,Ny_m-1
704        do I=2,Nx_m-1
705           Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
706           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
707
708           Chaldef_maj(I,J,Nz_m)=(Qdef2(I,J)+Qdef2(I+1,J+1))+(Qdef2(I,J+1)+Qdef2(I+1,J))
709           Chaldef_maj(I,J,Nz_m)=Chaldef_maj(I,J,Nz_m)*0.25
710
711           ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
712           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)
713
714        end do
715     end do
716
717
718
719
720     !
721     !Write(126,*)'Time=',Time
722     !J=41
723     !Do I=20,30
724     !   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) &
725     !       , Qdef2(I,J+1),Qdef2(I+1,J)
726     !End Do
727
728     ! Flux Total A Rajouter À La Base
729     therm_var_m%Phid(:,:)= Chaldef_maj(:,:,Nz_m)+Chalgliss_maj(:,:)
730
731  case (6)    ! Q_sia_stag  : concentre La Chaleur  Au Socle Mais Sur Les Mailles Staggered
732              ! Au Milieux Des Mailles Vitesses
733
734     !----------------------------------------------------------------------------------
735     ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
736     ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
737     ! Nom Du Sub Routine Avant Et Var Utilisé 
738     ! Q_sia_stag(T_m%T,T_m%Tpmp,Geom_m%H,Geom_m%Slop_x,Geom_m%Slop_y,therm_var_m%Phid)
739
740
741     Pente_stag = 0.
742     Vit_stag3   = 0.
743     Uslid_stag = 0.
744     Hmxy  = 0.
745     Qslid = 0.
746     Qdef3 = 0.
747
748     do J=2,Ny_m
749        do I=2,Nx_m
750
751           !     Variables 2d
752           Pente_stag(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I,J-1))**2 + &          ! Pente
753                (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I-1,J))**2
754
755           Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
756           Pente_stag(I,J)=Pente_stag(I,J)**0.5
757
758           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
759                (Ice_flow_m%Veloc_y(I,J,Nz_m)+Ice_flow_m%Veloc_y(I-1,J,Nz_m))**2
760           Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
761           Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
762
763           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
764
765           !     En Vertical : Calcul De La Deformation Par Differentiation Des Vitesses
766
767           do K=1,Nz_m
768
769              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
770                   (Ice_flow_m%Veloc_y(I,J,K)+Ice_flow_m%Veloc_y(I-1,J,K))**2            ! A Tous Niveaux K
771
772              Vit_stag3(I,J,K)=Vit_stag3(I,J,K)*0.25
773              Vit_stag3(I,J,K)=Vit_stag3(I,J,K)**0.5
774           end do
775
776           Qdef3(I,J,1)=0.
777           do K=2,Nz_m-1
778              Qdef3(I,J,K)=Vit_stag3(I,J,K-1)-Vit_stag3(I,J,K+1)           ! Difference Des Vitesses
779              Qdef3(I,J,K)=Ro*G*Pente_stag(I,J)*E(K)*Qdef3(I,J,K)/2./De    ! Gamma Tau
780           end do
781
782
783           Qdef3(I,J,Nz_m)=Vit_stag3(I,J,Nz_m-1)-Vit_stag3(I,J,Nz_m)              ! Pour Le Fond Differentiation
784           Qdef3(I,J,Nz_m)=Ro*G*Pente_stag(I,J)*Qdef3(I,J,K)/De               ! Sur Une Seule Maille
785        end do
786     end do
787
788
789     Qslid(:,:)=Ro*G*Hmxy(:,:)*Pente_stag(:,:)*Uslid_stag(:,:)
790
791     ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
792
793     do J=2,Ny_m-1
794        do I=2,Nx_m-1
795           Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
796           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
797
798           ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
799           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)
800
801           ! Chaleur Due A La Defromation
802           do K=1,Nz_m
803              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))
804              Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)*0.25
805              Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)/Cp(I,J,K)
806           end do
807
808
809        end do
810     end do
811
812     !Write(126,*)'Time=',Time
813     !J=41
814     !Do I=20,30
815     !   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) &
816     !       , Qdef2(I,J+1),Qdef2(I+1,J)
817     !End Do
818
819     ! Flux Total A Rajouter À La Base
820
821     therm_var_m%Phid(:,:)=Chalgliss_maj(:,:)+Chaldef_maj(:,:,Nz_m)*Fracq*Geom_m%Thickness(:,:)*Cp(:,:,Nz_m)
822
823  case (7)   !  Q_all_stag idem case 6 mais Le Calcul De Taub Tient Compte Des Noeuds Grzmx
824     !----------------------------------------------------------------------------------
825     ! Routine Qui Prend La Chaleur Remise Au Socle Mais Sur Les Mailles Staggered
826     ! Les Noeuds Stag Sont Ceux Au Milieux Des Mailles Vitesses
827     ! Le Calcul De Taub Tient Compte Des Noeuds Grzmx     
828     ! Nom Du Sub Routine Avant Et Var Utilisé 
829     ! 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)
830
831     where (mask_flot_m%Flgz_mx(:,:))
832        Tox(:,:)=Deform_m%Cisail_Basal_mx(:,:)
833     elsewhere
834        Tox(:,:)=Geom_m%Slop_x(:,:)*Geom_m%Thick_mx(:,:)*Ro*G
835     end where
836
837     where (mask_flot_m%Flgz_my(:,:))
838        Toy(:,:)=Deform_m%Cisail_Basal_my(:,:)
839     elsewhere
840        Toy(:,:)=Geom_m%Slop_y(:,:)*Geom_m%Thick_my(:,:)*Ro*G
841     end where
842
843     do J=2,Ny_m
844        do I=2,Nx_m
845
846           !     Variables 2d
847           Pente_stag(I,J)=(Geom_m%Slop_x(I,J)+Geom_m%Slop_x(I,J-1))**2 + &    ! Pente
848                (Geom_m%Slop_y(I,J)+Geom_m%Slop_y(I-1,J))**2
849
850           Pente_stag(I,J)=Pente_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
851           Pente_stag(I,J)=Pente_stag(I,J)**0.5
852
853
854           Tob_stag(I,J)=(Tox(I,J)+Tox(I,J-1))**2 + &          ! Pente
855                (Toy(I,J)+Toy(I-1,J))**2
856
857           Tob_stag(I,J)=Tob_stag(I,J)*0.25    ! 0.25pour La Moyenne Sur Sdx Et Sdy
858           Tob_stag(I,J)=Tob_stag(I,J)**0.5
859
860           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
861                (Ice_flow_m%Veloc_y(I,J,Nz_m)+Ice_flow_m%Veloc_y(I-1,J,Nz_m))**2
862           Uslid_stag(I,J)=Uslid_stag(I,J)*0.25
863           Uslid_stag(I,J)=Uslid_stag(I,J)**0.5
864
865           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
866
867           !     En Vertical : Calcul De La Deformation Par Differentiation Des Vitesses
868
869           do K=1,Nz_m
870
871              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
872                   (Ice_flow_m%Veloc_y(I,J,K)+Ice_flow_m%Veloc_y(I-1,J,K))**2            ! A Tous Niveaux K
873
874              Vit_stag3(I,J,K)=Vit_stag3(I,J,K)*0.25
875              Vit_stag3(I,J,K)=Vit_stag3(I,J,K)**0.5
876           end do
877
878           Qdef3(I,J,1)=0.
879           do K=2,Nz_m-1
880              Qdef3(I,J,K)=Vit_stag3(I,J,K-1)-Vit_stag3(I,J,K+1)            ! Difference Des Vitesses
881              Qdef3(I,J,K)=Ro*G*Pente_stag(I,J)*E(K)*Qdef3(I,J,K)/2./De    ! Gamma Tau
882           end do
883
884
885           Qdef3(I,J,Nz_m)=Vit_stag3(I,J,Nz_m-1)-Vit_stag3(I,J,Nz_m)               ! Pour Le Fond Differentiation
886           Qdef3(I,J,Nz_m)=Ro*G*Pente_stag(I,J)*Qdef3(I,J,K)/De               ! Sur Une Seule Maille
887        end do
888     end do
889
890     Qslid(:,:)=Tob_stag*Uslid_stag(:,:)
891
892     ! On Fait La Moyenne De La Chaleur Produite Sur Les Mailles Stag
893
894     do J=2,Ny_m-1
895        do I=2,Nx_m-1
896           Chalgliss_maj(I,J)=(Qslid(I,J)+Qslid(I+1,J+1))+(Qslid(I,J+1)+Qslid(I+1,J))
897           Chalgliss_maj(I,J)=Chalgliss_maj(I,J)*0.25
898
899           ! Plus La Base Est Loin Du Point De Fusion Moins La Chaleur De Glissement Est Prise En Compte
900           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)
901
902           ! Chaleur Due A La Defromation
903           do K=1,Nz_m         
904              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))
905              Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)*0.25
906              Chaldef_maj(I,J,K)=Chaldef_maj(I,J,K)/Cp(I,J,K)
907           end do
908
909        end do
910     end do
911
912     !
913     !Write(126,*)'Time=',Time
914     !J=41
915     !Do I=20,30
916     !   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) &
917     !       , Qdef2(I,J+1),Qdef2(I+1,J)
918     !End Do
919
920     ! Flux Total A Rajouter À La Base
921     therm_var_m%Phid(:,:)=Chalgliss_maj(:,:)+Chaldef_maj(:,:,Nz_m)*Fracq*Geom_m%Thickness(:,:)*Cp(:,:,Nz_m)
922
923  end select
924end subroutine Qprod
925
926
927!///////////////////////////////////////////////////////////////////////
928!///////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.