Skip to content
Snippets Groups Projects
ice_adjust.f90 32 KiB
Newer Older
!MNH_LIC Copyright 1996-2022 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_ICE_ADJUST
!     ######################
!
INTERFACE
!
      SUBROUTINE ICE_ADJUST (KKA, KKU, KKL, KRR, HFRAC_ICE, HCONDENS, HLAMBDA3,&
                             HBUNAME, OSUBG_COND, OSIGMAS, HSUBG_MF_PDF,       &
                             PTSTEP, PSIGQSAT,                                 &
                             PRHODJ, PEXNREF, PRHODREF, PSIGS, PMFCONV,        &
                             PPABST, PZZ,                                      &
                             PEXN, PCF_MF, PRC_MF, PRI_MF,                     &
                             PRV, PRC, PRVS, PRCS, PTH, PTHS, PSRCS, PCLDFR,   &
                             PRR, PRI, PRIS, PRS, PRG, PRH,                    &
                             POUT_RV, POUT_RC, POUT_RI, POUT_TH,               &
                             PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF            )
!
INTEGER,                  INTENT(IN)    :: KKA   !near ground array index
INTEGER,                  INTENT(IN)    :: KKU   !uppest atmosphere array index
INTEGER,                  INTENT(IN)    :: KKL   !vert. levels type 1=MNH -1=ARO
INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
CHARACTER(len=1),         INTENT(IN)    :: HFRAC_ICE
CHARACTER(len=80),        INTENT(IN)    :: HCONDENS
CHARACTER(len=4),         INTENT(IN)    :: HLAMBDA3 ! formulation for lambda3 coeff
CHARACTER(len=*),         INTENT(IN)    :: HBUNAME  ! Name of the budget
LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
                                                    ! Condensation
LOGICAL                                 :: OSIGMAS  ! Switch for Sigma_s: 
                                                    ! use values computed in CONDENSATION
                                                    ! or that from turbulence scheme
CHARACTER(len=*),         INTENT(IN)   :: HSUBG_MF_PDF
REAL,                     INTENT(IN)   :: PTSTEP    ! Double Time step
                                                    ! (single if cold start)
REAL,                     INTENT(IN)   :: PSIGQSAT  ! coeff applied to qsat variance contribution
!
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODJ  ! Dry density * Jacobian
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PEXNREF ! Reference Exner function
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODREF
!
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PSIGS   ! Sigma_s at time t
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PMFCONV ! convective mass flux
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PPABST  ! Absolute Pressure at t
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PZZ     ! height of model layer
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PEXN    ! Exner function
!
REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PCF_MF! Convective Mass Flux Cloud fraction
REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRI_MF! Convective Mass Flux ice mixing ratio
REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRC_MF! Convective Mass Flux liquid mixing ratio
!
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRV     ! Water vapor m.r. to adjust
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRC     ! Cloud water m.r. to adjust
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRVS    ! Water vapor m.r. source
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRCS    ! Cloud water m.r. source
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTH     ! Theta to adjust
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PSRCS   ! Second-order flux
                                                   ! s'rc'/2Sigma_s2 at time t+1
                                                   ! multiplied by Lambda_3
REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PCLDFR  ! Cloud fraction
REAL, DIMENSION(:,:,:),   INTENT(INOUT) ::  PRIS   ! Cloud ice  m.r. at t+1
REAL, DIMENSION(:,:,:),   INTENT(IN)    ::  PRR    ! Rain water m.r. to adjust
REAL, DIMENSION(:,:,:),   INTENT(IN)    ::  PRI    ! Cloud ice  m.r. to adjust
REAL, DIMENSION(:,:,:),   INTENT(IN)    ::  PRS    ! Aggregate  m.r. to adjust
REAL, DIMENSION(:,:,:),   INTENT(IN)    ::  PRG    ! Graupel    m.r. to adjust
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN)   ::  PRH  ! Hail       m.r. to adjust
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  POUT_RV ! Adjusted value
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  POUT_RC ! Adjusted value
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  POUT_RI ! Adjusted value
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  POUT_TH ! Adjusted value
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  PHLC_HRC
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  PHLC_HCF
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  PHLI_HRI
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  PHLI_HCF
!
!
END SUBROUTINE ICE_ADJUST
!
END INTERFACE
!
END MODULE MODI_ICE_ADJUST

!     ##########################################################################
      SUBROUTINE ICE_ADJUST (KKA, KKU, KKL, KRR, HFRAC_ICE, HCONDENS, HLAMBDA3,&
                             HBUNAME, OSUBG_COND, OSIGMAS, HSUBG_MF_PDF,       &
                             PTSTEP, PSIGQSAT,                                 &
                             PRHODJ, PEXNREF, PRHODREF, PSIGS, PMFCONV,        &
                             PPABST, PZZ,                                      &
                             PEXN, PCF_MF, PRC_MF, PRI_MF,                     &
                             PRV, PRC, PRVS, PRCS, PTH, PTHS, PSRCS, PCLDFR,   &
                             PRR, PRI, PRIS, PRS, PRG, PRH,                    &
                             POUT_RV, POUT_RC, POUT_RI, POUT_TH,               &
                             PHLC_HRC, PHLC_HCF, PHLI_HRI, PHLI_HCF            )
