  MODULE pressure

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

#include "toolcpp.h"

   IMPLICIT NONE
   PRIVATE

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


 CONTAINS

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

#ifdef key_profil_back
  SUBROUTINE pres_gradt(boz,xec,pr,ipx,ipy,boz2)
#else
  SUBROUTINE pres_gradt(boz,xec,pr,ipx,ipy)
#endif

   !&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
#if key_dyn_pg_djcs
   USE comcolig3d, ONLY : zr
#endif
#if defined key_siggen || defined key_gencoord
   USE comsiggen,  ONLY : b_sig,theta_sig,hc,hc_sig,csu_sig
#endif
   USE toolgeom,  ONLY : f_zw,f_zu,f_dzu,f_hzu,f_hzw
   USE_MPI toolmpi,ONLY : ex_i_rsh,ex_j_rsh

#ifdef key_profil_back
   REAL(KIND=rsh), OPTIONAL, DIMENSION(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in) :: boz2
#endif
   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
#ifdef key_dyn_pg_djcs
   REAL(KIND=rsh), PARAMETER :: demi   = 0.5_rsh
   REAL(KIND=rsh), PARAMETER :: uncinq = 0.2_rsh
   REAL(KIND=rsh), PARAMETER :: undouz = 1.0_rsh/12.0_rsh
   REAL(KIND=rsh)            :: cff,tbc,dxzr,dxboz, &
                                dxzrp1,dxbozp1,     &
                                dxzrp2,dxbozp2,     &
                                dyboz,dybozp1,      &
                                dybozp2,dyzr,dyzrp1,&
                                dyzrp2
#elif key_dyn_pg_pint
   REAL(KIND=rsh),DIMENSION(kmax)   :: p_primu_i,p_primu_ip1,sec_gp, &
                                       p_primu_j,p_primu_jp1
   REAL(KIND=rsh),DIMENSION(0:kmax) :: p_prim_interface_i,  &
                                       p_prim_interface_ip1,&
                                       p_prim_interface_j,  &
                                       p_prim_interface_jp1
   REAL(KIND=rsh)                   :: hz_i_kp1,hz_i,hz_ip1_kp1,hz_ip1, &
                                       zk_i,zk_ip1,htxee,htxeo,dxb,     &
                                       hz_j_kp1,hz_j,hz_jp1_kp1,hz_jp1, &
                                       zk_j,zk_jp1,htxes,htxen,dyb
#elif key_dyn_pg_dj
   REAL(KIND=rsh)    :: dxb,jacobian,dxbs,dszx,dsbx,dxzs,z1,z2,z3,z4,dzb,b2surf,b4surf
   REAL(KIND=rsh)    :: dyb,dybs,dszy,dsby,dyzs
   REAL(KIND=rsh)    :: f31,f43,f42,f21
#else
#if defined key_siggen || defined key_gencoord
   REAL(KIND=rlg)    :: bms,bm,pri,tbc,dxb,dhxx,dhcx,dxex,bhze,bhzo, & !dyn3dx_ipg_pj_sig.h
                        dyb,dhyy,dhcy,dxey,bhzn,bhzs       !dyn3dy_ipg_pj_sig.h
#else
   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
#endif
#endif
!#else
   REAL(KIND=rsh)    :: zukp1,zuk,zwk
!#endif
   REAL(KIND=rsh)    :: verif,dzu


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

   PRINT_DBG*, 'enter PRESSURE routine'

#ifdef key_dyn_pg_dj

   DO j=ljmin,ljmaxp1
     DO i=limin,limaxp1
       zuk  = f_zu(h0(i,j),xec(i,j),kmax,i,j)
       pr(kmax,i,j)=boz(kmax,i,j)*(0.0_rsh-zuk)
       DO k=kmax-1,1,-1
         zukp1= f_zu(h0(i,j),xec(i,j),k+1,i,j)
         zuk  = f_zu(h0(i,j),xec(i,j),k  ,i,j)
         zwk  = f_zw(h0(i,j),xec(i,j),k,i,j)
         pr(k,i,j) = pr(k+1,i,j)+boz(k+1,i,j)*(zukp1-zwk)
         pr(k,i,j) = pr(k,i,j)+boz(k,i,j)*(zwk-zuk)
       END DO
     END DO
   END DO

#elif key_dyn_pg_pint

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

     dxb=dx_u(i,j)
     htxeo=h0(i,j)+xec(i,j)

   ! surface - estimate of the second pressure gradient term : boz*dz/dx

     k = kmax

     p_prim_interface_i  (k)=0.0_rlg

#if defined key_siggen || defined key_gencoord
     hz_i   = f_hzu(h0(i,j),xec(i  ,j),k,i  ,j)
#else
     hz_i   = htxeo
#endif
     zk_i   = f_zu(h0(i  ,j),xec(i  ,j),k,i  ,j)

     DO k=kmax-1,1,-1
       hz_i_kp1 = hz_i
#if defined key_siggen || defined key_gencoord
       hz_i   = f_hzu(h0(i,j),xec(i  ,j),k,i  ,j)
#else
       hz_i   = htxeo
#endif
       zk_i   = f_zu(h0(i  ,j),xec(i  ,j),k,i  ,j)

       ! horizontal gradient of pressure (sigma = cst)
       p_prim_interface_i(k)   = p_prim_interface_i(k+1)   + boz(k+1,i,j)*hz_i_kp1  *dsigu(k+1)
       p_primu_i(k+1)  =0.5_rlg*(p_prim_interface_i(k)  +p_prim_interface_i(k+1)  )

