source: branches/GRISLIv3/SOURCES/neffect-0.4.f90 @ 465

!> \file neffect-0.4.f90
!! Calcule la pression effective
!<

!> SUBROUTINE: Neffect()
!! \author Catherine 
!! \date 18 janvier 2005
!! @note Cette routine permet de calculer la pression effective
!! @note Used modules:
!! @note    - use module3D_phy
!! @note    - use module_choix
!<
subroutine Neffect

  !     ===============================================================
  !     *    Pression effective                 Cat   18 janvier 2005
  !     ===============================================================

  use param_phy_mod, only: ro,row,g
  use runparam, only: itracebug
  use module3D_phy, only: flotmx,flotmy,coefmxbmelt,coefmybmelt,hwater,neffmx,neffmy,hmx,hmy
  use geography, only: nx,ny
!  use module_choix, only: eaubasale <- afq, not clean, I prefer to use  explicitly eau_basale
  use eau_basale, only: hwatermax,eaubasale

  implicit none

  integer :: i,j
  !cdc      real :: Nefmin=1.e5                  ! Pression effective minimum (~ 10 m de glace)
  !afq  - no longer used - real :: Nefmin=0.  !cdc 22/01/2019 on autorise une pression effective nulle
  real :: nefminfrac = 0. ! put 0.04 to have something similar to PISM

  if (itracebug.eq.1)  call tracebug(' Entree dans routine neffect')




  ! le calcul de pente est fait dans diffusiv


  ! calcul de la hauteur d'eau dans eaubasale


  call eaubasale

  if (itracebug.eq.1)  call tracebug(' Apres eau basale')
  ! calcul de la presion effective Neffmy

  do j=2,ny
     do i=2,nx-1
        pose_y :    if (.not.flotmy(i,j)) then

           !     coefmybmelt est la hauteur d'eau au point considere
           !     moyenne pondérée des points voisins.
           !     coefmybmelt est utilise dans la vitesse de glissement 

           !    ancienne version : moyenne large
           !             COEFMYBMELT(I,J)=(HWATER(I,J)+HWATER(I,J-1))/4.+
           !     &               ((HWATER(I+1,J)+HWATER(I+1,J-1))+
           !     &                (HWATER(I-1,J)+HWATER(I-1,J-1)))/8.

           !    nouvelle version, moyenne sur les deux plus proches voisins pour mieux 
           !    localiser les fleuves de glace

           coefmybmelt(i,j)=(hwater(i,j)+HWATER(i,j-1))*0.5
           Neffmy(i,j)=ro*g*hmy(i,j)-row*g*coefmybmelt(i,j)
!!!!!        Neffmy(i,j)=RO*G*HMY(I,J)-pwater(i,j)

           !             la pression effective mini est Nefmin (en Pa) pour la glace posée 
! afq --           Neffmy(i,j)=max(Neffmy(i,j),Nefmin)
           Neffmy(i,j)=max(Neffmy(i,j),nefminfrac*RO*G*HMY(I,J)) !PISM: nefminfrac=4%

           !     coefmybmelt est utilise dans la vitesse de glissement 
           !     et pour cela est limité par hwatermax
           coefmybmelt(i,j)=min(hwatermax,coefmybmelt(i,j))



           ! on calcule tobmy dans diagnoshelf apres etre passe dans remplimat
           ! pour utiliser les nouvelles vitesses.

        else
           !           points completement flottants
           Neffmy(i,j)=0.

        endif  pose_y
     end do
  end do

  ! calcul de la presion effective Neffmx

  do j=2,Ny-1
     do i=2,nx
        pose_x :       if (.not.flotmx(i,j)) then

           !     coefmybmelt est la hauteur d'eau au point considere
           !     moyenne pondérée des points voisins.
           !     coefmybmelt est utilise dans la vitesse de glissement 

           !    ancienne version : moyenne large
           !             COEFMYBMELT(I,J)=(HWATER(I,J)+HWATER(I-1,J))/4.+
           !     &               ((HWATER(I,J+1)+HWATER(I-1,J+1))+
           !     &                (HWATER(I,J-1)+HWATER(I-1,J-1)))/8.

           !    nouvelle version, moyenne sur les deux plus proches voisins pour mieux 
           !    localiser les fleuves de glace

           coefmxbmelt(i,j)=(hwater(i,j)+hwater(i-1,j))*0.5
           Neffmx(i,j)=ro*g*hmx(I,J)-row*g*coefmxbmelt(i,j)
!!!!!        Neffmx(i,j)=RO*G*HMX(I,J)-pwater(i,j)

           !             la pression effective mini est 1 Pascal pour la glace posée 
! afq --           Neffmx(i,j)=max(Neffmx(i,j),Nefmin)
           Neffmx(i,j)=max(Neffmx(i,j),nefminfrac*RO*G*HMX(I,J)) !PISM: nefminfrac=4%

           !     coefmybmelt est utilise dans la vitesse de glissement 
           !     et pour cela est limité par hwatermax
           coefmxbmelt(i,j)=min(hwatermax,coefmxbmelt(i,j))




           ! on calcule tobmy dans diagnoshelf apres etre passe dans remplimat
           ! pour utiliser les nouvelles vitesses.

        else
           !           points completement flottants
           Neffmx(i,j)=0.

        endif  pose_x
     end do
  end do

  ! nom_table='pwater'
  ! table_interm=ROW*G*COEFMXBMELT/10e5      
  ! call printtable_r(table_interm,nom_table)         

  if (itracebug.eq.1)  call tracebug(' Sortie routine neffect')
  return
end subroutine Neffect