!     #########################################################################
!
!!****  *ICE_ADJUST* -  compute the ajustment of water vapor in mixed-phase 
!!                      clouds
!!
!!    PURPOSE
!!    -------
!!    The purpose of this routine is to compute the fast microphysical sources
!!    through a saturation ajustement procedure in case of mixed-phase clouds.
!!
!!
!!**  METHOD
!!    ------
!!    Langlois, Tellus, 1973 for the cloudless version.
!!    When cloud water is taken into account, refer to book 1 of the
!!    documentation.
!!
!!     
!!
!!    EXTERNAL
!!    --------
!!      None
!!     
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!      Module MODD_CST
!!         XP00               ! Reference pressure
!!         XMD,XMV            ! Molar mass of dry air and molar mass of vapor
!!         XRD,XRV            ! Gaz constant for dry air, gaz constant for vapor
!!         XCPD,XCPV          ! Cpd (dry air), Cpv (vapor)
!!         XCL                ! Cl (liquid)
!!         XCI                ! Ci (ice)
!!         XTT                ! Triple point temperature
!!         XLVTT              ! Vaporization heat constant
!!         XLSTT              ! Sublimation  heat constant
!!         XALPW,XBETAW,XGAMW ! Constants for saturation vapor over liquid
!!                            !  pressure  function 
!!         XALPI,XBETAI,XGAMI ! Constants for saturation vapor over ice
!!                            !  pressure  function 
!!      Module  MODD_CONF 
!!         CCONF
!!      Module MODD_BUDGET:
!!         NBUMOD 
!!         CBUTYPE
!!         LBU_RTH    
!!         LBU_RRV  
!!         LBU_RRC  
!!         LBU_RRI  
!!
!!
!!    REFERENCE
!!    ---------
!!      Book 1 and Book2 of documentation ( routine ICE_ADJUST )
!!      Langlois, Tellus, 1973
!!
!!    AUTHOR
!!    ------
!!      J.-P. Pinty    * Laboratoire d'Aerologie*
!!   
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    06/12/96 
!!      M. Tomasini 27/11/00 Change CND and DEP fct of the T instead of rc and ri
!!                           Avoid the sub- and super-saturation before the ajustment
!!                           Avoid rc>0 if T<T00 before the ajustment
!!      P Bechtold 12/02/02  change subgrid condensation
!!      JP Pinty   29/11/02  add ICE2 and IC4 cases
!!      (P. Jabouille) 27/05/04 safety test for case where esw/i(T)> pabs (~Z>40km)
!!      J.Pergaud and S.Malardel Add EDKF case
!!      S. Riette ice for EDKF
!!      2012-02 Y. Seity,  add possibility to run with reversed vertical levels
!!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
!!      2016-07 S. Riette: adjustement is now realized on state variables (PRV, PRC, PRI, PTH)
!!                         whereas tendencies are still applied on S variables.
!!                         This modification allows to call ice_adjust on T variable
!!                         or to call it on S variables
!!      2016-11 S. Riette: all-or-nothing adjustment now uses condensation
!  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
!  P. Wautelet    02/2020: use the new data structures and subroutines for budgets
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
use modd_budget,       only: lbudget_th, lbudget_rv, lbudget_rc, lbudget_ri, &
                             NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, NBUDGET_RI, &
                             tbudgets
USE MODD_CONF
USE MODD_CST
USE MODD_PARAMETERS
USE MODD_RAIN_ICE_PARAM, ONLY : XCRIAUTC, XCRIAUTI, XACRIAUTI, XBCRIAUTI

use mode_budget,       only: Budget_store_init, Budget_store_end
#ifdef MNH_OPENACC
USE MODE_MNH_ZWORK,    ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
#endif
USE MODE_MPPDB
#ifdef MNH_OPENACC
use mode_msg
#endif
use mode_tools_ll,        only: GET_INDICE_ll

#ifdef MNH_BITREP
USE MODI_BITREP
#endif
USE MODI_CONDENSATION
USE MODI_GET_HALO

IMPLICIT NONE
!
!
!*       0.1   Declarations of dummy arguments :
!
!
INTEGER,                  INTENT(IN)    :: KKA  !near ground array index
INTEGER,                  INTENT(IN)    :: KKU  !uppest atmosphere array index
INTEGER,                  INTENT(IN)    :: KKL  !vert. levels type 1=MNH -1=ARO
INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
CHARACTER(len=1),         INTENT(IN)    :: HFRAC_ICE
CHARACTER(len=80),        INTENT(IN)    :: HCONDENS
CHARACTER(len=4),         INTENT(IN)    :: HLAMBDA3 ! formulation for lambda3 coeff
CHARACTER(len=*),         INTENT(IN)    :: HBUNAME  ! Name of the budget
LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
                                                    ! Condensation
LOGICAL                                 :: OSIGMAS  ! Switch for Sigma_s: 
                                                    ! use values computed in CONDENSATION
                                                    ! or that from turbulence scheme
CHARACTER(len=*),         INTENT(IN)    :: HSUBG_MF_PDF
REAL,                     INTENT(IN)   :: PTSTEP    ! Double Time step
                                                    ! (single if cold start)
