!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.
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source: /srv/cvsroot/MNH-VX-Y-Z/src/MNH/ini_cturb.f90,v $ $Revision: 1.2.4.1.2.2.2.1.2.2 $ $Date: 2014/01/09 15:01:55 $
!-----------------------------------------------------------------
!-----------------------------------------------------------------
!       #######################
        MODULE MODI_IBM_FORCING2
!       #######################
!
INTERFACE
!
SUBROUTINE IBM_FORCING2(PRUS,PRVS,PRWS,PTHS,PRRS,PSVS,PTKS)
!
REAL, DIMENSION(:,:,:)     ,INTENT(INOUT) :: PRUS,PRVS,PRWS
REAL, DIMENSION(:,:,:)     ,INTENT(INOUT) :: PTHS
REAL, DIMENSION(:,:,:,:)   ,INTENT(INOUT) :: PRRS
REAL, DIMENSION(:,:,:,:)   ,INTENT(INOUT) :: PSVS
REAL, DIMENSION(:,:,:)     ,INTENT(INOUT) :: PTKS

!
END SUBROUTINE IBM_FORCING2
!
END INTERFACE
!
END MODULE MODI_IBM_FORCING2
!
!
!
!       ##########################################################
        SUBROUTINE IBM_FORCING2(PRUS,PRVS,PRWS,PTHS,PRRS,PSVS,PTKS)
!       ##########################################################
!
!!****     *IBM_FORCING2*  - routine to force all desired fields 
!!
!!      PURPOSE
!!      -------
!         The purpose of this routine is to compute variables in the virtual
!         embedded solid region in regard of variables computed in the real
!         fluid region when the classical GCT induces artefacts
!         (if necessary, see the output file messages)
!
!!      METHOD
!!      ------
!!
!!      EXTERNAL
!!      --------
!!        NONE
!!
!!      IMPLICIT ARGUMENTS
!!      ------------------
!!
!!      REFERENCE
!!      ---------
!!
!!      AUTHOR
!!      ------
!!        Franck Auguste       * CERFACS(AE) *
!!
!!      MODIFICATIONS
!!      -------------
!!        Original          01/01/2015
!!
!------------------------------------------------------------------------------
!       
!**** 0. DECLARATIONS
!     ---------------
!
! module
USE MODE_POS
USE MODE_ll
USE MODE_IO_ll
USE MODD_ARGSLIST_ll,     ONLY : LIST_ll
!
! declaration
USE MODD_REF_n,           ONLY: XRHODJ,XRHODREF,XTHVREF
USE MODD_PARAMETERS,      ONLY: JPVEXT,JPHEXT
USE MODD_IBM_PARAM_n
USE MODD_LBC_n
USE MODD_CONF
USE MODD_CONF_n
USE MODD_NSV
USE MODD_CTURB
USE MODD_PARAM_n
!
! interface
!
IMPLICIT NONE
!
!------------------------------------------------------------------------------
!
!       0.1  declarations of arguments                                                               
REAL, DIMENSION(:,:,:)     ,INTENT(INOUT) :: PRUS,PRVS,PRWS
REAL, DIMENSION(:,:,:)     ,INTENT(INOUT) :: PTHS
REAL, DIMENSION(:,:,:,:)   ,INTENT(INOUT) :: PRRS
REAL, DIMENSION(:,:,:,:)   ,INTENT(INOUT) :: PSVS
REAL, DIMENSION(:,:,:)     ,INTENT(INOUT) :: PTKS

