source: trunk/SOURCES/relaxation_water_diffusion.f90 @ 334

!> \file relaxation_water_diffusion.f90
!! Module pour la resolution de l'equation de relaxation et de l equation de diffusion.
!<

!> \namespace relaxation_waterdif_mod
!! Module pour la resolution de l'equation de relaxation et de l equation de diffusion .
!! \author ...
!! \date ...
!! @note dhwat/dt = bmelt-infiltr-d/dx(Kond*dhwat/dx)+d/dx(Kond*pgx)
!<
module relaxation_waterdif_mod




CONTAINS


subroutine relaxation_waterdif(NXX,NYY,DT,DX,vieuxHWATER,limit_hw,klimit,BMELT,INFILTR,PGMX,PGMY,KOND,HWATER)

  !$ USE OMP_LIB

  implicit none


    ! declaration des variables en entree
    !------------------------------------------------
    INTEGER, intent(in) :: NXX, NYY       !< defini la taille des tableaux
    REAL,    intent(in) ::  DT            !< pas de temps court
    REAL,    intent(in) ::  DX            !< pas en x
    REAL,    intent(in) :: INFILTR        !< basal infiltration (lose of water) 

    REAL,dimension(NXX,NYY), intent(in) :: limit_hw    !< conditions aux limites
    integer,dimension(NXX,NYY), intent(in) :: klimit    !< ou appliquer les conditions 
    REAL,dimension(NXX,NYY), intent(in) :: vieuxHWATER    !< H au pas de temps precedent 'o'
    REAL,dimension(NXX,NYY), intent(in) :: BMELT     !< basal melting  'o'
    REAL,dimension(NXX,NYY), intent(in) :: PGMX     !< hydaulic potential gratient '>' 
    REAL,dimension(NXX,NYY), intent(in) :: PGMY     !< hydaulic potential gratient '^'
    REAL,dimension(NXX,NYY), intent(in) :: KOND     !< hydaulic conductivity 'o'

    ! declaration des variables en sortie
    !-------------------------------------
    REAL,dimension(NXX,NYY), intent(out):: HWATER      !< basal water thickness  'o' (pressure equivalent)


    ! declaration des variables locales
    !----------------------------------
    INTEGER :: I,J
    REAL  :: TESTH 
    REAL,dimension(NXX,NYY) :: ARELAX,BRELAX,CRELAX,DRELAX,ERELAX,FRELAX
    REAL,dimension(NXX,NYY) :: DELTAH
    REAL :: RESTE,DELH,VH
    INTEGER :: ntour
    REAL  :: DTSRGDX,dtwdx2
    LOGICAL :: STOPP
    REAL,dimension(NXX,NYY) :: KMX, KMY
    !    REAL  :: RHO,RHOW,RHOG  !,SECYEAR!! ice density, water density, density*acceleration
!    write(6,*) 'entree relaxation'

! write(166,*)' entree  relaxation waterdif'
!$OMP PARALLEL
!$OMP WORKSHARE
    HWATER(:,:)= vieuxHWATER(:,:)
!$OMP END WORKSHARE
    !   calcul de kmx et kmx a partir de KOND
    !   conductivite hyrdraulique sur les noeuds mineurs
    !   moyenne harmonique
    !   ----------------------------------------
!$OMP DO
    do j=2,nyy
       do i=2,nxx

          if  ((kond(i,j).lt.1.e-20).or.(kond(i-1,j).lt.1.e-20)) then
             kmx(i,j)=0. ! to avoid division by o
          else
             kmx(i,j)=2*(kond(i,j)*kond(i-1,j))/(kond(i,j)+kond(i-1,j))
          endif

       end do
    end do
!$OMP END DO

!$OMP DO
    do j=2,nyy
       do i=2,nxx
          if  ((kond(i,j).lt.1.e-20).or.(kond(i,j-1).lt.1.e-20)) then
             kmy(i,j)=0. ! to avoid division by o
          else
             kmy(i,j)=2*(kond(i,j)*kond(i,j-1))/(kond(i,j)+kond(i,j-1))
          endif

       enddo
    enddo
!$OMP END DO

    !   attribution des coefficients  arelax ....
    !   ----------------------------------------
    !    SECYEAR=365.*24.*3600.
    !   rho=910.
    !    rhow=1000.
    !    rhog=rhow*9.81
    !    dtsrgdx=dt/(rhog*DX) a mon avis c'est rhow qu'il fallait utiliser. Maintenant cette 
    !    division est faite dans eaubasale 

    dtsrgdx=dt/DX
    dtwdx2=dt/dx/dx

