diff --git a/src/MNH/conjgrad.f90 b/src/MNH/conjgrad.f90 index f80d9042eb4db6acafd3b1e15a98c7265bc6a682..4db68c47a5b83636ae2d08755db811caf50562ad 100644 --- a/src/MNH/conjgrad.f90 +++ b/src/MNH/conjgrad.f90 @@ -117,7 +117,7 @@ END MODULE MODI_CONJGRAD !! !! AUTHOR !! ------ -!! P. HÅreil and J. Stein * Meteo France * +!! P. Hareil and J. Stein * Meteo France * !! !! MODIFICATIONS !! ------------- @@ -290,11 +290,7 @@ DO JM = 1,KITR ! -1 ! compute the vector DELTA = F * ( Y - Q ( PHI ) ) ! -#ifndef MNH_OPENACC - ZWORK(:,:,:) = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) ! Q (PHI) -#else - CALL QLAP_DEVICE( ZWORK, HLBCX, HLBCY, PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, PPHI ) ! Q (PHI) -#endif + CALL QLAP( ZWORK, HLBCX, HLBCY, PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, PPHI ) ! Q (PHI) ! !$acc kernels ZWORK(:,:,:) = PY(:,:,:) - ZWORK(:,:,:) ! Y - Q (PHI) @@ -313,11 +309,7 @@ DO JM = 1,KITR ZP(:,:,:) = ZDELTA(:,:,:) ! P = DELTA at the first solver iteration !$acc end kernels ELSE -#ifndef MNH_OPENACC - ZWORK(:,:,:) = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZDELTA) ! Q ( DELTA ) -#else - CALL QLAP_DEVICE( ZWORK, HLBCX, HLBCY, PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, ZDELTA ) ! Q ( DELTA ) -#endif + CALL QLAP( ZWORK, HLBCX, HLBCY, PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, ZDELTA ) ! Q ( DELTA ) CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! -1 PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZWORKD) ! F * Q ( DELTA ) ! @@ -331,11 +323,7 @@ DO JM = 1,KITR !* 2.4 compute LAMBDA ! ! -#ifndef MNH_OPENACC - ZWORK(:,:,:) = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) ! Q ( P ) -#else - CALL QLAP_DEVICE( ZWORK, HLBCX, HLBCY, PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, ZP ) ! Q ( P ) -#endif + CALL QLAP( ZWORK, HLBCX, HLBCY, PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, ZP ) ! Q ( P ) CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,& ! -1 PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZWORKP) ! F * Q ( P ) ! diff --git a/src/MNH/conresol.f90 b/src/MNH/conresol.f90 index d7493a23e800f84cb8b5613a24dab5683d55c33e..51bb0bb575b1c383a27e4ad0a4df895d78be071c 100644 --- a/src/MNH/conresol.f90 +++ b/src/MNH/conresol.f90 @@ -258,14 +258,10 @@ CALL MNH_MEM_GET( ZRESIDUE, JIU, JJU, JKU ) ! !* 1.1 compute the vector: r^(0) = Q(PHI) - Y ! -#ifndef MNH_OPENACC -ZRESIDUE(:,:,:) = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) - PY(:,:,:) -#else -CALL QLAP_DEVICE(ZRESIDUE,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) +CALL QLAP(ZRESIDUE,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) !$acc kernels present( PY, ZRESIDUE ) ZRESIDUE(:,:,:) = ZRESIDUE(:,:,:) - PY(:,:,:) !$acc end kernels -#endif ! !* 1.2 compute the vector: p^(0) = F^(-1)*( Q(PHI) - Y ) ! @@ -274,11 +270,7 @@ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! !* 1.3 compute the vector: delta^(0) = Q ( p^(0) ) ! -#ifndef MNH_OPENACC -ZDELTA = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) -#else -CALL QLAP_DEVICE(ZDELTA,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) -#endif +CALL QLAP(ZDELTA,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) ! !------------------------------------------------------------------------------- ! @@ -315,11 +307,7 @@ DO JM = 1,KITR ! !* 2.5 compute the auxiliary field: ksi = Q ( q ) ! -#ifndef MNH_OPENACC - ZKSI= QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) -#else - CALL QLAP_DEVICE(ZKSI,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) -#endif + CALL QLAP(ZKSI,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) ! -1 !* 2.6 compute the step ALPHA ! diff --git a/src/MNH/conresolz.f90 b/src/MNH/conresolz.f90 index fb1dc51bdd99c4550ac05f760a3b9ef8ebcd11df..7e7f5ac60070e6b7363981857913b4e7046b31d0 100644 --- a/src/MNH/conresolz.f90 +++ b/src/MNH/conresolz.f90 @@ -268,14 +268,10 @@ CALL MNH_MEM_GET( ZRESIDUE , JIU, JJU, JKU ) ! !* 1.1 compute the vector: r^(0) = Q(PHI) - Y ! -#ifndef MNH_OPENACC -ZRESIDUE(:,:,:) = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) - PY(:,:,:) -#else -CALL QLAP_DEVICE(ZRESIDUE,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) +CALL QLAP(ZRESIDUE,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) !$acc kernels present( PY, ZRESIDUE ) ZRESIDUE(:,:,:) = ZRESIDUE(:,:,:) - PY(:,:,:) !$acc end kernels -#endif ! !* 1.2 compute the vector: p^(0) = F^(-1)*( Q(PHI) - Y ) ! @@ -304,11 +300,7 @@ END IF ! !* 1.3 compute the vector: delta^(0) = Q ( p^(0) ) ! -#ifndef MNH_OPENACC -ZDELTA = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) -#else -CALL QLAP_DEVICE(ZDELTA,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) -#endif +CALL QLAP(ZDELTA,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) ! !------------------------------------------------------------------------------- ! @@ -347,11 +339,7 @@ CALL FLAT_INVZ(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! !* 2.5 compute the auxiliary field: ksi = Q ( q ) ! -#ifndef MNH_OPENACC - ZKSI= QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) -#else - CALL QLAP_DEVICE(ZKSI,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) -#endif + CALL QLAP(ZKSI,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) ! -1 !* 2.6 compute the step ALPHA ! diff --git a/src/MNH/qlap.f90 b/src/MNH/qlap.f90 index f7f68d8db60045b3699d1a5bb817c0b6348ef359..5e6d15843ec5b72395cf12f929940871597ecec7 100644 --- a/src/MNH/qlap.f90 +++ b/src/MNH/qlap.f90 @@ -9,37 +9,11 @@ ! INTERFACE ! - FUNCTION QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PY) & - RESULT(PQLAP) -! -IMPLICIT NONE -! -CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type -CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type -! - ! Metric coefficients: -REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx -REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDZX ! d*zx -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDZY ! d*zy -REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz -! -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J -REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at time t -! -REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! field components -! -REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PQLAP ! final divergence -! -END FUNCTION QLAP -! -! -#ifdef MNH_OPENACC -SUBROUTINE QLAP_DEVICE(PQLAP,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PY) +SUBROUTINE QLAP(PQLAP,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PY) ! IMPLICIT NONE ! -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQLAP ! final divergence +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PQLAP ! final divergence ! CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type @@ -58,8 +32,8 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! field components ! ! -END SUBROUTINE QLAP_DEVICE -#endif +END SUBROUTINE QLAP + END INTERFACE ! END MODULE MODI_QLAP @@ -67,272 +41,7 @@ END MODULE MODI_QLAP ! ! ! ######################################################################### - FUNCTION QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PY) & - RESULT(PQLAP) -! ######################################################################### -! -!!**** *QLAP * - compute the complete quasi-laplacien QLAP of a field P -!! -!! PURPOSE -!! ------- -! This function computes the quasi-laplacien QLAP of the scalar field P -! localized at a mass point, with non-vanishing orography. -! The result is localized at a mass point and defined by: -! for Durran and MAE anelastic equations -! ( ( GX_M_U (PY) ) ) -! PQLAP = GDIV ( rho * CPd * Thetav * J ( GX_M_V (PY) ) ) -! ( ( GX_M_W (PY) ) ) -! or for Lipps and Hemler -! ( ( GX_M_U (PY) ) ) -! PQLAP = GDIV ( rho * J ( GX_M_V (PY) ) ) -! ( ( GX_M_W (PY) ) ) -! Where GX_M_.. are the cartesian components of the gradient of PY and -! GDIV is the operator acting on a vector AA: -! -! GDIV ( AA ) = J * divergence (1/J AA ) -! -!!** METHOD -!! ------ -!! First, we compute the gradients along x, y , z of the P field. The -!! result is multiplied by rhod * CPd * Thetav or rhod depending on the -!! chosen anelastic system where the suffixes indicate -!! d dry and v for virtual. -!! Then, the pseudo-divergence ( J * DIV (1/J o ) ) is computed by the -!! subroutine GDIV. The result is localized at a mass point. -!! -!! EXTERNAL -!! -------- -!! Function GX_M_U : compute the gradient along x -!! Function GY_M_V : compute the gradient along y -!! Function GZ_M_W : compute the gradient along z -!! FUNCTION MXM: compute an average in the x direction for a variable -!! at a mass localization -!! FUNCTION MYM: compute an average in the y direction for a variable -!! at a mass localization -!! FUNCTION MZM: compute an average in the z direction for a variable -!! at a mass localization -!! Subroutine GDIV : compute J times the divergence of 1/J times a vector -!! -!! IMPLICIT ARGUMENTS -!! ------------------ -!! Module MODD_PARAMETERS: JPHEXT, JPVEXT -!! Module MODD_CONF: L2D,CEQNSYS -!! Module MODD_CST : XCPD -!! -!! REFERENCE -!! --------- -!! Pressure solver documentation ( Scientific documentation ) -!! -!! AUTHOR -!! ------ -!! P. Hereil and J. Stein * Meteo France * -!! -!! MODIFICATIONS -!! ------------- -!! Original 11/07/94 -!! Modification 16/03/95 change the argument list of the gradient -!! 14/01/97 New anelastic equation ( Stein ) -!! 17/12/97 include the case of non-vanishing orography -!! at the lbc ( Stein ) -!! 06/12 V.Masson : update_halo due to CONTRAV changes -!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 -! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine -! F. Auguste 02/21: add IBM -!------------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -USE MODE_ll -! -USE MODD_PARAMETERS -USE MODD_CONF -USE MODD_CST -USE MODI_GDIV -USE MODI_GRADIENT_M -USE MODI_SHUMAN -! -USE MODE_MPPDB -! -USE MODD_IBM_PARAM_n, ONLY: XIBM_LS, LIBM, XIBM_SU -USE MODI_IBM_BALANCE -! -IMPLICIT NONE -! -!* 0.1 declarations of arguments -! -! -CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type -CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type -! - ! Metric coefficients: -REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX ! d*xx -REAL, DIMENSION(:,:,:), INTENT(IN) :: PDYY ! d*yy -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDZX ! d*zx -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDZY ! d*zy -REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ ! d*zz -! -REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference * J -REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHETAV ! virtual potential temp. at time t -! -REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! field components -! -REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: PQLAP ! final divergence -! -!* 0.2 declarations of local variables -! -REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: ZU ! rho*J*gradient along x -! -REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: ZV ! rho*J*gradient along y -! -REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: ZW ! rho*J*gradient along z -! -INTEGER :: IIU,IJU,IKU ! I,J,K array sizes -INTEGER :: JK,JJ,JI ! vertical loop index -TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange -INTEGER :: IINFO_ll -INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE -!------------------------------------------------------------------------------- -IF (MPPDB_INITIALIZED) THEN - !Check all IN arrays - CALL MPPDB_CHECK(PDXX,"QLAP beg:PDXX") - CALL MPPDB_CHECK(PDYY,"QLAP beg:PDYY") - CALL MPPDB_CHECK(PDZX,"QLAP beg:PDZX") - CALL MPPDB_CHECK(PDZY,"QLAP beg:PDZY") - CALL MPPDB_CHECK(PDZZ,"QLAP beg:PDZZ") - CALL MPPDB_CHECK(PRHODJ,"QLAP beg:PRHODJ") - CALL MPPDB_CHECK(PTHETAV,"QLAP beg:PTHETAV") - CALL MPPDB_CHECK(PY,"QLAP beg:PY") -END IF -! -! -!* 1. COMPUTE THE GRADIENT COMPONENTS -! ------------------------------- -! -CALL GET_DIM_EXT_ll('B',IIU,IJU) -CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) -IKU=SIZE(PY,3) -IKE = IKU - JPVEXT -IKB = 1 + JPVEXT -! -ZU(:,:,:) = GX_M_U( 1, IKU, 1, PY(:,:,:), PDXX(:,:,:), PDZZ(:,:,:), PDZX(:,:,:) ) -CALL MPPDB_CHECK3D(ZU,'QLAP::ZU',PRECISION) -! -IF ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() ) THEN - DO JK=2,IKU-1 - DO JJ=1,IJU - ZU(IIB,JJ,JK)= (PY(IIB,JJ,JK) - PY(IIB-1,JJ,JK) - 0.5 * ( & - PDZX(IIB,JJ,JK) * (PY(IIB,JJ,JK)-PY(IIB,JJ,JK-1)) / PDZZ(IIB,JJ,JK) & - +PDZX(IIB,JJ,JK+1) * (PY(IIB,JJ,JK+1)-PY(IIB,JJ,JK)) / PDZZ(IIB,JJ,JK+1) & - ) ) / PDXX(IIB,JJ,JK) - END DO - END DO -END IF -CALL MPPDB_CHECK3D(ZU,'QLAP::ZU/W',PRECISION) -! -IF ( HLBCX(1) /= 'CYCL' .AND. LEAST_ll() ) THEN - DO JK=2,IKU-1 - DO JJ=1,IJU - ZU(IIE+1,JJ,JK)= (PY(IIE+1,JJ,JK) - PY(IIE+1-1,JJ,JK) - 0.5 * ( & - PDZX(IIE+1,JJ,JK) * (PY(IIE+1-1,JJ,JK)-PY(IIE+1-1,JJ,JK-1)) / PDZZ(IIE+1-1,JJ,JK) & - +PDZX(IIE+1,JJ,JK+1) * (PY(IIE+1-1,JJ,JK+1)-PY(IIE+1-1,JJ,JK)) / PDZZ(IIE+1-1,JJ,JK+1)& - ) ) / PDXX(IIE+1,JJ,JK) - END DO - END DO -END IF -CALL MPPDB_CHECK3D(ZU,'QLAP::ZU/E',PRECISION) -! -IF(.NOT. L2D) THEN -! - ZV(:,:,:) = GY_M_V( 1, IKU, 1, PY(:,:,:), PDYY(:,:,:), PDZZ(:,:,:), PDZY(:,:,:) ) - CALL MPPDB_CHECK3D(ZV,'QLAP::ZV',PRECISION) -! - IF ( HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() ) THEN - DO JK=2,IKU-1 - DO JI=1,IIU - ZV(JI,IJB,JK)= (PY(JI,IJB,JK) - PY(JI,IJB-1,JK) - 0.5 * ( & - PDZY(JI,IJB,JK) * (PY(JI,IJB,JK)-PY(JI,IJB,JK-1)) / PDZZ(JI,IJB,JK) & - +PDZY(JI,IJB,JK+1) * (PY(JI,IJB,JK+1)-PY(JI,IJB,JK)) / PDZZ(JI,IJB,JK+1) & - ) ) / PDYY(JI,IJB,JK) - END DO - END DO - END IF - CALL MPPDB_CHECK3D(ZV,'QLAP::ZV/S',PRECISION) - IF ( HLBCY(1) /= 'CYCL' .AND. LNORTH_ll() ) THEN -! - DO JK=2,IKU-1 - DO JI=1,IIU - ZV(JI,IJE+1,JK)= (PY(JI,IJE+1,JK) - PY(JI,IJE+1-1,JK) - 0.5 * ( & - PDZY(JI,IJE+1,JK) * (PY(JI,IJE+1-1,JK)-PY(JI,IJE+1-1,JK-1)) / PDZZ(JI,IJE+1-1,JK) & - +PDZY(JI,IJE+1,JK+1) * (PY(JI,IJE+1-1,JK+1)-PY(JI,IJE+1-1,JK)) / PDZZ(JI,IJE+1-1,JK+1)& - ) ) / PDYY(JI,IJE+1,JK) - END DO - END DO - END IF - CALL MPPDB_CHECK3D(ZV,'QLAP::ZV/N',PRECISION) -! -ELSE - ZV(:,:,:) = 0. -ENDIF -! -IF ( CEQNSYS == 'DUR' .OR. CEQNSYS == 'MAE' ) THEN - ZU(:,:,:) = MXM( PRHODJ(:,:,:) * XCPD * PTHETAV(:,:,:) ) * ZU(:,:,:) - IF(.NOT. L2D) THEN - ZV(:,:,:) = MYM( PRHODJ(:,:,:) * XCPD * PTHETAV(:,:,:) ) * ZV(:,:,:) - END IF - ZW(:,:,:) = MZM( PRHODJ(:,:,:) * XCPD * PTHETAV(:,:,:) ) * GZ_M_W( 1, IKU, 1, PY(:,:,:), PDZZ(:,:,:) ) -ELSEIF ( CEQNSYS == 'LHE' ) THEN - ZU(:,:,:) = MXM( PRHODJ(:,:,:) ) * ZU(:,:,:) - IF(.NOT. L2D) THEN - ZV(:,:,:) = MYM( PRHODJ(:,:,:) ) * ZV(:,:,:) - ENDIF - ZW(:,:,:) = MZM( PRHODJ(:,:,:) ) * GZ_M_W( 1, IKU, 1, PY(:,:,:), PDZZ(:,:,:) ) -END IF -! -!------------------------------------------------------------------------------- -! -!* 2. COMPUTE THE DIVERGENCE -! ---------------------- -! -NULLIFY(TZFIELDS_ll) -CALL ADD3DFIELD_ll( TZFIELDS_ll, ZU, 'QLAP::ZU' ) -CALL ADD3DFIELD_ll( TZFIELDS_ll, ZV, 'QLAP::ZV' ) -CALL ADD3DFIELD_ll( TZFIELDS_ll, ZW, 'QLAP::ZW' ) -CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll) -CALL CLEANLIST_ll(TZFIELDS_ll) -! -CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PQLAP) -! -IF (LIBM) THEN - ! - CALL IBM_BALANCE(XIBM_LS,XIBM_SU,ZU,ZV,ZW,PQLAP) - ! - PQLAP(:,:,IKB-1) = PQLAP(:,:,IKB-1)*XIBM_SU(:,:,IKB,1) - PQLAP(:,:,IKE+1) = PQLAP(:,:,IKE+1)*XIBM_SU(:,:,IKE,1) - ! - IF ( HLBCX(1) /= 'CYCL' ) THEN - IF(LWEST_ll()) PQLAP(IIB-1,:,:) = PQLAP(IIB-1,:,:)*XIBM_SU(IIB,:,:,1) - IF(LEAST_ll()) PQLAP(IIE+1,:,:) = PQLAP(IIE+1,:,:)*XIBM_SU(IIE,:,:,1) - ENDIF - ! - IF ( HLBCY(1) /= 'CYCL' ) THEN - IF (LSOUTH_ll()) PQLAP(:,IJB-1,:) = PQLAP(:,IJB-1,:)*XIBM_SU(:,IJB,:,1) - IF (LNORTH_ll()) PQLAP(:,IJE+1,:) = PQLAP(:,IJE+1,:)*XIBM_SU(:,IJE,:,1) - ENDIF - ! -ENDIF -! -IF (MPPDB_INITIALIZED) THEN - !Check all OUT arrays - CALL MPPDB_CHECK(PQLAP,"QLAP end:PQLAP") -END IF -!------------------------------------------------------------------------------- -! -END FUNCTION QLAP - -#ifdef MNH_OPENACC -! ######################################################################### - SUBROUTINE QLAP_DEVICE(PQLAP,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PY) + SUBROUTINE QLAP(PQLAP,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PY) ! ######################################################################### ! !!**** *QLAP * - compute the complete quasi-laplacien QLAP of a field P @@ -414,12 +123,16 @@ USE MODD_CONF USE MODD_CST USE MODI_GDIV USE MODI_GRADIENT_M +#ifndef MNH_OPENACC USE MODI_SHUMAN +#else USE MODI_SHUMAN_DEVICE +#endif ! USE MODE_MPPDB -! +#ifdef MNH_OPENACC USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE +#endif ! USE MODI_GET_HALO USE MODD_IBM_PARAM_n, ONLY: XIBM_LS, LIBM, XIBM_SU @@ -430,7 +143,7 @@ IMPLICIT NONE !* 0.1 declarations of arguments ! ! -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PQLAP ! final divergence +REAL, DIMENSION(:,:,:), INTENT(OUT) :: PQLAP ! final divergence ! CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! x-direction LBC type CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! y-direction LBC type @@ -449,21 +162,29 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! field components ! !* 0.2 declarations of local variables ! -REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZU ! rho*J*gradient along x -! -REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZV ! rho*J*gradient along y -! -REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZW ! rho*J*gradient along z +#ifndef MNH_OPENACC +REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZU ! rho*J*gradient along x +REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZV ! rho*J*gradient along y +REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZW ! rho*J*gradient along z +#else +REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZU ! rho*J*gradient along x +REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZV ! rho*J*gradient along y +REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZW ! rho*J*gradient along z +#endif ! INTEGER :: IIU,IJU,IKU ! I,J,K array sizes INTEGER :: JK,JJ,JI ! vertical loop index +#ifndef MNH_OPENACC TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange INTEGER :: IINFO_ll +#endif INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE ! -REAL, DIMENSION(:,:,:), pointer , contiguous :: ZMXM,ZMYM,ZMZM,ZRHODJ,ZGZMW +#ifdef MNH_OPENACC +REAL, DIMENSION(:,:,:), pointer, contiguous :: ZMXM, ZMYM, ZMZM, ZRHODJ, ZGZMW +#endif ! -LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH +LOGICAL :: GWEST, GEAST, GSOUTH, GNORTH !------------------------------------------------------------------------------- IF (MPPDB_INITIALIZED) THEN !Check all IN arrays @@ -492,7 +213,12 @@ GEAST = ( HLBCX(2) /= 'CYCL' .AND. LEAST_ll() ) GSOUTH = ( HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() ) GNORTH = ( HLBCY(2) /= 'CYCL' .AND. LNORTH_ll() ) ! -CALL MNH_MEM_POSITION_PIN() +#ifndef MNH_OPENACC +ALLOCATE( ZU(IIU, IJU, IKU) ) +ALLOCATE( ZV(IIU, IJU, IKU) ) +ALLOCATE( ZW(IIU, IJU, IKU) ) +#else +CALL MNH_MEM_POSITION_PIN( 'QLAP' ) CALL MNH_MEM_GET( ZU, IIU, IJU, IKU ) CALL MNH_MEM_GET( ZV, IIU, IJU, IKU ) CALL MNH_MEM_GET( ZW, IIU, IJU, IKU ) @@ -502,12 +228,21 @@ CALL MNH_MEM_GET( ZMYM, IIU, IJU, IKU ) CALL MNH_MEM_GET( ZMZM, IIU, IJU, IKU ) CALL MNH_MEM_GET( ZRHODJ, IIU, IJU, IKU ) CALL MNH_MEM_GET( ZGZMW, IIU, IJU, IKU ) +#endif + +!$acc data present( PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, PY, & +!$acc & ZU, ZV, ZW, ZMXM, ZMYM, ZMZM, ZRHODJ, ZGZMW ) + ! +#ifndef MNH_OPENACC +ZU(:,:,:) = GX_M_U( 1, IKU, 1, PY(:,:,:), PDXX(:,:,:), PDZZ(:,:,:), PDZX(:,:,:) ) +#else CALL GX_M_U_DEVICE( 1, IKU, 1, PY(:,:,:), PDXX(:,:,:), PDZZ(:,:,:), PDZX(:,:,:), ZU(:,:,:) ) +#endif CALL MPPDB_CHECK3D(ZU,'QLAP::ZU',PRECISION) ! IF ( GWEST ) THEN - !$acc kernels async present_cr(ZU) + !$acc kernels async DO JK=2,IKU-1 DO JJ=1,IJU ZU(IIB,JJ,JK)= (PY(IIB,JJ,JK) - PY(IIB-1,JJ,JK) - 0.5 * ( & @@ -520,7 +255,7 @@ IF ( GWEST ) THEN END IF ! IF ( GEAST ) THEN - !$acc kernels async present_cr(ZU) + !$acc kernels async DO JK=2,IKU-1 DO JJ=1,IJU ZU(IIE+1,JJ,JK)= (PY(IIE+1,JJ,JK) - PY(IIE+1-1,JJ,JK) - 0.5 * ( & @@ -533,16 +268,20 @@ IF ( GEAST ) THEN END IF !$acc wait ! -CALL MPPDB_CHECK3D(ZU,'QLAP::ZU/W',PRECISION) -CALL MPPDB_CHECK3D(ZU,'QLAP::ZU/E',PRECISION) +!CALL MPPDB_CHECK3D(ZU,'QLAP::ZU/W',PRECISION) +!CALL MPPDB_CHECK3D(ZU,'QLAP::ZU/E',PRECISION) ! IF(.NOT. L2D) THEN ! - CALL GY_M_V_DEVICE( 1, IKU, 1, PY(:,:,:), PDYY(:,:,:), PDZZ(:,:,:), PDZY(:,:,:), ZV(:,:,:) ) - CALL MPPDB_CHECK3D(ZV,'QLAP::ZV',PRECISION) +#ifndef MNH_OPENACC + ZV(:,:,:) = GY_M_V( 1, IKU, 1, PY(:,:,:), PDYY(:,:,:), PDZZ(:,:,:), PDZY(:,:,:) ) +#else + CALL GY_M_V_DEVICE( 1, IKU, 1, PY(:,:,:), PDYY(:,:,:), PDZZ(:,:,:), PDZY(:,:,:), ZV(:,:,:) ) +#endif +! CALL MPPDB_CHECK3D(ZV,'QLAP::ZV',PRECISION) ! IF ( GSOUTH ) THEN - !$acc kernels async present_cr(ZV) + !$acc kernels async DO JK=2,IKU-1 DO JI=1,IIU ZV(JI,IJB,JK)= (PY(JI,IJB,JK) - PY(JI,IJB-1,JK) - 0.5 * ( & @@ -556,7 +295,7 @@ IF(.NOT. L2D) THEN IF ( GNORTH ) THEN - !$acc kernels async present_cr(ZV) + !$acc kernels async DO JK=2,IKU-1 DO JI=1,IIU ZV(JI,IJE+1,JK)= (PY(JI,IJE+1,JK) - PY(JI,IJE+1-1,JK) - 0.5 * ( & @@ -570,38 +309,72 @@ IF(.NOT. L2D) THEN !$acc end kernels END IF !$acc wait - CALL MPPDB_CHECK3D(ZV,'QLAP::ZV/S',PRECISION) - CALL MPPDB_CHECK3D(ZV,'QLAP::ZV/N',PRECISION) +! CALL MPPDB_CHECK3D(ZV,'QLAP::ZV/S',PRECISION) +! CALL MPPDB_CHECK3D(ZV,'QLAP::ZV/N',PRECISION) ! ELSE ZV(:,:,:) = 0. ENDIF ! IF ( CEQNSYS == 'DUR' .OR. CEQNSYS == 'MAE' ) THEN - !$acc kernels present_cr(ZRHODJ) +#ifndef MNH_OPENACC +ZU(:,:,:) = MXM( PRHODJ(:,:,:) * XCPD * PTHETAV(:,:,:) ) * ZU(:,:,:) +#else + !$acc kernels ZRHODJ(:,:,:) = PRHODJ(:,:,:) * XCPD * PTHETAV(:,:,:) !$acc end kernels CALL MXM_DEVICE( ZRHODJ(:,:,:), ZMXM(:,:,:) ) !$acc kernels present(ZU) ZU(:,:,:) = ZMXM(:,:,:) * ZU(:,:,:) !$acc end kernels +#endif IF(.NOT. L2D) THEN +#ifndef MNH_OPENACC + ZV(:,:,:) = MYM( PRHODJ(:,:,:) * XCPD * PTHETAV(:,:,:) ) * ZV(:,:,:) +#else CALL MYM_DEVICE( ZRHODJ(:,:,:), ZMYM(:,:,:) ) !$acc kernels present(ZV) ZV(:,:,:) = ZMYM(:,:,:) * ZV(:,:,:) !$acc end kernels +#endif END IF - CALL MZM_DEVICE( ZRHODJ(:,:,:), ZMZM(:,:,:) ) +#ifndef MNH_OPENACC + ZW(:,:,:) = MZM( PRHODJ(:,:,:) * XCPD * PTHETAV(:,:,:) ) * GZ_M_W( 1, IKU, 1, PY(:,:,:), PDZZ(:,:,:) ) +#else + CALL MZM_DEVICE( ZRHODJ(:,:,:), ZMZM(:,:,:) ) CALL GZ_M_W_DEVICE( 1, IKU, 1, PY(:,:,:), PDZZ(:,:,:), ZGZMW(:,:,:) ) !$acc kernels present(ZW) ZW(:,:,:) = ZMZM(:,:,:) * ZGZMW(:,:,:) !$acc end kernels +#endif ELSE IF ( CEQNSYS == 'LHE' ) THEN +#ifndef MNH_OPENACC ZU(:,:,:) = MXM( PRHODJ(:,:,:) ) * ZU(:,:,:) +#else + CALL MXM_DEVICE( PRHODJ(:,:,:), ZMXM(:,:,:) ) + !$acc kernels present(ZU) + ZU(:,:,:) = ZMXM(:,:,:) * ZU(:,:,:) + !$acc end kernels +#endif IF(.NOT. L2D) THEN +#ifndef MNH_OPENACC ZV(:,:,:) = MYM( PRHODJ(:,:,:) ) * ZV(:,:,:) +#else + CALL MYM_DEVICE( PRHODJ(:,:,:), ZMYM(:,:,:) ) + !$acc kernels present(ZV) + ZV(:,:,:) = ZMYM(:,:,:) * ZV(:,:,:) + !$acc end kernels +#endif ENDIF +#ifndef MNH_OPENACC ZW(:,:,:) = MZM( PRHODJ(:,:,:) ) * GZ_M_W( 1, IKU, 1, PY(:,:,:), PDZZ(:,:,:) ) +#else + CALL MZM_DEVICE( PRHODJ(:,:,:), ZMZM(:,:,:) ) + CALL GZ_M_W_DEVICE( 1, IKU, 1, PY(:,:,:), PDZZ(:,:,:), ZGZMW(:,:,:) ) + !$acc kernels present(ZW) + ZW(:,:,:) = ZMZM(:,:,:) * ZGZMW(:,:,:) + !$acc end kernels +#endif END IF ! !------------------------------------------------------------------------------- @@ -609,9 +382,18 @@ END IF !* 2. COMPUTE THE DIVERGENCE ! ---------------------- ! +#ifndef MNH_OPENACC +NULLIFY(TZFIELDS_ll) +CALL ADD3DFIELD_ll( TZFIELDS_ll, ZU, 'QLAP::ZU' ) +CALL ADD3DFIELD_ll( TZFIELDS_ll, ZV, 'QLAP::ZV' ) +CALL ADD3DFIELD_ll( TZFIELDS_ll, ZW, 'QLAP::ZW' ) +CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll) +CALL CLEANLIST_ll(TZFIELDS_ll) +#else CALL GET_HALO_D(ZU,HNAME='QLAP::ZU') CALL GET_HALO_D(ZV,HNAME='QLAP::ZV') CALL GET_HALO_D(ZW,HNAME='QLAP::ZW') +#endif ! #ifndef MNH_OPENACC CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PQLAP) @@ -637,9 +419,15 @@ IF (LIBM) THEN ENDIF ! ENDIF -! + +!$acc end data + +#ifndef MNH_OPENACC +DEALLOCATE( ZU, ZV, ZW ) +#else !Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN -CALL MNH_MEM_RELEASE() +CALL MNH_MEM_RELEASE( 'QLAP' ) +#endif IF (MPPDB_INITIALIZED) THEN !Check all OUT arrays @@ -647,5 +435,4 @@ IF (MPPDB_INITIALIZED) THEN END IF !------------------------------------------------------------------------------- ! -END SUBROUTINE QLAP_DEVICE -#endif +END SUBROUTINE QLAP diff --git a/src/MNH/richardson.f90 b/src/MNH/richardson.f90 index b2a1d2f55266eda6f7002bd6a399d7ec97551702..753d4d54ad67f80d844513cf3d135b6cd06f31f2 100644 --- a/src/MNH/richardson.f90 +++ b/src/MNH/richardson.f90 @@ -118,7 +118,7 @@ END MODULE MODI_RICHARDSON !! !! AUTHOR !! ------ -!! P. HÅreil and J. Stein * Meteo France * +!! P. Hareil and J. Stein * Meteo France * !! !! MODIFICATIONS !! ------------- @@ -257,11 +257,7 @@ END IF ! DO JM = 1,KITR ! -#ifndef MNH_OPENACC - ZWORK(:,:,:) = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) ! Q (PHI) -#else - CALL QLAP_DEVICE( ZWORK, HLBCX, HLBCY, PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, PPHI ) ! Q (PHI) -#endif + CALL QLAP( ZWORK, HLBCX, HLBCY, PDXX, PDYY, PDZX, PDZY, PDZZ, PRHODJ, PTHETAV, PPHI ) ! Q (PHI) ! !$acc kernels ZCORREC(:,:,:) = 0. diff --git a/src/MNH/zconjgrad.f90 b/src/MNH/zconjgrad.f90 index ed927d553c471b034f9f631c5ecb89534f899be6..afa5c46ef543e416445ebc9edca7c0725264d565 100644 --- a/src/MNH/zconjgrad.f90 +++ b/src/MNH/zconjgrad.f90 @@ -120,7 +120,7 @@ END MODULE MODI_ZCONJGRAD !! !! AUTHOR !! ------ -!! P. HÅreil and J. Stein * Meteo France * +!! P. Hareil and J. Stein * Meteo France * !! !! MODIFICATIONS !! ------------- @@ -233,7 +233,7 @@ DO JM = 1,KITR ! -1 ! compute the vector DELTA = F * ( Y - Q ( PHI ) ) ! - ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) + CALL QLAP(ZWORK,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) ! Q (PHI) ! ZWORK = PY - ZWORK ! Y - Q (PHI) @@ -249,7 +249,7 @@ DO JM = 1,KITR IF (JM == 1) THEN ZP = ZDELTA ! P = DELTA at the first solver iteration ELSE - ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY, & + CALL QLAP(ZWORK,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY, & PDZZ,PRHODJ,PTHETAV,ZDELTA) ! Q ( DELTA ) CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! -1 PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZWORKD) ! F * Q ( DELTA ) @@ -263,7 +263,7 @@ DO JM = 1,KITR !* 2.4 compute LAMBDA ! ! - ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,& + CALL QLAP(ZWORK,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,& PDZZ,PRHODJ,PTHETAV,ZP) ! Q ( P ) CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,& ! -1 PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZWORKP) ! F * Q ( P ) diff --git a/src/MNH/zsolver.f90 b/src/MNH/zsolver.f90 index 2c618266e0fe1f016a886da6591295daf1753e04..f744786678305d3f9831293be93d01c8291e3836 100644 --- a/src/MNH/zsolver.f90 +++ b/src/MNH/zsolver.f90 @@ -265,16 +265,12 @@ CALL MNH_MEM_GET( ZRESIDUE , JIU, JJU, JKU ) !* 1.1 compute the vector: r^(0) = Q(PHI) - Y ! ! -#ifndef MNH_OPENACC -ZRESIDUE = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) - PY -#else !$acc data copyin( PTHETAV, PPHI ) -CALL QLAP_DEVICE(ZRESIDUE,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) +CALL QLAP(ZRESIDUE,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) !$acc end data !$acc kernels present(PY,ZRESIDUE) ZRESIDUE = ZRESIDUE - PY !$acc end kernels -#endif ! !* 1.2 compute the vector: p^(0) = F^(-1)*( Q(PHI) - Y ) ! @@ -286,13 +282,9 @@ CALL ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! !* 1.3 compute the vector: delta^(0) = Q ( p^(0) ) ! -#ifndef MNH_OPENACC -ZDELTA = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) -#else !$acc data copyin( PTHETAV, ZP ) -CALL QLAP_DEVICE(ZDELTA,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) +CALL QLAP(ZDELTA,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZP) !$acc end data -#endif ! !------------------------------------------------------------------------------- ! @@ -332,13 +324,9 @@ DO JM = 1,KITR ! !* 2.5 compute the auxiliary field: ksi = Q ( q ) ! -#ifndef MNH_OPENACC - ZKSI= QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) -#else !$acc data copyin( PTHETAV, ZQ ) - CALL QLAP_DEVICE(ZKSI,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) + CALL QLAP(ZKSI,HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,ZQ) !$acc end data -#endif ! -1 !* 2.6 compute the step ALPHA !