REAL,                     INTENT(IN)   :: PSIGQSAT  ! coeff applied to qsat variance contribution
!
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODJ  ! Dry density * Jacobian
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PEXNREF ! Reference Exner function
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODREF
!
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PSIGS   ! Sigma_s at time t
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PMFCONV ! convective mass flux
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PPABST  ! Absolute Pressure at t
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PZZ     ! height of model layer
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PEXN    ! Exner function
!
REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PCF_MF! Convective Mass Flux Cloud fraction
REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRC_MF! Convective Mass Flux liquid mixing ratio
REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRI_MF! Convective Mass Flux ice mixing ratio
!
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRV     ! Water vapor m.r. to adjust
REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRC     ! Cloud water m.r. to adjust
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRVS    ! Water vapor m.r. source
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRCS    ! Cloud water m.r. source
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTH     ! Theta to adjust
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PSRCS   ! Second-order flux
                                                   ! s'rc'/2Sigma_s2 at time t+1
                                                   ! multiplied by Lambda_3
REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PCLDFR  ! Cloud fraction
!
REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRIS ! Cloud ice  m.r. at t+1
REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRR  ! Rain water m.r. to adjust
REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRI  ! Cloud ice  m.r. to adjust
REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRS  ! Aggregate  m.r. to adjust
REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRG  ! Graupel    m.r. to adjust
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN)   ::  PRH  ! Hail       m.r. to adjust
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  POUT_RV ! Adjusted value
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  POUT_RC ! Adjusted value
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  POUT_RI ! Adjusted value
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  POUT_TH ! Adjusted value
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  PHLC_HRC
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  PHLC_HCF
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  PHLI_HRI
REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(OUT)  ::  PHLI_HCF
!
!*       0.2   Declarations of local variables :
!
!
INTEGER             :: IIU,IJU,IKU! dimensions of dummy arrays
INTEGER             :: IIB,IJB    ! Horz index values of the first inner mass points
INTEGER             :: IIE,IJE    ! Horz index values of the last inner mass points
INTEGER             :: IKB        ! K index value of the first inner mass point
INTEGER             :: IKE        ! K index value of the last inner mass point
INTEGER             :: JITER,ITERMAX ! iterative loop for first order adjustment
INTEGER             :: JI,JJ,JK
!
#ifndef MNH_OPENACC
LOGICAL,DIMENSION(:,:,:), allocatable :: GTEMP
!
REAL, DIMENSION(:,:,:), allocatable :: ZSIGS,ZSRCS
REAL, DIMENSION(:,:,:), allocatable &
                         :: ZT,   &  ! adjusted temperature
                   ZRV, ZRC, ZRI, &  ! adjusted state
                            ZCPH, &  ! guess of the CPh for the mixing
                            ZLV,  &  ! guess of the Lv at t+1
                            ZLS,  &  ! guess of the Ls at t+1
                         ZW1,ZW2, &  ! Work arrays for intermediate fields
                         ZCRIAUT, &  ! Autoconversion thresholds
                       ZHCF, ZHR
#else
LOGICAL,DIMENSION(:,:,:), POINTER, CONTIGUOUS :: GTEMP
!
REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZSIGS,ZSRCS
REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS &
                         :: ZT,   &  ! adjusted temperature
                   ZRV, ZRC, ZRI, &  ! adjusted state
                            ZCPH, &  ! guess of the CPh for the mixing
                            ZLV,  &  ! guess of the Lv at t+1
                            ZLS,  &  ! guess of the Ls at t+1
                         ZW1,ZW2, &  ! Work arrays for intermediate fields
                         ZCRIAUT, &  ! Autoconversion thresholds
                         ZHCF, ZHR
REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTEMP_BUD
LOGICAL :: GPOUT_RV,GPOUT_RC,GPOUT_RI,GPOUT_TH
!-------------------------------------------------------------------------------
!
GPOUT_RV = PRESENT(POUT_RV)
GPOUT_RC = PRESENT(POUT_RC)
GPOUT_RI = PRESENT(POUT_RI)
GPOUT_TH = PRESENT(POUT_TH)
!
! IN variables
!
!$acc data present( PRHODJ, PEXNREF, PSIGS, PMFCONV, PPABST, PZZ,     &
!$acc &             PEXN, PCF_MF, PRC_MF, PRI_MF, PRV, PRC,           &
!$acc &             PTH, PRR, PRI, PRS, PRG, PRH,                     &
!
! INOUT variables
!
!$acc &             PRVS, PRCS, PRIS, PTHS,                           &
!
! OUT variables
!
!$acc &             PSRCS, PCLDFR, POUT_RV, POUT_RC, POUT_RI, POUT_TH )
!
!
!*       1.     PRELIMINARIES
!               -------------
!
!Check all IN arrays
CALL MPPDB_CHECK3D(PRHODJ,"ICE_ADJUST beg:PRHODJ",PRECISION)
CALL MPPDB_CHECK3D(PEXNREF,"ICE_ADJUST beg:PEXNREF",PRECISION)
CALL MPPDB_CHECK3D(PRHODREF,"ICE_ADJUST beg:PRHODREF")
CALL MPPDB_CHECK3D(PSIGS,"ICE_ADJUST beg:PSIGS",PRECISION)
CALL MPPDB_CHECK3D(PMFCONV,"ICE_ADJUST beg:PMFCONV",PRECISION)
CALL MPPDB_CHECK3D(PPABST,"ICE_ADJUST beg:PPABST",PRECISION)
CALL MPPDB_CHECK3D(PZZ,"ICE_ADJUST beg:PZZ",PRECISION)
CALL MPPDB_CHECK3D(PEXN,"ICE_ADJUST beg:PEXN",PRECISION)
CALL MPPDB_CHECK3D(PCF_MF,"ICE_ADJUST beg:PCF_MF",PRECISION)
CALL MPPDB_CHECK3D(PRC_MF,"ICE_ADJUST beg:PRC_MF",PRECISION)
CALL MPPDB_CHECK3D(PRI_MF,"ICE_ADJUST begPRI_MF:",PRECISION)
CALL MPPDB_CHECK3D(PRV,"ICE_ADJUST beg:PRV",PRECISION)
CALL MPPDB_CHECK3D(PRC,"ICE_ADJUST beg:PRC",PRECISION)
CALL MPPDB_CHECK3D(PTH,"ICE_ADJUST beg:PTH",PRECISION)
CALL MPPDB_CHECK3D(PRR,"ICE_ADJUST beg:PRR",PRECISION)
CALL MPPDB_CHECK3D(PRI,"ICE_ADJUST beg:PRI",PRECISION)
CALL MPPDB_CHECK3D(PRS,"ICE_ADJUST beg:PRS",PRECISION)
CALL MPPDB_CHECK3D(PRG,"ICE_ADJUST beg:PRG",PRECISION)
IF (PRESENT(PRH)) CALL MPPDB_CHECK3D(PRH,"ICE_ADJUST beg:PRH",PRECISION)
!Check all INOUT arrays
CALL MPPDB_CHECK3D(PRVS,"ICE_ADJUST beg:PRVS",PRECISION)
CALL MPPDB_CHECK3D(PRCS,"ICE_ADJUST beg:PRCS",PRECISION)
CALL MPPDB_CHECK3D(PTHS,"ICE_ADJUST beg:PTHS",PRECISION)
CALL MPPDB_CHECK3D(PRIS,"ICE_ADJUST beg:PRIS",PRECISION)

