Context Navigation

source: trunk/SOURCES/relaxation_water_diffusion.f90 @ 25

Last change on this file since 25 was 4, checked in by dumas, 10 years ago
initial import GRISLI trunk
File size: 6.1 KB

Line
1	!> \file relaxation_water_diffusion.f90
2	!! Module pour la resolution de l'equation de relaxation et de l equation de diffusion.
3	!<
4
5	!> \namespace relaxation_waterdif_mod
6	!! Module pour la resolution de l'equation de relaxation et de l equation de diffusion .
7	!! \author ...
8	!! \date ...
9	!! @note dhwat/dt = bmelt-infiltr-d/dx(Konddhwat/dx)+d/dx(Kondpgx)
10	!<
11	module relaxation_waterdif_mod
12
13
14	CONTAINS
15
16
17	subroutine relaxation_waterdif(NXX,NYY,DT,DX,vieuxHWATER,limit_hw,klimit,BMELT,INFILTR,PGMX,PGMY,KOND,KONDMAX,HWATER)
18	implicit none
19
20
21	! declaration des variables en entree
22	!------------------------------------------------
23	INTEGER, intent(in) :: NXX, NYY !< defini la taille des tableaux
24	REAL, intent(in) :: DT !< pas de temps court
25	REAL, intent(in) :: DX !< pas en x
26	REAL, intent(in) :: INFILTR !< basal infiltration (lose of water)
27	REAL, intent(in) :: KONDMAX !< maximum hydaulic conductivity (outside ice sheet)
28
29	REAL,dimension(NXX,NYY), intent(in) :: limit_hw !< conditions aux limites
30	integer,dimension(NXX,NYY), intent(in) :: klimit !< ou appliquer les conditions
31	REAL,dimension(NXX,NYY), intent(in) :: vieuxHWATER !< H au pas de temps precedent 'o'
32	REAL,dimension(NXX,NYY), intent(in) :: BMELT !< basal melting 'o'
33	REAL,dimension(NXX,NYY), intent(in) :: PGMX !< hydaulic potential gratient '>'
34	REAL,dimension(NXX,NYY), intent(in) :: PGMY !< hydaulic potential gratient '^'
35	REAL,dimension(NXX,NYY), intent(in) :: KOND !< hydaulic conductivity 'o'
36
37	! declaration des variables en sortie
38	!-------------------------------------
39	REAL,dimension(NXX,NYY), intent(out):: HWATER !< basal water thickness 'o' (pressure equivalent)
40
41
42	! declaration des variables locales
43	!----------------------------------
44	INTEGER :: I,J
45	REAL :: TESTH
46	REAL,dimension(NXX,NYY) :: ARELAX,BRELAX,CRELAX,DRELAX,ERELAX,FRELAX
47	REAL,dimension(NXX,NYY) :: DELTAH
48	REAL :: RESTE,DELH,VH
49	INTEGER :: ntour
50	INTEGER :: mbord
51	REAL :: DTSRGDX,dtwdx2
52	LOGICAL :: STOPP
53	REAL,dimension(NXX,NYY) :: KMX, KMY
54	! REAL :: RHO,RHOW,RHOG !,SECYEAR!! ice density, water density, density*acceleration
55	! write(6,*) 'entree relaxation'
56
57	! write(166,*)' entree relaxation waterdif'
58	HWATER(:,:)= vieuxHWATER(:,:)
59
60	! calcul de kmx et kmx a partir de KOND
61	! conductivite hyrdraulique sur les noeuds mineurs
62	! moyenne harmonique
63	! ----------------------------------------
64
65	do j=2,nyy
66	do i=2,nxx
67
68	if ((kond(i,j).lt.1.e-20).or.(kond(i-1,j).lt.1.e-20)) then
69	kmx(i,j)=0. ! to avoid division by o
70	else
71	kmx(i,j)=2(kond(i,j)kond(i-1,j))/(kond(i,j)+kond(i-1,j))
72	endif
73
74	end do
75	end do
76
77	do j=2,nyy
78	do i=2,nxx
79	if ((kond(i,j).lt.1.e-20).or.(kond(i,j-1).lt.1.e-20)) then
80	kmy(i,j)=0. ! to avoid division by o
81	else
82	kmy(i,j)=2(kond(i,j)kond(i,j-1))/(kond(i,j)+kond(i,j-1))
83	endif
84
85	enddo
86	enddo
87
88
89	! attribution des coefficients arelax ....
90	! ----------------------------------------
91	! SECYEAR=365.24.3600.
92	! rho=910.
93	! rhow=1000.
94	! rhog=rhow*9.81
95	! dtsrgdx=dt/(rhog*DX) a mon avis c'est rhow qu'il fallait utiliser. Maintenant cette
96	! division est faite dans eaubasale
97
98	dtsrgdx=dt/DX
99	dtwdx2=dt/dx/dx
100
101	arelax(:,:)=0.
102	brelax(:,:)=0.
103	crelax(:,:)=1.
104	drelax(:,:)=0.
105	erelax(:,:)=0.
106	frelax(:,:)=limit_hw(:,:)
107
108
109
110	do J=2,NYY-1
111	do I=2,NXX-1
112
113	if (klimit(i,j).eq.0) then
114
115	! calcul du vecteur
116
117	FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT
118	frelax(i,j)=frelax(i,j)+(kmx(i,j)pgmx(i,j)-kmx(i+1,j)pgmx(i+1,j))*dtsrgdx
119	frelax(i,j)=frelax(i,j)+(kmy(i,j)pgmy(i,j)-kmy(i,j+1)pgmy(i,j+1))*dtsrgdx
120
121	! calcul des diagonales
122	arelax(i,j)=-kmx(i,j)*dtwdx2 ! arelax : diagonale i-1,j
123
124	brelax(i,j)=-kmx(i+1,j)*dtwdx2 ! brelax : diagonale i+1,j
125
126	drelax(i,j)=-kmy(i,j)*dtwdx2 ! drelax : diagonale i,j-1
127
128	erelax(i,j)=-kmy(i,j+1)*dtwdx2 ! drelax : diagonale i,j+1
129
130	crelax(i,j)=1.+((kmx(i,j)+kmx(i+1,j))+(kmy(i,j+1)+kmy(i,j+1)))*dtwdx2
131	!crelax : diagonale i,j
132
133	else if (klimit(i,j).eq.1) then
134	hwater(i,j)=limit_hw(i,j)
135	! write(6,*) i,j,hwater(i,j),crelax(i,j),frelax(i,j),arelax(i,j)
136
137	endif
138
139	end do
140	end do
141
142
143
144	! Boucle de relaxation :
145	! ----------------------
146	STOPP = .false.
147
148	ntour=0
149	stopp=.false.
150	Do WHILE(.NOT.STOPP)
151	ntour=ntour+1
152	! write(6,*) 'boucle de relaxation numero',ntour
153
154	do j=2,NYY-1
155	do i=2,NXX-1
156
157	RESTE = (((ARELAX(I,J)HWATER(I-1,J) + BRELAX(I,J)HWATER(I+1,J)) &
158	+ (DRELAX(I,J)HWATER(I,J-1) + ERELAX(I,J)HWATER(I,J+1))) &
159	+ CRELAX(I,J)*HWATER(I,J))- FRELAX(I,J)
160
161	DELTAH(I,J) = RESTE/CRELAX(I,J)
162
163	end do
164	end do
165
166	deltah(:,:)=min(deltah(:,:),10.)
167	deltah(:,:)=max(deltah(:,:),-10.)
168
169
170
171	! il faut faire le calcul suivant dans une autre boucle car RESTE est fonction
172	! de hwater sur les points voisins.
173	do j=2,NYY-1
174	do i=2,NXX-1
175	HWATER(I,J) = HWATER(I,J) - DELTAH(I,J)
176
177	end do
178	end do
179
180
181	! critere d'arret:
182
183
184
185	Delh=0
186	Vh=0
187
188	DO j=2,NYY-1
189	DO i=2,NXX-1
190
191	! write(166,*) I,J,delh,deltah(i,j)
192	Delh=Delh+deltah(i,j)**2
193	! Vh=Vh+h(i,j)**2.
194	END DO
195	END DO
196
197	! write(6,*) delh,maxval(deltah),minval(deltah)
198	! testh=SQRT(Delh/Vh)
199	if (delh.gt.0.) then
200	testh=SQRT(Delh)/((NXX-2)*(NYY-2))
201	else
202	testh=0.
203	endif
204	STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000)
205	! write(6,*) ntour,testh
206
207
208	362 continue
209	end do
210	end subroutine relaxation_waterdif
211	end module relaxation_waterdif_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats: