turb.f90

   CALL GZ_M_W_PHY(D,ZMWTH,PDZZ,ZWORK1)    ! -d(w'2th' )/dz
   CALL GZ_W_M_PHY(D,ZMTH2,PDZZ,ZWORK2)    ! -d(w'th'2 )/dz
   !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
   ZFWTH(:,:) = -ZWORK1(:,:)
   ZFTH2(:,:) = -ZWORK2(:,:)
   !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!
  ZFWTH(:,IKTE:) = 0.
  ZFWTH(:,:IKTB) = 0.
  ZFWR(:,:)  = 0.
  ZFTH2(:,IKTE:) = 0.
  ZFTH2(:,:IKTB) = 0.
  ZFR2(:,:)  = 0.
  ZFTHR(:,:) = 0.
ELSE
!$acc kernels present_cr(ZFWTH,ZFWR,ZFTH2,ZFR2,ZFTHR)
  ZFWTH(:,:) = 0.
  ZFWR(:,:)  = 0.
  ZFTH2(:,:) = 0.
  ZFR2(:,:)  = 0.
  ZFTHR(:,:) = 0.
!$acc end kernels
ENDIF
!
!----------------------------------------------------------------------------
!
!*      5. TURBULENT SOURCES
!          -----------------
!
IF( BUCONF%LBUDGET_U )  CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_U ), 'VTURB', PRUS(:,:)    )
IF( BUCONF%LBUDGET_V )  CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_V ), 'VTURB', PRVS(:,:)    )
IF( BUCONF%LBUDGET_W )  CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_W ), 'VTURB', PRWS(:,:)    )

IF( BUCONF%LBUDGET_TH ) THEN
  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLVOCPEXNM(:,:) * PRRS(:,:, 2) + ZLSOCPEXNM(:,:) * PRRS(:,:, 4) 
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', ZTEMP_BUD(:,:) )
  ELSE IF( KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD, ZLOCPEXNM)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:, 2)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', ZTEMP_BUD(:,:) )
  ELSE
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:,:) )
  END IF
END IF

IF( BUCONF%LBUDGET_RV ) THEN
  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRRS(:,:, 1) - PRRS(:,:, 2) - PRRS(:,:, 4) 
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', ZTEMP_BUD(:,:) )
  ELSE IF( KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRRS(:,:, 1) - PRRS(:,:, 2)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', ZTEMP_BUD(:,:) )
  ELSE
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:,:, 1) )
  END IF
END IF

IF( BUCONF%LBUDGET_RC ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RC), 'VTURB', PRRS  (:,:, 2) )
IF( BUCONF%LBUDGET_RI ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RI), 'VTURB', PRRS  (:,:, 4) )

IF( BUCONF%LBUDGET_SV ) THEN
  DO JSV = 1, KSV
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'VTURB', PRSVS(:,:, JSV) )
  END DO
END IF
!$acc update device(PRHODJ)
!$acc update device(PRUS,PRVS,PRWS,PRSVS)
CALL TURB_VER(D,CST,CSTURB,TURBN,TLES,                   &
          KRR,KRRL,KRRI,KGRADIENTS,                      &
          GOCEAN, ODEEPOC, OCOMPUTE_SRC,                 &
          KSV,KSV_LGBEG,KSV_LGEND,                       &
          ZEXPL,HPROGRAM, O2D, ONOMIXLG, OFLAT,          &
          OCOUPLES,OBLOWSNOW,OFLYER, PRSNOW,             &
          PTSTEP,TPFILE,                                 &
          PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ,        &
          PCOSSLOPE,PSINSLOPE,                           &
          PRHODJ,PTHVREF,PSFU,PSFV,                      &
          PSFTH,PSFRV,PSFSV,PSFTH,PSFRV,PSFSV,           &
          ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU33M,               &
          PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT,    &
          PTKET,ZLM,PLENGTHM,PLENGTHH,ZLEPS,MFMOIST,     &
          ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE,     &
          ZFWTH,ZFWR,ZFTH2,ZFR2,ZFTHR,PBL_DEPTH,         &
          PSBL_DEPTH,ZLMO,PHGRAD,PZS,                    &
          PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,              &
          PDP,PTP,PSIGS,PWTH,PWRC,PWSV,                  &
          PSSTFL, PSSTFL_C, PSSRFL_C,PSSUFL_C,PSSVFL_C,  &
          PSSUFL,PSSVFL                                  )
!$acc update self(PWTH,PWRC,PWSV)
!IF (HCLOUD == 'LIMA') THEN
!   IF (KSV_LIMA_NR.GT.0) PRSVS(:,:,KSV_LIMA_NR) = ZRSVS(:,:,KSV_LIMA_NR) 
!   IF (KSV_LIMA_NS.GT.0) PRSVS(:,:,KSV_LIMA_NS) = ZRSVS(:,:,KSV_LIMA_NS)
!   IF (KSV_LIMA_NG.GT.0) PRSVS(:,:,KSV_LIMA_NG) = ZRSVS(:,:,KSV_LIMA_NG) 
!   IF (KSV_LIMA_NH.GT.0) PRSVS(:,:,KSV_LIMA_NH) = ZRSVS(:,:,KSV_LIMA_NH)
!END IF

