!> \file tab-litho-0.3.f90
!! Repartition des enfoncement en fonction de la rigidite de la lithosphere.
!<

!> SUBROUTINE:  tab_litho 
!! \author ...
!! \date ...
!!  @note Subroutine qui donne la repartition des enfoncements
!!        en fonction de la rigidite de la lithosphere.
!!  @note  variables en entree
!!       - ROM masse volumique du manteau
!!       - RL  rayon de rigidite relative
!!       - DL  rigidite flexural
!!  
!!   @note variables en sortie
!!       - WE  deflection due a une charge unitaire      
!!
!! @note Used modules:
!! @note    - use module3D_phy
!! @note    - use param_phy_mod
!! @note    - use ISO_DECLAR
!<

subroutine tab_litho 

  USE module3D_phy
  USE param_phy_mod
  USE ISO_DECLAR ! module de declaration pour l'isostasie


  implicit none


  INTEGER, PARAMETER :: nk=1000
  REAL, dimension(0:nk) :: kei
  REAL :: stepk,AL,XL,DIST
  REAL :: som

  ! pour la lithosphere
  STEPK=100.
  AL=-RL*RL/(2.*PI*DL)*DX*DY      

  !      fonction de kelvin
  !      lecture de la table kei qui est tous les 0.01 entre 0 et 10
  !      STEPK=100=1/ecart 

  open(num_kelvin,file=trim(dirnameinp)//'../kelvin.res')
  read(num_kelvin,*)
  do K=0,NK
     read(num_kelvin,*) XL,kei(K)
  end do
  close(num_kelvin)


  do I=-LBLOC,LBLOC    
     do J=-LBLOC,LBLOC
        DIST=dx*sqrt(1.*(I*I+J*J))                                
        XL=DIST/RL*STEPK                                          
        K=int(XL)
        if ((K.gt.834).or.(DIST.gt.dx*LBLOC)) then                 
           We(I,J)=0.                                              
        else 
           We(I,J)=kei(k)+(kei(k+1)-kei(k))*(XL-K)
           if (K.eq.834) We(I,J)=max(We(I,J),0.)                   
        endif
        We(I,J)=We(I,J)*Al   
     end do
  end do



  !open(55,file='iso-cat')
  !write(55,*) 'WE avec LBLOC sur 480 km'
  !do i=-lbloc,lbloc
  !   write(55,'25(i4,1x)') (int(WE(i,j)*1.e08),j=-lbloc,lbloc)
  !end do
  !close(55)

  ! sommation des elements de WE dans som pour normaliser
  ! la somme des enfoncements a 1
  som = SUM(WE(:,:))


  ! normalisation
  WE(:,:) = WE(:,:)/(som*ROM*g)



  ! ATTENTION DESSOUS VALABLE SEULEMENT EN 40 km

  ! la partie ci-dessous doit etre mise dans une routine specifique Philippe
  !       pour faire exactement la meme isostasie que Philippe
  !if ((icouple.eq.2).or.(icouple.eq.4)) then

  !   LBLOC=10

  !   do I=-LBLOC,LBLOC
  !      do J=-LBLOC,LBLOC
  !         We(i,j)=0.   ! mise a 0 partout
  !      end do
  !   end do
  !   open(num_phil,file='../INPUT/iso-philippe')

  !   read(num_phil,*)
  !   read(num_phil,*)
  !   read(num_phil,*)
  !   read(num_phil,*)
  !   do J=-10,10
  !      read(num_phil,*)(We(i,j),i=-10,10)
  !   end do
  !   close(num_phil)
  !   do i=-LBLOC,LBLOC
  !      do j=-LBLOC,LBLOC
  !         we(i,j)=we(i,j)/11030/g/rom
  !      end do
  !   end do
  !endif

  return
end subroutine tab_litho