source: trunk/SOURCES/litho-0.4.f90 @ 4

!> \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 module3D_phy
      USE param_phy_mod
      USE ISO_DECLAR ! module de declaration des variables specifiques a l'isostasie

      implicit none

      INTEGER :: IP,JP,LPX,LPY,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

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
       
!     write(6,*)'enfoncement total avec rigidite de la lithosphere:',
!    &  som1
!     write(6,*)'enfoncement total avec isostasie locale :',
!    &  som2


    end subroutine litho