!
!------------------------------------------------------------------------------
!
!       0.2  declaration of local variables
INTEGER                                 :: JI,JJ,JK,JI2,JJ2,JK2,IIU,IJU,IKU
INTEGER                                 :: JIM1,JJM1,JKM1,JIP1,JJP1,JKP1
INTEGER                                 :: IIE,IIB,IJE,IJB,IKB,IKE
REAL                                    :: ZSUM1,ZSUM2,ZSUM3,ZSUM4
TYPE(LIST_ll), POINTER                  :: TZFIELDS_ll   
INTEGER                                 :: IINFO_ll
!
!------------------------------------------------------------------------------
!
RETURN
!**** 0. ALLOCATIONS
!     --------------
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IIU = SIZE(PRUS,1)
IJU = SIZE(PRUS,2)
IKU = SIZE(PRUS,3)
IKB =   1 + JPVEXT
IKE = SIZE(PRUS,3) - JPVEXT
!
!------------------------------------------------------------------------------
! 
  DO JI=IIB,IIE
  DO JJ=IJB,IJE
  DO JK=IKB,IKE

  IF (XIBM_SU(JI,JJ,JK,6).LT.0.5) THEN

    JIM1 = JI-1
    JJM1 = JJ-1
    JKM1 = JK-1
    JIP1 = JI+1
    JJP1 = JJ+1
    JKP1 = JK+1

    ZSUM1 = 0.
    ZSUM2 = 0.
    ZSUM3 = 0.
    ZSUM4 = 0.

    DO JI2=JIM1,JIP1
    DO JJ2=JJM1,JJP1
    DO JK2=JKM1,JKP1

                         ZSUM1 = ZSUM1 + (XIBM_SU(JI2,JJ2,JK2,6))
                         ZSUM2 = ZSUM2 + (XIBM_SU(JI2,JJ2,JK2,6))*PTHS(JI2,JJ2,JK2)
             IF (NSV==1) ZSUM3 = ZSUM3 + (XIBM_SU(JI2,JJ2,JK2,6))*PSVS(JI2,JJ2,JK2,1)
      IF (CTURB/='NONE') ZSUM4 = ZSUM4 + (XIBM_SU(JI2,JJ2,JK2,6))*PTKS(JI2,JJ2,JK2)

    ENDDO
    ENDDO
    ENDDO

    PTHS(JI,JJ,JK) = XIBM_THT(JK)
    IF (NSV==1) PSVS(JI,JJ,JK,1) = 0.
    IF (CTURB/='NONE') PTKS(JI,JJ,JK) = XTKEMIN
    IF (ZSUM1.GT.XIBM_EPSI)     PTHS(JI,JJ,JK) =  ZSUM2/ZSUM1
    IF (ZSUM1.GT.XIBM_EPSI.AND.NSV==1) PSVS(JI,JJ,JK,1) = ZSUM3/ZSUM1
    IF (ZSUM1.GT.XIBM_EPSI.AND.(CTURB/='NONE'))   PTKS(JI,JJ,JK) = ZSUM4/ZSUM1

  ENDIF

  IF (XIBM_SU(JI,JJ,JK,7).LT.0.5) THEN

    JIM1 = JI-1
    JJM1 = JJ-1
    JKM1 = JK-1
    JIP1 = JI+1
    JJP1 = JJ+1
    JKP1 = JK+1

    ZSUM1 = 0.
    ZSUM2 = 0.

    DO JI2=JIM1,JIP1
    DO JJ2=JJM1,JJP1
    DO JK2=JKM1,JKP1

                         ZSUM1 = ZSUM1 + (XIBM_SU(JI2,JJ2,JK2,7))
                         ZSUM2 = ZSUM2 + (XIBM_SU(JI2,JJ2,JK2,7))*PRUS(JI2,JJ2,JK2)

    ENDDO
    ENDDO
    ENDDO

    PRUS(JI,JJ,JK) = 0.
    IF (ZSUM1.GT.XIBM_EPSI)     PRUS(JI,JJ,JK) =  ZSUM2/ZSUM1

  ENDIF

  IF (XIBM_SU(JI,JJ,JK,8).LT.0.5) THEN

    JIM1 = JI-1
    JJM1 = JJ-1
    JKM1 = JK-1
    JIP1 = JI+1
    JJP1 = JJ+1
    JKP1 = JK+1

    ZSUM1 = 0.
    ZSUM2 = 0.

    DO JI2=JIM1,JIP1
    DO JJ2=JJM1,JJP1
    DO JK2=JKM1,JKP1

                         ZSUM1 = ZSUM1 + (XIBM_SU(JI2,JJ2,JK2,8))
                         ZSUM2 = ZSUM2 + (XIBM_SU(JI2,JJ2,JK2,8))*PRVS(JI2,JJ2,JK2)

    ENDDO
    ENDDO
    ENDDO

    PRVS(JI,JJ,JK) = 0.
    IF (ZSUM1.GT.XIBM_EPSI)     PRVS(JI,JJ,JK) =  ZSUM2/ZSUM1

  ENDIF

  IF (XIBM_SU(JI,JJ,JK,9).LT.0.5) THEN

    JIM1 = JI-1
    JJM1 = JJ-1
    JKM1 = JK-1
    JIP1 = JI+1
    JJP1 = JJ+1
    JKP1 = JK+1

    ZSUM1 = 0.
    ZSUM2 = 0.

    DO JI2=JIM1,JIP1
    DO JJ2=JJM1,JJP1
    DO JK2=JKM1,JKP1

                         ZSUM1 = ZSUM1 + (XIBM_SU(JI2,JJ2,JK2,9))
                         ZSUM2 = ZSUM2 + (XIBM_SU(JI2,JJ2,JK2,9))*PRWS(JI2,JJ2,JK2)

    ENDDO
    ENDDO
    ENDDO

    PRWS(JI,JJ,JK) = 0.
    IF (ZSUM1.GT.XIBM_EPSI)     PRWS(JI,JJ,JK) =  ZSUM2/ZSUM1

  ENDIF

  ENDDO
  ENDDO
  ENDDO

  PTHS(:,:,IKB-1)=PTHS(:,:,IKB)
  PTHS(:,:,IKE+1)=PTHS(:,:,IKE)
  IF (CTURB/='NONE') PTKS(:,:,IKB-1)=PTKS(:,:,IKB)
  IF (CTURB/='NONE')  PTKS(:,:,IKE+1)=PTKS(:,:,IKE)
  IF (NSV==1) PSVS(:,:,IKB-1,1)=PSVS(:,:,IKB,1)
  IF (NSV==1) PSVS(:,:,IKE+1,1)=PSVS(:,:,IKE,1)
  PRUS(:,:,IKB-1)=PRUS(:,:,IKB)
  PRUS(:,:,IKE+1)=PRUS(:,:,IKE)
  PRVS(:,:,IKB-1)=PRVS(:,:,IKB)
  PRVS(:,:,IKE+1)=PRVS(:,:,IKE)
  PRWS(:,:,IKB-1)=-PRWS(:,:,IKB+1)
  PRWS(:,:,IKB)  =0.
  PRWS(:,:,IKE+1)=0.

  !
  NULLIFY(TZFIELDS_ll)
  CALL ADD3DFIELD_ll(TZFIELDS_ll,PTHS(:,:,:))
  IF (CTURB/='NONE') CALL ADD3DFIELD_ll(TZFIELDS_ll,PTKS(:,:,:))
  CALL ADD3DFIELD_ll(TZFIELDS_ll,PRUS(:,:,:))
  CALL ADD3DFIELD_ll(TZFIELDS_ll,PRVS(:,:,:))
  CALL ADD3DFIELD_ll(TZFIELDS_ll,PRWS(:,:,:))
  IF (NSV==1) CALL ADD3DFIELD_ll(TZFIELDS_ll,PSVS(:,:,:,1))
  CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
  CALL CLEANLIST_ll(TZFIELDS_ll)