!cval 0.5_rlg*(p_prim_interface_jp1(k)+p_prim_interface_jp1(k+1)) alors que u n est pas
!cval necessairement au milieu de la couche ?

     END DO

   ! bottom estimate of the first pressure gradient term : d/dx integral(boz) (k=1 level)

     k = 0
     hz_i_kp1 = hz_i
     p_prim_interface_i(k)   = p_prim_interface_i(k+1)   + boz(k+1,i,j)*hz_i_kp1  *dsigu(k+1)
     p_primu_i(k+1)  =0.5_rlg*(p_prim_interface_i(k)  +p_prim_interface_i(k+1))

    do k=1,kmax
#ifdef key_ktest_iniv3d
     pr(k,i,j)=p_primu_i(k)
#endif
    end do

     END DO
   END DO


#elif key_dyn_pg_djcs

   !DO j = ljminp2, ljmaxm1
    ! DO i = MAX0(limin,ig(j)), MIN0(limax,id(j))
   DO j=ljmin,ljmax
     DO i=limin,limax

       IF (h0(i,j) > -valmanq) THEN
         zr(kmax,i,j) = f_zu(h0(i,j),xec(i,j),kmax,i,j)
         DO k = kmax-1, 1, -1
           zr(k,i,j) = f_zu(h0(i,j),xec(i,j),k,i,j)
           dzr(k) = zr(k+1,i,j)-zr(k,i,j)
           dboz(k) = boz(k+1,i,j)-boz(k,i,j)
         END DO
         dzr(kmax) = dzr(kmax-1)
         dboz(kmax) = dboz(kmax-1)
         dzr(0) = dzr(1)
         dboz(0) = dboz(1)

         DO k = kmax,1,-1
           cff = 2.0_rsh*dboz(k)*dboz(k-1)
           IF (cff > epsilon) THEN
             dboz(k) = cff / (dboz(k)+dboz(k-1))
           ELSE
             dboz(k) = 0.0_rsh
           END  IF
           dzr(k) = 2.0_rsh*dzr(k)*dzr(k-1) / (dzr(k)+dzr(k-1))
         END DO
         pr(kmax,i,j) = (  boz(kmax,i,j) + &
                  0.5_rsh * (xec(i,j)-zr(kmax,i,j)) * (boz(kmax,i,j)-boz(kmax-1,i,j)) / (zr(kmax,i,j)-zr(kmax-1,i,j))  ) &
                  * (xec(i,j)-zr(kmax,i,j))
         DO k = kmax-1, 1, -1
           pr(k,i,j) = pr(k+1,i,j) + &
                       demi*(  (boz(k+1,i,j)+boz(k,i,j))*(zr(k+1,i,j)-zr(k,i,j)) &
                       - uncinq*( (dboz(k+1)-dboz(k))*(zr(k+1,i,j)-zr(k,i,j)   &
                                           -undouz*(dzr(k+1)+dzr(k)))          &
                                 -(dzr(k+1)-dzr(k))*(boz(k+1,i,j)-boz(k,i,j)   &
                                           -undouz*(dboz(k+1)+dboz(k))) )  )
         END DO

       ELSE
         pr(:,i,j)=valmanq
       END IF

     END DO
   END DO

   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,zr(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,zr(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_i_rsh( 0,1,kmax,limin  ,limaxp1,ljmin  ,ljmaxp1,pr(1:kmax,limin  :limaxp1,ljmin  :ljmaxp1))
   CALL_MPI ex_j_rsh( 0,1,kmax,limin  ,limaxp1,ljmin  ,ljmaxp1,pr(1:kmax,limin  :limaxp1,ljmin  :ljmaxp1))

#else

   DO j=ljmin,ljmaxp1
     DO i=limin,limaxp1
#if defined key_siggen || defined key_gencoord
    k=kmax
    bhzo=f_hzw(h0(i  ,j),xec(i  ,j),k,i  ,j)*boz(k,i  ,j)
    pr(k,i,j)=bhzo*(0.0_rlg-sig(k))
    DO k=kmax-1,1,-1
      bhzo=(f_hzw(h0(i  ,j),xec(i  ,j),k,i  ,j))*0.5_rlg*(boz(k,i  ,j)+boz(k+1,i  ,j))
      pr(k,i,j)=pr(k+1,i,j)+bhzo*dsigw(k)
    ENDDO
#else
    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
#endif
     END DO
   END DO

#endif

   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
#ifdef key_dyn_pg_djcs
#          include "../DYN3DX/dyn3dx_ipg_cs.h"
#elif key_dyn_pg_pint
#          include "../DYN3DX/dyn3dx_ipg_pint.h"
#elif key_dyn_pg_dj
#          include "../DYN3DX/dyn3dx_ipg_dj.h"
#else
#if defined key_siggen || defined key_gencoord
#          include "../DYN3DX/dyn3dx_ipg_pj_sig.h"
#else
#          include "../DYN3DX/dyn3dx_ipg_pj.h"
#endif
#endif
       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
#ifdef key_dyn_pg_djcs
#          include "../DYN3DY/dyn3dy_ipg_cs.h"
#elif key_dyn_pg_pint
#          include "../DYN3DY/dyn3dy_ipg_pint.h"
#elif key_dyn_pg_dj
#          include "../DYN3DY/dyn3dy_ipg_dj.h"
#else
#if defined key_siggen || defined key_gencoord
#          include "../DYN3DY/dyn3dy_ipg_pj_sig.h"
#else
#          include "../DYN3DY/dyn3dy_ipg_pj.h"
#endif
#endif
       END  IF                             ! fin test igdv
     END DO
   END DO

   END IF

   PRINT_DBG*, 'exit PRESSURE routine'

  END SUBROUTINE pres_gradt

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

  END MODULE pressure