IF( BUCONF%LBUDGET_U ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_U), 'VTURB', PRUS(:,:) )
IF( BUCONF%LBUDGET_V ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_V), 'VTURB', PRVS(:,:) )
IF( BUCONF%LBUDGET_W ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_W), 'VTURB', PRWS(:,:) )

IF( BUCONF%LBUDGET_TH ) THEN
  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLVOCPEXNM(:,:) * PRRS(:,:, 2) + ZLSOCPEXNM(:,:) * PRRS(:,:, 4) 
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', ZTEMP_BUD(:,:) )
  ELSE IF( KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD, ZLOCPEXNM)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:, 2)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', ZTEMP_BUD(:,:) )
  ELSE
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'VTURB', PRTHLS(:,:) )
  END IF
END IF

IF( BUCONF%LBUDGET_RV ) THEN
  IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRRS(:,:, 1) - PRRS(:,:, 2) - PRRS(:,:, 4) 
     !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', ZTEMP_BUD(:,:) )
  ELSE IF( KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRRS(:,:, 1) - PRRS(:,:, 2) 
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', ZTEMP_BUD(:,:))
  ELSE
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'VTURB', PRRS(:,:, 1) )
  END IF
END IF

IF( BUCONF%LBUDGET_RC ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RC), 'VTURB', PRRS(:,:, 2) )
IF( BUCONF%LBUDGET_RI ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RI), 'VTURB', PRRS(:,:, 4) )

IF( BUCONF%LBUDGET_SV )  THEN
  DO JSV = 1, KSV
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'VTURB', PRSVS(:,:, JSV) )
  END DO
END IF
!
!Les budgets des termes horizontaux de la turb sont présents dans AROME
! alors que ces termes ne sont pas calculés
#ifdef REPRO48
#else
IF( TURBN%CTURBDIM == '3DIM' ) THEN
#endif
  IF( BUCONF%LBUDGET_U  ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_U ), 'HTURB', PRUS  (:,:) )
  IF( BUCONF%LBUDGET_V  ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_V ), 'HTURB', PRVS  (:,:) )
  IF( BUCONF%LBUDGET_W  ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_W ), 'HTURB', PRWS  (:,:) )

  IF(BUCONF%LBUDGET_TH)  THEN
    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLVOCPEXNM(:,:) * PRRS(:,:, 2) + ZLSOCPEXNM(:,:) * PRRS(:,:, 4)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', ZTEMP_BUD(:,:) )
    ELSE IF( KRRL >= 1 ) THEN
      !$acc kernels present_cr(ZTEMP_BUD)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:, 2)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', ZTEMP_BUD(:,:)  )
    ELSE
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:,:) )
    END IF
  END IF

  IF( BUCONF%LBUDGET_RV ) THEN
    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
      !$acc kernels present_cr(ZTEMP_BUD)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZTEMP_BUD(:,:) =  PRRS(:,:, 1) - PRRS(:,:, 2) - PRRS(:,:, 4)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', ZTEMP_BUD(:,:) )
    ELSE IF( KRRL >= 1 ) THEN
      !$acc kernels present_cr(ZTEMP_BUD)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZTEMP_BUD(:,:) =  PRRS(:,:, 1) - PRRS(:,:, 2)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', ZTEMP_BUD(:,:) )
    ELSE
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:,:, 1) )
    END IF
  END IF

  IF( BUCONF%LBUDGET_RC ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RC), 'HTURB', PRRS(:,:, 2) )
  IF( BUCONF%LBUDGET_RI ) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_RI), 'HTURB', PRRS(:,:, 4) )

  IF( BUCONF%LBUDGET_SV )  THEN
    DO JSV = 1, KSV
      CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:,:, JSV) )
    END DO
  END IF
!à supprimer une fois le précédent ifdef REPRO48 validé
#ifdef REPRO48
#else
    CALL TURB_HOR_SPLT(D,CST,CSTURB, TURBN, TLES,              &
          KSPLIT, KRR, KRRL, KRRI, KSV,KSV_LGBEG,KSV_LGEND,    & 
          PTSTEP,HLBCX,HLBCY, OFLAT,O2D, ONOMIXLG,             & 
          GOCEAN,OCOMPUTE_SRC,OBLOWSNOW,PRSNOW,                &
          TPFILE, HPROGRAM, KHALO,                             &
          PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ,                        &
          PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,                       &
          PCOSSLOPE,PSINSLOPE,                                 &
          PRHODJ,PTHVREF,                                      &
          PSFTH,PSFRV,PSFSV,                                   &
          ZCDUEFF,ZTAU11M,ZTAU12M,ZTAU22M,ZTAU33M,             &
          PUT,PVT,PWT,ZUSLOPE,ZVSLOPE,PTHLT,PRT,PSVT,          &
          PTKET,ZLM,ZLEPS,                                     &
          ZLOCPEXNM,ZATHETA,ZAMOIST,PSRCT,ZFRAC_ICE,           &
          PDP,PTP,PSIGS,                                       &
          ZTRH,                                                &
          PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS                     )