IIU = SIZE(PEXNREF,1)
IJU = SIZE(PEXNREF,2)
IKU = SIZE(PEXNREF,3)
CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
IKB=KKA+JPVEXT*KKL
IKE=KKU-JPVEXT*KKL

#ifndef MNH_OPENACC
allocate( gtemp  (IIU, IJU, IKU ) )
allocate( zsigs  (IIU, IJU, IKU ) )
allocate( zsrcs  (IIU, IJU, IKU ) )
allocate( zt     (IIU, IJU, IKU ) )
allocate( zrv    (IIU, IJU, IKU ) )
allocate( zrc    (IIU, IJU, IKU ) )
allocate( zri    (IIU, IJU, IKU ) )
allocate( zcph   (IIU, IJU, IKU ) )
allocate( zlv    (IIU, IJU, IKU ) )
allocate( zls    (IIU, IJU, IKU ) )
allocate( zw1    (IIU, IJU, IKU ) )
allocate( zw2    (IIU, IJU, IKU ) )
allocate( zcriaut(IIU, IJU, IKU ) )
allocate( zhcf   (IIU, IJU, IKU ) )
allocate( zhr    (IIU, IJU, IKU ) )
allocate( ZTEMP_BUD (IIU, IJU, IKU ) )
#else
!Pin positions in the pools of MNH memory
CALL MNH_MEM_POSITION_PIN()

CALL MNH_MEM_GET( gtemp  , IIU, IJU, IKU )
CALL MNH_MEM_GET( zsigs  , IIU, IJU, IKU )
CALL MNH_MEM_GET( zsrcs  , IIU, IJU, IKU )
CALL MNH_MEM_GET( zt     , IIU, IJU, IKU )
CALL MNH_MEM_GET( zrv    , IIU, IJU, IKU )
CALL MNH_MEM_GET( zrc    , IIU, IJU, IKU )
CALL MNH_MEM_GET( zri    , IIU, IJU, IKU )
CALL MNH_MEM_GET( zcph   , IIU, IJU, IKU )
CALL MNH_MEM_GET( zlv    , IIU, IJU, IKU )
CALL MNH_MEM_GET( zls    , IIU, IJU, IKU )
CALL MNH_MEM_GET( zw1    , IIU, IJU, IKU )
CALL MNH_MEM_GET( zw2    , IIU, IJU, IKU )
CALL MNH_MEM_GET( zcriaut, IIU, IJU, IKU )
CALL MNH_MEM_GET( zhcf   , IIU, IJU, IKU )
CALL MNH_MEM_GET( zhr    , IIU, IJU, IKU )
CALL MNH_MEM_GET( ZTEMP_BUD    , IIU, IJU, IKU )

!$acc data present( gtemp, zsigs, zsrcs, zt, zrv, zrc, zri, zcph, zlv, zls, zw1, zw2, zcriaut, zhcf, zhr )
#endif

if ( lbudget_th ) then
   !$acc kernels
   ZTEMP_BUD(:,:,:) = pths(:, :, :) * prhodj(:, :, :)
   !$acc end kernels
   call Budget_store_init( tbudgets(NBUDGET_TH), trim( hbuname ), ZTEMP_BUD(:,:,:) )
end if
if ( lbudget_rv ) then
   !$acc kernels
   ZTEMP_BUD(:,:,:) = prvs(:, :, :) * prhodj(:, :, :)
   !$acc end kernels
   call Budget_store_init( tbudgets(NBUDGET_RV), trim( hbuname ), ZTEMP_BUD(:,:,:) )
