source: trunk/SOURCES/eaubasale-0.5_mod.f90 @ 23

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

initial import GRISLI trunk

File size: 13.1 KB
Line 
1!> \file eaubasale-0.5_mod.f90
2!! TOOOOOOOOOOO DOOOOOOOOOOOOO
3!<
4
5!> \namespace eau_basale
6!! TOOOOOOOOOOO DOOOOOOOOOOOOO
7!! \author ...
8!! \date ...
9!! @note Used modules
10!! @note   - use module3D_phy
11!! @note   - use eau_basale
12!! @note   - use param_phy_mod
13!! @note   - use relaxation_waterdif_mod
14!<
15module eau_basale
16
17use module3d_phy
18use param_phy_mod
19use relaxation_waterdif_mod ! module qui contient la routine relaxation
20                           ! pour l'eau avec interface explicite
21
22
23implicit none
24
25! dimensionnements
26!------------------------------------------------------------------------
27    LOGICAL :: ecoulement_eau
28
29    REAL :: KONDMAX
30    real :: kond0
31    REAL :: INFILTR
32    REAL :: DTWAT               !< pas de temps pour l'eau
33    REAL :: compress_w
34    REAL :: hmax_till    !< épaisseur de la couche de till
35    REAL :: hmax_wat     !< épaisseur de la couche d'eau dans le till
36    REAL :: poro_till    !< porosité du till
37    REAL :: Zflot        !< hauteur d'eau donnant la flottaison
38    real :: keffmax      !< kondmax*hmax_wat
39    real :: nefflocal
40
41    REAL,dimension(NX,NY) :: limit_hw    !< conditions aux limites
42    integer,dimension(NX,NY) :: klimit    !< ou appliquer les conditions
43    REAL,dimension(NX,NY) :: pot_w             !< pour calculer le gradient de pression
44    REAL,dimension(NX,NY) :: pot_f             !< pour les points flottants
45    REAL,dimension(NX,NY) :: hw                !< hauteur d'eau dans le till, limite a hmax_wat
46                                                 !< c'est la hauteur sur laquelle peut se faire
47                                                 !< l'ecoulement de l'eau
48    REAL,dimension(nx,ny) :: keff              !< conductivite effective : Kond*hw
49
50
51    REAL,dimension(NX,NY) :: bmelt_w          !< fusion (terme source) exprimé en m d'eau
52    REAL,dimension(NX,NY) :: vieuxhwater      !< valeur de hwater au debut de l'appel
53    REAL,dimension(NX,NY) :: tetar            !< pour le routage de l'eau dans l'ocean
54
55
56    INTEGER :: NXlocal,NYlocal !< pour passer NX et NY en parametre a la routine relaxation
57
58
59contains
60!> SUBROUTINE: init_eaubasale
61!!Initialisation du block eaubasale
62!>
63subroutine init_eaubasale
64
65namelist/eaubasale1/ecoulement_eau,hwatermax,infiltr 
66namelist/param_hydr/hmax_till,poro_till,kond0
67
68! formats pour les ecritures dans 42
69428 format(A)
70
71! lecture des parametres du run                      block eaubasale1
72!--------------------------------------------------------------------
73
74
75rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
76read(num_param,eaubasale1)
77write(num_rep_42,428)'!___________________________________________________________' 
78write(num_rep_42,428) '&eaubasale1              ! nom du premier bloc eau basale '
79write(num_rep_42,*)
80write(num_rep_42,*) 'ecoulement_eau = ',ecoulement_eau
81write(num_rep_42,*) 'hwatermax      = ',hwatermax
82write(num_rep_42,*) 'infiltr        = ', infiltr
83write(num_rep_42,*)'/'                           
84write(num_rep_42,428) '! ecoulement eau : .false. -> modele bucket, sinon equ. diffusion'
85write(num_rep_42,428) '! hwatermax :  hauteur d eau basale maximum dans le sediment (m)'
86write(num_rep_42,428) '! infiltr est la quantite d eau qui peut s infiltrer dans le sol (m/an)'
87write(num_rep_42,*)
88
89
90!valeurs numeriques des parametres hydrauliques
91
92rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
93read(num_param,param_hydr)
94
95write(num_rep_42,428)'!___________________________________________________________' 
96write(num_rep_42,428) '&param_hydr             ! nom du  bloc parametres hydrauliques '
97write(num_rep_42,*)
98write(num_rep_42,*) 'hmax_till      = ',hmax_till
99write(num_rep_42,*) 'poro_till      = ',poro_till
100write(num_rep_42,*) 'kond0          = ',kond0 
101write(num_rep_42,*)'/'                           
102write(num_rep_42,428) '! hmax_till (m) : epaisseur max du sediment '
103write(num_rep_42,428) '! poro_till : porosite du sediment '
104write(num_rep_42,428) '! conductivite du sediment :  kond0 (m/s)'
105write(num_rep_42,*)
106
107
108
109
110hmax_wat=hmax_till*poro_till ! hauteur maxi que peut atteindre l'eau dans la couche de till
111
112
113
114kond(:,:)=kond0
115
116
117! Conductivite hydraulique : cond passée en m/an ( car le dt est en années)
118kond(:,:)=kond0
119kond(:,:) = kond(:,:)*SECYEAR 
120kondmax   = 1.*SECYEAR 
121keffmax=kondmax*hmax_wat
122hdotwater(:,:)=0.
123
124
125NXlocal=NX
126NYlocal=NY
127
128end subroutine init_eaubasale
129
130!> SUBROUTINE: eaubasale
131!! to do
132!>
133subroutine eaubasale !(pwater)   version correspondant à la thèse de Vincent
134  ! A terme, il faudrait en faire un module
135
136
137
138
139
140
141  if (itracebug.eq.1)  call tracebug(' Entree dans routine eau basale')
142
143  vieuxhwater(:,:) = hwater(:,:)
144  kond(:,:) = kond0*SECYEAR 
145
146
147  ! conditions aux limites
148  klimit(:,:)=0
149  limit_hw(:,:)=-9999.
150  do j=1,ny
151     do i=1,nx
152
153
154        if (flot(i,j))then  ! points flottants
155           klimit(i,j)=1
156           limit_hw(i,j)=(sealevel-Bsoc(i,j))*rowg/rofreshg
157
158        else if (IBASE(I,J).eq.1) then ! base froide
159           klimit(i,j)=1
160           limit_hw(i,j)=0.
161
162        else if ((.not.flot(i,j)).and.(H(i,j).lt.1.)) then  ! bord de la calotte
163           klimit(i,j)=1
164           limit_hw(i,j)=10.    ! riviere de 10 m de profondeur
165        endif
166
167     end do
168  end do
169
170
171
172  do i=1,nx
173     do j=2,ny-1
174        bmelt_w(i,j)=bmelt(i,j)*rofresh/ro   
175        hwater(i,j)=max(hwater(i,j),0.)
176        hw(i,j)=min(hwater(i,j),hmax_wat) 
177
178     enddo
179  enddo
180
181
182  ecoul:  if (ecoulement_eau) then
183     !  print*,'===dans eaubasale t, dt',time,dt,'kond 1.0e-5'
184     !   write(6,*) 'entree ecoulement_eau'
185
186     do J=1,NY
187        do I=1,NX
188           tetar(i,j)=(xlong(i,j))*PI/180. ! pourrait etre fait une fois pour toute
189           PGX(I,J)=101*sin(tetar(i,j))*1.e-2                 
190           PGY(I,J)=101*cos(tetar(i,j))*1.e-2                 
191
192           !   calcul des potentiels
193           pot_w(I,J)=rofreshg*B(I,J)  ! potentiel de gravite
194
195           ! on ajoute pression glace mais pas la pression d'eau qui est traitée comme diffusion
196
197           pot_w(I,J)=pot_w(I,J)+rog*H(I,J)
198
199           pot_f(I,J)=rowg*(sealevel-S(i,j)+H(I,J))  ! pression a la base de l'ice shelf
200        enddo
201     enddo
202
203     ! sorties debug 17 juillet 2007
204     debug_3D(:,:,5)=pot_w(:,:)
205     debug_3D(:,:,6)=pot_f(:,:)
206     debug_3D(:,:,7)=pot_w(:,:)+rofreshg*hwater(:,:)
207     debug_3D(:,:,8)=hwater(:,:)
208
209     do j=2,ny
210        do i=2,nx
211           if (H(i,j).gt.25.) then
212
213              !           calcul du gradient de pression
214              ! _______________________________________
215
216              if (flotmx(i,j)) then
217                 pgx(i,j)=(pot_f(i-1,j)-pot_f(i,j))/dx
218              else
219                 pgx(i,j)=(pot_w(i-1,j)-pot_w(i,j))/dx
220              endif
221
222              if (flotmy(i,j)) then
223                 pgy(i,j)=(pot_f(i,j-1)-pot_f(i,j))/dy
224              else
225                 pgy(i,j)=(pot_w(i,j-1)-pot_w(i,j))/dy
226              endif
227
228
229           endif
230
231
232           pgx(i,j)=pgx(i,j)/rofreshg   ! pour passer pgx sans unité
233           pgy(i,j)=pgy(i,j)/rofreshg     
234        end do
235     end do
236
237
238
239     if (nt.gt.0) then
240        if (dt.gt.0.)  then !!!!!!!!!!!!!!!!!!!!!!relax_water si dt>0
241
242
243           !-------------------------------------------------------------------------
244           ! les points de la grounding line ont une conductivité hydraulique élevée
245           ! même si la base est froide.  modif catritz 18 janvier 2005
246           !open(166,file='erreur_eau')
247
248           do J=2,NY-1
249              do I=2,NX-1
250
251                 ! cond=0 si glace froide et pas sur la grounding line
252
253                 if ((IBASE(I,J).eq.1).and.   &
254                      (.not.flot(i,j-1)).and.(.not.flot(i,j+1)).and.  &
255                      (.not.flot(i-1,j)).and.(.not.flot(i+1,j)))  KOND(I,J)=0.! 1.0e-5
256
257                 ! cond infinie quand epaisseur faible et glace flottante
258                 if (flot(i,j).or.H(i,j).le.1.5)  kond(i,j)= kondmax
259
260                 ! conductivite forte lorsque N est faible (croit à partir de 100 bars)
261                 nefflocal=0.91*H(i,j)-hwater(i,j)
262                 nefflocal=max(100.,nefflocal)
263                 if (nefflocal.le.1000.) then
264                    kond(i,j)=kond(i,j)*1000./nefflocal
265                 endif
266                 kond(i,j)=min(kondmax,kond(i,j)) 
267
268                 ! conductivite effective (conductivité * taille du tuyau en m2/an)
269
270                 keff(i,j)=kond(i,j)*hw(i,j)
271              end do
272           end do
273
274           ! condition sur les bords de la grille
275
276           kond(1,:)=kondmax
277           kond(nx,:)=kondmax
278           kond(:,1)=kondmax
279           kond(:,ny)=kondmax
280
281           ! fin de la modif 18 janvier 2005---------------------------------------
282
283           vieuxhwater(:,:) = hwater(:,:)
284
285
286           call relaxation_waterdif(nxlocal,nylocal,dt,dx,vieuxhwater,limit_hw,klimit,bmelt_w,infiltr,pgx,pgy,keff,keffmax,hwater)
287
288
289
290        else
291           !         print*,'dt=',dt,'pas de relax_water'
292        endif!!!!!!!!!!!!!!!!!! fin relax_water si dt>0
293
294     endif
295
296
297
298
299     ! Apres relaxation, boundary conditions, extremum values
300     ! ================, ===================, ===============
301
302     !   ------------variation d'epaisseur entre 2 pas de temps ------------
303
304     ! on le fait avant les butoirs pour justement pouvoir les estimer
305
306     if (dt.gt.0.) then
307        !         print*,'dt=',dt,'pas de relax_water'
308        do j=1,ny
309           do i=1,nx
310              hdotwater(i,j)=(hwater(i,j)-vieuxhwater(i,j))/dt 
311           end do
312        end do
313     endif
314
315
316     do i=1,nx
317        do j=1,ny
318
319
320           !  ______________ valeurs de hwater et pwater _____________________________
321           !
322           if (flot(I,J).or.H(I,J).le.1.5) then ! we are outside of the ice sheet
323
324              if (flot(i,j)) then  ! if flot > hwater=hwater in ocean
325                 hwater(i,j)= sealevel - bsoc(i,j)
326                 !             hwater(i,j)= max(0.,hwater(i,j))
327                 if (hwater(i,j).lt.0.) hwater(i,j)=0.
328                 pwater(i,j)= hwater(i,j)*rowg
329              else
330                 hwater(i,j) = 0. ! bare grounded land > no hwater
331                 pwater(i,j)=0.
332              endif
333
334           else          !            sous la  calotte ----------------------
335
336              ! Attention le bloc suivant est pour le run 20
337              !           Zflot=row/ro*(sealevel-Bsoc(i,j))-10.
338              Zflot=H(i,j)*rog/rofreshg
339
340              if (hwater(i,j).le.0.) then
341                 hwater(i,j)=0.
342
343              else if (hwater(i,j).gt.zflot) then
344                 hwater(i,j)=zflot
345                 hw(i,j)=min(hwater(i,j),hmax_wat)
346                 pwater(i,j)=rofreshg*Hwater(i,j)
347              endif
348              !            if ((i.eq.60).and.(j.eq.60)) write(6,*) hw(i,j),hwater(i,j)
349
350              ! sinon
351
352              !          hw(i,j)=min(hw(i,j),hmax_wat) 
353              !          Hwater(i,j)=hw(i,j)
354              !          pwater(i,j)=rog*H(I,J)*(hw(I,J)/hmax_wat)**(3.5)
355
356
357
358
359           endif
360
361
362           ! bloc qui pourrait servir pour mettre l'eau encore plus sous pression
363           ! -----------------------------------------------------------------------------
364           !            if (HWATER(i,j).gt.poro_till*hmax_till) then
365           !             pwater(i,j)=pwater(i,j)+(HWATER(i,j)-poro_till*hmax_till)/(compress_w*hmax_till)
366           !            endif
367
368
369
370
371        end do
372     end do
373
374
375     !   ************ valeurs de HWATER pour les coins ********
376
377     hwater(1,1)=(hwater(1,2)+hwater(2,1))/2.
378     hwater(1,ny)=(hwater(1,ny-1)+hwater(2,ny))/2.
379     hwater(nx,1)=(hwater(nx,2)+hwater(nx-1,1))/2.
380     hwater(nx,ny)=(hwater(nx,ny-1)+hwater(nx-1,ny))/2.
381
382     ! pour les sorties de flux d'eau
383     do j=2,ny
384        do i=2,nx
385           if  (keff(i,j)==0..or.keff(i-1,j)==0.) then
386              phiwx(i,j)=0. ! to avoid division by o             
387           else
388              phiwx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx)
389
390              phiwx(i,j)=phiwx(i,j)*2*(keff(i,j)*keff(i-1,j))/(keff(i,j)+keff(i-1,j))
391           endif
392
393           ! ligne pour sortir les pgx
394
395           pgx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx)
396        end do
397     end do
398
399     do j=2,ny
400        do i=2,nx
401           if  (keff(i,j)==0..or.keff(i,j-1)==0.) then
402              phiWy(i,j)=0. ! to avoid division by o
403           else
404              phiWy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx)
405              phiWy(i,j)=phiWy(i,j)*2*(keff(i,j)*keff(i,j-1))/(keff(i,j)+keff(i,j-1))
406           endif
407           pgy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx)
408
409        enddo
410     enddo
411
412
413
414  else  ! partie originellement dans icetemp à mettre dans un autre module
415     ! (système module choix)     hauteur d'eau cumulee
416
417     if (ISYNCHRO.eq.1) then
418        do j=1,ny
419           do i=1,nx
420              if (.not.flot(i,j)) then
421                 hwater(i,j)=hwater(i,j)+(bmelt(i,j)*dtt)
422                 hwater(i,j)=hwater(i,j)-(infiltr*dtt)
423                 hwater(i,j)=max(hwater(i,j),0.)
424                 hwater(i,j)=min(hwater(i,j),hwatermax)
425              else
426                 hwater(i,j)=hwatermax
427              endif
428              !        if (IBASE(I,J).eq.1) HWATER(I,J)=0.
429           end do
430        end do
431     end if
432  endif ecoul
433
434  return
435end subroutine eaubasale
436 end module eau_basale
Note: See TracBrowser for help on using the repository browser.