#endif
  !
!  IF (HCLOUD == 'LIMA') THEN
!     IF (KSV_LIMA_NR.GT.0) PRSVS(:,:,KSV_LIMA_NR) = ZRSVS(:,:,KSV_LIMA_NR) 
!     IF (KSV_LIMA_NS.GT.0) PRSVS(:,:,KSV_LIMA_NS) = ZRSVS(:,:,KSV_LIMA_NS)
!     IF (KSV_LIMA_NG.GT.0) PRSVS(:,:,KSV_LIMA_NG) = ZRSVS(:,:,KSV_LIMA_NG) 
!     IF (KSV_LIMA_NH.GT.0) PRSVS(:,:,KSV_LIMA_NH) = ZRSVS(:,:,KSV_LIMA_NH)
!  END IF
  !
  IF( BUCONF%LBUDGET_U ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_U), 'HTURB', PRUS(:,:) )
  IF( BUCONF%LBUDGET_V ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_V), 'HTURB', PRVS(:,:) )
  IF( BUCONF%LBUDGET_W ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_W), 'HTURB', PRWS(:,:) )

  IF( BUCONF%LBUDGET_TH ) THEN
    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
      !$acc kernels present_cr(ZTEMP_BUD)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLVOCPEXNM(:,:) * PRRS(:,:, 2) + ZLSOCPEXNM(:,:) * PRRS(:,:, 4)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', ZTEMP_BUD(:,:) )
    ELSE IF( KRRL >= 1 ) THEN
      !$acc kernels present_cr(ZTEMP_BUD)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:, 2)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', ZTEMP_BUD(:,:) )
    ELSE
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'HTURB', PRTHLS(:,:) )
    END IF
  END IF

  IF( BUCONF%LBUDGET_RV ) THEN
    IF( KRRI >= 1 .AND. KRRL >= 1 ) THEN
      !$acc kernels present_cr(ZTEMP_BUD)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZTEMP_BUD(:,:) =  PRRS(:,:, 1) - PRRS(:,:, 2) - PRRS(:,:, 4)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', ZTEMP_BUD(:,:) )
    ELSE IF( KRRL >= 1 ) THEN
      !$acc kernels present_cr(ZTEMP_BUD)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZTEMP_BUD(:,:) =  PRRS(:,:, 1) - PRRS(:,:, 2)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', ZTEMP_BUD(:,:) )
    ELSE
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RV), 'HTURB', PRRS(:,:, 1) )
    END IF
  END IF

  IF( BUCONF%LBUDGET_RC ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RC), 'HTURB', PRRS(:,:, 2) )
  IF( BUCONF%LBUDGET_RI ) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_RI), 'HTURB', PRRS(:,:, 4) )

  IF( BUCONF%LBUDGET_SV )  THEN
    DO JSV = 1, KSV
      CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:,:, JSV) )
    END DO
  END IF
#ifdef REPRO48
#else
END IF
#endif
!$acc update self(PSIGS,PRUS,PRVS,PRWS,PRSVS)
!----------------------------------------------------------------------------
!
!*      6. EVOLUTION OF THE TKE AND ITS DISSIPATION
!          ----------------------------------------
!
!  6.1 Contribution of mass-flux in the TKE buoyancy production if
!      cloud computation is not statistical
CALL MZF_PHY(D,PFLXZTHVMF,ZWORK1)
!$acc kernels present_cr(PTHP)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
PTP(:,:) = PTP(:,:) &
                             + CST%XG / PTHVREF(:,:) * ZWORK1(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
!
IF(PRESENT(PTPMF))  THEN
  !$acc kernels
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  PTPMF(:,:)=CST%XG / PTHVREF(:,:) * ZWORK1(:,:)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  !$acc end kernels
END IF
!  6.2 TKE evolution equation

IF (.NOT. TURBN%LHARAT) THEN
!
IF (BUCONF%LBUDGET_TH)  THEN
  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:)+ ZLVOCPEXNM(:,:) * PRRS(:,:,2) + ZLSOCPEXNM(:,:) * PRRS(:,:,4)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', ZTEMP_BUD(:,:) )
  ELSE IF ( KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:,2)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', ZTEMP_BUD(:,:) )
  ELSE
    CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:,:) )
  END IF
END IF
!
IF(PRESENT(PRTKEMS)) THEN
!$acc kernels
  ZRTKEMS(:,:)=PRTKEMS(:,:)
!$acc end kernels
ELSE
!$acc kernels
  ZRTKEMS(:,:)=0.
!$acc end kernels
END IF
!
!$acc update device(PRTKES)
CALL TKE_EPS_SOURCES(D,CST,CSTURB,BUCONF,TURBN,TLES,HPROGRAM,           &
                   & KMI,PTKET,ZLM,ZLEPS,PDP,ZTRH,                      &
                   & PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ,               &
                   & PTSTEP,ZEXPL,                                      &
                   & TPFILE,ODIAG_IN_RUN,GOCEAN,                        &
                   & PSFU,PSFV,                                         &
                   & PTP,PRTKES,PRTHLS,ZCOEF_DISS,PTDIFF,PTDISS,ZRTKEMS,&
                   & TBUDGETS,KBUDGETS, PEDR=PEDR, PTR=PTR,PDISS=PDISS, &
                   & PCURRENT_TKE_DISS=PCURRENT_TKE_DISS                )
                   !
