diff --git a/src/ZSOLVER/gdiv.f90 b/src/ZSOLVER/gdiv.f90 deleted file mode 100644 index 7ebb32eb6b655f3c657dd90c0f3de9bdac0adaeb..0000000000000000000000000000000000000000 --- a/src/ZSOLVER/gdiv.f90 +++ /dev/null @@ -1,373 +0,0 @@ -!MNH_LIC Copyright 1994-2022 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 MODI_GDIV -! ################ -! -INTERFACE -! - SUBROUTINE GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PU,PV,PW,PGDIV) -! -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 -! - ! Field components -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PU ! along x -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PV ! along y -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PW ! along z -! -REAL, DIMENSION(:,:,:), INTENT(OUT) :: PGDIV ! divergence at - ! a mass point -! -END SUBROUTINE GDIV -! -END INTERFACE -! -END MODULE MODI_GDIV -! -! #################################################################### - SUBROUTINE GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PU,PV,PW,PGDIV) -! #################################################################### -! -!!**** *GDIV * - Compute J times the divergence of 1/J times a vector defined -!! by its cartesian components -!! -!! PURPOSE -!! ------- -! The purpose of this function is to compute J times the divergence of -! 1/J times a vector which cartesian components are (U, V, W). The result -! is localized at a mass point: -! -! GDIV = dxf (UC) + dyf (VC) + dzf (WC) -! -! where UC, VC, WC are the contravariant components of the vector. -! The array is completed outside the physical domain by the value of the -! normal component at the boundary. -! -!!** METHOD -!! ------ -!! First, we compute the contravariant components by using -!! the suboutine CONTRAV (The metric coefficients are dummy arguments). Then -!! we use the Shuman finite difference operators DXF, DYF, DZF to compute -!! the divergence. The result is localized at a mass point. -!! -!! EXTERNAL -!! -------- -!! SUBROUTINE CONTRAV : compute the contavariant components -!! -!! IMPLICIT ARGUMENTS -!! ------------------ -!! Module MODD_PARAMETERS: declaration of parameter variables -!! JPHEXT, JPVEXT: define the number of marginal points out of the -!! physical domain along horizontal and vertical directions respectively -!! Module MODD_CONF: model configurations -!! L2D: logical for 2D model version -!! Module MODI_SHUMAN : interface for the Shuman operators -!! Module MODI_CONTRAV : interface for the contravariant components -!! -!! REFERENCE -!! --------- -!! Book2 of documentation (routine GDIV) -!! -!! -!! AUTHOR -!! ------ -!! P. Hereil and J. Stein * Meteo France * -!! -!! MODIFICATIONS -!! ------------- -!! Original 11/07/94 -!! 17/07/97 ( J. Stein and V. Masson) initialize the corner -!! verticals -!! 30/09/97 ( J. Stein ) bug correction for the case of -!! non-vanishing oro. at the open lbc -!! Modification 15/06/98 (D.Lugato, R.Guivarch) Parallelisation -!! 22/08/02 (P Jabouille) simplification of parallel coding -!------------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -USE MODD_PARAMETERS -USE MODD_CONF -USE MODI_CONTRAV -! -USE MODE_ll -! -#ifdef MNH_OPENACC -USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE -#endif -! -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 -! - ! Field components -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PU ! along x -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PV ! along y -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PW ! along z -! -REAL, DIMENSION(:,:,:), INTENT(OUT) :: PGDIV ! divergence at - ! a mass point -! -!* 0.2 declarations of local variables -! - ! Contravariant components along: -REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZUC ! x -REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZVC ! y -REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZWC ! z -REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: Z1,Z2,Z3 !work arrays -! -INTEGER :: IIB ! indice I for the first inner mass point along x -INTEGER :: IIE ! indice I for the last inner mass point along x -INTEGER :: IJB ! indice J for the first inner mass point along y -INTEGER :: IJE ! indice J for the last inner mass point along y -INTEGER :: IKB ! indice K for the first inner mass point along z -INTEGER :: IKE ! indice K for the last inner mass point along z -! -INTEGER :: JI,JJ,JK ! loop indexes -! -INTEGER :: IIU,IJU,IKU -! -LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH -! -!------------------------------------------------------------------------------- -! -!* 1. COMPUTE LOOP BOUNDS -! ------------------- -! -! -CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) -IKB=1+JPVEXT -IKE=SIZE(PU,3) - JPVEXT -! -! -IIU=SIZE(PU,1) -IJU=SIZE(PU,2) -IKU=SIZE(PU,3) -! -#ifndef MNH_OPENACC -ALLOCATE(ZUC(IIU,IJU,IKU),ZVC(IIU,IJU,IKU),ZWC(IIU,IJU,IKU)) -ALLOCATE(Z1(IIU,IJU,IKU),Z2(IIU,IJU,IKU),Z3(IIU,IJU,IKU)) -#else -!Pin positions in the pools of MNH memory -CALL MNH_MEM_POSITION_PIN() - -CALL MNH_MEM_GET( ZUC, IIU, IJU, IKU ) -CALL MNH_MEM_GET( ZVC, IIU, IJU, IKU ) -CALL MNH_MEM_GET( ZWC, IIU, IJU, IKU ) -CALL MNH_MEM_GET( Z1, IIU, IJU, IKU ) -CALL MNH_MEM_GET( Z2, IIU, IJU, IKU ) -CALL MNH_MEM_GET( Z3, IIU, IJU, IKU ) -!$acc data present( ZUC, ZVC, ZWC, Z1, Z2, Z3 ) copyout( PGDIV ) -#endif -! -GWEST = ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() ) -GEAST = ( HLBCX(2) /= 'CYCL' .AND. LEAST_ll() ) -GSOUTH = ( .NOT. L2D .AND. HLBCY(1) /= 'CYCL' .AND. LSOUTH_ll() ) -GNORTH = ( .NOT. L2D .AND. HLBCY(2) /= 'CYCL' .AND. LNORTH_ll() ) -! -!------------------------------------------------------------------------------- -! -!* 2. COMPUTE THE CONTRAVARIANT COMPONENTS -! ------------------------------------ -! -!* 2.1 prepare the boundary conditions -! - -! -!$acc kernels - DO CONCURRENT ( JI=1:IIU,JJ=1:IJU ) - PU(JI,JJ,IKB-1)=PU(JI,JJ,IKB) - PU(JI,JJ,IKE+1)=PU(JI,JJ,IKE) - PV(JI,JJ,IKB-1)=PV(JI,JJ,IKB) - PV(JI,JJ,IKE+1)=PV(JI,JJ,IKE) - END DO -!$acc end kernels -! -! -!* 2.1 compute the contravariant components -! -#ifndef MNH_OPENACC -CALL CONTRAV(HLBCX,HLBCY,PU,PV,PW,PDXX,PDYY,PDZZ,PDZX,PDZY,ZUC,ZVC,ZWC,4) -#else -CALL CONTRAV_DEVICE(HLBCX,HLBCY,PU,PV,PW,PDXX,PDYY,PDZZ,PDZX,PDZY,ZUC,ZVC,ZWC,4, & - ODATA_ON_DEVICE=.TRUE.) -#endif -! -!------------------------------------------------------------------------------- -! -!* 3. COMPUTE THE DIVERGENCE -! ---------------------- -! -!$acc kernels -#ifdef MNH_COMPILER_NVHPC -!$acc loop independent collapse(3) -#endif -DO CONCURRENT (JI=1:IIU,JJ=1:IJU,JK=1:IKU) - PGDIV(JI,JJ,JK)=0. !usefull for the four corners and halo zones -ENDDO -! -#ifdef MNH_COMPILER_NVHPC -!$acc loop independent collapse(3) -#endif -DO CONCURRENT (JI=IIB:IIE,JJ=1:IJU,JK=1:IKU) - Z1(JI,JJ,JK)=ZUC(JI+IIB+1-(IIB) ,JJ,JK)-ZUC(JI,JJ,JK) -ENDDO -#ifdef MNH_COMPILER_NVHPC -!$acc loop independent collapse(3) -#endif -DO CONCURRENT (JI=1:IIU,JJ=IJB:IJE,JK=1:IKU) - Z2(JI,JJ,JK)=ZVC(JI,JJ+IJB+1-(IJB) ,JK)-ZVC(JI,JJ,JK) -ENDDO -#ifdef MNH_COMPILER_NVHPC -!$acc loop independent collapse(3) -#endif -DO CONCURRENT (JI=1:IIU,JJ=1:IJU,JK=IKB:IKE) - Z3(JI,JJ,JK)=ZWC(JI,JJ,JK+IKB+1-(IKB) )-ZWC(JI,JJ,JK) -ENDDO -! -PGDIV(IIB:IIE,IJB:IJE,IKB:IKE)= Z1(IIB:IIE,IJB:IJE,IKB:IKE) + & - Z2(IIB:IIE,IJB:IJE,IKB:IKE) + & - Z3(IIB:IIE,IJB:IJE,IKB:IKE) - ! only the divergences computed - ! in the inner mass points are meaningful -!$acc end kernels -! -!------------------------------------------------------------------------------- -! -!* 4. SET DIVERGENCE AT THE OUTER POINTS -! ---------------------------------- -! -!* 4.1 set divergence at the upper and lower boundary -! -! we set the divergence equal to the vertical contravariant component above -! and under the physical domain -!$acc kernels async -DO JJ=IJB,IJE - DO JI=IIB,IIE - PGDIV(JI,JJ,IKB-1)=ZWC(JI,JJ,IKB) - PGDIV(JI,JJ,IKE+1)=ZWC(JI,JJ,IKE+1) - END DO -END DO -!$acc end kernels -! -!* 4.2 set divergence at the lateral boundaries -! -! we set the divergence equal to the horizontal contravariant component at -! the right and the left of the physical domain in both horizontal directions -! for non-periodic cases -! -IF( GWEST ) THEN - !$acc kernels async - DO JK=IKB,IKE - DO JJ=IJB,IJE - PGDIV(IIB-1,JJ,JK)=ZUC(IIB,JJ,JK) - END DO - END DO - !$acc end kernels -END IF -! -IF( GEAST ) THEN - !$acc kernels async - DO JK=IKB,IKE - DO JJ=IJB,IJE - PGDIV(IIE+1,JJ,JK)=ZUC(IIE+1,JJ,JK) - END DO - END DO - !$acc end kernels -END IF -! -! -IF ( GSOUTH ) THEN - !$acc kernels async - DO JK=IKB,IKE - DO JI=IIB,IIE - PGDIV(JI,IJB-1,JK)=ZVC(JI,IJB,JK) - END DO - END DO - !$acc end kernels -END IF -! -IF ( GNORTH ) THEN - !$acc kernels async - DO JK=IKB,IKE - DO JI=IIB,IIE - PGDIV(JI,IJE+1,JK)=ZVC(JI,IJE+1,JK) - END DO - END DO - !$acc end kernels -END IF -! -! wait on GPU for all boundary condition update -!$acc wait -! -!* 4.3 set divergence at the corner points -! -! it is the following of the condition of copy the horizontal component -! under the bottom of the model -! -IF( GWEST ) THEN - !$acc kernels async - PGDIV(IIB-1,IJB:IJE,IKB-1)=PGDIV(IIB-1,IJB:IJE,IKB) - PGDIV(IIB-1,IJB:IJE,IKE+1)=PGDIV(IIB-1,IJB:IJE,IKE) - !$acc end kernels -END IF -! -IF ( GEAST ) THEN - !$acc kernels async - PGDIV(IIE+1,IJB:IJE,IKB-1)=PGDIV(IIE+1,IJB:IJE,IKB) - PGDIV(IIE+1,IJB:IJE,IKE+1)=PGDIV(IIE+1,IJB:IJE,IKE) - !$acc end kernels -END IF -! -IF ( GSOUTH ) THEN - !$acc kernels async - PGDIV(IIB:IIE,IJB-1,IKB-1)=PGDIV(IIB:IIE,IJB-1,IKB) - PGDIV(IIB:IIE,IJB-1,IKE+1)=PGDIV(IIB:IIE,IJB-1,IKE) - !$acc end kernels -END IF -! -IF ( GNORTH ) THEN - !$acc kernels async - PGDIV(IIB:IIE,IJE+1,IKB-1)=PGDIV(IIB:IIE,IJE+1,IKB) - PGDIV(IIB:IIE,IJE+1,IKE+1)=PGDIV(IIB:IIE,IJE+1,IKE) - !$acc end kernels -END IF -! -! wait on GPU for all corner update -!$acc wait -! -#ifndef MNH_OPENACC -DEALLOCATE( ZUC, ZVC, ZWC, Z1, Z2, Z3 ) -#else -!$acc end data -!Release all memory allocated with MNH_MEM_GET calls since last call to MNH_MEM_POSITION_PIN -CALL MNH_MEM_RELEASE() -#endif -! -!------------------------------------------------------------------------------- -! -END SUBROUTINE GDIV