!> \file eaubasale-0.5_mod.f90 !! TOOOOOOOOOOO DOOOOOOOOOOOOO !< !> \namespace eau_basale !! TOOOOOOOOOOO DOOOOOOOOOOOOO !! \author ... !! \date ... !! @note Used modules !! @note - use module3D_phy !! @note - use eau_basale !! @note - use param_phy_mod !! @note - use relaxation_waterdif_mod !< module eau_basale use module3d_phy use param_phy_mod use relaxation_waterdif_mod ! module qui contient la routine relaxation ! pour l'eau avec interface explicite implicit none ! dimensionnements !------------------------------------------------------------------------ LOGICAL :: ecoulement_eau REAL :: KONDMAX real :: kond0 REAL :: INFILTR REAL :: hmax_till !< épaisseur de la couche de till REAL :: hmax_wat !< épaisseur de la couche d'eau dans le till REAL :: poro_till !< porosité du till REAL :: Zflot !< hauteur d'eau donnant la flottaison real :: keffmax !< kondmax*hmax_wat real :: nefflocal REAL,dimension(NX,NY) :: limit_hw !< conditions aux limites integer,dimension(NX,NY) :: klimit !< ou appliquer les conditions REAL,dimension(NX,NY) :: pot_w !< pour calculer le gradient de pression REAL,dimension(NX,NY) :: pot_f !< pour les points flottants REAL,dimension(NX,NY) :: hw !< hauteur d'eau dans le till, limite a hmax_wat !< c'est la hauteur sur laquelle peut se faire !< l'ecoulement de l'eau REAL,dimension(nx,ny) :: keff !< conductivite effective : Kond*hw REAL,dimension(NX,NY) :: bmelt_w !< fusion (terme source) exprimé en m d'eau REAL,dimension(NX,NY) :: vieuxhwater !< valeur de hwater au debut de l'appel INTEGER :: NXlocal,NYlocal !< pour passer NX et NY en parametre a la routine relaxation contains !> SUBROUTINE: init_eaubasale !!Initialisation du block eaubasale !> subroutine init_eaubasale namelist/eaubasale1/ecoulement_eau,hwatermax,infiltr namelist/param_hydr/hmax_till,poro_till,kond0 ! formats pour les ecritures dans 42 428 format(A) ! lecture des parametres du run block eaubasale1 !-------------------------------------------------------------------- rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,eaubasale1) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&eaubasale1 ! nom du premier bloc eau basale ' write(num_rep_42,*) write(num_rep_42,*) 'ecoulement_eau = ',ecoulement_eau write(num_rep_42,*) 'hwatermax = ',hwatermax write(num_rep_42,*) 'infiltr = ', infiltr write(num_rep_42,*)'/' write(num_rep_42,428) '! ecoulement eau : .false. -> modele bucket, sinon equ. diffusion' write(num_rep_42,428) '! hwatermax : hauteur d eau basale maximum dans le sediment (m)' write(num_rep_42,428) '! infiltr est la quantite d eau qui peut s infiltrer dans le sol (m/an)' write(num_rep_42,*) !valeurs numeriques des parametres hydrauliques rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,param_hydr) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '¶m_hydr ! nom du bloc parametres hydrauliques ' write(num_rep_42,*) write(num_rep_42,*) 'hmax_till = ',hmax_till write(num_rep_42,*) 'poro_till = ',poro_till write(num_rep_42,*) 'kond0 = ',kond0 write(num_rep_42,*)'/' write(num_rep_42,428) '! hmax_till (m) : epaisseur max du sediment ' write(num_rep_42,428) '! poro_till : porosite du sediment ' write(num_rep_42,428) '! conductivite du sediment : kond0 (m/s)' write(num_rep_42,*) hmax_wat=hmax_till*poro_till ! hauteur maxi que peut atteindre l'eau dans la couche de till ! Conductivite hydraulique : cond passée en m/an ( car le dt est en années) kond(:,:)=kond0 kond(:,:) = kond(:,:)*SECYEAR kondmax = 1.*SECYEAR hdotwater(:,:)=0. NXlocal=NX NYlocal=NY end subroutine init_eaubasale !> SUBROUTINE: eaubasale !! to do !> subroutine eaubasale !(pwater) version correspondant à la thèse de Vincent ! A terme, il faudrait en faire un module !$ USE OMP_LIB !~ integer :: t1,t2,ir !~ real :: temps, t_cpu_0, t_cpu_1, t_cpu, norme !~ !~ ! Temps CPU de calcul initial. !~ call cpu_time(t_cpu_0) !~ ! Temps elapsed de reference. !~ call system_clock(count=t1, count_rate=ir) if (itracebug.eq.1) call tracebug(' Entree dans routine eau basale') !$OMP PARALLEL !$OMP WORKSHARE vieuxhwater(:,:) = hwater(:,:) kond(:,:) = kond0*SECYEAR ! conditions aux limites klimit(:,:)=0 limit_hw(:,:)=-9999. !$OMP END WORKSHARE !$OMP DO do j=1,ny do i=1,nx if (flot(i,j))then ! points flottants klimit(i,j)=1 limit_hw(i,j)=(sealevel-Bsoc(i,j))*rowg/rofreshg else if (IBASE(I,J).eq.1) then ! base froide klimit(i,j)=1 limit_hw(i,j)=0. else if ((.not.flot(i,j)).and.(H(i,j).lt.1.)) then ! bord de la calotte klimit(i,j)=1 limit_hw(i,j)=10. ! riviere de 10 m de profondeur endif end do end do !$OMP END DO !$OMP DO ! do j=2,ny-1 do j=1,ny do i=1,nx bmelt_w(i,j)=bmelt(i,j)*rofresh/ro hwater(i,j)=max(hwater(i,j),0.) hw(i,j)=min(hwater(i,j),hmax_wat) enddo enddo !$OMP END DO !$OMP END PARALLEL ecoul: if (ecoulement_eau) then ! print*,'===dans eaubasale t, dt',time,dt,'kond 1.0e-5' ! write(6,*) 'entree ecoulement_eau' !$OMP PARALLEL !$OMP DO do j=1,ny do i=1,nx ! calcul des potentiels pot_w(I,J)=rofreshg*B(I,J) ! potentiel de gravite ! on ajoute pression glace mais pas la pression d'eau qui est traitée comme diffusion pot_w(I,J)=pot_w(I,J)+rog*H(I,J) pot_f(I,J)=rowg*(sealevel-S(i,j)+H(I,J)) ! pression a la base de l'ice shelf enddo enddo !$OMP END DO ! sorties debug 17 juillet 2007 debug_3D(:,:,5)=pot_w(:,:) debug_3D(:,:,6)=pot_f(:,:) debug_3D(:,:,7)=pot_w(:,:)+rofreshg*hwater(:,:) debug_3D(:,:,8)=hwater(:,:) !$OMP DO do j=2,ny do i=2,nx if (H(i,j).gt.25.) then ! calcul du gradient de pression if (flotmx(i,j)) then pgx(i,j)=(pot_f(i-1,j)-pot_f(i,j))/dx else pgx(i,j)=(pot_w(i-1,j)-pot_w(i,j))/dx endif if (flotmy(i,j)) then pgy(i,j)=(pot_f(i,j-1)-pot_f(i,j))/dy else pgy(i,j)=(pot_w(i,j-1)-pot_w(i,j))/dy endif endif pgx(i,j)=pgx(i,j)/rofreshg ! pour passer pgx sans unité pgy(i,j)=pgy(i,j)/rofreshg end do end do !$OMP END DO !$OMP END PARALLEL if (nt.gt.0) then if (dt.gt.0.) then !!!!!!!!!!!!!!!!!!!!!!relax_water si dt>0 !------------------------------------------------------------------------- ! les points de la grounding line ont une conductivité hydraulique élevée ! même si la base est froide. modif catritz 18 janvier 2005 !open(166,file='erreur_eau') !$OMP PARALLEL !$OMP DO do j=2,Ny-1 do i=2,Nx-1 ! cond=0 si glace froide et pas sur la grounding line if ((IBASE(i,j).eq.1).and. & (.not.flot(i,j-1)).and.(.not.flot(i,j+1)).and. & (.not.flot(i-1,j)).and.(.not.flot(i+1,j))) KOND(i,j)=0.! 1.0e-5 ! cond infinie quand epaisseur faible et glace flottante if (flot(i,j).or.H(i,j).le.1.5) kond(i,j)= kondmax ! conductivite forte lorsque N est faible (croit à partir de 100 bars) nefflocal=0.91*H(i,j)-hwater(i,j) nefflocal=max(100.,nefflocal) if (nefflocal.le.1000.) then kond(i,j)=kond(i,j)*1000./nefflocal endif kond(i,j)=min(kondmax,kond(i,j)) ! conductivite effective (conductivité * taille du tuyau en m2/an) keff(i,j)=kond(i,j)*hw(i,j) end do end do !$OMP END DO ! condition sur les bords de la grille !$OMP WORKSHARE kond(1,:)=kondmax kond(nx,:)=kondmax kond(:,1)=kondmax kond(:,ny)=kondmax vieuxhwater(:,:) = hwater(:,:) !$OMP END WORKSHARE !$OMP END PARALLEL !!$OMP SINGLE call relaxation_waterdif(nxlocal,nylocal,dt,dx,vieuxhwater,limit_hw,klimit,bmelt_w,infiltr,pgx,pgy,keff,hwater) !!$OMP END SINGLE else ! print*,'dt=',dt,'pas de relax_water' endif!!!!!!!!!!!!!!!!!! fin relax_water si dt>0 endif ! Apres relaxation, boundary conditions, extremum values ! ================, ===================, =============== ! ------------variation d'epaisseur entre 2 pas de temps ------------ ! on le fait avant les butoirs pour justement pouvoir les estimer !$OMP PARALLEL if (dt.gt.0.) then ! print*,'dt=',dt,'pas de relax_water' !$OMP DO do j=1,ny do i=1,nx hdotwater(i,j)=(hwater(i,j)-vieuxhwater(i,j))/dt end do end do !$OMP END DO endif !$OMP DO PRIVATE(Zflot) do i=1,nx do j=1,ny ! ______________ valeurs de hwater et pwater _____________________________ ! if (flot(I,J).or.H(I,J).le.1.5) then ! we are outside of the ice sheet if (flot(i,j)) then ! if flot > hwater=hwater in ocean hwater(i,j)= sealevel - bsoc(i,j) ! hwater(i,j)= max(0.,hwater(i,j)) if (hwater(i,j).lt.0.) hwater(i,j)=0. pwater(i,j)= hwater(i,j)*rowg else hwater(i,j) = 0. ! bare grounded land > no hwater pwater(i,j)=0. endif else ! sous la calotte ---------------------- ! Attention le bloc suivant est pour le run 20 ! Zflot=row/ro*(sealevel-Bsoc(i,j))-10. Zflot=H(i,j)*rog/rofreshg if (hwater(i,j).le.0.) then hwater(i,j)=0. else if (hwater(i,j).gt.zflot) then hwater(i,j)=zflot hw(i,j)=min(hwater(i,j),hmax_wat) pwater(i,j)=rofreshg*Hwater(i,j) endif ! if ((i.eq.60).and.(j.eq.60)) write(6,*) hw(i,j),hwater(i,j) ! sinon ! hw(i,j)=min(hw(i,j),hmax_wat) ! Hwater(i,j)=hw(i,j) ! pwater(i,j)=rog*H(I,J)*(hw(I,J)/hmax_wat)**(3.5) endif ! bloc qui pourrait servir pour mettre l'eau encore plus sous pression ! ----------------------------------------------------------------------------- ! if (HWATER(i,j).gt.poro_till*hmax_till) then ! pwater(i,j)=pwater(i,j)+(HWATER(i,j)-poro_till*hmax_till)/(compress_w*hmax_till) ! endif end do end do !$OMP END DO ! ************ valeurs de HWATER pour les coins ******** hwater(1,1)=(hwater(1,2)+hwater(2,1))/2. hwater(1,ny)=(hwater(1,ny-1)+hwater(2,ny))/2. hwater(nx,1)=(hwater(nx,2)+hwater(nx-1,1))/2. hwater(nx,ny)=(hwater(nx,ny-1)+hwater(nx-1,ny))/2. ! pour les sorties de flux d'eau !$OMP DO do j=2,ny do i=2,nx if (keff(i,j)==0..or.keff(i-1,j)==0.) then phiwx(i,j)=0. ! to avoid division by o else phiwx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx) phiwx(i,j)=phiwx(i,j)*2*(keff(i,j)*keff(i-1,j))/(keff(i,j)+keff(i-1,j)) endif ! ligne pour sortir les pgx pgx(i,j)=(pgx(i,j)+(hwater(i-1,j)-hwater(i,j))/dx) end do end do !$OMP END DO !$OMP DO do j=2,ny do i=2,nx if (keff(i,j)==0..or.keff(i,j-1)==0.) then phiWy(i,j)=0. ! to avoid division by o else phiWy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx) phiWy(i,j)=phiWy(i,j)*2*(keff(i,j)*keff(i,j-1))/(keff(i,j)+keff(i,j-1)) endif pgy(i,j)=(pgy(i,j)+(hwater(i,j-1)-hwater(i,j))/dx) enddo enddo !$OMP END DO !$OMP END PARALLEL else ! partie originellement dans icetemp à mettre dans un autre module ! (système module choix) hauteur d'eau cumulee if (ISYNCHRO.eq.1) then !$OMP PARALLEL !$OMP DO do j=1,ny do i=1,nx if (.not.flot(i,j)) then hwater(i,j)=hwater(i,j)+(bmelt(i,j)*dtt) hwater(i,j)=hwater(i,j)-(infiltr*dtt) hwater(i,j)=max(hwater(i,j),0.) hwater(i,j)=min(hwater(i,j),hwatermax) else hwater(i,j)=hwatermax endif ! if (IBASE(I,J).eq.1) HWATER(I,J)=0. end do end do !$OMP END DO !$OMP END PARALLEL end if endif ecoul !~ ! Temps elapsed final !~ call system_clock(count=t2, count_rate=ir) !~ temps=real(t2 - t1,kind=4)/real(ir,kind=4) !~ ! Temps CPU de calcul final !~ call cpu_time(t_cpu_1) !~ t_cpu = t_cpu_1 - t_cpu_0 !~ ! Impression du resultat. !~ print '(//,3X,"Valeurs de nx et ny : ",I5,I5/, & !~ & 3X,"Temps elapsed : ",1PE10.3," sec.",/, & !~ & 3X,"Temps CPU : ",1PE10.3," sec.",/, & !~ & 3X,"Norme (PB si /= 0) : ",1PE10.3,//)', & !~ nx,ny,temps,t_cpu,norme return end subroutine eaubasale end module eau_basale