!$acc update self(PTR,PDISS)
!
!$acc update self(PRTKES)
IF (BUCONF%LBUDGET_TH)  THEN
  IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:)+ ZLVOCPEXNM(:,:) * PRRS(:,:,2) + ZLSOCPEXNM(:,:) * PRRS(:,:,4)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', ZTEMP_BUD(:,:) )
  ELSE IF ( KRRL >= 1 ) THEN
    !$acc kernels present_cr(ZTEMP_BUD)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZTEMP_BUD(:,:) =  PRTHLS(:,:) + ZLOCPEXNM(:,:) * PRRS(:,:,2)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', ZTEMP_BUD(:,:) )
  ELSE
    CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TH), 'DISSH', PRTHLS(:,:) )
  END IF
END IF
!
ENDIF
!
!----------------------------------------------------------------------------
!
!*      7. STORES SOME INFORMATIONS RELATED TO THE TURBULENCE SCHEME
!          ---------------------------------------------------------
!
!$acc update self(PLEM)
IF ( TURBN%LTURB_DIAG .AND. TPFILE%LOPENED ) THEN
!
! stores the mixing length
!
  TZFIELD = TFIELDMETADATA(       &
    CMNHNAME   = 'LM',            &
    CSTDNAME   = '',              &
    CLONGNAME  = 'LM',            &
    CUNITS     = 'm',             &
    CDIR       = 'XY',            &
    CCOMMENT   = 'Mixing length', &
    NGRID      = 1,               &
    NTYPE      = TYPEREAL,        &
    NDIMS      = 3,               &
    LTIMEDEP   = .TRUE.           )
  CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZLM)
!
  IF (KRR /= 0) THEN
!
! stores the conservative potential temperature
!
    TZFIELD = TFIELDMETADATA(                          &
    CMNHNAME   = 'THLM',                               &
    CSTDNAME   = '',                                   &
    CLONGNAME  = 'THLM',                               &
    CUNITS     = 'K',                                  &
    CDIR       = 'XY',                                 &
    CCOMMENT   = 'Conservative potential temperature', &
    NGRID      = 1,                                    &
    NTYPE      = TYPEREAL,                             &
    NDIMS      = 3,                                    &
    LTIMEDEP   = .TRUE.                                )
!$acc update self(PTHLT)
    CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,PTHLT)
!
! stores the conservative mixing ratio
!
    TZFIELD = TFIELDMETADATA(                &
    CMNHNAME   = 'RNPM',                     &
    CSTDNAME   = '',                         &
    CLONGNAME  = 'RNPM',                     &
    CUNITS     = 'kg kg-1',                  &
    CDIR       = 'XY',                       &
    CCOMMENT   = 'Conservative mixing ratio',&
    NGRID      = 1,                          &
    NTYPE      = TYPEREAL,                   &
    NDIMS      = 3,                          &
    LTIMEDEP   = .TRUE.                      )
!$acc update self(PRT)
    CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,PRT(:,:,1))
   END IF
END IF
!
!* stores value of conservative variables & wind before turbulence tendency (AROME only)
IF(PRESENT(PDRUS_TURB)) THEN
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  PDRUS_TURB(:,:)   = PRUS(:,:) - PDRUS_TURB(:,:)
  PDRVS_TURB(:,:)   = PRVS(:,:) - PDRVS_TURB(:,:)
  PDRTHLS_TURB(:,:) = PRTHLS(:,:) - PDRTHLS_TURB(:,:)
  PDRRTS_TURB(:,:)  = PRRS(:,:,1) - PDRRTS_TURB(:,:)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT,JSV=1:KSV)  
  PDRSVS_TURB(:,:,:)  = PRSVS(:,:,:) - PDRSVS_TURB(:,:,:)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT,JSV=1:KSV)
END IF
!----------------------------------------------------------------------------
!
!*      8. RETRIEVE NON-CONSERVATIVE VARIABLES
!          -----------------------------------
!
IF ( KRRL >= 1 ) THEN
  IF ( KRRI >= 1 ) THEN
!$acc kernels present_cr(PRT,PRRS,PTHLT,PRTHLS)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    PRT(:,:,1)  = PRT(:,:,1)  - PRT(:,:,2)  &
                                    - PRT(:,:,4)
    PRRS(:,:,1) = PRRS(:,:,1) - PRRS(:,:,2) &
                                    - PRRS(:,:,4)
    PTHLT(:,:)  = PTHLT(:,:)  + ZLVOCPEXNM(:,:) &
                                    * PRT(:,:,2) &
                                    + ZLSOCPEXNM(:,:) * PRT(:,:,4)
    PRTHLS(:,:) = PRTHLS(:,:) + ZLVOCPEXNM(:,:) &
                                    * PRRS(:,:,2) &
                                    + ZLSOCPEXNM(:,:) * PRRS(:,:,4)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