end if
if ( lbudget_rc ) then
   !$acc kernels
   ZTEMP_BUD(:,:,:) = prcs(:, :, :) * prhodj(:, :, :)
   !$acc end kernels
   call Budget_store_init( tbudgets(NBUDGET_RC), trim( hbuname ), ZTEMP_BUD(:,:,:) )
end if
if ( lbudget_ri ) then
   !$acc kernels
   ZTEMP_BUD(:,:,:) = pris(:, :, :) * prhodj(:, :, :)
   !$acc end kernels
   call Budget_store_init( tbudgets(NBUDGET_RI), trim( hbuname ), ZTEMP_BUD(:,:,:) )
end if
!
ITERMAX=1
!
!-------------------------------------------------------------------------------
!
!*       2.     COMPUTE QUANTITIES WITH THE GUESS OF THE FUTURE INSTANT
!               -------------------------------------------------------
!
!
!    beginning of the iterative loop (to compute the adjusted state)
!$acc kernels
ZRV(:,:,:)=PRV(:,:,:)
ZRC(:,:,:)=PRC(:,:,:)
ZRI(:,:,:)=PRI(:,:,:)
ZT(:,:,:)=PTH(:,:,:) * PEXN(:,:,:)
!$acc end kernels
!
DO JITER =1,ITERMAX
  !
  !*       2.3    compute the latent heat of vaporization Lv(T*) at t+1
  !                   and the latent heat of sublimation  Ls(T*) at t+1
  !
  ZLV(:,:,:) = XLVTT + ( XCPV - XCL ) * ( ZT(:,:,:) -XTT )
  ZLS(:,:,:) = XLSTT + ( XCPV - XCI ) * ( ZT(:,:,:) -XTT )
  !
  !*       2.4    compute the specific heat for moist air (Cph) at t+1
  !
  IF     ( KRR == 7 ) THEN
    ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:)                             &
                       + XCL  * (ZRC(:,:,:) + PRR(:,:,:))             &
                       + XCI  * (ZRI(:,:,:) + PRS(:,:,:) + PRG(:,:,:) + PRH(:,:,:))
  ELSE IF( KRR == 6 ) THEN
    ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:)                             &
                       + XCL  * (ZRC(:,:,:) + PRR(:,:,:))             &
                       + XCI  * (ZRI(:,:,:) + PRS(:,:,:) + PRG(:,:,:))
  ELSE IF( KRR == 5 ) THEN
    ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:)                             &
                       + XCL  * (ZRC(:,:,:) + PRR(:,:,:))             &
                       + XCI  * (ZRI(:,:,:) + PRS(:,:,:))
  ELSE IF( KRR == 3 ) THEN
    ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:)               &
                       + XCL  * (ZRC(:,:,:) + PRR(:,:,:))
  ELSE IF( KRR == 2 ) THEN
    ZCPH(:,:,:) = XCPD + XCPV * ZRV(:,:,:) &
                       + XCL  * ZRC(:,:,:)
  END IF
!$acc end kernels
  !
  IF ( OSUBG_COND ) THEN
    !
    !*       3.     SUBGRID CONDENSATION SCHEME
    !               ---------------------------
    !
    !   PSRC= s'rci'/Sigma_s^2
    !   ZT, ZRV, ZRC and ZRI are INOUT
    CALL CONDENSATION(IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, KKL, &
         HFRAC_ICE, HCONDENS, HLAMBDA3,                                    &
         PPABST, PZZ, PRHODREF, ZT, ZRV, ZRC, ZRI, PRS, PRG, PSIGS, PMFCONV, PCLDFR, &
         PSRCS, .TRUE., OSIGMAS,                                           &
         PSIGQSAT, PLV=ZLV, PLS=ZLS, PCPH=ZCPH, PHLC_HRC=PHLC_HRC, PHLC_HCF=PHLC_HCF, PHLI_HRI=PHLI_HRI, PHLI_HCF=PHLI_HCF)
  ELSE
    !
    !*       4.     ALL OR NOTHING CONDENSATION SCHEME
    !                            FOR MIXED-PHASE CLOUD
    !               -----------------------------------------------
    !
    !
    !   ZT, ZRV, ZRC and ZRI are INOUT
    !
    !CALL ADJUST_LANGLOIS(IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, KKL,                &
    !                     PPABST, ZT, ZRV, ZRC, ZRI, ZLV, ZLS, ZCPH) HFRAC_ICE must be implemented in Langlois before using it again
!PW: TODO: ZSRCS not used !!! (OUT variable from CONDENSATION)
!PW:TODO: autres optis: ZSIGS=0 and PSIGQSAT=0 -> less computation...
!$acc kernels
    ZSIGS=0.
!$acc end kernels
    CALL CONDENSATION(IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, KKL,                                            &
         HFRAC_ICE, HCONDENS, HLAMBDA3,                                                                            &
         PPABST, PZZ, PRHODREF, ZT, ZRV, ZRC, ZRI, PRS, PRG, ZSIGS, PMFCONV, PCLDFR,                               &
         ZSRCS, .TRUE., OSIGMAS=.TRUE.,                                                                            &
         PSIGQSAT=0., PLV=ZLV, PLS=ZLS, PHLC_HRC=PHLC_HRC, PHLC_HCF=PHLC_HCF, PHLI_HRI=PHLI_HRI, PHLI_HCF=PHLI_HCF )
  END IF
