!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!MNH_LIC for details. version 1.
!------------------------------------------------------------------------------
!    #######################
     MODULE MODI_IBM_BALANCE
!    ####################### 
!
INTERFACE
!
SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
!
REAL, DIMENSION(:,:,:,:)     ,INTENT(IN)    :: PPHI
REAL, DIMENSION(:,:,:,:)     ,INTENT(IN)    :: PVOL
REAL, DIMENSION(:,:,:)       ,INTENT(INOUT) :: PRUS,PRVS,PRWS
REAL, DIMENSION(:,:,:)       ,INTENT(INOUT) :: PBAL
!
!------------------------------------------------------------------------------
!
END SUBROUTINE IBM_BALANCE
!
END INTERFACE
!
END MODULE MODI_IBM_BALANCE
!
!     ################################################
      SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
!     ################################################
!
!!
!!****  IBM_BALANCE computes the velocity divergence using a volumic approach 
!!                                
!!    PURPOSE
!!    -------
!!****  The purpose of this routine is to compute div(U)=1/V*int_S(u.n)dS
!!      S is the modified surface and is estimated before MNH
!!      U is approximated using adjacents points  
!
!!    METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!      SUBROUTINE ?
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!       MODD_?   
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!      Franck Auguste (CERFACS-AE)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original         01/01/2015
!!
!------------------------------------------------------------------------------
!       
!**** 0. DECLARATIONS
!     ---------------
! module
USE MODE_ll
USE MODE_IO_ll
!
! declaration
USE MODD_IBM_PARAM_n        
USE MODD_GRID_n,          ONLY: XXHAT,XYHAT,XZHAT,XZZ 
USE MODD_PARAMETERS,      ONLY: JPVEXT,JPHEXT
USE MODD_LBC_n
USE MODD_REF_n
!
! interface
USE MODI_IBM_WRITE
USE MODI_SHUMAN
!
IMPLICIT NONE
!
!------------------------------------------------------------------------------
!
!       0.1  declarations of arguments
!
REAL, DIMENSION(:,:,:,:)     ,INTENT(IN)    :: PPHI
REAL, DIMENSION(:,:,:,:)     ,INTENT(IN)    :: PVOL
REAL, DIMENSION(:,:,:)       ,INTENT(INOUT) :: PRUS,PRVS,PRWS
REAL, DIMENSION(:,:,:)       ,INTENT(INOUT) :: PBAL
!
!------------------------------------------------------------------------------
!
!       0.2  declaration of local variables
!
INTEGER                               :: IIU,IJU,IKU
INTEGER                               :: IIE,IIB,IJE,IJB,IKE,IKB
INTEGER                               :: JI,JJ,JK,JL,JI2,JJ2,JK2,JM
REAL                                  :: ZPH0,ZPH1,ZPH2,ZDEL,ZBAR,ZRAY,ZCOE,ZCO2
REAL                                  :: ZVIT1,ZVIT2,ZVIT0,ZSIG0,ZSIG1,ZSIG2
REAL, DIMENSION(:,:,:,:) ,ALLOCATABLE :: ZIBM_FLUX
REAL, DIMENSION(:,:,:)   ,ALLOCATABLE :: ZIBM_TEST,ZFLU
REAL                                  :: ZADD1,ZADD2,ZADD3,ZADD4,ZADD5,ZADD6
REAL                                  :: ZTOTO
!
!------------------------------------------------------------------------------
!
!       0.3  Allocation
!
CALL GET_DIM_EXT_ll('B',IIU,IJU)
IKU = SIZE(PPHI,3)
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IKE = IKU - JPVEXT
IKB =   1 + JPVEXT
ALLOCATE(ZIBM_FLUX(IIU,IJU,IKU,3))
ALLOCATE(ZFLU(IIU,IJU,IKU))
! 
!------------------------------------------------------------------------------
!
!**** 1. PRELIMINARIES
!     ----------------
!
ZIBM_FLUX = 0.
ZFLU = 0.
ZTOTO = 1.0
!
!------------------------------------------------------------------------------
!
!**** 2. EXECUTIONS
!     -------------
!
 !
  DO JK=IKB,IKE
    DO JJ=IJB,IJE
      DO JI=IIB,IIE

       IF (PVOL(JI,JJ,JK,4).gt.XIBM_EPSI) THEN

        ! Flux Ouest
        JL = 2
        JI2 = JI
        ZIBM_FLUX(JI2,JJ,JK,JL-1) = 0.

        ZDEL = SQRT((XYHAT(JJ+1)-XYHAT(JJ))*0.5*(XZZ(JI2,JJ,JK+1)-XZZ(JI2,JJ,JK)+XZZ(JI2-1,JJ,JK+1)-XZZ(JI2-1,JJ,JK)))
        ZPH1 = PPHI(JI2 ,JJ  ,JK  ,JL)

        ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
        ZVIT1 = ZSIG1*PRUS(JI2,JJ ,JK )

        DO JM=1,8

        IF (JM==1) THEN
        JJ2 = JJ-1
        JK2 = JK-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==2) THEN
        JJ2 = JJ-1
        JK2 = JK
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==3) THEN
        JJ2 = JJ-1
        JK2 = JK+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==4) THEN
        JJ2 = JJ+1
        JK2 = JK-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==5) THEN
        JJ2 = JJ+1
        JK2 = JK
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==6) THEN
        JJ2 = JJ+1
        JK2 = JK+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==7) THEN
        JJ2 = JJ
        JK2 = JK-1
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==8) THEN
        JJ2 = JJ
        JK2 = JK+1
        ZCOE = sqrt(1.)
        ENDIF

        ZPH2 = PPHI(JI2,JJ2,JK2,JL)

        ZSIG2 = max(0.,-ZPH2/abs(ZPH2))
        ZSIG0 = max(0.,-ZPH1*ZPH2/abs(ZPH1*ZPH2))

        ZVIT2 = ZSIG2*PRUS(JI2,JJ2,JK2)

        ZRAY = ZDEL*0.5641895*ZTOTO

        ZBAR = 0.!ZRAY
        IF (ABS(ZPH2-ZPH1).GT.XIBM_EPSI) ZBAR=-ZPH1/(ZPH2-ZPH1)*ZDEL*ZCOE
        ZBAR = MIN(ZBAR,ZRAY)
        ZBAR = MAX(ZBAR,  0.)

        ZIBM_FLUX(JI2,JJ,JK,JL-1) = ZIBM_FLUX(JI2,JJ,JK,JL-1) + &
                                  (ZSIG1*ZSIG2*ZVIT1+ZSIG0*(ZVIT1+ZVIT2)*abs(ZSIG2-(ZBAR/ZRAY)**2.))/8.*ZDEL**2.*ZTOTO**-2.


        ENDDO

        ! Flux Est
        JL = 2
        JI2 = JI+1
        ZIBM_FLUX(JI2,JJ,JK,JL-1) = 0.

        ZDEL = SQRT((XYHAT(JJ+1)-XYHAT(JJ))*0.5*(XZZ(JI2,JJ,JK+1)-XZZ(JI2,JJ,JK)+XZZ(JI2-1,JJ,JK+1)-XZZ(JI2-1,JJ,JK)))
        ZPH1 = PPHI(JI2  ,JJ  ,JK  ,JL)

        ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
        ZVIT1 = ZSIG1*PRUS(JI2,JJ ,JK )

        DO JM=1,8

        IF (JM==1) THEN
        JJ2 = JJ-1
        JK2 = JK-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==2) THEN
        JJ2 = JJ-1
        JK2 = JK
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==3) THEN
        JJ2 = JJ-1
        JK2 = JK+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==4) THEN
        JJ2 = JJ+1
        JK2 = JK-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==5) THEN
        JJ2 = JJ+1
        JK2 = JK
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==6) THEN
        JJ2 = JJ+1
        JK2 = JK+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==7) THEN
        JJ2 = JJ
        JK2 = JK-1
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==8) THEN
        JJ2 = JJ
        JK2 = JK+1
        ZCOE = sqrt(1.)
        ENDIF

        ZPH2 = PPHI(JI2,JJ2,JK2,JL)

        ZSIG2 = max(0.,-ZPH2/abs(ZPH2))
        ZSIG0 = max(0.,-ZPH1*ZPH2/abs(ZPH1*ZPH2))

        ZVIT2 = ZSIG2*PRUS(JI2,JJ2,JK2)

        ZRAY = ZDEL*0.5641895*ZTOTO

        ZBAR = 0.!ZRAY
        IF (ABS(ZPH2-ZPH1).GT.XIBM_EPSI) ZBAR=-ZPH1/(ZPH2-ZPH1)*ZDEL*ZCOE
        ZBAR = MIN(ZBAR,ZRAY)
        ZBAR = MAX(ZBAR,  0.)


        ZIBM_FLUX(JI2,JJ,JK,JL-1) = ZIBM_FLUX(JI2,JJ,JK,JL-1) + &
                                  (ZSIG1*ZSIG2*ZVIT1+ZSIG0*(ZVIT1+ZVIT2)*abs(ZSIG2-(ZBAR/ZRAY)**2.))/8.*ZDEL**2.*ZTOTO**-2.

        ENDDO

        ! Flux Sud
        JL = 3 
        JJ2 = JJ
        ZIBM_FLUX(JI,JJ2,JK,JL-1) = 0.

        ZDEL = SQRT((XXHAT(JI+1)-XXHAT(JI))*0.5*(XZZ(JI,JJ2,JK+1)-XZZ(JI,JJ2,JK)+XZZ(JI,JJ2-1,JK+1)-XZZ(JI,JJ2-1,JK)))
        ZPH1 = PPHI(JI  ,JJ2  ,JK  ,JL)

        ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
        ZVIT1 = ZSIG1*PRVS(JI ,JJ2,JK )

        DO JM=1,8

        IF (JM==1) THEN
        JI2 = JI-1
        JK2 = JK-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==2) THEN
        JI2 = JI-1
        JK2 = JK
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==3) THEN
        JI2 = JI-1
        JK2 = JK+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==4) THEN
        JI2 = JI+1
        JK2 = JK-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==5) THEN
        JI2 = JI+1
        JK2 = JK
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==6) THEN
        JI2 = JI+1
        JK2 = JK+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==7) THEN
        JI2 = JI
        JK2 = JK-1
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==8) THEN
        JI2 = JI
        JK2 = JK+1
        ZCOE = sqrt(1.)
        ENDIF

        ZPH2 = PPHI(JI2,JJ2,JK2,JL)

        ZSIG2 = max(0.,-ZPH2/abs(ZPH2))
        ZSIG0 = max(0.,-ZPH1*ZPH2/abs(ZPH1*ZPH2))

        ZVIT2 = ZSIG2*PRVS(JI2,JJ2,JK2)

        ZRAY = ZDEL*0.5641895*ZTOTO

        ZBAR = 0.!ZRAY
        IF (ABS(ZPH2-ZPH1).GT.XIBM_EPSI) ZBAR = -ZPH1/(ZPH2-ZPH1)*ZDEL*ZCOE
        ZBAR = MIN(ZBAR,ZRAY)
        ZBAR = MAX(ZBAR,  0.)

        ZIBM_FLUX(JI,JJ2,JK,JL-1) = ZIBM_FLUX(JI,JJ2,JK,JL-1) + &
                                  (ZSIG1*ZSIG2*ZVIT1+ZSIG0*(ZVIT1+ZVIT2)*abs(ZSIG2-(ZBAR/ZRAY)**2.))/8.*ZDEL**2.*ZTOTO**-2.
        ENDDO

        ! Flux Nord
        JL = 3 
        JJ2 = JJ+1
        ZIBM_FLUX(JI,JJ2,JK,JL-1) = 0.

        ZDEL = SQRT((XXHAT(JI+1)-XXHAT(JI))*0.5*(XZZ(JI,JJ2,JK+1)-XZZ(JI,JJ2,JK)+XZZ(JI,JJ2-1,JK+1)-XZZ(JI,JJ2-1,JK)))
        ZPH1 = PPHI(JI  ,JJ2  ,JK  ,JL)

        ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
        ZVIT1 = ZSIG1*PRVS(JI ,JJ2,JK )

        DO JM=1,8

        IF (JM==1) THEN
        JI2 = JI-1
        JK2 = JK-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==2) THEN
        JI2 = JI-1
        JK2 = JK
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==3) THEN
        JI2 = JI-1
        JK2 = JK+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==4) THEN
        JI2 = JI+1
        JK2 = JK-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==5) THEN
        JI2 = JI+1
        JK2 = JK
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==6) THEN
        JI2 = JI+1
        JK2 = JK+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==7) THEN
        JI2 = JI
        JK2 = JK-1
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==8) THEN
        JI2 = JI
        JK2 = JK+1
        ZCOE = sqrt(1.)
        ENDIF

        ZPH2 = PPHI(JI2,JJ2,JK2,JL)

        ZSIG2 = max(0.,-ZPH2/abs(ZPH2))
        ZSIG0 = max(0.,-ZPH1*ZPH2/abs(ZPH1*ZPH2))

        ZVIT2 = ZSIG2*PRVS(JI2,JJ2,JK2)

        ZRAY = ZDEL*0.5641895*ZTOTO

        ZBAR = 0.!ZRAY
        IF (ABS(ZPH2-ZPH1).GT.XIBM_EPSI) ZBAR = -ZPH1/(ZPH2-ZPH1)*ZDEL*ZCOE
        ZBAR = MIN(ZBAR,ZRAY)
        ZBAR = MAX(ZBAR,  0.)

        ZIBM_FLUX(JI,JJ2,JK,JL-1) = ZIBM_FLUX(JI,JJ2,JK,JL-1) + &
                                  (ZSIG1*ZSIG2*ZVIT1+ZSIG0*(ZVIT1+ZVIT2)*abs(ZSIG2-(ZBAR/ZRAY)**2.))/8.*ZDEL**2.*ZTOTO**-2.


        ENDDO

        ! Flux bas
        JL = 4 
        JK2 = JK
        ZIBM_FLUX(JI,JJ,JK2,JL-1) = 0.

        ZDEL = SQRT((XXHAT(JI+1)-XXHAT(JI))*(XYHAT(JJ+1)-XYHAT(JJ)))
        ZPH1 = PPHI(JI  ,JJ  ,JK2 ,JL)

        ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
        ZVIT1 = ZSIG1*PRWS(JI ,JJ ,JK2)

        DO JM=1,8

        IF (JM==1) THEN
        JJ2 = JJ-1
        JI2 = JI-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==2) THEN
        JJ2 = JJ-1
        JI2 = JI
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==3) THEN
        JJ2 = JJ-1
        JI2 = JI+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==4) THEN
        JJ2 = JJ+1
        JI2 = JI-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==5) THEN
        JJ2 = JJ+1
        JI2 = JI
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==6) THEN
        JJ2 = JJ+1
        JI2 = JI+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==7) THEN
        JJ2 = JJ
        JI2 = JI-1
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==8) THEN
        JJ2 = JJ
        JI2 = JI+1
        ZCOE = sqrt(1.)
        ENDIF

        ZPH2 = PPHI(JI2,JJ2,JK2,JL)

        ZSIG2 = max(0.,-ZPH2/abs(ZPH2))
        ZSIG0 = max(0.,-ZPH1*ZPH2/abs(ZPH1*ZPH2))

        ZVIT2 = ZSIG2*PRWS(JI2,JJ2,JK2)

        ZRAY = ZDEL*0.5641895*ZTOTO

        ZBAR = 0.!ZRAY
        IF (ABS(ZPH2-ZPH1).GT.XIBM_EPSI) ZBAR = -ZPH1/(ZPH2-ZPH1)*ZDEL*ZCOE
        ZBAR = MIN(ZBAR,ZRAY)
        ZBAR = MAX(ZBAR,  0.)

        ZIBM_FLUX(JI,JJ,JK2,JL-1) = ZIBM_FLUX(JI,JJ,JK2,JL-1) + &
                                  (ZSIG1*ZSIG2*ZVIT1+ZSIG0*(ZVIT1+ZVIT2)*abs(ZSIG2-(ZBAR/ZRAY)**2.))/8.*ZDEL**2.*ZTOTO**-2.

        ENDDO

        ! Flux haut
        JL = 4 
        JK2 = JK+1
        ZIBM_FLUX(JI,JJ,JK2,JL-1) = 0.

        ZDEL = SQRT((XXHAT(JI+1)-XXHAT(JI))*(XYHAT(JJ+1)-XYHAT(JJ)))
        ZPH1 = PPHI(JI  ,JJ  ,JK2 ,JL)

        ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
        ZVIT1 = ZSIG1*PRWS(JI ,JJ ,JK2)

        DO JM=1,8

        IF (JM==1) THEN
        JJ2 = JJ-1
        JI2 = JI-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==2) THEN
        JJ2 = JJ-1
        JI2 = JI
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==3) THEN
        JJ2 = JJ-1
        JI2 = JI+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==4) THEN
        JJ2 = JJ+1
        JI2 = JI-1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==5) THEN
        JJ2 = JJ+1
        JI2 = JI
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==6) THEN
        JJ2 = JJ+1
        JI2 = JI+1
        ZCOE = sqrt(2.)
        ENDIF
        IF (JM==7) THEN
        JJ2 = JJ
        JI2 = JI-1
        ZCOE = sqrt(1.)
        ENDIF
        IF (JM==8) THEN
        JJ2 = JJ
        JI2 = JI+1
        ZCOE = sqrt(1.)
        ENDIF

        ZPH2 = PPHI(JI2,JJ2,JK2,JL)

        ZSIG2 = max(0.,-ZPH2/abs(ZPH2))
        ZSIG0 = max(0.,-ZPH1*ZPH2/abs(ZPH1*ZPH2))

        ZVIT2 = ZSIG2*PRWS(JI2,JJ2,JK2)

        ZRAY = ZDEL*0.5641895*ZTOTO

        ZBAR = 0.!ZRAY
        IF (ABS(ZPH2-ZPH1).GT.XIBM_EPSI) ZBAR = -ZPH1/(ZPH2-ZPH1)*ZDEL*ZCOE
        ZBAR = MIN(ZBAR,ZRAY)
        ZBAR = MAX(ZBAR,  0.)

        ZIBM_FLUX(JI,JJ,JK2,JL-1) = ZIBM_FLUX(JI,JJ,JK2,JL-1) + &
                                  (ZSIG1*ZSIG2*ZVIT1+ZSIG0*(ZVIT1+ZVIT2)*abs(ZSIG2-(ZBAR/ZRAY)**2.))/8.*ZDEL**2.*ZTOTO**-2.

        ENDDO

        ENDIF

      ENDDO
    ENDDO
  ENDDO
  !

  ZFLU(IIB:IIE,IJB:IJE,IKB:IKE) = (ZIBM_FLUX(IIB+1:IIE+1,IJB  :IJE  ,IKB  :IKE  ,1)-ZIBM_FLUX(IIB:IIE,IJB:IJE,IKB:IKE,1) +&
                                   ZIBM_FLUX(IIB  :IIE  ,IJB+1:IJE+1,IKB  :IKE  ,2)-ZIBM_FLUX(IIB:IIE,IJB:IJE,IKB:IKE,2) +&
                                   ZIBM_FLUX(IIB  :IIE  ,IJB  :IJE  ,IKB+1:IKE+1,3)-ZIBM_FLUX(IIB:IIE,IJB:IJE,IKB:IKE,3))*&
                                   XRHODREF(IIB:IIE,IJB:IJE,IKB:IKE)/XRHODJ(IIB:IIE,IJB:IJE,IKB:IKE)

  PBAL(IIB-1:IIE+1,IJB-1:IJE+1,IKB-1:IKE+1) = PBAL(IIB-1:IIE+1,IJB-1:IJE+1,IKB-1:IKE+1)*PVOL(IIB-1:IIE+1,IJB-1:IJE+1,IKB-1:IKE+1,2)+&
                                              ZFLU(IIB-1:IIE+1,IJB-1:IJE+1,IKB-1:IKE+1)*PVOL(IIB-1:IIE+1,IJB-1:IJE+1,IKB-1:IKE+1,4)
!
!**** X. DEALLOCATIONS/CLOSES
!     -----------------------
DEALLOCATE(ZIBM_FLUX,ZFLU)
! 
RETURN
!
END SUBROUTINE IBM_BALANCE