!> \file taubed-0.3.f90
!! Calcul la charge en chaque point de la grille puis appel la routine litho pour calculer la contribution 
!! de chaque point a la deflexion de la lithosphere    
!<

!> SUBROUTINE: taubed()
!! \author ...
!! \date  16 Novembre 1999
!! @note  Routine qui calcul la charge en chaque point de la grille
!! puis appel la routine litho pour calculer la contribution 
!! de chaque point a la deflexion de la lithosphere    
!! @note En entree
!!        - H, 
!!        - BSOC, 
!!        - SEALEVEL, 
!!        - RO, 
!!        - ROW
!! @note En sortie
!!        - CHARGE(1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC) : poids par unite de surface
!!               (unite ?)   Elle est calculee initialement dans initial2
!!               Poids de la colonne d'eau ou de la colonne de glace.
!!               a l'exterieur du domaine : 1-LBLOC:1 et NX+1:NX+LBLOC
!!               on donne les valeurs des bords de la grille
!!               CHARGE est utilise par litho uniquement
!!        - W1(NX,NY) est l'enfoncement courant, c'est le resultat 
!!               de la routine litho
!!               W1 peut etre calcule de plusieurs facons selon le modele 
!!               d'isostasie utilise
!! @note Used module:
!! @note   - use module3D_phy
!! @note   - use param_phy_mod
!! @note   - use ISO_DECLAR
!!
!<

!!     ****************************************************************
!!     *           BEDROCK ADJUSTMENT avec temps de reaction          *
!!     *  changement de nom B -> Bsoc                                 * 
!!     ****************************************************************



subroutine taubed()

  !$USE OMP_LIB 
  USE module3D_phy
  USE param_phy_mod
  USE ISO_DECLAR ! module de declaration des variables de l'isostasie

  implicit none

  !    ********* calcul de W1 l'enfoncement d'equilibre au temps t
  ! NLITH est defini dans isostasie et permet le choix du modele d'isostasie

  
  if (NLITH.eq.1) then
     !       avec rigidite de la lithosphere
     !$OMP PARALLEL
     !$OMP DO
     do J=1,NY 
        do I=1,NX
           if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then
              !           glace ou terre 
              CHARGE(I,J)=ROG*H(I,J)
           else
              !           ocean
              CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL_2D(i,j))
           endif
        end do
     end do
     !$OMP END DO

     ! il faut remplir CHARGE dans les parties a l'exterieur de la grille :
     ! a l'exterieur de la grille CHARGE est egale a la valeur sur le bord
     !$OMP DO
     do J=1,NY
        CHARGE(1-LBLOC:0,J)=CHARGE(1,J)      ! bord W
        CHARGE(NX+1:NX+LBLOC,J)=CHARGE(NX,J) ! bord E
     end do
     !$OMP END DO
     !$OMP DO
     do I=1,NX
        CHARGE(I,1-LBLOC:0)=CHARGE(I,1)      ! bord S
        CHARGE(I,NY+1:NY+LBLOC)=CHARGE(I,NY) ! bord N
     end do
     !$OMP END DO

     ! valeur dans les quatres coins exterieurs au domaine       
	  !$OMP WORKSHARE
     CHARGE(1-LBLOC:0,1-LBLOC:0)=CHARGE(1,1)           ! coin SW
     CHARGE(1-LBLOC:0,NY+1:NY+LBLOC)=CHARGE(1,NY)      ! coin NW
     CHARGE(NX+1:NX+LBLOC,1-LBLOC:0)=CHARGE(NX,1)      ! coin SE
     CHARGE(NX+1:NX+LBLOC,NY+1:NY+LBLOC)=CHARGE(NX,NY) ! coin NE
     !$OMP END WORKSHARE
     !$OMP END PARALLEL 
     call litho

  else
     !     enfoncement local
     !$OMP PARALLEL
     !$OMP DO
     do J=1,NY
        do I=1,NX
           if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then
              !           glace ou terre 
              W1(I,J)=RO/ROM*H(I,J)
           else
              !           ocean
              W1(I,J)=-ROW/ROM*(BSOC(I,J)-SEALEVEL_2D(i,j))
           endif
        end do
     end do
     !$OMP END DO
     !$OMP END PARALLEL
  endif

  !     decroissance exponentielle de l'enfoncement
  !$OMP PARALLEL
  !$OMP DO
  do J=1,NY
     do I=1,NX
        BDOT(I,J)=((Bsoc0(I,J)-BSOC(I,J))-(W1(I,J)-W0(I,J)))/TAUSOC
        BSOC(I,J)=BSOC(I,J)+BDOT(I,J)*DT_ISO
     end do
  end do
  !$OMP END DO
  !$OMP END PARALLEL

end subroutine taubed