ENDDO         ! end of the iterative loop
!
!*       5.     COMPUTE THE SOURCES AND STORES THE CLOUD FRACTION
!               -------------------------------------------------
!
!$acc kernels
#ifdef MNH_COMPILER_NVHPC
!$acc loop independent collapse(3)
#endif
DO CONCURRENT (JI=1:IIU,JJ=1:IJU,JK=1:IKU)

!
!*       5.0    compute the variation of mixing ratio
!
                                                      !         Rc - Rc*
ZW1(JI,JJ,JK) = (ZRC(JI,JJ,JK) - PRC(JI,JJ,JK)) / PTSTEP       ! Pcon = ----------
                                                      !         2 Delta t

ZW2(JI,JJ,JK) = (ZRI(JI,JJ,JK) - PRI(JI,JJ,JK)) / PTSTEP       ! idem ZW1 but for Ri
!
!*       5.1    compute the sources
!

GTEMP(JI,JJ,JK) = ZW1(JI,JJ,JK) < 0.0
IF( GTEMP(JI,JJ,JK) )THEN
  ZW1(JI,JJ,JK) = MAX ( ZW1(JI,JJ,JK), -PRCS(JI,JJ,JK) )
ELSE
  ZW1(JI,JJ,JK) = MIN ( ZW1(JI,JJ,JK),  PRVS(JI,JJ,JK) )
ENDIF
PRVS(JI,JJ,JK) = PRVS(JI,JJ,JK) - ZW1(JI,JJ,JK)
PRCS(JI,JJ,JK) = PRCS(JI,JJ,JK) + ZW1(JI,JJ,JK)
PTHS(JI,JJ,JK) = PTHS(JI,JJ,JK) +        &
                ZW1(JI,JJ,JK) * ZLV(JI,JJ,JK) / (ZCPH(JI,JJ,JK) * PEXNREF(JI,JJ,JK))
!
GTEMP(JI,JJ,JK) = ZW2(JI,JJ,JK) < 0.0
IF( GTEMP(JI,JJ,JK) )THEN
  ZW2(JI,JJ,JK) = MAX ( ZW2(JI,JJ,JK), -PRIS(JI,JJ,JK) )
ELSE
  ZW2(JI,JJ,JK) = MIN ( ZW2(JI,JJ,JK),  PRVS(JI,JJ,JK) )
ENDIF
PRVS(JI,JJ,JK) = PRVS(JI,JJ,JK) - ZW2(JI,JJ,JK)
PRIS(JI,JJ,JK) = PRIS(JI,JJ,JK) + ZW2(JI,JJ,JK)
PTHS(JI,JJ,JK) = PTHS(JI,JJ,JK) +        &
     ZW2(JI,JJ,JK) * ZLS(JI,JJ,JK) / (ZCPH(JI,JJ,JK) * PEXNREF(JI,JJ,JK))
ENDDO

!$acc end kernels
!
!
!*       5.2    compute the cloud fraction PCLDFR
!
IF ( .NOT. OSUBG_COND ) THEN
!$acc kernels present_cr(GTEMP)
DO CONCURRENT (JI=1:IIU,JJ=1:IJU,JK=1:IKU)   
  GTEMP(JI,JJ,JK) = PRCS(JI,JJ,JK) + PRIS(JI,JJ,JK) > 1.E-12 / PTSTEP
  IF ( GTEMP(JI,JJ,JK) )THEN
    PCLDFR(JI,JJ,JK)  = 1.
  ELSE
    PCLDFR(JI,JJ,JK)  = 0.
 ENDIF
ENDDO
  IF ( SIZE(PSRCS,3) /= 0 ) THEN
    PSRCS(:,:,:) = PCLDFR(:,:,:)
  END IF
!$acc end kernels
ELSE
#ifdef MNH_OPENACC
  call Print_msg( NVERB_FATAL, 'GEN', 'ICE_ADJUST', 'OpenACC: OSUBG_COND=.T. not yet implemented' )
