source: trunk/SOURCES/Temperature-routines/Qprod_icetemp.f90 @ 111

Last change on this file since 111 was 71, checked in by dumas, 8 years ago

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