From b6dd5b9be68b799bc629e06d1e52413258ff989a Mon Sep 17 00:00:00 2001
From: Juan Escobar <escj@aero.obs-mip.fr>
Date: Thu, 22 Nov 2018 14:41:26 +0100
Subject: [PATCH] Juan 22/11/2018: add modif for full *_ZS array without
average on Z , not really working yet !
---
anel_balancen.f90 | 23 ++++++-
ini_dynamics.f90 | 23 +++++--
ini_modeln.f90 | 12 +++-
ini_spectren.f90 | 4 +-
modd_dynn.f90 | 24 ++++++++
modeln.f90 | 6 +-
pressurez.f90 | 23 +++++--
tridz.f90 | 97 +++++++++++++++++++++--------
zconjgrad.f90 | 6 +-
zsolver.f90 | 28 +++++++--
zsolver_inv.f90 | 153 +++++++++++++++++++++++++++++++++++++++-------
11 files changed, 328 insertions(+), 71 deletions(-)
diff --git a/anel_balancen.f90 b/anel_balancen.f90
index 6c64978f6..2d39a7aab 100644
--- a/anel_balancen.f90
+++ b/anel_balancen.f90
@@ -179,7 +179,10 @@ INTEGER :: IMI ! model index
!JUAN
INTEGER :: IIU_B,IJU_B,IKU
INTEGER :: IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll
-REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZBFB,ZBF_SXP2_YP1_Z
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZBFB,ZBF_SXP2_YP1_Z
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZAF_ZS,ZBF_ZS,ZCF_ZS
+REAL, DIMENSION(:,:) , ALLOCATABLE :: ZDXATH_ZS,ZDYATH_ZS
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRHO_ZS
!JUAN
!
INTEGER :: IINFO_ll
@@ -202,8 +205,17 @@ ALLOCATE(ZTRIGSY(3*(NJMAX_ll+2*JPHEXT)))
IKU=SIZE(XRHODJ,3)
CALL GET_DIM_EXT_ll('B',IIU_B,IJU_B)
ALLOCATE(ZBFB(IIU_B,IJU_B,IKU))
+!
+ALLOCATE(ZAF_ZS(IIU_B,IJU_B,IKU))
+ALLOCATE(ZBF_ZS(IIU_B,IJU_B,IKU))
+ALLOCATE(ZCF_ZS(IIU_B,IJU_B,IKU))
+ALLOCATE(ZDXATH_ZS(IIU_B,IJU_B))
+ALLOCATE(ZDYATH_ZS(IIU_B,IJU_B))
+ALLOCATE(ZRHO_ZS(IIU_B,IJU_B,IKU))
+!
CALL GET_DIM_EXTZ_ll('SXP2_YP1_Z',IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll)
ALLOCATE(ZBF_SXP2_YP1_Z(IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll))
+!
!JUAN Z_SPLITING
CALL MPPDB_CHECK3D(XRHODJ,"anel_balancen1-::XRHODJ",PRECISION)
CALL MPPDB_CHECK3D(XUT,"anel_balancen1-::XUT",PRECISION)
@@ -216,7 +228,9 @@ CALL MPPDB_CHECK3D(XUT,"anel_balancen1-::XUT",PRECISION)
!
CALL TRIDZ(CLBCX,CLBCY,XMAP,XDXHAT,XDYHAT,ZDXHATM,ZDYHATM,ZRHOM, &
ZAF,ZCF,ZTRIGSX,ZTRIGSY,IIFAXX,IIFAXY,XRHODJ,XTHVREF,XZZ,ZBFY,&
- ZBFB,ZBF_SXP2_YP1_Z)
+ ZBFB,ZBF_SXP2_YP1_Z, &
+ ZAF_ZS,ZBF_ZS,ZCF_ZS, &
+ ZDXATH_ZS,ZDYATH_ZS,ZRHO_ZS) !JUAN FULL ZSOLVER
CALL MPPDB_CHECK3D(XRHODJ,"anel_balancen1-after TRIDZ::XRHODJ",PRECISION)
!
!-------------------------------------------------------------------------------
@@ -280,7 +294,10 @@ CALL PRESSUREZ(CLBCX,CLBCY,CPRESOPT,NITR,LITRADJ,ITCOUNT,XRELAX,IMI, &
IRR,IRRL,IRRI,ZDRYMASST,ZREFMASS,ZMASS_O_PHI0, &
ZTH,ZRR,XRHODREF,XTHVREF,XRVREF,XEXNREF, XLINMASS, &
ZRU,ZRV,ZRW,ZPABST, &
- ZBFB,ZBF_SXP2_YP1_Z,PRESIDUAL )
+ ZBFB,ZBF_SXP2_YP1_Z, &
+ XAF_ZS,XBF_ZS,XCF_ZS, &
+ XDXATH_ZS,XDYATH_ZS,XRHO_ZS, &
+ PRESIDUAL )
!
CALL MPPDB_CHECK3D(XRHODJ,"anel_balancen3.2-after pressurez halo::XRHODJ",PRECISION)
CALL MPPDB_CHECK3D(ZRU,"anel_balancen3.2-after pressurez::ZRU",PRECISION)
diff --git a/ini_dynamics.f90 b/ini_dynamics.f90
index e4d00f5bb..4e2d0ff8c 100644
--- a/ini_dynamics.f90
+++ b/ini_dynamics.f90
@@ -26,8 +26,10 @@ SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ, &
PALK,PALKW,KALBOT,PALKBAS,PALKWBAS,KALBAS, &
OMASK_RELAX,PKURELAX, PKVRELAX, PKWRELAX, &
PDK2U,PDK4U,PDK2TH,PDK4TH,PDK2SV,PDK4SV,OZDIFFU,PZDIFFU_HALO2,&
- PBFB,&
- PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
+ PBFB, &
+ PBF_SXP2_YP1_Z, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS) !JUAN FULL ZSOLVER
! intent in arguments
!
USE MODE_TYPE_ZDIFFU
@@ -172,6 +174,10 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB ! elements of the tri-diag matrix
! on an b-slice of global physical domain
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide
! matrix in the pressure eq.
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(OUT) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRHO_ZS
+
END SUBROUTINE INI_DYNAMICS
!
END INTERFACE
@@ -197,8 +203,10 @@ SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ, &
PALK,PALKW,KALBOT,PALKBAS,PALKWBAS,KALBAS, &
OMASK_RELAX,PKURELAX, PKVRELAX, PKWRELAX, &
PDK2U,PDK4U,PDK2TH,PDK4TH,PDK2SV,PDK4SV,OZDIFFU,PZDIFFU_HALO2,&
- PBFB,&
- PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
+ PBFB, &
+ PBF_SXP2_YP1_Z, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS) !JUAN FULL ZSOLVER
! ######################################################################
!
!!**** *INI_DYNAMICS* - routine to initialize the parameters for the dynamics
@@ -440,6 +448,9 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB ! elements of the tri-diag matrix
! on an b-slice of global physical domain
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide
! matrix in the pressure eq.
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(OUT) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRHO_ZS
!
!* 0.2 declarations of local variables
!
@@ -495,7 +506,9 @@ IF (.NOT.L1D) THEN
PMAP,PDXHAT,PDYHAT,PDXHATM,PDYHATM,PRHOM,PAF, &
PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
PRHODJ,PTHVREF,PZZ,PBFY,PBFB, &
- PBF_SXP2_YP1_Z)
+ PBF_SXP2_YP1_Z, &
+ PAF_ZS,PBF_ZS,PCF_ZS,&
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS) !JUAN FULL ZSOLVER
END IF
!
!
diff --git a/ini_modeln.f90 b/ini_modeln.f90
index c16588e4a..4697cf476 100644
--- a/ini_modeln.f90
+++ b/ini_modeln.f90
@@ -899,6 +899,14 @@ ELSE
END IF
CALL GET_DIM_EXT_ll('B',IIU_B,IJU_B)
ALLOCATE(XBFB(IIU_B,IJU_B,IKU))
+
+ALLOCATE(XAF_ZS(IIU_B,IJU_B,IKU))
+ALLOCATE(XBF_ZS(IIU_B,IJU_B,IKU))
+ALLOCATE(XCF_ZS(IIU_B,IJU_B,IKU))
+ALLOCATE(XDXATH_ZS(IIU_B,IJU_B))
+ALLOCATE(XDYATH_ZS(IIU_B,IJU_B))
+ALLOCATE(XRHO_ZS(IIU_B,IJU_B,IKU))
+
CALL GET_DIM_EXTZ_ll('SXP2_YP1_Z',IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll)
ALLOCATE(XBF_SXP2_YP1_Z(IIU_SXP2_YP1_Z_ll,IJU_SXP2_YP1_Z_ll,IKU_SXP2_YP1_Z_ll))
ALLOCATE(XAF(IKU),XCF(IKU))
@@ -1967,7 +1975,9 @@ CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT, &
LMASK_RELAX,XKURELAX,XKVRELAX,XKWRELAX, &
XDK2U,XDK4U,XDK2TH,XDK4TH,XDK2SV,XDK4SV, &
LZDIFFU,XZDIFFU_HALO2, &
- XBFB,XBF_SXP2_YP1_Z )
+ XBFB,XBF_SXP2_YP1_Z, &
+ XAF_ZS,XBF_ZS,XCF_ZS, &
+ XDXATH_ZS,XDYATH_ZS,XRHO_ZS )
!
!
!* 16.1 Initialize the XDRAG array
diff --git a/ini_spectren.f90 b/ini_spectren.f90
index d2564f410..ee95c73f5 100644
--- a/ini_spectren.f90
+++ b/ini_spectren.f90
@@ -938,7 +938,9 @@ CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT, &
LMASK_RELAX,XKURELAX,XKVRELAX,XKWRELAX, &
XDK2U,XDK4U,XDK2TH,XDK4TH,XDK2SV,XDK4SV, &
LZDIFFU,XZDIFFU_HALO2, &
- XBFB,XBF_SXP2_YP1_Z )
+ XBFB,XBF_SXP2_YP1_Z, &
+ XAF_ZS,XBF_ZS,XCF_ZS, &
+ XDXATH_ZS,XDYATH_ZS,XRHO_ZS)
!
!-------------------------------------------------------------------------------
!
diff --git a/modd_dynn.f90 b/modd_dynn.f90
index c3b64335c..251e8cbd6 100644
--- a/modd_dynn.f90
+++ b/modd_dynn.f90
@@ -65,6 +65,9 @@ TYPE DYN_t
REAL, DIMENSION(:,:,:), POINTER :: XBF_SXP2_YP1_Z=>NULL() ! Vectors giving the non
REAL, DIMENSION(:,:,:), POINTER :: XBF=>NULL() ! vanishing elements of the
REAL, DIMENSION(:), POINTER :: XAF=>NULL(),XCF=>NULL() ! tri-diag matrix in the pressure equation
+ REAL, DIMENSION(:,:,:), POINTER :: XAF_ZS=>NULL(),XBF_ZS=>NULL(),XCF_ZS=>NULL() ! coef for Zsolver
+ REAL, DIMENSION(:,:) , POINTER :: XDXATH_ZS=>NULL(),XDYATH_ZS=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XRHO_ZS=>NULL()
!
! Arrays of sinus or cosinus
! values for the FFT
@@ -186,6 +189,9 @@ REAL, DIMENSION(:,:,:), POINTER :: XBFB=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XBF_SXP2_YP1_Z=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XBF=>NULL()
REAL, DIMENSION(:), POINTER :: XAF=>NULL(),XCF=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XAF_ZS=>NULL(),XBF_ZS=>NULL(),XCF_ZS=>NULL()
+REAL, DIMENSION(:,:) , POINTER :: XDXATH_ZS=>NULL(),XDYATH_ZS=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XRHO_ZS=>NULL()
REAL, DIMENSION(:), POINTER :: XTRIGSX=>NULL()
REAL, DIMENSION(:), POINTER :: XTRIGSY=>NULL()
INTEGER, DIMENSION(:), POINTER :: NIFAXX=>NULL()
@@ -268,6 +274,15 @@ DYN_MODEL(KFROM)%XBF_SXP2_YP1_Z=>XBF_SXP2_YP1_Z
DYN_MODEL(KFROM)%XBF=>XBF
DYN_MODEL(KFROM)%XAF=>XAF
DYN_MODEL(KFROM)%XCF=>XCF
+
+DYN_MODEL(KFROM)%XAF_ZS=>XAF_ZS
+DYN_MODEL(KFROM)%XBF_ZS=>XBF_ZS
+DYN_MODEL(KFROM)%XCF_ZS=>XCF_ZS
+
+DYN_MODEL(KFROM)%XDXATH_ZS=>XDXATH_ZS
+DYN_MODEL(KFROM)%XDYATH_ZS=>XDYATH_ZS
+DYN_MODEL(KFROM)%XRHO_ZS=>XRHO_ZS
+
DYN_MODEL(KFROM)%XTRIGSX=>XTRIGSX
DYN_MODEL(KFROM)%XTRIGSY=>XTRIGSY
DYN_MODEL(KFROM)%XRHOM=>XRHOM
@@ -289,6 +304,15 @@ XBF_SXP2_YP1_Z=>DYN_MODEL(KTO)%XBF_SXP2_YP1_Z
XBF=>DYN_MODEL(KTO)%XBF
XAF=>DYN_MODEL(KTO)%XAF
XCF=>DYN_MODEL(KTO)%XCF
+
+XAF_ZS=>DYN_MODEL(KTO)%XAF_ZS
+XBF_ZS=>DYN_MODEL(KTO)%XBF_ZS
+XCF_ZS=>DYN_MODEL(KTO)%XCF_ZS
+
+XDXATH_ZS=>DYN_MODEL(KFROM)%XDXATH_ZS
+XDYATH_ZS=>DYN_MODEL(KFROM)%XDYATH_ZS
+XRHO_ZS=>DYN_MODEL(KFROM)%XRHO_ZS
+
XTRIGSX=>DYN_MODEL(KTO)%XTRIGSX
XTRIGSY=>DYN_MODEL(KTO)%XTRIGSY
NIFAXX=>DYN_MODEL(KTO)%NIFAXX
diff --git a/modeln.f90 b/modeln.f90
index 6836f74fa..823724640 100644
--- a/modeln.f90
+++ b/modeln.f90
@@ -1671,8 +1671,10 @@ CALL MPPDB_CHECK3DM("before pressurez:XRU/V/WS",PRECISION,XRUS,XRVS,XRWS)
NRR,NRRL,NRRI,XDRYMASST,XREFMASS,XMASS_O_PHI0, &
XTHT,XRT,XRHODREF,XTHVREF,XRVREF,XEXNREF, XLINMASS, &
XRUS, XRVS, XRWS, XPABST, &
- XBFB,&
- XBF_SXP2_YP1_Z) !JUAN Z_SPLITING
+ XBFB, &
+ XBF_SXP2_YP1_Z, &
+ XAF_ZS,XBF_ZS,XCF_ZS, &
+ XDXATH_ZS,XDYATH_ZS,XRHO_ZS) !JUAN FULL ZSOLVER
!
XRUS_PRES = XRUS - XRUS_PRES + ZRUS
XRVS_PRES = XRVS - XRVS_PRES + ZRVS
diff --git a/pressurez.f90 b/pressurez.f90
index 7ef9fb3e7..12071e18c 100644
--- a/pressurez.f90
+++ b/pressurez.f90
@@ -18,7 +18,9 @@ INTERFACE
PRUS,PRVS,PRWS,PPABST, &
PBFB, &
PBF_SXP2_YP1_Z, &
- PRESIDUAL ) !JUAN Z_SPLITING
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS, &
+ PRESIDUAL) !JUAN FULL ZSOLVER
!
IMPLICIT NONE
!
@@ -96,6 +98,11 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PBFB ! elements of the tri-diag b-sl
! matrix in the pressure eq.
REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide
! matrix in the pressure eq.
+
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(IN) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_ZS
+
REAL, OPTIONAL :: PRESIDUAL
!JUAN Z_SPLITING
END SUBROUTINE PRESSUREZ
@@ -113,7 +120,9 @@ END MODULE MODI_PRESSUREZ
PRUS,PRVS,PRWS,PPABST, &
PBFB, &
PBF_SXP2_YP1_Z, &
- PRESIDUAL ) !JUAN Z_SPLITING
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS, &
+ PRESIDUAL) !JUAN FULL ZSOLVER
! ######################################################################
!
!!**** *PRESSUREZ * - solve the pressure equation and add the pressure term
@@ -330,6 +339,10 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PBFB ! elements of the tri-diag b-sl
! matrix in the pressure eq.
REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide
! matrix in the pressure eq.
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(IN) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_ZS
+
REAL, OPTIONAL :: PRESIDUAL
!JUAN Z_SPLITING
!
@@ -513,8 +526,10 @@ ELSE
!
CASE('ZSOLV') ! Conjugate Residual method
CALL ZSOLVER(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV, &
- PDXHATM,PDYHATM,PRHOT,PAF,PBFB,PCF,PTRIGSX,PTRIGSY, &
- KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
+ PDXHATM,PDYHATM,PRHOT,PAF,PBFB,PCF,PTRIGSX,PTRIGSY, &
+ KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS)
!
CASE('ZRESI') ! Conjugate Residual method
CALL CONRESOLZ(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV, &
diff --git a/tridz.f90 b/tridz.f90
index b5e0161d4..8c27aaf58 100644
--- a/tridz.f90
+++ b/tridz.f90
@@ -12,8 +12,10 @@ INTERFACE
SUBROUTINE TRIDZ(HLBCX,HLBCY, &
PMAP,PDXHAT,PDYHAT,PDXHATM,PDYHATM,PRHOM, &
PAF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
- PRHODJ,PTHVREF,PZZ,PBFY,PBFB,&
- PBF_SXP2_YP1_Z)!JUAN Z_SPLITING
+ PRHODJ,PTHVREF,PZZ,PBFY,PBFB, &
+ PBF_SXP2_YP1_Z, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS) !JUAN FULL ZSOLVER
!
IMPLICIT NONE
!
@@ -48,6 +50,9 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB ! elements (bsplit slide) of the
! matrix in the pressure eq.
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide
! matrix in the pressure eq
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(OUT) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRHO_ZS
!JUAN
!
! arrays of sin or cos values
@@ -71,8 +76,10 @@ END MODULE MODI_TRIDZ
SUBROUTINE TRIDZ(HLBCX,HLBCY, &
PMAP,PDXHAT,PDYHAT,PDXHATM,PDYHATM,PRHOM, &
PAF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
- PRHODJ,PTHVREF,PZZ,PBFY,PBFB,&
- PBF_SXP2_YP1_Z) !JUAN Z_SPLITING
+ PRHODJ,PTHVREF,PZZ,PBFY,PBFB, &
+ PBF_SXP2_YP1_Z, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS) !JUAN FULL ZSOLVER
! ####################################################################
!
!!**** *TRIDZ * - Compute coefficients for the flat operator
@@ -229,6 +236,9 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBFB ! elements (bsplit slide) of the
! matrix in the pressure eq.
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PBF_SXP2_YP1_Z ! elements of the tri-diag. SXP2_YP1_Z-slide
! matrix in the pressure eq.
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(OUT) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRHO_ZS
!JUAN
!
! arrays of sin or cos values
@@ -292,6 +302,7 @@ INTEGER :: IXMODE_SXP2_YP1_Z_ll,IYMODE_SXP2_YP1_Z_ll ! number of modes in the x
!JUAN16
!TYPE(DOUBLE_DOUBLE) , DIMENSION (SIZE(PZZ,3)) :: ZRHOM_ll , ZDZM_ll
REAL, ALLOCATABLE, DIMENSION(:,:) :: ZRHOM_2D , ZDZM_2D
+REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ZDZM_ZS
!JUAN16
!
!
@@ -362,63 +373,59 @@ ZGWNY = XPI/REAL(IJMAX_ll)
!
ZINVMEAN = 1./REAL(IIMAX_ll*IJMAX_ll)
!JUAN16
-ALLOCATE(ZRHOM_2D(IIB:IIE, IJB:IJE))
+ALLOCATE(ZRHOM_2D(IIB:IIE, IJB:IJE))
+PRHO_ZS = 1.0
!
DO JK = 1,SIZE(PZZ,3)
IF ( CEQNSYS == 'DUR' .OR. CEQNSYS == 'MAE' ) THEN
DO JJ = IJB,IJE
DO JI = IIB,IIE
- ZRHOM_2D(JI,JJ) = PRHODJ(JI,JJ,JK)*XCPD*PTHVREF(JI,JJ,JK)*ZINVMEAN
+ ZRHOM_2D(JI,JJ) = PRHODJ(JI,JJ,JK)*XCPD*PTHVREF(JI,JJ,JK)
+ PRHO_ZS(JI,JJ,JK) = ZRHOM_2D(JI,JJ)
END DO
END DO
ELSEIF ( CEQNSYS == 'LHE' ) THEN
DO JJ = IJB,IJE
DO JI = IIB,IIE
- ZRHOM_2D(JI,JJ) = PRHODJ(JI,JJ,JK)*ZINVMEAN
+ ZRHOM_2D(JI,JJ) = PRHODJ(JI,JJ,JK)
+ PRHO_ZS(JI,JJ,JK) = ZRHOM_2D(JI,JJ)
END DO
END DO
END IF
! global sum
- PRHOM(JK) = SUM_DD_R2_ll(ZRHOM_2D)
+ PRHOM(JK) = SUM_DD_R2_ll(ZRHOM_2D) * ZINVMEAN
END DO
-!
-! global sum
-!CALL REDUCESUM_ll(ZRHOM_ll,IINFO_ll)
-!PRHOM = ZRHOM_ll
-!JUAN16
!
!------------------------------------------------------------------------------
!
!* 3. COMPUTE THE MEAN INCREMENT BETWEEN Z LEVELS
! -------------------------------------------
!
-!JUAN16
-!ZDZM_ll = 0.
ALLOCATE(ZDZM_2D(IIB:IIE, IJB:IJE))
+ALLOCATE(ZDZM_ZS(IIB:IIE, IJB:IJE,IKU))
+ZDZM_ZS = 1.0
!
DO JK = IKB-1,IKE
DO JJ = IJB,IJE
DO JI = IIB,IIE
- ZDZM_2D(JI,JJ) = (PZZ(JI,JJ,JK+1)-PZZ(JI,JJ,JK))*ZINVMEAN
+ ZDZM_2D(JI,JJ) = (PZZ(JI,JJ,JK+1)-PZZ(JI,JJ,JK))
+ ZDZM_ZS(JI,JJ,JK) = ZDZM_2D(JI,JJ)
END DO
END DO
- ZDZM(JK) = SUM_DD_R2_ll(ZDZM_2D)
+ ZDZM(JK) = SUM_DD_R2_ll(ZDZM_2D) * ZINVMEAN
END DO
ZDZM(IKE+1) = ZDZM(IKE)
-!
-! global sum
-!CALL REDUCESUM_ll(ZDZM_ll,IINFO_ll)
-!ZDZM = ZDZM_ll
-!JUAN16
-!
+ZDZM_ZS(:,:,IKE+1) = ZDZM_ZS(:,:,IKE)
!
! vertical average to arrive at a w-level
DO JK = IKE+1,IKB,-1
ZDZM(JK) = (ZDZM(JK) + ZDZM(JK-1))*0.5
+ ZDZM_ZS(IIB:IIE,IJB:IJE,JK) = (ZDZM_ZS(IIB:IIE,IJB:IJE,JK) + ZDZM_ZS(IIB:IIE,IJB:IJE,JK-1)) * 0.5
END DO
!
ZDZM(IKB-1) = -999.
+ZDZM_ZS(IIB:IIE,IJB:IJE,IKB-1) = -999.
!
!------------------------------------------------------------------------------
!
@@ -428,7 +435,10 @@ ZDZM(IKB-1) = -999.
PDXHATM =0.
! . local sum
IF (LCARTESIAN) THEN
- PDXHATM = SUM_1DFIELD_ll ( PDXHAT,'X',IIB_ll,IIE_ll,IINFO_ll)
+ PDXHATM = SUM_1DFIELD_ll ( PDXHAT,'X',IIB_ll,IIE_ll,IINFO_ll)
+ DO JJ=1,SIZE(PDXATH_ZS,2)
+ PDXATH_ZS(:,JJ) = PDXHAT(:)
+ END DO
ELSE
! Extraction of x-slice PMAP at j=(IJB_ll+IJE_ll)/2
CALL GET_SLICE_ll (PMAP,'X',(IJB_ll+IJE_ll)/2,ZXMAP(IIB:IIE) &
@@ -436,6 +446,9 @@ ELSE
! initialize the work array = PDXHAT/ZXMAP
ZWORKX(IIB:IIE) = PDXHAT(IIB:IIE)/ ZXMAP (IIB:IIE)
PDXHATM = SUM_1DFIELD_ll ( ZWORKX,'X',IIB_ll,IIE_ll,IINFO_ll)
+ DO JJ=1,SIZE(PDXATH_ZS,2)
+ PDXATH_ZS(:,JJ) = PDXHAT(:) / PMAP(:,JJ)
+ END DO
END IF
! . division to complete sum
PDXHATM = PDXHATM / REAL(IIMAX_ll)
@@ -448,6 +461,9 @@ PDXHATM = PDXHATM / REAL(IIMAX_ll)
PDYHATM = 0.
IF (LCARTESIAN) THEN
PDYHATM = SUM_1DFIELD_ll ( PDYHAT,'Y',IJB_ll,IJE_ll,IINFO_ll)
+ DO JI=1,SIZE(PDYATH_ZS,1)
+ PDYATH_ZS(JI,:) = PDYHAT(:)
+ END DO
ELSE
! Extraction of y-slice PMAP at i=IIB_ll+IIE_ll/2
CALL GET_SLICE_ll (PMAP,'Y',(IIB_ll+IIE_ll)/2,ZYMAP(IJB:IJE) &
@@ -455,6 +471,9 @@ ELSE
! initialize the work array = PDYHAT / ZYMAP
ZWORKY(IJB:IJE) = PDYHAT(IJB:IJE) / ZYMAP (IJB:IJE)
PDYHATM = SUM_1DFIELD_ll ( ZWORKY,'Y',IJB_ll,IJE_ll,IINFO_ll)
+ DO JI=1,SIZE(PDYATH_ZS,1)
+ PDYATH_ZS(JI,:) = PDYHAT(:) / PMAP(JI,:)
+ END DO
END IF
! . division to complete sum
PDYHATM= PDYHATM / REAL(IJMAX_ll)
@@ -464,9 +483,16 @@ PDYHATM= PDYHATM / REAL(IJMAX_ll)
!* 6. COMPUTE THE OUT-DIAGONAL ELEMENTS A AND C OF THE MATRIX
! -------------------------------------------------------
!
+PAF_ZS = 1.0
+PCF_ZS = 1.0
DO JK = IKB,IKE
PAF(JK) = 0.5 * ( PRHOM(JK-1) + PRHOM(JK) ) / ZDZM(JK) **2
PCF(JK) = 0.5 * ( PRHOM(JK) + PRHOM(JK+1) ) / ZDZM(JK+1) **2
+
+ PAF_ZS(IIB:IIE,IJB:IJE,JK) = 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,JK-1) + PRHO_ZS(IIB:IIE,IJB:IJE,JK) ) &
+ / ZDZM_ZS(IIB:IIE,IJB:IJE,JK) **2
+ PCF_ZS(IIB:IIE,IJB:IJE,JK) = 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,JK) + PRHO_ZS(IIB:IIE,IJB:IJE,JK+1) ) &
+ / ZDZM_ZS(IIB:IIE,IJB:IJE,JK+1) **2
END DO
!
! at the upper and lower levels PAF and PCF are computed using the Neumann
@@ -477,6 +503,14 @@ PCF(IKB-1) = 0.5 * ( PRHOM(IKB-1) + PRHOM(IKB) ) / ZDZM(IKB) **2
!
PAF(IKB-1) = 0.0 ! 0.5 * ( PRHOM(IKB-1) + PRHOM(IKB) ) / ZDZM(IKB) **2
PCF(IKE+1) = 0.0 ! 0.5 * ( PRHOM(IKE) + PRHOM(IKE+1) ) / ZDZM(IKE+1) **2
+
+PAF_ZS(IIB:IIE,IJB:IJE,IKE+1) = -0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKE) + PRHO_ZS(IIB:IIE,IJB:IJE,IKE+1) ) &
+ / ZDZM_ZS(IIB:IIE,IJB:IJE,IKE+1) **2
+PCF_ZS(IIB:IIE,IJB:IJE,IKB-1) = 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKB-1) + PRHO_ZS(IIB:IIE,IJB:IJE,IKB) ) &
+ / ZDZM_ZS(IIB:IIE,IJB:IJE,IKB) **2
+
+PAF_ZS(IIB:IIE,IJB:IJE,IKB-1) = 0.0
+PCF_ZS(IIB:IIE,IJB:IJE,IKE+1) = 0.0
!------------------------------------------------------------------------------
!* 7. COMPUTE THE DIAGONAL ELEMENTS B OF THE MATRIX
! ---------------------------------------------
@@ -619,15 +653,22 @@ ELSE
!JUAN Z_SPLITTING
!JUAN for Z splitting we need to do the vertical substitution
!JUAN in Bsplitting slides so need PBFB
+PBF_ZS = 1.0
DO JK= IKB,IKE
!!$ DO JJ= IJB-1,IJE+1
!!$ DO JI= IIB-1,IIE+1
!!$ PBFB(JI,JJ,JK) = PRHOM(JK)* ZEIGEN_ll(JI+IORXB_ll-IIB_ll,JJ+IORYB_ll-IJB_ll) - 0.5 * &
DO JJ= IJB,IJE
DO JI= IIB,IIE
+
PBFB(JI,JJ,JK) = PRHOM(JK)* ( -2.0 / PDXHATM**2 -2.0 /PDYHATM**2 ) - 0.5 * &
- ( ( PRHOM(JK-1) + PRHOM(JK) ) / ZDZM(JK) **2 &
+ ( ( PRHOM(JK-1) + PRHOM(JK) ) / ZDZM(JK) **2 &
+( PRHOM(JK) + PRHOM(JK+1) ) / ZDZM(JK+1)**2 )
+
+ PBF_ZS(JI,JJ,JK) = PRHO_ZS(JI,JJ,JK)* ( -2.0 / PDXATH_ZS(JI,JJ)**2 -2.0 /PDYATH_ZS(JI,JJ)**2 ) - 0.5 * &
+ ( ( PRHO_ZS(JI,JJ,JK-1) + PRHO_ZS(JI,JJ,JK) ) / ZDZM_ZS(JI,JJ,JK) **2 &
+ +( PRHO_ZS(JI,JJ,JK) + PRHO_ZS(JI,JJ,JK+1) ) / ZDZM_ZS(JI,JJ,JK+1)**2 )
+
END DO
END DO
END DO
@@ -640,6 +681,12 @@ ELSE
PBFB(IIB:IIE,IJB:IJE,IKE+1) = 0.5 * ( PRHOM(IKE) + PRHOM(IKE+1) ) / &
ZDZM(IKE+1) **2
+ PBF_ZS(IIB:IIE,IJB:IJE,IKB-1) = - 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKB-1) + PRHO_ZS(IIB:IIE,IJB:IJE,IKB) ) &
+ / ZDZM_ZS(IIB:IIE,IJB:IJE,IKB)**2
+
+ PBF_ZS(IIB:IIE,IJB:IJE,IKE+1) = 0.5 * ( PRHO_ZS(IIB:IIE,IJB:IJE,IKE) + PRHO_ZS(IIB:IIE,IJB:IJE,IKE+1) ) &
+ / ZDZM_ZS(IIB:IIE,IJB:IJE,IKE+1)**2
+
!!$ PBFB(IIB-1:IIE+1,IJB-1:IJE+1,IKB-1) = -0.5 * ( PRHOM(IKB-1) + PRHOM(IKB) ) / &
!!$ ZDZM(IKB) **2
!!$ !
diff --git a/zconjgrad.f90 b/zconjgrad.f90
index c3e977aae..35d82862f 100644
--- a/zconjgrad.f90
+++ b/zconjgrad.f90
@@ -242,7 +242,7 @@ DO JM = 1,KITR
!
ZWORK = PY - ZWORK ! Y - Q (PHI)
!
- CALL ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &! -1
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &! -1
PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZDELTA) ! F (Y - Q (PHI)))
!
!* 2.3 compute the auxiliary field P
@@ -255,7 +255,7 @@ DO JM = 1,KITR
ELSE
ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY, &
PDZZ,PRHODJ,PTHETAV,ZDELTA) ! Q ( DELTA )
- CALL ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! -1
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! -1
PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZWORKD) ! F * Q ( DELTA )
!
ZALPHA = - DOTPROD(ZWORKD,ZWORKP,HLBCX,HLBCY)/ZDOTPP ! ZWORKP,ZDOTPP come
@@ -269,7 +269,7 @@ DO JM = 1,KITR
!
ZWORK = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,&
PDZZ,PRHODJ,PTHETAV,ZP) ! Q ( P )
- CALL ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,& ! -1
+ CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,& ! -1
PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZWORK,ZWORKP) ! F * Q ( P )
!
! store the scalar product to compute lambda and next P
diff --git a/zsolver.f90 b/zsolver.f90
index f218af070..943eff22d 100644
--- a/zsolver.f90
+++ b/zsolver.f90
@@ -15,8 +15,10 @@
INTERFACE
!
SUBROUTINE ZSOLVER(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV, &
- PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
- KITR,PY,PPHI)
+ PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ KITR,PY,PPHI, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS) !JUAN FULL ZSOLVER
!
IMPLICIT NONE
!
@@ -61,6 +63,10 @@ INTEGER, INTENT(IN) :: KITR ! number of iterations for the
REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(IN) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_ZS
!
END SUBROUTINE ZSOLVER
!
@@ -72,8 +78,10 @@ END MODULE MODI_ZSOLVER
!
! #########################################################################
SUBROUTINE ZSOLVER(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV, &
- PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
- KITR,PY,PPHI)
+ PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY, &
+ KITR,PY,PPHI, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS)
! #########################################################################
!
!!**** *CONRESOL * - solve an elliptic equation by the conjugate residual
@@ -188,6 +196,10 @@ INTEGER, INTENT(IN) :: KITR ! number of iterations for the
REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
!
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHI ! solution of the equation
+
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(IN) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_ZS
!
!* 0.2 declarations of local variables
!
@@ -217,7 +229,9 @@ ZRESIDUE = QLAP(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,PTHETAV,PPHI) - PY
!* 1.2 compute the vector: p^(0) = F^(-1)*( Q(PHI) - Y )
!
CALL ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
- PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZP)
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZP, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS)
!
!* 1.3 compute the vector: delta^(0) = Q ( p^(0) )
!
@@ -250,7 +264,9 @@ DO JM = 1,KITR
!* 2.4 compute the vector: q = F^(-1)*( Q(PHI) - Y )
!
CALL ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
- PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZQ)
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZRESIDUE,ZQ, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS)
!
!* 2.5 compute the auxiliary field: ksi = Q ( q )
!
diff --git a/zsolver_inv.f90 b/zsolver_inv.f90
index f736a409c..ad850b8e5 100644
--- a/zsolver_inv.f90
+++ b/zsolver_inv.f90
@@ -15,7 +15,9 @@
INTERFACE
!
SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
- PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PY,PF_1_Y)
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PY,PF_1_Y, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS)
!
!
IMPLICIT NONE
@@ -49,14 +51,21 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
!
REAL, DIMENSION(:,:,:), INTENT(OUT):: PF_1_Y ! solution of the equation
!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(IN) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_ZS
+!
+!
END SUBROUTINE ZSOLVER_INV
!
END INTERFACE
!
END MODULE MODI_ZSOLVER_INV
! ######################################################################
- SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
- PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PY,PF_1_Y)
+ SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
+ PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PY,PF_1_Y, &
+ PAF_ZS,PBF_ZS,PCF_ZS, &
+ PDXATH_ZS,PDYATH_ZS,PRHO_ZS)
! ######################################################################
!
!!**** *FLAT_INV * - Invert the flat quasi-laplacian operator
@@ -168,6 +177,10 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PY ! RHS of the equation
!
REAL, DIMENSION(:,:,:), INTENT(OUT):: PF_1_Y ! solution of the equation
!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PAF_ZS,PBF_ZS,PCF_ZS
+REAL, DIMENSION(:,:) , INTENT(IN) :: PDXATH_ZS,PDYATH_ZS
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_ZS
+!
!* 0.2 declaration of local variables
!
REAL, DIMENSION(SIZE(PY,1),SIZE(PY,2),SIZE(PY,3)) :: ZY ! work array to store
@@ -229,6 +242,7 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZBAND_YRT ! array in Y slices distributio
!
INTEGER :: IIU,IJU
REAL :: ZMEAN
+INTEGER :: IT , NT = 1
!-------------------------------------------------------------------------------
!
!* 1. COMPUTE LOOP BOUNDS
@@ -346,21 +360,16 @@ IF (L2D) THEN
CALL FAST_SUBSTITUTION_2D(ZBAND_YR,ZBETX,PBF,ZGAM,PCF,ZAF &
,ZBAND_Y,IIY,IJY,IKU)
ELSE
+
CALL FAST_SUBSTITUTION_3D(PF_1_Y,ZBETX,PBF,ZGAM,PCF,PAF &
,ZY,IIU,IJU,IKU)
-END IF
-!
ZY = 0.0
DO JK=IKB,IKE
ZY(IIB:IIE,IJB:IJE,JK) = PAF(JK) * PF_1_Y(IIB:IIE,IJB:IJE,JK-1) &
+ ( PBF(IIB:IIE,IJB:IJE,JK) ) * PF_1_Y(IIB:IIE,IJB:IJE,JK) &
+ PCF(JK ) * PF_1_Y(IIB:IIE,IJB:IJE,JK+1) &
- - PY(IIB:IIE,IJB:IJE,JK)
-
-!!$ - ( PAF(JK) + PCF(JK) ) * PF_1_Y(IIB:IIE,IJB:IJE,JK) &
-
-
+ - PY(IIB:IIE,IJB:IJE,JK)
END DO
ZY(IIB:IIE,IJB:IJE,IKB-1) = - ( + PCF(IKB-1) ) * PF_1_Y(IIB:IIE,IJB:IJE,IKB-1) &
+ PCF(IKB-1) * PF_1_Y(IIB:IIE,IJB:IJE,IKB) &
@@ -370,9 +379,81 @@ ZY(IIB:IIE,IJB:IJE,IKE+1) = PAF(IKE+1) * PF_1_Y(IIB:IIE,IJB
- ( PAF(IKE+1) ) * PF_1_Y(IIB:IIE,IJB:IJE,IKE+1) &
- PY(IIB:IIE,IJB:IJE,IKE+1)
-!!$ZMEAN = SUM( PF_1_Y(IIB:IIE,IJB:IJE,IKE+1)) / FLOAT ( (IIE-IIB)*(IJE-IJB) )
-!!$PF_1_Y(:,:,:) = PF_1_Y(:,:,:) - ZMEAN
+ZDXM2 = PDXHATM*PDXHATM
+ZDYM2 = PDYHATM*PDYHATM
+DO JK=IKB,IKE
+ PF_1_Y(IIB-1,IJB:IJE,JK) = PF_1_Y(IIB,IJB:IJE,JK) - PY(IIB-1,IJB:IJE,JK)*ZDXM2/PRHOM(JK)
+ PF_1_Y(IIE+1,IJB:IJE,JK) = PF_1_Y(IIE,IJB:IJE,JK) + PY(IIE+1,IJB:IJE,JK)*ZDXM2/PRHOM(JK)
+ PF_1_Y(IIB:IIE,IJB-1,JK) = PF_1_Y(IIB:IIE,IJB,JK) - PY(IIB:IIE,IJB-1,JK)*ZDYM2/PRHOM(JK)
+ PF_1_Y(IIB:IIE,IJE+1,JK) = PF_1_Y(IIB:IIE,IJE,JK) + PY(IIB:IIE,IJE+1,JK)*ZDYM2/PRHOM(JK)
+
+END DO
+PF_1_Y(IIB-1,IJB:IJE,IKE+1) = PF_1_Y(IIB,IJB:IJE,IKE+1)
+PF_1_Y(IIE+1,IJB:IJE,IKB-1) = PF_1_Y(IIE,IJB:IJE,IKB-1)
+PF_1_Y(IIE+1,IJB:IJE,IKE+1) = PF_1_Y(IIE,IJB:IJE,IKE+1)
+PF_1_Y(IIB-1,IJB:IJE,IKB-1) = PF_1_Y(IIB,IJB:IJE,IKB-1)
+
+PF_1_Y(IIB:IIE,IJB-1,IKB-1) = PF_1_Y(IIB:IIE,IJB,IKB-1)
+PF_1_Y(IIB:IIE,IJB-1,IKE+1) = PF_1_Y(IIB:IIE,IJB,IKE+1)
+PF_1_Y(IIB:IIE,IJE+1,IKB-1) = PF_1_Y(IIB:IIE,IJE,IKB-1)
+PF_1_Y(IIB:IIE,IJE+1,IKE+1) = PF_1_Y(IIB:IIE,IJE,IKE+1)
+
+ZY = 0.0
+DO JK=IKB,IKE
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ ZY(JI,JJ,JK) = &
+!!$ PAF(JK) * PF_1_Y(JI,JJ,JK-1) &
+!!$ + ( PBF(JI,JJ,JK) ) * PF_1_Y(JI,JJ,JK) &
+!!$ + PCF(JK ) * PF_1_Y(JI,JJ,JK+1) + &
+ PRHOM(JK)* ( ( PF_1_Y(JI+1,JJ,JK) + PF_1_Y(JI-1,JJ,JK) ) / (PDXHATM*PDXHATM) &
+ + ( PF_1_Y(JI,JJ+1,JK) + PF_1_Y(JI,JJ-1,JK) ) / (PDYHATM*PDYHATM) )
+!!$ - PY(JI,JJ,JK)
+ END DO
+ END DO
+END DO
+DO IT=1,NT
+ZY = 0.0
+DO JK=IKB,IKE
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ ZY(JI,JJ,JK) = ( PY(JI,JJ,JK) &
+ - PAF(JK) * PF_1_Y(JI,JJ,JK-1) &
+ - PCF(JK) * PF_1_Y(JI,JJ,JK+1) &
+ - PRHOM(JK)* ( ( PF_1_Y(JI+1,JJ,JK) + PF_1_Y(JI-1,JJ,JK) ) / (PDXHATM*PDXHATM) &
+ + ( PF_1_Y(JI,JJ+1,JK) + PF_1_Y(JI,JJ-1,JK) ) / (PDYHATM*PDYHATM) ) ) &
+ / PBF(JI,JJ,JK)
+ END DO
+ END DO
+END DO
+
+PF_1_Y = ZY
+
+END DO
+
+!!$
+!!$ CALL FAST_SUBSTITUTION_3D_ZS(PF_1_Y(IIB:IIE,IJB:IJE,:),ZBETX,PBF_ZS(IIB:IIE,IJB:IJE,:),&
+!!$ ZGAM,PCF_ZS(IIB:IIE,IJB:IJE,:),PAF_ZS(IIB:IIE,IJB:IJE,:)&
+!!$ ,ZY(IIB:IIE,IJB:IJE,:),(IIE-IIB+1),(IJE-IJB+1),IKU)
+!!$
+!!$ZY = 0.0
+!!$DO JK=IKB,IKE
+!!$ ZY(IIB:IIE,IJB:IJE,JK) = PAF_ZS(IIB:IIE,IJB:IJE,JK) * PF_1_Y(IIB:IIE,IJB:IJE,JK-1) &
+!!$ + ( PBF_ZS(IIB:IIE,IJB:IJE,JK) ) * PF_1_Y(IIB:IIE,IJB:IJE,JK) &
+!!$ + PCF_ZS(IIB:IIE,IJB:IJE,JK ) * PF_1_Y(IIB:IIE,IJB:IJE,JK+1) &
+!!$ - PY(IIB:IIE,IJB:IJE,JK)
+!!$END DO
+!!$ZY(IIB:IIE,IJB:IJE,IKB-1) = - ( + PCF_ZS(IIB:IIE,IJB:IJE,IKB-1) ) * PF_1_Y(IIB:IIE,IJB:IJE,IKB-1) &
+!!$ + PCF_ZS(IIB:IIE,IJB:IJE,IKB-1) * PF_1_Y(IIB:IIE,IJB:IJE,IKB) &
+!!$ - PY(IIB:IIE,IJB:IJE,IKB-1)
+!!$
+!!$ZY(IIB:IIE,IJB:IJE,IKE+1) = PAF_ZS(IIB:IIE,IJB:IJE,IKE+1) * PF_1_Y(IIB:IIE,IJB:IJE,IKE) &
+!!$ - ( PAF_ZS(IIB:IIE,IJB:IJE,IKE+1) ) * PF_1_Y(IIB:IIE,IJB:IJE,IKE+1) &
+!!$ - PY(IIB:IIE,IJB:IJE,IKE+1)
+
+END IF
+!
!
!-------------------------------------------------------------------------------
!
@@ -513,6 +594,45 @@ CALL GET_HALO(PF_1_Y)
!
CONTAINS
+SUBROUTINE FAST_SUBSTITUTION_3D_ZS(PBAND_YR,PBETX,PPBF,PGAM,PPCF,PAF &
+ ,PBAND_Y,KIY,KJY,KKU)
+INTEGER :: KIY,KJY,KKU
+REAL, DIMENSION (KIY*KJY,KKU) :: PBAND_YR,PBAND_Y,PPBF,PGAM
+REAL, DIMENSION (KIY*KJY) :: PBETX
+REAL, DIMENSION (KIY*KJY,KKU) :: PPCF,PAF
+INTEGER :: JK
+!
+!
+! initialization
+!
+!
+PBAND_YR = 0.0
+PBETX(:) = PPBF(:,IKB-1)
+
+PBAND_YR(:,IKB-1) = PBAND_Y(:,IKB-1) &
+ / PBETX(:)
+!
+! 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
+!
+!
+END SUBROUTINE FAST_SUBSTITUTION_3D_ZS
SUBROUTINE FAST_SUBSTITUTION_3D(PBAND_YR,PBETX,PPBF,PGAM,PPCF,PAF &
,PBAND_Y,KIY,KJY,KKU)
@@ -527,11 +647,7 @@ INTEGER :: JK
!
!
PBAND_YR = 0.0
-!PAF(IKB-1) = 0.0 ! - PAF(IKB) ! 0.0
-!PPCF(IKE+1) = - PPCF(IKE) ! - PPCF(IKE)
-
PBETX(:) = PPBF(:,IKB-1)
-!!$PBETX(:) = - ( PPCF(IKB-1) + PAF(IKB-1) )
PBAND_YR(:,IKB-1) = PBAND_Y(:,IKB-1) &
/ PBETX(:)
@@ -541,15 +657,10 @@ PBAND_YR(:,IKB-1) = PBAND_Y(:,IKB-1) &
DO JK = IKB,IKE+1
PGAM(:,JK) = PPCF(JK-1) / PBETX(:)
-
!
-
PBETX(:) = PPBF(:,JK) - PAF(JK)*PGAM(:,JK)
-!!$ PBETX(:) = - ( PPCF(JK) + PAF(JK) ) - PAF(JK)*PGAM(:,JK)
!
-
PBAND_YR(:,JK) = ( PBAND_Y(:,JK) - PAF(JK)*PBAND_YR(:,JK- 1) ) /PBETX(:)
-
!
END DO
!
--
GitLab