!$acc update self(PRT(:,:,1))
!
  ELSE
!$acc kernels present_cr(PRT,PRRS,PTHLT,PRTHLS, zlocpexnm )
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    PRT(:,:,1)  = PRT(:,:,1)  - PRT(:,:,2)
    PRRS(:,:,1) = PRRS(:,:,1) - PRRS(:,:,2)
    PTHLT(:,:)  = PTHLT(:,:)  + ZLOCPEXNM(:,:) &
                                    * PRT(:,:,2)
    PRTHLS(:,:) = PRTHLS(:,:) + ZLOCPEXNM(:,:) &
                                    * PRRS(:,:,2)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
!$acc update self(PRT(:,:,1))
  END IF
END IF

! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets
CALL SOURCES_NEG_CORRECT_PHY(D,KSV,HCLOUD, 'NETUR',KRR,PTSTEP,PPABST,PTHLT,PRT,PRTHLS,PRRS,PRSVS)
!$acc update self( PTHLT ) !PTHLT not modified in Sources_neg_correct
!$acc update self( PRTHLS, PRRS )
!----------------------------------------------------------------------------
!
!*      9. LES averaged surface fluxes
!          ---------------------------
!
IF (TLES%LLES_CALL) THEN
  CALL SECOND_MNH(ZTIME1)
  !$acc data copy(TLES%X_LES_Q0,TLES%X_LES_E0,TLES%X_LES_SV0,TLES%X_LES_UW0,TLES%X_LES_VW0,TLES%X_LES_USTAR)
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,PSFTH,TLES%X_LES_Q0)
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,PSFRV,TLES%X_LES_E0)
  DO JSV=1,KSV
    CALL LES_MEAN_SUBGRID_PHY(D,TLES,PSFSV(:,JSV),TLES%X_LES_SV0(:,JSV))
  END DO
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,PSFU,TLES%X_LES_UW0)
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,PSFV,TLES%X_LES_VW0)
  !
  !$acc kernels present_cr(ZWORK2D)
  !$mnh_expand_array(JIJ=IIJB:IIJE)
  ZWORK2D(:) = (PSFU(:)*PSFU(:)+PSFV(:)*PSFV(:))**0.25
  !$mnh_end_expand_array(JIJ=IIJB:IIJE)
  !$acc end kernels
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZWORK2D,TLES%X_LES_USTAR)
!$acc end data
!----------------------------------------------------------------------------
!
!*     10. LES for 3rd order moments
!          -------------------------
!
!$acc data copy(TLES%X_LES_SUBGRID_W2Thl,TLES%X_LES_SUBGRID_WThl2)
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZMWTH,TLES%X_LES_SUBGRID_W2Thl)
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZMTH2,TLES%X_LES_SUBGRID_WThl2)
!$acc end data
  IF (KRR>0) THEN
!$acc data copy(TLES%X_LES_SUBGRID_W2Rt,TLES%X_LES_SUBGRID_WThlRt,TLES%X_LES_SUBGRID_WRt2)
    CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZMWR,TLES%X_LES_SUBGRID_W2Rt)
    CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZMTHR,TLES%X_LES_SUBGRID_WThlRt)
    CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZMR2,TLES%X_LES_SUBGRID_WRt2)
!$acc end data
  END IF
!
!----------------------------------------------------------------------------
!
!*     11. LES quantities depending on <w'2> in "1DIM" mode
!          ------------------------------------------------
!
  IF (TURBN%CTURBDIM=="1DIM") THEN
!$acc data copy(TLES%X_LES_SUBGRID_U2,TLES%X_LES_SUBGRID_V2,TLES%X_LES_SUBGRID_W2,TLES%X_LES_RES_ddz_Thl_SBG_W2)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZWORK1(:,:) = 2./3.*PTKET(:,:)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZWORK1,TLES%X_LES_SUBGRID_U2)
    TLES%X_LES_SUBGRID_V2(:,:,:) = TLES%X_LES_SUBGRID_U2(:,:,:)
    TLES%X_LES_SUBGRID_W2(:,:,:) = TLES%X_LES_SUBGRID_U2(:,:,:)
    !
    CALL GZ_M_W_PHY(D,PTHLT,PDZZ,ZWORK1)
    CALL MZF_PHY(D,ZWORK1,ZWORK2)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    ZWORK2(:,:)  = 2./3.*PTKET(:,:) *ZWORK2(:,:)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZWORK2,TLES%X_LES_RES_ddz_Thl_SBG_W2)
!$acc end data
    !
    IF (KRR>=1) THEN
!$acc data copy(TLES%X_LES_RES_ddz_Rt_SBG_W2)
      CALL GZ_M_W_PHY(D,PRT(:,:,1),PDZZ,ZWORK1)
      CALL MZF_PHY(D,ZWORK1,ZWORK2)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZWORK2(:,:)  = 2./3.*PTKET(:,:) *ZWORK2(:,:)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZWORK2,TLES%X_LES_RES_ddz_Rt_SBG_W2)
