From a90e04873986214515c8e4eea672b05fd312e1d1 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Wed, 7 Sep 2022 14:39:23 +0200
Subject: [PATCH] Philippe 07/09/2022: deduplication of PX/Y/ZHATM
 interpolations + update_halo for horizontal values

---
 src/MNH/ini_modeln.f90      |   2 +-
 src/MNH/ini_spectren.f90    |   2 +-
 src/MNH/prep_ideal_case.f90 |   7 +-
 src/MNH/read_hgridn.f90     |   7 +-
 src/MNH/read_ver_grid.f90   |   4 +-
 src/MNH/set_grid.f90        | 178 +++++++++++++++++++-----------------
 src/MNH/spawn_grid2.f90     |   9 +-
 src/MNH/write_surf_mnh.f90  |   9 +-
 8 files changed, 106 insertions(+), 112 deletions(-)

diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 282936099..9ef32c017 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -412,6 +412,7 @@ USE MODE_ll
 USE MODE_MODELN_HANDLER
 USE MODE_MPPDB
 USE MODE_MSG
+USE MODE_SET_GRID
 USE MODE_SPLITTINGZ_ll,     only: GET_DIM_EXTZ_ll
 USE MODE_TYPE_ZDIFFU
 
@@ -460,7 +461,6 @@ USE MODI_MNHGET_SURF_PARAM_n
 USE MODI_MNHREAD_ZS_DUMMY_n
 USE MODI_READ_FIELD
 USE MODI_SET_DIRCOS
-USE MODI_SET_GRID
 USE MODI_SET_REF
 #ifdef CPLOASIS
 USE MODI_SFX_OASIS_READ_NAM
diff --git a/src/MNH/ini_spectren.f90 b/src/MNH/ini_spectren.f90
index 7f62014dd..c17715ec6 100644
--- a/src/MNH/ini_spectren.f90
+++ b/src/MNH/ini_spectren.f90
@@ -105,6 +105,7 @@ USE MODE_IO_FIELD_READ, only: IO_Field_read
 USE MODE_ll
 USE MODE_MODELN_HANDLER
 USE MODE_MSG
+USE MODE_SET_GRID
 USE MODE_SPLITTINGZ_ll, ONLY: GET_DIM_EXTZ_ll
 USE MODE_TYPE_ZDIFFU
 !
@@ -114,7 +115,6 @@ USE MODI_INI_DYNAMICS
 USE MODI_INI_SPAWN_LS_n
 USE MODI_GET_SIZEX_LB
 USE MODI_GET_SIZEY_LB
-USE MODI_SET_GRID
 USE MODI_METRICS
 USE MODI_SET_REF
 USE MODI_UPDATE_METRICS
diff --git a/src/MNH/prep_ideal_case.f90 b/src/MNH/prep_ideal_case.f90
index d1604d2f1..6be3241c1 100644
--- a/src/MNH/prep_ideal_case.f90
+++ b/src/MNH/prep_ideal_case.f90
@@ -379,6 +379,7 @@ USE MODE_ll
 USE MODE_MODELN_HANDLER
 use mode_field,            only: Alloc_field_scalars, Ini_field_list, Ini_field_scalars
 USE MODE_MSG
+USE MODE_SET_GRID,         only: INTERP_HORGRID_TO_MASSPOINTS
 !
 USE MODI_DEFAULT_DESFM_n    ! Interface modules
 USE MODI_DEFAULT_EXPRE
@@ -1286,11 +1287,7 @@ ELSE
   END IF
 
   ! Interpolations of positions to mass points
-  XXHATM(1:NIU-1) = 0.5 * XXHAT(1:NIU-1) + 0.5 * XXHAT(2:NIU)
-  XXHATM(  NIU)   = 1.5 * XXHAT(  NIU)   - 0.5 * XXHAT(NIU-1)
-
-  XYHATM(1:NJU-1) = 0.5 * XYHAT(1:NJU-1) + 0.5 * XYHAT(2:NJU)
-  XYHATM(  NJU)   = 1.5 * XYHAT(  NJU)   - 0.5 * XYHAT(NJU-1)
+  CALL INTERP_HORGRID_TO_MASSPOINTS( XXHAT, XYHAT, XXHATM, XYHATM )
 END IF
 !
 !*       5.1.2  Orography and Gal-Chen Sommerville transformation :
