source: trunk/SOURCES/Hudson_files/eaubasale-0.5_hudson_mod.f90 @ 237

Last change on this file since 237 was 237, checked in by aquiquet, 6 years ago

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