  MODULE pressure

!!======================================================================
!!                   ***  MODULE  pressure  ***
!! Get the pressure and the pressure gradient to initialize the test case
!!======================================================================





























































   IMPLICIT NONE
   PRIVATE

!! * Accessibility
   PUBLIC pres_gradt      ! function called by cas*.F90


 CONTAINS

!!======================================================================


  SUBROUTINE pres_gradt(boz,xec,pr,ipx,ipy)


!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE pressure  ***
!&E
!&E ** Purpose : estimate the pressure in each cell
!&E
!&E ** Description : the estimate is done according to the choice of the
!&E                  pressure gradient
!&E                  Usefull for the initialization
!&E
!&E ** Called by : cas_acc
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables : pr : pressure at u level
!&E
!&E ** Reference :
!&E                XXXXXXXXXi et al, JPO 2007
!&E
!&E ** History :
!&E       !  2013     (X. Couvelard) Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : rsh,rlg,kmax,ljminm1,ljmaxp2,liminm1,limaxp2,valmanq, &
                          limin,limax,liminp2,limaxp1,limaxm1,liminm2,ljminm2, &
                          ljminp2,ljmaxm1,ljminp1,ljmin,ljmax,ljmaxp1,epsilon,grav
   USE comvars2d, ONLY : h0,ig,id,jb,jh,dy_v,dx_u
   USE comvars3d, ONLY : dsigu,dsigw,sig
   USE comvarp3d, ONLY : uz,vz
   USE comcolig2d,ONLY : igdu,igdv
   USE comcolig3d,ONLY : yz_exp


   USE toolgeom,  ONLY : f_zw,f_zu,f_dzu,f_hzu,f_hzw
   USE toolmpi,ONLY : ex_i_rsh,ex_j_rsh


   REAL(KIND=rsh), DIMENSION(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in) :: boz
   REAL(KIND=rsh), DIMENSION(1:kmax,limin:limaxp1,ljmin:ljmaxp1), INTENT(inout) :: pr
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in) :: xec
   REAL(KIND=rsh), DIMENSION(1:kmax,limin:limax,ljmin:ljmax), INTENT(inout), OPTIONAL :: ipx
   REAL(KIND=rsh), DIMENSION(1:kmax,limin:limax,ljmin:ljmax), INTENT(inout), OPTIONAL :: ipy

!! * Arguments

!! * Local declarations
   INTEGER           :: i,j,k
   REAL(KIND=rsh),DIMENSION(0:kmax) :: dboz,dzr


   REAL(KIND=rsh)    :: bms,bm,pri,tbc,dxb,dhxx,dhtx,htxee,htxeo, & !dyn3dx_ipg_pj.h
                        dyb,dhyy,dhty,htxen,htxes !dyn3dy_ipg_pj.h


!#else
   REAL(KIND=rsh)    :: zukp1,zuk,zwk
!#endif
   REAL(KIND=rsh)    :: verif,dzu


!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter PRESSURE routine'



   DO j=ljmin,ljmaxp1
     DO i=limin,limaxp1

    htxeo=h0(i,j)+xec(i,j)
    k=kmax
    bms=htxeo*boz(kmax,i,j)
    pr(k,i,j)=-bms*sig(kmax)
    DO k=kmax-1,1,-1
      bm=htxeo*boz(k,i,j)
      pr(k,i,j)=pr(k+1,i,j)+(bm+bms)*0.5_rlg*dsigw(k)
      bms=bm
    END DO

     END DO
   END DO



   IF (PRESENT(ipx)) THEN
   ipx(:,:,:)=0.0_rsh
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(ig(j),limin),MIN0(limax,id(j)-1)

       IF(igdu(i,j) == 1) THEN


