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

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

Cleaning-up: unused variables removed

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