#endif
  !We limit PRC_MF+PRI_MF to PRVS*PTSTEP to avoid negative humidity
  ZW1(:,:,:)=PRC_MF(:,:,:)/PTSTEP
  ZW2(:,:,:)=PRI_MF(:,:,:)/PTSTEP
  GTEMP(:,:,:) = ZW1(:,:,:)+ZW2(:,:,:)>PRVS(:,:,:)
  WHERE(GTEMP(:,:,:))
    ZW1(:,:,:)=ZW1(:,:,:)*PRVS(:,:,:)/(ZW1(:,:,:)+ZW2(:,:,:))
    ZW2(:,:,:)=PRVS(:,:,:)-ZW1(:,:,:)
  ENDWHERE
  IF(PRESENT(PHLC_HRC) .AND. PRESENT(PHLC_HCF)) THEN
    ZCRIAUT(:,:,:)=XCRIAUTC/PRHODREF
    IF(HSUBG_MF_PDF=='NONE')THEN
      WHERE(ZW1(:,:,:)*PTSTEP>PCF_MF * ZCRIAUT)
        PHLC_HRC(:,:,:)=PHLC_HRC(:,:,:)+ZW1(:,:,:)*PTSTEP
        PHLC_HCF(:,:,:)=MIN(1.,PHLC_HCF(:,:,:)+PCF_MF(:,:,:))
      ENDWHERE
    ELSEIF(HSUBG_MF_PDF=='TRIANGLE')THEN
      !ZHCF is the precipitating part of the *cloud* and not of the grid cell
      WHERE(ZW1(:,:,:)*PTSTEP>PCF_MF*ZCRIAUT(:,:,:))
        ZHCF(:,:,:)=1.-.5*(ZCRIAUT(:,:,:)*PCF_MF(:,:,:) / MAX(1.E-20, ZW1(:,:,:)*PTSTEP))**2
        ZHR(:,:,:)=ZW1(:,:,:)*PTSTEP-(ZCRIAUT(:,:,:)*PCF_MF(:,:,:))**3 / &
                                    &(3*MAX(1.E-20, ZW1(:,:,:)*PTSTEP)**2)
      ELSEWHERE(2.*ZW1(:,:,:)*PTSTEP<=PCF_MF * ZCRIAUT(:,:,:))
        ZHCF(:,:,:)=0.
        ZHR(:,:,:)=0.
      ELSEWHERE
        ZHCF(:,:,:)=(2.*ZW1(:,:,:)*PTSTEP-ZCRIAUT(:,:,:)*PCF_MF(:,:,:))**2 / &
                   &(2.*MAX(1.E-20, ZW1(:,:,:)*PTSTEP)**2)
        ZHR(:,:,:)=(4.*(ZW1(:,:,:)*PTSTEP)**3-3.*ZW1(:,:,:)*PTSTEP*(ZCRIAUT(:,:,:)*PCF_MF(:,:,:))**2+&
                    (ZCRIAUT(:,:,:)*PCF_MF(:,:,:))**3) / &
                  &(3*MAX(1.E-20, ZW1(:,:,:)*PTSTEP)**2)
      ENDWHERE
      ZHCF(:,:,:)=ZHCF(:,:,:)*PCF_MF(:,:,:) !to retrieve the part of the grid cell
      PHLC_HCF(:,:,:)=MIN(1.,PHLC_HCF(:,:,:)+ZHCF(:,:,:)) !total part of the grid cell that is precipitating
      PHLC_HRC(:,:,:)=PHLC_HRC(:,:,:)+ZHR(:,:,:)
    ENDIF
  ENDIF
  IF(PRESENT(PHLI_HRI) .AND. PRESENT(PHLI_HCF)) THEN
    ZCRIAUT(:,:,:)=MIN(XCRIAUTI,10**(XACRIAUTI*(ZT(:,:,:)-XTT)+XBCRIAUTI))
    IF(HSUBG_MF_PDF=='NONE')THEN
      WHERE(ZW2(:,:,:)*PTSTEP>PCF_MF * ZCRIAUT(:,:,:))
        PHLI_HRI(:,:,:)=PHLI_HRI(:,:,:)+ZW2(:,:,:)*PTSTEP
        PHLI_HCF(:,:,:)=MIN(1.,PHLI_HCF(:,:,:)+PCF_MF(:,:,:))
      ENDWHERE
    ELSEIF(HSUBG_MF_PDF=='TRIANGLE')THEN
      !ZHCF is the precipitating part of the *cloud* and not of the grid cell
      WHERE(ZW2(:,:,:)*PTSTEP>PCF_MF*ZCRIAUT)
        ZHCF(:,:,:)=1.-.5*(ZCRIAUT*PCF_MF(:,:,:) / (ZW2(:,:,:)*PTSTEP))**2
        ZHR(:,:,:)=ZW2(:,:,:)*PTSTEP-(ZCRIAUT*PCF_MF(:,:,:))**3/(3*(ZW2(:,:,:)*PTSTEP)**2)
      ELSEWHERE(2.*ZW2(:,:,:)*PTSTEP<=PCF_MF * ZCRIAUT)
        ZHCF(:,:,:)=0.
        ZHR(:,:,:)=0.
      ELSEWHERE
        ZHCF(:,:,:)=(2.*ZW2(:,:,:)*PTSTEP-ZCRIAUT*PCF_MF(:,:,:))**2 / (2.*(ZW2(:,:,:)*PTSTEP)**2)
        ZHR(:,:,:)=(4.*(ZW2(:,:,:)*PTSTEP)**3-3.*ZW2(:,:,:)*PTSTEP*(ZCRIAUT*PCF_MF(:,:,:))**2+&
                    (ZCRIAUT*PCF_MF(:,:,:))**3)/(3*(ZW2(:,:,:)*PTSTEP)**2)
      ENDWHERE
      ZHCF(:,:,:)=ZHCF(:,:,:)*PCF_MF(:,:,:) !to retrieve the part of the grid cell
      PHLI_HCF(:,:,:)=MIN(1.,PHLI_HCF(:,:,:)+ZHCF(:,:,:)) !total part of the grid cell that is precipitating
      PHLI_HRI(:,:,:)=PHLI_HRI(:,:,:)+ZHR(:,:,:)
    ENDIF
  ENDIF
  PCLDFR(:,:,:)=MIN(1.,PCLDFR(:,:,:)+PCF_MF(:,:,:))
  PRCS(:,:,:)=PRCS(:,:,:)+ZW1(:,:,:)
  PRIS(:,:,:)=PRIS(:,:,:)+ZW2(:,:,:)
  PRVS(:,:,:)=PRVS(:,:,:)-(ZW1(:,:,:)+ZW2(:,:,:))
  PTHS(:,:,:) = PTHS(:,:,:) + &
                (ZW1 * ZLV(:,:,:) + ZW2 * ZLS(:,:,:)) / ZCPH(:,:,:)     &
                /  PEXNREF(:,:,:)
  IF(GPOUT_RV .OR. GPOUT_RC .OR. &
    &GPOUT_RI .OR. GPOUT_TH) THEN
    ZW1(:,:,:)=PRC_MF(:,:,:)
    ZW2(:,:,:)=PRI_MF(:,:,:)
    GTEMP(:,:,:) = ZW1(:,:,:)+ZW2(:,:,:)>ZRV(:,:,:)
    WHERE(GTEMP(:,:,:))
      ZW1(:,:,:)=ZW1(:,:,:)*ZRV(:,:,:)/(ZW1(:,:,:)+ZW2(:,:,:))
      ZW2(:,:,:)=ZRV(:,:,:)-ZW1(:,:,:)
    ENDWHERE
    ZRC(:,:,:)=ZRC(:,:,:)+ZW1(:,:,:)
    ZRI(:,:,:)=ZRI(:,:,:)+ZW2(:,:,:)
    ZRV(:,:,:)=ZRV(:,:,:)-(ZW1(:,:,:)+ZW2(:,:,:))
    ZT(:,:,:) = ZT(:,:,:) + &
                (ZW1 * ZLV(:,:,:) + ZW2 * ZLS(:,:,:)) / ZCPH(:,:,:)
  ENDIF
