Quentin 02/10/2023: ice4_nucleation written as ELEMENTAL included file in...

Quentin 02/10/2023: ice4_nucleation written as ELEMENTAL included file in order to be called either on full arrays (in rain_ice) or packed arrays (in mode_ice4_tendencies)

src/common/micro/mode_ice4_nucleation.F90 → src/common/micro/ice4_nucleation.func.h

@@ -3,10 +3,7 @@
 !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(CST, PARAMI, ICEP, ICED, KSIZE, ODCOMPUTE, &
+ELEMENTAL SUBROUTINE ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, ODCOMPUTE, &
                            PTHT, PPABST, PRHODREF, PEXN, PLSFACT, PT, &
                            PRVT, &
                            PCIT, PRVHENI_MR)
@@ -32,7 +29,6 @@ USE MODD_CST,            ONLY: CST_t
 USE MODD_PARAM_ICE_n,      ONLY: PARAM_ICE_t
 USE MODD_RAIN_ICE_DESCR_n, ONLY: RAIN_ICE_DESCR_t
 USE MODD_RAIN_ICE_PARAM_n, ONLY: RAIN_ICE_PARAM_t
-USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
 !
 IMPLICIT NONE
 !
@@ -42,114 +38,86 @@ TYPE(CST_t),              INTENT(IN)    :: CST
 TYPE(PARAM_ICE_t),        INTENT(IN)    :: PARAMI
 TYPE(RAIN_ICE_PARAM_t),   INTENT(IN)    :: ICEP
 TYPE(RAIN_ICE_DESCR_t),   INTENT(IN)    :: ICED
-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
+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, DIMENSION(KSIZE) :: ZW ! work array
-REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
-LOGICAL, DIMENSION(KSIZE) :: GNEGT  ! Test where to compute the HEN process
-REAL, DIMENSION(KSIZE)  :: ZZW,      & ! Work array
+REAL  :: ZW ! work array
+LOGICAL  :: GNEGT  ! Test where to compute the HEN process
+REAL   :: ZZW,      & ! Work array
                            ZUSW,     & ! Undersaturation over water
                            ZSSI        ! Supersaturation over ice
-INTEGER :: JL
 !-------------------------------------------------------------------------------
 !
-IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION', 0, ZHOOK_HANDLE)!
 !
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(ODCOMPUTE(:))
-  GNEGT(:)=PT(:)<CST%XTT .AND. PRVT(:)>ICED%XRTMIN(1)
-ELSEWHERE
-  GNEGT(:)=.FALSE.
-ENDWHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
+  IF (ODCOMPUTE) THEN
+GNEGT=PT<CST%XTT .AND. PRVT>ICED%XRTMIN(1)
+ELSE
+GNEGT=.FALSE.
+END IF
 
-ZUSW(:)=0.
-ZZW(:)=0.
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  ZZW(:)=ALOG(PT(:))
-  ZUSW(:)=EXP(CST%XALPW - CST%XBETAW/PT(:) - CST%XGAMW*ZZW(:))          ! es_w
-  ZZW(:)=EXP(CST%XALPI - CST%XBETAI/PT(:) - CST%XGAMI*ZZW(:))           ! es_i
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
+ZUSW=0.
+ZZW=0.
 
-ZSSI(:)=0.
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  ZZW(:)=MIN(PPABST(:)/2., ZZW(:))             ! safety limitation
-  ZSSI(:)=PRVT(:)*(PPABST(:)-ZZW(:)) / (CST%XEPSILO*ZZW(:)) - 1.0
+IF (GNEGT) THEN
+ ZZW=ALOG(PT)
+ ZUSW=EXP(CST%XALPW - CST%XBETAW/PT - CST%XGAMW*ZZW)          ! es_w
+ ZZW=EXP(CST%XALPI - CST%XBETAI/PT - CST%XGAMI*ZZW)           ! es_i
+END IF
+  
+ZSSI=0.
+IF (GNEGT) THEN
+  ZZW=MIN(PPABST/2., ZZW)             ! safety limitation
+  ZSSI=PRVT*(PPABST-ZZW) / (CST%XEPSILO*ZZW) - 1.0
                                                ! Supersaturation over ice