!
  DO JI=IIB,IIE
  DO JJ=IJB,IJE
  DO JK=IKB,IKE

  IF (XIBM_LS(JI,JJ,JK,2).GT.XIBM_EPSI) THEN

    ZSUM1 = (XIBM_SU(JI,JJ,JK,5)+XIBM_SU(JI-1,JJ,JK,5))/2.
    ZSUM1 = ABS(ZSUM1)
    ZSUM1 = MIN(1.,ZSUM1)
   
    JIM1 = JI-1
    JJM1 = JJ-1
    JKM1 = JK-1
    JIP1 = JI+1
    JJP1 = JJ+1
    JKP1 = JK+1

    ZSUM2 = 0.

    DO JI2=JIM1,JIP1
    DO JJ2=JJM1,JJP1
    DO JK2=JKM1,JKP1

                         ZSUM2 = ZSUM2 + PRUS(JI2,JJ2,JK2)

    ENDDO
    ENDDO
    ENDDO

    PRUS(JI,JJ,JK) =   (1.-ZSUM1)*PRUS(JI,JJ,JK)+ZSUM1*ZSUM2/27.

  ENDIF

  IF (XIBM_LS(JI,JJ,JK,3).GT.XIBM_EPSI) THEN

    ZSUM1 = (XIBM_SU(JI,JJ,JK,5)+XIBM_SU(JI,JJ-1,JK,5))/2.
    ZSUM1 = ABS(ZSUM1)
    ZSUM1 = MIN(1.,ZSUM1)
   
    JIM1 = JI-1
    JJM1 = JJ-1
    JKM1 = JK-1
    JIP1 = JI+1
    JJP1 = JJ+1
    JKP1 = JK+1

    ZSUM2 = 0.

    DO JI2=JIM1,JIP1
    DO JJ2=JJM1,JJP1
    DO JK2=JKM1,JKP1

                         ZSUM2 = ZSUM2 + PRVS(JI2,JJ2,JK2)

    ENDDO
    ENDDO
    ENDDO

    PRVS(JI,JJ,JK) =   (1.-ZSUM1)*PRVS(JI,JJ,JK)+ZSUM1*ZSUM2/27.

  ENDIF

  IF (XIBM_LS(JI,JJ,JK,4).GT.XIBM_EPSI) THEN

    ZSUM1 = (XIBM_SU(JI,JJ,JK,5)+XIBM_SU(JI,JJ,JK-1,5))/2.
    ZSUM1 = ABS(ZSUM1)
    ZSUM1 = MIN(1.,ZSUM1)
   
    JIM1 = JI-1
    JJM1 = JJ-1
    JKM1 = JK-1
    JIP1 = JI+1
    JJP1 = JJ+1
    JKP1 = JK+1

    ZSUM2 = 0.

    DO JI2=JIM1,JIP1
    DO JJ2=JJM1,JJP1
    DO JK2=JKM1,JKP1

                         ZSUM2 = ZSUM2 + PRWS(JI2,JJ2,JK2)

    ENDDO
    ENDDO
    ENDDO

    PRWS(JI,JJ,JK) =   (1.-ZSUM1)*PRWS(JI,JJ,JK)+ZSUM1*ZSUM2/27.

  ENDIF
  ENDDO
  ENDDO
  ENDDO

  PTHS(:,:,IKB-1)=PTHS(:,:,IKB)
  PTHS(:,:,IKE+1)=PTHS(:,:,IKE)
  IF (CTURB/='NONE') PTKS(:,:,IKB-1)=PTKS(:,:,IKB)
  IF (CTURB/='NONE')  PTKS(:,:,IKE+1)=PTKS(:,:,IKE)
  IF (NSV==1) PSVS(:,:,IKB-1,1)=PSVS(:,:,IKB,1)
  IF (NSV==1) PSVS(:,:,IKE+1,1)=PSVS(:,:,IKE,1)
  PRUS(:,:,IKB-1)=PRUS(:,:,IKB)
  PRUS(:,:,IKE+1)=PRUS(:,:,IKE)
  PRVS(:,:,IKB-1)=PRVS(:,:,IKB)
  PRVS(:,:,IKE+1)=PRVS(:,:,IKE)
  PRWS(:,:,IKB-1)=0.
  PRWS(:,:,IKB)  =0.
  PRWS(:,:,IKE+1)=0.

  !
  NULLIFY(TZFIELDS_ll)
  CALL ADD3DFIELD_ll(TZFIELDS_ll,PTHS(:,:,:))
  IF (CTURB/='NONE') CALL ADD3DFIELD_ll(TZFIELDS_ll,PTKS(:,:,:))
  CALL ADD3DFIELD_ll(TZFIELDS_ll,PRUS(:,:,:))
  CALL ADD3DFIELD_ll(TZFIELDS_ll,PRVS(:,:,:))
  CALL ADD3DFIELD_ll(TZFIELDS_ll,PRWS(:,:,:))
  IF (NSV==1) CALL ADD3DFIELD_ll(TZFIELDS_ll,PSVS(:,:,:,1))
  CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
  CALL CLEANLIST_ll(TZFIELDS_ll)

RETURN

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END SUBROUTINE IBM_FORCING2