!$acc end data
    END IF
!$acc data copy(TLES%X_LES_RES_ddz_Sv_SBG_W2(:,:,:,1:KSV))    
    DO JSV=1,KSV
      CALL GZ_M_W_PHY(D,PSVT(:,:,JSV),PDZZ,ZWORK1)
      CALL MZF_PHY(D,ZWORK1,ZWORK2)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      ZWORK2(:,:)  = 2./3.*PTKET(:,:) *ZWORK2(:,:)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZWORK2, TLES%X_LES_RES_ddz_Sv_SBG_W2(:,:,:,JSV))
    END DO
!$acc end data
  END IF

!----------------------------------------------------------------------------
!
!*     12. LES mixing end dissipative lengths, presso-correlations
!          -------------------------------------------------------
!
!$acc data copy(TLES%X_LES_SUBGRID_LMix,TLES%X_LES_SUBGRID_LDiss,TLES%X_LES_SUBGRID_WP)
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZLM,TLES%X_LES_SUBGRID_LMix)
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZLEPS,TLES%X_LES_SUBGRID_LDiss)
!
!* presso-correlations for subgrid Tke are equal to zero.
!
!$acc kernels present_cr(ZLEPS)
  ZLEPS(:,:) = 0. !ZLEPS is used as a work array (not used anymore)
!$acc end kernels
  CALL LES_MEAN_SUBGRID_PHY(D,TLES,ZLEPS,TLES%X_LES_SUBGRID_WP)
!$acc end data
!
  CALL SECOND_MNH(ZTIME2)
  TLES%XTIME_LES = TLES%XTIME_LES + ZTIME2 - ZTIME1
END IF
!
IF(PRESENT(PLEM)) PLEM(:,:) = ZLM(:,:)
!----------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('TURB',1,ZHOOK_HANDLE)
CONTAINS
!
!     ########################################################################
      SUBROUTINE COMPUTE_FUNCTION_THERMO(PALP,PBETA,PGAM,PLTT,PC,PT,PEXN,PCP,&
                                         PLOCPEXN,PAMOIST,PATHETA            )
!     ########################################################################
!!
!!****  *COMPUTE_FUNCTION_THERMO* routine to compute several thermo functions
!
!!    AUTHOR
!!    ------
!!
!!     JP Pinty      *LA*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   24/02/03
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments
!
REAL,                   INTENT(IN)    :: PALP,PBETA,PGAM,PLTT,PC
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN)    :: PT,PEXN,PCP
!
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PLOCPEXN
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PAMOIST,PATHETA
!
!-------------------------------------------------------------------------------
!
  IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',0,ZHOOK_HANDLE2)
  ZEPS = CST%XMV / CST%XMD
!
!*       1.1 Lv/Cph at  t
!
!$acc kernels present_cr(PLOCPEXN)  ! present(ZRVSAT,ZDRVSATDT) ! present(PLOCPEXN) ! present ZDRVSATDT)
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  PLOCPEXN(:,:) = ( PLTT + (CST%XCPV-PC) *  (PT(:,:)-CST%XTT) ) &
                                     / PCP(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)                                     
!$acc end kernels
!
!$acc kernels present_cr(ZRVSAT,ZDRVSATDT)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!*      1.2 Saturation vapor pressure at t
!
  ZRVSAT(:,:) =  EXP( PALP - PBETA/PT(:,:) - PGAM*ALOG( PT(:,:) ) )
!
!*      1.3 saturation  mixing ratio at t
!
  ZRVSAT(:,:) =  ZRVSAT(:,:) &
                                    * ZEPS / ( PPABST(:,:) - ZRVSAT(:,:) )
