From a5462c37e1184346adeb37a7d09ef5d999b9e3ca Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Riette?= <sebastien.riette@meteo.fr>
Date: Mon, 23 May 2022 18:24:50 +0200
Subject: [PATCH] S. Riette 23/05/2022 Revert ice4_nucleation
For GPU, ice4_nucleation cannot be put in an include file
as included subroutine must not contain local variables.
---
src/common/micro/ice4_nucleation_elem.func.h | 110 ----------------
src/common/micro/mode_ice4_nucleation.F90 | 127 +++++++++++++++++++
src/common/micro/mode_ice4_tendencies.F90 | 14 +-
src/common/micro/rain_ice.F90 | 27 ++--
4 files changed, 147 insertions(+), 131 deletions(-)
delete mode 100644 src/common/micro/ice4_nucleation_elem.func.h
create mode 100644 src/common/micro/mode_ice4_nucleation.F90
diff --git a/src/common/micro/ice4_nucleation_elem.func.h b/src/common/micro/ice4_nucleation_elem.func.h
deleted file mode 100644
index f9ef73991..000000000
--- a/src/common/micro/ice4_nucleation_elem.func.h
+++ /dev/null
@@ -1,110 +0,0 @@
-!MNH_LIC Copyright 1994-2021 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.
-ELEMENTAL SUBROUTINE ICE4_NUCLEATION_ELEM(ODCOMPUTE, &
- PTHT, PPABST, PRHODREF, PEXN, PLSFACT, PT, &
- PRVT, &
- PCIT, PRVHENI_MR)
-! ******* TO BE INCLUDED IN THE *CONTAINS* OF A SUBROUTINE, IN ORDER TO EASE AUTOMATIC INLINING ******
-! => Don't use drHook !!!
-!
-
-!!
-!!** PURPOSE
-!! -------
-!! Computes the nucleation
-!!
-!! AUTHOR
-!! ------
-!! S. Riette from the splitting of rain_ice source code (nov. 2014)
-!!
-!! MODIFICATIONS
-!! -------------
-!!
-!! R. El Khatib 24-Aug-2021 Optimizations
-!! S. Riette Feb 2022: as an include file
-!
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_CST, ONLY: XALPI, XALPW, XBETAI, XBETAW, XGAMI, XGAMW, XMD, XMV, XTT, XEPSILO
-USE MODD_PARAM_ICE, ONLY: LFEEDBACKT
-USE MODD_RAIN_ICE_PARAM, ONLY: XALPHA1, XALPHA2, XBETA1, XBETA2, XMNU0, XNU10, XNU20
-USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN
-!
-IMPLICIT NONE
-!
-!* 0.1 Declarations of dummy arguments :
-!
-LOGICAL, INTENT(IN) :: ODCOMPUTE
-REAL, INTENT(IN) :: PTHT ! Theta at t
-REAL, INTENT(IN) :: PPABST ! absolute pressure at t
-REAL, INTENT(IN) :: PRHODREF! Reference density
-REAL, INTENT(IN) :: PEXN ! Exner function
-REAL, INTENT(IN) :: PLSFACT
-REAL, INTENT(IN) :: PT ! Temperature at time t
-REAL, INTENT(IN) :: PRVT ! Water vapor m.r. at t
-REAL, INTENT(INOUT) :: PCIT ! Pristine ice n.c. at t
-REAL, INTENT(OUT) :: PRVHENI_MR ! Mixing ratio change due to the heterogeneous nucleation
-!
-!* 0.2 declaration of local variables
-!
-REAL :: ZW ! work array
-LOGICAL :: GNEGT ! Test where to compute the HEN process
-REAL :: ZZW, & ! Work scalar
- ZUSW, & ! Undersaturation over water
- ZSSI ! Supersaturation over ice
-!-------------------------------------------------------------------------------
-!
-GNEGT=PT<XTT .AND. PRVT>XRTMIN(1) .AND. ODCOMPUTE
-
-PRVHENI_MR=0.
-IF(GNEGT) THEN
- ZZW=ALOG(PT)
- ZUSW=EXP(XALPW - XBETAW/PT - XGAMW*ZZW) ! es_w
- ZZW=EXP(XALPI - XBETAI/PT - XGAMI*ZZW) ! es_i
-
- ZZW=MIN(PPABST/2., ZZW) ! safety limitation
- ZSSI=PRVT*(PPABST-ZZW) / (XEPSILO*ZZW) - 1.0 ! Supersaturation over ice
- ZUSW=MIN(PPABST/2., ZUSW) ! safety limitation
- ZUSW=(ZUSW/ZZW)*((PPABST-ZZW)/(PPABST-ZUSW)) - 1.0
- ! Supersaturation of saturated water vapor over ice
- !
- !* 3.1 compute the heterogeneous nucleation source RVHENI
- !
- !* 3.1.1 compute the cloud ice concentration
- !
- ZSSI=MIN(ZSSI, ZUSW) ! limitation of SSi according to SSw=0
-
- IF(PT<XTT-5. .AND. ZSSI>0.) THEN
- ZZW=XNU20*EXP(XALPHA2*ZSSI-XBETA2)
- ELSEIF(PT<=XTT-2. .AND. PT>=XTT-5. .AND. ZSSI>0.) THEN
- ZZW=MAX(XNU20*EXP(-XBETA2 ), &
- XNU10*EXP(-XBETA1*(PT-XTT))*(ZSSI/ZUSW)**XALPHA1)
- ELSE
- ZZW=0.
- ENDIF
-
- ZZW=ZZW-PCIT
- ZZW=MIN(ZZW, 50.E3) ! limitation provisoire a 50 l^-1
- !
- !* 3.1.2 update the r_i and r_v mixing ratios
- !
- PRVHENI_MR=MAX(ZZW, 0.0)*XMNU0/PRHODREF
- PRVHENI_MR=MIN(PRVT, PRVHENI_MR)
- !
- !Limitation due to 0 crossing of temperature
- !
- IF(LFEEDBACKT) THEN
- ZW=MIN(PRVHENI_MR, MAX(0., (XTT/PEXN-PTHT)/PLSFACT)) / &
- MAX(PRVHENI_MR, 1.E-20)
- PRVHENI_MR=PRVHENI_MR*ZW
- ZZW=ZZW*ZW
- ENDIF
- !
- PCIT=MAX(ZZW+PCIT, PCIT)
-ENDIF
-!
-END SUBROUTINE ICE4_NUCLEATION_ELEM
diff --git a/src/common/micro/mode_ice4_nucleation.F90 b/src/common/micro/mode_ice4_nucleation.F90
new file mode 100644
index 000000000..7d5423327
--- /dev/null
+++ b/src/common/micro/mode_ice4_nucleation.F90
@@ -0,0 +1,127 @@
+!MNH_LIC Copyright 1994-2021 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 MODE_ICE4_NUCLEATION
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE ICE4_NUCLEATION(KSIZE, ODCOMPUTE, &
+ PTHT, PPABST, PRHODREF, PEXN, PLSFACT, PT, &
+ PRVT, &
+ PCIT, PRVHENI_MR)
+!!
+!!** PURPOSE
+!! -------
+!! Computes the nucleation
+!!
+!! AUTHOR
+!! ------
+!! S. Riette from the splitting of rain_ice source code (nov. 2014)
+!!
+!! MODIFICATIONS
+!! -------------
+!!
+!! R. El Khatib 24-Aug-2021 Optimizations
+!
+!
+!* 0. DECLARATIONS
+! ------------
+!
+USE MODD_CST, ONLY: XALPI, XALPW, XBETAI, XBETAW, XGAMI, XGAMW, XMD, XMV, XTT, XEPSILO
+USE MODD_PARAM_ICE, ONLY: LFEEDBACKT
+USE MODD_RAIN_ICE_PARAM, ONLY: XALPHA1, XALPHA2, XBETA1, XBETA2, XMNU0, XNU10, XNU20
+USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+INTEGER, INTENT(IN) :: KSIZE
+LOGICAL, DIMENSION(KSIZE),INTENT(IN) :: ODCOMPUTE
+REAL, DIMENSION(KSIZE), INTENT(IN) :: PTHT ! Theta at t
+REAL, DIMENSION(KSIZE), INTENT(IN) :: PPABST ! absolute pressure at t
+REAL, DIMENSION(KSIZE), INTENT(IN) :: PRHODREF! Reference density
+REAL, DIMENSION(KSIZE), INTENT(IN) :: PEXN ! Exner function
+REAL, DIMENSION(KSIZE), INTENT(IN) :: PLSFACT
+REAL, DIMENSION(KSIZE), INTENT(IN) :: PT ! Temperature at time t
+REAL, DIMENSION(KSIZE), INTENT(IN) :: PRVT ! Water vapor m.r. at t
+REAL, DIMENSION(KSIZE), INTENT(INOUT) :: PCIT ! Pristine ice n.c. at t
+REAL, DIMENSION(KSIZE), INTENT(OUT) :: PRVHENI_MR ! Mixing ratio change due to the heterogeneous nucleation
+!
+!* 0.2 declaration of local variables
+!
+REAL, DIMENSION(KSIZE) :: ZW ! work array
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+LOGICAL, DIMENSION(KSIZE) :: GNEGT ! Test where to compute the HEN process
+REAL, DIMENSION(KSIZE) :: ZZW, & ! Work array
+ ZUSW, & ! Undersaturation over water
+ ZSSI ! Supersaturation over ice
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION', 0, ZHOOK_HANDLE)!
+!
+GNEGT(:)=PT(:)<XTT .AND. PRVT(:)>XRTMIN(1) .AND. ODCOMPUTE(:)
+
+ZUSW(:)=0.
+ZZW(:)=0.
+WHERE(GNEGT(:))
+ ZZW(:)=ALOG(PT(:))
+ ZUSW(:)=EXP(XALPW - XBETAW/PT(:) - XGAMW*ZZW(:)) ! es_w
+ ZZW(:)=EXP(XALPI - XBETAI/PT(:) - XGAMI*ZZW(:)) ! es_i
+END WHERE
+
+ZSSI(:)=0.
+WHERE(GNEGT(:))
+ ZZW(:)=MIN(PPABST(:)/2., ZZW(:)) ! safety limitation
+ ZSSI(:)=PRVT(:)*(PPABST(:)-ZZW(:)) / (XEPSILO*ZZW(:)) - 1.0
+ ! Supersaturation over ice
+ ZUSW(:)=MIN(PPABST(:)/2., ZUSW(:)) ! safety limitation
+ ZUSW(:)=(ZUSW(:)/ZZW(:))*((PPABST(:)-ZZW(:))/(PPABST(:)-ZUSW(:))) - 1.0
+ ! Supersaturation of saturated water vapor over ice
+ !
+ !* 3.1 compute the heterogeneous nucleation source RVHENI
+ !
+ !* 3.1.1 compute the cloud ice concentration
+ !
+ ZSSI(:)=MIN(ZSSI(:), ZUSW(:)) ! limitation of SSi according to SSw=0
+END WHERE
+
+ZZW(:)=0.
+WHERE(GNEGT(:) .AND. PT(:)<XTT-5.0 .AND. ZSSI(:)>0.0 )
+ ZZW(:)=XNU20*EXP(XALPHA2*ZSSI(:)-XBETA2)
+ELSEWHERE(GNEGT(:) .AND. PT(:)<=XTT-2.0 .AND. PT(:)>=XTT-5.0 .AND. ZSSI(:)>0.0)
+ ZZW(:)=MAX(XNU20*EXP(-XBETA2 ), &
+ XNU10*EXP(-XBETA1*(PT(:)-XTT))*(ZSSI(:)/ZUSW(:))**XALPHA1)
+END WHERE
+WHERE(GNEGT(:))
+ ZZW(:)=ZZW(:)-PCIT(:)
+ ZZW(:)=MIN(ZZW(:), 50.E3) ! limitation provisoire a 50 l^-1
+END WHERE
+
+PRVHENI_MR(:)=0.
+WHERE(GNEGT(:))
+ !
+ !* 3.1.2 update the r_i and r_v mixing ratios
+ !
+ PRVHENI_MR(:)=MAX(ZZW(:), 0.0)*XMNU0/PRHODREF(:)
+ PRVHENI_MR(:)=MIN(PRVT(:), PRVHENI_MR(:))
+END WHERE
+!Limitation due to 0 crossing of temperature
+IF(LFEEDBACKT) THEN
+ ZW(:)=0.
+ WHERE(GNEGT(:))
+ ZW(:)=MIN(PRVHENI_MR(:), &
+ MAX(0., (XTT/PEXN(:)-PTHT(:))/PLSFACT(:))) / &
+ MAX(PRVHENI_MR(:), 1.E-20)
+ END WHERE
+ PRVHENI_MR(:)=PRVHENI_MR(:)*ZW(:)
+ ZZW(:)=ZZW(:)*ZW(:)
+ENDIF
+PCIT(:)=MAX(ZZW(:)+PCIT(:), PCIT(:))
+!
+IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION', 1, ZHOOK_HANDLE)
+END SUBROUTINE ICE4_NUCLEATION
+END MODULE MODE_ICE4_NUCLEATION
diff --git a/src/common/micro/mode_ice4_tendencies.F90 b/src/common/micro/mode_ice4_tendencies.F90
index 047e481c5..e00c28266 100644
--- a/src/common/micro/mode_ice4_tendencies.F90
+++ b/src/common/micro/mode_ice4_tendencies.F90
@@ -75,6 +75,7 @@ USE MODE_ICE4_FAST_RS, ONLY: ICE4_FAST_RS
USE MODE_ICE4_FAST_RG, ONLY: ICE4_FAST_RG
USE MODE_ICE4_FAST_RH, ONLY: ICE4_FAST_RH
USE MODE_ICE4_FAST_RI, ONLY: ICE4_FAST_RI
+USE MODE_ICE4_NUCLEATION, ONLY: ICE4_NUCLEATION
!
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
@@ -206,13 +207,10 @@ ELSE
!
!* 2. COMPUTES THE SLOW COLD PROCESS SOURCES
! --------------------------------------
-!DIR$ VECTOR ALWAYS
- DO CONCURRENT (JL=1:KSIZE)
- CALL ICE4_NUCLEATION_ELEM(LLCOMPUTE(JL), &
- ZVART(JL,ITH), PPRES(JL), PRHODREF(JL), PEXN(JL), PLSFACT(JL), ZT(JL), &
- ZVART(JL,IRV), &
- PCIT(JL), PRVHENI_MR(JL))
- ENDDO
+ CALL ICE4_NUCLEATION(KSIZE, LLCOMPUTE(:), &
+ ZVART(:,ITH), PPRES(:), PRHODREF(:), PEXN(:), PLSFACT(:), ZT(:), &
+ ZVART(:,IRV), &
+ PCIT(:), PRVHENI_MR(:))
DO JL=1, KSIZE
ZVART(JL,ITH)=ZVART(JL,ITH) + PRVHENI_MR(JL)*PLSFACT(JL)
ZT(JL) = ZVART(JL,ITH) * PEXN(JL)
@@ -492,7 +490,5 @@ CALL ICE4_FAST_RI(KSIZE, ODSOFT, PCOMPUTE, &
!
IF (LHOOK) CALL DR_HOOK('ICE4_TENDENCIES', 1, ZHOOK_HANDLE)
!
-CONTAINS
-INCLUDE "ice4_nucleation_elem.func.h"
END SUBROUTINE ICE4_TENDENCIES
END MODULE MODE_ICE4_TENDENCIES
diff --git a/src/common/micro/rain_ice.F90 b/src/common/micro/rain_ice.F90
index 9894bc0bf..3f2e20f3b 100644
--- a/src/common/micro/rain_ice.F90
+++ b/src/common/micro/rain_ice.F90
@@ -207,6 +207,7 @@ USE MODE_ICE4_SEDIMENTATION_STAT, ONLY: ICE4_SEDIMENTATION_STAT
USE MODE_ICE4_SEDIMENTATION_SPLIT, ONLY: ICE4_SEDIMENTATION_SPLIT
USE MODE_ICE4_SEDIMENTATION_SPLIT_MOMENTUM, ONLY: ICE4_SEDIMENTATION_SPLIT_MOMENTUM
USE MODE_ICE4_TENDENCIES, ONLY: ICE4_TENDENCIES
+USE MODE_ICE4_NUCLEATION, ONLY: ICE4_NUCLEATION
!
IMPLICIT NONE
!
@@ -434,7 +435,9 @@ REAL :: ZDEVIDE, ZX, ZRICE
!
INTEGER :: IC, JMICRO
LOGICAL :: LLSIGMA_RC, LL_ANY_ITER, LL_AUCV_ADJU
-
+!
+REAL, DIMENSION(KIT,KJT,KKT) :: ZW3D
+LOGICAL, DIMENSION(KIT,KJT,KKT) :: LLW3D
!
!-------------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('RAIN_ICE', 0, ZHOOK_HANDLE)
@@ -1115,21 +1118,23 @@ PCIT(:,:,:)=ZCITOUT(:,:,:)
!* 6. COMPUTES THE SLOW COLD PROCESS SOURCES OUTSIDE OF ODMICRO POINTS
! ----------------------------------------------------------------
!
-DO JK=1, KKT
+DO JK=1, KKT
DO JJ=1, KJT
-!DIR$ VECTOR ALWAYS
- DO CONCURRENT (JI=1:KIT)
+ DO JI=1, KIT
IF (.NOT. ODMICRO(JI, JJ, JK)) THEN
- ZW0D=ZZ_LSFACT(JI, JJ, JK)/PEXN(JI, JJ, JK)
+ LLW3D(JI, JJ, JK)=.TRUE.
+ ZW3D(JI, JJ, JK)=ZZ_LSFACT(JI, JJ, JK)/PEXN(JI, JJ, JK)
+ ELSE
+ LLW3D(JI, JJ, JK)=.FALSE.
ENDIF
- CALL ICE4_NUCLEATION_ELEM(.NOT. ODMICRO(JI, JJ, JK), &
- PTHT(JI, JJ, JK), PPABST(JI, JJ, JK), PRHODREF(JI, JJ, JK), &
- PEXN(JI, JJ, JK), ZW0D, ZT(JI, JJ, JK), &
- PRVT(JI, JJ, JK), &
- PCIT(JI, JJ, JK), ZZ_RVHENI_MR(JI, JJ, JK))
ENDDO
ENDDO
ENDDO
+CALL ICE4_NUCLEATION(KIT*KJT*KKT, LLW3D(:,:,:), &
+ PTHT(:, :, :), PPABST(:, :, :), PRHODREF(:, :, :), &
+ PEXN(:, :, :), ZW3D(:, :, :), ZT(:, :, :), &
+ PRVT(:, :, :), &
+ PCIT(:, :, :), ZZ_RVHENI_MR(:, :, :))
!
!-------------------------------------------------------------------------------
!
@@ -1808,6 +1813,4 @@ CONTAINS
!
END SUBROUTINE CORRECT_NEGATIVITIES
!
-INCLUDE "ice4_nucleation_elem.func.h"
-!
END SUBROUTINE RAIN_ICE
--
GitLab