!$OMP WORKSHARE
  deltah(:,:)=0.
  arelax(:,:)=0.
  brelax(:,:)=0.
  crelax(:,:)=1.
  drelax(:,:)=0.
  erelax(:,:)=0.
  frelax(:,:)=limit_hw(:,:)
!$OMP END WORKSHARE
  reste      =0.

!$OMP DO
    do J=2,NYY-1
       do I=2,NXX-1

          if (klimit(i,j).eq.0) then
! calcul du vecteur
          FRELAX(I,J)= VIEUXHWATER(I,J)+(BMELT(I,J)-INFILTR)*DT
          frelax(i,j)=frelax(i,j)+(kmx(i,j)*pgmx(i,j)-kmx(i+1,j)*pgmx(i+1,j))*dtsrgdx
          frelax(i,j)=frelax(i,j)+(kmy(i,j)*pgmy(i,j)-kmy(i,j+1)*pgmy(i,j+1))*dtsrgdx
! calcul des diagonales      
          arelax(i,j)=-kmx(i,j)*dtwdx2        ! arelax : diagonale i-1,j 

          brelax(i,j)=-kmx(i+1,j)*dtwdx2      ! brelax : diagonale i+1,j 

          drelax(i,j)=-kmy(i,j)*dtwdx2        ! drelax : diagonale i,j-1 
          
          erelax(i,j)=-kmy(i,j+1)*dtwdx2      ! drelax : diagonale i,j+1 

          crelax(i,j)=1.+((kmx(i,j)+kmx(i+1,j))+(kmy(i,j)+kmy(i,j+1)))*dtwdx2 
                                              !crelax : diagonale i,j 
          else if (klimit(i,j).eq.1) then
             hwater(i,j)=limit_hw(i,j)
!             write(6,*) i,j,hwater(i,j),crelax(i,j),frelax(i,j),arelax(i,j)
          endif
       end do
    end do
!$OMP END DO
!$OMP END PARALLEL

    ! Boucle de relaxation :
    ! ----------------------
    STOPP = .false.

    ntour=0
    stopp=.false.
    Do  WHILE(.NOT.STOPP)
       ntour=ntour+1
!       write(6,*) 'boucle de relaxation numero',ntour
       !$OMP PARALLEL
       !$OMP DO PRIVATE(reste)
       do j=2,NYY-1
          do i=2,NXX-1

             RESTE = (((ARELAX(I,J)*HWATER(I-1,J) + BRELAX(I,J)*HWATER(I+1,J)) &
                  + (DRELAX(I,J)*HWATER(I,J-1) + ERELAX(I,J)*HWATER(I,J+1))) &
                  + CRELAX(I,J)*HWATER(I,J))- FRELAX(I,J)

             DELTAH(I,J) = RESTE/CRELAX(I,J)             
          end do
       end do
       !$OMP END DO

       !$OMP WORKSHARE
       deltah(:,:)=min(deltah(:,:),10.)
       deltah(:,:)=max(deltah(:,:),-10.)
       !$OMP END WORKSHARE
       ! il faut faire le calcul suivant dans une autre boucle car RESTE est fonction
       ! de hwater sur les points voisins.
       !$OMP DO
       do j=2,NYY-1
          do i=2,NXX-1
             HWATER(I,J) = HWATER(I,J) - DELTAH(I,J)
          end do
       end do
       !$OMP END DO

       ! critere d'arret:
       Delh=0
       Vh=0
       
       !$OMP DO REDUCTION(+:Delh)
       DO j=2,NYY-1
          DO i=2,NXX-1
             
             !   write(166,*) I,J,delh,deltah(i,j)
             Delh=Delh+deltah(i,j)**2
             !            Vh=Vh+h(i,j)**2.
          END DO
       END DO
       !$OMP END DO
       !$OMP END PARALLEL
!       write(6,*) delh,maxval(deltah),minval(deltah)
       !      testh=SQRT(Delh/Vh)
       if (delh.gt.0.) then
          testh=SQRT(Delh)/((NXX-2)*(NYY-2))
       else
          testh=0.
       endif
       STOPP = (testh.lt.1.E-3).or.(ntour.gt.1000)
!       write(6,*) ntour,testh


362   continue  
end do
  end subroutine relaxation_waterdif
end module relaxation_waterdif_mod