diff --git a/src/MNH/read_hgridn.f90 b/src/MNH/read_hgridn.f90
index b60277ad2..39586f979 100644
--- a/src/MNH/read_hgridn.f90
+++ b/src/MNH/read_hgridn.f90
@@ -93,6 +93,7 @@ USE MODE_IO,            only: IO_Pack_set
 USE MODE_IO_FIELD_READ, only: IO_Field_read
 USE MODE_MSG
 USE MODE_MODELN_HANDLER
+USE MODE_SET_GRID,       only: INTERP_HORGRID_TO_MASSPOINTS
 use MODE_TOOLS_ll,       only: GET_DIM_EXT_ll, GET_DIM_PHYS_ll, GET_INDICE_ll
 !
 IMPLICIT NONE
@@ -255,11 +256,7 @@ IF ( .NOT. ASSOCIATED(XXHATM) ) ALLOCATE( XXHATM(SIZE( XXHAT )) )
 IF ( .NOT. ASSOCIATED(XYHATM) ) ALLOCATE( XYHATM(SIZE( XYHAT )) )
 
 ! Interpolations of positions to mass points
-XXHATM( : UBOUND(XXHATM,1)-1 ) = 0.5 * XXHAT( : UBOUND(XXHAT,1)-1 ) + 0.5 * XXHAT( LBOUND(XXHAT,1)+1 : UBOUND(XXHAT,1) )
-XXHATM( UBOUND(XXHATM,1)     ) = 1.5 * XXHAT( UBOUND(XXHAT,1)     ) - 0.5 * XXHAT( UBOUND(XXHAT,1)-1 )
-
-XYHATM( : UBOUND(XYHATM,1)-1 ) = 0.5 * XYHAT( : UBOUND(XYHAT,1)-1 ) + 0.5 * XYHAT( LBOUND(XYHAT,1)+1 : UBOUND(XYHAT,1) )
-XYHATM( UBOUND(XYHATM,1)     ) = 1.5 * XYHAT( UBOUND(XYHAT,1)     ) - 0.5 * XYHAT( UBOUND(XYHAT,1)-1 )
+CALL INTERP_HORGRID_TO_MASSPOINTS( XXHAT, XYHAT, XXHATM, XYHATM )
 
 !JUAN REALZ
 IF ( CPROGRAM .EQ. "REAL  " ) THEN
diff --git a/src/MNH/read_ver_grid.f90 b/src/MNH/read_ver_grid.f90
index bd1de11a9..653bbec02 100644
--- a/src/MNH/read_ver_grid.f90
+++ b/src/MNH/read_ver_grid.f90
@@ -111,6 +111,7 @@ USE MODD_PARAMETERS
 !
 USE MODE_MSG
 USE MODE_POS
+USE MODE_SET_GRID, ONLY: INTERP_VERGRID_TO_MASSPOINTS
 !
 USE MODI_DEFAULT_SLEVE
 !
@@ -328,8 +329,7 @@ END SELECT
 XZTOP = XZHAT(IKU)
 
 ! Interpolations of positions to mass points
-XZHATM(1:IKU-1 ) = 0.5 * XZHAT(1:IKU-1) + 0.5 * XZHAT(2:IKU  )
-XZHATM(  IKU   ) = 1.5 * XZHAT(  IKU  ) - 0.5 * XZHAT(  IKU-1)
+CALL INTERP_VERGRID_TO_MASSPOINTS( XZHAT, XZHATM )
 
 !-------------------------------------------------------------------------------
 !
diff --git a/src/MNH/set_grid.f90 b/src/MNH/set_grid.f90
index dd7daa678..1f972fd21 100644
--- a/src/MNH/set_grid.f90
+++ b/src/MNH/set_grid.f90
@@ -3,86 +3,24 @@
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
+! Modifications:
+!  P. Wautelet 31/08/2022: add PXHATM and PYHATM variables
+!  P. Wautelet 07/09/2022: add INTERP_HORGRID_1DIR_TO_MASSPOINTS, INTERP_HORGRID_TO_MASSPOINTS, INTERP_VERGRID_TO_MASSPOINTS
+!-----------------------------------------------------------------
 !     ####################
-      MODULE MODI_SET_GRID
+      MODULE MODE_SET_GRID
 !     ####################
