From 815b94a5b39bf39fe72f337ab41c80f542917958 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Mon, 22 May 2023 15:48:35 +0200
Subject: [PATCH] Philippe 22/05/2023: FFT: deduplicate FAST_TRANSPOSE,
FAST_SUBSTITUTION_2D, FAST_SUBSTITUTION_3D and FAST_SPREAD
---
src/MNH/compute_spectre.f90 | 36 +-------
src/MNH/fft_tools.f90 | 169 +++++++++++++++++++++++++++++++++++
src/MNH/flat_inv.f90 | 164 +---------------------------------
src/MNH/flat_invz.f90 | 171 ++----------------------------------
4 files changed, 180 insertions(+), 360 deletions(-)
create mode 100644 src/MNH/fft_tools.f90
diff --git a/src/MNH/compute_spectre.f90 b/src/MNH/compute_spectre.f90
index c900a114d..f40f291a5 100644
--- a/src/MNH/compute_spectre.f90
+++ b/src/MNH/compute_spectre.f90
@@ -64,6 +64,7 @@ USE MODD_SPECTRE, ONLY: LSMOOTH, LSTAT
USE MODE_FFT, ONLY: SET99
USE MODE_FFT55, ONLY: FFT55
+USE MODE_FFT_TOOLS, ONLY: FAST_TRANSPOSE
USE MODE_ll
USE MODE_MSG
USE MODE_SPLITTINGZ_ll
@@ -482,39 +483,6 @@ DEALLOCATE(ZVARP)
DEALLOCATE(ZAP)
DEALLOCATE(ZTRIGSX)
DEALLOCATE(ZTRIGSY)
-
-!
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-!
-CONTAINS
- SUBROUTINE FAST_TRANSPOSE(PX,PXT,KNI,KNJ,KNK)
- INTEGER :: KNI,KNJ,KNK ! 3D dimension of X and XT
- REAL, DIMENSION(KNI*KNJ,KNK) :: PX
- REAL, DIMENSION(KNJ*KNI,KNK) :: PXT
- !
- INTEGER :: IJI,II,IJ,IIJ ! index in array X and XT
- INTEGER :: JK
-!
- DO JK=1,KNK
- ! PERMUTATION(PX,PXT)
- !CDIR NODEP
- !OCL NOVREC
- DO IJI = 1, KNJ*KNI
- ! I,J Indice in XT array from linearised index IJI
- II = 1 + (IJI-1)/KNJ
- IJ = IJI - (II-1)*KNJ
- ! linearised index in X
- IIJ = II + (IJ-1)*KNI
- ! transposition
- PXT(IJI,JK) = PX(IIJ,JK)
-!
- END DO
- END DO
-!
-END SUBROUTINE FAST_TRANSPOSE
!
-!------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
END SUBROUTINE COMPUTE_SPECTRE
diff --git a/src/MNH/fft_tools.f90 b/src/MNH/fft_tools.f90
new file mode 100644
index 000000000..2fa896ad8
--- /dev/null
+++ b/src/MNH/fft_tools.f90
@@ -0,0 +1,169 @@
+!MNH_LIC Copyright 1994-2023 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.
+!-----------------------------------------------------------------
+! Creation: 22/05/2023
+! codes comes from flat_inv.f90 and flat_invz.f90 (deduplication of code)
+!-----------------------------------------------------------------
+! ##############
+MODULE MODE_FFT_TOOLS
+! ##############
+
+ IMPLICIT NONE
+
+ CONTAINS
+
+ SUBROUTINE FAST_TRANSPOSE( PX, PXT, KNI, KNJ, KNK )
+ INTEGER, INTENT(IN) :: KNI, KNJ, KNK ! 3D dimension of X and XT
+ REAL, DIMENSION(KNI*KNJ,KNK), INTENT(IN) :: PX
+ REAL, DIMENSION(KNJ*KNI,KNK), INTENT(OUT) :: PXT
+
+ INTEGER :: IJI,II,IJ,IIJ ! index in array X and XT
+ INTEGER :: JK
+
+ !$acc data present( PX, PXT )
+
+ !$acc kernels
+ DO JK=1,KNK
+ ! PERMUTATION(PX,PXT)
+ !CDIR NODEP
+ !OCL NOVREC
+ DO IJI = 1, KNJ*KNI
+ ! I,J Indice in XT array from linearised index IJI
+ II = 1 + (IJI-1)/KNJ
+ IJ = IJI - (II-1)*KNJ
+ ! linearised index in X
+ IIJ = II + (IJ-1)*KNI
+ ! transposition
+ PXT(IJI,JK) = PX(IIJ,JK)
+
+ END DO
+ END DO
+ !$acc end kernels
+
+ !$acc end data
+
+ END SUBROUTINE FAST_TRANSPOSE
+
+
+ SUBROUTINE FAST_SUBSTITUTION_3D( PBAND_YR, PBETX, PPBF, PGAM, PPCF, PAF, PBAND_Y, KIY, KJY, KKB, KKE, KKU )
+ INTEGER, INTENT(IN) :: KIY, KJY, KKB, KKE, KKU
+ REAL, DIMENSION (KIY*KJY,KKU),INTENT(OUT) :: PBAND_YR, PGAM
+ REAL, DIMENSION (KIY*KJY,KKU),INTENT(IN) :: PBAND_Y, PPBF, PAF
+ REAL, DIMENSION (KIY*KJY), INTENT(OUT) :: PBETX
+ REAL, DIMENSION (KKU), INTENT(IN) :: PPCF
+
+ INTEGER :: JK
+
+ !$acc data present( PBAND_YR, PGAM, PBAND_Y, PPBF, PAF, PBETX, PPCF )
+
+ !$acc kernels
+ !
+ ! initialization
+ !
+ PBAND_YR = 0.0
+ PBETX(:) = PPBF(:,KKB-1)
+ PBAND_YR(:,KKB-1) = PBAND_Y(:,KKB-1) &
+ / PBETX(:)
+ PGAM(:,1:KKB-1) = 0.
+ !
+ ! decomposition and forward substitution
+ !
+ DO JK = KKB,KKE+1
+ PGAM(:,JK) = PPCF(JK-1) / PBETX(:)
+ !
+ PBETX(:) = PPBF(:,JK) - &
+ PAF(:,JK)*PGAM(:,JK)
+ !
+ PBAND_YR(:,JK) = ( PBAND_Y(:,JK) - &
+ PAF(:,JK)*PBAND_YR(:,JK- 1) ) &
+ /PBETX(:)
+ !
+ END DO
+ !
+ ! backsubstitution
+ !
+ DO JK = KKE,KKB-1,-1
+ PBAND_YR(:,JK) = PBAND_YR(:,JK) - &
+ PGAM(:,JK+1)*PBAND_YR(:,JK+1)
+ END DO
+ !$acc end kernels
+
+ !$acc end data
+
+ END SUBROUTINE FAST_SUBSTITUTION_3D
+
+
+ SUBROUTINE FAST_SUBSTITUTION_2D( PBAND_YR, PBETX, PPBF, PGAM, PPCF, PAF, PBAND_Y, KIY, KJY, KKB, KKE, KKU )
+ INTEGER, INTENT(IN) :: KIY, KJY, KKB, KKE, KKU
+ REAL, DIMENSION(KIY,KJY,KKU), INTENT(OUT) :: PBAND_YR
+ REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PBAND_Y
+ REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PPBF
+ REAL, DIMENSION(KIY,KJY,KKU), INTENT(INOUT) :: PGAM
+ REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PAF
+ REAL, DIMENSION(KIY,KJY), INTENT(INOUT) :: PBETX
+ REAL, DIMENSION(KKU), INTENT(IN) :: PPCF
+
+ INTEGER :: JK
+
+ !$acc data present( PBAND_YR, PGAM, PBAND_Y, PPBF, PAF, PBETX, PPCF )
+
+ !$acc kernels
+ !
+ ! initialization
+ !
+ PBAND_YR = 0.0
+ PBETX(:,1) = PPBF(:,1,KKB-1)
+ PBAND_YR(:,1,KKB-1) = PBAND_Y(:,1,KKB-1) &
+ / PBETX(:,1)
+ !
+ ! decomposition and forward substitution
+ !
+ DO JK = KKB,KKE+1
+ PGAM(:,1,JK) = PPCF(JK-1) / PBETX(:,1)
+ !
+ PBETX(:,1) = PPBF(:,1,JK) - &
+ PAF(:,1,JK)*PGAM(:,1,JK)
+ !
+ PBAND_YR(:,1,JK) = ( PBAND_Y(:,1,JK) - &
+ PAF(:,1,JK)*PBAND_YR(:,1,JK- 1) ) &
+ /PBETX(:,1)
+ !
+ END DO
+ !
+ ! backsubstitution
+ !
+ DO JK = KKE,KKB-1,-1
+ PBAND_YR(:,1,JK) = PBAND_YR(:,1,JK) - &
+ PGAM(:,1,JK+1)*PBAND_YR(:,1,JK+1)
+ END DO
+ !$acc end kernels
+
+ !$acc end data
+
+ END SUBROUTINE FAST_SUBSTITUTION_2D
+
+
+ SUBROUTINE FAST_SPREAD( PTAB1D, PTAB3D, KIY, KJY, KKU )
+ INTEGER, INTENT(IN) :: KIY, KJY, KKU
+ REAL, DIMENSION(KKU), INTENT(IN) :: PTAB1D
+ REAL, DIMENSION(KIY*KJY,KKU), INTENT(OUT) :: PTAB3D
+
+ INTEGER :: JIJ,JK
+
+ !$acc data present( PTAB1D, PTAB3D )
+
+ !$acc kernels
+ DO JK=1,KKU
+ DO JIJ=1,KIY*KJY
+ PTAB3D(JIJ,JK) = PTAB1D(JK)
+ END DO
+ END DO
+ !$acc end kernels
+
+ !$acc end data
+
+ END SUBROUTINE FAST_SPREAD
+
+END MODULE MODE_FFT_TOOLS
diff --git a/src/MNH/flat_inv.f90 b/src/MNH/flat_inv.f90
index e9cb8d81a..80e4c312b 100644
--- a/src/MNH/flat_inv.f90
+++ b/src/MNH/flat_inv.f90
@@ -128,6 +128,7 @@ USE MODD_PARAMETERS
!
USE MODE_FFT, ONLY: FFT991
USE MODE_FFT55, ONLY: FFT55
+USE MODE_FFT_TOOLS
USE MODE_ll
#ifdef MNH_OPENACC
USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
@@ -473,11 +474,9 @@ IF (L2D) THEN
#ifdef MNH_OPENACC
CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'FLAT_INV', 'L2D=T not tested' )
#endif
- CALL FAST_SUBSTITUTION_2D(ZBAND_YR,ZBETX,PBF,ZGAM,PCF,ZAF &
- ,ZBAND_Y,IIY,IJY,IKU)
+ CALL FAST_SUBSTITUTION_2D( ZBAND_YR, ZBETX, PBF, ZGAM, PCF, ZAF, ZBAND_Y, IIY, IJY, IKB, IKE, IKU )
ELSE
- CALL FAST_SUBSTITUTION_3D(ZBAND_YRT,ZBETX,PBF,ZGAM,PCF,ZAF &
- ,ZBAND_YT,IIY,IJY,IKU)
+ CALL FAST_SUBSTITUTION_3D( ZBAND_YRT, ZBETX, PBF, ZGAM, PCF, ZAF, ZBAND_YT, IIY, IJY, IKB, IKE, IKU )
END IF
!
!
@@ -689,161 +688,4 @@ END IF
!
!-------------------------------------------------------------------------------
!
-CONTAINS
- SUBROUTINE FAST_TRANSPOSE(PX,PXT,KNI,KNJ,KNK)
- INTEGER, INTENT(IN) :: KNI,KNJ,KNK ! 3D dimension of X and XT
- REAL, DIMENSION(KNI*KNJ,KNK), INTENT(IN) :: PX
- REAL, DIMENSION(KNJ*KNI,KNK), INTENT(OUT) :: PXT
- !
- INTEGER :: IJI,II,IJ,IIJ ! index in array X and XT
- INTEGER :: JK
-!
-!$acc data present( PX, PXT )
-
-!$acc kernels
- DO JK=1,KNK
- ! PERMUTATION(PX,PXT)
- !CDIR NODEP
- !OCL NOVREC
- DO IJI = 1, KNJ*KNI
- ! I,J Indice in XT array from linearised index IJI
- II = 1 + (IJI-1)/KNJ
- IJ = IJI - (II-1)*KNJ
- ! linearised index in X
- IIJ = II + (IJ-1)*KNI
- ! transposition
- PXT(IJI,JK) = PX(IIJ,JK)
-
- END DO
- END DO
-!$acc end kernels
-
-!$acc end data
-!
-END SUBROUTINE FAST_TRANSPOSE
-
-SUBROUTINE FAST_SUBSTITUTION_3D(PBAND_YR,PBETX,PPBF,PGAM,PPCF,PAF &
- ,PBAND_Y,KIY,KJY,KKU)
-INTEGER, INTENT(IN) :: KIY,KJY,KKU
-REAL, DIMENSION (KIY*KJY,KKU),INTENT(OUT) :: PBAND_YR,PGAM
-REAL, DIMENSION (KIY*KJY,KKU),INTENT(IN) :: PBAND_Y,PPBF,PAF
-REAL, DIMENSION (KIY*KJY), INTENT(OUT) :: PBETX
-REAL, DIMENSION (KKU), INTENT(IN) :: PPCF
-!
-INTEGER :: JK
-!
-!
-! initialization
-!
-!
-!$acc data present( PBAND_YR, PGAM, PBAND_Y, PPBF, PAF, PBETX, PPCF )
-
-!$acc kernels
-PBAND_YR = 0.0
-PBETX(:) = PPBF(:,IKB-1)
-PBAND_YR(:,IKB-1) = PBAND_Y(:,IKB-1) &
- / PBETX(:)
-PGAM(:,1:IKB-1) = 0.
-!
-! decomposition and forward substitution
-!
-DO JK = IKB,IKE+1
- PGAM(:,JK) = PPCF(JK-1) / PBETX(:)
-!
- PBETX(:) = PPBF(:,JK) - &
- PAF(:,JK)*PGAM(:,JK)
-!
- PBAND_YR(:,JK) = ( PBAND_Y(:,JK) - &
- PAF(:,JK)*PBAND_YR(:,JK- 1) ) &
- /PBETX(:)
-!
-END DO
-!
-! backsubstitution
-!
-DO JK = IKE,IKB-1,-1
- PBAND_YR(:,JK) = PBAND_YR(:,JK) - &
- PGAM(:,JK+1)*PBAND_YR(:,JK+1)
-END DO
-!$acc end kernels
-
-!$acc end data
-!
-!
-END SUBROUTINE FAST_SUBSTITUTION_3D
-!
-SUBROUTINE FAST_SUBSTITUTION_2D(PBAND_YR,PBETX,PPBF,PGAM,PPCF,PAF &
- ,PBAND_Y,KIY,KJY,KKU)
-INTEGER, INTENT(IN) :: KIY, KJY, KKU
-REAL, DIMENSION(KIY,KJY,KKU), INTENT(OUT) :: PBAND_YR
-REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PBAND_Y
-REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PPBF
-REAL, DIMENSION(KIY,KJY,KKU), INTENT(INOUT) :: PGAM
-REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PAF
-REAL, DIMENSION(KIY,KJY), INTENT(INOUT) :: PBETX
-REAL, DIMENSION(KKU), INTENT(IN) :: PPCF
-!
-INTEGER :: JK
-!
-!
-! initialization
-!
-!
-!$acc data present( PBAND_YR, PGAM, PBAND_Y, PPBF, PAF, PBETX, PPCF )
-
-!$acc kernels
-PBAND_YR = 0.0
-PBETX(:,1) = PPBF(:,1,IKB-1)
-PBAND_YR(:,1,IKB-1) = PBAND_Y(:,1,IKB-1) &
- / PBETX(:,1)
-!
-! decomposition and forward substitution
-!
-DO JK = IKB,IKE+1
- PGAM(:,1,JK) = PPCF(JK-1) / PBETX(:,1)
-!
- PBETX(:,1) = PPBF(:,1,JK) - &
- PAF(:,1,JK)*PGAM(:,1,JK)
-!
- PBAND_YR(:,1,JK) = ( PBAND_Y(:,1,JK) - &
- PAF(:,1,JK)*PBAND_YR(:,1,JK- 1) ) &
- /PBETX(:,1)
-!
-END DO
-!
-! backsubstitution
-!
-DO JK = IKE,IKB-1,-1
- PBAND_YR(:,1,JK) = PBAND_YR(:,1,JK) - &
- PGAM(:,1,JK+1)*PBAND_YR(:,1,JK+1)
-END DO
-!$acc end kernels
-
-!$acc end data
-!
-!
-END SUBROUTINE FAST_SUBSTITUTION_2D
-
-SUBROUTINE FAST_SPREAD(PTAB1D,PTAB3D,KIY,KJY,KKU)
-INTEGER, INTENT(IN) :: KIY,KJY,KKU
-REAL, DIMENSION(KKU), INTENT(IN) :: PTAB1D
-REAL, DIMENSION(KIY*KJY,KKU), INTENT(OUT) :: PTAB3D
-
-INTEGER :: JIJ,JK
-!
-!$acc data present( PTAB1D, PTAB3D )
-
-!$acc kernels
-DO JK=1,KKU
- DO JIJ=1,KIY*KJY
- PTAB3D(JIJ,JK) = PTAB1D(JK)
- ENDDO
-ENDDO
-!$acc end kernels
-
-!$acc end data
-!
-END SUBROUTINE FAST_SPREAD
-!
-!------------------------------------------------------------------------------
END SUBROUTINE FLAT_INV
diff --git a/src/MNH/flat_invz.f90 b/src/MNH/flat_invz.f90
index ea7b6b3cd..fe6016718 100644
--- a/src/MNH/flat_invz.f90
+++ b/src/MNH/flat_invz.f90
@@ -151,6 +151,7 @@ SUBROUTINE FLAT_INVZ(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
USE MODE_FFT, ONLY: FFT991
USE MODE_FFT55, ONLY: FFT55
+ USE MODE_FFT_TOOLS
#ifdef MNH_OPENACC
USE MODE_MNH_ZWORK, ONLY: MNH_MEM_GET, MNH_MEM_POSITION_PIN, MNH_MEM_RELEASE
#endif
@@ -779,11 +780,9 @@ CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'FLAT_INVZ', 'OpenACC: HLBCY(1)==CYCL not po
#ifdef MNH_OPENACC
CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'FLAT_INVZ', 'L2D=T not tested' )
#endif
- CALL FAST_SUBSTITUTION_2D(ZBAND_YR,ZBETX,PBF,ZGAM,PCF,ZAF &
- ,ZBAND_Y,IIY,IJY,IKU)
+ CALL FAST_SUBSTITUTION_2D( ZBAND_YR, ZBETX, PBF, ZGAM, PCF, ZAF, ZBAND_Y, IIY, IJY, IKB, IKE, IKU )
ELSE
- CALL FAST_SUBSTITUTION_3D(ZBAND_YRT,ZBETX,PBF,ZGAM,PCF,ZAF &
- ,ZBAND_YT,IIY,IJY,IKU)
+ CALL FAST_SUBSTITUTION_3D( ZBAND_YRT, ZBETX, PBF, ZGAM, PCF, ZAF, ZBAND_YT, IIY, IJY, IKB, IKE, IKU )
END IF
END IF ! NZ_SPLITTING
@@ -856,9 +855,9 @@ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'FLAT_INVZ', 'L2D=T not tested' )
!Juan a faire CALL FAST_SUBSTITUTION_2D(ZBAND_YR,ZBETX,PBF,ZGAM,PCF,ZAF &
! ,ZBAND_Y,IIY,IJY,IKU)
ELSE
- CALL FAST_SUBSTITUTION_3D(ZBAND_SXP2_YP1_ZR,ZBETX_SXP2_YP1_Z,PBF_SXP2_YP1_Z,&
- ZGAM_SXP2_YP1_Z,PCF,ZAF_SXP2_YP1_Z,ZBAND_SXP2_YP1_Z,&
- II_SXP2_YP1_Z,IJ_SXP2_YP1_Z,IK_SXP2_YP1_Z)
+ CALL FAST_SUBSTITUTION_3D( ZBAND_SXP2_YP1_ZR, ZBETX_SXP2_YP1_Z, PBF_SXP2_YP1_Z, &
+ ZGAM_SXP2_YP1_Z, PCF, ZAF_SXP2_YP1_Z, ZBAND_SXP2_YP1_Z, &
+ II_SXP2_YP1_Z, IJ_SXP2_YP1_Z, IKB, IKE, IK_SXP2_YP1_Z )
END IF
END IF ! NZ_SPLITTING
!
@@ -1204,162 +1203,4 @@ CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'FLAT_INVZ', 'L2D=T not tested' )
!JUAN SCALASCA TEST
!
!-------------------------------------------------------------------------------
- !
-CONTAINS
- SUBROUTINE FAST_TRANSPOSE(PX,PXT,KNI,KNJ,KNK)
- INTEGER, INTENT(IN) :: KNI,KNJ,KNK ! 3D dimension of X and XT
- REAL, DIMENSION(KNI*KNJ,KNK), INTENT(IN) :: PX
- REAL, DIMENSION(KNJ*KNI,KNK), INTENT(OUT) :: PXT
- !
- INTEGER :: IJI,II,IJ,IIJ ! index in array X and XT
- INTEGER :: JK
- !
-!$acc data present( PX, PXT )
-
-!$acc kernels
- DO JK=1,KNK
- ! PERMUTATION(PX,PXT)
- !CDIR NODEP
- !OCL NOVREC
- DO IJI = 1, KNJ*KNI
- ! I,J Indice in XT array from linearised index IJI
- II = 1 + (IJI-1)/KNJ
- IJ = IJI - (II-1)*KNJ
- ! linearised index in X
- IIJ = II + (IJ-1)*KNI
- ! transposition
- PXT(IJI,JK) = PX(IIJ,JK)
-
- END DO
- END DO
-!$acc end kernels
-
-!$acc end data
- !
- END SUBROUTINE FAST_TRANSPOSE
-
- SUBROUTINE FAST_SUBSTITUTION_3D(PBAND_YR,PBETX,PPBF,PGAM,PPCF,PAF &
- ,PBAND_Y,KIY,KJY,KKU)
- INTEGER, INTENT(IN) :: KIY,KJY,KKU
- REAL, DIMENSION (KIY*KJY,KKU),INTENT(OUT) :: PBAND_YR,PGAM
- REAL, DIMENSION (KIY*KJY,KKU),INTENT(IN) :: PBAND_Y,PPBF,PAF
- REAL, DIMENSION (KIY*KJY), INTENT(OUT) :: PBETX
- REAL, DIMENSION (KKU), INTENT(IN) :: PPCF
- !
- INTEGER :: JK
- !
- !
- ! initialization
- !
- !
-!$acc data present( PBAND_YR, PGAM, PBAND_Y, PPBF, PAF, PBETX, PPCF )
-
-!$acc kernels
- PBAND_YR = 0.0
- PBETX(:) = PPBF(:,IKB-1)
- PBAND_YR(:,IKB-1) = PBAND_Y(:,IKB-1) &
- / PBETX(:)
- PGAM(:,1:IKB-1) = 0.
- !
- ! decomposition and forward substitution
- !
- DO JK = IKB,IKE+1
- PGAM(:,JK) = PPCF(JK-1) / PBETX(:)
- !
- PBETX(:) = PPBF(:,JK) - &
- PAF(:,JK)*PGAM(:,JK)
- !
- PBAND_YR(:,JK) = ( PBAND_Y(:,JK) - &
- PAF(:,JK)*PBAND_YR(:,JK- 1) ) &
- /PBETX(:)
- !
- END DO
- !
- ! backsubstitution
- !
- DO JK = IKE,IKB-1,-1
- PBAND_YR(:,JK) = PBAND_YR(:,JK) - &
- PGAM(:,JK+1)*PBAND_YR(:,JK+1)
- END DO
-!$acc end kernels
-
-!$acc end data
- !
- !
- END SUBROUTINE FAST_SUBSTITUTION_3D
- !
- SUBROUTINE FAST_SUBSTITUTION_2D(PBAND_YR,PBETX,PPBF,PGAM,PPCF,PAF &
- ,PBAND_Y,KIY,KJY,KKU)
- INTEGER, INTENT(IN) :: KIY, KJY, KKU
- REAL, DIMENSION(KIY,KJY,KKU), INTENT(OUT) :: PBAND_YR
- REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PBAND_Y
- REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PPBF
- REAL, DIMENSION(KIY,KJY,KKU), INTENT(INOUT) :: PGAM
- REAL, DIMENSION(KIY,KJY,KKU), INTENT(IN) :: PAF
- REAL, DIMENSION(KIY,KJY), INTENT(INOUT) :: PBETX
- REAL, DIMENSION(KKU), INTENT(IN) :: PPCF
- !
- INTEGER :: JK
- !
- !
- ! initialization
- !
- !
-!$acc data present( PBAND_YR, PGAM, PBAND_Y, PPBF, PAF, PBETX, PPCF )
-
-!$acc kernels
- PBAND_YR = 0.0
- PBETX(:,1) = PPBF(:,1,IKB-1)
- PBAND_YR(:,1,IKB-1) = PBAND_Y(:,1,IKB-1) &
- / PBETX(:,1)
- !
- ! decomposition and forward substitution
- !
- DO JK = IKB,IKE+1
- PGAM(:,1,JK) = PPCF(JK-1) / PBETX(:,1)
- !
- PBETX(:,1) = PPBF(:,1,JK) - &
- PAF(:,1,JK)*PGAM(:,1,JK)
- !
- PBAND_YR(:,1,JK) = ( PBAND_Y(:,1,JK) - &
- PAF(:,1,JK)*PBAND_YR(:,1,JK- 1) ) &
- /PBETX(:,1)
- !
- END DO
- !
- ! backsubstitution
- !
- DO JK = IKE,IKB-1,-1
- PBAND_YR(:,1,JK) = PBAND_YR(:,1,JK) - &
- PGAM(:,1,JK+1)*PBAND_YR(:,1,JK+1)
- END DO
-!$acc end kernels
-
-!$acc end data
- !
- !
- END SUBROUTINE FAST_SUBSTITUTION_2D
-
- SUBROUTINE FAST_SPREAD(PTAB1D,PTAB3D,KIY,KJY,KKU)
- INTEGER, INTENT(IN) :: KIY,KJY,KKU
- REAL, DIMENSION(KKU), INTENT(IN) :: PTAB1D
- REAL, DIMENSION(KIY*KJY,KKU), INTENT(OUT) :: PTAB3D
-
- INTEGER :: JIJ,JK
- !
-!$acc data present( PTAB1D, PTAB3D )
-
-!$acc kernels
- DO JK=1,KKU
- DO JIJ=1,KIY*KJY
- PTAB3D(JIJ,JK) = PTAB1D(JK)
- ENDDO
- ENDDO
-!$acc end kernels
-
-!$acc end data
- !
- END SUBROUTINE FAST_SPREAD
- !
- !------------------------------------------------------------------------------
END SUBROUTINE FLAT_INVZ
--
GitLab