!> \file isostasie_mod-0.3.f90
!!Package isostasie 
!<

!> \namespace isostasie_mod
!! Module pour le calcul de l'isostasie
!! \author ...
!! \date ...
!! @note Used modules
!! @note    - use iso_declar
!! @note    - use module3D_phy
!! @note    - use param_phy_mod
!<



module isostasie_mod

  use iso_declar
  
contains

!> SUBROUTINE: init_iso
!! Routine qui permet d'initialiser les variables 
!>
  subroutine init_iso ! routine qui permet d'initialiser les variables 

    use module3D_phy
    use param_phy_mod

    implicit none

    !     nbed=0 pas d'isostasie
    !     nbed=1 temps de reaction


    if (GEOPLACE.eq.'eismint') then
       NBED=0
       NLITH=0
    else if (GEOPLACE(1:6).eq.'marine') then
       NBED=1
       NLITH=1
    else if ((geoplace.eq.'anteis1').or.(geoplace.eq.'ant20km')) then

       if (icouple.eq.2) then
          NBED=0
          NLITH=0
          if (marine) then
             NBED=1
             NLITH=1
          endif
       else
          NBED=1
          NLITH=1
       endif
    else
       NBED=1
       !       switch sur la lithosphere
       NLITH=1
    endif

    !      temps de reaction en annees
    if (NBED.eq.1) tausoc=3000

    if (NLITH.eq.1) then
       !        DL lithosphere flexural rigidity (N.m)
       DL=9.87E24
       !        radius of relative stiffness  (metre)
       RL=131910.
       !        LBLOC, 400 km de part et d'autre
       LBLOC=int((400000.-0.1)/DX)+1
       !         LBLOC=int((480000.-0.1)/DX)+1
       !         LBLOC=int(400000./DX)+1


       !******** allocation des tableaux en fonction de la valeur de LBLOC *****

       if (.not.allocated(WE)) then
          allocate(WE(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err)
          if (err/=0) then
             print *,"Erreur a l'allocation du tableau WE",err
             stop 4
          end if
       end if
       !    PRINT*,'shape(we)',SHAPE(we),'lbloc',lbloc

       if (.not.allocated(CHARGE)) then
          allocate(CHARGE(1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC),stat=err)
          if (err/=0) then
             print *,"Erreur a l'allocation du tableau CHARGE",err
             stop 4
          end if
       end if
       !********* FIN DE L'ALLOCATION DES TABLEAUX ***********
       call tab_litho
    end if

    !******** initialisation de CHARGE ***********

    !       pour calcul du socle initial on calcule la charge
    !       avec l'etat initial suppose en equilibre S0, H0, Bsoc0

    do i=1,nx
       do j=1,ny
          if (ro*h0(i,j).ge.-row*(BSOC0(i,j)-sealevel)) then

             !        glace ou terre
             charge(i,j)=rog*h0(i,j)

          else
             !        ocean
             charge(i,j)=-rowg*(Bsoc0(i,j)-sealevel)
          endif
       end do
    end do

    do j=1,ny  ! parties de charge a l'exterieure de la grille
       charge(1-lbloc:0,j)=charge(1,j)
       charge(nx+1:nx+lbloc,j)=charge(nx,j)
    end do
    do i=1,nx
       charge(i,1-lbloc:0)=charge(i,1)
       charge(i,ny+1:ny+lbloc)=charge(i,ny)
    end do

    do j=1,ny  ! parties de charge a l'exterieure de la grille
       charge(1-lbloc:0,j)=charge(1,j)
       charge(nx+1:nx+lbloc,j)=charge(nx,j)
    end do
    do i=1,nx
       charge(i,1-lbloc:0)=charge(i,1)
       charge(i,ny+1:ny+lbloc)=charge(i,ny)
    end do
    charge(1-lbloc:0,1-lbloc:0)=charge(1,1) !valeurs aux quatres coins
    charge(1-lbloc:0,ny+1:ny+lbloc)=charge(1,ny) ! exterieurs au domaine
    charge(nx+1:nx+lbloc,1-lbloc:0)=charge(nx,1)
    charge(nx+1:nx+lbloc,ny+1:ny+lbloc)=charge(nx,ny)

    !     deflection de la lithosphere
    if (nlith.eq.1) then

       call litho

       do j=1,ny
          do i=1,nx
             w0(i,j)=w1(i,j)
          end do
       end do
    else
       do j=1,ny
          do i=1,nx
             w0(i,j)=charge(i,j)/romg
          end do
       end do
    end if

  end subroutine init_iso



!> SUBROUTINE: bedrock
!! Routine qui calcule l'isostasie 
!>

subroutine bedrock
  call taubed  ! routine qui calcule l'isostasie
end subroutine bedrock

end module isostasie_mod