!
!*      1.4 compute the saturation mixing ratio derivative (rvs')
!
  ZDRVSATDT(:,:) = ( PBETA / PT(:,:)  - PGAM ) / PT(:,:)   &
                 * ZRVSAT(:,:) * ( 1. + ZRVSAT(:,:) / ZEPS )
!
!*      1.5 compute Amoist
!
  PAMOIST(:,:)=  0.5 / ( 1.0 + ZDRVSATDT(:,:) * PLOCPEXN(:,:) )
!
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
!*      1.6 compute Atheta
!
!$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  PATHETA(:,:)= PAMOIST(:,:) * PEXN(:,:) *               &
        ( ( ZRVSAT(:,:) - PRT(:,:,1) ) * PLOCPEXN(:,:) / &
          ( 1. + ZDRVSATDT(:,:) * PLOCPEXN(:,:) )        *               &
          (                                                                  &
           ZRVSAT(:,:) * (1. + ZRVSAT(:,:)/ZEPS)                         &
                        * ( -2.*PBETA/PT(:,:) + PGAM ) / PT(:,:)**2      &
          +ZDRVSATDT(:,:) * (1. + 2. * ZRVSAT(:,:)/ZEPS)                 &
                        * ( PBETA/PT(:,:) - PGAM ) / PT(:,:)             &
          )                                                                  &
         - ZDRVSATDT(:,:)                                                  &
        )
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
!
!*      1.7 Lv/Cph/Exner at t-1
!
!$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  PLOCPEXN(:,:) = PLOCPEXN(:,:) / PEXN(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
!
IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',1,ZHOOK_HANDLE2)
END SUBROUTINE COMPUTE_FUNCTION_THERMO

!     ########################################################################
      SUBROUTINE COMPUTE_FUNCTION_THERMO_NEW_STAT(PALP,PBETA,PGAM,PLTT,PC,PT,PEXN,PCP,&
                                         PLOCPEXN,PAMOIST,PATHETA            )
!     ########################################################################
!!
!!****  *COMPUTE_FUNCTION_THERMO* routine to compute several thermo functions
!
!!    AUTHOR
!!    ------
!!
!!     JP Pinty      *LA*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   24/02/03
!!     Modified: Wim de Rooy 06-02-2019
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
USE MODD_CST
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments
!
REAL, INTENT(IN)                      :: PALP,PBETA,PGAM,PLTT,PC
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN)    :: PT,PEXN,PCP
!
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PLOCPEXN
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PAMOIST,PATHETA
!
!-------------------------------------------------------------------------------
!
  IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO_NEW_STAT',0,ZHOOK_HANDLE2)
  ZEPS = CST%XMV / CST%XMD
!
!*       1.1 Lv/Cph at  t
!
!$acc kernels
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  PLOCPEXN(:,:) = ( PLTT + (CST%XCPV-PC) *  (PT(:,:)-CST%XTT) ) / PCP(:,:)
!
!*      1.2 Saturation vapor pressure at t
!
  ZRVSAT(:,:) =  EXP( PALP - PBETA/PT(:,:) - PGAM*ALOG( PT(:,:) ) )
!
!*      1.3 saturation  mixing ratio at t
!
  ZRVSAT(:,:) =  ZRVSAT(:,:) * ZEPS / ( PPABST(:,:) - ZRVSAT(:,:) )
!
!*      1.4 compute the saturation mixing ratio derivative (rvs')
!
  ZDRVSATDT(:,:) = ( PBETA / PT(:,:)  - PGAM ) / PT(:,:)   &
                 * ZRVSAT(:,:) * ( 1. + ZRVSAT(:,:) / ZEPS )
!
!*      1.5 compute Amoist
!
  PAMOIST(:,:)=  1.0 / ( 1.0 + ZDRVSATDT(:,:) * PLOCPEXN(:,:) )
!
!*      1.6 compute Atheta
!
  PATHETA(:,:)= PAMOIST(:,:) * PEXN(:,:) * ZDRVSATDT(:,:)
!
!*      1.7 Lv/Cph/Exner at t-1
!
  PLOCPEXN(:,:) = PLOCPEXN(:,:) / PEXN(:,:)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
!
IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO_NEW_STAT',1,ZHOOK_HANDLE2)
END SUBROUTINE COMPUTE_FUNCTION_THERMO_NEW_STAT

!
!     ####################
      SUBROUTINE DELT(PLM,ODZ)
!     ####################
!!
!!****  *DELT* routine to compute mixing length for DELT case
!
!!    AUTHOR
!!    ------
!!
!!     M Tomasini      *Meteo-France
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   01/05
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
!*       0.1   Declarations of dummy arguments
!
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PLM
LOGICAL,                INTENT(IN)    :: ODZ
!-------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('TURB:DELT',0,ZHOOK_HANDLE2)
!
CALL MXF_PHY(D,PDXX,ZWORK1)
IF (.NOT. O2D) THEN
  CALL MYF_PHY(D,PDYY,ZWORK2)
END IF
!
IF (ODZ) THEN
!$acc kernels present_cr(PLM)
  ! Dz is take into account in the computation
  DO JK = IKTB,IKTE ! 1D turbulence scheme
    !$mnh_expand_array(JIJ=IIJB:IIJE)
    PLM(:,JK) = PZZ(:,JK+IKL) - PZZ(:,JK)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE)
  END DO
  !$mnh_expand_array(JIJ=IIJB:IIJE)
  PLM(:,IKU) = PLM(:,IKE)
  PLM(:,IKA) = PZZ(:,IKB) - PZZ(:,IKA)
  !$mnh_end_expand_array(JIJ=IIJB:IIJE)
!$acc end kernels
  IF ( TURBN%CTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
    IF ( O2D) THEN
!$acc kernels present_cr(PLM)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      PLM(:,:) = SQRT( PLM(:,:)*ZWORK1(:,:) )
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
    ELSE
!$acc kernels present_cr(PLM)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      PLM(:,:) = (PLM(:,:)*ZWORK1(:,:) &
                                   * ZWORK2(:,:) ) ** (1./3.)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
    END IF
  END IF
ELSE
  ! Dz not taken into account in computation to assure invariability with vertical grid mesh
!$acc kernels present_cr(PLM)
  !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
  PLM(:,:)=1.E10
  !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
!$acc end kernels
  IF ( TURBN%CTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
    IF ( O2D) THEN
      !$acc kernels present_cr(PLM)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      PLM(:,:) = ZWORK1(:,:)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
    ELSE
      !$acc kernels present_cr(PLM)
      !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      PLM(:,:) = (ZWORK1(:,:)*ZWORK2(:,:) ) ** (1./2.)
      !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
      !$acc end kernels
    END IF
  END IF
END IF
!
!  mixing length limited by the distance normal to the surface
!  (with the same factor as for BL89)
!
IF (.NOT. TURBN%LRMC01) THEN
  ZALPHA=0.5**(-1.5)
  !
!$acc kernels
  DO JIJ=IIJB,IIJE
    IF (GOCEAN) THEN
      DO JK=IKTE,IKTB,-1
        ZD=ZALPHA*(PZZ(JIJ,IKTE+1)-PZZ(JIJ,JK))
        IF ( PLM(JIJ,JK)>ZD) THEN
          PLM(JIJ,JK)=ZD
        ELSE
          EXIT
        ENDIF
      END DO
    ELSE
      DO JK=IKTB,IKTE
        ZD=ZALPHA*(0.5*(PZZ(JIJ,JK)+PZZ(JIJ,JK+IKL))&
        -PZZ(JIJ,IKB)) *PDIRCOSZW(JIJ)
        IF ( PLM(JIJ,JK)>ZD) THEN
          PLM(JIJ,JK)=ZD
        ELSE
          EXIT
        ENDIF
      END DO
    ENDIF
  END DO
!$acc end kernels
END IF
!
!$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE)
PLM(:,IKA) = PLM(:,IKB)
PLM(:,IKU) = PLM(:,IKE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!$acc end kernels
!
IF (LHOOK) CALL DR_HOOK('TURB:DELT',1,ZHOOK_HANDLE2)
END SUBROUTINE DELT
!
!     ####################
      SUBROUTINE DEAR(PLM)
!     ####################
!!
!!****  *DEAR* routine to compute mixing length for DEARdorff case
!
!!    AUTHOR
!!    ------
!!
!!     M Tomasini      *Meteo-France
!!
!!    MODIFICATIONS
!!    -------------
!!      Original   01/05
!!      I.Sandu (Sept.2006) : Modification of the stability criterion
!!                            (theta_v -> theta_l)
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
!*       0.1   Declarations of dummy arguments
!
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(OUT)   :: PLM
!
!-------------------------------------------------------------------------------
!
!   initialize the mixing length with the mesh grid
IF (LHOOK) CALL DR_HOOK('TURB:DEAR',0,ZHOOK_HANDLE2)
IF ( TURBN%CTURBDIM /= '1DIM' ) THEN
  CALL MXF_PHY(D,PDXX,ZWORK1)
  IF (.NOT. O2D) THEN
    CALL MYF_PHY(D,PDYY,ZWORK2)
  END IF
END IF
! 1D turbulence scheme
!$acc kernels present_cr(PLM)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=IKTB:IKTE)
PLM(:,IKTB:IKTE) = PZZ(:,IKTB+IKL:IKTE+IKL) - PZZ(:,IKTB:IKTE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=IKTB:IKTE)
!$mnh_expand_array(JIJ=IIJB:IIJE)
PLM(:,IKU) = PLM(:,IKE)
PLM(:,IKA) = PZZ(:,IKB) - PZZ(:,IKA)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
!$acc end kernels
!
IF ( TURBN%CTURBDIM /= '1DIM' ) THEN  ! 3D turbulence scheme
  IF ( O2D) THEN
    !$acc kernels present_cr(PLM)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    PLM(:,:) = SQRT( PLM(:,:)*ZWORK1(:,:) )
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
  ELSE
    !$acc kernels present_cr(PLM)
    !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    PLM(:,:) = (PLM(:,:)*ZWORK1(:,:) &
                                 * ZWORK2(:,:) ) ** (1./3.)
    !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
    !$acc end kernels
  END IF
END IF
!   compute a mixing length limited by the stability
!
CALL ETHETA(D,CST,KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZATHETA,PSRCT,GOCEAN,OCOMPUTE_SRC,ZETHETA)
CALL EMOIST(D,CST,KRR,KRRI,PTHLT,PRT,ZLOCPEXNM,ZAMOIST,PSRCT,GOCEAN,ZEMOIST)
!
IF (KRR>0) THEN
!$acc kernels
!$acc loop independent collapse(2)
  DO CONCURRENT (JK=IKTB+1:IKTE-1,JIJ=IIJB:IIJE)
      ZDTHLDZ(JIJ,JK)= 0.5*((PTHLT(JIJ,JK+IKL)-PTHLT(JIJ,JK    ))/PDZZ(JIJ,JK+IKL)+ &
                              (PTHLT(JIJ,JK    )-PTHLT(JIJ,JK-IKL))/PDZZ(JIJ,JK    ))
      ZDRTDZ(JIJ,JK) = 0.5*((PRT(JIJ,JK+IKL,1)-PRT(JIJ,JK    ,1))/PDZZ(JIJ,JK+IKL)+ &
                              (PRT(JIJ,JK    ,1)-PRT(JIJ,JK-IKL,1))/PDZZ(JIJ,JK    ))