source: trunk/SOURCES/relaxation_water_mod-0.4.f90 @ 37

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

initial import GRISLI trunk
File size: 7.8 KB
Line 
1!> \file relaxation_water_mod-0.4.f90
2!! Module pour le calcule de l'ecoulement de l'eau sous la glace 
3!! selon une equation de Darcy avec un shema amont.     
4!<
5
6!> \namespace RELAXATION_WATER_MOD
7!! Module pour le calcule de l'ecoulement de l'eau sous la glace 
8!! selon une equation de Darcy avec un shema amont.   
9!! \author Vincent Peyaud
10!! \date septembre 2004
11!! @note  reprise de relaxation_mod (pour le calcul de l'epaisseur)           
12!!     et adaptation pour le calcul de l'ecoulement de l'eau sous la glace 
13!!     selon une equation de Darcy avec un shema amont.                     
14!<
15
16
17! ************************************************************************ !
18!                                                                          !                                     
19!                       Vincent Peyaud septembre 2004                      !
20!     reprise de relaxation_mod (pour le calcul de l'epaisseur)            !
21!     et adaptation pour le calcul de l'ecoulement de l'eau sous la glace  !
22!     selon une equation de Darcy avec un shema amont.                     !
23!                                                                          !
24! ************************************************************************ !
25
26MODULE RELAXATION_WATER_MOD
27
28CONTAINS
29
30
31  subroutine RELAXATION_WAT(NXX,NYY,DT,DX,vieuxHWATER,BMELT,INFILTR,PGMX,PGMY,KOND,KONDMAX,HWATER)
32    implicit none
33
34
35    ! declaration des variables en entree
36    !------------------------------------------------
37    INTEGER, intent(in) :: NXX, NYY       !< defini la taille des tableaux
38    REAL,    intent(in) ::  DT            !< pas de temps court
39    REAL,    intent(in) ::  DX            !< pas en x
40    REAL,    intent(in) :: INFILTR        !< basal infiltration (lose of water)
41    REAL,    intent(in) :: KONDMAX        !< maximum hydaulic conductivity (outside ice sheet)
42
43    REAL,dimension(NXX,NYY), intent(in) :: vieuxHWATER    !< H au pas de temps precedent 'o'
44    REAL,dimension(NXX,NYY), intent(in) :: BMELT     !< basal melting  'o'
45    REAL,dimension(NXX,NYY), intent(in) :: PGMX     !< hydaulic potential gratient '>'
46    REAL,dimension(NXX,NYY), intent(in) :: PGMY     !< hydaulic potential gratient '^'
47    REAL,dimension(NXX,NYY), intent(in) :: KOND     !< hydaulic conductivity 'o'
48
49    ! declaration des variables en sortie
50    !-------------------------------------
51    REAL,dimension(NXX,NYY), intent(out):: HWATER      !< basal water thickness  'o' (pressure equivalent)
52
53
54    ! declaration des variables locales
55    !----------------------------------
56    INTEGER :: I,J
57    REAL  :: TESTH 
58    REAL,dimension(NXX,NYY) :: ARELAX,BRELAX,CRELAX,DRELAX,ERELAX,FRELAX
59    REAL,dimension(NXX,NYY) :: DELTAH
60    REAL :: RESTE,DELH,VH
61    INTEGER :: ntour
62    INTEGER :: mbord
63    REAL  :: DTSRGDX 
64    LOGICAL :: STOPP
65    REAL,dimension(NXX,NYY) :: KMX, KMY
66    !    REAL  :: RHO,RHOW,RHOG  !,SECYEAR!! ice density, water density, density*acceleration
67
68    HWATER(:,:)= vieuxHWATER(:,:)
69
70    !   calcul de kmx et kmx a partir de KOND
71    !   conductivite hyrdraulique sur les noeuds mineurs
72    !   moyenne harmonique
73    !   ----------------------------------------
74
75    do j=1,nyy
76       do i=2,nxx
77          if  (kond(i,j)==0.or.kond(i-1,j)==0) then
78             kmx(i,j)=0. ! to avoid division by o
79          else
80             kmx(i,j)=2*(kond(i,j)*kond(i-1,j))/(kond(i,j)+kond(i-1,j))
81          endif
82
83       end do
84    end do
85
86    do j=2,nyy
87       do i=2,nxx
88          if  (kond(i,j)==0.or.kond(i,j-1)==0) then
89             kmy(i,j)=0. ! to avoid division by o
90          else
91             kmy(i,j)=2*(kond(i,j)*kond(i,j-1))/(kond(i,j)+kond(i,j-1))
92          endif
93
94       enddo
95    enddo
96
97
98    !   attribution des coefficients  arelax ....
99    !   ----------------------------------------
100    !    SECYEAR=365.*24.*3600.
101    !   rho=910.
102    !    rhow=1000.
103    !    rhog=rhow*9.81
104    !    dtsrgdx=dt/(rhog*DX) a mon avis c'est rhow qu'il fallait utiliser. Maintenant cette
105    !    division est faite dans eaubasale
106
107    dtsrgdx=dt/DX
108
109
110    do J=2,NYY-1
111       do I=2,NXX-1
112
113          FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT !
114
115          !          selon x
116          !                                     ecoulement vers la droite (sens axe x)
117          if ((PGMX(I,J).ge.0.).and.(PGMX(I+1,J).ge.0.)) then
118             ARELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
119             BRELAX(I,J)=0.
120             CRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
121
122             !                                     ecoulement vers la gauche
123          else if ((PGMX(I,J).le.0.).and.(PGMX(I+1,J).le.0.)) then
124             ARELAX(I,J)=0.
125             BRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
126             CRELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
127
128             !                                     ecoulement divergent
129          else if ((PGMX(I,J).le.0.).and.(PGMX(I+1,J).ge.0.)) then
130             ARELAX(I,J)=0.
131             BRELAX(I,J)=0.
132             CRELAX(I,J)=(KMX(I+1,J)*PGMX(I+1,J)-KMX(I,J)*PGMX(I,J))*dtsrgdx
133
134             !                                     ecoulement convergent
135          else if ((PGMX(I,J).ge.0.).and.(PGMX(I+1,J).le.0.)) then
136             ARELAX(I,J)=-KMX(I,J)*PGMX(I,J)*dtsrgdx
137             BRELAX(I,J)=KMX(I+1,J)*PGMX(I+1,J)*dtsrgdx
138             CRELAX(I,J)=0.
139          endif
140 
141
142
143          !          selon y
144          !                                 ecoulement vers le haut (sens axe y)
145          if ((PGMY(I,J).ge.0.).and.(PGMY(I,J+1).ge.0.)) then
146             DRELAX(I,J)=-KMY(I,J)*PGMY(I,J)*dtsrgdx
147             ERELAX(I,J)=0.
148             CRELAX(I,J)=CRELAX(I,J)+1. +(KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx)
149
150             !                                     ecoulement vers le bas
151          else if ((PGMY(I,J).le.0.).and.(PGMY(I,J+1).le.0.)) then
152             DRELAX(I,J)=0.
153             ERELAX(I,J)=KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx
154             CRELAX(I,J)=CRELAX(I,J)+1. -(KMY(I,J)*PGMY(I,J)*dtsrgdx)
155
156             !                                     ecoulement divergent
157          else if ((PGMY(I,J).le.0.).and.(PGMY(I,J+1).ge.0.)) then
158             DRELAX(I,J)=0.
159             ERELAX(I,J)=0.
160             CRELAX(I,J)=CRELAX(I,J)+1.+((KMY(I,J+1)*PGMY(I,J+1)-KMY(I,J)*PGMY(I,J))*dtsrgdx)
161
162             !                                     ecoulement convergent
163          else if ((PGMY(I,J).ge.0.).and.(PGMY(I,J+1).le.0.)) then
164             DRELAX(I,J)=-KMY(I,J)*PGMY(I,J)*dtsrgdx
165             ERELAX(I,J)=KMY(I,J+1)*PGMY(I,J+1)*dtsrgdx
166             CRELAX(I,J)=CRELAX(I,J)+1.
167          endif
168
169       end do
170    end do
171
172    ! Boucle de relaxation :
173    ! ----------------------
174    STOPP = .false.
175
176    ntour=0
177
178    Do 362 WHILE(.NOT.STOPP)
179       ntour=ntour+1
180
181       do j=2,NYY-1
182          do i=2,NXX-1
183
184             RESTE = (((ARELAX(I,J)*HWATER(I-1,J) + BRELAX(I,J)*HWATER(I+1,J)) &
185                  + (DRELAX(I,J)*HWATER(I,J-1) + ERELAX(I,J)*HWATER(I,J+1))) &
186                  + CRELAX(I,J)*HWATER(I,J))- FRELAX(I,J)
187
188             DELTAH(I,J) = RESTE/CRELAX(I,J)             
189
190          end do
191       end do
192
193     
194
195       ! il faut faire le calcul suivant dans une autre boucle car RESTE est fonction
196       ! de hwater sur les points voisins.
197       do j=2,NYY-1
198          do i=2,NXX-1
199             HWATER(I,J) = HWATER(I,J) - DELTAH(I,J)
200          end do
201       end do
202
203       ! critere d'arret:
204
205     
206
207       Delh=0
208       Vh=0
209
210       DO j=2,NYY-1
211          DO i=2,NXX-1
212
213             !   write(6,*) I,J,delh,deltah(i,j)
214             Delh=Delh+deltah(i,j)**2
215             !            Vh=Vh+h(i,j)**2.
216          END DO
217       END DO
218
219!       write(6,*) delh,maxval(deltah),minval(deltah)
220       !      testh=SQRT(Delh/Vh)
221       if (delh.gt.0.) then
222          testh=SQRT(Delh)/((NXX-2)*(NYY-2))
223       else
224          testh=0.
225       endif
226       STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000)
227
228
229362    Continue
230
231     end subroutine relaxation_wat
232
233end module relaxation_water_mod
Note: See TracBrowser for help on using the repository browser.