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

Last change on this file since 29 was 24, checked in by dumas, 9 years ago

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