!*********    including dyn3dx_ipg_pj.h      ****************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dx_ipg_pj.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the internal
!&E              pressure gradient for the zonal momentum forcing
!&E              (3D equation)
!&E
!&E ** Description : explicit
!&E                  vertical integration of the pressure jacobian
!&E
!&E ** Called by : included into dyn3dxtlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xec,h0,boz,dx_u
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  2007-09 (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E
!&E---------------------------------------------------------------------

    dxb=dx_u(i,j)

!    dhxx=(h0(i+1,j)-h0(i,j))/dxb
    dhxx=(xec(i+1,j)-xec(i,j))/dxb
    htxee=h0(i+1,j)+xec(i+1,j)
    htxeo=h0(i,j)+xec(i,j)
    dhtx=(htxee-htxeo)/dxb
    pri=0.0_rlg

    k=kmax
! internal pressure gradient : dP/dx part

    bms=(htxee*boz(kmax,i+1,j)-htxeo*boz(kmax,i,j))/dxb
    pri=pri-bms*sig(kmax)
    tbc=pri+(boz(k,i,j)+boz(k,i+1,j))*0.5_rlg*(sig(k)*dhtx+dhxx)


    yz_exp(k,i,j) = yz_exp(k,i,j)+tbc

    DO k=kmax-1,1,-1

      bm=(htxee*boz(k,i+1,j)-htxeo*boz(k,i,j))/dxb
      pri=pri+(bm+bms)*0.5_rlg*dsigw(k)
      tbc=pri+(boz(k,i,j)+boz(k,i+1,j))*0.5_rlg*(sig(k)*dhtx+dhxx)


      yz_exp(k,i,j) = yz_exp(k,i,j)+tbc
      bms=bm
    END DO

!********* end  including dyn3dx_ipg_pj.h      ****************************


       END IF

     END DO
   END DO
   END IF

   IF (PRESENT(ipy)) THEN
   ipy(:,:,:)=0.0_rsh
   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF(igdv(i,j) == 1) THEN


!*********    including dyn3dy_ipg_pj.h      *******************************

!&E---------------------------------------------------------------------
!&E                 ***  dyn3dy_ipg_pj.h  ***
!&E
!&E ** Purpose : estimate the trend term relative to the internal
!&E              pressure gradient for the meridional momentum forcing
!&E              (3D equation)
!&E
!&E ** Description : explicit
!&E                  vertical integration of the pressure jacobian
!&E
!&E ** Called by : included into dyn3dytlscoef.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : xec,h0,boz,dx_v
!&E
!&E ** Modified variables : yz_exp
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  1987-02  (P. Lazure)  Original code
!&E       !  2007-09 (V. Garnier) homogeneneity of dx dy...
!&E       !  2008-07  (L. Debreu, M. Honnorat, T. Duhaut) time scheme
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E
!&E---------------------------------------------------------------------

    dyb=dy_v(i,j)

!    dhyy=(h0(i,j+1)-h0(i,j))/dyb
    dhyy=(xec(i,j+1)-xec(i,j))/dyb
    htxen=h0(i,j+1)+xec(i,j+1)
    htxes=h0(i,j)+xec(i,j)
    dhty=(htxen-htxes)/dyb
    pri=0.0_rlg

    k=kmax
! internal pressure gradient : dP/dy part

    bms=(htxen*boz(k,i,j+1)-htxes*boz(k,i,j))/dyb
    pri=pri-bms*sig(kmax)
    tbc=pri+(boz(k,i,j)+boz(k,i,j+1))*0.5_rlg*(sig(k)*dhty+dhyy)


    yz_exp(k,i,j)=yz_exp(k,i,j)+tbc
!
    DO k=kmax-1,1,-1

      bm=(htxen*boz(k,i,j+1)-htxes*boz(k,i,j))/dyb
      pri=pri+(bm+bms)*0.5_rlg*dsigw(k)
      tbc=pri+(boz(k,i,j)+boz(k,i,j+1))*0.5_rlg*(sig(k)*dhty+dhyy)

      bms=bm

      yz_exp(k,i,j)=yz_exp(k,i,j)+tbc
    END DO

!********* end  including dyn3dy_ipg_pj.h      ***********************


       END  IF                             ! fin test igdv
     END DO
   END DO

   END IF

   !PRINT*, 'exit PRESSURE routine'

  END SUBROUTINE pres_gradt

!!======================================================================

  END MODULE pressure