-  ZUSW(:)=MIN(PPABST(:)/2., ZUSW(:))            ! safety limitation
-  ZUSW(:)=(ZUSW(:)/ZZW(:))*((PPABST(:)-ZZW(:))/(PPABST(:)-ZUSW(:))) - 1.0
+  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
-!$mnh_end_expand_where(JL=1:KSIZE)
+  ZSSI=MIN(ZSSI, ZUSW) ! limitation of SSi according to SSw=0
+END IF
+  
+ZZW=0.
 
-ZZW(:)=0.
-DO JL=1,KSIZE
-  IF(GNEGT(JL)) THEN
-    IF(PT(JL)<CST%XTT-5.0 .AND. ZSSI(JL)>0.0) THEN
-      ZZW(JL)=ICEP%XNU20*EXP(ICEP%XALPHA2*ZSSI(JL)-ICEP%XBETA2)
-    ELSEIF(PT(JL)<=CST%XTT-2.0 .AND. PT(JL)>=CST%XTT-5.0 .AND. ZSSI(JL)>0.0) THEN
-      ZZW(JL)=MAX(ICEP%XNU20*EXP(-ICEP%XBETA2 ), &                                                                                       
-                  ICEP%XNU10*EXP(-ICEP%XBETA1*(PT(JL)-CST%XTT))*(ZSSI(JL)/ZUSW(JL))**ICEP%XALPHA1)
-    ENDIF
+IF(GNEGT) THEN
+  IF(PT<CST%XTT-5.0 .AND. ZSSI>0.0) THEN
+    ZZW=ICEP%XNU20*EXP(ICEP%XALPHA2*ZSSI-ICEP%XBETA2)
+  ELSEIF(PT<=CST%XTT-2.0 .AND. PT>=CST%XTT-5.0 .AND. ZSSI>0.0) THEN
+    ZZW=MAX(ICEP%XNU20*EXP(-ICEP%XBETA2 ), &                                                                                       
+                ICEP%XNU10*EXP(-ICEP%XBETA1*(PT-CST%XTT))*(ZSSI/ZUSW)**ICEP%XALPHA1)
   ENDIF
-ENDDO
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  ZZW(:)=ZZW(:)-PCIT(:)
-  ZZW(:)=MIN(ZZW(:), 50.E3) ! limitation provisoire a 50 l^-1
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
+ENDIF
+IF (GNEGT) THEN
+  ZZW=ZZW-PCIT
+  ZZW=MIN(ZZW, 50.E3) ! limitation provisoire a 50 l^-1
+END IF
 
-PRVHENI_MR(:)=0.
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  !
-  !*       3.1.2   update the r_i and r_v mixing ratios
-  !
-  PRVHENI_MR(:)=MAX(ZZW(:), 0.0)*ICEP%XMNU0/PRHODREF(:)
-  PRVHENI_MR(:)=MIN(PRVT(:), PRVHENI_MR(:))
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
+PRVHENI_MR=0.
+
+IF (GNEGT) THEN
+  PRVHENI_MR=MAX(ZZW, 0.0)*ICEP%XMNU0/PRHODREF
+  PRVHENI_MR=MIN(PRVT, PRVHENI_MR)
+END IF
 !Limitation due to 0 crossing of temperature
 IF(PARAMI%LFEEDBACKT) THEN
-  ZW(:)=0.
-  !$mnh_expand_where(JL=1:KSIZE)
-  WHERE(GNEGT(:))
-    ZW(:)=MIN(PRVHENI_MR(:), &
-              MAX(0., (CST%XTT/PEXN(:)-PTHT(:))/PLSFACT(:))) / &
-              MAX(PRVHENI_MR(:), 1.E-20)
-  END WHERE
-  PRVHENI_MR(:)=PRVHENI_MR(:)*ZW(:)
-  ZZW(:)=ZZW(:)*ZW(:)
-  !$mnh_end_expand_where(JL=1:KSIZE)
+  ZW=0.
+  IF (GNEGT) THEN
+    ZW=MIN(PRVHENI_MR, &
+              MAX(0., (CST%XTT/PEXN-PTHT)/PLSFACT)) / &
+              MAX(PRVHENI_MR, 1.E-20)
+  END IF
+  PRVHENI_MR=PRVHENI_MR*ZW
+  ZZW=ZZW*ZW
 ENDIF
