!> \file litho-0.3.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 Novembre 1999         
!! @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
!<
subroutine litho

  USE module3D_phy
  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
  REAL, dimension(:), allocatable :: croix

  !----- 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

  if (.not.allocated(croix)) THEN
     allocate(croix(0:LBLOC),stat=err)
     if (err/=0) then
        print *,"Erreur a l'allocation du tableau croix",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

        !          IP,JP compteurs de la sommation dans le repere 
        !                  (1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC)

        !          LPX,LPY compteurs dans le repere -LBLOC:LBLOC

        !      do IP=I-LBLOC,I+LBLOC
        !         do JP=J-LBLOC,J+LBLOC
        !            LPX = IP-I
        !            LPY = JP-J
        !            WLOC(LPX,LPY)=We(LPX,LPY)*CHARGE(IP,JP)
        !            ii=ii+1
        !         end do
        !      end do
        ! ecriture tableau f90 de la boucle ci-dessus :

        WLOC(:,:) = WE(:,:) * CHARGE(I-LBLOC:I+LBLOC,J-LBLOC:J+LBLOC)



        !         fin de l'attribution des wloc


        !           sommation de tous les effets (tentative pour
        !           eviter les erreurs d'arrondi)
        W1(I,J)=WLOC(0,0)

        !  dans croix on calcul la somme des effets WLOC puis on met le resultat dans W1   
        do LPX=1,LBLOC
           LPY=0
           !              print*,i,j,LPX,WLOC(LPX,LPY),WLOC(-LPX,LPY),WLOC(LPY,LPX),WLOC(LPY,-LPX)
           !               print*,(WLOC(LPX,LPY)+WLOC(-LPX,LPY))
           !                print*,(WLOC(LPY,LPX)+WLOC(LPY,-LPX))
           !                 CROIX(LPY)=(WLOC(LPX,LPY)+WLOC(-LPX,LPY))

           !                 CROIX(LPY)=(WLOC(LPY,LPX)+WLOC(LPY,-LPX)) + CROIX(LPY) 


           !                 write(6,*) i,j,lpy,lpx,lbloc,WLOC(LPY,LPX)
           !                 write(6,*) WLOC(LPY,-LPX),WLOC(LPX,LPY),WLOC(-LPX,LPY)
           !                 write(6,*) (WLOC(LPY,LPX)+WLOC(LPY,-LPX))+(WLOC(LPX,LPY)+WLOC(-LPX,LPY))
           CROIX(LPY)=(WLOC(LPY,LPX)+WLOC(LPY,-LPX))+(WLOC(LPX,LPY)+WLOC(-LPX,LPY))


           !                 print*,'croix ',croix(lpy)
           LPY=LPX
           CROIX(LPY)= ((WLOC(LPX,LPY)+WLOC(LPX,-LPY)) &
                +(WLOC(-LPX,LPY)+WLOC(-LPX,-LPY))) 
           do LPY=1,LPX-1
              CROIX(LPY)= (((WLOC(LPX,LPY)+WLOC(LPX,-LPY)) &
                   +(WLOC(-LPX,LPY)+WLOC(-LPX,-LPY))) &
                   +((WLOC(LPY,LPX)+WLOC(LPY,-LPX)) &
                   +(WLOC(-LPY,LPX)+WLOC(-LPY,-LPX))))
           end do
           do LPY=0,LPX 
              W1(I,J)=W1(I,J)+CROIX(LPY) ! sommation dans W1
           end do

        end do


        !      if (ii.ne.(2*LBLOC+1)*(2*LBLOC+1)) then
        !         write(6,*)'attention mauvais compteur pour i,j',I,J,II
        !      endif
        !    ------------------------ 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

  deallocate(WLOC,stat=err)
  if (err/=0) then
     print *,"Erreur a l'allocation du tableau WLOC",err
     stop 4
  end if

  deallocate(croix,stat=err)
  if (err/=0) then
     print *,"Erreur a l'allocation du tableau croix",err
     stop 4
  end if
  return
end subroutine litho