From 169770cbe796b3c9a1d226137846c831e322f499 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Thu, 10 Jan 2019 12:14:10 +0100
Subject: [PATCH] Quentin 10/01/2019: add a local filtering for very high
slopes in PGD step (for a better chance of pressure solver convergence). In
namelist NAM_ZSFILTER : possibility to activate or not (LHSLOP), choose the
number of filtering iterations (NLOCZSFILTER) and the slope value (XHSLOP=
Dz/Dx and Dz/Dy) where the filter is applied
---
src/MNH/prep_pgd.f90 | 11 ++-
src/MNH/zsmt_pgd.f90 | 205 +++++++++++++++++++++++++++++++++++++++----
2 files changed, 196 insertions(+), 20 deletions(-)
diff --git a/src/MNH/prep_pgd.f90 b/src/MNH/prep_pgd.f90
index 45a14fcf0..001b1e5c8 100644
--- a/src/MNH/prep_pgd.f90
+++ b/src/MNH/prep_pgd.f90
@@ -72,6 +72,8 @@
!! 10/2016 (S.Faroux S.Bielli) correction for NHALO=0
!! 01/2018 (G.Delautier) SURFEX 8.1
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Q. Rodier 01/2019 : add a new filtering for very high slopes in NAM_ZSFILTER
+!!
!----------------------------------------------------------------------------
!
!* 0. DECLARATION
@@ -132,8 +134,9 @@ CHARACTER(LEN=28) :: YDAD =' ' ! name of dad of input FM file
CHARACTER(LEN=28) :: CPGDFILE ='PGDFILE' ! name of the output file
CHARACTER(LEN=100) :: YMSG
INTEGER :: NZSFILTER=1 ! number of iteration for filter for fine orography
-LOGICAL :: LHSLOP=.FALSE. ! filtering of slopes higher than XHSLOP
-REAL :: XHSLOP=1.2 ! if LHSLOP filtering of slopes higher than XHSLOP
+INTEGER :: NLOCZSFILTER=3 ! number of iteration for filter of local fine orography
+LOGICAL :: LHSLOP=.FALSE. ! filtering of local slopes higher than XHSLOP
+REAL :: XHSLOP=1.0 ! slopes where the local fine filtering is applied
INTEGER :: NSLEVE =12 ! number of iteration for filter for smooth orography
REAL :: XSMOOTH_ZS = XUNDEF ! optional uniform smooth orography for SLEVE coordinate
REAL, DIMENSION(:,:),ALLOCATABLE :: ZWORK ! work array for lat and lon reshape
@@ -145,7 +148,7 @@ TYPE(TFILEDATA),POINTER :: TZFILE => NULL()
TYPE(TFILEDATA),POINTER :: TZNMLFILE => NULL() ! Namelist file
!
NAMELIST/NAM_PGDFILE/CPGDFILE, NHALO
-NAMELIST/NAM_ZSFILTER/NZSFILTER,LHSLOP,XHSLOP
+NAMELIST/NAM_ZSFILTER/NZSFILTER,NLOCZSFILTER,LHSLOP,XHSLOP
NAMELIST/NAM_SLEVE/NSLEVE, XSMOOTH_ZS
NAMELIST/NAM_CONF_PGD/JPHEXT, NHALO_MNH
!------------------------------------------------------------------------------
@@ -279,7 +282,7 @@ NULLIFY(TFILE_SURFEX) !Probably not necessary
!
!* 4. Computes and writes smooth orography for SLEVE coordinate
! ---------------------------------------------------------
-CALL ZSMT_PGD(TZFILE,NZSFILTER,NSLEVE,XSMOOTH_ZS)
+CALL ZSMT_PGD(TZFILE,NZSFILTER,NSLEVE,NLOCZSFILTER,LHSLOP,XHSLOP,XSMOOTH_ZS)
!
IF (.NOT.LCARTESIAN) THEN
!!!! WRITE LAT and LON
diff --git a/src/MNH/zsmt_pgd.f90 b/src/MNH/zsmt_pgd.f90
index 0de2cd8a2..14a4be23b 100644
--- a/src/MNH/zsmt_pgd.f90
+++ b/src/MNH/zsmt_pgd.f90
@@ -8,14 +8,17 @@
!
INTERFACE
!
- SUBROUTINE ZSMT_PGD(TPFILE,KZSFILTER,KSLEVE,PSMOOTH_ZS)
+ SUBROUTINE ZSMT_PGD(TPFILE,KZSFILTER,KSLEVE,KLOCZSFILTER,OHSLOP,PHSLOP,PSMOOTH_ZS)
!
USE MODD_IO_ll, ONLY : TFILEDATA
!
-TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! File characteristics
-INTEGER, INTENT(IN) :: KZSFILTER ! number of iterations for fine orography
-INTEGER, INTENT(IN) :: KSLEVE ! number of iterations
-REAL, INTENT(IN) :: PSMOOTH_ZS ! optional uniform smooth orography for SLEVE coordinate
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! File characteristics
+INTEGER, INTENT(IN) :: KZSFILTER ! number of iterations for fine orography
+INTEGER, INTENT(IN) :: KSLEVE ! number of iterations
+INTEGER, INTENT(IN) :: KLOCZSFILTER ! number of iteration for filter of local fine orography
+LOGICAL, INTENT(IN) :: OHSLOP ! filtering of local slopes higher than XHSLOP
+REAL, INTENT(IN) :: PHSLOP ! slopes where the local fine filtering is applied
+REAL, INTENT(IN) :: PSMOOTH_ZS ! optional uniform smooth orography for SLEVE coordinate
!
END SUBROUTINE ZSMT_PGD
!
@@ -26,7 +29,7 @@ END MODULE MODI_ZSMT_PGD
!
!
! #############################
- SUBROUTINE ZSMT_PGD(TPFILE,KZSFILTER,KSLEVE,PSMOOTH_ZS)
+ SUBROUTINE ZSMT_PGD(TPFILE,KZSFILTER,KSLEVE,KLOCZSFILTER,OHSLOP,PHSLOP,PSMOOTH_ZS)
! #############################
!
!!**** *ZSMT_PGD* computes smoothed orography for SLEVE coordinate
@@ -56,7 +59,8 @@ END MODULE MODI_ZSMT_PGD
!! -------------
!! Original nov 2005
!! J.Escobar 23/06/2015 : correction for JPHEXT<>1
-!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Q. Rodier 01/2019 : add a new filtering for very high slopes (applied locally)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -70,15 +74,20 @@ USE MODE_FMREAD
USE MODE_FMWRIT
USE MODE_ll , ONLY : GET_DIM_EXT_ll , ADD2DFIELD_ll , CLEANLIST_ll , UPDATE_HALO_ll
USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+USE MODI_SUM_ll
+USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL
!
IMPLICIT NONE
!
!* 0.1 declarations of arguments
!
-TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! File characteristics
-INTEGER, INTENT(IN) :: KZSFILTER ! number of iterations for fine orography
-INTEGER, INTENT(IN) :: KSLEVE ! number of iterations
-REAL, INTENT(IN) :: PSMOOTH_ZS ! optional uniform smooth orography for SLEVE coordinate
+TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! File characteristics
+INTEGER, INTENT(IN) :: KZSFILTER ! number of iterations for fine orography
+INTEGER, INTENT(IN) :: KSLEVE ! number of iterations
+INTEGER, INTENT(IN) :: KLOCZSFILTER ! number of iteration for filter of local fine orography
+LOGICAL, INTENT(IN) :: OHSLOP ! filtering of local slopes higher than XHSLOP
+REAL, INTENT(IN) :: PHSLOP ! slopes where the local fine filtering is applied
+REAL, INTENT(IN) :: PSMOOTH_ZS ! optional uniform smooth orography for SLEVE coordinate
!
!
!* 0.2 declarations of local variables
@@ -88,8 +97,6 @@ INTEGER :: JN ! loop counter on iterations
INTEGER :: JI ! loop counter on X coordinate
INTEGER :: JJ ! loop counter on Y coordinate
!
-INTEGER :: IIMAX ! number of physical points in X direction
-INTEGER :: IJMAX ! number of physical points in Y direction
INTEGER :: IIU ! number of points in X direction
INTEGER :: IJU ! number of points in Y direction
!
@@ -98,11 +105,22 @@ REAL, DIMENSION(:,:), ALLOCATABLE :: ZFINE_ZS ! smoothed fine orography
REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLEVE_ZS ! smooth orography for sleve coordinate
REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS ! orography at previous iteration
REAL, DIMENSION(:,:), ALLOCATABLE :: ZMASK ! sea mask
+!
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLOPEX ! terrain slope along x (flux pt)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLOPEY ! terrain slope along y (flux pt)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLOPEMX ! terrain slope along x (max. of flux, at mass pt)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSLOPEMY ! terrain slope along y (max. of flux, at mass pt)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCOFIL ! coefficient for filtering high slopes
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSMOOTH_ZSINI ! initial orography before slope filtering
+!
INTEGER :: JIM, JIP, JJM, JJP
TYPE(LIST_ll) , POINTER :: THALO_ll => NULL() ! halo
INTEGER :: INFO_ll ! error return code
INTEGER :: IIB,IIE,IJB,IJE
+REAL, DIMENSION(:), ALLOCATABLE :: ZXHAT
+REAL, DIMENSION(:), ALLOCATABLE :: ZYHAT
+TYPE(TFIELDDATA) :: TZFIELD
!-------------------------------------------------------------------------------
!
!* 1. Read orography in the file
@@ -123,7 +141,18 @@ ALLOCATE(ZSMOOTH_ZS(IIU,IJU))
ALLOCATE(ZFINE_ZS(IIU,IJU))
ALLOCATE(ZSLEVE_ZS(IIU,IJU))
ALLOCATE(ZMASK(IIU,IJU))
-!
+ALLOCATE(ZSLOPEX(IIU,IJU))
+ALLOCATE(ZSLOPEY(IIU,IJU))
+ALLOCATE(ZSLOPEMX(IIU,IJU))
+ALLOCATE(ZSLOPEMY(IIU,IJU))
+ALLOCATE(ZCOFIL(IIU,IJU))
+ALLOCATE(ZSMOOTH_ZSINI(IIU,IJU))
+ALLOCATE(ZXHAT(IIU))
+ALLOCATE(ZYHAT(IJU))
+!
+CALL IO_READ_FIELD(TPFILE,'XHAT',ZXHAT)
+CALL IO_READ_FIELD(TPFILE,'YHAT',ZYHAT)
+
!PW: bug/TODO: read a field in a file opened in WRITE mode
!There is a test in IO_READ_FIELD_BYFIELD_X2 to allow this if TPFILE%CMODE='LFICDF4'
CALL IO_READ_FIELD(TPFILE,'ZS',ZZS)
@@ -162,6 +191,8 @@ ZSMOOTH_ZS = ZZS
!
CALL ADD2DFIELD_ll(THALO_ll,ZZS)
CALL ADD2DFIELD_ll(THALO_ll,ZSMOOTH_ZS)
+CALL ADD2DFIELD_ll(THALO_ll,ZSLOPEX)
+CALL ADD2DFIELD_ll(THALO_ll,ZSLOPEY)
!
CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
@@ -192,7 +223,6 @@ CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
IF (JN==KSLEVE) ZSLEVE_ZS= ZSMOOTH_ZS
ENDDO
!
-CALL CLEANLIST_ll(THALO_ll)
!-------------------------------------------------------------------------------
!
!* 3. Case of uniform smooth orography prescribed
@@ -205,7 +235,142 @@ IF (PSMOOTH_ZS /= XUNDEF) THEN
END IF
!-------------------------------------------------------------------------------
!
-!* 4. writes smoothed orographies in the file
+!* 4. Filtering very high slopes for steep orography
+! -------------------------------------------
+!
+IF(OHSLOP) THEN
+ ZSMOOTH_ZSINI=ZFINE_ZS
+ ZZS=ZFINE_ZS
+ ZSLOPEMX=0.
+ ZSLOPEMY=0.
+ ZSLOPEX=0.
+ ZSLOPEY=0.
+ !
+ DO JN=1, KLOCZSFILTER
+ ZCOFIL=0.
+ ZSLOPEX=0.
+ ZSLOPEY=0.
+ ZSLOPEMX=0.
+ ZSLOPEMY=0.
+ !
+ !Slope calculation along Y at flux and mass point
+ DO JI=1,IIU
+ DO JJ=2,IJU-1
+ ZSLOPEY(JI,JJ) = (ZZS(JI,JJ)-ZZS(JI,JJ-1))/(0.5*(ZYHAT(JJ+1)-ZYHAT(JJ-1)))
+ END DO
+ END DO
+ !
+ CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
+ !
+ DO JI=1,IIU
+ DO JJ=1,IJU-1
+ ZSLOPEMY(JI,JJ) = MAX(ABS(ZSLOPEY(JI,JJ)),ABS(ZSLOPEY(JI,JJ+1)))
+ END DO
+ END DO
+ !
+ !Slope calculation along X at flux and mass point
+ DO JJ=1,IJU
+ DO JI=2,IIU-1
+ ZSLOPEX(JI,JJ) = (ZZS(JI,JJ)-ZZS(JI-1,JJ))/(0.5*(ZXHAT(JI+1)-ZXHAT(JI-1)))
+ END DO
+ END DO
+ !
+ CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
+ !
+ DO JJ=1,IJU
+ DO JI=1,IIU-1
+ ZSLOPEMX(JI,JJ) = MAX(ABS(ZSLOPEX(JI+1,JJ)),ABS(ZSLOPEX(JI,JJ)))
+ END DO
+ END DO
+ !
+ !Filtering coefficient
+ DO JI=1,IIU
+ DO JJ=1,IJU
+ IF(ZSLOPEMX(JI,JJ)>=PHSLOP .OR. ZSLOPEMY(JI,JJ)>=PHSLOP) THEN
+ ZCOFIL(JI,JJ) = 0.25 !Fixed coefficient
+ END IF
+ END DO
+ END DO
+ !
+ !Filtering
+ DO JJ = 1,IJU
+ DO JI = 1,IIU
+ JIP = MIN(JI+1,IIU)
+ JIM = MAX(JI-1,1 )
+ JJP = MIN(JJ+1,IJU)
+ JJM = MAX(JJ-1,1 )
+ ZSMOOTH_ZS(JI,JJ) = ZZS(JI,JJ) &
+ + ZCOFIL(JI,JJ)* ZMASK(JI,JJ) &
+ * ( ZMASK(JIM,JJ) * ZZS(JIM,JJ) &
+ + ZMASK(JIP,JJ) * ZZS(JIP,JJ) &
+ + ZMASK(JI,JJM) * ZZS(JI,JJM) &
+ + ZMASK(JI,JJP) * ZZS(JI,JJP) &
+ - ( ZMASK(JIM,JJ) &
+ +ZMASK(JIP,JJ) &
+ +ZMASK(JI,JJM) &
+ +ZMASK(JI,JJP) ) * ZZS(JI,JJ) )
+ ENDDO
+ ENDDO
+ !
+ CALL UPDATE_HALO_ll(THALO_ll,INFO_ll)
+ ZZS=ZSMOOTH_ZS
+ END DO
+ !
+ ZFINE_ZS=ZSMOOTH_ZS
+ !
+ !Filtered slopes computed again for output
+ DO JI=1,IIU
+ DO JJ=2,IJU-1
+ ZSLOPEY(JI,JJ) = (ZZS(JI,JJ)-ZZS(JI,JJ-1))/(0.5*(ZYHAT(JJ+1)-ZYHAT(JJ-1)))
+ END DO
+ END DO
+ !
+ DO JJ=1,IJU
+ DO JI=2,IIU-1
+ ZSLOPEX(JI,JJ) = (ZZS(JI,JJ)-ZZS(JI-1,JJ))/(0.5*(ZXHAT(JI+1)-ZXHAT(JI-1)))
+ END DO
+ END DO
+ !
+ ! Writes filtred orography and slopes along i and j
+ TZFIELD%CMNHNAME = 'ZSLOPEX'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ZSLOPEX'
+ TZFIELD%CUNITS = ''
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'orography slope along x'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZSLOPEX)
+ !
+ TZFIELD%CMNHNAME = 'ZSLOPEY'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ZSLOPEY'
+ TZFIELD%CUNITS = ''
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'orography slope along y'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZSLOPEY)
+ !
+ TZFIELD%CMNHNAME = 'ZS_FILTR'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'ZS_FILTR'
+ TZFIELD%CUNITS = 'm'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'filtred orography'
+ TZFIELD%NGRID = 4
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .FALSE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZSMOOTH_ZSINI-ZFINE_ZS)
+END IF
+!-------------------------------------------------------------------------------
+!
+!* 5. writes smoothed orographies in the file
! ---------------------------------------
!
!
@@ -217,6 +382,14 @@ DEALLOCATE(ZFINE_ZS)
DEALLOCATE(ZSMOOTH_ZS)
DEALLOCATE(ZSLEVE_ZS)
DEALLOCATE(ZMASK)
+DEALLOCATE(ZSLOPEX)
+DEALLOCATE(ZSLOPEY)
+DEALLOCATE(ZSLOPEMX)
+DEALLOCATE(ZSLOPEMY)
+DEALLOCATE(ZCOFIL)
+DEALLOCATE(ZSMOOTH_ZSINI)
+!
+CALL CLEANLIST_ll(THALO_ll)
!
!-------------------------------------------------------------------------------
!
--
GitLab