!> \file litho-0.4.f90      
!! Donne la repartition des enfoncements
!! en fonction de la rigidite de la lithosphere.
!<

!> SUBROUTINE: litho()
!!Petit routine qui donne la repartition des enfoncements
!!en fonction de la rigidite de la lithosphere.
!! \author ....
!! \date  10 septembre 2009     
!! @note Cette version utilise la syntaxe f90 pour calculer W1
!!     - avantage : plus facile a lire
!!         pas de pb d'execution avac ifor 10
!!     - inconvenient : pourrait perdre la symetrie.
!! @note  - En entree 
!! @note      WE(-LBLOC:LBLOC,-LBLOC:LBLOC)  : deflection due a une charge unitaire 
!! defini dans tab-litho
!! @note      LBLOC : relie a la distance : distance en noeud autour de laquelle 
!! la flexure de la lithosphere est calculee
!! @note      CHARGE(1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC) : poids par unite de surface
!! (unite ?) au temps time, calcule avant  'appel a litho 
!! dans taubed ou initial2 
!! @note  - En sortie
!! @note      W1(nx,ny)   : enfoncement courant
!! @note  - tableaux de travail.
!! @note      WLOC(-LBLOC:LBLOC,-LBLOC:LBLOC) : enfoncement au point i,j du a la charge 
!! au point IP,JP 
!! @note      LPX et LPY sont les coordonnees entre (-LBLOC:LBLOC,-LBLOC:LBLOC) autour 
!! du point I,J pour calculer tous les WLOC
!! @note      IP et JP sont les coordonnees de CHARGE dans un repere :
!! (1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC)
!! @note Used modules:
!! @note    - use module3D_phy
!! @note    - use ISO_DECLAR
!! @note    - use param_phy_mod
!<

 
subroutine litho
  !$ USE OMP_LIB
  USE module3D_phy
  USE param_phy_mod
  USE ISO_DECLAR ! module de declaration des variables specifiques a l'isostasie

  implicit none

  INTEGER :: II,SOM1,SOM2
  REAL, dimension(:,:), allocatable :: WLOC


  !----- allocation de WLOC  et de croix -----------

  if (.not.allocated(WLOC)) THEN
     allocate(WLOC(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err)
     if (err/=0) then
	print *,"Erreur a l'allocation du tableau WLOC",err
	stop 4
     end if
  end if

  !----- fin de l'allocation --------------

  !     calcul de la deflexion par sommation des voisins appartenant
  !     au bloc de taille 2*LBLOC
  som1=0.
  som2=0.

  !      On somme aussi les contributions des points exterieurs au domaine
  !      lorsque la charge est due a l'ocean. On suppose alors  que
  !      ces points ont la meme charge que les limites
  !$OMP PARALLEL
  !$OMP DO PRIVATE(Wloc)
  boucleij : do J=1,NY
     do I=1,NX


	W1(I,J)=0. 
 !   ii=0

 ! Wloc : enfoncement centre sur ij du a chaque charges autour
	Wloc(:,:) = We(:,:) * charge(i-lbloc:i+lbloc,j-lbloc:j+lbloc)


	W1(i,j) = sum (Wloc(:,:))   ! effet en ij de toutes les charges autour (distance < lbloc)

 !    ------------------------ FIN DE L'INTEGRATION SUR LE PAVE LBLOC
 !   som1=som1+W1(I,J)
 !   som2=som2-charge(I,J)/ROMG

     end do
  end do  boucleij
  !$OMP END DO
  !$OMP END PARALLEL

  !     write(6,*)'enfoncement total avec rigidite de la lithosphere:',
  !    &  som1
  !     write(6,*)'enfoncement total avec isostasie locale :',
  !    &  som2


end subroutine litho