ENDIF
!
!$acc kernels
IF(GPOUT_RV) POUT_RV=ZRV
IF(GPOUT_RC) POUT_RC=ZRC
IF(GPOUT_RI) POUT_RI=ZRI
IF(GPOUT_TH) POUT_TH=ZT / PEXN(:,:,:)
!$acc end kernels
!
!
!*       6.  STORE THE BUDGET TERMS
!            ----------------------
!
if ( lbudget_th ) then
   !$acc kernels
   ZTEMP_BUD(:,:,:) = pths(:, :, :) * prhodj(:, :, :)
   !$acc end kernels
   call Budget_store_end( tbudgets(NBUDGET_TH), trim( hbuname ), ZTEMP_BUD(:,:,:) )
   !$acc kernels
   ZTEMP_BUD(:,:,:) = prvs(:, :, :) * prhodj(:, :, :)
   !$acc end kernels
   call Budget_store_end( tbudgets(NBUDGET_RV), trim( hbuname ), ZTEMP_BUD(:,:,:) )
   !$acc kernels
   ZTEMP_BUD(:,:,:) = prcs(:, :, :) * prhodj(:, :, :)
   !$acc end kernels
   call Budget_store_end( tbudgets(NBUDGET_RC), trim( hbuname ), ZTEMP_BUD(:,:,:) )
   !$acc kernels
   ZTEMP_BUD(:,:,:) = pris(:, :, :) * prhodj(:, :, :)
   !$acc end kernels
   call Budget_store_end( tbudgets(NBUDGET_RI), trim( hbuname ), ZTEMP_BUD(:,:,:) )
deallocate( gtemp, zsigs, zsrcs, zt, zrv, zrc, zri, zcph, zlv, zls, zw1, zw2, zcriaut, zhcf, zhr , ZTEMP_BUD )
#else
!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN
CALL MNH_MEM_RELEASE()
#endif
!------------------------------------------------------------------------------
!
!Check all INOUT arrays
CALL MPPDB_CHECK3D(PRVS,"ICE_ADJUST end:PRVS",PRECISION)
CALL MPPDB_CHECK3D(PRCS,"ICE_ADJUST end:PRCS",PRECISION)
CALL MPPDB_CHECK3D(PTHS,"ICE_ADJUST end:PTHS",PRECISION)
CALL MPPDB_CHECK3D(PRIS,"ICE_ADJUST end:PRIS",PRECISION)
!Check all OUT arrays
CALL MPPDB_CHECK3D(PSRCS,"ICE_ADJUST end:PSRCS",PRECISION)
CALL MPPDB_CHECK3D(PCLDFR,"ICE_ADJUST end:PCLDFR",PRECISION)
IF (GPOUT_RV) CALL MPPDB_CHECK3D(POUT_RV,"ICE_ADJUST end:POUT_RV")
IF (GPOUT_RC) CALL MPPDB_CHECK3D(POUT_RC,"ICE_ADJUST end:POUT_RC")
IF (GPOUT_RI) CALL MPPDB_CHECK3D(POUT_RI,"ICE_ADJUST end:POUT_RI")
IF (GPOUT_TH) CALL MPPDB_CHECK3D(POUT_TH,"ICE_ADJUST end:POUT_TH")
IF (PRESENT(PHLC_HRC)) CALL MPPDB_CHECK3D(PHLC_HRC,"ICE_ADJUST end:PHLC_HRC")
IF (PRESENT(PHLC_HCF)) CALL MPPDB_CHECK3D(PHLC_HCF,"ICE_ADJUST end:PHLC_HCF")
IF (PRESENT(PHLI_HRI)) CALL MPPDB_CHECK3D(PHLI_HRI,"ICE_ADJUST end:PHLI_HRI")
IF (PRESENT(PHLI_HCF)) CALL MPPDB_CHECK3D(PHLI_HCF,"ICE_ADJUST end:PHLI_HCF")

!$acc end data

END SUBROUTINE ICE_ADJUST