-!$mnh_expand_where(JL=1:KSIZE)
-WHERE(GNEGT(:))
-  PCIT(:)=MAX(ZZW(:)+PCIT(:), PCIT(:))
-END WHERE
-!$mnh_end_expand_where(JL=1:KSIZE)
-!
-IF (LHOOK) CALL DR_HOOK('ICE4_NUCLEATION', 1, ZHOOK_HANDLE)
+IF (GNEGT) THEN
+  PCIT=MAX(ZZW+PCIT, PCIT)
+END IF
 END SUBROUTINE ICE4_NUCLEATION
-END MODULE MODE_ICE4_NUCLEATION

src/common/micro/mode_ice4_tendencies.F90

@@ -54,7 +54,6 @@ 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 YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
 !
@@ -135,10 +134,10 @@ ELSE
   !
   !*       2.     COMPUTES THE SLOW COLD PROCESS SOURCES
   !               --------------------------------------
-  CALL ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, KSIZE, LDCOMPUTE(:), &
-                   ZVART(:,ITH), PPRES(:), PRHODREF(:), PEXN(:), PLSFACT(:), ZT(:), &
-                   ZVART(:,IRV), &
-                   PCIT(:), PBU_INST(:, IRVHENI_MR))
+  CALL ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, LDCOMPUTE(1:KSIZE), &
+                   ZVART(1:KSIZE,ITH), PPRES(1:KSIZE), PRHODREF(1:KSIZE), PEXN(1:KSIZE), PLSFACT(1:KSIZE), ZT(1:KSIZE), &
+                   ZVART(1:KSIZE,IRV), &
+                   PCIT(1:KSIZE), PBU_INST(1:KSIZE, IRVHENI_MR))
   DO JL=1, KSIZE
     ZVART(JL,ITH)=ZVART(JL,ITH) + PBU_INST(JL, IRVHENI_MR)*PLSFACT(JL)
     ZT(JL) = ZVART(JL,ITH) * PEXN(JL)
@@ -549,5 +548,7 @@ ENDDO
 !
 IF (LHOOK) CALL DR_HOOK('ICE4_TENDENCIES', 1, ZHOOK_HANDLE)
 !
+CONTAINS
+INCLUDE "ice4_nucleation.func.h"
 END SUBROUTINE ICE4_TENDENCIES
 END MODULE MODE_ICE4_TENDENCIES

src/common/micro/rain_ice.F90

@@ -203,7 +203,6 @@ USE MODE_ICE4_RAINFR_VERT,   ONLY: ICE4_RAINFR_VERT
 USE MODE_ICE4_COMPUTE_PDF,   ONLY: ICE4_COMPUTE_PDF
 USE MODE_ICE4_SEDIMENTATION, ONLY: ICE4_SEDIMENTATION
 USE MODE_ICE4_PACK, ONLY: ICE4_PACK
-USE MODE_ICE4_NUCLEATION, ONLY: ICE4_NUCLEATION
 USE MODE_ICE4_CORRECT_NEGATIVITIES, ONLY: ICE4_CORRECT_NEGATIVITIES
 !
 IMPLICIT NONE
@@ -418,7 +417,7 @@ DO JK=IKTB,IKTE
     ENDIF
   ENDDO
 ENDDO
-CALL ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, D%NIJT*D%NKT, LLW3D(:,:), &
+CALL ICE4_NUCLEATION(CST, PARAMI, ICEP, ICED, LLW3D(:,:), &
                      PTHT(:, :), PPABST(:, :), PRHODREF(:, :), &                                       
                      PEXN(:, :), ZW3D(:, :), ZT(:, :), &                                                           
                      PRVT(:, :), &                                                                                 
@@ -669,4 +668,6 @@ ENDIF
 
 IF (LHOOK) CALL DR_HOOK('RAIN_ICE', 1, ZHOOK_HANDLE)
 !
+CONTAINS
+INCLUDE "ice4_nucleation.func.h"
 END SUBROUTINE RAIN_ICE