source: trunk/SOURCES/Old-sources/icetemp-0.5_mod.F90 @ 111

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

initial import GRISLI trunk
File size: 45.8 KB
Line 
1
2module ice_temp_mod                 ! mis en module le 13 mars 2007 (Cat)
3
4! la production de chaleur est mise dans une routine du module
5
6use module3d_phy
7use tridiagmod     ! module dans lequel est la routine tridiag
8use module_choix
9
10
11implicit none
12
13integer :: iec,jec,ngr
14integer :: nfracq   ! exposant fracq
15
16integer :: iq       ! choix du type de routine chaleur
17
18real :: sx,sy,sx2,sy2,deh22,tss
19real ::  dou,dah,duu,dzz,dzi,chalbed,ct_bas,ct_haut
20real :: acof1,bcof1,ccof1,s0mer ! pour la température de l'eau de mer
21                                ! à la base de l'ice shelf
22real :: advecx,advecy
23
24!integer :: ifail             ! permet de detecter les erreurs
25real :: chalk_1               ! utilise pour le calcul de chalk : glace posée
26real :: chalk_2               ! utilise pour le calcul de chalk : ice streams et ice shelves
27real :: chalbed_1             ! utilise pour le calcul de chalbed
28real :: coefadv               ! pour limiter le flux de chaleur horiz.
29
30real :: ecart_phid            ! pour limiter l'impact de la chaleur de glissement si base froide
31                              ! ancienne valeur 0.5.
32
33! tableaux 1D
34! _______________
35real,dimension(nzm+1) :: abis,bbis,cbis,rbis,hbis
36
37real,dimension(nn) :: aa=0.              ! work arrays for tridiag
38real,dimension(nn) :: bb=0.              ! work arrays for tridiag
39real,dimension(nn) :: cc=1.              ! work arrays for tridiag
40real,dimension(nn) :: rr=0.              ! work arrays for tridiag
41real,dimension(nn) :: hh=0.              ! work arrays for tridiag
42
43
44! tableaux de travail 2D
45! ___________________________
46real, dimension(nx,ny):: tbmer   ! temperature de la mer a la base de l'ice shelf
47real,dimension(nx,ny) :: alpha   ! pente locale sur noeud majeur
48real,dimension(nx,ny) :: ubar    ! vitesse sur noeud majeur
49
50real,dimension(nx,ny) :: chalglissx,chalglissy ! chaleur de glissement
51
52
53
54! Tableaux de travail 3D
55! __________________________
56
57real, dimension(nx,ny,nz,n1poly:n2poly) :: chalx, chaly
58real, dimension(nx,ny,nz) :: chal2_x, chal2_y, chal2_z, chal2_xy
59! utilisé pour calculer la chaleur de deformation selon xx yy zz et xy
60real, dimension(nx,ny,n1poly:n2poly) :: ffx,ffy
61real, dimension(nx,ny,nz+nzm) :: T3D_new
62
63
64
65
66contains
67!------------------------------------------------------------------------------------------
68subroutine init_icetemp
69
70!     variables dependant du pas de temps dtt
71
72  ctm=dtt*cm/rom/cpm/dzm/dzm
73  dttdx=dtt/dx
74
75
76!________________________________________________________________
77! La temperature du point de congelation de l'eau de mer
78! est tiree de Jenkins (1991), Tb est en degres C
79! _______________________________________________________________
80
81  acof1 = - 0.0575
82  bcof1 = 0.0901
83  ccof1 = 7.61E-4
84  s0mer = 34.75
85
86
87!!     nouvelle formulation de fracq
88!  nfracq=15  ! ancienne formulation
89nfracq=1
90fracq=(1.-(1.-de/2.)**nfracq)/nfracq
91
92write(num_rep_42,*) 'calcul des temperatures'
93write(num_rep_42,*) 'nfracq=',nfracq
94
95iq=4
96
97
98if (iq.eq.1) write(num_rep_42,*) 'chaleur demi maille, iq=',iq
99if (iq.eq.2) write(num_rep_42,*) 'chaleur centree, iq=',iq
100if (iq.eq.3) write(num_rep_42,*) 'chaleur pente, iq=',iq
101if (iq.eq.4) write(num_rep_42,*) 'chaleur u_taub, iq=',iq
102if (iq.eq.5) write(num_rep_42,*) 'chaleur u_taub_stag, iq=',iq
103if (iq.eq.6) write(num_rep_42,*) 'chaleur Q_SIA_stag,iq=',iq
104if (iq.eq.7) write(num_rep_42,*) 'chaleur Q_all_stag,iq=',iq
105write(num_rep_42,*)
106
107ecart_phid=0.5
108write(num_rep_42,*) ' si base froide, la chaleur de glissement est en exp(delta T * ecart_phid)'
109write(num_rep_42,*) 'ecart_phid=',ecart_phid
110
111write(num_rep_42,*) '---------------------------------------------------------------------------'
112
113end subroutine init_icetemp
114
115
116!------------------------------------------------------------------------------------------
117subroutine icetemp()
118
119  if (itracebug.eq.1)  call tracebug(' Entree dans routine icetemp')
120
121  !________________________________________________________________
122  ! La temperature du point de congelation de l'eau de mer
123  ! est tiree de Jenkins (1991), Tb est en degres C
124  ! _______________________________________________________________
125
126  do j=1,ny
127     do i=1,nx
128        if (flot(i,j)) &
129             tbmer(i,j) = acof1*s0mer + bcof1 + ccof1*H(i,j)*ro/row
130     end do
131  end do
132
133
134  ! Appel aux proprietes thermiques de la glace : conductivite, capacite claorifique -tpmp
135
136  call thermal_prop
137
138
139
140
141  if (iq.eq.1) call  Q_prod_demi
142  if (iq.eq.2) call  Q_prod_centre
143  if (iq.eq.3) call  Q_prod_pente
144  if (iq.eq.4) call  Q_U_taub
145  if (iq.eq.5) call  Q_U_taub_stag
146  if (iq.eq.6) call  Q_SIA_stag
147  if (iq.eq.7) call  Q_all_stag
148
149  !     dans le socle : les éléments de la matrice tridiag sont invariants dans l'espace
150  do k=nz+1,nz+nzm-1
151     aa(k)=-ctm
152     bb(k)=1.+2.*ctm
153     cc(k)=-ctm
154  end do
155
156  !     conditions aux limites pour le socle
157  aa(nz+nzm)=-1.
158  bb(nz+nzm)=1.
159  cc(nz+nzm)=0.
160
161
162  !     boucle de resolution des temperatures colonne par colonne
163  do j=2,ny-1
164     do i=2,nx-1
165        rr(nz+nzm)=-dzm*ghf(i,j)/cm  ! dépend de i,j pour un flux géothermique variable
166        !      -----------------------------cas general (H>10m)
167
168        if (H(i,j).gt.10.) then
169
170           T(i,j,1)=min(0.,Ts(i,j))
171
172           aa(1)=0.
173           bb(1)=1.
174           cc(1)=0.
175           rr(1)=T(i,j,1)
176
177
178           !        variables de calcul dans la glace
179
180           deh22=2.*de*h(i,j)**2
181           dou=dtt/de/de/h(i,j)/h(i,j)
182           dah=dtt/h(i,j)/de
183           duu=dtt/2./de
184           ct_bas=2.*(Ct(i,j,1)*Ct(i,j,2))/(Ct(i,j,1)+Ct(i,j,2))
185
186
187           !        indice upwind pour l'advection horizontale
188           !        Le test se fait maintenant sur les vitesses
189
190
191           !           ecoulement vers l'OUEST
192           if ((uxbar(i+1,j).lt.0.).and.(uxbar(i,j).le.0.)) then
193              i1=i+1
194              i2=i
195              iec=0
196
197              !           ecoulement vers l'EST
198           else if ((uxbar(i+1,j).ge.0.).and.(uxbar(i,j).gt.0.)) then
199              i1=i
200              i2=i-1
201              iec=0
202              !           ecoulement  convergent
203           else if ((uxbar(i+1,j).lt.0.).and.(uxbar(i,j).gt.0.)) then
204              i1=i+1
205              i2=i-1
206              iec=1
207              !           ecoulement divergent
208           else if ((uxbar(i+1,j).ge.0.).and.(uxbar(i,j).le.0.)) then
209              i1=i
210              i2=i
211              iec=-1
212           endif
213
214           !           ecoulement vers lE SUD
215           if ((uybar(i,j+1).lt.0.).and.(uybar(i,j).le.0.)) then
216              j1=j+1
217              j2=j
218              jec=0
219              !           ecoulement vers lE NORD
220           else if ((uybar(i,j+1).ge.0.).and.(uybar(i,j).gt.0.)) then
221              j1=j
222              j2=j-1
223              jec=0
224              !           ecoulement  convergent
225           else if ((uybar(i,j+1).lt.0.).and.(uybar(i,j).gt.0.)) then
226              j1=j+1
227              j2=j-1
228              jec=1
229              !           ecoulement  divergent
230           else if ((uybar(i,j+1).ge.0.).and.(uybar(i,j).le.0.)) then
231              j1=j
232              j2=j
233              jec=-1
234           endif
235
236
237           do k=2,nz-1
238
239
240              dzz=dou/cp(i,j,k)
241
242              ! conductivité au milieu des mailles
243              Ct_haut=ct_bas
244              Ct_bas=2.*(Ct(i,j,k)*Ct(i,j,k+1))/(Ct(i,j,k)+Ct(i,j,k+1))
245
246              !         advection verticale centree
247              aa(k)=-dzz*ct_haut-uzr(i,j,k)*dah/2.
248              bb(k)=1.+dzz*(ct_haut+ct_bas)
249              cc(k)=-dzz*ct_bas+uzr(i,j,k)*dah/2.
250
251              !         avection selon x
252              !         ----------------
253              if (iec.eq.0) then
254
255                 !    coefadv limite l'advection afin de ne pas amener plus de
256                 !    chaleur qu'il n'y en a (limite a delta T/2)
257                 !  Attention : avril 2010 :
258                 !  on enleve cette limitation qui a l'air de modifier sensiblement les resultats
259                 ! (comparaison avec Hassine)
260
261                 !                 coefadv=max(1.,abs(uxbar(I1,J)))
262                 !                 coefadv=min(0.5/dttdx/coefadv,1.)
263                 coefadv=1.                                    !correction adv avril 2010
264
265                 advecx = coefadv*(UX(I1,J,K)*(T(I1,J,K)-T(I2,J,K)))
266
267              else
268                 !                 coefadv=max(1.,abs(uxbar(i,j)))
269                 !                 coefadv=min(0.5/dttdx/coefadv,1.)
270                 coefadv=1.                                    !correction adv avril 2010
271                 advecx=coefadv*ux(i,j,k)*(T(i,j,k)-T(i-1,j,k))
272                 coefadv=max(1.,abs(uxbar(i+1,j)))
273                 coefadv=min(0.5/dttdx/coefadv,1.)
274                 coefadv=1.                                    !correction adv avril 2010
275                 advecx=iec*(advecx+ &
276                      coefadv*ux(i+1,j,k)*(T(i+1,j,k)-T(i,j,k)))
277              endif
278
279              !         avection selon y
280              !         ----------------
281              if (jec.eq.0) then
282                 !                 coefadv=max(1.,abs(uybar(i,j1)))
283                 !                 coefadv=min(0.5/dttdx/coefadv,1.)
284                 coefadv=1.                                    !correction adv avril 2010
285                 advecy = coefadv*uy(i,j1,k)*(T(i,j1,k)-T(i,j2,k))
286
287              else
288                 !                 coefadv=max(1.,abs(uybar(i,j)))
289                 !                 coefadv=min(0.5/dttdx/coefadv,1.)
290                 coefadv=1.                                    !correction adv avril 2010
291                 advecy=coefadv*uy(i,j,k)*(T(i,j,k)-T(i,j-1,k))
292                 coefadv=max(1.,abs(uybar(i+1,j)))
293                 coefadv=min(0.5/dttdx/coefadv,1.)
294                 advecy=jec*(advecy+ &
295                      (uy(i,j+1,k)*(T(i,j+1,k)-T(i,j,k))))
296                 advecy=jec*((uy(i,j,k)*(T(i,j,k)-T(i,j-1,k)))+ &
297                      coefadv*uy(i,j+1,k)*(T(i,j+1,k)-T(i,j,k)))
298              endif
299
300              rr(k)=T(i,j,k)+chaldef_maj(i,j,k)*dtt
301              rr(k)=rr(k)-dttdx*(advecx+advecy)
302
303
304           end do
305           !        fin boucle k, ecriture matrice tridiag dans la glace
306
307
308
309
310
311           ! Conditions aux limites à la base de la glace
312           !________________________________________________
313
314           !           base froide
315           if ( ((ibase(i,j).eq.1).or.(ibase(i,j).eq.4) &
316                .or. ((ibase(i,j).eq.5).and.(T(i,j,nz).lt.TPMP(i,j,nz)))) &
317                .and..not.flot(i,j)) then
318
319              ibase(i,j)=1
320              bb(nz)=1.
321
322              if (ncond.eq.1) then ! avec socle
323                 dzi=h(i,j)*de*cm/ct(i,j,nz)
324                 aa(nz)=-dzm/(dzm+dzi)
325                 cc(nz)=-dzi/(dzm+dzi)
326                 rr(nz)=h(i,j)*de*phid(i,j)*dzm/ct(i,j,nz)/(dzm+dzi)
327
328              else                 ! sans socle
329                 aa(nz)=-1.
330                 cc(nz)=0.
331                 rr(nz)=-(ghf(i,j)-phid(i,j))/ct(i,j,nz)*h(i,j)*de
332              endif
333
334           else
335
336              !          base temperee ou Shelf
337              !          ----------------------               
338              if (ibase(i,j).eq.5)  ibase(i,j)=2
339              ibase(i,j)=max(ibase(i,j),2)
340              aa(nz)=0.
341              bb(nz)=1.
342              cc(nz)=0.
343
344              if (.not.flot(i,j)) then
345                 rr(nz)=TPMP(i,j,nz)
346              else
347                 rr(nz)=Tbmer(i,j)
348              end if
349
350           endif
351
352           if (ncond.eq.1) then ! avec socle
353              do k=nz+1,nz+nzm-1
354                 rr(k)=T(i,j,k)
355              end do
356              call TRIDIAG (aa,bb,cc,rr,hh,nz+nzm,ifail)
357
358           else                 ! sans socle
359
360
361              call TRIDIAG (AA,BB,CC,RR,HH,NZ,iFAIL)
362           endif
363
364           ! test permettant d'ecrire les variables quand il y a une erreur
365           ! iFAIL = 1 : erreur dans tridiag
366
367           if (ifail.eq.1) then
368              write(6,*)'erreur dans TRIDIAG'
369              write(6,*)'i = ',i,' j = ',j
370              write(6,*) 'a, b, c, r, u'
371              do k=1,nz+nzm
372                 write(6,*)'k = ',k
373                 write(6,*) aa(k),bb(k),cc(k),rr(k),hh(k)
374              end do
375              stop
376           end if
377
378           !         temperature dans le socle si ncond=0, lineaire avec le gradient ghf
379           if (ncond.eq.0) then
380              do k=nz+1,nz+nzm
381                 hh(k)=hh(nz)-dzm*(k-nz)*ghf(i,j)/cm
382              end do
383           endif
384
385
386
387           do k=1,nz+nzm
388              T3D_new(i,j,k)=hh(k)
389           end do
390           tbdot(i,j)=(T3D_new(i,j,nz)-t(i,j,nz))/dtt
391
392
393           !      ------------------------------------------------- glace tres fine (H<10m ou H=0)
394        else if (H(i,j).le.10.) then
395
396           !       Pour eviter des problemes lorsque H est inferieur a 10 m
397           !       profil lineaire, pente=flux geothermique
398           !       Le cas sans glace est traite de la meme facon  (DOU=0)
399
400           ! ........................................ avec calcul dans le socle
401           if (NCOND.eq.1.AND..NOT.FLOT(I,J)) then
402
403              if (H(I,J).gt.0.) then
404                 !           gradient dans le socle
405                 DOU=(T(I,J,NZ+1)-T(I,J,NZ))/DZM*CM
406                 DOU=DOU/CT(I,J,NZ)*DE*H(I,J)
407
408                 TSS=min(0.,TS(I,J))
409                 do K=1,NZ
410                    T3D_new(I,J,K)=TSS+DOU*(K-1.)
411                 end do
412
413                 !         pas de glace
414              else
415                 TSS=Ts(I,J)
416                 do K=1,NZ
417                    T3D_new(I,J,K)=TSS
418                 end do
419              endif
420
421              !         calcul dans le socle meme s'il n'y a pas de glace
422              RR(NZ)=T3D_new(I,J,NZ)
423              AA(NZ)=0.
424              BB(NZ)=1.
425              CC(NZ)=0.
426
427              do K=NZ+1,NZ+NZM-1
428                 RR(K)=T(I,J,K)
429              end do
430
431              !         creation de nouveaux tableaux juste pour tridiag     
432              do k=1,NZM+1
433                 ABIS(k)=AA(NZ-1+k)
434                 BBIS(k)=BB(NZ-1+k)
435                 CBIS(k)=CC(NZ-1+k)
436                 RBIS(k)=RR(NZ-1+k)
437              end do
438
439              call TRIDIAG (ABIS,BBIS,CBIS,RBIS,HBIS,NZM+1,iFAIL)
440              ! test permettant d'ecrire les variables quand il y a une erreur
441              ! iFAIL = 1 : erreur dans tridiag
442              IF (iFAIL.EQ.1) THEN
443                 write(6,*)'Erreur dans TRIDIAG'
444                 write(6,*)'I = ',I,' J = ',J
445                 write(6,*) 'A, B, C, R, U'
446                 do k=1,NZ+NZM
447                    write(6,*)'k = ',k
448                    write(6,*) ABIS(k),BBIS(k),CBIS(k),RBIS(k),HBIS(k)
449                 end do
450                 STOP
451              END IF
452
453              do k=1,NZM+1
454                 HH(NZ-1+k)=HBIS(k)
455              end do
456
457              ! ........................................ sans calcul dans le socle
458           else
459              !                        calotte posee                   
460              if ((H(I,J).gt.0.).AND..NOT.FLOT(I,J)) then
461                 DOU=-GHF(i,j)/CT(I,J,NZ)*DE*H(I,J)
462                 TSS=min(0.,TS(I,J))
463                 do K=1,NZ
464                    T3D_new(I,J,K)=TSS+DOU*(K-1.)
465                 end do
466
467                 !                        shelf
468              else if ((H(I,J).gt.0.).AND.FLOT(I,J)) then
469                 TSS=min(0.,TS(I,J))
470                 DOU=(TBMER(I,J)-TSS)*DE
471                 do K=1,NZ
472                    T3D_new(I,J,K)=TSS+DOU*(K-1.)
473                 end do
474
475                 !                       pas de glace
476              else
477                 TSS=Ts(I,J)
478                 do K=1,NZ
479                    T3D_new(I,J,K)=TSS
480                 end do
481              endif
482
483              !         temperature dans le socle, lineaire avec le gradient ghf
484
485              IF (.NOT.FLOT(I,J)) THEN
486                 do K=NZ+1,NZ+NZM
487                    HH(K)=T3D_new(I,J,NZ)-DZM*(K-NZ)*GHF(i,j)/CM
488                 end do
489
490              ELSE
491
492                 DO K=NZ+1,NZ+NZM
493                    HH(K)=Tbmer(i,j)-DZM*(K-NZ)*GHF(i,j)/CM
494                 END DO
495              ENDIF
496
497
498           endif
499           !         fin du test socle
500
501           do K=NZ+1,NZ+NZM
502              T3D_new(I,J,K)=HH(K)
503           end do
504
505
506           TBDOT(I,J)=(T3D_new(I,J,NZ)-T(I,J,NZ))/DTT
507           BMELT(I,J)=0.
508           IBASE(I,J)=5
509           phid(i,j)=0.
510
511        endif
512        !       fin du test sur l'epaisseur de glace
513
514     end do
515  end do
516
517
518  do j=1,ny
519     do i=1,nx
520        H1(I,J)=H(I,J)
521        B1(I,J)=B(I,J)
522        TPMP(I,J,1)=0.
523     end do
524  end do
525
526
527  !     attribution finale de la temperature
528  ! limitation à 3deg de variation par dtt et pas plus froid que -70
529
530  do k=1,nz+nzm
531     do j=1,ny
532        do i=1,nx
533           if (H(i,j).gt.10.) then
534
535              Tdot(k)=T3D_new(i,j,k)-T(i,j,k)
536
537              if ((k.gt.1).and.(Tdot(k).lt.-3.)) Tdot(k)=-3.
538              if ((k.gt.1).and.(Tdot(k).gt.3.)) Tdot(k)=3.
539
540              T(i,j,k)=T(i,j,k)+Tdot(k)
541
542              if (T(i,j,k).lt.-70.) then
543                 T(i,j,k)=-70.
544              endif
545           endif
546
547        end do
548     end do
549  end do
550
551  do k=2,nz
552     do j=1,ny
553        do i=1,nx
554           if (T(i,j,k).gt.TPMP(i,j,k)) then
555              T(i,j,k)=TPMP(i,j,k)
556              ibase(i,j)=2
557           endif
558        end do
559     end do
560  end do
561
562  !     temperature sur les bords :
563  !     -------------------------
564
565
566  do k=1,nz
567     do j=1,ny
568        !        T(nx,j,k)= T(nx-1,j,k)
569        !        T(1,j,k)= T(2,j,k)
570        T(nx,j,k)= Ts(nx,j)
571        T(1,j,k)= Ts(1,j)
572     end do
573     !
574     do i=1,nx
575        !        T(i,1,k)=T(i,2,k)
576        !        T(i,ny,k)=T(i,ny-1,k)
577        T(i,1,k)=Ts(i,1)
578        T(i,ny,k)=Ts(i,ny)
579     end do
580  end do
581
582
583  do k=1,nz
584     do j=1,ny
585        do i=1,nx
586           if (H(i,j).le.1.) T(i,j,k)=Ts(i,j)
587        end do
588     end do
589  end do
590
591
592  return
593end subroutine icetemp
594!-------------------------------------------------------------------------------------
595! calcul avec l'ancienne méthode sur les demi-mailles
596
597subroutine Q_prod_demi
598
599  do k=2,nz
600     do j=2,ny-1
601        do i=2,nx-1
602
603
604! calcul de la chaleur de deformation selon xx yy zz et xy
605! pour les ice-streams et ice-shelves
606! les divers eps sont calculés dans strain-rate pour l'ensemble de la grille
607
608        if ((flot(i,j).or.flgzmx(i,j).or.flgzmx(i+1,j)).or.  &
609                        (flgzmy(i,j).or.flgzmy(i,j+1))) then
610
611              chal2_x(i,j,k)=2.*visc(i,j,k)*epsxx(i,j)**2
612              chal2_y(i,j,k)=2.*visc(i,j,k)*epsyy(i,j)**2
613              chal2_z(i,j,k)=2.*visc(i,j,k)*(-epsxx(i,j)-epsyy(i,j))**2
614
615!  epsxy est calcule sur les noeuds mineur 1/2,1/2, faire la moyenne
616              chal2_xy(i,j,k)=(epsxy(i,j)+epsxy(i+1,j)+epsxy(i+1,j+1)+epsxy(i,j+1))
617              chal2_xy(i,j,k)=2.*visc(i,j,k)*(chal2_xy(i,j,k)*0.25)**2
618
619        else   ! glace posée
620              chal2_x(i,j,k)=0.
621              chal2_y(i,j,k)=0.
622              chal2_z(i,j,k)=0.
623              chal2_xy(i,j,k)=0.
624        endif
625        end do
626     end do
627  end do
628
629!  partie SIA   calcul de la chaleur produite sur chaque demi maille
630  do l=n1poly,n2poly
631     do j=2,ny
632        do i=2,nx
633
634!          ffx a 3 dimensions !
635           ffx(i,j,l)=ddx(i,j,l)*sdx(i,j)*sdx(i,j)
636           ffy(i,j,l)=ddy(i,j,l)*sdy(i,j)*sdy(i,j)
637        end do
638     end do
639  end do
640
641
642
643  do l=n1poly,n2poly
644     do k=2,nz
645        do j=2,ny
646           do i=2,nx
647              if ((.not.flotmx(i,j)).and.(.not.gzmx(i,j))) then
648                 chalx(i,j,k,l)=(btt(i-1,j,k,l)+btt(i,j,k,l))*ffx(i,j,l) !&
649!                      *ro*g*e(k)**(glen(l)+1)/cp(i,j,k)
650
651              else if (gzmx(i,j)) then ! ice streams
652                  chalx(i,j,k,l)=0.
653
654              else                     ! ice shelves 
655                  chalx(i,j,k,l)=0.
656
657              endif
658
659              if ((.not.flotmy(i,j)).and.(.not.gzmy(i,j))) then
660                 chaly(i,j,k,l)=(btt(i,j-1,k,l)+btt(i,j,k,l))*ffy(i,j,l) !&
661!                      *ro*g*e(k)**(glen(l)+1)/cp(i,j,k)
662
663              else if (gzmy(i,j)) then  ! ice streams
664                  chaly(i,j,k,l)=0.
665
666              else                      ! ice shelves
667                  chaly(i,j,k,l)=0.
668              endif
669
670           end do
671        end do
672     end do
673  end do
674
675
676 
677
678
679!         nouvelle formulation de chaldef_maj(i,j,k), le 4 vient des moyennes
680!         Btt et gauche et droite (ou haut et bas) mais il faut sommer
681!         les productions x et y
682!         ancienne formulation CHAL=(RO*G*H(I,J))**4*(SX2+SY2)*(SX*SX+SY*SY)
683
684do k=2,nz
685   do j=1,ny-1
686      do i=1,nx-1
687
688! modif christophe mars 2000 : chalx et chaly sont a 4 dim
689              chaldef_maj(i,j,k)= 0.
690
691! chalk_2 pour ice shelves et ice streams
692              chalk_2=(chal2_x(i,j,k)+chal2_y(i,j,k) + &
693                   chal2_z(i,j,k)+chal2_xy(i,j,k))/cp(i,j,k)
694
695! chalk_1 pour la partie posée
696              do l=n1poly,n2poly
697
698! on somme la chaleur due aux diverses lois de déformation et celles en x et y
699                 chalk_1=(chalx(i,j,k,l)+chalx(i+1,j,k,l))+ &
700                      (chaly(i,j,k,l)+chaly(i,j+1,k,l))
701
702                 chalk_1=ro*g*chalk_1/4.*(e(k)**(glen(l)+1))/cp(i,j,k)
703                 chaldef_maj(i,j,k)= chalk_1 + chaldef_maj(i,j,k) 
704              enddo
705
706! Pour shelves et streams,  on ajoute chalk_2
707                chaldef_maj(i,j,k) = chaldef_maj(i,j,k) + chalk_2
708             end do
709          end do
710       end do
711
712! chaleur produite a la base par le glissement
713
714  do j=2,ny
715     do i=2,nx
716       
717        if (gzmx(i,j)) then
718           chalglissx(i,j)= abs(uxbar(i,j)*tobmx(i,j))
719        else
720           chalglissx(i,j)=ddbx(i,j)*sdx(i,j)**2*ro*g*Hmx(i,j)
721        endif
722
723        if (gzmy(i,j)) then
724           chalglissy(i,j)= abs(uybar(i,j)*tobmy(i,j))
725        else
726           chalglissy(i,j)=ddby(i,j)*sdy(i,j)**2*ro*g*Hmy(i,j)
727        endif
728
729     end do
730  end do
731
732
733!  BOUNDARY CONDITION ICE-ROCK interface
734
735k=nz
736
737
738do j=2,ny-1
739   do i=2,nx-1
740
741! l'ancien chalbed correspond a chaldef_maj(i,j,nz)  : pas la peine de recalculer
742!----------------------------------------------------
743!!$      chalbed=0.
744!!$      do l=n1poly,n2poly
745!!$         chalbed_1=(chalx(i,j,k,l)+chalx(i+1,j,k,l)) &
746!!$              + (chaly(i,j,k,l)+chaly(i,j+1,k,l))
747!!$         chalbed_1=ro*g*chalbed_1/4.*h(i,j)
748!!$         chalbed= chalbed_1 +chalbed
749!!$      enddo
750
751
752!        rajouter un flux de chaleur pour la production par deformation
753!        dans la derniere 1/2 maille et par le glissement
754!        attention phid est >0 et ghf est <0
755           if (icouple.ge.1.and..not.flot(i,j)) then
756
757!           phid avec fonte sous les streams
758!              phid(i,j)=0.25*(chalglissx(i,j)+chalglissx(i+1,j)+ &
759!                        chalglissy(i,j)+chalglissy(i,j+1)) + &
760!                        chalbed*fracq
761
762
763! moyenne phid sur 4 points : attention 0.5 pour moyenne gauche - droite
764! mais les chaleurs x,y s'ajoutent
765
766!              phid(i,j)=0.5*(chalglissx(i,j)+chalglissx(i+1,j)+ &
767!                        chalglissy(i,j)+chalglissy(i,j+1))
768
769! moyenne phid sur 12 points. On moyenne chaque chaleur et on somme x et y
770
771             phid(i,j)=0.25*(chalglissx(i,j)+chalglissx(i+1,j)+ &
772                        chalglissy(i,j)+chalglissy(i,j+1)) 
773             phid(i,j)=phid(i,j)+0.125*( &
774                  ((chalglissx(i,j-1)+chalglissx(i+1,j-1))+ &
775                  (chalglissx(i,j+1)+chalglissx(i+1,j+1))) + &
776                  ((chalglissy(i-1,j)+chalglissy(i-1,j+1)) +  &
777                  (chalglissy(i+1,j)+chalglissy(i+1,j+1))))
778                       
779! plus la base est loin du point de fusion moins la chaleur de glissement est prise en compte
780              phid(i,j)=phid(i,j)*exp((T(i,j,nz)-tpmp(i,j,nz))*ecart_phid)
781
782! pour sorties Heino
783              chalgliss_maj(i,j)=phid(i,j)
784
785              phid(i,j)=phid(i,j)+ chaldef_maj(i,j,nz)*fracq*H(i,j)*cp(i,j,nz)
786
787
788           else
789              phid(i,j)=0.
790           endif
791        end do
792     end do
793return
794end subroutine Q_prod_demi
795
796!-------------------------------------------------------------------------------------
797
798! calcul avec la somme des carres
799
800subroutine Q_prod_centre
801
802  do k=2,nz 
803     do j=2,ny-1
804        do i=2,nx-1
805
806
807! calcul de la chaleur de deformation selon xx yy zz et xy
808! pour les ice-streams et ice-shelves
809! les divers eps sont calculés dans strain-rate pour l'ensemble de la grille
810
811        if ((flot(i,j).or.flgzmx(i,j).or.flgzmx(i+1,j)).or.  &
812                        (flgzmy(i,j).or.flgzmy(i,j+1))) then
813
814              chal2_x(i,j,k)=2.*visc(i,j,k)*epsxx(i,j)**2
815              chal2_y(i,j,k)=2.*visc(i,j,k)*epsyy(i,j)**2
816              chal2_z(i,j,k)=2.*visc(i,j,k)*(-epsxx(i,j)-epsyy(i,j))**2
817
818!  epsxy est calcule sur les noeuds mineur 1/2,1/2, faire la moyenne
819              chal2_xy(i,j,k)=(epsxy(i,j)+epsxy(i+1,j)+epsxy(i+1,j+1)+epsxy(i,j+1))
820              chal2_xy(i,j,k)=2.*visc(i,j,k)*(chal2_xy(i,j,k)*0.25)**2
821
822        else   ! glace posée
823              chal2_x(i,j,k)=0.
824              chal2_y(i,j,k)=0.
825              chal2_z(i,j,k)=0.
826              chal2_xy(i,j,k)=0.
827        endif
828        end do
829     end do
830  end do
831
832!  partie SIA   calcul de la chaleur produite sur chaque demi maille
833  do l=n1poly,n2poly
834     do j=2,ny
835        do i=2,nx
836
837!          ffx a 3 dimensions !
838           ffx(i,j,l)=ddx(i,j,l)*sdx(i,j)*sdx(i,j)
839           ffy(i,j,l)=ddy(i,j,l)*sdy(i,j)*sdy(i,j)
840        end do
841     end do
842  end do
843
844
845
846  do l=n1poly,n2poly
847     do k=2,nz
848        do j=2,ny
849           do i=2,nx
850              if ((.not.flotmx(i,j)).and.(.not.gzmx(i,j))) then
851                 chalx(i,j,k,l)=(btt(i-1,j,k,l)+btt(i,j,k,l))*ffx(i,j,l) !&
852!                      *ro*g*e(k)**(glen(l)+1)/cp(i,j,k)
853
854              else if (gzmx(i,j)) then ! ice streams
855                  chalx(i,j,k,l)=0.
856
857              else                     ! ice shelves 
858                  chalx(i,j,k,l)=0.
859
860              endif
861
862              if ((.not.flotmy(i,j)).and.(.not.gzmy(i,j))) then
863                 chaly(i,j,k,l)=(btt(i,j-1,k,l)+btt(i,j,k,l))*ffy(i,j,l) !&
864!                      *ro*g*e(k)**(glen(l)+1)/cp(i,j,k)
865
866              else if (gzmy(i,j)) then  ! ice streams
867                  chaly(i,j,k,l)=0.
868
869              else                      ! ice shelves
870                  chaly(i,j,k,l)=0.
871              endif
872
873           end do
874        end do
875     end do
876  end do
877
878
879 
880
881
882!         nouvelle formulation de chaldef_maj(i,j,k), le 4 vient des moyennes
883!         Btt et gauche et droite (ou haut et bas) mais il faut sommer
884!         les productions x et y
885!         ancienne formulation CHAL=(RO*G*H(I,J))**4*(SX2+SY2)*(SX*SX+SY*SY)
886
887do k=2,nz
888   do j=1,ny-1
889      do i=1,nx-1
890
891! modif christophe mars 2000 : chalx et chaly sont a 4 dim
892              chaldef_maj(i,j,k)= 0.
893
894! chalk_2 pour ice shelves et ice streams
895              chalk_2=(chal2_x(i,j,k)+chal2_y(i,j,k) + &
896                   chal2_z(i,j,k)+chal2_xy(i,j,k))/cp(i,j,k)
897
898! chalk_1 pour la partie posée
899              do l=n1poly,n2poly
900
901! on somme la chaleur due aux diverses lois de déformation et celles en x et y
902! la somme se fait par (Cx2+Cy2)** 0.5. Le 4 vient des moyennes Btt + des moyennes gauche-droite
903
904
905                 chalk_1=(chalx(i,j,k,l)+chalx(i+1,j,k,l))**2+ &
906                      (chaly(i,j,k,l)+chaly(i,j+1,k,l))**2
907                 chalk_1=chalk_1**0.5
908
909                 chalk_1=ro*g*chalk_1/4.*(e(k)**(glen(l)+1))/cp(i,j,k)
910                 chaldef_maj(i,j,k)= chalk_1 + chaldef_maj(i,j,k) 
911              enddo
912
913! Pour shelves et streams,  on ajoute chalk_2
914                chaldef_maj(i,j,k) = chaldef_maj(i,j,k) + chalk_2
915             end do
916          end do
917       end do
918
919! chaleur produite a la base par le glissement
920
921  do j=2,ny
922     do i=2,nx
923       
924        if (gzmx(i,j)) then
925           chalglissx(i,j)= abs(uxbar(i,j)*tobmx(i,j))
926        else
927           chalglissx(i,j)=ddbx(i,j)*sdx(i,j)**2*ro*g*Hmx(i,j)
928        endif
929
930        if (gzmy(i,j)) then
931           chalglissy(i,j)= abs(uybar(i,j)*tobmy(i,j))
932        else
933           chalglissy(i,j)=ddby(i,j)*sdy(i,j)**2*ro*g*Hmy(i,j)
934        endif
935
936     end do
937  end do
938
939
940!  BOUNDARY CONDITION ICE-ROCK interface
941
942k=nz
943
944
945do j=2,ny-1
946   do i=2,nx-1
947
948
949
950!        rajouter un flux de chaleur pour la production par deformation
951!        dans la derniere 1/2 maille et par le glissement
952!        attention phid est >0 et ghf est <0
953           if (icouple.ge.1.and..not.flot(i,j)) then
954
955!           phid avec fonte sous les streams
956!              phid(i,j)=0.25*(chalglissx(i,j)+chalglissx(i+1,j)+ &
957!                        chalglissy(i,j)+chalglissy(i,j+1)) + &
958!                        chalbed*fracq
959
960
961! moyenne phid sur 4 points formulation (A2+B2)**0.5
962!
963              chalgliss_maj(i,j)=(chalglissx(i,j)+chalglissx(i+1,j))**2+  &
964                   (chalglissy(i,j)+chalglissy(i,j+1))**2
965
966              chalgliss_maj(i,j)=chalgliss_maj(i,j)**0.5
967
968              chalgliss_maj(i,j)=chalgliss_maj(i,j)*0.5 ! les moyennes droite gauche
969
970! plus la base est loin du point de fusion moins la chaleur de glissement est prise en compte
971              chalgliss_maj(i,j)=chalgliss_maj(i,j)*exp((T(i,j,nz)-tpmp(i,j,nz))*ecart_phid)
972
973
974! flux total a rajouter à la base
975              phid(i,j)=chalgliss_maj(i,j)+chaldef_maj(i,j,nz)*fracq*H(i,j)*cp(i,j,nz)
976
977
978           else
979              phid(i,j)=0.
980           endif
981        end do
982     end do
983return
984end subroutine Q_prod_centre
985
986!-------------------------------------------------------------------------------------
987
988
989! calcul pour essayer d'avoir la chaleur en alpha4
990
991subroutine Q_prod_pente
992
993implicit none
994real, dimension(nx,ny) :: pente2_maj   ! pente au carre sur les noeuds majeurs
995
996do j=2,ny
997   do i=2,nx
998      pente2_maj(i,j)=(sdx(i,j)+sdx(i+1,j))**2 + &
999                      (sdy(i,j)+sdy(i,j+1))**2
1000
1001      pente2_maj(i,j)=pente2_maj(i,j)*0.25    !pour la moyenne sur sdx et sdy
1002   end do
1003end do
1004
1005  do k=2,nz 
1006     do j=2,ny-1
1007        do i=2,nx-1
1008
1009
1010! calcul de la chaleur de deformation selon xx yy zz et xy
1011! pour les ice-streams et ice-shelves
1012! les divers eps sont calculés dans strain-rate pour l'ensemble de la grille
1013
1014        if ((flot(i,j).or.flgzmx(i,j).or.flgzmx(i+1,j)).or.  &
1015                        (flgzmy(i,j).or.flgzmy(i,j+1))) then
1016
1017              chal2_x(i,j,k)=2.*visc(i,j,k)*epsxx(i,j)**2
1018              chal2_y(i,j,k)=2.*visc(i,j,k)*epsyy(i,j)**2
1019              chal2_z(i,j,k)=2.*visc(i,j,k)*(-epsxx(i,j)-epsyy(i,j))**2
1020
1021!  epsxy est calcule sur les noeuds mineur 1/2,1/2, faire la moyenne
1022              chal2_xy(i,j,k)=(epsxy(i,j)+epsxy(i+1,j)+epsxy(i+1,j+1)+epsxy(i,j+1))
1023              chal2_xy(i,j,k)=2.*visc(i,j,k)*(chal2_xy(i,j,k)*0.25)**2
1024
1025        else   ! glace posée
1026              chal2_x(i,j,k)=0.
1027              chal2_y(i,j,k)=0.
1028              chal2_z(i,j,k)=0.
1029              chal2_xy(i,j,k)=0.
1030        endif
1031        end do
1032     end do
1033  end do
1034
1035!  partie SIA   calcul de la chaleur produite sur chaque demi maille
1036!!$  do l=n1poly,n2poly
1037!!$     do j=2,ny
1038!!$        do i=2,nx
1039!!$
1040!!$!          ffx a 3 dimensions ! 
1041!!$           ffx(i,j,l)=ddx(i,j,l)*sdx(i,j)*sdx(i,j)
1042!!$           ffy(i,j,l)=ddy(i,j,l)*sdy(i,j)*sdy(i,j)
1043!!$        end do
1044!!$     end do
1045!!$  end do
1046
1047
1048
1049  do l=n1poly,n2poly
1050     do k=2,nz
1051        do j=2,ny
1052           do i=2,nx
1053              if ((.not.flotmx(i,j)).and.(.not.gzmx(i,j))) then
1054                 chalx(i,j,k,l)=(btt(i-1,j,k,l)+btt(i,j,k,l))*ddx(i,j,l)
1055
1056              else if (gzmx(i,j)) then ! ice streams
1057                  chalx(i,j,k,l)=0.
1058
1059              else                     ! ice shelves 
1060                  chalx(i,j,k,l)=0.
1061
1062              endif
1063
1064              if ((.not.flotmy(i,j)).and.(.not.gzmy(i,j))) then
1065                 chaly(i,j,k,l)=(btt(i,j-1,k,l)+btt(i,j,k,l))*ddy(i,j,l)
1066
1067              else if (gzmy(i,j)) then  ! ice streams
1068                  chaly(i,j,k,l)=0.
1069
1070              else                      ! ice shelves
1071                  chaly(i,j,k,l)=0.
1072              endif
1073
1074           end do
1075        end do
1076     end do
1077  end do
1078
1079
1080 
1081
1082
1083!         nouvelle formulation de chaldef_maj(i,j,k), le 4 vient des moyennes
1084!         Btt et gauche et droite (ou haut et bas) mais il faut sommer
1085!         les productions x et y
1086!         ancienne formulation CHAL=(RO*G*H(I,J))**4*(SX2+SY2)*(SX*SX+SY*SY)
1087
1088do k=2,nz
1089   do j=1,ny-1
1090      do i=1,nx-1
1091
1092! modif christophe mars 2000 : chalx et chaly sont a 4 dim
1093              chaldef_maj(i,j,k)= 0.
1094
1095! chalk_2 pour ice shelves et ice streams
1096              chalk_2=(chal2_x(i,j,k)+chal2_y(i,j,k) + &
1097                   chal2_z(i,j,k)+chal2_xy(i,j,k))/cp(i,j,k)
1098
1099! chalk_1 pour la partie posée
1100              do l=n1poly,n2poly
1101
1102! on somme la chaleur due aux diverses lois de déformation et celles en x et y
1103! la somme se fait par (Cx2+Cy2)** 0.5. Le 4 vient des moyennes Btt + des moyennes gauche-droite
1104
1105! on fait la moyenne des termes ddx*btt (*0.25 pour cette moyenne, le 0.5 est pour les btt)
1106                 chalk_1=(chalx(i,j,k,l)+chalx(i+1,j,k,l))+ &
1107                      (chaly(i,j,k,l)+chaly(i,j+1,k,l))
1108
1109                 chalk_1=chalk_1*0.25*0.5
1110
1111! on multiplie par la pente moyenne au carre sur le noeud majeur
1112                 chalk_1=chalk_1*pente2_maj(i,j)
1113
1114                 chalk_1=ro*g*chalk_1*(e(k)**(glen(l)+1))/cp(i,j,k)  ! attention plus de /4
1115                 chaldef_maj(i,j,k)= chalk_1 + chaldef_maj(i,j,k) 
1116              enddo
1117
1118! Pour shelves et streams,  on ajoute chalk_2
1119                chaldef_maj(i,j,k) = chaldef_maj(i,j,k) + chalk_2
1120             end do
1121          end do
1122       end do
1123
1124! chaleur produite a la base par le glissement
1125
1126  do j=2,ny
1127     do i=2,nx
1128       
1129        if (gzmx(i,j)) then
1130           chalglissx(i,j)= abs(uxbar(i,j)*tobmx(i,j))
1131        else
1132           chalglissx(i,j)=ddbx(i,j)*sdx(i,j)**2*ro*g*Hmx(i,j)
1133        endif
1134
1135        if (gzmy(i,j)) then
1136           chalglissy(i,j)= abs(uybar(i,j)*tobmy(i,j))
1137        else
1138           chalglissy(i,j)=ddby(i,j)*sdy(i,j)**2*ro*g*Hmy(i,j)
1139        endif
1140
1141     end do
1142  end do
1143
1144
1145!  BOUNDARY CONDITION ICE-ROCK interface
1146
1147k=nz
1148
1149
1150do j=2,ny-1
1151   do i=2,nx-1
1152
1153
1154
1155!        rajouter un flux de chaleur pour la production par deformation
1156!        dans la derniere 1/2 maille et par le glissement
1157!        attention phid est >0 et ghf est <0
1158           if (icouple.ge.1.and..not.flot(i,j)) then
1159
1160!           phid avec fonte sous les streams
1161!              phid(i,j)=0.25*(chalglissx(i,j)+chalglissx(i+1,j)+ &
1162!                        chalglissy(i,j)+chalglissy(i,j+1)) + &
1163!                        chalbed*fracq
1164
1165
1166! moyenne phid sur 4 points formulation (A2+B2)**0.5  ! je garde pour l'instant a revoir
1167!
1168              chalgliss_maj(i,j)=(chalglissx(i,j)+chalglissx(i+1,j))**2+  &
1169                   (chalglissy(i,j)+chalglissy(i,j+1))**2
1170
1171              chalgliss_maj(i,j)=chalgliss_maj(i,j)**0.5
1172
1173              chalgliss_maj(i,j)=chalgliss_maj(i,j)*0.5 ! les moyennes droite gauche
1174
1175! plus la base est loin du point de fusion moins la chaleur de glissement est prise en compte
1176              chalgliss_maj(i,j)=chalgliss_maj(i,j)*exp((T(i,j,nz)-tpmp(i,j,nz))*ecart_phid)
1177
1178
1179! flux total a rajouter à la base
1180
1181              phid(i,j)=chalgliss_maj(i,j)+chaldef_maj(i,j,nz)*fracq*H(i,j)*cp(i,j,nz)
1182
1183
1184           else
1185              phid(i,j)=0.
1186           endif
1187        end do
1188     end do
1189return
1190end subroutine Q_prod_pente
1191
1192!-------------------------------------------------------------------------------------
1193
1194! routine qui prend la chaleur remise au socle
1195
1196subroutine Q_U_taub
1197
1198implicit none
1199real, dimension(nx,ny) :: pente_maj   ! pente au carre sur les noeuds majeurs
1200real,dimension(nx,ny) :: vit_maj      ! Ubar moyennee sur les noeuds majeurs
1201real,dimension(nx,ny) :: uslid_maj      ! Uslid moyennee sur les noeuds majeurs
1202
1203do j=2,ny
1204   do i=2,nx
1205
1206
1207      pente_maj(i,j)=(sdx(i,j)+sdx(i+1,j))**2 + &          ! pente
1208                      (sdy(i,j)+sdy(i,j+1))**2
1209
1210      pente_maj(i,j)=pente_maj(i,j)*0.25    ! 0.25pour la moyenne sur sdx et sdy
1211      pente_maj(i,j)=pente_maj(i,j)**0.5
1212
1213      vit_maj(i,j)=(uxbar(i,j)+uxbar(i+1,j))**2 + &        ! vitesse de bilan
1214                      (uybar(i,j)+uybar(i,j+1))**2
1215      vit_maj(i,j)=vit_maj(i,j)*0.25
1216      vit_maj(i,j)=vit_maj(i,j)**0.5
1217
1218      uslid_maj(i,j)=(ux(i,j,nz)+ux(i+1,j,nz))**2 + &        ! vitesse de bilan
1219                      (uy(i,j,nz)+uy(i,j+1,nz))**2
1220      uslid_maj(i,j)=uslid_maj(i,j)*0.25
1221      uslid_maj(i,j)=uslid_maj(i,j)**0.5
1222
1223   end do
1224end do
1225
1226
1227
1228chaldef_maj(:,:,:)=0.
1229chalgliss_maj(:,:)=ro*g*H(:,:)*pente_maj(:,:)*uslid_maj(:,:)
1230chaldef_maj(:,:,nz)=ro*g*H(:,:)*pente_maj(:,:)*vit_maj(:,:)-chalgliss_maj(:,:)
1231
1232! plus la base est loin du point de fusion moins la chaleur de glissement est prise en compte
1233
1234if (ispinup.eq.0) then
1235   chalgliss_maj(:,:)=chalgliss_maj(:,:)*exp((T(:,:,nz)-tpmp(:,:,nz))*ecart_phid)
1236end if
1237
1238! flux total a rajouter à la base
1239phid(:,:)=chaldef_maj(:,:,nz)+chalgliss_maj(:,:)
1240
1241
1242if (ispinup.eq.3) then
1243   phid(:,:)=ro*g*H(:,:)*pente_maj(:,:)*vit_maj(:,:)
1244   chaldef_maj(:,:,nz)=0
1245end if
1246
1247
1248return
1249end subroutine Q_U_taub
1250
1251!----------------------------------------------------------------------------------
1252! routine qui prend la chaleur remise au socle mais sur les mailles staggered
1253! les noeuds stag sont ceux au milieux des mailles vitesses
1254
1255subroutine Q_U_taub_stag
1256
1257implicit none
1258real, dimension(nx,ny) :: pente_stag   ! pente au carre sur les noeuds stag
1259real,dimension(nx,ny) :: vit_stag      ! Ubar moyennee sur les noeuds stag
1260real,dimension(nx,ny) :: uslid_stag    ! Uslid moyennee sur les noeuds stag
1261real,dimension(nx,ny) :: Hmxy          ! Epaisseur moyenne sur les noeuds stag
1262real,dimension(nx,ny) :: Qslid         ! chaleur produite par glissement sur les noeuds stag
1263real,dimension(nx,ny) :: Qdef          ! chaleur produite par deformation sur les noeuds stag
1264
1265
1266do j=2,ny
1267   do i=2,nx
1268
1269
1270      pente_stag(i,j)=(sdx(i,j)+sdx(i,j-1))**2 + &          ! pente
1271                      (sdy(i,j)+sdy(i-1,j))**2
1272
1273      pente_stag(i,j)=pente_stag(i,j)*0.25    ! 0.25pour la moyenne sur sdx et sdy
1274      pente_stag(i,j)=pente_stag(i,j)**0.5
1275
1276      vit_stag(i,j)=(uxbar(i,j)+uxbar(i,j-1))**2 + &        ! vitesse de bilan
1277                      (uybar(i,j)+uybar(i-1,j))**2
1278      vit_stag(i,j)=vit_stag(i,j)*0.25
1279      vit_stag(i,j)=vit_stag(i,j)**0.5
1280
1281      uslid_stag(i,j)=(ux(i,j,nz)+ux(i,j-1,nz))**2 + &        ! vitesse de bilan
1282                      (uy(i,j,nz)+uy(i-1,j,nz))**2
1283      uslid_stag(i,j)=uslid_stag(i,j)*0.25
1284      uslid_stag(i,j)=uslid_stag(i,j)**0.5
1285
1286      Hmxy(i,j)=((H(i,j)+H(i-1,j-1))+(H(i,j-1)+H(i-1,j)))*0.25
1287
1288   end do
1289end do
1290
1291
1292Qslid(:,:)=ro*g*Hmxy(:,:)*pente_stag(:,:)*uslid_stag(:,:)
1293Qdef(:,:)=ro*g*Hmxy(:,:)*pente_stag(:,:)*vit_stag(:,:)-Qslid(:,:)
1294chaldef_maj(:,:,:)=0.
1295
1296
1297! On fait la moyenne de la chaleur produite sur les mailles stag
1298
1299do j=2,ny-1
1300   do i=2,nx-1
1301      chalgliss_maj(i,j)=(Qslid(i,j)+Qslid(i+1,j+1))+(Qslid(i,j+1)+Qslid(i+1,j))
1302      chalgliss_maj(i,j)=chalgliss_maj(i,j)*0.25
1303
1304      chaldef_maj(i,j,nz)=(Qdef(i,j)+Qdef(i+1,j+1))+(Qdef(i,j+1)+Qdef(i+1,j))
1305      chaldef_maj(i,j,nz)=chaldef_maj(i,j,nz)*0.25
1306
1307! plus la base est loin du point de fusion moins la chaleur de glissement est prise en compte
1308      chalgliss_maj(i,j)=chalgliss_maj(i,j)*exp((T(i,j,nz)-tpmp(i,j,nz))*ecart_phid)
1309
1310   end do
1311end do
1312
1313!!$
1314!!$write(126,*)'time=',time
1315!!$j=41
1316!!$do i=20,30
1317!!$   write(126,'(i3,6(e14.4,1x))') i,chalgliss_maj(i,j),chaldef_maj(i,j,nz),Qdef(i,j),Qdef(i+1,j+1) &
1318!!$       , Qdef(i,j+1),Qdef(i+1,j)
1319!!$end do
1320
1321! flux total a rajouter à la base
1322phid(:,:)=chaldef_maj(:,:,nz)+chalgliss_maj(:,:)
1323
1324return
1325end subroutine Q_U_taub_stag
1326
1327!----------------------------------------------------------------------------------
1328! routine qui prend la chaleur remise au socle mais sur les mailles staggered
1329! les noeuds stag sont ceux au milieux des mailles vitesses
1330
1331subroutine Q_SIA_stag
1332
1333implicit none
1334real,dimension(nx,ny) :: pente_stag   ! pente au carre sur les noeuds stag
1335real,dimension(nx,ny,nz) :: vit_stag  ! Ubar moyennee sur les noeuds stag
1336real,dimension(nx,ny) :: uslid_stag   ! Uslid moyennee sur les noeuds stag
1337real,dimension(nx,ny) :: Hmxy         ! Epaisseur moyenne sur les noeuds stag
1338real,dimension(nx,ny) :: Qslid        ! chaleur produite par glissement sur les noeuds stag
1339real,dimension(nx,ny,nz) :: Qdef      ! chaleur produite par deformation sur les noeuds stag
1340
1341pente_stag = 0.
1342vit_stag   = 0.
1343uslid_stag = 0.
1344Hmxy  = 0.
1345Qslid = 0.
1346Qdef = 0.
1347
1348do j=2,ny
1349   do i=2,nx
1350
1351!     variables 2D
1352      pente_stag(i,j)=(sdx(i,j)+sdx(i,j-1))**2 + &          ! pente
1353                      (sdy(i,j)+sdy(i-1,j))**2
1354
1355      pente_stag(i,j)=pente_stag(i,j)*0.25    ! 0.25pour la moyenne sur sdx et sdy
1356      pente_stag(i,j)=pente_stag(i,j)**0.5
1357
1358      uslid_stag(i,j)=(ux(i,j,nz)+ux(i,j-1,nz))**2 + &        ! vitesse de bilan
1359                      (uy(i,j,nz)+uy(i-1,j,nz))**2
1360      uslid_stag(i,j)=uslid_stag(i,j)*0.25
1361      uslid_stag(i,j)=uslid_stag(i,j)**0.5
1362
1363      Hmxy(i,j)=((H(i,j)+H(i-1,j-1))+(H(i,j-1)+H(i-1,j)))*0.25
1364
1365!     en vertical : calcul de la deformation par differentiation des vitesses
1366
1367      do k=1,nz
1368
1369         vit_stag(i,j,k)=(ux(i,j,k)+ux(i,j-1,k))**2 + &        ! magnitude vitesse
1370                         (uy(i,j,k)+uy(i-1,j,k))**2            ! a tous niveaux k
1371
1372         vit_stag(i,j,k)=vit_stag(i,j,k)*0.25
1373         vit_stag(i,j,k)=vit_stag(i,j,k)**0.5
1374      end do
1375
1376      Qdef(i,j,1)=0.
1377      do k=2,nz-1
1378         Qdef(i,j,k)=vit_stag(i,j,k-1)-vit_stag(i,j,k+1)            ! difference des vitesses
1379         Qdef(i,j,k)=ro*g*pente_stag(i,j)*e(k)*Qdef(i,j,k)/2./de    ! gamma tau
1380      end do
1381
1382
1383      Qdef(i,j,nz)=vit_stag(i,j,nz-1)-vit_stag(i,j,nz)               ! pour le fond differentiation
1384      Qdef(i,j,nz)=ro*g*pente_stag(i,j)*Qdef(i,j,k)/de               ! sur une seule maille
1385   end do
1386end do
1387
1388
1389Qslid(:,:)=ro*g*Hmxy(:,:)*pente_stag(:,:)*uslid_stag(:,:)
1390
1391
1392
1393! On fait la moyenne de la chaleur produite sur les mailles stag
1394
1395do j=2,ny-1
1396   do i=2,nx-1
1397         chalgliss_maj(i,j)=(Qslid(i,j)+Qslid(i+1,j+1))+(Qslid(i,j+1)+Qslid(i+1,j))
1398         chalgliss_maj(i,j)=chalgliss_maj(i,j)*0.25
1399
1400! plus la base est loin du point de fusion moins la chaleur de glissement est prise en compte
1401         chalgliss_maj(i,j)=chalgliss_maj(i,j)*exp((T(i,j,nz)-tpmp(i,j,nz))*ecart_phid)
1402
1403! chaleur due a la defromation
1404      do k=1,nz
1405         chaldef_maj(i,j,k)=(Qdef(i,j,k)+Qdef(i+1,j+1,k))+(Qdef(i,j+1,k)+Qdef(i+1,j,k))
1406         chaldef_maj(i,j,k)=chaldef_maj(i,j,k)*0.25
1407         chaldef_maj(i,j,k)=chaldef_maj(i,j,k)/cp(i,j,k)
1408      end do
1409
1410
1411   end do
1412end do
1413
1414
1415!!$write(126,*)'time=',time
1416!!$j=41
1417!!$do i=20,30
1418!!$   write(126,'(i3,6(e14.4,1x))') i,chalgliss_maj(i,j),chaldef_maj(i,j,nz),Qdef(i,j),Qdef(i+1,j+1) &
1419!!$       , Qdef(i,j+1),Qdef(i+1,j)
1420!!$end do
1421
1422! flux total a rajouter à la base
1423
1424phid(:,:)=chalgliss_maj(:,:)+chaldef_maj(:,:,nz)*fracq*H(:,:)*cp(:,:,nz)
1425
1426return
1427end subroutine Q_SIA_stag
1428
1429!----------------------------------------------------------------------------------
1430! routine qui prend la chaleur remise au socle mais sur les mailles staggered
1431! les noeuds stag sont ceux au milieux des mailles vitesses
1432! le calcul de taub tient compte des noeuds grzmx
1433
1434subroutine Q_all_stag
1435
1436implicit none
1437
1438real, dimension(nx,ny) :: pente_stag   ! pente au carre sur les noeuds stag
1439real, dimension(nx,ny) :: tob_stag   ! contrainte basale
1440real, dimension(nx,ny) :: tox       ! contraintes sur maille mx
1441real, dimension(nx,ny) :: toy       ! contraintes sur maille mx
1442real,dimension(nx,ny,nz) :: vit_stag      ! Ubar moyennee sur les noeuds stag
1443real,dimension(nx,ny) :: uslid_stag    ! Uslid moyennee sur les noeuds stag
1444real,dimension(nx,ny) :: Hmxy          ! Epaisseur moyenne sur les noeuds stag
1445real,dimension(nx,ny) :: Qslid         ! chaleur produite par glissement sur les noeuds stag
1446real,dimension(nx,ny,nz) :: Qdef          ! chaleur produite par deformation sur les noeuds stag
1447
1448
1449where (flgzmx(:,:))
1450   tox(:,:)=tobmx(:,:)
1451elsewhere
1452   tox(:,:)=sdx(:,:)*hmx(:,:)*ro*g
1453end where
1454
1455where (flgzmy(:,:))
1456   toy(:,:)=tobmy(:,:)
1457elsewhere
1458   toy(:,:)=sdy(:,:)*hmy(:,:)*ro*g
1459end where
1460
1461
1462do j=2,ny
1463   do i=2,nx
1464
1465!     variables 2D
1466      pente_stag(i,j)=(sdx(i,j)+sdx(i,j-1))**2 + &          ! pente
1467                      (sdy(i,j)+sdy(i-1,j))**2
1468
1469      pente_stag(i,j)=pente_stag(i,j)*0.25    ! 0.25pour la moyenne sur sdx et sdy
1470      pente_stag(i,j)=pente_stag(i,j)**0.5
1471
1472
1473      tob_stag(i,j)=(tox(i,j)+tox(i,j-1))**2 + &          ! pente
1474                      (toy(i,j)+toy(i-1,j))**2
1475
1476      tob_stag(i,j)=tob_stag(i,j)*0.25    ! 0.25pour la moyenne sur sdx et sdy
1477      tob_stag(i,j)=tob_stag(i,j)**0.5
1478
1479      uslid_stag(i,j)=(ux(i,j,nz)+ux(i,j-1,nz))**2 + &        ! vitesse de bilan
1480                      (uy(i,j,nz)+uy(i-1,j,nz))**2
1481      uslid_stag(i,j)=uslid_stag(i,j)*0.25
1482      uslid_stag(i,j)=uslid_stag(i,j)**0.5
1483
1484      Hmxy(i,j)=((H(i,j)+H(i-1,j-1))+(H(i,j-1)+H(i-1,j)))*0.25
1485
1486!     en vertical : calcul de la deformation par differentiation des vitesses
1487
1488      do k=1,nz
1489
1490         vit_stag(i,j,k)=(ux(i,j,k)+ux(i,j-1,k))**2 + &        ! magnitude vitesse
1491                         (uy(i,j,k)+uy(i-1,j,k))**2            ! a tous niveaux k
1492
1493         vit_stag(i,j,k)=vit_stag(i,j,k)*0.25
1494         vit_stag(i,j,k)=vit_stag(i,j,k)**0.5
1495      end do
1496
1497      Qdef(i,j,1)=0.
1498      do k=2,nz-1
1499         Qdef(i,j,k)=vit_stag(i,j,k-1)-vit_stag(i,j,k+1)            ! difference des vitesses
1500         Qdef(i,j,k)=ro*g*pente_stag(i,j)*e(k)*Qdef(i,j,k)/2./de    ! gamma tau
1501      end do
1502
1503
1504      Qdef(i,j,nz)=vit_stag(i,j,nz-1)-vit_stag(i,j,nz)               ! pour le fond differentiation
1505      Qdef(i,j,nz)=ro*g*pente_stag(i,j)*Qdef(i,j,k)/de               ! sur une seule maille
1506   end do
1507end do
1508
1509
1510Qslid(:,:)=tob_stag*uslid_stag(:,:)
1511
1512
1513
1514! On fait la moyenne de la chaleur produite sur les mailles stag
1515
1516do j=2,ny-1
1517   do i=2,nx-1
1518         chalgliss_maj(i,j)=(Qslid(i,j)+Qslid(i+1,j+1))+(Qslid(i,j+1)+Qslid(i+1,j))
1519         chalgliss_maj(i,j)=chalgliss_maj(i,j)*0.25
1520
1521! plus la base est loin du point de fusion moins la chaleur de glissement est prise en compte
1522         chalgliss_maj(i,j)=chalgliss_maj(i,j)*exp((T(i,j,nz)-tpmp(i,j,nz))*ecart_phid)
1523
1524! chaleur due a la defromation
1525      do k=1,nz         
1526         chaldef_maj(i,j,k)=(Qdef(i,j,k)+Qdef(i+1,j+1,k))+(Qdef(i,j+1,k)+Qdef(i+1,j,k))
1527         chaldef_maj(i,j,k)=chaldef_maj(i,j,k)*0.25
1528         chaldef_maj(i,j,k)=chaldef_maj(i,j,k)/cp(i,j,k)
1529      end do
1530
1531
1532   end do
1533end do
1534
1535!!$
1536!!$write(126,*)'time=',time
1537!!$j=41
1538!!$do i=20,30
1539!!$   write(126,'(i3,6(e14.4,1x))') i,chalgliss_maj(i,j),chaldef_maj(i,j,nz),Qdef(i,j),Qdef(i+1,j+1) &
1540!!$       , Qdef(i,j+1),Qdef(i+1,j)
1541!!$end do
1542
1543! flux total a rajouter à la base
1544phid(:,:)=chalgliss_maj(:,:)+chaldef_maj(:,:,nz)*fracq*H(:,:)*cp(:,:,nz)
1545 
1546return
1547end subroutine Q_all_stag
1548
1549
1550
1551end module ice_temp_mod
Note: See TracBrowser for help on using the repository browser.