-!
-INTERFACE
-!
-      SUBROUTINE SET_GRID( KMI, TPINIFILE,                                &
-                           KKU, KIMAX_ll, KJMAX_ll,                       &
-                           PTSTEP, PSEGLEN,                               &
-                           PLONORI, PLATORI, PLON, PLAT,                  &
-                           PXHAT, PYHAT, PDXHAT, PDYHAT, PXHATM, PYHATM,  &
-                           PMAP, PZS, PZZ, PZHAT, PZHATM, PZTOP, OSLEVE,  &
-                           PLEN1, PLEN2, PZSMT, PJ,                       &
-                           TPDTMOD, TPDTCUR, KSTOP,                       &
-                           KBAK_NUMB, KOUT_NUMB, TPBACKUPN, TPOUTPUTN     )
-!
-USE MODD_TYPE_DATE
-USE MODD_IO, ONLY: TFILEDATA,TOUTBAK
-!
-INTEGER,                INTENT(IN)  :: KMI       ! Model index
-TYPE(TFILEDATA),        INTENT(IN)  :: TPINIFILE !Initial file
-INTEGER,                INTENT(IN)  :: KKU       ! Upper dimension in z direction
-                                                 ! for domain arrays
-INTEGER,                INTENT(IN)  :: KIMAX_ll  !  Dimensions  in x direction
-                                                 ! of the physical domain,
-INTEGER,                INTENT(IN)  :: KJMAX_ll  !  Dimensions  in y direction
-                                                 ! of the physical domain,
-!
-REAL,                   INTENT(IN)  :: PTSTEP    ! time step of model KMI
-REAL,                   INTENT(INOUT) :: PSEGLEN ! segment duration (in seconds)
-!
-REAL,                   INTENT(OUT) :: PLONORI   ! Longitude  of the
-                                                 ! Origine point of
-                                                 ! conformal projection
-REAL,                   INTENT(OUT) :: PLATORI   ! Latitude of the
-                                                 ! Origine point of
-                                                 ! conformal projection
-REAL, DIMENSION(:,:),   INTENT(OUT) :: PLON,PLAT ! Longitude and latitude
-REAL, DIMENSION(:),     INTENT(OUT) :: PXHAT     ! Position x in the conformal
-                                                 ! plane or on the cartesian plane
-REAL, DIMENSION(:),     INTENT(OUT) :: PYHAT     ! Position y in the conformal
-                                                 ! plane or on the cartesian plane
-REAL, DIMENSION(:),     INTENT(OUT) :: PDXHAT    ! horizontal stretching in x
-REAL, DIMENSION(:),     INTENT(OUT) :: PDYHAT    ! horizontal stretching in y
-REAL, DIMENSION(:),     INTENT(OUT) :: PXHATM    ! Position x in the conformal plane or on the cartesian plane at mass points
-REAL, DIMENSION(:),     INTENT(OUT) :: PYHATM    ! Position y in the conformal plane or on the cartesian plane at mass points
-REAL, DIMENSION(:,:),   INTENT(OUT) :: PMAP      ! Map factor
-!
-REAL, DIMENSION(:,:),   INTENT(OUT) :: PZS       ! orography
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PZZ       ! Height z
-REAL, DIMENSION(:),     INTENT(OUT) :: PZHAT     ! Height level
-REAL, DIMENSION(:),     INTENT(OUT) :: PZHATM    ! Height level at mass points
-REAL,                   INTENT(OUT) :: PZTOP     ! Model top
-LOGICAL,                INTENT(OUT) :: OSLEVE    ! flag for SLEVE coordinate
-REAL,                   INTENT(OUT) :: PLEN1     ! Decay scale for smooth topography
-REAL,                   INTENT(OUT) :: PLEN2     ! Decay scale for small-scale topography deviation
-REAL, DIMENSION(:,:),   INTENT(OUT) :: PZSMT     ! smooth-orography
-!
-TYPE (DATE_TIME),       INTENT(OUT) :: TPDTMOD   ! date and time of the model
-                                                 ! beginning
-TYPE (DATE_TIME),       INTENT(OUT) :: TPDTCUR   ! Current date and time
-INTEGER,                INTENT(OUT) :: KSTOP     ! number of time steps for
-                                                 ! current segment
-INTEGER,POINTER,        INTENT(OUT) :: KBAK_NUMB ! number of backups
-INTEGER,POINTER,        INTENT(OUT) :: KOUT_NUMB ! number of outputs
-TYPE(TOUTBAK),DIMENSION(:),POINTER,INTENT(OUT) :: TPBACKUPN ! List of backups
-TYPE(TOUTBAK),DIMENSION(:),POINTER,INTENT(OUT) :: TPOUTPUTN ! List of outputs
-!
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PJ        ! Jacobian
-!
-END SUBROUTINE SET_GRID
-!
-END INTERFACE
-!
-END MODULE MODI_SET_GRID
-!
-!
-!
-!
-!
+
+IMPLICIT NONE
+
+PRIVATE
+
+PUBLIC :: SET_GRID
+
+PUBLIC :: INTERP_HORGRID_1DIR_TO_MASSPOINTS, INTERP_HORGRID_TO_MASSPOINTS, INTERP_VERGRID_TO_MASSPOINTS
+
+CONTAINS
+
 !     #####################################################################
       SUBROUTINE SET_GRID( KMI, TPINIFILE,                                &
                            KKU, KIMAX_ll, KJMAX_ll,                       &
@@ -405,14 +343,8 @@ CALL IO_Field_read(TPINIFILE,'DTSEG',TDTSEG)
 !
 
 ! Interpolations of positions to mass points
-PXHATM( : UBOUND(PXHATM,1)-1 ) = 0.5 * PXHAT( : UBOUND(PXHAT,1)-1 ) + 0.5 * PXHAT( LBOUND(PXHAT,1)+1 : UBOUND(PXHAT,1) )
-PXHATM( UBOUND(PXHATM,1)     ) = 1.5 * PXHAT( UBOUND(PXHAT,1)     ) - 0.5 * PXHAT( UBOUND(PXHAT,1)-1 )
-
-PYHATM( : UBOUND(PYHATM,1)-1 ) = 0.5 * PYHAT( : UBOUND(PYHAT,1)-1 ) + 0.5 * PYHAT( LBOUND(PYHAT,1)+1 : UBOUND(PYHAT,1) )
-PYHATM( UBOUND(PYHATM,1)     ) = 1.5 * PYHAT( UBOUND(PYHAT,1)     ) - 0.5 * PYHAT( UBOUND(PYHAT,1)-1 )
-
-PZHATM( : UBOUND(PZHATM,1)-1 ) = 0.5 * PZHAT( : UBOUND(PZHAT,1)-1 ) + 0.5 * PZHAT( LBOUND(PZHAT,1)+1 : UBOUND(PZHAT,1) )
-PZHATM( UBOUND(PZHATM,1)     ) = 1.5 * PZHAT( UBOUND(PZHAT,1)     ) - 0.5 * PZHAT( UBOUND(PZHAT,1)-1 )
+CALL INTERP_HORGRID_TO_MASSPOINTS( PXHAT, PYHAT, PXHATM, PYHATM )
+CALL INTERP_VERGRID_TO_MASSPOINTS( PZHAT, PZHATM )
 
 IF (LCARTESIAN) THEN
   CALL SM_GRIDCART(PXHAT,PYHAT,PZHAT,PZS,OSLEVE,PLEN1,PLEN2,PZSMT,PDXHAT,PDYHAT,PZZ,PJ)
@@ -496,3 +428,77 @@ CALL SM_PRINT_TIME(TDTSEG,TLUOUT,YTITLE)
 !-------------------------------------------------------------------------------
 !
 END SUBROUTINE SET_GRID
+
+
+SUBROUTINE INTERP_HORGRID_1DIR_TO_MASSPOINTS( HDIR, PHAT, PHATM )
+! Interpolate 1 direction of horizontal grid to mass points
+
+USE MODD_ARGSLIST_ll, ONLY: LIST1D_ll
+
+USE MODE_ARGSLIST_ll, ONLY: ADD1DFIELD_ll, CLEANLIST1D_ll
+USE MODE_EXCHANGE_ll, ONLY: UPDATE_1DHALO_ll
+USE MODE_MSG
+
+IMPLICIT NONE
+
+CHARACTER(LEN=1),   INTENT(IN)  :: HDIR  ! Direction ('X' or 'Y')
+REAL, DIMENSION(:), INTENT(IN)  :: PHAT  ! Position x or y in the conformal or cartesian plane
+REAL, DIMENSION(:), INTENT(OUT) :: PHATM ! Position x or y in the conformal or cartesian plane at mass points
+
+CHARACTER(LEN=:), ALLOCATABLE :: YNAME
+INTEGER                       :: IINFO_ll ! return code
+TYPE(LIST1D_ll),  POINTER     :: TZLIST   ! pointer for the list of 1D fields to be communicated
+
+
+! Interpolate inside subdomain
+PHATM( : UBOUND(PHATM,1)-1 ) = 0.5 * PHAT( : UBOUND(PHAT,1)-1 ) + 0.5 * PHAT( LBOUND(PHAT,1)+1 : UBOUND(PHAT,1) )
+PHATM( UBOUND(PHATM,1)     ) = 1.5 * PHAT( UBOUND(PHAT,1)     ) - 0.5 * PHAT( UBOUND(PHAT,1)-1 )
+
+! Update data between subdomains
+NULLIFY( TZLIST )
+
+IF ( HDIR == 'X' ) THEN
+  YNAME = 'XHATM'
+ELSE IF ( HDIR == 'Y' ) THEN
+  YNAME = 'YHATM'
+ELSE
+  CALL PRINT_MSG( NVERB_ERROR, 'GEN', 'INTERP_HORGRID_1DIR_TO_MASSPOINTS', 'invalid direction (valid: X or Y)' )
+END IF
+
+CALL ADD1DFIELD_ll( HDIR, TZLIST, PHATM, YNAME )
+CALL UPDATE_1DHALO_ll( TZLIST, IINFO_ll )
+CALL CLEANLIST1D_ll( TZLIST )
+
+END SUBROUTINE INTERP_HORGRID_1DIR_TO_MASSPOINTS
+
+
+SUBROUTINE INTERP_HORGRID_TO_MASSPOINTS( PXHAT, PYHAT, PXHATM, PYHATM )
+
+IMPLICIT NONE
+
+REAL, DIMENSION(:), INTENT(IN)  :: PXHAT  ! Position x in the conformal or cartesian plane
+REAL, DIMENSION(:), INTENT(IN)  :: PYHAT  ! Position y in the conformal or cartesian plane
+REAL, DIMENSION(:), INTENT(OUT) :: PXHATM ! Position x in the conformal or cartesian plane at mass points
+REAL, DIMENSION(:), INTENT(OUT) :: PYHATM ! Position y in the conformal or cartesian plane at mass points
+
+CALL INTERP_HORGRID_1DIR_TO_MASSPOINTS( 'X', PXHAT, PXHATM )
+CALL INTERP_HORGRID_1DIR_TO_MASSPOINTS( 'Y', PYHAT, PYHATM )
+
+END SUBROUTINE INTERP_HORGRID_TO_MASSPOINTS
+
+
+PURE SUBROUTINE INTERP_VERGRID_TO_MASSPOINTS( PZHAT, PZHATM )
+! Interpolate vertical grid to mass points
+
+IMPLICIT NONE
+
+REAL, DIMENSION(:), INTENT(IN)  :: PZHAT  ! Position z in the conformal or cartesian plane
+REAL, DIMENSION(:), INTENT(OUT) :: PZHATM ! Position z in the conformal or cartesian plane at mass points
+
+PZHATM( : UBOUND(PZHATM,1)-1 ) = 0.5 * PZHAT( : UBOUND(PZHAT,1)-1 ) + 0.5 * PZHAT( LBOUND(PZHAT,1)+1 : UBOUND(PZHAT,1) )
+PZHATM( UBOUND(PZHATM,1)     ) = 1.5 * PZHAT( UBOUND(PZHAT,1)     ) - 0.5 * PZHAT( UBOUND(PZHAT,1)-1 )
+
+END SUBROUTINE INTERP_VERGRID_TO_MASSPOINTS
+
+
+END MODULE MODE_SET_GRID
diff --git a/src/MNH/spawn_grid2.f90 b/src/MNH/spawn_grid2.f90
index 7131b7d5a..6909c693e 100644
--- a/src/MNH/spawn_grid2.f90
+++ b/src/MNH/spawn_grid2.f90
@@ -167,6 +167,7 @@ USE MODD_BIKHARDT_n
 USE MODD_VAR_ll
 use mode_bikhardt
 USE MODE_ll
+USE MODE_SET_GRID,   only: INTERP_HORGRID_TO_MASSPOINTS
 USE MODE_TIME
 USE MODE_GRIDPROJ
 !
@@ -394,11 +395,6 @@ PLEN2    = XLEN21
   DEALLOCATE(ZXHAT_2D_F)
   DEALLOCATE(ZXHAT_EXTENDED_C)
   DEALLOCATE(ZXHAT_2D_C)
-
-  ! Interpolations of positions to mass points
-  PXHATM(1:IIU_C-1) = 0.5 * PXHAT(1:IIU_C-1) + 0.5 * PXHAT(2:IIU_C)
-  PXHATM(  IIU_C)   = 1.5 * PXHAT(  IIU_C)   - 0.5 * PXHAT(IIU_C-1)
-
 !
 !     YHAT
 !
@@ -459,8 +455,7 @@ PLEN2    = XLEN21
   DEALLOCATE(ZYHAT_2D_C)
 
   ! Interpolations of positions to mass points
-  PYHATM(1:IJU_C-1) = 0.5 * PYHAT(1:IJU_C-1) + 0.5 * PYHAT(2:IJU_C)
-  PYHATM(  IJU_C)   = 1.5 * PYHAT(  IJU_C)   - 0.5 * PYHAT(IJU_C-1)
+  CALL INTERP_HORGRID_TO_MASSPOINTS( PXHAT, PYHAT, PXHATM, PYHATM )
 
 !!$=======
 !!$  IXSIZE1=SIZE(XXHAT1)
diff --git a/src/MNH/write_surf_mnh.f90 b/src/MNH/write_surf_mnh.f90
index ef001e661..ca3d3d124 100644
--- a/src/MNH/write_surf_mnh.f90
+++ b/src/MNH/write_surf_mnh.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1997-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1997-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.
@@ -297,6 +297,7 @@ USE MODD_PARAMETERS,    ONLY: XUNDEF, JPHEXT
 use mode_field,          only: Find_field_id_from_mnhname
 use MODE_IO_FIELD_WRITE, only: IO_Field_write
 USE MODE_MSG
+USE MODE_SET_GRID,       only: INTERP_HORGRID_1DIR_TO_MASSPOINTS
 USE MODE_TOOLS_ll
 USE MODE_WRITE_SURF_MNH_TOOLS
 
@@ -451,8 +452,7 @@ IF (      (CSTORAGE_TYPE=='PG' .OR. CSTORAGE_TYPE=='SU')  &
       XXHAT(:) = ZW1D(1+NHALO:IIU-NHALO)
 
       ! Interpolations of positions to mass points
-      XXHATM( : UBOUND(XXHATM,1)-1 ) = 0.5 * XXHAT( : UBOUND(XXHAT,1)-1 ) + 0.5 * XXHAT( LBOUND(XXHAT,1)+1 : UBOUND(XXHAT,1) )
-      XXHATM( UBOUND(XXHATM,1)     ) = 1.5 * XXHAT( UBOUND(XXHAT,1)     ) - 0.5 * XXHAT( UBOUND(XXHAT,1)-1 )
+      CALL INTERP_HORGRID_1DIR_TO_MASSPOINTS( 'X', XXHAT, XXHATM )
     END IF
   END IF
   DEALLOCATE(ZW1D)
@@ -487,8 +487,7 @@ ELSE IF ( (CSTORAGE_TYPE=='PG' .OR. CSTORAGE_TYPE=='SU')  &
       XYHAT(:) = ZW1D(1+NHALO:IJU-NHALO)
 
       ! Interpolations of positions to mass points
-      XYHATM( : UBOUND(XYHATM,1)-1 ) = 0.5 * XYHAT( : UBOUND(XYHAT,1)-1 ) + 0.5 * XYHAT( LBOUND(XYHAT,1)+1 : UBOUND(XYHAT,1) )
-      XYHATM( UBOUND(XYHATM,1)     ) = 1.5 * XYHAT( UBOUND(XYHAT,1)     ) - 0.5 * XYHAT( UBOUND(XYHAT,1)-1 )
+      CALL INTERP_HORGRID_1DIR_TO_MASSPOINTS( 'Y', XYHAT, XYHATM )
     END IF
   END IF
   DEALLOCATE(ZW1D)
-- 
GitLab