Skip to content
Snippets Groups Projects
Commit 8911db7f authored by WAUTELET Philippe's avatar WAUTELET Philippe
Browse files

Philippe 31/08/2018: OpenACC: condensation and compute_frac_ice3d_device

parent 3cc5ca8a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/MNH/compute_frac_ice.f90
+ 95
− 2
View file @ 8911db7f
!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC Copyright 2006-2018 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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
! ######spl
MODULE MODI_COMPUTE_FRAC_ICE
......@@ -34,6 +34,17 @@ REAL, DIMENSION(:), INTENT(INOUT) :: PFRAC_ICE
END INTERFACE
!
!
CONTAINS
!
SUBROUTINE COMPUTE_FRAC_ICE3D_DEVICE(HFRAC_ICE,PFRAC_ICE,PT)
!
CHARACTER*1 , INTENT(IN) :: HFRAC_ICE
REAL, DIMENSION(:,:,:), INTENT(IN) :: PT
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFRAC_ICE
!
END SUBROUTINE COMPUTE_FRAC_ICE3D_DEVICE
!
END MODULE MODI_COMPUTE_FRAC_ICE
!
! ##############################
......@@ -284,3 +295,85 @@ PFRAC_ICE(:) = MAX( 0., MIN(1., PFRAC_ICE(:) ) )
END SUBROUTINE COMPUTE_FRAC_ICE1D
#ifdef _OPENACC
! #############################################################
SUBROUTINE COMPUTE_FRAC_ICE3D_DEVICE(HFRAC_ICE,PFRAC_ICE,PT)
! #############################################################
!
!
!!**** *COMPUTE_FRAC_ICE* - computes ice fraction
!!
!! PURPOSE
!! -------
!!
!!** METHOD
!! ------
!!
!!
!! EXTERNAL
!! --------
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!!
!!
!! REFERENCE
!! ---------
!!
!! AUTHOR
!! ------
!! Julien PERGAUD * Meteo-France *
!!
!! MODIFICATIONS
!! -------------
!! Original 13/03/06
!! S. Riette April 2011 optimisation
!! S. Riette 08/2016 add option O
!!
!! --------------------------------------------------------------------------
! 0. DECLARATIONS
! ------------
!
USE MODD_NEB, ONLY : XTMINMIX, XTMAXMIX
USE MODD_CST, ONLY : XTT
USE MODE_MSG
!
IMPLICIT NONE
!
!
!* 0.1 declarations of arguments
!
CHARACTER(LEN=1), INTENT(IN) :: HFRAC_ICE ! scheme to use
REAL, DIMENSION(:,:,:), INTENT(IN) :: PT ! temperature
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFRAC_ICE ! Ice fraction (1 for ice only, 0 for liquid only)
!$acc declare present(PT,PFRAC_ICE)
!
! 0.2 declaration of local variables
!
!
! 0.2 initialisation
!
!
!------------------------------------------------------------------------
! 1. Compute FRAC_ICE
!
!$acc kernels present(PT,PFRAC_ICE)
SELECT CASE(HFRAC_ICE)
CASE ('T') !using Temperature
PFRAC_ICE(:,:,:) = ( XTMAXMIX - PT(:,:,:) ) / ( XTMAXMIX - XTMINMIX ) ! freezing interval
CASE ('O') !using Temperature with old formulae
PFRAC_ICE(:,:,:) = ( XTT - PT(:,:,:) ) / 40. ! freezing interval
CASE ('N') !No ice
PFRAC_ICE(:,:,:) = 0.
CASE ('S') !Same as previous
!nothing to do
CASE DEFAULT
CALL PRINT_MSG(NVERB_FATAL,'GEN','COMPUTE_FRAC_ICE3D_DEVICE','invalid option')
END SELECT
!
PFRAC_ICE(:,:,:) = MAX( 0., MIN(1., PFRAC_ICE(:,:,:) ) )
!$acc end kernels
!
END SUBROUTINE COMPUTE_FRAC_ICE3D_DEVICE
#endif
src/MNH/condensation.f90
+ 62
− 31
View file @ 8911db7f
......@@ -182,17 +182,21 @@ REAL, DIMENSION(KIU,KJU,KKU) :: ZLV, ZLS, ZCPD
REAL :: ZTEMP, ZPV, ZQSL, ZPIV, ZQSI, ZCOND, ZLVS ! thermodynamics
REAL :: ZLL, DZZ, ZZZ ! used for length scales
REAL :: ZAH, ZA, ZB, ZSBAR, ZQ1, ZSIGMA, ZDRW, ZDTL, ZSIG_CONV ! related to computation of Sig_s
INTEGER :: INQ1
REAL :: ZINC
LOGICAL :: GPRESENT_PLV, GPRESENT_PLS, GPRESENT_PCPH
LOGICAL, DIMENSION(KIU,KJU,KKU) :: GWORK
!
!* 0.3 Definition of constants :
!
!-------------------------------------------------------------------------------
!
REAL :: ZL0 = 600. ! tropospheric length scale
REAL :: ZCSIGMA = 0.2 ! constant in sigma_s parameterization
REAL :: ZCSIG_CONV = 0.30E-2 ! scaling factor for ZSIG_CONV as function of mass flux
REAL,PARAMETER :: ZL0 = 600. ! tropospheric length scale
REAL,PARAMETER :: ZCSIGMA = 0.2 ! constant in sigma_s parameterization
REAL,PARAMETER :: ZCSIG_CONV = 0.30E-2 ! scaling factor for ZSIG_CONV as function of mass flux
!
REAL, DIMENSION(-22:11) :: ZSRC_1D =(/ &
REAL, DIMENSION(-22:11),PARAMETER :: ZSRC_1D =(/ &
0. , 0. , 2.0094444E-04, 0.316670E-03, &
4.9965648E-04, 0.785956E-03 , 1.2341294E-03, 0.193327E-02, &
3.0190963E-03, 0.470144E-02 , 7.2950651E-03, 0.112759E-01, &
......@@ -202,8 +206,8 @@ REAL, DIMENSION(-22:11) :: ZSRC_1D =(/ &
0.8413813 , 0.933222E+00 , 0.9772662 , 0.993797E+00, &
0.9986521 , 0.999768E+00 , 0.9999684 , 0.999997E+00, &
1.0000000 , 1.000000 /)
INTEGER :: INQ1
REAL :: ZINC
!
!$acc declare create(ZFRAC) device_resident(ZTLK,ZRT,ZL,ITPL,ZTMIN,ZLV,ZLS,ZCPD,GWORK)
!
!-------------------------------------------------------------------------------
!
......@@ -211,28 +215,18 @@ REAL :: ZINC
IKTB=1+JPVEXT
IKTE=KKU-JPVEXT
!$acc update self(PZZ,PRS,PRG)
!PW:TODO: the following update can be removed once CONDENSATION is fully openacc
!PT,PSIGS,PLV,PLS,PCPH,PRV,PRC,PRIare OK on device in entry of this SUBROUTINE
!$acc update self(PT,PSIGS,PRV,PRC,PRI)
IF(PRESENT(PLV)) THEN
!$acc update self(PLV)
END IF
IF(PRESENT(PLS)) THEN
!$acc update self(PLS)
END IF
IF(PRESENT(PCPH)) THEN
!$acc update self(PCPH)
END IF
!
GPRESENT_PLV = .FALSE.
GPRESENT_PLS = .FALSE.
GPRESENT_PCPH = .FALSE.
IF (PRESENT(PLV)) GPRESENT_PLV = .TRUE.
IF (PRESENT(PLS)) GPRESENT_PLS = .TRUE.
IF (PRESENT(PCPH)) GPRESENT_PCPH = .TRUE.
!
!$acc kernels
PCLDFR(:,:,:) = 0. ! Initialize values
!$acc end kernels
!$acc update self(PCLDFR)
PSIGRC(:,:,:) = 0. ! Initialize values
!
!
!-------------------------------------------------------------------------------
! store total water mixing ratio
DO JK=IKTB,IKTE
......@@ -245,7 +239,7 @@ END DO
!-------------------------------------------------------------------------------
! Preliminary calculations
! latent heat of vaporisation/sublimation
IF(PRESENT(PLV) .AND. PRESENT(PLS)) THEN
IF(GPRESENT_PLV .AND. GPRESENT_PLS) THEN
ZLV(:,:,:)=PLV(:,:,:)
ZLS(:,:,:)=PLS(:,:,:)
ELSE
......@@ -260,7 +254,7 @@ ELSE
ENDDO
ENDDO
ENDIF
IF(PRESENT(PCPH)) THEN
IF(GPRESENT_PCPH) THEN
ZCPD(:,:,:)=PCPH(:,:,:)
ELSE
DO JK=IKTB,IKTE
......@@ -321,11 +315,20 @@ END IF
!
!Ice fraction
ZFRAC(:,:,:) = 0.
WHERE(PRC(:,:,:)+PRI(:,:,:) > 1.E-20)
ZFRAC(:,:,:) = PRI(:,:,:) / (PRC(:,:,:)+PRI(:,:,:))
ENDWHERE
CALL COMPUTE_FRAC_ICE(HFRAC_ICE, ZFRAC, PT)
IF(.NOT. OUSERI) ZFRAC(:,:,:)=0.
IF (OUSERI) THEN
GWORK(:,:,:) = PRC(:,:,:)+PRI(:,:,:) > 1.E-20
WHERE( GWORK(:,:,:) )
ZFRAC(:,:,:) = PRI(:,:,:) / (PRC(:,:,:)+PRI(:,:,:))
ENDWHERE
ENDIF
!$acc end kernels
!
#ifndef _OPENACC
IF (OUSERI) CALL COMPUTE_FRAC_ICE(HFRAC_ICE, ZFRAC, PT)
#else
IF (OUSERI) CALL COMPUTE_FRAC_ICE3D_DEVICE(HFRAC_ICE, ZFRAC, PT)
#endif
!$acc kernels
!
DO JK=IKTB,IKTE
JKP=MAX(MIN(JK+KKL,IKTE),IKTB)
......@@ -335,11 +338,19 @@ DO JK=IKTB,IKTE
! latent heats
ZTEMP = PT(JI,JJ,JK)
! saturated water vapor mixing ratio over liquid water
#ifndef MNH_BITREP
ZPV = MIN(EXP( XALPW - XBETAW / ZTEMP - XGAMW * LOG( ZTEMP ) ), .99*PPABS(JI,JJ,JK))
#else
ZPV = MIN(BR_EXP( XALPW - XBETAW / ZTEMP - XGAMW * BR_LOG( ZTEMP ) ), .99*PPABS(JI,JJ,JK))
#endif
ZQSL = XRD / XRV * ZPV / ( PPABS(JI,JJ,JK) - ZPV )
! saturated water vapor mixing ratio over ice
#ifndef MNH_BITREP
ZPIV = MIN(EXP( XALPI - XBETAI / ZTEMP - XGAMI * LOG( ZTEMP ) ), .99*PPABS(JI,JJ,JK))
#else
ZPIV = MIN(BR_EXP( XALPI - XBETAI / ZTEMP - XGAMI * BR_LOG( ZTEMP ) ), .99*PPABS(JI,JJ,JK))
#endif
ZQSI = XRD / XRV * ZPIV / ( PPABS(JI,JJ,JK) - ZPIV )
! interpolate between liquid and solid as function of temperature
......@@ -348,7 +359,11 @@ DO JK=IKTB,IKTE
& ZFRAC(JI,JJ,JK) * ZLS(JI,JJ,JK)
! coefficients a and b
#ifndef MNH_BITREP
ZAH = ZLVS * ZQSL / ( XRV * ZTEMP**2 ) * (XRV * ZQSL / XRD + 1.)
#else
ZAH = ZLVS * ZQSL / ( XRV * BR_P2(ZTEMP) ) * (XRV * ZQSL / XRD + 1.)
#endif
ZA = 1. / ( 1. + ZLVS/ZCPD(JI,JJ,JK) * ZAH )
ZB = ZAH * ZA
......@@ -359,7 +374,11 @@ DO JK=IKTB,IKTE
! or that of Meso-NH turbulence scheme
IF ( OSIGMAS ) THEN
IF (PSIGQSAT/=0.) THEN
#ifndef MNH_BITREP
ZSIGMA = SQRT((2*PSIGS(JI,JJ,JK))**2 + (PSIGQSAT*ZQSL*ZA)**2)
#else
ZSIGMA = BR_POW(BR_P2(2*PSIGS(JI,JJ,JK)) + BR_P2(PSIGQSAT*ZQSL*ZA) , 0.5)
#endif
ELSE
ZSIGMA = 2*PSIGS(JI,JJ,JK)
END IF
......@@ -385,15 +404,27 @@ DO JK=IKTB,IKTE
ZQ1 = ZSBAR/ZSIGMA
! cloud fraction
#ifndef MNH_BITREP
PCLDFR(JI,JJ,JK) = MAX( 0., MIN(1.,0.5+0.36*ATAN(1.55*ZQ1)) )
#else
PCLDFR(JI,JJ,JK) = MAX( 0., MIN(1.,0.5+0.36*BR_ATAN(1.55*ZQ1)) )
#endif
! total condensate
IF (ZQ1 > 0. .AND. ZQ1 <= 2 ) THEN
#ifndef MNH_BITREP
ZCOND = MIN(EXP(-1.)+.66*ZQ1+.086*ZQ1**2, 2.) ! We use the MIN function for continuity
#else
ZCOND = MIN(BR_EXP(-1.)+.66*ZQ1+.086*BR_P2(ZQ1), 2.) ! We use the MIN function for continuity
#endif
ELSE IF (ZQ1 > 2.) THEN
ZCOND = ZQ1
ELSE
#ifndef MNH_BITREP
ZCOND = EXP( 1.2*ZQ1-1. )
#else
ZCOND = BR_EXP( 1.2*ZQ1-1. )
#endif
ENDIF
ZCOND = ZCOND * ZSIGMA
......@@ -432,7 +463,7 @@ DO JK=IKTB,IKTE
END DO
END DO
END DO
!$acc end kernels
!
!$acc update device(PCLDFR,PT,PRV,PRC,PRI,PSIGRC)
!
END SUBROUTINE CONDENSATION
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

Documentation | Mentions légales | CGU | Accessibilité : non conforme