diff --git a/src/LIB/SURCOUCHE/src/mode_field.f90 b/src/LIB/SURCOUCHE/src/mode_field.f90
index b6f5a6ab733fd7c13cbeeed96a4b0e09cb984dd5..a0df6d83f884d7c1e18f0f206f64e06440115d28 100644
--- a/src/LIB/SURCOUCHE/src/mode_field.f90
+++ b/src/LIB/SURCOUCHE/src/mode_field.f90
@@ -440,7 +440,6 @@ call Add_field2list( TFIELDDATA( &
 
 call Add_field2list( TFIELDDATA( &
   CMNHNAME   = 'XHAT',           &
-!TODO: check stdname
   CSTDNAME   = 'projection_x_coordinate', &
   CLONGNAME  = 'XHAT',           &
   CUNITS     = 'm',              &
@@ -453,7 +452,6 @@ call Add_field2list( TFIELDDATA( &
 
 call Add_field2list( TFIELDDATA( &
   CMNHNAME   = 'YHAT',           &
-!TODO: check stdname
   CSTDNAME   = 'projection_y_coordinate', &
   CLONGNAME  = 'YHAT',           &
   CUNITS     = 'm',              &
@@ -464,6 +462,30 @@ call Add_field2list( TFIELDDATA( &
   NDIMS      = 1,                &
   LTIMEDEP   = .FALSE.           ) )
 
+call Add_field2list( TFIELDDATA( &
+  CMNHNAME   = 'XHATM',          &
+  CSTDNAME   = 'projection_x_coordinate', &
+  CLONGNAME  = 'XHATM',          &
+  CUNITS     = 'm',              &
+  CDIR       = 'XX',             &
+  CCOMMENT   = 'Position x in the conformal or cartesian plane at mass point', &
+  NGRID      = 1,                &
+  NTYPE      = TYPEREAL,         &
+  NDIMS      = 1,                &
+  LTIMEDEP   = .FALSE.           ) )
+
+call Add_field2list( TFIELDDATA( &
+  CMNHNAME   = 'YHATM',          &
+  CSTDNAME   = 'projection_y_coordinate', &
+  CLONGNAME  = 'YHATM',          &
+  CUNITS     = 'm',              &
+  CDIR       = 'YY',             &
+  CCOMMENT   = 'Position y in the conformal or cartesian plane at mass point', &
+  NGRID      = 1,                &
+  NTYPE      = TYPEREAL,         &
+  NDIMS      = 1,                &
+  LTIMEDEP   = .FALSE.           ) )
+
 call Add_field2list( TFIELDDATA( &
   CMNHNAME   = 'ZHAT',           &
 !TODO: check stdname
@@ -477,6 +499,18 @@ call Add_field2list( TFIELDDATA( &
   NDIMS      = 1,                &
   LTIMEDEP   = .FALSE.           ) )
 
+call Add_field2list( TFIELDDATA( &
+  CMNHNAME   = 'ZHATM',          &
+  CSTDNAME   = '',               &
+  CLONGNAME  = 'ZHATM',          &
+  CUNITS     = 'm',              &
+  CDIR       = 'ZZ',             &
+  CCOMMENT   = 'Height level without orography at mass point', &
+  NGRID      = 4,                &
+  NTYPE      = TYPEREAL,         &
+  NDIMS      = 1,                &
+  LTIMEDEP   = .FALSE.           ) )
+
 call Add_field2list( TFIELDDATA(      &
   CMNHNAME   = 'ZTOP',                &
   CSTDNAME   = 'altitude_at_top_of_atmosphere_model', &
@@ -3594,7 +3628,10 @@ call Goto_model_1field( 'ZS'  ,  kfrom, kto, XZS    )
 call Goto_model_1field( 'ZSMT',  kfrom, kto, XZSMT  )
 call Goto_model_1field( 'XHAT',  kfrom, kto, XXHAT  )
 call Goto_model_1field( 'YHAT',  kfrom, kto, XYHAT  )
+call Goto_model_1field( 'XHATM', kfrom, kto, XXHATM )
+call Goto_model_1field( 'YHATM', kfrom, kto, XYHATM )
 call Goto_model_1field( 'ZHAT',  kfrom, kto, XZHAT  )
+call Goto_model_1field( 'ZHATM', kfrom, kto, XZHATM )
 call Goto_model_1field( 'ZTOP',  kfrom, kto, XZTOP  )
 call Goto_model_1field( 'DXHAT', kfrom, kto, XDXHAT )
 call Goto_model_1field( 'DYHAT', kfrom, kto, XDYHAT )
diff --git a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90 b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
index 8ca884014dae686642a85c7892796d175d604a3d..65ffbdc8e20b5ad6b6494ecbc97899a40c76829f 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_write_nc4.f90
@@ -1450,7 +1450,7 @@ use modd_field,      only: NMNHDIM_NI, NMNHDIM_NJ, NMNHDIM_NI_U, NMNHDIM_NJ_U, N
                            NMNHDIM_PROFILER_TIME, NMNHDIM_STATION_TIME,                                    &
                            tfieldlist
 use modd_grid,       only: xlatori, xlonori
-use modd_grid_n,     only: lsleve, xxhat, xyhat, xzhat
+use modd_grid_n,     only: lsleve, xxhat, xxhatm, xyhat, xyhatm, xzhat, xzhatm
 use modd_les,        only: cles_level_type, cspectra_level_type, nlesn_iinf, nlesn_isup, nlesn_jinf, nlesn_jsup, &
                            nles_k, nles_levels, nspectra_k, nspectra_levels,                                     &
                            xles_altitudes, xspectra_altitudes
@@ -1475,7 +1475,7 @@ character(len=*), optional, intent(in) :: hprogram_orig !To emulate a file comin
 
 character(len=:),                         allocatable :: ystdnameprefix
 character(len=:),                         allocatable :: yprogram
-integer                                               :: iiu, iju, iku
+integer                                               :: iiu, iju
 integer                                               :: id, iid, iresp
 integer                                               :: imi
 integer                                               :: ji
@@ -1486,7 +1486,7 @@ logical                                               :: gchangemodel
 logical                                               :: gdealloc
 logical,                         pointer              :: gsleve
 real,            dimension(:),   pointer              :: zxhat, zyhat, zzhat
-real,            dimension(:),            allocatable :: zxhatm, zyhatm, zzhatm !Coordinates at mass points in the transformed space
+real,            dimension(:),   pointer              :: zxhatm, zyhatm, zzhatm !Coordinates at mass points in the transformed space
 real,            dimension(:),            allocatable :: zles_levels
 real,            dimension(:),            allocatable :: zspectra_levels
 real,            dimension(:,:), pointer              :: zlat, zlon
@@ -1505,9 +1505,12 @@ real, dimension(:,:), pointer, save :: zlatf_glob  => null(), zlonf_glob  => nul
 
 call Print_msg( NVERB_DEBUG, 'IO', 'IO_Coordvar_write_nc4', 'called for ' // Trim( tpfile%cname ) )
 
-zxhat => null()
-zyhat => null()
-zzhat => null()
+zxhat  => null()
+zyhat  => null()
+zzhat  => null()
+zxhatm => null()
+zyhatm => null()
+zzhatm => null()
 
 gchangemodel = .false.
 
@@ -1527,8 +1530,14 @@ if ( tpfile%nmodel > 0 ) then
   zxhat => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
   call Find_field_id_from_mnhname( 'YHAT', iid, iresp )
   zyhat => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  call Find_field_id_from_mnhname( 'XHATM', iid, iresp )
+  zxhatm => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  call Find_field_id_from_mnhname( 'YHATM', iid, iresp )
+  zyhatm => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
   call Find_field_id_from_mnhname( 'ZHAT', iid, iresp )
   zzhat => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
+  call Find_field_id_from_mnhname( 'ZHATM', iid, iresp )
+  zzhatm => tfieldlist(iid)%tfield_x1d(tpfile%nmodel)%data
   call Find_field_id_from_mnhname( 'SLEVE', iid, iresp )
   gsleve => tfieldlist(iid)%tfield_l0d(tpfile%nmodel)%data
 
@@ -1538,21 +1547,17 @@ if ( tpfile%nmodel > 0 ) then
     gchangemodel = .true.
   end if
 else
-  zxhat => xxhat
-  zyhat => xyhat
-  zzhat => xzhat
+  zxhat  => xxhat
+  zyhat  => xyhat
+  zzhat  => xzhat
+  zxhatm => xxhatm
+  zyhatm => xyhatm
+  zzhatm => xzhatm
   gsleve => lsleve
 end if
 
 iiu = Size( zxhat )
 iju = Size( zyhat )
-Allocate( zxhatm(iiu), zyhatm(iju) )
-!zxhatm(iiu) and zyhatm(iju) are correct only on some processes
-!but it is OK due to the way Gather_xxfield is done
-zxhatm(1 : iiu - 1) = 0.5 * ( zxhat(1 : iiu - 1) + zxhat(2 : iiu) )
-zxhatm(iiu)         = 2. * zxhat(iiu) - zxhatm(iiu - 1)
-zyhatm(1 : iju - 1) = 0.5 * ( zyhat(1 : iju - 1) + zyhat(2 : iju) )
-zyhatm(iju)         = 2. * zyhat(iju) - zyhatm(iju - 1)
 
 if ( lcartesian ) then
   ystdnameprefix = 'plane'
@@ -1634,17 +1639,10 @@ if ( .not. lcartesian ) then
   if ( gdealloc ) Deallocate( zlatm_glob, zlonm_glob, zlatu_glob, zlonu_glob, zlatv_glob, zlonv_glob, zlatf_glob, zlonf_glob )
 end if
 
-Deallocate( zxhatm, zyhatm )
-
 if ( tpfile%lmaster ) then !vertical coordinates in the transformed space are the same on all processes
   if ( Trim( yprogram ) /= 'PGD' .and. Trim( yprogram ) /= 'NESPGD' .and. Trim( yprogram ) /= 'ZOOMPG' &
       .and. .not. ( Trim( yprogram ) == 'REAL' .and. cstorage_type == 'SU') ) then !condition to detect prep_surfex
 
-    iku = Size( zzhat )
-    Allocate( zzhatm(iku) )
-    zzhatm(1 : iku - 1) = 0.5 * ( zzhat(2 : iku) + zzhat(1 : iku - 1) )
-    zzhatm(iku)         = 2.* zzhat(iku) - zzhatm(iku - 1)
-
     call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_LEVEL),  'position z in the transformed space',              '', &
                           'altitude',                0.,  JPVEXT, JPVEXT, ZZHATM )
     call Write_ver_coord( tpfile%tncdims%tdims(NMNHDIM_LEVEL_W),'position z in the transformed space at w location','', &
diff --git a/src/MNH/aircraft_balloon_evol.f90 b/src/MNH/aircraft_balloon_evol.f90
index c79631fa0e07aa38c1e42c2e2c6e3bf70ce27821..e6cfc21bad4548fe49bd153d979a04d2a7a1739d 100644
--- a/src/MNH/aircraft_balloon_evol.f90
+++ b/src/MNH/aircraft_balloon_evol.f90
@@ -138,6 +138,7 @@ USE MODD_CONF
 USE MODD_CST
 USE MODD_DIAG_IN_RUN
 USE MODD_GRID
+USE MODD_GRID_n,           ONLY: XXHATM, XYHATM
 USE MODD_LUNIT_n,          ONLY: TLUOUT
 USE MODD_NESTING
 USE MODD_NSV,              ONLY : NSV_LIMA_NI,NSV_LIMA_NR,NSV_LIMA_NC
@@ -229,9 +230,6 @@ INTEGER :: IKU
 !
 INTEGER :: JK         ! loop index
 !
-REAL, DIMENSION(SIZE(PXHAT))        :: ZXHATM ! mass point coordinates
-REAL, DIMENSION(SIZE(PYHAT))        :: ZYHATM ! mass point coordinates
-!
 REAL, DIMENSION(2,2,SIZE(PZ,3))     :: ZZM    ! mass point coordinates
 REAL, DIMENSION(2,2,SIZE(PZ,3))     :: ZZU    ! U points z coordinates
 REAL, DIMENSION(2,2,SIZE(PZ,3))     :: ZZV    ! V points z coordinates
@@ -384,11 +382,6 @@ IKU = SIZE(PZ,3)
 IIU=SIZE(PXHAT)
 IJU=SIZE(PYHAT)
 !
-ZXHATM(1:IIU-1)=0.5*PXHAT(1:IIU-1)+0.5*PXHAT(2:IIU  )
-ZXHATM(  IIU  )=1.5*PXHAT(  IIU  )-0.5*PXHAT(  IIU-1)
-!
-ZYHATM(1:IJU-1)=0.5*PYHAT(1:IJU-1)+0.5*PYHAT(2:IJU  )
-ZYHATM(  IJU  )=1.5*PYHAT(  IJU  )-0.5*PYHAT(  IJU-1)
 !----------------------------------------------------------------------------
 !
 !*      2.3  Compute time until launch by comparison of dates and times
@@ -504,7 +497,7 @@ IF ( TPFLYER%LFLY ) THEN
 !            ----------
 !
   IU=COUNT( PXHAT (:)<=TPFLYER%XX_CUR )
-  II=COUNT( ZXHATM(:)<=TPFLYER%XX_CUR )
+  II=COUNT( XXHATM(:)<=TPFLYER%XX_CUR )
 !
   IF ( IU < IIB .AND. LWEST_ll() ) THEN
     IF ( TPFLYER%CMODEL == 'FIX' .OR. TPFLYER%NMODEL == 1 ) THEN
@@ -528,7 +521,7 @@ IF ( TPFLYER%LFLY ) THEN
 !            ----------
 !
   IV=COUNT( PYHAT (:)<=TPFLYER%XY_CUR )
-  IJ=COUNT( ZYHATM(:)<=TPFLYER%XY_CUR )
+  IJ=COUNT( XYHATM(:)<=TPFLYER%XY_CUR )
 !
   IF ( IV < IJB .AND. LSOUTH_ll() ) THEN
     IF ( TPFLYER%CMODEL == 'FIX'  .OR. TPFLYER%NMODEL == 1 ) THEN
@@ -632,9 +625,9 @@ IF ( TPFLYER%LFLY ) THEN
 !*      5.2.1 Iso-density balloon
 !
             CASE ( 'ISODEN' )
-              ZXCOEF = (TPFLYER%XX_CUR - ZXHATM(II)) / (ZXHATM(II+1) - ZXHATM(II))
+              ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II)) / (XXHATM(II+1) - XXHATM(II))
               ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
-              ZYCOEF = (TPFLYER%XY_CUR - ZYHATM(IJ)) / (ZYHATM(IJ+1) - ZYHATM(IJ))
+              ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ)) / (XYHATM(IJ+1) - XYHATM(IJ))
               ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
               IF ( TPFLYER%XALTLAUNCH /= XNEGUNDEF ) THEN
                 IK00 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(1,1,:)), 1)
@@ -676,9 +669,9 @@ IF ( TPFLYER%LFLY ) THEN
 !*      5.2.4 Constant Volume Balloon
 !
             CASE ( 'CVBALL' )
-              ZXCOEF = (TPFLYER%XX_CUR - ZXHATM(II)) / (ZXHATM(II+1) - ZXHATM(II))
+              ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II)) / (XXHATM(II+1) - XXHATM(II))
               ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
-              ZYCOEF = (TPFLYER%XY_CUR - ZYHATM(IJ)) / (ZYHATM(IJ+1) - ZYHATM(IJ))
+              ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ)) / (XYHATM(IJ+1) - XYHATM(IJ))
               ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
               IF ( TPFLYER%XALTLAUNCH /= XNEGUNDEF ) THEN
                 IK00 = MAX ( COUNT (TPFLYER%XALTLAUNCH >= ZZM(1,1,:)), 1)
@@ -837,14 +830,14 @@ IF ( TPFLYER%LFLY ) THEN
 !*      6.1  Interpolation coefficient for X
 !            -------------------------------
 !
-      ZXCOEF = (TPFLYER%XX_CUR - ZXHATM(II)) / (ZXHATM(II+1) - ZXHATM(II))
+      ZXCOEF = (TPFLYER%XX_CUR - XXHATM(II)) / (XXHATM(II+1) - XXHATM(II))
       ZXCOEF = MAX (0.,MIN(ZXCOEF,1.))
 !
 !
 !*      6.2  Interpolation coefficient for y
 !            -------------------------------
 !
-      ZYCOEF = (TPFLYER%XY_CUR - ZYHATM(IJ)) / (ZYHATM(IJ+1) - ZYHATM(IJ))
+      ZYCOEF = (TPFLYER%XY_CUR - XYHATM(IJ)) / (XYHATM(IJ+1) - XYHATM(IJ))
       ZYCOEF = MAX (0.,MIN(ZYCOEF,1.))
 !
 !
diff --git a/src/MNH/eddy_fluxn.f90 b/src/MNH/eddy_fluxn.f90
index ed52def268d0ea872591a158078514d3fef3316e..dc20500dcd9390c6bd54a6847fbdf54ea661b953 100644
--- a/src/MNH/eddy_fluxn.f90
+++ b/src/MNH/eddy_fluxn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2004-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2004-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.
@@ -302,8 +302,8 @@ DO JJ=IJB,IJE
       !
         ZKYY(JI,JJ,JK) =  &
             (ZA0(JI) * XG*ZNDW(JI,JJ)/(ZTHM_W(JI,JJ)*(ZCORIOZ(JI,JJ,JK)**2)) ) &
-            * (ZDW(JI,JJ)**2) * (EXP(-0.5*(XZHAT(JK)+XZHAT(JK+1))/ZDW(JI,JJ))) &
-            * ZADTDXW(JI,JJ)   
+            * ZDW(JI,JJ)**2 * EXP( - XZHATM(JK) / ZDW(JI,JJ) ) &
+            * ZADTDXW(JI,JJ)
       ENDDO       
       !
     ! CASE WHERE NO INSTABILITY IS DETECTED
@@ -323,8 +323,7 @@ ENDDO
 
 !
 DO JK=IKB,IKE
-   ZVTH_FLUX(:,:,JK) = - 0.5 * ZKYY(:,:,JK)*ZDTHM_DY(:,:,JK) * &
-                    (1-EXP(-0.5*(XZHAT(JK)+XZHAT(JK+1))/ZDELTAZ))
+   ZVTH_FLUX(:,:,JK) = - 0.5 * ZKYY(:,:,JK) * ZDTHM_DY(:,:,JK) * ( 1 - EXP( -XZHATM(JK) / ZDELTAZ ) )
 END DO
 !
 !       2.1 Smoothing in equatorial region
diff --git a/src/MNH/eol_maths.f90 b/src/MNH/eol_maths.f90
index bda3aa6949d919d1616eb68b28bef118ab7802fd..81e2a0a885d7012caff260547722e3829d0e7aee 100644
--- a/src/MNH/eol_maths.f90
+++ b/src/MNH/eol_maths.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2018-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2018-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.
@@ -283,7 +283,7 @@ END FUNCTION INTERP_SPLCUB
 !#########################################################
 FUNCTION INTERP_LIN8NB(PPOS, KI, KJ, KK, PVAR, PZH)
 !
-USE MODD_GRID_n, ONLY: XXHAT,XYHAT
+USE MODD_GRID_n, ONLY: XDXHAT, XXHATM, XDYHAT, XYHATM
 !
 REAL                               :: INTERP_LIN8NB  ! Return
 REAL, DIMENSION(3),     INTENT(IN) :: PPOS           ! Position where we want to evaluate
@@ -309,7 +309,7 @@ REAL     :: ZUX                           ! Interpolated variable (VAR) in Z pla
 !
 ! FINDING 8 NEIGHBOORS 
 ! -- X axis
-IF (PPOS(1) <= 0.5*(XXHAT(KI)+XXHAT(KI+1))) THEN
+IF (PPOS(1) <= XXHATM(KI)) THEN
  IIP = KI - 1
  IIN = KI
 ELSE   
@@ -317,7 +317,7 @@ ELSE
  IIN = KI + 1
 END IF
 ! -- Y axis
-IF (PPOS(2) <= 0.5*((XYHAT(KJ)+XYHAT(KJ+1)))) THEN
+IF (PPOS(2) <= XYHATM(KJ)) THEN
  IJP = KJ - 1
  IJN = KJ
 ELSE   
@@ -336,7 +336,7 @@ END IF
 ! INTERPOLATION 
 ! -- Along X
 ! -- -- Alpha
-ZALPHAX = (PPOS(1) -  0.5*(XXHAT(IIP)+XXHAT(IIN))) / (XXHAT(IIN) - XXHAT(IIP))
+ZALPHAX = (PPOS(1) -  XXHATM(IIP)) / XDXHAT(IIP)
 !!PRINT*, "ZALPHAX = ", ZALPHAX
 ! -- -- -- Wind
 ! -- -- Interpolated variable in temporary plane X
@@ -353,7 +353,7 @@ ZHXPN = (1-ZALPHAX)*PZH(IIP,IJP,IKN) + ZALPHAX*PZH(IIN,IJP,IKN)
 !
 ! -- Along Y
 ! -- -- Alpha
-ZALPHAY = (PPOS(2) -  0.5*(XYHAT(IJP)+XYHAT(IJN))) / (XYHAT(IJN) - XYHAT(IJP))
+ZALPHAY = (PPOS(2) -  XYHATM(IJP)) / XDYHAT(IJP)
 !PRINT*, "ZALPHAY = ", ZALPHAY
 ! -- -- Interpolated variable in temporary plane Y
 ! -- -- -- Wind
diff --git a/src/MNH/extract_vortex.f90 b/src/MNH/extract_vortex.f90
index df1bb103b63398cd752ce8f6d5b2dbc5adeda879..311fcecd19f99b461609155dfd3ae1df6ee84159 100644
--- a/src/MNH/extract_vortex.f90
+++ b/src/MNH/extract_vortex.f90
@@ -1,12 +1,8 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2001-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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-!-----------------------------------------------------------------
 !     ##########################
       MODULE MODI_EXTRACT_VORTEX
 !     ##########################
@@ -102,7 +98,7 @@ END MODULE MODI_EXTRACT_VORTEX
 USE MODD_CST, ONLY: XPI
 USE MODD_PARAMETERS, ONLY: XUNDEF
 USE MODD_DIM_n,      ONLY: NIMAX,NJMAX
-USE MODD_GRID_n,     ONLY: XXHAT,XYHAT
+USE MODD_GRID_n,     ONLY: XDXHAT, XXHAT, XDYHAT, XYHAT
 USE MODE_ll
 !
 IMPLICIT NONE
@@ -142,8 +138,8 @@ IPHI=SIZE(PR0,1)
 !
 CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
 !
-ZDELTAX = XXHAT(3) - XXHAT(2)
-ZDELTAY = XYHAT(3) - XYHAT(2)
+ZDELTAX = XDXHAT(2)
+ZDELTAY = XDYHAT(2)
 ZDELTAR = MAX(ZDELTAX,ZDELTAY)
 ZDPHI = 2.*XPI / IPHI
 !
diff --git a/src/MNH/flash_geom_elec.f90 b/src/MNH/flash_geom_elec.f90
index 8ae281d266727272811789065a02bce9cad1ec3c..14265cba17427a97b008a4ae1fef61a785116e5b 100644
--- a/src/MNH/flash_geom_elec.f90
+++ b/src/MNH/flash_geom_elec.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2010-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2010-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.
@@ -101,6 +101,7 @@ END MODULE MODI_FLASH_GEOM_ELEC_n
 !  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
 !  P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
 !  P. Wautelet 18/09/2019: correct support of 64bit integers (MNH_INT=8)
+!  P. Wautelet 31/08/2022: remove ZXMASS and ZYMASS (use XXHATM and XYHATM instead)
 !-------------------------------------------------------------------------------
 !
 !*      0.      DECLARATIONS
@@ -119,7 +120,7 @@ USE MODD_ELEC_PARAM,     ONLY: XFQLIGHTR, XEXQLIGHTR, &
                                XFQLIGHTH, XEXQLIGHTH, &
                                XFQLIGHTC
 USE MODD_GRID,           ONLY: XLATORI,XLONORI
-USE MODD_GRID_n,         ONLY: XXHAT, XYHAT, XZHAT
+USE MODD_GRID_n,         ONLY: XXHATM, XYHATM, XZHAT
 USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_LMA_SIMULATOR
 USE MODD_METRICS_n,      ONLY: XDXX, XDYY, XDZZ ! in linox_production
@@ -352,8 +353,6 @@ INBFTS_MAX = ANINT(1000 * PTSTEP / 60)
 !
 IF (GEFIRSTCALL) THEN
   GEFIRSTCALL = .FALSE.
-  ALLOCATE (ZXMASS(SIZE(XXHAT)))
-  ALLOCATE (ZYMASS(SIZE(XYHAT)))
   ALLOCATE (ZZMASS(SIZE(PZZ,1), SIZE(PZZ,2), SIZE(PZZ,3)))
   ALLOCATE (ZPRES_COEF(SIZE(PZZ,1), SIZE(PZZ,2), SIZE(PZZ,3)))
   IF(LLMA) THEN
@@ -379,8 +378,6 @@ IF (GEFIRSTCALL) THEN
   ALLOCATE (ZSNEUT_POS(NFLASH_WRITE))
   ALLOCATE (ZSNEUT_NEG(NFLASH_WRITE))
 !
-  ZXMASS(IIB:IIE) = 0.5 * (XXHAT(IIB:IIE) + XXHAT(IIB+1:IIE+1))
-  ZYMASS(IJB:IJE) = 0.5 * (XYHAT(IJB:IJE) + XYHAT(IJB+1:IJE+1))
   ZZMASS = MZF(PZZ)
   ZPRES_COEF = EXP(ZZMASS/8400.)
   ZSCOORD_SEG(:,:,:) = 0.0
@@ -588,9 +585,8 @@ DEALLOCATE (ZEMAX)
 IF (INB_CELL .GE. 1) THEN
 !
 ! mean mesh size
-  ZMEAN_GRID = (XDXHATM**2 + XDYHATM**2 +                            &
-               (SUM(XZHAT(2:SIZE(PRT,3)) - XZHAT(1:SIZE(PRT,3)-1)) / &
-               (SIZE(PRT,3)-1.))**2)**0.5
+  ZMEAN_GRID = (XDXHATM**2 + XDYHATM**2 +                                       &
+                ( ( XZHAT(UBOUND(XZHAT,1)) - XZHAT(1) ) / (SIZE(PRT,3)-1.) )**2 )**0.5
 ! chaque proc calcule son propre zmean_grid
 ! mais cette valeur peut etre differente sur chaque proc (ex: relief)
 ! laisse tel quel pour le moment
@@ -913,8 +909,8 @@ ENDIF
             DO IJ = IJB, IJE
               DO IK = IKB, IKE
                 IF (GPROP(II,IJ,IK,IL)) THEN
-                  ZDIST(II,IJ,IK) = ((ZXMASS(II) - ZCOORD_TRIG(1,IL))**2 + &
-                                     (ZYMASS(IJ) - ZCOORD_TRIG(2,IL))**2 + &
+                  ZDIST(II,IJ,IK) = ((XXHATM(II) - ZCOORD_TRIG(1,IL))**2 + &
+                                     (XYHATM(IJ) - ZCOORD_TRIG(2,IL))**2 + &
                                      (ZZMASS(II,IJ,IK) - ZCOORD_TRIG(3,IL))**2)**0.5
                 END IF
               END DO
@@ -1664,8 +1660,8 @@ DO IL = 1, INB_CELL
           ISEG_GLOB(2,IL) = IJ_TRIG_GLOB
           ISEG_GLOB(3,IL) = IK_TRIG
 !
-          ZCOORD_TRIG(1,IL) = ZXMASS(II_TRIG_LOC)
-          ZCOORD_TRIG(2,IL) = ZYMASS(IJ_TRIG_LOC)
+          ZCOORD_TRIG(1,IL) = XXHATM(II_TRIG_LOC)
+          ZCOORD_TRIG(2,IL) = XYHATM(IJ_TRIG_LOC)
           ZCOORD_TRIG(3,IL) = ZZMASS(II_TRIG_LOC, IJ_TRIG_LOC, IK_TRIG)
 !
           ZCOORD_SEG(1:3,IL) = ZCOORD_TRIG(1:3,IL)
@@ -1767,8 +1763,8 @@ IF (IPROC .EQ. IPROC_TRIG(IL)) THEN
     ISEG_GLOB(IIDECAL+2,IL) = IJBL_LOC + IYOR - 1
     ISEG_GLOB(IIDECAL+3,IL) = IKBL
 !
-    ZCOORD_SEG(IIDECAL+1,IL) = ZXMASS(IIBL_LOC)
-    ZCOORD_SEG(IIDECAL+2,IL) = ZYMASS(IJBL_LOC)
+    ZCOORD_SEG(IIDECAL+1,IL) = XXHATM(IIBL_LOC)
+    ZCOORD_SEG(IIDECAL+2,IL) = XYHATM(IJBL_LOC)
     ZCOORD_SEG(IIDECAL+3,IL) = ZZMASS(IIBL_LOC, IJBL_LOC, IKBL)
 !
     INBSEG(IL) = INBSEG(IL) + 1
@@ -2239,8 +2235,8 @@ IF (INB_SEG_AFT .GT. INB_SEG_BEF) THEN
           ISEG_GLOB(IIDECALB+2,IL) = IJ + IYOR - 1
           ISEG_GLOB(IIDECALB+3,IL) = IK
 !
-          ZCOORD_SEG(IIDECALB+1,IL) = ZXMASS(II)
-          ZCOORD_SEG(IIDECALB+2,IL) = ZYMASS(IJ)
+          ZCOORD_SEG(IIDECALB+1,IL) = XXHATM(II)
+          ZCOORD_SEG(IIDECALB+2,IL) = XYHATM(IJ)
           ZCOORD_SEG(IIDECALB+3,IL) = ZZMASS(II,IJ,IK)
           INBSEG(IL) = INBSEG(IL) + 1
         END IF
diff --git a/src/MNH/fun.f90 b/src/MNH/fun.f90
index 46db8b201c4b7d5352b188bf5489b2e580bdbeb1..64c54caa6b08b2254b3819cf2c932ab21dd224af 100644
--- a/src/MNH/fun.f90
+++ b/src/MNH/fun.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2017 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -126,9 +126,9 @@ ZWIDTHY =ZYHAT_ll(IJ0+IJU_ll/5)-ZYHAT_ll(IJ0)
 ZWIDTHZ =PZHAT(IK0+KKU/5)-PZHAT(IK0)
 DO JJ = 1,KJU-1
   DO JK = 1,KKU
-    FUNUYZ(JJ,JK) = 1./COSH(                                    &
-     (( (XYHAT(JJ)+XYHAT(JJ+1))*0.5-ZYHAT_ll(IJ0))/ZWIDTHY) **2 &
-    +((  PZHAT(JK)                 -   PZHAT(IK0))/ZWIDTHZ) **2 )
+    FUNUYZ(JJ,JK) = 1./COSH(                   &
+     (( XYHATM(JJ)-ZYHAT_ll(IJ0))/ZWIDTHY) **2 &
+    +(( PZHAT(JK) -   PZHAT(IK0))/ZWIDTHZ) **2 )
   END DO
 END DO
 DEALLOCATE(ZYHAT_ll)
@@ -217,7 +217,7 @@ ALLOCATE(ZYHAT_ll(IJU_ll))
 CALL GATHERALL_FIELD_ll('YY',XYHAT,ZYHAT_ll,IINFO_ll)
 ZWIDTH=ZYHAT_ll(IJ0+IJU_ll/5)-ZYHAT_ll(IJ0)
 DO JJ = 1,KJU-1
- FUNUY(JJ) = 1./COSH(((XYHAT(JJ)+XYHAT(JJ+1))*0.5-ZYHAT_ll(IJ0))/ZWIDTH)
+ FUNUY(JJ) = 1./COSH((XYHATM(JJ)-ZYHAT_ll(IJ0))/ZWIDTH)
 END DO
 FUNUY(KJU)=2.*FUNUY(KJU-1)-FUNUY(KJU-2) !simple extrapolation for the last point
 !
@@ -314,9 +314,9 @@ ZWIDTHX=ZXHAT_ll(II0+IIU_ll/5)-ZXHAT_ll(II0)
 ZWIDTHZ=PZHAT(IK0+KKU/5)-PZHAT(IK0)
 DO JI = 1,KIU-1
   DO JK = 1,KKU
-    FUNVXZ(JI,JK) = 1./COSH(                                   &
-     (( (XXHAT(JI)+XXHAT(JI+1))*0.5-ZXHAT_ll(II0))/ZWIDTHX)**2 &
-    +((  PZHAT (JK)                 - PZHAT (IK0))/ZWIDTHZ)**2 )
+    FUNVXZ(JI,JK) = 1./COSH(                  &
+     (( XXHATM(JI)-ZXHAT_ll(II0))/ZWIDTHX)**2 &
+    +(( PZHAT (JK) - PZHAT (IK0))/ZWIDTHZ)**2 )
   END DO
 END DO
 FUNVXZ(KIU,:)=2.*FUNVXZ(KIU-1,:)-FUNVXZ(KIU-2,:) !simple extrapolation for the last point
@@ -403,7 +403,7 @@ ALLOCATE(ZXHAT_ll(IIU_ll))
 CALL GATHERALL_FIELD_ll('XX',XXHAT,ZXHAT_ll,IINFO_ll)
 ZWIDTH=ZXHAT_ll(II0+IIU_ll/5)-ZXHAT_ll(II0)
 DO JI = 1,KIU
- FUNVX(JI)=1./COSH(((XXHAT(JI)+XXHAT(JI))*0.5-ZXHAT_ll(II0))/ZWIDTH)
+ FUNVX(JI)=1./COSH((XXHATM(JI)-ZXHAT_ll(II0))/ZWIDTH)
 END DO
 DEALLOCATE(ZXHAT_ll)
 !-------------------------------------------------------------------------------
diff --git a/src/MNH/gps_zenith.f90 b/src/MNH/gps_zenith.f90
index a11a074d6f511b7909c55934712aa9294cd0a3c2..5fe9cd741b0508c9b21133516957926682ae0942 100644
--- a/src/MNH/gps_zenith.f90
+++ b/src/MNH/gps_zenith.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2004-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2004-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.
@@ -66,8 +66,8 @@ END MODULE MODI_GPS_ZENITH
 !!    -------------
 !!      Original    18/11/04
 !!      Modified    4/12/2007
-!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!!
+!  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
+!  P. Wautelet 31/08/2022: use XXHATM and XYHATM (remove ZXHATM and ZYHATM local variables)
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -120,7 +120,6 @@ INTEGER :: JI,JJ,JK               ! Loop variables of control
 INTEGER :: IIU,IJU,IKU      ! Loop variables of model 
 INTEGER :: ILUOUT0, IRESP   ! file unit and return code for output
 INTEGER :: JL               !
-REAL,  DIMENSION(:),ALLOCATABLE         :: ZXHATM,ZYHATM   ! mass-point positions 
 REAL,  DIMENSION(:,:,:),ALLOCATABLE     :: ZZHATM   ! mass level altitude  
 !-------- Physical parameters for the integration ----------------------------
 REAL, DIMENSION(:,:,:),ALLOCATABLE ::  ZE              !  Partial pressure of water vapor
@@ -166,8 +165,6 @@ IKU = SIZE (PTEMP,3)
 IKB = JPVEXT + 1
 IKE = IKU - JPVEXT
 !
-ALLOCATE(ZXHATM(IIU))
-ALLOCATE(ZYHATM(IJU))
 ALLOCATE(ZZHATM(IIU,IJU,IKU))
 ALLOCATE(ZE(IIU,IJU,IKU))
 ALLOCATE(ZTV(IIU,IJU,IKU))
@@ -203,10 +200,6 @@ ZRDSRV=XRD/XRV
 !              ------------------- 
 !
 ! 
-ZXHATM(1:IIU-1) = 0.5*(XXHAT(1:IIU-1)+XXHAT(2:IIU))
-ZXHATM(IIU)     = 2.*XXHAT(IIU)-ZXHATM(IIU-1)
-ZYHATM(1:IJU-1) = 0.5*(XYHAT(1:IJU-1)+XYHAT(2:IJU))
-ZYHATM(IJU)     = 2.*XYHAT(IJU)-ZYHATM(IJU-1)  
 ZZHATM(:,:,1:IKU-1)=0.5*(XZZ(:,:,1:IKU-1)+XZZ(:,:,2:IKU))
 ZZHATM(:,:,IKU)    = 2.*XZZ(:,:,IKU) -ZZHATM(:,:,IKU-1)
 ! 
@@ -304,7 +297,7 @@ IF (ISTATIONS >0 ) THEN
     CALL SM_XYHAT(XLATORI,XLONORI,   &
        XLAT_GPS(JL),XLON_GPS(JL),ZXHAT_GPS(JL),ZYHAT_GPS(JL))
 !
-    II(JL)=COUNT(ZXHATM(:)<=ZXHAT_GPS(JL))
+    II(JL)=COUNT(XXHATM(:)<=ZXHAT_GPS(JL))
     IX=COUNT(XXHAT(:)<=ZXHAT_GPS(JL))
     IF (IX<IIB .AND. LWEST_ll()) THEN
      ! station outside the MESO-NH domain 
@@ -314,7 +307,7 @@ IF (ISTATIONS >0 ) THEN
      ! station outside the MESO-NH domain 
       GSTATION(JL)=.FALSE.
     ENDIF
-    IJ(JL)=COUNT(ZYHATM(:)<=ZYHAT_GPS(JL))
+    IJ(JL)=COUNT(XYHATM(:)<=ZYHAT_GPS(JL))
     IY=COUNT(XYHAT(:)<=ZYHAT_GPS(JL))
     IF (IY<IJB .AND. LSOUTH_ll()) THEN
      ! stations outside MESO-NH domain 
@@ -337,7 +330,7 @@ IF (ISTATIONS >0 ) THEN
       II1=II(JL)
       IJ1=IJ(JL)
 !     interpolate Z at station position and check that the difference between model relief and station altitude is weaker than XDIFORO
-      CALL INTERPOL_STATION(XZZ(:,:,:),ZXHATM,ZYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZZ_STA(:,JL))
+      CALL INTERPOL_STATION(XZZ(:,:,:),XXHATM,XYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZZ_STA(:,JL))
       IF ( ABS( ZZ_STA(IKB,JL)-XZS_GPS(JL)) > XDIFFORO ) THEN
         WRITE(IFGRI,*) 'NO DATA, Difference between the model orography and the GPS station height too large for ',CNAM_GPS(JL)
         GSTATION(JL)=.FALSE.
@@ -347,10 +340,10 @@ IF (ISTATIONS >0 ) THEN
 !        6.3 Interpolate to the station positions 
 !
 ! interpolate model variables to obs point
-      CALL INTERPOL_STATION(PPABSM(:,:,:),ZXHATM,ZYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZP_STA(:))
-      CALL INTERPOL_STATION(ZE(:,:,:),ZXHATM,ZYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZE_STA(:))
-      CALL INTERPOL_STATION(PTEMP(:,:,:),ZXHATM,ZYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZT_STA(:))
-      CALL INTERPOL_STATION(ZTV(:,:,:),ZXHATM,ZYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZTV_STA(:))
+      CALL INTERPOL_STATION(PPABSM(:,:,:),XXHATM,XYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZP_STA(:))
+      CALL INTERPOL_STATION(ZE(:,:,:),XXHATM,XYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZE_STA(:))
+      CALL INTERPOL_STATION(PTEMP(:,:,:),XXHATM,XYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZT_STA(:))
+      CALL INTERPOL_STATION(ZTV(:,:,:),XXHATM,XYHATM,II1,IJ1,ZXHAT_GPS(JL),ZYHAT_GPS(JL),ZTV_STA(:))
 !  
 !            6.3.1 For stations above model orography
 !  
@@ -386,7 +379,7 @@ IF (ISTATIONS >0 ) THEN
              /(XRD* 0.5 *(ZTV_STAT+ZTV_STA(IKB))))
         
 ! assume constant rvn for Vapor pressure
-        CALL INTERPOL_STATION( PRM(:,:,IKB),ZXHATM,ZYHATM,II(JL),IJ(JL),ZXHAT_GPS(JL), &
+        CALL INTERPOL_STATION( PRM(:,:,IKB),XXHATM,XYHATM,II(JL),IJ(JL),ZXHAT_GPS(JL), &
              ZYHAT_GPS(JL),ZRV_STAT)
         ZEM_STAT =  ZPM_STAT *  ZRV_STAT  / ( ZRDSRV + ZRV_STAT )
 ! add contribution below the model orography        
@@ -408,11 +401,11 @@ IF (ISTATIONS >0 ) THEN
 !
 !            6.3.4 Add external  contribution for ZHD
 !
-      CALL INTERPOL_STATION( ZPTOP(:,:),ZXHATM,ZYHATM,II(JL),IJ(JL),ZXHAT_GPS(JL), &
+      CALL INTERPOL_STATION( ZPTOP(:,:),XXHATM,XYHATM,II(JL),IJ(JL),ZXHAT_GPS(JL), &
            ZYHAT_GPS(JL),ZPTOP_STAT)
-      CALL INTERPOL_STATION( ZGTOP(:,:),ZXHATM,ZYHATM,II(JL),IJ(JL),ZXHAT_GPS(JL), &
+      CALL INTERPOL_STATION( ZGTOP(:,:),XXHATM,XYHATM,II(JL),IJ(JL),ZXHAT_GPS(JL), &
            ZYHAT_GPS(JL),ZGTOP_STAT)
-      CALL INTERPOL_STATION( ZTEMPTOP(:,:),ZXHATM,ZYHATM,II(JL),IJ(JL),ZXHAT_GPS(JL), &
+      CALL INTERPOL_STATION( ZTEMPTOP(:,:),XXHATM,XYHATM,II(JL),IJ(JL),ZXHAT_GPS(JL), &
            ZYHAT_GPS(JL),ZTEMPTOP_STAT)
       ZGPS_ZHD(JL)=ZGPS_ZHD(JL)+ ( 1.E-6 * ZK1 * ZPTOP_STAT * XRD * ( 1. + 2. * XRD * ZTEMPTOP_STAT &
            / ( ( XRADIUS + ZZ_STA(IKE+1,JL) ) * ZGTOP_STAT )  +  2. * ( XRD * ZTEMPTOP_STAT &
@@ -435,8 +428,6 @@ IF (ISTATIONS >0 ) THEN
   CALL IO_File_close(TZFILE,IRESP)
   PRINT *,'File ',TRIM(HFGRI),' closed, IRESP= ',IRESP
 !
-  DEALLOCATE(ZXHATM)
-  DEALLOCATE(ZYHATM)
   DEALLOCATE(ZGPS_ZTD)
   DEALLOCATE(ZGPS_ZHD)
   DEALLOCATE(ZGPS_ZWD)
diff --git a/src/MNH/gps_zenith_grid.f90 b/src/MNH/gps_zenith_grid.f90
index 86b72414fd58c1ab1ec065bc60a2e8c1f585e3b0..f3859e62cffe105ca97de8c6c8418eaa08348d98 100644
--- a/src/MNH/gps_zenith_grid.f90
+++ b/src/MNH/gps_zenith_grid.f90
@@ -1,6 +1,6 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2004-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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 !##########################################
@@ -66,8 +66,8 @@ END MODULE MODI_GPS_ZENITH_GRID
 !!      Modified    4/12/2007
 !!              July, 2015 (O.Nuissier/F.Duffourg) Add microphysics diagnostic for
 !!                                      aircraft, ballon and profiler
-!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!!
+!  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
+!  P. Wautelet 31/08/2022: remove ZXHATM and ZYHATM (unused variables)
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -107,7 +107,6 @@ INTEGER :: IJB,IJE          ! Loop limits for coordinate Y
 INTEGER :: IKB,IKE          ! Loop limits for coordinate Z
 INTEGER :: JK               ! Loop variables of control
 INTEGER :: IIU,IJU,IKU      ! Loop variables of model 
-REAL,  DIMENSION(:),ALLOCATABLE         :: ZXHATM,ZYHATM   ! mass-point positions 
 REAL,  DIMENSION(:,:,:),ALLOCATABLE     :: ZZHATM   ! mass level altitude  
 !-------- Physical parameters for the integration ----------------------------
 REAL, DIMENSION(:,:,:),ALLOCATABLE ::  ZE              !  Partial pressure of water vapor
@@ -135,8 +134,6 @@ IKU = SIZE (PTEMP,3)
 IKB = JPVEXT + 1
 IKE = IKU - JPVEXT
 !
-ALLOCATE(ZXHATM(IIU))
-ALLOCATE(ZYHATM(IJU))
 ALLOCATE(ZZHATM(IIU,IJU,IKU))
 ALLOCATE(ZE(IIU,IJU,IKU))
 ALLOCATE(ZTV(IIU,IJU,IKU))
@@ -172,10 +169,6 @@ ZRDSRV=XRD/XRV
 !              ------------------- 
 !
 ! 
-ZXHATM(1:IIU-1) = 0.5*(XXHAT(1:IIU-1)+XXHAT(2:IIU))
-ZXHATM(IIU)     = 2.*XXHAT(IIU)-ZXHATM(IIU-1)
-ZYHATM(1:IJU-1) = 0.5*(XYHAT(1:IJU-1)+XYHAT(2:IJU))
-ZYHATM(IJU)     = 2.*XYHAT(IJU)-ZYHATM(IJU-1)  
 ZZHATM(:,:,1:IKU-1)=0.5*(XZZ(:,:,1:IKU-1)+XZZ(:,:,2:IKU))
 ZZHATM(:,:,IKU)    = 2.*XZZ(:,:,IKU) -ZZHATM(:,:,IKU-1)
 ! 
diff --git a/src/MNH/ibm_affectp.f90 b/src/MNH/ibm_affectp.f90
index b0c998744f7971f296e01b421b8f959ed890136d..4f6bd44514015955aa2c20885d9f2535b20a898e 100644
--- a/src/MNH/ibm_affectp.f90
+++ b/src/MNH/ibm_affectp.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2019-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.
@@ -121,7 +121,7 @@ SUBROUTINE IBM_AFFECTP(PVAR,KIBM_LAYER,PRADIUS,PPOWERS,&
   USE MODD_RADIATIONS_n
   USE MODD_DYN_n
   USE MODD_FIELD_n
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT
+  USE MODD_GRID_n, ONLY: XDXHAT, XDYHAT
   !
   IMPLICIT NONE
   !
@@ -225,7 +225,7 @@ SUBROUTINE IBM_AFFECTP(PVAR,KIBM_LAYER,PRADIUS,PPOWERS,&
         DO JN = 1,3
            !
            Z_LOCAT_IMAG(JN,:)= XIBM_IMAGE_P(JM,JMM,1  ,JN,:)
-           Z_DELTA_IMAG      = ((XXHAT(JI+1)-XXHAT(JI))*(XYHAT(JJ+1)-XYHAT(JJ)))**0.5       
+           Z_DELTA_IMAG      = ( XDXHAT(JI) * XDYHAT(JJ) ) ** 0.5
            I_INDEX_CORN(:)   = NIBM_IMAGE_P(JM,JMM,1,1,JN,:)
            IF (I_INDEX_CORN(1)==0.AND.JN==2) ZIBM_HALO=0.
            IF (I_INDEX_CORN(2)==0.AND.JN==2) ZIBM_HALO=0.
diff --git a/src/MNH/ibm_affectv.f90 b/src/MNH/ibm_affectv.f90
index 1a5711e10bf7388cb89e7c965f4b657e51f52a1e..dddd95e97ca1c9ad4992d08e3ba86cdd58a98e8f 100644
--- a/src/MNH/ibm_affectv.f90
+++ b/src/MNH/ibm_affectv.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2019-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.
@@ -129,7 +129,7 @@ SUBROUTINE IBM_AFFECTV(PVAR,PVAR2,PVAR3,HVAR,KIBM_LAYER,HIBM_MODE_INTE3,&
   USE MODD_IBM_PARAM_n
   USE MODD_FIELD_n
   USE MODD_PARAM_n, ONLY: CTURB
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ
+  USE MODD_GRID_n,  ONLY: XDXHAT, XDYHAT
   USE MODD_VAR_ll,  ONLY: IP
   USE MODD_LBC_n
   USE MODD_REF_n, ONLY: XRHODJ,XRHODREF
@@ -291,7 +291,7 @@ SUBROUTINE IBM_AFFECTV(PVAR,PVAR2,PVAR3,HVAR,KIBM_LAYER,HIBM_MODE_INTE3,&
         DO JN = 1,3
            !
            Z_LOCAT_IMAG(JN,:)= XIBM_IMAGE_V(JM,JMM,JH  ,JN,:)
-           Z_DELTA_IMAG      = ((XXHAT(JI+1)-XXHAT(JI))*(XYHAT(JJ+1)-XYHAT(JJ)))**0.5
+           Z_DELTA_IMAG      = ( XDXHAT(JI) * XDYHAT(JJ) ) ** 0.5
            !
            DO JLL=1,3
               I_INDEX_CORN(:)       = NIBM_IMAGE_V(JM,JMM,JH,JLL,JN,:)
diff --git a/src/MNH/ibm_balance.f90 b/src/MNH/ibm_balance.f90
index 2256cd097bc547fd789da11f5fc85507a808242d..e1ed43c51bd559873d0e4adb640a49a9c68c9d66 100644
--- a/src/MNH/ibm_balance.f90
+++ b/src/MNH/ibm_balance.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2019-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.
@@ -57,6 +57,7 @@ SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
   !    MODIFICATIONS
   !    -------------
   !      Original         01/01/2019
+  !  P. Wautelet 31/08/2022: use XDXHAT and XDYHAT instead of XXHAT and XYHAT
   !
   !------------------------------------------------------------------------------
   !       
@@ -69,7 +70,7 @@ SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
   ! declaration
   USE MODD_CST, ONLY: XPI
   USE MODD_IBM_PARAM_n
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZHAT,XZZ 
+  USE MODD_GRID_n,     ONLY: XDXHAT, XDYHAT, XZZ
   USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT
   USE MODD_LBC_n
   USE MODD_REF_n
@@ -140,7 +141,7 @@ SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
               JL = 2
               JI2 = JI
               ZIBM_FLUX(JI2,JJ,JK,JL-1) = 0.
-              ZDEL = SQRT((XYHAT(JJ+1)-XYHAT(JJ))*0.5*(XZZ(JI2,JJ,JK+1)-XZZ(JI2,JJ,JK)+XZZ(JI2-1,JJ,JK+1)-XZZ(JI2-1,JJ,JK)))
+              ZDEL = SQRT( XDYHAT(JJ) * 0.5 * (XZZ(JI2,JJ,JK+1)-XZZ(JI2,JJ,JK)+XZZ(JI2-1,JJ,JK+1)-XZZ(JI2-1,JJ,JK)) )
               ZPH1 = PPHI(JI2 ,JJ  ,JK  ,JL)
               ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
               ZVIT1 = ZSIG1*PRUS(JI2,JJ ,JK )
@@ -206,7 +207,7 @@ SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
               JL = 2
               JI2 = JI+1
               ZIBM_FLUX(JI2,JJ,JK,JL-1) = 0.
-              ZDEL = SQRT((XYHAT(JJ+1)-XYHAT(JJ))*0.5*(XZZ(JI2,JJ,JK+1)-XZZ(JI2,JJ,JK)+XZZ(JI2-1,JJ,JK+1)-XZZ(JI2-1,JJ,JK)))
+              ZDEL = SQRT( XDYHAT(JJ) * 0.5 * (XZZ(JI2,JJ,JK+1)-XZZ(JI2,JJ,JK)+XZZ(JI2-1,JJ,JK+1)-XZZ(JI2-1,JJ,JK)) )
               ZPH1 = PPHI(JI2  ,JJ  ,JK  ,JL)
               ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
               ZVIT1 = ZSIG1*PRUS(JI2,JJ ,JK )
@@ -270,7 +271,7 @@ SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
               JL = 3 
               JJ2 = JJ
               ZIBM_FLUX(JI,JJ2,JK,JL-1) = 0.
-              ZDEL = SQRT((XXHAT(JI+1)-XXHAT(JI))*0.5*(XZZ(JI,JJ2,JK+1)-XZZ(JI,JJ2,JK)+XZZ(JI,JJ2-1,JK+1)-XZZ(JI,JJ2-1,JK)))
+              ZDEL = SQRT( XDXHAT(JI) * 0.5 * (XZZ(JI,JJ2,JK+1)-XZZ(JI,JJ2,JK)+XZZ(JI,JJ2-1,JK+1)-XZZ(JI,JJ2-1,JK)) )
               ZPH1 = PPHI(JI  ,JJ2  ,JK  ,JL)
               ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
               ZVIT1 = ZSIG1*PRVS(JI ,JJ2,JK )
@@ -335,7 +336,7 @@ SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
               JL = 3 
               JJ2 = JJ+1
               ZIBM_FLUX(JI,JJ2,JK,JL-1) = 0.
-              ZDEL = SQRT((XXHAT(JI+1)-XXHAT(JI))*0.5*(XZZ(JI,JJ2,JK+1)-XZZ(JI,JJ2,JK)+XZZ(JI,JJ2-1,JK+1)-XZZ(JI,JJ2-1,JK)))
+              ZDEL = SQRT( XDXHAT(JI) * 0.5 * (XZZ(JI,JJ2,JK+1)-XZZ(JI,JJ2,JK)+XZZ(JI,JJ2-1,JK+1)-XZZ(JI,JJ2-1,JK)) )
               ZPH1 = PPHI(JI  ,JJ2  ,JK  ,JL)
               ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
               ZVIT1 = ZSIG1*PRVS(JI ,JJ2,JK )
@@ -401,7 +402,7 @@ SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
               JL = 4 
               JK2 = JK
               ZIBM_FLUX(JI,JJ,JK2,JL-1) = 0.
-              ZDEL = SQRT((XXHAT(JI+1)-XXHAT(JI))*(XYHAT(JJ+1)-XYHAT(JJ)))
+              ZDEL = SQRT( XDXHAT(JI) * XDYHAT(JJ) )
               ZPH1 = PPHI(JI  ,JJ  ,JK2 ,JL)
               ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
               ZVIT1 = ZSIG1*PRWS(JI ,JJ ,JK2)
@@ -467,7 +468,7 @@ SUBROUTINE IBM_BALANCE(PPHI,PVOL,PRUS,PRVS,PRWS,PBAL)
               JL = 4 
               JK2 = JK+1
               ZIBM_FLUX(JI,JJ,JK2,JL-1) = 0.
-              ZDEL = SQRT((XXHAT(JI+1)-XXHAT(JI))*(XYHAT(JJ+1)-XYHAT(JJ)))
+              ZDEL = SQRT( XDXHAT(JI) * XDYHAT(JJ) )
               ZPH1 = PPHI(JI  ,JJ  ,JK2 ,JL)
               ZSIG1 = max(0.,-ZPH1/abs(ZPH1))
               ZVIT1 = ZSIG1*PRWS(JI ,JJ ,JK2)
diff --git a/src/MNH/ibm_detect.f90 b/src/MNH/ibm_detect.f90
index ca4530964ff6f617309a5df187c9c689cbf73d53..b88a80445aad736413f45e64431eae8f3a3c3db1 100644
--- a/src/MNH/ibm_detect.f90
+++ b/src/MNH/ibm_detect.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2019-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.
@@ -75,7 +75,6 @@ SUBROUTINE IBM_DETECT(PPHI)
   ! declaration
   USE MODD_IBM_PARAM_n
   USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZHAT,XZZ
   USE MODD_METRICS_n, ONLY: XDXX,XDYY,XDZZ,XDZX,XDZY
   USE MODD_LBC_n
   USE MODD_CONF, ONLY: NHALO
diff --git a/src/MNH/ibm_generls.f90 b/src/MNH/ibm_generls.f90
index a129d210930de85d7ade5ed783ce6c57c7714ad7..1b768365207207d3f7c07875bc3c44ccdfb61a0d 100644
--- a/src/MNH/ibm_generls.f90
+++ b/src/MNH/ibm_generls.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2021-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2021-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.
@@ -76,7 +76,6 @@ SUBROUTINE IBM_GENERLS(PIBM_FACES,PNORM_FACES,PV1,PV2,PV3,PX_MIN,PY_MIN,PX_MAX,P
   USE MODD_IBM_LSF
   USE MODD_DIM_n, ONLY: NIMAX,NJMAX,NKMAX
   USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT,XUNDEF
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZHAT,XZZ
   USE MODD_METRICS_n, ONLY: XDXX,XDYY,XDZZ
   USE MODD_VAR_ll, ONLY: IP
   USE MODD_CST, ONLY: XMNH_EPSILON 
diff --git a/src/MNH/ibm_idealee.f90 b/src/MNH/ibm_idealee.f90
index e08be780d96d538d079536dc127c04074e860dd8..1c2d6449fec35679e66e1b036ac9773fddce7504 100644
--- a/src/MNH/ibm_idealee.f90
+++ b/src/MNH/ibm_idealee.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2019-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.
@@ -69,7 +69,6 @@ SUBROUTINE IBM_IDEALEE(KNUMB_OBS,PIBM_XYZ,PPHI)
   ! declaration
   USE MODD_IBM_PARAM_n
   USE MODD_DIM_n, ONLY: NIMAX,NJMAX,NKMAX
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ
   USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT
   !
   ! interface
diff --git a/src/MNH/ibm_init_ls.f90 b/src/MNH/ibm_init_ls.f90
index 2d881e1fd564803d1d7f786993fa5efc2fe14042..bdfbff5ad10d7304d4ad3efa355e280305fc8250 100644
--- a/src/MNH/ibm_init_ls.f90
+++ b/src/MNH/ibm_init_ls.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2019-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.
@@ -70,7 +70,7 @@ SUBROUTINE IBM_INIT_LS(PPHI)
   USE MODD_IBM_PARAM_n, ONLY: XIBM_EPSI,XIBM_IEPS
   USE MODD_IBM_LSF, ONLY: LIBM_LSF,CIBM_TYPE,NIBM_SMOOTH,XIBM_SMOOTH
   USE MODD_VAR_ll, ONLY: IP
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZHAT,XZZ
+  USE MODD_GRID_n,      ONLY: XZZ
   USE MODD_PARAMETERS, ONLY: XUNDEF,JPHEXT,JPVEXT
   !
   ! interface
diff --git a/src/MNH/ibm_mixinglength.f90 b/src/MNH/ibm_mixinglength.f90
index 14bb0dd89b6effcd7d00516aea80b62a9b222b78..c4c1a2c10bd3d50e77e484a3b2a0d833e537b05b 100644
--- a/src/MNH/ibm_mixinglength.f90
+++ b/src/MNH/ibm_mixinglength.f90
@@ -75,7 +75,7 @@ SUBROUTINE IBM_MIXINGLENGTH(PLM,PLEPS,PMU,PHI,PTKE)
   USE MODD_REF_n, ONLY: XRHODJ,XRHODREF
   USE MODD_CTURB
   USE MODD_CST
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ
+  USE MODD_GRID_n, ONLY: XZZ
   !
   ! interface
   !
diff --git a/src/MNH/ibm_smooth_ls.f90 b/src/MNH/ibm_smooth_ls.f90
index 96144123454422ed14e94f1ae592098171937ca9..7e06c354df84832fcb6a2ffaf27490fecd3d71b2 100644
--- a/src/MNH/ibm_smooth_ls.f90
+++ b/src/MNH/ibm_smooth_ls.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2019-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.
@@ -74,7 +74,7 @@ SUBROUTINE IBM_SMOOTH_LS(KIBM_SMOOTH,PIBM_SMOOTH,PPHI)
   USE MODD_IBM_PARAM_n    
   USE MODD_IBM_LSF   
   USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT   
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ 
+  USE MODD_GRID_n,      ONLY: XDXHAT, XDYHAT
   USE MODD_METRICS_n, ONLY: XDXX,XDYY,XDZZ,XDZX,XDZY
   USE MODD_ARGSLIST_ll, ONLY: LIST_ll
   USE MODD_VAR_ll, ONLY: IP
@@ -131,8 +131,8 @@ SUBROUTINE IBM_SMOOTH_LS(KIBM_SMOOTH,PIBM_SMOOTH,PPHI)
   !
   IKE = IKU - JPVEXT
   IKB =   1 + JPVEXT
-  ZREF =(1.e-2)*((XXHAT(IIB+1)-XXHAT(IIB))*(XYHAT(IJB+1)-XYHAT(IJB)))**0.5
-  ZREF3=((XXHAT(IIB+1)-XXHAT(IIB))*(XYHAT(IJB+1)-XYHAT(IJB)))**0.5
+  ZREF =(1.e-2)*( XDXHAT(IIB) * XDYHAT(IJB) )**0.5
+  ZREF3=        ( XDXHAT(IIB) * XDYHAT(IJB) )**0.5
   !
   ! Boundary symmetry
   !
diff --git a/src/MNH/ibm_volume.f90 b/src/MNH/ibm_volume.f90
index af4012a42b9cba352d527c22179f6de698ff59dc..ec734278ffa65c8cff5eef79aa8247ebe816a099 100644
--- a/src/MNH/ibm_volume.f90
+++ b/src/MNH/ibm_volume.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2019-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.
@@ -66,7 +66,6 @@ SUBROUTINE IBM_VOLUME(PPHI,PVOL)
   !
   ! declaration
   USE MODD_IBM_PARAM_n        
-  USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZHAT,XZZ 
   USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT
   USE MODD_LBC_n
   USE MODD_LUNIT_n, ONLY: TLUOUT
diff --git a/src/MNH/ini_dynamics.f90 b/src/MNH/ini_dynamics.f90
index e4d00f5bb6aa29d04e9ce3cdb3cf10e0ae4187cb..add1cc9e389d189ec736e64cc2c1c4f87e6d4758 100644
--- a/src/MNH/ini_dynamics.f90
+++ b/src/MNH/ini_dynamics.f90
@@ -1,6 +1,6 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 !     ########################
@@ -8,7 +8,7 @@
 !     ########################
 INTERFACE
 SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ,                   &
-               PDXHAT,PDYHAT,PZHAT,HLBCX,HLBCY,PTSTEP,                       &
+               PDXHAT,PDYHAT,PZHAT,PZHATM,HLBCX,HLBCY,PTSTEP,                &
                OVE_RELAX,OVE_RELAX_GRD,OHORELAX_UVWTH,OHORELAX_RV,           &
                OHORELAX_RC,OHORELAX_RR,OHORELAX_RI,OHORELAX_RS,OHORELAX_RG,  &
                OHORELAX_RH,OHORELAX_TKE,OHORELAX_SV,                         &
@@ -41,7 +41,8 @@ REAL, DIMENSION(:,:),   INTENT(IN)        :: PMAP      ! Map factor
 REAL, DIMENSION(:,:,:), INTENT(IN)        :: PZZ       ! height 
 REAL, DIMENSION(:),     INTENT(IN)        :: PDXHAT     ! Stretching in x direction
 REAL, DIMENSION(:),     INTENT(IN)        :: PDYHAT     ! Stretching in y direction
-REAL, DIMENSION(:),     INTENT(IN)        :: PZHAT     ! Gal-Chen Height   
+REAL, DIMENSION(:),     INTENT(IN)        :: PZHAT     ! Gal-Chen Height
+REAL, DIMENSION(:),     INTENT(IN)        :: PZHATM    ! ... at mass points
 CHARACTER(LEN=4), DIMENSION(:), INTENT(IN) :: HLBCX    ! x-direction LBC type
 CHARACTER(LEN=4), DIMENSION(:), INTENT(IN) :: HLBCY    ! y-direction LBC type
 LOGICAL,                INTENT(IN)        :: OVE_RELAX ! logical
@@ -179,7 +180,7 @@ END INTERFACE
 END MODULE MODI_INI_DYNAMICS
 !     ######################################################################
 SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ,                   &
-               PDXHAT,PDYHAT,PZHAT,HLBCX,HLBCY,PTSTEP,                       &
+               PDXHAT,PDYHAT,PZHAT,PZHATM,HLBCX,HLBCY,PTSTEP,                &
                OVE_RELAX,OVE_RELAX_GRD,OHORELAX_UVWTH,OHORELAX_RV,           &
                OHORELAX_RC,OHORELAX_RR,OHORELAX_RI,OHORELAX_RS,OHORELAX_RG,  &
                OHORELAX_RH,OHORELAX_TKE,OHORELAX_SV,                         &
@@ -312,6 +313,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN)        :: PZZ       ! height
 REAL, DIMENSION(:),     INTENT(IN)        :: PDXHAT     ! Stretching in x direction
 REAL, DIMENSION(:),     INTENT(IN)        :: PDYHAT     ! Stretching in y direction
 REAL, DIMENSION(:),     INTENT(IN)        :: PZHAT     ! Gal-Chen Height   
+REAL, DIMENSION(:),     INTENT(IN)        :: PZHATM    ! ... at mass points
 CHARACTER(LEN=4), DIMENSION(:), INTENT(IN) :: HLBCX    ! x-direction LBC type
 CHARACTER(LEN=4), DIMENSION(:), INTENT(IN) :: HLBCY    ! y-direction LBC type
 LOGICAL,                INTENT(IN)        :: OVE_RELAX ! logical
@@ -540,7 +542,7 @@ IF (GHORELAX .OR. OVE_RELAX.OR.OVE_RELAX_GRD) THEN
      OHORELAX_SVCHEM, OHORELAX_SVAER, OHORELAX_SVDST, OHORELAX_SVSLT, &
      OHORELAX_SVPP, OHORELAX_SVCS, OHORELAX_SVCHIC,OHORELAX_SVSNW,    &
      PALKTOP,PALKGRD, PALZBOT,PALZBAS,                                &
-     PZZ, PZHAT, PTSTEP,                                              &
+     PZZ, PZHAT, PZHATM, PTSTEP,                                      &
      PRIMKMAX,KRIMX,KRIMY,                                            &
      PALK, PALKW, KALBOT,                                             &
      PALKBAS, PALKWBAS, KALBAS,                                       &
diff --git a/src/MNH/ini_lg.f90 b/src/MNH/ini_lg.f90
index 0ce0d7b1b535d813371881bfc2daf5d6b968aae2..8f5428cae1f9421a238f4c35a74b69ff17b81544 100644
--- a/src/MNH/ini_lg.f90
+++ b/src/MNH/ini_lg.f90
@@ -1,24 +1,19 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2001-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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 init 2006/05/18 13:07:25
-!-----------------------------------------------------------------
 !     ##################
       MODULE MODI_INI_LG
 !     ##################
 INTERFACE
 !
-      SUBROUTINE INI_LG(PXHAT,PYHAT,PZZ,PSVT,PLBXSVM,PLBYSVM)
+      SUBROUTINE INI_LG( PXHATM, PYHATM, PZZ, PSVT, PLBXSVM, PLBYSVM )
 !
-REAL,DIMENSION(:),      INTENT(IN) :: PXHAT,PYHAT ! Positions x,y in the cartesian plane
-REAL,DIMENSION(:,:,:), INTENT(IN)  :: PZZ         ! True altitude of the w grid-point
-REAL,DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT        ! scalar var. at t
-REAL,DIMENSION(:,:,:,:), INTENT(INOUT) :: PLBXSVM,PLBYSVM  ! LB in x and y-dir.
+REAL, DIMENSION(:),       INTENT(IN)    :: PXHATM, PYHATM    ! Positions x,y in the cartesian plane at mass points
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PZZ               ! True altitude of the w grid-point
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT              ! scalar var. at t
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PLBXSVM, PLBYSVM  ! LB in x and y-dir.
 !
 END SUBROUTINE INI_LG
 !
@@ -28,9 +23,9 @@ END MODULE MODI_INI_LG
 !
 !
 !
-!     ############################################################
-      SUBROUTINE INI_LG(PXHAT,PYHAT,PZZ,PSVT,PLBXSVM,PLBYSVM)
-!     ############################################################
+!     ################################################################
+      SUBROUTINE INI_LG( PXHATM, PYHATM, PZZ, PSVT, PLBXSVM, PLBYSVM )
+!     ################################################################
 !
 !!****  *INI_LG* - routine to initialize lagrangian variables
 !!
@@ -78,10 +73,10 @@ IMPLICIT NONE
 !
 !*       0.1   declarations of argument
 !
-REAL,DIMENSION(:),      INTENT(IN) :: PXHAT,PYHAT ! Positions x,y in the cartesian plane
-REAL,DIMENSION(:,:,:), INTENT(IN)  :: PZZ         ! True altitude of the w grid-point
-REAL,DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT        ! scalar var. at t
-REAL,DIMENSION(:,:,:,:), INTENT(INOUT) :: PLBXSVM,PLBYSVM  ! LB in x and y-dir.
+REAL, DIMENSION(:),       INTENT(IN)    :: PXHATM, PYHATM    ! Positions x,y in the cartesian plane at mass points
+REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PZZ               ! True altitude of the w grid-point
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVT              ! scalar var. at t
+REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PLBXSVM, PLBYSVM  ! LB in x and y-dir.
 !
 !
 !*       0.2   declarations of local variables
@@ -103,7 +98,7 @@ IKU=SIZE(PZZ,3)
 DO JK=1,IKU
   DO JJ=1,IJU
     DO JI=1,IIU-1
-      PSVT(JI,JJ,JK,NSV_LGBEG)=0.5*(PXHAT(JI)+PXHAT(JI+1))
+      PSVT(JI,JJ,JK,NSV_LGBEG)=PXHATM(JI)
     END DO
     PSVT(IIU,JJ,JK,NSV_LGBEG)=2.*PSVT(IIU-1,JJ,JK,NSV_LGBEG)-PSVT(IIU-2,JJ,JK,NSV_LGBEG)
   END DO
@@ -112,7 +107,7 @@ END DO
 DO JK=1,IKU
   DO JI=1,IIU
     DO JJ=1,IJU-1
-      PSVT(JI,JJ,JK,NSV_LGBEG+1)=0.5*(PYHAT(JJ)+PYHAT(JJ+1))
+      PSVT(JI,JJ,JK,NSV_LGBEG+1)=PYHATM(JJ)
     END DO
     PSVT(JI,IJU,JK,NSV_LGBEG+1)=2.*PSVT(JI,IJU-1,JK,NSV_LGBEG+1)-PSVT(JI,IJU-2,JK,NSV_LGBEG+1)
   END DO
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 11dad29e3b7aae197cde18308ab71eaa3e3a58b3..282936099e6da5e136f458cf944751555a23c8d4 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -1001,10 +1001,13 @@ ALLOCATE(XXHAT(IIU))
 ALLOCATE(XDXHAT(IIU))
 ALLOCATE(XYHAT(IJU))
 ALLOCATE(XDYHAT(IJU))
+ALLOCATE(XXHATM(IIU))
+ALLOCATE(XYHATM(IJU))
 ALLOCATE(XZS(IIU,IJU))
 ALLOCATE(XZSMT(IIU,IJU))
 ALLOCATE(XZZ(IIU,IJU,IKU))
 ALLOCATE(XZHAT(IKU))
+ALLOCATE(XZHATM(IKU))
 ALLOCATE(XDIRCOSZW(IIU,IJU))
 ALLOCATE(XDIRCOSXW(IIU,IJU))
 ALLOCATE(XDIRCOSYW(IIU,IJU))
@@ -1852,13 +1855,13 @@ END IF
 !*       7.    INITIALIZE GRIDS AND METRIC COEFFICIENTS
 !              ----------------------------------------
 !
-CALL SET_GRID(KMI,TPINIFILE,IKU,NIMAX_ll,NJMAX_ll,                       &
-              XTSTEP,XSEGLEN,                                            &
-              XLONORI,XLATORI,XLON,XLAT,                                 &
-              XXHAT,XYHAT,XDXHAT,XDYHAT, XMAP,                           &
-              XZS,XZZ,XZHAT,XZTOP,LSLEVE,XLEN1,XLEN2,XZSMT,              &
-              ZJ,                                                        &
-              TDTMOD,TDTCUR,NSTOP,NBAK_NUMB,NOUT_NUMB,TBACKUPN,TOUTPUTN)
+CALL SET_GRID( KMI, TPINIFILE, IKU, NIMAX_ll, NJMAX_ll,                        &
+               XTSTEP, XSEGLEN,                                                &
+               XLONORI, XLATORI, XLON, XLAT,                                   &
+               XXHAT, XYHAT, XDXHAT, XDYHAT, XXHATM, XYHATM, XMAP,             &
+               XZS, XZZ, XZHAT, XZHATM, XZTOP, LSLEVE, XLEN1, XLEN2, XZSMT,    &
+               ZJ,                                                             &
+               TDTMOD, TDTCUR, NSTOP, NBAK_NUMB, NOUT_NUMB, TBACKUPN, TOUTPUTN )
 !
 CALL METRICS(XMAP,XDXHAT,XDYHAT,XZZ,XDXX,XDYY,XDZX,XDZY,XDZZ)
 !
@@ -1938,10 +1941,10 @@ CALL READ_FIELD(KMI,TPINIFILE,IIU,IJU,IKU,                                    &
 !              ---------------------------
 !
 !
-CALL SET_REF(KMI,TPINIFILE,                         &
-             XZZ,XZHAT,ZJ,XDXX,XDYY,CLBCX,CLBCY,    &
-             XREFMASS,XMASS_O_PHI0,XLINMASS,        &
-             XRHODREF,XTHVREF,XRVREF,XEXNREF,XRHODJ )
+CALL SET_REF( KMI, TPINIFILE,                            &
+              XZZ, XZHATM, ZJ, XDXX, XDYY, CLBCX, CLBCY, &
+              XREFMASS, XMASS_O_PHI0, XLINMASS,          &
+              XRHODREF, XTHVREF, XRVREF, XEXNREF, XRHODJ )
 !
 !-------------------------------------------------------------------------------
 !
@@ -2187,7 +2190,7 @@ END IF
 !               -------------------------------
 !
 IF (LLG .AND. LINIT_LG .AND. CPROGRAM=='MESONH') &
-  CALL INI_LG(XXHAT,XYHAT,XZZ,XSVT,XLBXSVM,XLBYSVM)
+  CALL INI_LG( XXHATM, XYHATM, XZZ, XSVT, XLBXSVM, XLBYSVM )
 
 !
 !-------------------------------------------------------------------------------
@@ -2196,7 +2199,7 @@ IF (LLG .AND. LINIT_LG .AND. CPROGRAM=='MESONH') &
 !               ------------------------------------------
 !
 CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT,            &
-             XZHAT,CLBCX,CLBCY,XTSTEP,                                        &
+             XZHAT,XZHATM,CLBCX,CLBCY,XTSTEP,                                 &
              LVE_RELAX,LVE_RELAX_GRD,LHORELAX_UVWTH,LHORELAX_RV,              &
              LHORELAX_RC,LHORELAX_RR,LHORELAX_RI,LHORELAX_RS,LHORELAX_RG,     &
              LHORELAX_RH,LHORELAX_TKE,LHORELAX_SV,                            &
diff --git a/src/MNH/ini_posprofilern.f90 b/src/MNH/ini_posprofilern.f90
index ad634021133d7a805054226b2ab90c1c4387bf78..966c343cf93cf3cc28790b6ac2926a90754a1434 100644
--- a/src/MNH/ini_posprofilern.f90
+++ b/src/MNH/ini_posprofilern.f90
@@ -62,7 +62,7 @@ USE MODD_ALLPROFILER_n
 USE MODD_CONF,           ONLY: LCARTESIAN
 USE MODD_DYN,            ONLY: XSEGLEN
 USE MODD_DYN_n,          ONLY: DYN_MODEL, XTSTEP
-USE MODD_GRID_n,         ONLY: XXHAT, XYHAT
+USE MODD_GRID_n,         ONLY: XXHAT, XXHATM, XYHAT, XYHATM
 USE MODD_PARAMETERS,     ONLY: JPHEXT, JPVEXT
 USE MODD_PROFILER_n,     ONLY: LPROFILER, NUMBPROFILER_LOC, TPROFILERS, TPROFILERS_TIME
 USE MODD_TYPE_STATPROF,  ONLY: TPROFILERDATA
@@ -95,8 +95,6 @@ LOGICAL :: GINSIDE                          ! True if profiler is inside physica
 LOGICAL :: GPRESENT                         ! True if profiler is present on the current process
 REAL    :: ZXHATM_PHYS_MIN, ZYHATM_PHYS_MIN ! Minimum X coordinate of mass points in the physical domain
 REAL    :: ZXHATM_PHYS_MAX, ZYHATM_PHYS_MAX ! Minimum X coordinate of mass points in the physical domain
-REAL, DIMENSION(SIZE(XXHAT)) :: ZXHATM      ! mass point coordinates
-REAL, DIMENSION(SIZE(XYHAT)) :: ZYHATM      ! mass point coordinates
 REAL, DIMENSION(:), POINTER  :: ZXHAT_GLOB
 REAL, DIMENSION(:), POINTER  :: ZYHAT_GLOB
 TYPE(TPROFILERDATA)          :: TZPROFILER
@@ -126,13 +124,6 @@ IF ( CFILE_PROF /= "NO_INPUT_CSV" .OR. NNUMB_PROF > 0 ) THEN
   CALL GATHERALL_FIELD_ll( 'XX', XXHAT, ZXHAT_GLOB, IERR )
   CALL GATHERALL_FIELD_ll( 'YY', XYHAT, ZYHAT_GLOB, IERR )
 
-  ! Interpolations of model variables to mass points
-  ZXHATM(1:IIU-1) = 0.5 * XXHAT(1:IIU-1) + 0.5 * XXHAT(2:IIU  )
-  ZXHATM(  IIU  ) = 1.5 * XXHAT(  IIU  ) - 0.5 * XXHAT(  IIU-1)
-
-  ZYHATM(1:IJU-1) = 0.5 * XYHAT(1:IJU-1) + 0.5 * XYHAT(2:IJU  )
-  ZYHATM(  IJU  ) = 1.5 * XYHAT(  IJU  ) - 0.5 * XYHAT(  IJU-1)
-
   ZXHATM_PHYS_MIN = 0.5 * ( ZXHAT_GLOB(1+JPHEXT) + ZXHAT_GLOB(2+JPHEXT) )
   ZXHATM_PHYS_MAX = 0.5 * ( ZXHAT_GLOB(UBOUND(ZXHAT_GLOB,1)-JPHEXT) + ZXHAT_GLOB(UBOUND(ZXHAT_GLOB,1)-JPHEXT+1) )
   ZYHATM_PHYS_MIN = 0.5 * ( ZYHAT_GLOB(1+JPHEXT) + ZYHAT_GLOB(2+JPHEXT) )
@@ -159,7 +150,7 @@ IF (CFILE_PROF=="NO_INPUT_CSV") THEN
       TZPROFILER%XZ    = XZ_PROF(JI)
       TZPROFILER%CNAME = CNAME_PROF(JI)
 
-      CALL STATPROF_POSITION( TZPROFILER, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,                 &
+      CALL STATPROF_POSITION( TZPROFILER, ZXHAT_GLOB, ZYHAT_GLOB, XXHATM, XYHATM,                 &
                               ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX, ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
                               GINSIDE, GPRESENT                                                   )
 
@@ -173,7 +164,7 @@ IF (CFILE_PROF=="NO_INPUT_CSV") THEN
   END IF
 ELSE
   !Treat CSV datafile
-  CALL STATPROF_CSV_READ( TZPROFILER, CFILE_PROF, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,    &
+  CALL STATPROF_CSV_READ( TZPROFILER, CFILE_PROF, ZXHAT_GLOB, ZYHAT_GLOB, XXHATM, XYHATM,    &
                           ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX,ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
                           INUMBPROF                                                          )
 END IF
diff --git a/src/MNH/ini_size_spawn.f90 b/src/MNH/ini_size_spawn.f90
index abab91578153c475a65cd826bb4f559cad97d426..3695fe140491677c653e6c4af6d6909873b28bb0 100644
--- a/src/MNH/ini_size_spawn.f90
+++ b/src/MNH/ini_size_spawn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1999-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1999-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.
@@ -287,7 +287,7 @@ IF (LEN_TRIM(CDOMAIN)>0) THEN
   CALL READ_HGRID(2,TZDOMAIN,YMY_NAME,YDAD_NAME,YSTORAGE_TYPE)
   CALL IO_File_close(TZDOMAIN)
   CALL RETRIEVE1_NEST_INFO_n(1,2,NXOR,NYOR,NXSIZE,NYSIZE,NDXRATIO,NDYRATIO)
-  DEALLOCATE(XZS,XZSMT,XXHAT,XYHAT)
+  DEALLOCATE( XZS, XZSMT, XXHAT, XYHAT, XXHATM, XYHATM)
 !
 END IF
 !
diff --git a/src/MNH/ini_spectren.f90 b/src/MNH/ini_spectren.f90
index 1067f2ceffdbf2de2b8de02dcf94039a76a01378..7f62014dd039f0f6aa81cde9cd70307df5fbb601 100644
--- a/src/MNH/ini_spectren.f90
+++ b/src/MNH/ini_spectren.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2015-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2015-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.
@@ -262,10 +262,13 @@ ALLOCATE(XXHAT(IIU))
 ALLOCATE(XDXHAT(IIU))
 ALLOCATE(XYHAT(IJU))
 ALLOCATE(XDYHAT(IJU))
+ALLOCATE(XXHATM(IIU))
+ALLOCATE(XYHATM(IJU))
 ALLOCATE(XZS(IIU,IJU))
 ALLOCATE(XZSMT(IIU,IJU))
 ALLOCATE(XZZ(IIU,IJU,IKU))
 ALLOCATE(XZHAT(IKU))
+ALLOCATE(XZHATM(IKU))
 !
 ALLOCATE(XDXX(IIU,IJU,IKU))
 ALLOCATE(XDYY(IIU,IJU,IKU))
@@ -682,13 +685,13 @@ CALL INI_BIKHARDT_n (NDXRATIO_ALL(KMI),NDYRATIO_ALL(KMI),KMI)
 !*       6.    INITIALIZE GRIDS AND METRIC COEFFICIENTS
 !              ----------------------------------------
 !
-CALL SET_GRID(KMI,TPINIFILE,IKU,NIMAX_ll,NJMAX_ll,                       &
-              XTSTEP,XSEGLEN,                                            &
-              XLONORI,XLATORI,XLON,XLAT,                                 &
-              XXHAT,XYHAT,XDXHAT,XDYHAT, XMAP,                           &
-              XZS,XZZ,XZHAT,XZTOP,LSLEVE,XLEN1,XLEN2,XZSMT,              &
-              ZJ,                                                        &
-              TDTMOD,TDTCUR,NSTOP,NBAK_NUMB,NOUT_NUMB,TBACKUPN,TOUTPUTN)
+CALL SET_GRID( KMI, TPINIFILE, IKU, NIMAX_ll, NJMAX_ll,                        &
+               XTSTEP, XSEGLEN,                                                &
+               XLONORI, XLATORI, XLON, XLAT,                                   &
+               XXHAT, XYHAT, XDXHAT, XDYHAT, XXHATM, XYHATM, XMAP,             &
+               XZS, XZZ, XZHAT, XZHATM, XZTOP, LSLEVE, XLEN1, XLEN2, XZSMT,    &
+               ZJ,                                                             &
+               TDTMOD, TDTCUR, NSTOP, NBAK_NUMB, NOUT_NUMB, TBACKUPN, TOUTPUTN )
 !
 CALL METRICS(XMAP,XDXHAT,XDYHAT,XZZ,XDXX,XDYY,XDZX,XDZY,XDZZ)
 !
@@ -743,10 +746,10 @@ END IF
 !              ---------------------------
 !
 !
-CALL SET_REF(KMI,TPINIFILE,                        &
-             XZZ,XZHAT,ZJ,XDXX,XDYY,CLBCX,CLBCY,   &
-             XREFMASS,XMASS_O_PHI0,XLINMASS,       &
-             XRHODREF,XTHVREF,XRVREF,XEXNREF,XRHODJ)
+CALL SET_REF( KMI, TPINIFILE,                            &
+              XZZ, XZHATM, ZJ, XDXX, XDYY, CLBCX, CLBCY, &
+              XREFMASS, XMASS_O_PHI0, XLINMASS,          &
+              XRHODREF, XTHVREF, XRVREF, XEXNREF, XRHODJ )
 !-------------------------------------------------------------------------------
 !
 !*       11.    INITIALIZE THE SOURCE OF TOTAL DRY MASS Md
@@ -894,7 +897,7 @@ ALLOCATE(XALKBAS(0))
 ALLOCATE(XALKWBAS(0))
 !
 CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT,            &
-             XZHAT,CLBCX,CLBCY,XTSTEP,                                        &
+             XZHAT,XZHATM,CLBCX,CLBCY,XTSTEP,                                 &
              LVE_RELAX,LVE_RELAX_GRD,LHORELAX_UVWTH,LHORELAX_RV,              &
              LHORELAX_RC,LHORELAX_RR,LHORELAX_RI,LHORELAX_RS,LHORELAX_RG,     &
              LHORELAX_RH,LHORELAX_TKE,LHORELAX_SV,                            &
diff --git a/src/MNH/ini_surfstationn.f90 b/src/MNH/ini_surfstationn.f90
index 5452a5f0e005ca50e8f0007a3d342b8a81aaa614..94332d77c09dcd643ea5da119bf0e1bf37ff1943 100644
--- a/src/MNH/ini_surfstationn.f90
+++ b/src/MNH/ini_surfstationn.f90
@@ -63,7 +63,7 @@ USE MODD_ALLSTATION_n
 USE MODD_CONF,           ONLY: LCARTESIAN
 USE MODD_DYN,            ONLY: XSEGLEN
 USE MODD_DYN_n,          ONLY: DYN_MODEL, XTSTEP
-USE MODD_GRID_n,         ONLY: XXHAT, XYHAT
+USE MODD_GRID_n,         ONLY: XXHAT, XXHATM, XYHAT, XYHATM
 USE MODD_PARAMETERS,     ONLY: JPHEXT, JPVEXT
 USE MODD_STATION_n
 USE MODD_TYPE_STATPROF
@@ -96,8 +96,6 @@ LOGICAL :: GINSIDE                          ! True if station is inside physical
 LOGICAL :: GPRESENT                         ! True if station is present on the current process
 REAL    :: ZXHATM_PHYS_MIN, ZYHATM_PHYS_MIN ! Minimum X coordinate of mass points in the physical domain
 REAL    :: ZXHATM_PHYS_MAX, ZYHATM_PHYS_MAX ! Minimum X coordinate of mass points in the physical domain
-REAL, DIMENSION(SIZE(XXHAT)) :: ZXHATM      ! mass point coordinates
-REAL, DIMENSION(SIZE(XYHAT)) :: ZYHATM      ! mass point coordinates
 REAL, DIMENSION(:), POINTER  :: ZXHAT_GLOB
 REAL, DIMENSION(:), POINTER  :: ZYHAT_GLOB
 TYPE(TSTATIONDATA)           :: TZSTATION
@@ -127,13 +125,6 @@ IF ( CFILE_STAT /= "NO_INPUT_CSV" .OR. NNUMB_STAT > 0 ) THEN
   CALL GATHERALL_FIELD_ll( 'XX', XXHAT, ZXHAT_GLOB, IERR )
   CALL GATHERALL_FIELD_ll( 'YY', XYHAT, ZYHAT_GLOB, IERR )
 
-  ! Interpolations of model variables to mass points
-  ZXHATM(1:IIU-1) = 0.5 * XXHAT(1:IIU-1) + 0.5 * XXHAT(2:IIU  )
-  ZXHATM(  IIU  ) = 1.5 * XXHAT(  IIU  ) - 0.5 * XXHAT(  IIU-1)
-
-  ZYHATM(1:IJU-1) = 0.5 * XYHAT(1:IJU-1) + 0.5 * XYHAT(2:IJU  )
-  ZYHATM(  IJU  ) = 1.5 * XYHAT(  IJU  ) - 0.5 * XYHAT(  IJU-1)
-
   ZXHATM_PHYS_MIN = 0.5 * ( ZXHAT_GLOB(1+JPHEXT) + ZXHAT_GLOB(2+JPHEXT) )
   ZXHATM_PHYS_MAX = 0.5 * ( ZXHAT_GLOB(UBOUND(ZXHAT_GLOB,1)-JPHEXT) + ZXHAT_GLOB(UBOUND(ZXHAT_GLOB,1)-JPHEXT+1) )
   ZYHATM_PHYS_MIN = 0.5 * ( ZYHAT_GLOB(1+JPHEXT) + ZYHAT_GLOB(2+JPHEXT) )
@@ -160,7 +151,7 @@ IF (CFILE_STAT=="NO_INPUT_CSV") THEN
       TZSTATION%XZ    = XZ_STAT(JI)
       TZSTATION%CNAME = CNAME_STAT(JI)
 
-      CALL STATPROF_POSITION( TZSTATION, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,                 &
+      CALL STATPROF_POSITION( TZSTATION, ZXHAT_GLOB, ZYHAT_GLOB, XXHATM, XYHATM,                 &
                              ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX, ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
                              GINSIDE, GPRESENT                                                   )
 
@@ -174,7 +165,7 @@ IF (CFILE_STAT=="NO_INPUT_CSV") THEN
   END IF
 ELSE
   !Treat CSV datafile
-  CALL STATPROF_CSV_READ( TZSTATION, CFILE_STAT, ZXHAT_GLOB, ZYHAT_GLOB, ZXHATM, ZYHATM,     &
+  CALL STATPROF_CSV_READ( TZSTATION, CFILE_STAT, ZXHAT_GLOB, ZYHAT_GLOB, XXHATM, XYHATM,     &
                           ZXHATM_PHYS_MIN, ZXHATM_PHYS_MAX,ZYHATM_PHYS_MIN, ZYHATM_PHYS_MAX, &
                           INUMBSTAT                                                          )
 END IF
diff --git a/src/MNH/modd_gridn.f90 b/src/MNH/modd_gridn.f90
index 055d3c88f76b4634b987607db83365a12e58710f..6437cfc8bcff5cdbfcfc63c283990283b7c8a06b 100644
--- a/src/MNH/modd_gridn.f90
+++ b/src/MNH/modd_gridn.f90
@@ -1,6 +1,6 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 !     ##################
@@ -36,6 +36,7 @@
 !!      V. Ducrocq   13/08/98  // : add XLATOR_ll and XLONOR_ll       
 !!      V. Masson   nov 2004  supress XLATOR,XLONOR,XLATOR_ll,XLONOR_ll
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!  P. Wautelet 02/09/2022: add XXHATM, XYHATM and XZHATM
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -47,6 +48,8 @@ IMPLICIT NONE
 REAL, DIMENSION(:,:),  POINTER :: XLON=>NULL(),XLAT=>NULL() ! Longitude and latitude  
 REAL, DIMENSION(:),    POINTER :: XXHAT=>NULL()             ! Position x in the conformal or cartesian plane
 REAL, DIMENSION(:),    POINTER :: XYHAT=>NULL()             ! Position y in the conformal or cartesian plane
+REAL, DIMENSION(:),    POINTER :: XXHATM=>NULL()            ! Position x in the conformal or cartesian plane at mass points
+REAL, DIMENSION(:),    POINTER :: XYHATM=>NULL()            ! Position y in the conformal or cartesian plane at mass points
 REAL, DIMENSION(:),    POINTER :: XDXHAT=>NULL()            ! horizontal stretching in x
 REAL, DIMENSION(:),    POINTER :: XDYHAT=>NULL()            ! horizontal stretching in y
 REAL, DIMENSION(:,:),  POINTER :: XMAP=>NULL()              ! Map factor 
@@ -54,6 +57,7 @@ REAL, DIMENSION(:,:),  POINTER :: XZS=>NULL()               ! orography
 REAL, DIMENSION(:,:,:),POINTER :: XZZ=>NULL()               ! height z 
 REAL,                  POINTER :: XZTOP=>NULL()             ! model top (m)
 REAL, DIMENSION(:),    POINTER :: XZHAT=>NULL()             ! height level without orography
+REAL, DIMENSION(:),    POINTER :: XZHATM=>NULL()            ! height level without orography at mass points
 REAL, DIMENSION(:,:),  POINTER :: XDIRCOSXW=>NULL(),XDIRCOSYW=>NULL(),XDIRCOSZW=>NULL() ! director cosinus of the normal
                                                                                         ! to the ground surface
 REAL, DIMENSION(:,:),  POINTER  :: XCOSSLOPE=>NULL()         ! cosinus of the angle between i and the slope vector
diff --git a/src/MNH/modd_spawn.f90 b/src/MNH/modd_spawn.f90
index 8d432e588f5fbce395ee2385fb50883f8eb5f41e..6efedc643dbbef8fdf2a9e5e52c8f72c68afcaf7 100644
--- a/src/MNH/modd_spawn.f90
+++ b/src/MNH/modd_spawn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1999-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.
@@ -66,6 +66,7 @@ CHARACTER (LEN=28) :: CDADSPAFILE ! DAD fm-file for spawning file
 REAL,DIMENSION(:),    SAVE,POINTER :: XXHAT1  => NULL()
 REAL,DIMENSION(:),    SAVE,POINTER :: XYHAT1  => NULL()
 REAL,DIMENSION(:),    SAVE,POINTER :: XZHAT1  => NULL()
+REAL,DIMENSION(:),    SAVE,POINTER :: XZHATM1 => NULL()
 REAL,                 SAVE,POINTER :: XZTOP1  => NULL()
 REAL,DIMENSION(:,:),  SAVE,POINTER :: XZS1    => NULL()
 REAL,DIMENSION(:,:),  SAVE,POINTER :: XZSMT1  => NULL()
diff --git a/src/MNH/modd_sub_elecn.f90 b/src/MNH/modd_sub_elecn.f90
index d25df3084adf326595431c9d9e122baa7a90aa1e..1f5b6b9405c83f091bcb41bcf5bb6f5270b1d820 100644
--- a/src/MNH/modd_sub_elecn.f90
+++ b/src/MNH/modd_sub_elecn.f90
@@ -1,13 +1,8 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2011-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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 modd 2006/06/27 12:43:28
-!-----------------------------------------------------------------
 !!    ############################ 
       MODULE MODD_SUB_ELEC_n
 !!    ############################ 
@@ -32,6 +27,7 @@
 !!    MODIFICATIONS
 !!    -------------
 !!      Original    07/11
+!  P. Wautelet 31/08/2022: remove ZXMASS and ZYMASS (use XXHATM and XYHATM instead)
 !-------------------------------------------------------------------------------
 !
 !*       0.   DECLARATIONS
@@ -50,8 +46,6 @@ TYPE SUB_ELEC_t
   INTEGER , DIMENSION(:), POINTER :: ISNBSEG=>NULL() ! Number of flash segments
   INTEGER , DIMENSION(:), POINTER :: ISTCOUNT_NUMBER=>NULL() ! Temporal loop number of the flash
   INTEGER , DIMENSION(:), POINTER :: ISTYPE=>NULL() ! flash type :IC, CGN or CGP
-  REAL    , DIMENSION(:), POINTER :: ZXMASS=>NULL() ! Coord. at mass points
-  REAL    , DIMENSION(:), POINTER :: ZYMASS=>NULL() ! Coord. at mass points
   REAL    , DIMENSION(:,:,:), POINTER :: ZZMASS=>NULL() ! Coord. at mass points
   REAL    , DIMENSION(:,:,:), POINTER :: ZPRES_COEF=>NULL() ! Pressure effect for E
   REAL    , DIMENSION(:,:,:), POINTER :: ZSCOORD_SEG=>NULL() ! Global coordinates of segments
@@ -71,8 +65,6 @@ TYPE(SUB_ELEC_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: SUB_ELEC_MODEL
   INTEGER , DIMENSION(:), POINTER :: ISNBSEG=>NULL() 
   INTEGER , DIMENSION(:), POINTER :: ISTCOUNT_NUMBER=>NULL() 
   INTEGER , DIMENSION(:), POINTER :: ISTYPE=>NULL() 
-  REAL    , DIMENSION(:), POINTER :: ZXMASS=>NULL() 
-  REAL    , DIMENSION(:), POINTER :: ZYMASS=>NULL() 
   REAL    , DIMENSION(:,:,:), POINTER :: ZZMASS=>NULL() 
   REAL    , DIMENSION(:,:,:), POINTER :: ZPRES_COEF=>NULL() 
   REAL    , DIMENSION(:,:,:), POINTER :: ZSCOORD_SEG=>NULL()
@@ -92,8 +84,6 @@ SUB_ELEC_MODEL(KFROM)%ISCELL_NUMBER=>ISCELL_NUMBER
 SUB_ELEC_MODEL(KFROM)%ISNBSEG=>ISNBSEG            
 SUB_ELEC_MODEL(KFROM)%ISTCOUNT_NUMBER=>ISTCOUNT_NUMBER
 SUB_ELEC_MODEL(KFROM)%ISTYPE=>ISTYPE           
-SUB_ELEC_MODEL(KFROM)%ZXMASS=>ZXMASS
-SUB_ELEC_MODEL(KFROM)%ZYMASS=>ZYMASS
 SUB_ELEC_MODEL(KFROM)%ZZMASS=>ZZMASS
 SUB_ELEC_MODEL(KFROM)%ZPRES_COEF=>ZPRES_COEF
 SUB_ELEC_MODEL(KFROM)%ZSCOORD_SEG=>ZSCOORD_SEG
@@ -109,8 +99,6 @@ ISCELL_NUMBER=>SUB_ELEC_MODEL(KTO)%ISCELL_NUMBER
 ISNBSEG=>SUB_ELEC_MODEL(KTO)%ISNBSEG        
 ISTCOUNT_NUMBER=>SUB_ELEC_MODEL(KTO)%ISTCOUNT_NUMBER 
 ISTYPE=>SUB_ELEC_MODEL(KTO)%ISTYPE              
-ZXMASS=>SUB_ELEC_MODEL(KTO)%ZXMASS 
-ZYMASS=>SUB_ELEC_MODEL(KTO)%ZYMASS 
 ZZMASS=>SUB_ELEC_MODEL(KTO)%ZZMASS 
 ZPRES_COEF=>SUB_ELEC_MODEL(KTO)%ZPRES_COEF 
 ZSCOORD_SEG=>SUB_ELEC_MODEL(KTO)%ZSCOORD_SEG
diff --git a/src/MNH/mode_gridproj.f90 b/src/MNH/mode_gridproj.f90
index d907e68015248c67d0c3c1bcbc229c7a3cf6d1c5..c055e7c064aa36041dd96176ad866f847c904e10 100644
--- a/src/MNH/mode_gridproj.f90
+++ b/src/MNH/mode_gridproj.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -71,11 +71,11 @@ CONTAINS
 !*                1.  ROUTINE  SM_GRIDPROJ   
 !                     --------------------
 !-------------------------------------------------------------------------------
-!      ####################################################################
-       SUBROUTINE SM_GRIDPROJ(PXHAT,PYHAT,PZHAT,PZS,                  &
-                              OSLEVE,PLEN1,PLEN2,PZSMT,PLATOR,PLONOR, &
-                              PMAP,PLAT,PLON,PDXHAT,PDYHAT,PZZ,PJ     )
-!      ####################################################################
+!      ######################################################################
+       SUBROUTINE SM_GRIDPROJ( PXHAT, PYHAT, PZHAT, PXHATM, PYHATM, PZS,    &
+                               OSLEVE, PLEN1, PLEN2, PZSMT, PLATOR, PLONOR, &
+                               PMAP, PLAT, PLON, PDXHAT, PDYHAT, PZZ, PJ    )
+!      ######################################################################
 !
 !!*****  *SM_GRIDPROJ * - Computes Jacobian J, map factor M,
 !!    horizontal grid-meshes, latitude and longitude  at the 
@@ -199,8 +199,8 @@ IMPLICIT NONE
 !
 !*       0.1   Declarations of arguments
 !
-REAL, DIMENSION(:),      INTENT(IN) :: PXHAT,PYHAT,PZHAT ! Positions x,y,z in 
-                                                         ! the cartesian plane
+REAL, DIMENSION(:),      INTENT(IN) :: PXHAT,PYHAT,PZHAT ! Positions x,y,z in the cartesian plane
+REAL, DIMENSION(:),      INTENT(IN) :: PXHATM, PYHATM    ! Positions x,y,z in the cartesian plane at mass points
 REAL, DIMENSION(:,:),    INTENT(IN) :: PZS               ! Orography
 LOGICAL,                INTENT(IN)  :: OSLEVE            ! flag for SLEVE coordinate
 REAL,                   INTENT(IN)  :: PLEN1             ! Decay scale for smooth topography
@@ -332,21 +332,13 @@ END IF
 !*       4.   COMPUTE ZXHAT AND ZYHAT AT MASS POINTS
 !              -------------------------------------
 !
-ZXHATM(:,:) = 0.
-ZYHATM(:,:) = 0.
-ZXHATM(1:IIU-1,1) = .5*(PXHAT(1:IIU-1)+PXHAT(2:IIU))
-ZXHATM(IIU,1)     = 2.*PXHAT(IIU)-ZXHATM(IIU-1,1)
-ZXHATM(:,2:IJU)   = SPREAD(ZXHATM(:,1),2,IJU-1)
+ZXHATM(:,:) = SPREAD( PXHATM(:), 2 , IJU )
+ZYHATM(:,:) = SPREAD( PYHATM(:), 1 , IIU )
 ! cancel MPPDB_CHECK if cprog=='SPAWN  '
-IF(CPROGRAM/='SPAWN ')&
-CALL MPPDB_CHECK2D(ZXHATM,"GRIDPROJ:ZXHATM",PRECISION)
-!
-ZYHATM(1,1:IJU-1) = .5*(PYHAT(1:IJU-1)+PYHAT(2:IJU))
-ZYHATM(1,IJU)     = 2.*PYHAT(IJU)-ZYHATM(1,IJU-1)
-ZYHATM(2:IIU,:)   = SPREAD(ZYHATM(1,:),1,IIU-1)
-! cancel MPPDB_CHECK if cprog=='SPAWN  '
-IF(CPROGRAM/='SPAWN ')&
-CALL MPPDB_CHECK2D(ZYHATM,"GRIDPROJ:ZYHATM",PRECISION)
+IF( CPROGRAM /= 'SPAWN ') THEN
+  CALL MPPDB_CHECK( ZXHATM, "GRIDPROJ:ZXHATM" )
+  CALL MPPDB_CHECK( ZYHATM, "GRIDPROJ:ZXHATM" )
+END IF
 ! ZXHATM and ZXHATM have to be updated
 CALL ADD2DFIELD_ll( TZHALO_ll, ZXHATM, 'SM_GRIDPROJ::ZXHATM' )
 CALL ADD2DFIELD_ll( TZHALO_ll, ZYHATM, 'SM_GRIDPROJ::ZYHATM' )
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index bf87ea0529b690e562e513feed23a747c2633638..8183346f19420ed036dcb5d83a11d227c90b6a8d 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -989,7 +989,7 @@ IF ( nfile_backup_current < NBAK_NUMB ) THEN
     !
     ! Reinitialize Lagragian variables at every model backup
     IF (LLG .AND. LINIT_LG .AND. CINIT_LG=='FMOUT') THEN
-      CALL INI_LG(XXHAT,XYHAT,XZZ,XSVT,XLBXSVM,XLBYSVM)
+      CALL INI_LG( XXHATM, XYHATM, XZZ, XSVT, XLBXSVM, XLBYSVM )
       IF (IVERB>=5) THEN
         WRITE(UNIT=ILUOUT,FMT=*) '************************************'
         WRITE(UNIT=ILUOUT,FMT=*) '*** Lagrangian variables refreshed after ',TRIM(TZBAKFILE%CNAME),' backup'
@@ -1204,8 +1204,9 @@ IF (LCARTESIAN) THEN
   CALL SM_GRIDCART(XXHAT,XYHAT,XZHAT,XZS,LSLEVE,XLEN1,XLEN2,XZSMT,XDXHAT,XDYHAT,XZZ,ZJ)
   XMAP=1.
 ELSE
-  CALL SM_GRIDPROJ(XXHAT,XYHAT,XZHAT,XZS,LSLEVE,XLEN1,XLEN2,XZSMT,XLATORI,XLONORI, &
-                   XMAP,XLAT,XLON,XDXHAT,XDYHAT,XZZ,ZJ)
+  CALL SM_GRIDPROJ( XXHAT, XYHAT, XZHAT, XXHATM, XYHATM, XZS,      &
+                    LSLEVE, XLEN1, XLEN2, XZSMT, XLATORI, XLONORI, &
+                    XMAP, XLAT, XLON, XDXHAT, XDYHAT, XZZ, ZJ      )
 END IF
 !
 IF ( LFORCING ) THEN
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
index be5d47e97df870f165baa47c3b3efb1d37f95961..bb929123903ef945648d36f621ffdec8f1831c5d 100644
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-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.
@@ -781,9 +781,9 @@ CALL SUNPOS_n   ( XZENITH, ZCOSZEN, ZSINZEN, ZAZIMSOL )
 !
   ZTIME1 = ZTIME2
 !
-  CALL SURF_RAD_MODIF (XMAP, XXHAT, XYHAT,                 &
-                  ZCOSZEN, ZSINZEN, ZAZIMSOL, XZS, XZS_XY, &
-                  XDIRFLASWD, XDIRSRFSWD                   )
+  CALL SURF_RAD_MODIF (XMAP, XDXHAT, XDYHAT, XXHATM, XYHATM, &
+                  ZCOSZEN, ZSINZEN, ZAZIMSOL, XZS, XZS_XY,   &
+                  XDIRFLASWD, XDIRSRFSWD                     )
 !
 !* Azimuthal angle to be sent later to surface processes
 !  Defined in radian, clockwise, from North
diff --git a/src/MNH/prep_ideal_case.f90 b/src/MNH/prep_ideal_case.f90
index 5169942d63cb80d8a8dafb21279d94def8a78ee8..d1604d2f123033ad145ea19dd031652f5925a3a0 100644
--- a/src/MNH/prep_ideal_case.f90
+++ b/src/MNH/prep_ideal_case.f90
@@ -1239,7 +1239,8 @@ ELSE
 ! the MESONH horizontal grid is built from the PRE_IDEA1.nam informations
 !------------------------------------------------------------------------
 !
-  ALLOCATE(XXHAT(NIU),XYHAT(NJU))
+  ALLOCATE( XXHAT(NIU),  XYHAT(NJU)  )
+  ALLOCATE( XXHATM(NIU), XYHATM(NJU) )
 !
 ! define the grid localization at the earth surface by the central point
 ! coordinates
@@ -1283,6 +1284,13 @@ ELSE
     XXHAT(:) = (/ (REAL(JLOOP-NIB)*XDELTAX, JLOOP=1,NIU) /)
     XYHAT(:) = (/ (REAL(JLOOP-NJB)*XDELTAY, JLOOP=1,NJU) /)
   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)
 END IF
 !
 !*       5.1.2  Orography and Gal-Chen Sommerville transformation :
@@ -1425,8 +1433,9 @@ IF (LCARTESIAN) THEN
   CALL SM_GRIDCART(XXHAT,XYHAT,XZHAT,XZS,LSLEVE,XLEN1,XLEN2,XZSMT,XDXHAT,XDYHAT,XZZ,XJ)
   XMAP=1.
 ELSE
-  CALL SM_GRIDPROJ(XXHAT,XYHAT,XZHAT,XZS,LSLEVE,XLEN1,XLEN2,XZSMT,XLATORI,XLONORI, &
-                   XMAP,XLAT,XLON,XDXHAT,XDYHAT,XZZ,XJ)
+  CALL SM_GRIDPROJ( XXHAT, XYHAT, XZHAT, XXHATM, XYHATM, XZS,      &
+                    LSLEVE, XLEN1, XLEN2, XZSMT, XLATORI, XLONORI, &
+                    XMAP, XLAT, XLON, XDXHAT, XDYHAT, XZZ, XJ      )
 END IF
 !*       5.4.1  metrics coefficients and update halos:
 !
@@ -1550,10 +1559,10 @@ CALL UPDATE_METRICS(CLBCX,CLBCY,XDXX,XDYY,XDZX,XDZY,XDZZ)
 !
 !*       5.4.2  3D reference state :
 !
-CALL SET_REF(0,TFILE_DUMMY,                        &
-             XZZ,XZHAT,XJ,XDXX,XDYY,CLBCX,CLBCY,   &
-             XREFMASS,XMASS_O_PHI0,XLINMASS,       &
-             XRHODREF,XTHVREF,XRVREF,XEXNREF,XRHODJ)
+CALL SET_REF( 0, TFILE_DUMMY,                            &
+              XZZ, XZHATM, XJ, XDXX, XDYY, CLBCX, CLBCY, &
+              XREFMASS, XMASS_O_PHI0, XLINMASS,          &
+              XRHODREF, XTHVREF, XRVREF, XEXNREF, XRHODJ )
 !
 !
 !*       5.5.1  Absolute pressure :
diff --git a/src/MNH/prep_real_case.f90 b/src/MNH/prep_real_case.f90
index 7014fc3b2d2233902125599097de6e2c5eefcb0f..f73bb09193149fe5afb291745a8abc138900635f 100644
--- a/src/MNH/prep_real_case.f90
+++ b/src/MNH/prep_real_case.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-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.
@@ -858,9 +858,9 @@ IF (LCARTESIAN) THEN
   CALL SM_GRIDCART(XXHAT,XYHAT,XZHAT,XZS,LSLEVE,XLEN1,XLEN2,XZSMT,XDXHAT,XDYHAT,XZZ,ZJ)
   XMAP=1.
 ELSE
-  CALL SM_GRIDPROJ(XXHAT,XYHAT,XZHAT,XZS,                    &
-                   LSLEVE,XLEN1,XLEN2,XZSMT,XLATORI,XLONORI, &
-                   XMAP,XLAT,XLON,XDXHAT,XDYHAT,XZZ,ZJ       )
+  CALL SM_GRIDPROJ( XXHAT, XYHAT, XZHAT, XXHATM, XYHATM, XZS,      &
+                    LSLEVE, XLEN1, XLEN2, XZSMT, XLATORI, XLONORI, &
+                    XMAP, XLAT, XLON, XDXHAT, XDYHAT, XZZ, ZJ      )
 END IF
 !
 CALL MPPDB_CHECK2D(XZS,"prep_real_case8:XZS",PRECISION)
diff --git a/src/MNH/radar_simulator.f90 b/src/MNH/radar_simulator.f90
index b855afc924ed0d8a391fc7d21a9d572158c6fd7f..d03bfb6571bbe080898d4b8396b266d560c68036 100644
--- a/src/MNH/radar_simulator.f90
+++ b/src/MNH/radar_simulator.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2004-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2004-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.
@@ -204,8 +204,6 @@ REAL                                     :: ZCLAT0,ZSLAT0          ! cos and sin
 REAL                                     :: ZMAP     ! Map factor
 REAL                                     :: ZGAMMA,ZCOSG,ZSING  ! angle of projection and its cos and sin values
 !
-REAL, DIMENSION(:),          ALLOCATABLE :: ZXHATM ! X values of the mass points
-REAL, DIMENSION(:),          ALLOCATABLE :: ZYHATM ! Y values of the mass points
 REAL, DIMENSION(:,:,:),      ALLOCATABLE :: ZZM    ! Z values of the mass points
 REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZT_RAY  ! temperature interpolated along the rays
 REAL, DIMENSION(:,:,:,:,:,:),ALLOCATABLE, TARGET :: ZR_RAY  ! rain            mixing ratio interpolated along the rays
@@ -519,16 +517,8 @@ ZSLAT0 = SIN(ZRDSDG*ZLAT0)
 !
 !  Positions of the mass points in the MESO-NH conformal projection 
 !
-ALLOCATE(ZXHATM(IIU))
-ALLOCATE(ZYHATM(IJU)) 
 ALLOCATE(ZZM(IIU,IJU,IKU))
 !
-ZXHATM(1:IIU-1) = .5*(XXHAT(1:IIU-1)+XXHAT(2:IIU))
-ZXHATM(IIU)     = 2.*XXHAT(IIU)-ZXHATM(IIU-1)
-!
-ZYHATM(1:IJU-1) = .5*(XYHAT(1:IJU-1)+XYHAT(2:IJU))
-ZYHATM(IJU)     = 2.*XYHAT(IJU)-ZYHATM(IJU-1)
-!
 ZZM(:,:,1:IKU-1)= .5*(XZZ(:,:,1:IKU-1)+XZZ(:,:,2:IKU))
 ZZM(:,:,IKU)= 2. * XZZ(:,:,IKU) - ZZM(:,:,IKU-1)
 !
@@ -626,7 +616,7 @@ DO JI=1,NBRAD
               ! first compute vertical position (height)          
               ! compute the index of refraction at the radar gate boundaries
               CALL INTERPOL_BEAM(ZN(:,:,:),ZN1,ZX_RAY(JI,JEL,JAZ,JL,JH,JV),&
-              ZY_RAY(JI,JEL,JAZ,JL,JH,JV),ZZ_RAY(JI,JEL,JAZ,JL,JH,JV),ZXHATM(:),ZYHATM(:),ZZM(:,:,:))
+              ZY_RAY(JI,JEL,JAZ,JL,JH,JV),ZZ_RAY(JI,JEL,JAZ,JL,JH,JV),XXHATM(:),XYHATM(:),ZZM(:,:,:))
               IF(LREFR) ZN_RAY(JI,JEL,JAZ,JL,JH,JV)=(ZN1-1.)*1.E6 !LREFR: if true writes out refractivity (N ≡ (n − 1) × 106)
               IF(LDNDZ) THEN !LDNDZ: if true writes out vertical gradient of refractivity
                 IF(JL==1) THEN
@@ -654,7 +644,7 @@ DO JI=1,NBRAD
                   ZDNDZ1=(ZN1-ZN0)/(ZZ_RAY(JI,JEL,JAZ,JL,1,1)-ZZ_RAY(JI,JEL,JAZ,JL-1,1,1))
                 ELSE ! for first gate DNDZ1 is the local value at radar
                   CALL INTERPOL_BEAM(ZDNDZ(:,:,:),ZDNDZ1,ZX_RAY(JI,JEL,JAZ,JL,JH,JV),&
-                            ZY_RAY(JI,JEL,JAZ,JL,JH,JV),ZZ_RAY(JI,JEL,JAZ,JL,JH,JV),ZXHATM(:),ZYHATM(:),ZZM(:,:,:))
+                            ZY_RAY(JI,JEL,JAZ,JL,JH,JV),ZZ_RAY(JI,JEL,JAZ,JL,JH,JV),XXHATM(:),XYHATM(:),ZZM(:,:,:))
                 END IF
                 IF(ZDNDZ1>-ZN1/XRADIUS/COS(ZELEV(JI,JEL,JL,JV))) THEN
                   ZKE=1./(1.+XRADIUS/ZN1*ZDNDZ1*COS(ZELEV(JI,JEL,JL,JV)))
@@ -802,10 +792,10 @@ ENDIF
 !interpolation from TVARMOD to TVARRAD of the model variables in the radar projection, using the position (ZX_RAY, ZY_RAY, ZZ_RAY)
 !of the beam in the model grid 
 CALL INTERPOL_BEAM(TVARMOD,TVARRAD,ZX_RAY(:,:,:,:,:,:),&
-     ZY_RAY(:,:,:,:,:,:),ZZ_RAY(:,:,:,:,:,:),ZXHATM(:),ZYHATM(:),ZZM(:,:,:)) 
+     ZY_RAY(:,:,:,:,:,:),ZZ_RAY(:,:,:,:,:,:),XXHATM(:),XYHATM(:),ZZM(:,:,:))
 !
 DEALLOCATE(ZBU_MASK)
-DEALLOCATE(ZXHATM,ZYHATM,ZZM)
+DEALLOCATE(ZZM)
 DEALLOCATE(ZX_RAY,ZY_RAY)
 DEALLOCATE(TVARMOD,TVARRAD)
 !
diff --git a/src/MNH/read_all_data_grib_case.f90 b/src/MNH/read_all_data_grib_case.f90
index 8a7335846cd898b6ee788269bdc88a11d8ae9d25..6682759c54f679d7c22d338c79281b6c1077fad3 100644
--- a/src/MNH/read_all_data_grib_case.f90
+++ b/src/MNH/read_all_data_grib_case.f90
@@ -370,12 +370,8 @@ ALLOCATE (ZXM(IIU,IJU))
 ALLOCATE (ZYM(IIU,IJU))
 ALLOCATE (ZLONM(IIU,IJU))
 ALLOCATE (ZLATM(IIU,IJU))
-ZXM(1:IIU-1,1) = (XXHAT(1:IIU-1) + XXHAT(2:IIU) ) / 2.
-ZXM(IIU,1)     = XXHAT(IIU) - XXHAT(IIU-1) + ZXM(IIU-1,1)
-ZXM(:,2:IJU)   = SPREAD(ZXM(:,1),2,IJU-1)
-ZYM(1,1:IJU-1) = (XYHAT(1:IJU-1) + XYHAT(2:IJU)) / 2.
-ZYM(1,IJU)     = XYHAT(IJU) - XYHAT(IJU-1) + ZYM(1,IJU-1)
-ZYM(2:IIU,:)   = SPREAD(ZYM(1,:),1,IIU-1)
+ZXM(:,:) = SPREAD(XXHATM(:),2,IJU)
+ZYM(:,:) = SPREAD(XYHATM(:),1,IIU)
 CALL SM_XYTOLATLON_A (XLAT0,XLON0,XRPK,XLATORI,XLONORI,ZXM,ZYM,ZLATM,ZLONM, &
                       IIU,IJU)
 ALLOCATE (ZLONOUT(INO))
diff --git a/src/MNH/read_cams_data_netcdf_case.f90 b/src/MNH/read_cams_data_netcdf_case.f90
index a4023d0917c837950a915e822d7b6ddcdc99c86e..ec6421713a1541e9b1f4f2f32460b345398177e1 100644
--- a/src/MNH/read_cams_data_netcdf_case.f90
+++ b/src/MNH/read_cams_data_netcdf_case.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2012-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2012-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.
@@ -207,12 +207,8 @@ ALLOCATE (ZXM(IIU,IJU))
 ALLOCATE (ZYM(IIU,IJU))
 ALLOCATE (ZLONM(IIU,IJU))
 ALLOCATE (ZLATM(IIU,IJU))
-ZXM(1:IIU-1,1) = (XXHAT(1:IIU-1) + XXHAT(2:IIU) ) / 2.
-ZXM(IIU,1)     = XXHAT(IIU) - XXHAT(IIU-1) + ZXM(IIU-1,1)
-ZXM(:,2:IJU)   = SPREAD(ZXM(:,1),2,IJU-1)
-ZYM(1,1:IJU-1) = (XYHAT(1:IJU-1) + XYHAT(2:IJU)) / 2.
-ZYM(1,IJU)     = XYHAT(IJU) - XYHAT(IJU-1) + ZYM(1,IJU-1)
-ZYM(2:IIU,:)   = SPREAD(ZYM(1,:),1,IIU-1)
+ZXM(:,:) = SPREAD(XXHATM(:),2,IJU)
+ZYM(:,:) = SPREAD(XYHATM(:),1,IIU)
 CALL SM_XYTOLATLON_A (XLAT0,XLON0,XRPK,XLATORI,XLONORI,ZXM,ZYM,ZLATM,ZLONM, &
                       IIU,IJU)
 ALLOCATE (ZLONOUT(INO))
diff --git a/src/MNH/read_chem_data_netcdf_case.f90 b/src/MNH/read_chem_data_netcdf_case.f90
index c1511ef3b4c008d386ba05dee7d3f3542032e7bc..2648709e0c93aed6ecd7ca7c21a4727a94d9881b 100644
--- a/src/MNH/read_chem_data_netcdf_case.f90
+++ b/src/MNH/read_chem_data_netcdf_case.f90
@@ -238,12 +238,8 @@ ALLOCATE (ZXM(IIU,IJU))
 ALLOCATE (ZYM(IIU,IJU))
 ALLOCATE (ZLONM(IIU,IJU))
 ALLOCATE (ZLATM(IIU,IJU))
-ZXM(1:IIU-1,1) = (XXHAT(1:IIU-1) + XXHAT(2:IIU) ) / 2.
-ZXM(IIU,1)     = XXHAT(IIU) - XXHAT(IIU-1) + ZXM(IIU-1,1)
-ZXM(:,2:IJU)   = SPREAD(ZXM(:,1),2,IJU-1)
-ZYM(1,1:IJU-1) = (XYHAT(1:IJU-1) + XYHAT(2:IJU)) / 2.
-ZYM(1,IJU)     = XYHAT(IJU) - XYHAT(IJU-1) + ZYM(1,IJU-1)
-ZYM(2:IIU,:)   = SPREAD(ZYM(1,:),1,IIU-1)
+ZXM(:,:) = SPREAD(XXHATM(:),2,IJU)
+ZYM(:,:) = SPREAD(XYHATM(:),1,IIU)
 CALL SM_XYTOLATLON_A (XLAT0,XLON0,XRPK,XLATORI,XLONORI,ZXM,ZYM,ZLATM,ZLONM, &
                       IIU,IJU)
 ALLOCATE (ZLONOUT(INO))
diff --git a/src/MNH/read_hgridn.f90 b/src/MNH/read_hgridn.f90
index d75ca9bbaf955f56de86e1245297c0b7cddb35ed..b60277ad25c36141ee4ce055ed3f31240dcf38b2 100644
--- a/src/MNH/read_hgridn.f90
+++ b/src/MNH/read_hgridn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1996-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1996-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.
@@ -250,7 +250,17 @@ ENDIF
 
 CALL IO_Field_read(TPFMFILE,'XHAT',XXHAT)
 CALL IO_Field_read(TPFMFILE,'YHAT',XYHAT)
-!
+
+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 )
+
 !JUAN REALZ
 IF ( CPROGRAM .EQ. "REAL  " ) THEN
 IF (.NOT. (ASSOCIATED(XZS))) ALLOCATE(XZS(IIU,IJU))
diff --git a/src/MNH/read_ver_grid.f90 b/src/MNH/read_ver_grid.f90
index 2f8b1fc47b98956bea3fd157a989b8b733546377..bd1de11a99bf5cd3326821565e9a8add17f5c12f 100644
--- a/src/MNH/read_ver_grid.f90
+++ b/src/MNH/read_ver_grid.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -7,7 +7,7 @@
       MODULE MODI_READ_VER_GRID
 !     #########################
 INTERFACE
-      SUBROUTINE READ_VER_GRID(TPPRE_REAL1,PZHAT,OSLEVE,PLEN1,PLEN2)
+      SUBROUTINE READ_VER_GRID( TPPRE_REAL1, PZHAT, OSLEVE, PLEN1, PLEN2 )
 !
 USE MODD_IO, ONLY : TFILEDATA
 !
@@ -21,9 +21,9 @@ END SUBROUTINE READ_VER_GRID
 END INTERFACE
 END MODULE MODI_READ_VER_GRID
 !
-!     ##############################################################
-      SUBROUTINE READ_VER_GRID(TPPRE_REAL1,PZHAT,OSLEVE,PLEN1,PLEN2)
-!     ##############################################################
+!     ####################################################################
+      SUBROUTINE READ_VER_GRID( TPPRE_REAL1, PZHAT, OSLEVE, PLEN1, PLEN2 )
+!     ####################################################################
 !
 !!****  *READ_VER_GRID* - reads namelist data in file PRE_REAL1 for the 
 !!                        initialization of the vertical grid for real cases.
@@ -121,10 +121,10 @@ IMPLICIT NONE
 !*       0.1   Declaration of arguments
 !              ------------------------
 TYPE(TFILEDATA),POINTER,      INTENT(IN) :: TPPRE_REAL1! namelist file
+REAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: PZHAT      ! vertival grid of input fmfile
 LOGICAL,            OPTIONAL, INTENT(IN) :: OSLEVE     ! flag for SLEVE coordinate
 REAL,               OPTIONAL, INTENT(IN) :: PLEN1      ! Decay scale for smooth topography
 REAL,               OPTIONAL, INTENT(IN) :: PLEN2      ! Decay scale for small-scale topography deviation
-REAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: PZHAT      ! vertival grid of input fmfile
 !
 !*       0.2   Declaration of local variables
 !              ------------------------------
@@ -193,7 +193,8 @@ XLEN1_n  = XLEN1
 XLEN2_n  = XLEN2
 !
 IF (CPROGRAM=='REAL  ') THEN
-  IF (ASSOCIATED (XZHAT) ) DEALLOCATE(XZHAT)
+  IF (ASSOCIATED (XZHAT)  ) DEALLOCATE(XZHAT)
+  IF (ASSOCIATED (XZHATM) ) DEALLOCATE(XZHATM)
   CALL POSNAM(IPRE_REAL1,'NAM_BLANKN',GFOUND,ILUOUT0)
   IF (GFOUND) THEN
     CALL INIT_NAM_BLANKn
@@ -214,9 +215,10 @@ SELECT CASE(YZGRID_TYPE)
 CASE('SAMEGR')
   IF (PRESENT(PZHAT) .AND. PRESENT(OSLEVE) .AND. PRESENT(PLEN1) .AND.  PRESENT(PLEN2)) THEN
     IF (NKMAX_n==0)  NKMAX_n=SIZE(PZHAT)-2*JPVEXT
-    ALLOCATE(XZHAT(NKMAX_n+2*JPVEXT))
+    ALLOCATE( XZHAT (IKU) )
+    ALLOCATE( XZHATM(IKU) )
 
-    IF ( (NKMAX_n+2*JPVEXT) > SIZE(PZHAT) ) THEN
+    IF ( (IKU) > SIZE(PZHAT) ) THEN
       WRITE(ILUOUT0,*) 'ERROR IN READ_VER_GRID :'
       WRITE(ILUOUT0,*) '  YOU WANT TO KEEP THE SAME VERTICAL GRID, BUT YOU ASK'
       WRITE(ILUOUT0,*) '  FOR MORE LEVELS THAN IN INPUT FM FILE.'
@@ -226,13 +228,13 @@ CASE('SAMEGR')
       CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_VER_GRID','')
     END IF
 
-    XZHAT(:)=PZHAT(1:NKMAX_n+2*JPVEXT)
+    XZHAT(:)  = PZHAT (1:IKU)
     LTHINSHELL = GTHINSHELL
     LSLEVE_n     = OSLEVE
     XLEN1_n      = PLEN1
     XLEN2_n      = PLEN2
 
-    IF ( (NKMAX_n+2*JPVEXT) == SIZE(PZHAT) ) THEN
+    IF ( (IKU) == SIZE(PZHAT) ) THEN
       WRITE(ILUOUT0,*) 'same vertical grid kept.'
     ELSE
       WRITE(ILUOUT0,*) NKMAX_n,' first levels in vertical grid kept.'
@@ -257,7 +259,8 @@ CASE('SAMEGR')
 !
 CASE('FUNCTN')
 !
-  IF (.NOT. ASSOCIATED(XZHAT)) ALLOCATE(XZHAT(IKU))
+  IF ( .NOT. ASSOCIATED(XZHAT)  ) ALLOCATE( XZHAT (IKU) )
+  IF ( .NOT. ASSOCIATED(XZHATM) ) ALLOCATE( XZHATM(IKU) )
 !
   IF (ABS(ZDZTOP-ZDZGRD) < 1.E-10) THEN
     XZHAT(:) = (/ (REAL(JK-IKB)*ZDZGRD, JK=1,IKU) /)
@@ -297,7 +300,8 @@ CASE('FUNCTN')
 !
 CASE('MANUAL')
 !
-  IF (.NOT. ASSOCIATED(XZHAT)) ALLOCATE(XZHAT(IKU))
+  IF ( .NOT. ASSOCIATED(XZHAT)  ) ALLOCATE( XZHAT (IKU) )
+  IF ( .NOT. ASSOCIATED(XZHATM) ) ALLOCATE( XZHATM(IKU) )
 !
   WRITE(ILUOUT0,FMT=*) 'YZGRID_TYPE="MANUAL", ATTEMPT TO READ VECTOR XZHAT(2,NKU)'
   CALL POSKEY(IPRE_REAL1,ILUOUT0,'ZHAT')
@@ -322,7 +326,11 @@ END SELECT
 !
 !Set model top
 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)
+
 !-------------------------------------------------------------------------------
 !
 !*       5.    TEST ON STRETCHING :
diff --git a/src/MNH/relaxdef.f90 b/src/MNH/relaxdef.f90
index 41665139b06cfe2cdd71187dbf66dba9ac34f17a..216dc038972b2c009df0b4b67e5ee19c2759440d 100644
--- a/src/MNH/relaxdef.f90
+++ b/src/MNH/relaxdef.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -17,7 +17,7 @@ INTERFACE
             OHORELAX_SVCHEM, OHORELAX_SVAER, OHORELAX_SVDST, OHORELAX_SVSLT, &
             OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC,OHORELAX_SVSNW,     &
             PALKTOP,PALKGRD, PALZBOT,PALZBAS,                                &
-            PZZ, PZHAT, PTSTEP,                                              &
+            PZZ, PZHAT, PZHATM, PTSTEP,                                      &
             PRIMKMAX,KRIMX, KRIMY,                                           &
             PALK, PALKW, KALBOT,PALKBAS,PALKWBAS,KALBAS,                     &  
             OMASK_RELAX,PKURELAX, PKVRELAX, PKWRELAX    )
@@ -79,8 +79,9 @@ REAL,                   INTENT(IN)        :: PALZBOT   ! Height of the abs.
 REAL,                   INTENT(IN)        :: PALZBAS   ! Height of the abs.
                                                        ! layer base            
 REAL, DIMENSION(:,:,:), INTENT(IN)        :: PZZ       ! Height 
-REAL, DIMENSION(:),     INTENT(IN)        :: PZHAT     ! Gal-Chen Height 
-REAL,                   INTENT(IN)        :: PTSTEP    ! Time step         
+REAL, DIMENSION(:),     INTENT(IN)        :: PZHAT     ! Gal-Chen Height
+REAL, DIMENSION(:),     INTENT(IN)        :: PZHATM    ! ... at mass points
+REAL,                   INTENT(IN)        :: PTSTEP    ! Time step
 REAL,                     INTENT(IN)    :: PRIMKMAX !Max. value of the horiz.
                                           ! relaxation coefficients
 INTEGER,                INTENT(IN)        :: KRIMX,KRIMY ! Number of points in 
@@ -122,7 +123,7 @@ END MODULE MODI_RELAXDEF
             OHORELAX_SVCHEM, OHORELAX_SVAER, OHORELAX_SVDST, OHORELAX_SVSLT, &
             OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC,OHORELAX_SVSNW,     &
             PALKTOP,PALKGRD, PALZBOT,PALZBAS,                                &
-            PZZ, PZHAT, PTSTEP,                                              &
+            PZZ, PZHAT, PZHATM, PTSTEP,                                      &
             PRIMKMAX,KRIMX, KRIMY,                                           &
             PALK, PALKW, KALBOT,PALKBAS,PALKWBAS,KALBAS,                     &
             OMASK_RELAX,PKURELAX, PKVRELAX, PKWRELAX    )
@@ -306,7 +307,8 @@ REAL,                   INTENT(IN)        :: PALZBAS   ! Height of the abs.
                                                        ! layer base                                                        
 REAL, DIMENSION(:,:,:), INTENT(IN)        :: PZZ       ! Height 
 REAL, DIMENSION(:),     INTENT(IN)        :: PZHAT     ! Gal-Chen Height 
-REAL,                   INTENT(IN)        :: PTSTEP    ! Time step         
+REAL, DIMENSION(:),     INTENT(IN)        :: PZHATM    ! ... at mass points
+REAL,                   INTENT(IN)        :: PTSTEP    ! Time step
 REAL,                     INTENT(IN)    :: PRIMKMAX !Max. value of the horiz.
                                           ! relaxation coefficients
 INTEGER,                INTENT(IN)        :: KRIMX,KRIMY ! Number of points in 
@@ -341,8 +343,6 @@ REAL                       :: ZZSHAT
 REAL                       :: ZZTOP       ! Height of the model top         
 REAL                       :: ZALZHAT     ! Gal-Chen height of the abs. layer 
                                           ! base
-REAL                       :: ZZHATK      ! Gal-Chen height of u(k), v(k), and 
-                                          ! theta(k)
 !
 INTEGER                    :: IKRIMAX     ! Maximum width of the rim zone
                                           !  (number of points)
@@ -445,8 +445,7 @@ IF(OVE_RELAX) THEN
   END DO
 !
   DO JK = KALBOT, IKE
-    ZZHATK = 0.5 * (PZHAT(JK) +PZHAT(JK+1))
-    PALK(JK) = PALKTOP * SIN ( ZWORK * (ZZHATK - PZHAT(KALBOT))) **2
+    PALK(JK) = PALKTOP * SIN ( ZWORK * (PZHATM(JK) - PZHAT(KALBOT))) **2
     PALK(JK) = PALK(JK) / (1. + 2. * PTSTEP * PALK(JK))
   END DO
 !
@@ -469,8 +468,7 @@ IF (OVE_RELAX_GRD) THEN
   END DO
 !
   DO JK = 1,KALBAS
-    ZZHATK = 0.5 * (PZHAT(JK) +PZHAT(JK+1))
-    PALKBAS(JK) = PALKGRD * SIN ( ZWORK * (-ZZHATK + PZHAT(KALBAS))) **2
+    PALKBAS(JK) = PALKGRD * SIN ( ZWORK * (-PZHATM(JK) + PZHAT(KALBAS))) **2
     PALKBAS(JK) = PALKBAS(JK) / (1. + 2. * PTSTEP * PALKBAS(JK))
   END DO
 END IF
diff --git a/src/MNH/series_cloud_elec.f90 b/src/MNH/series_cloud_elec.f90
index cb4d18b427e80960d528689416df35526711dc03..05ced2f2a1366d29bc2d76b3ae5dea0a280cd3ff 100644
--- a/src/MNH/series_cloud_elec.f90
+++ b/src/MNH/series_cloud_elec.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2010-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2010-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.
@@ -94,7 +94,6 @@ USE MODD_CST
 USE MODD_DYN_n,           ONLY: XDXHATM, XDYHATM
 USE MODD_ELEC_DESCR
 USE MODD_ELEC_PARAM
-USE MODD_GRID_n,          ONLY: XXHAT, XYHAT, XZHAT
 USE MODD_IO,              ONLY: TFILEDATA
 USE MODD_NSV,             ONLY: NSV_ELECBEG, NSV_ELECEND
 USE MODD_PARAMETERS
diff --git a/src/MNH/set_advfrc.f90 b/src/MNH/set_advfrc.f90
index 395660c8837f8278d7a8544a33a4f1208464dc9c..eac8ea7b8a8ca0741394ab1d1e84f03f1ee875e4 100644
--- a/src/MNH/set_advfrc.f90
+++ b/src/MNH/set_advfrc.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2013-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2013-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.
@@ -126,7 +126,6 @@ CHARACTER(LEN=48)   :: CFNAM_MEANVAR_ADV
 CHARACTER(LEN=48) :: CFNAM_ADV
 !
 REAL, DIMENSION(:), ALLOCATABLE::     ZHEIGHTMF,ZHEIGHTF,ZTHVUF
-REAL, DIMENSION(:), ALLOCATABLE::     ZZHATM
 REAL, DIMENSION(:), ALLOCATABLE::     ZTHDF,ZRVF,ZPRESS_ADV,ZTHVF
 REAL, DIMENSION(:,:,:), ALLOCATABLE:: ZTHFRC,ZRVFRC
 REAL, DIMENSION(:,:,:), ALLOCATABLE:: ZDRVFRC1D,ZDTHFRC1D,ZDVFRC1D
@@ -190,7 +189,6 @@ ALLOCATE(XDTHFRC(IIU,IJU,IKU,NADVFRC))
 ALLOCATE(XDRVFRC(IIU,IJU,IKU,NADVFRC))
 ! 
 ! For reading in PRE_IDEA1.nam
-ALLOCATE(ZZHATM(IKU))
 ALLOCATE(ZDRVFRC1D(IIU,IKU,NADVFRC))
 ALLOCATE(ZDTHFRC1D(IIU,IKU,NADVFRC))
 ALLOCATE(ZDVFRC1D(IIU,IKU,NADVFRC))
@@ -245,26 +243,24 @@ print*,"  ! 3.2  READ AND INTERPOLATE FORCING"
    !
   END IF
  !
-  ZZHATM(1:IKU-1) = 0.5*(XZHAT(2:IKU)+XZHAT(1:IKU-1))
-  ZZHATM(IKU) = 2.*XZHAT(IKU)-ZZHATM(IKU-1)
 
 print*,"  !! 3.2.2  Vertical interpolation"
   DO JK = 1,IKU
-    IF (ZZHATM(JK) <= ZHEIGHTF(1)) THEN
+    IF (XZHATM(JK) <= ZHEIGHTF(1)) THEN
 !
 print*,"! extrapolation below the first level"
 !   
    
-      ZDZSDH  = (ZZHATM(JK)-ZHEIGHTF(1)) / (ZHEIGHTF(2)-ZHEIGHTF(1))
+      ZDZSDH  = (XZHATM(JK)-ZHEIGHTF(1)) / (ZHEIGHTF(2)-ZHEIGHTF(1))
       ZDRVFRC1D(IIB:IIE,JK,JKT) = ZRVFRC(IIB:IIE,1,JKT) + & 
       (ZRVFRC(IIB:IIE,2,JKT) - ZRVFRC(IIB:IIE,1,JKT)) * ZDZSDH
       ZDTHFRC1D(IIB:IIE,JK,JKT) = ZTHFRC(IIB:IIE,1,JKT) + & 
       (ZTHFRC(IIB:IIE,2,JKT) - ZTHFRC(IIB:IIE,1,JKT)) * ZDZSDH
-    ELSE IF (ZZHATM(JK) > ZHEIGHTF(NPRESSLEV_ADV) ) THEN
+    ELSE IF (XZHATM(JK) > ZHEIGHTF(NPRESSLEV_ADV) ) THEN
 !
 print*,"! extrapolation above the last level"
 !
-      ZDZSDH  = (ZZHATM(JK) - ZHEIGHTF(NPRESSLEV_ADV)) /      &
+      ZDZSDH  = (XZHATM(JK) - ZHEIGHTF(NPRESSLEV_ADV)) /      &
                 (ZHEIGHTF(NPRESSLEV_ADV) - ZHEIGHTF(NPRESSLEV_ADV-1))
       ZDRVFRC1D(IIB:IIE,JK,JKT)  = ZRVFRC(IIB:IIE,NPRESSLEV_ADV,JKT) +  & 
       (ZRVFRC(IIB:IIE,NPRESSLEV_ADV,JKT)-ZRVFRC(IIB:IIE,NPRESSLEV_ADV-1,JKT)) * ZDZSDH
@@ -275,9 +271,9 @@ print*,"! extrapolation above the last level"
 print*,"! interpolation between first and last levels"
 !
       DO JKLEV = 1,NPRESSLEV_ADV-1
-        IF ( (ZZHATM(JK) > ZHEIGHTF(JKLEV)).AND. &
-             (ZZHATM(JK) <= ZHEIGHTF(JKLEV+1))   ) THEN
-          ZDZ1SDH = (ZZHATM(JK) - ZHEIGHTF(JKLEV)) /  &
+        IF ( (XZHATM(JK) > ZHEIGHTF(JKLEV)).AND. &
+             (XZHATM(JK) <= ZHEIGHTF(JKLEV+1))   ) THEN
+          ZDZ1SDH = (XZHATM(JK) - ZHEIGHTF(JKLEV)) /  &
                     (ZHEIGHTF(JKLEV+1)-ZHEIGHTF(JKLEV))
           ZDZ2SDH = 1.- ZDZ1SDH
           ZDRVFRC1D(IIB:IIE,JK,JKT)  = ZRVFRC(IIB:IIE,JKLEV,JKT)*ZDZ2SDH & 
@@ -360,7 +356,6 @@ DEALLOCATE(ZHEIGHTMF)
 DEALLOCATE(ZHEIGHTF)
 DEALLOCATE(ZTHVUF)
 ! pour lecture dans PREIDEA
-DEALLOCATE(ZZHATM)
 DEALLOCATE(ZDRVFRC1D)
 DEALLOCATE(ZDTHFRC1D)
 DEALLOCATE(ZDVFRC1D)
diff --git a/src/MNH/set_frc.f90 b/src/MNH/set_frc.f90
index 6c49fbbf288028ccd6c579aeddf88cc7988c05d8..3f82878d97a8f76f650a7aa1f5c919e2288507e9 100644
--- a/src/MNH/set_frc.f90
+++ b/src/MNH/set_frc.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-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.
@@ -153,8 +153,6 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZHEIGHTMF   ! Height at mass levels
 REAL, DIMENSION(:), ALLOCATABLE :: ZTHVF       ! Thetav at mass levels
 REAL, DIMENSION(:), ALLOCATABLE :: ZTHDF       ! Theta (dry) at mass levels
 REAL, DIMENSION(:), ALLOCATABLE :: ZMRF        ! Vapor mixing ratio at mass lev.
-REAL, DIMENSION(SIZE(XZHAT))    :: ZZHATM      ! Height of mass model grid
-                                               ! levels  without orography
 REAL, DIMENSION(SIZE(XZHAT))    :: ZSHEAR      ! vertical wind shear
 CHARACTER(LEN=4)                :: YZP         ! choice of zfrc or pfrc
 CHARACTER(LEN=100)              :: YMSG
@@ -385,12 +383,9 @@ DO JKT = 1,NFRC
 !
 ! Interpolate and extrapolate Ufrc on u-vertical-grid levels
 !           the other forcing variables.
-!
-  ZZHATM(1:IKU-1) = 0.5*(XZHAT(2:IKU)+XZHAT(1:IKU-1))
-  ZZHATM(IKU) = 2.*XZHAT(IKU)-ZZHATM(IKU-1)
 !
   DO JK = 1,IKU
-    IF (ZZHATM(JK) <= ZHEIGHTF(1)) THEN
+    IF (XZHATM(JK) <= ZHEIGHTF(1)) THEN
 !
 ! copy below the first level
 !
@@ -402,7 +397,7 @@ DO JKT = 1,NFRC
       XTENDRVFRC(JK,JKT) = ZGYRF(1)
       XTENDUFRC(JK,JKT) = ZTUF(1)
       XTENDVFRC(JK,JKT) = ZTVF(1)        
-    ELSE IF (ZZHATM(JK) > ZHEIGHTF(ILEVELF) ) THEN
+    ELSE IF (XZHATM(JK) > ZHEIGHTF(ILEVELF) ) THEN
 !
 ! copy above the last level
 !
@@ -419,9 +414,9 @@ DO JKT = 1,NFRC
 ! interpolation between first and last levels
 !
       DO JKLEV = 1,ILEVELF-1
-        IF ( (ZZHATM(JK) > ZHEIGHTF(JKLEV)).AND. &
-             (ZZHATM(JK) <= ZHEIGHTF(JKLEV+1))   ) THEN
-          ZDZ1SDH = (ZZHATM(JK) - ZHEIGHTF(JKLEV)) /  &
+        IF ( (XZHATM(JK) > ZHEIGHTF(JKLEV)).AND. &
+             (XZHATM(JK) <= ZHEIGHTF(JKLEV+1))   ) THEN
+          ZDZ1SDH = (XZHATM(JK) - ZHEIGHTF(JKLEV)) /  &
                     (ZHEIGHTF(JKLEV+1)-ZHEIGHTF(JKLEV))
           ZDZ2SDH = 1.- ZDZ1SDH
           XUFRC(JK,JKT)  = ZUF(JKLEV)*ZDZ2SDH  + ZUF(JKLEV+1)*ZDZ1SDH
diff --git a/src/MNH/set_geosbal.f90 b/src/MNH/set_geosbal.f90
index 28e528d8b1fbc1d269681cc40769c4f31e2ed273..2f8bd2e68e437c0c2158e26f8ac63d57380c7fca 100644
--- a/src/MNH/set_geosbal.f90
+++ b/src/MNH/set_geosbal.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2010-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2010-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.
@@ -49,7 +49,7 @@ END MODULE MODI_SET_GEOSBAL
 !     fields (at w-grid levels without ororgraphy i.e. at height XZHAT,
 !     with u along latitude circles and v along meridians.)
 !     The vertical profile of mass variable is given at mass-grid levels
-!     (without orography), i.e. at height ZZHATM.
+!     (without orography), i.e. at height XZHATM.
 !       The thermal wind balance is used to compute the potential virtual
 !     temperature from wind field and the vertical profile of thetav.
 !       The vapor mixing ratio is taken  uniform in horizontal plane, i.e.
@@ -289,9 +289,6 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTHV  !  potential virtual temperature
 !
 !*       0.2   Declarations of local variables :
 !
-REAL, DIMENSION(SIZE(XZHAT))    :: ZZHATM  ! Height of mass model grid levels
-                                             ! (without orography)
-
 INTEGER             :: IITRREF  !  iteration number for the determination of
                                  !  the thetav field
 REAL,DIMENSION(:,:),ALLOCATABLE :: ZGAMMA,ZGAMMA_ll ! K(lambda-lambda0) -beta
@@ -311,7 +308,7 @@ REAL, DIMENSION(SIZE(XXHAT),SIZE(XYHAT),SIZE(XZHAT)) :: ZDXX,ZDYY ! metric
                                        ! coefficients dxx,dyy for the grid
                                        ! without orography
 REAL, DIMENSION(SIZE(XZHAT))                         :: ZDZHATM,ZDZHAT  ! deltaz
-                                       ! for ZZHATM and XZHAT levels
+                                       ! for XZHATM and XZHAT levels
 REAL    :: ZRADSDG,ZRVSRD         !Pi/180,Rv/Rd
 INTEGER :: IKB,IKE                 ! useful area in z direction
 INTEGER :: IIU,IJU,IKU             ! Upper bounds in x,y,z directions
@@ -351,14 +348,7 @@ CALL GET_OR_ll('B',IXOR_ll,IYOR_ll)
 !               IN X,Y DIRECTION AND COMPUTE CORIOLIS PARAMETER, AND METRIC
 !               COEFFICIENTS :
 !               -------------------------------------------------------------
-  ZZHATM(1:IKU-1) = 0.5 * (XZHAT(2:IKU)+XZHAT(1:IKU-1))  ! ZHATm(k)=
-                                                     ! 0.5(ZHAT(k+1) +ZHAT(k)
-  ZZHATM(IKU) = 2.* XZHAT(IKU) - ZZHATM(IKU-1)    ! extrapolation for IKU 
-                                                  ! based on deltazhat(iku+1) =
-                                                  ! deltazhat(iku) and  Zhatm(k)
-                                                  ! is the middle point between 
-                                                  ! Zhat(k) and Zhat(k+1)  
-!                                                  
+!
 !*       2.1    Compute a first guess of the dry density
 !
 IF (LTHINSHELL .OR. LCARTESIAN) THEN
@@ -402,7 +392,7 @@ END IF
 !
 ZTHV3D(:,:,:) = SPREAD(SPREAD(PTHVM(:),1,IIU),2,IJU)  ! initialize  with
                                               ! (KILOC,KJLOC) vertical profile
-  ZDZHATM(2:IKU) = ZZHATM(2:IKU)-ZZHATM(1:IKU-1)
+  ZDZHATM(2:IKU) = XZHATM(2:IKU)-XZHATM(1:IKU-1)
   ZDZHAT(1:IKU-1) = XZHAT(2:IKU) - XZHAT(1:IKU-1)
 !
   IF (OBOUSS) THEN
@@ -652,28 +642,28 @@ ZZM(:,:,:)   = MZF(XZZ)                       ! compute height at mass level
                                               ! of grid with orography
 !
 ZZM(:,:,IKU) = 2. * XZZ(:,:,IKU) - ZZM(:,:,IKU-1) ! extrapolate on IKU mass level
-!  ZZM(:,:,1) is always greater than or equal to ZZHATM(1)
-!  ZZM(:,:,IKU) is always smaller than or equal to ZZHATM(IKU)
+!  ZZM(:,:,1) is always greater than or equal to XZHATM(1)
+!  ZZM(:,:,IKU) is always smaller than or equal to XZHATM(IKU)
 !
 DO JI = 1,IIU
   DO JJ = 1,IJU
 !
     DO JK = 1,IKU            ! loop on vertical levels of grid with orography
 !
-      IF (ZZM(JI,JJ,JK) >= ZZHATM(IKU)) THEN     ! copy out when
-        PTHV(JI,JJ,JK)  =  ZTHV3D (JI,JJ,IKU)    ! ZZM(IKU)= ZZHATM(IKU)
+      IF (ZZM(JI,JJ,JK) >= XZHATM(IKU)) THEN     ! copy out when
+        PTHV(JI,JJ,JK)  =  ZTHV3D (JI,JJ,IKU)    ! ZZM(IKU)= XZHATM(IKU)
         XRT(JI,JJ,JK,1) = PMRM(IKU)              ! (in case zs=0.)
 !
-      ELSEIF (ZZM(JI,JJ,JK) < ZZHATM(1)) THEN    ! copy out when  
-        PTHV(JI,JJ,JK)  =  ZTHV3D (JI,JJ,1)      ! ZZM(1)< ZZHATM(1)  
+      ELSEIF (ZZM(JI,JJ,JK) < XZHATM(1)) THEN    ! copy out when
+        PTHV(JI,JJ,JK)  =  ZTHV3D (JI,JJ,1)      ! ZZM(1)< XZHATM(1)
         XRT(JI,JJ,JK,1) = PMRM(1)                ! (in case zs=0.)
 !
       ELSE                  ! search levels on the mass grid without orography
         DO JKS = 2,IKU      ! that surrounded JK
-          IF((ZZM(JI,JJ,JK) >= ZZHATM(JKS-1)).AND.(ZZM(JI,JJ,JK) < ZZHATM(JKS))) &
+          IF((ZZM(JI,JJ,JK) >= XZHATM(JKS-1)).AND.(ZZM(JI,JJ,JK) < XZHATM(JKS))) &
           THEN              ! interpolation with the values on the grid without
                             ! orography
-            ZDZ1SDH = (ZZM(JI,JJ,JK)-ZZHATM(JKS-1))/ (ZZHATM(JKS)-ZZHATM(JKS-1))
+            ZDZ1SDH = (ZZM(JI,JJ,JK)-XZHATM(JKS-1))/ (XZHATM(JKS)-XZHATM(JKS-1))
             ZDZ2SDH = 1. - ZDZ1SDH
             PTHV(JI,JJ,JK)  = (ZDZ1SDH * ZTHV3D(JI,JJ,JKS)   )       &
                             + (ZDZ2SDH * ZTHV3D(JI,JJ,JKS-1) )
diff --git a/src/MNH/set_grid.f90 b/src/MNH/set_grid.f90
index 3945848c8adb12b18354c86efb9cf7dadbd8c13c..dd7daa678bf11bb0a0b0fa413a91b4999a0c8569 100644
--- a/src/MNH/set_grid.f90
+++ b/src/MNH/set_grid.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -9,15 +9,15 @@
 !
 INTERFACE
 !
-      SUBROUTINE SET_GRID(KMI,TPINIFILE,                                      &
-                          KKU,KIMAX_ll,KJMAX_ll,                              &
-                          PTSTEP,PSEGLEN,                                     &
-                          PLONORI,PLATORI,PLON,PLAT,                          &
-                          PXHAT,PYHAT,PDXHAT,PDYHAT, PMAP,                    &
-                          PZS,PZZ,PZHAT,PZTOP,OSLEVE,PLEN1,PLEN2,PZSMT,       &
-                          PJ,                                                 &
-                          TPDTMOD,TPDTCUR,KSTOP,                              &
-                          KBAK_NUMB,KOUT_NUMB,TPBACKUPN,TPOUTPUTN             )
+      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
@@ -47,11 +47,14 @@ 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) :: 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
@@ -80,17 +83,17 @@ END MODULE MODI_SET_GRID
 !
 !
 !
-!     #########################################################################
-      SUBROUTINE SET_GRID(KMI,TPINIFILE,                                      &
-                          KKU,KIMAX_ll,KJMAX_ll,                              &
-                          PTSTEP,PSEGLEN,                                     &
-                          PLONORI,PLATORI,PLON,PLAT,                          &
-                          PXHAT,PYHAT,PDXHAT,PDYHAT, PMAP,                    &
-                          PZS,PZZ,PZHAT,PZTOP,OSLEVE,PLEN1,PLEN2,PZSMT,       &
-                          PJ,                                                 &
-                          TPDTMOD,TPDTCUR,KSTOP,                              &
-                          KBAK_NUMB,KOUT_NUMB,TPBACKUPN,TPOUTPUTN             )
-!     #########################################################################
+!     #####################################################################
+      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     )
+!     #####################################################################
 !
 !!****  *SET_GRID* - routine to set grid variables
 !!
@@ -207,6 +210,7 @@ END MODULE MODI_SET_GRID
 !!     V.MASSON 12/10/00 read of the orography in all cases, even if LFLAT=T
 !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 19/04/2019: removed unused dummy arguments and variables
+!  P. Wautelet 31/08/2022: add PXHATM and PYHATM variables
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -260,11 +264,14 @@ 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) :: 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
@@ -396,11 +403,23 @@ CALL IO_Field_read(TPINIFILE,'DTSEG',TDTSEG)
 !
 !*       2.1    Spatial grid
 !
+
+! 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 )
+
 IF (LCARTESIAN) THEN
   CALL SM_GRIDCART(PXHAT,PYHAT,PZHAT,PZS,OSLEVE,PLEN1,PLEN2,PZSMT,PDXHAT,PDYHAT,PZZ,PJ)
 ELSE
-  CALL SM_GRIDPROJ(PXHAT,PYHAT,PZHAT,PZS,OSLEVE,PLEN1,PLEN2,PZSMT,PLATORI,PLONORI, &
-                   PMAP,PLAT,PLON,PDXHAT,PDYHAT,PZZ,PJ)
+  CALL SM_GRIDPROJ( PXHAT, PYHAT, PZHAT, PXHATM, PYHATM, PZS,      &
+                    OSLEVE, PLEN1, PLEN2, PZSMT, PLATORI, PLONORI, &
+                    PMAP, PLAT, PLON, PDXHAT, PDYHAT, PZZ, PJ      )
 END IF
 !
 !*       2.2    Temporal grid - segment length
diff --git a/src/MNH/set_perturb.f90 b/src/MNH/set_perturb.f90
index 42e384c2a602a5eff3e16415e897edff93ae7922..7e9e436875c30ea4b0dc79ae48fefce303c8079f 100644
--- a/src/MNH/set_perturb.f90
+++ b/src/MNH/set_perturb.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -275,11 +275,10 @@ SELECT CASE(CPERT_KIND)
     DO JK =IKB,IKE
       DO JJ = IJB,IJE
         DO JI = IIB,IIE
-          ZDIST(JI,JJ,JK) = SQRT(                         &
-          (( (XXHAT(JI)+XXHAT(JI+1))*0.5 - ZCENTERX ) / XRADX)**2 + &
-          (( (XYHAT(JJ)+XYHAT(JJ+1))*0.5 - ZCENTERY ) / XRADY)**2 + &
-          (( (XZHAT(JK)+XZHAT(JK+1))*0.5 - XCENTERZ ) / XRADZ)**2   &
-                                )
+          ZDIST(JI,JJ,JK) = SQRT(                                           &
+                                  ( ( XXHATM(JI) - ZCENTERX ) / XRADX)**2 + &
+                                  ( ( XYHATM(JJ) - ZCENTERY ) / XRADY)**2 + &
+                                  ( ( XZHATM(JK) - XCENTERZ ) / XRADZ)**2   )
         END DO
       END DO
     END DO
diff --git a/src/MNH/set_ref.f90 b/src/MNH/set_ref.f90
index 3fbd530b720dad8acead062a8a9854cfbc3950a0..02496e9208ff2a06760bf878f1626b76a2616898 100644
--- a/src/MNH/set_ref.f90
+++ b/src/MNH/set_ref.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -10,7 +10,7 @@ MODULE MODI_SET_REF
 INTERFACE
 !
       SUBROUTINE SET_REF(KMI,TPINIFILE,                                    &
-                         PZZ,PZHAT,PJ,PDXX,PDYY,HLBCX,HLBCY,               &
+                         PZZ,PZHATM,PJ,PDXX,PDYY,HLBCX,HLBCY,              &
                          PREFMASS,PMASS_O_PHI0,PLINMASS,                   &
                          PRHODREF,PTHVREF,PRVREF,PEXNREF,PRHODJ            )
 !
@@ -20,7 +20,7 @@ INTEGER,                INTENT(IN)  :: KMI       ! Model index
 TYPE(TFILEDATA),        INTENT(IN)  :: TPINIFILE ! Initial file
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PZZ       ! Height  of the w levels
                                                  ! with orography                                          
-REAL, DIMENSION(:),     INTENT(IN)  :: PZHAT     ! Height of the w levels
+REAL, DIMENSION(:),     INTENT(IN)  :: PZHATM    ! Height of the w levels at mass points
            ! in the transformed space (GCS transf.) or without orography                                         
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PJ        ! Jacobian 
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDXX      ! metric coefficient dxx
@@ -54,7 +54,7 @@ END MODULE MODI_SET_REF
 !
 !     #########################################################################
       SUBROUTINE SET_REF(KMI,TPINIFILE,                                &
-                         PZZ,PZHAT,PJ,PDXX,PDYY,HLBCX,HLBCY,           &
+                         PZZ,PZHATM,PJ,PDXX,PDYY,HLBCX,HLBCY,          &
                          PREFMASS,PMASS_O_PHI0,PLINMASS,               &
                          PRHODREF,PTHVREF,PRVREF,PEXNREF,PRHODJ        )
 !     #########################################################################
@@ -176,7 +176,7 @@ INTEGER,                INTENT(IN)  :: KMI       ! Model index
 TYPE(TFILEDATA),        INTENT(IN)  :: TPINIFILE ! Initial file
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PZZ       ! Height  of the w levels
                                                  ! with orography                                          
-REAL, DIMENSION(:),     INTENT(IN)  :: PZHAT     ! Height of the w levels
+REAL, DIMENSION(:),     INTENT(IN)  :: PZHATM    ! Height of the w levels at mass points
            ! in the transformed space (GCS transf.) or without orography                                         
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PJ        ! Jacobian 
 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PDXX      ! metric coefficient dxx
@@ -210,8 +210,6 @@ REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZZM
                                                  ! with orography
 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZRHOREF
                                                  ! Reference density
-REAL, DIMENSION(SIZE(PZZ,3))    :: ZZHATM        ! height of the mass levels
-               ! in the transformed space (GCS transf.) or without orography 
 REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)) :: ZDENSOC,ZPFLUX,ZPMASS
 !
 INTEGER             :: IIU        ! Upper dimension in x direction
@@ -282,11 +280,9 @@ END IF
 ! 
 DO JK = 1,IKU-1
   ZZM(:,:,JK) = 0.5*(PZZ(:,:,JK) + PZZ(:,:,JK+1))
-  ZZHATM(JK)  = 0.5*(PZHAT(JK)+PZHAT(JK+1))
 END DO
-ZZHATM(IKU) = 2.* PZHAT(IKU) -ZZHATM(IKU-1)
 ZZM(:,:,IKU)    = 2.* PZZ(:,:,IKU)   -ZZM(:,:,IKU-1)
-! ZZM(:,:,IKU) is always smaller than or equal ZZHATM(IKU)
+! ZZM(:,:,IKU) is always smaller than or equal PZHATM(IKU)
 !
 !
 CALL MPPDB_CHECK3D(ZZM,"SET_REF::ZZM",PRECISION)
@@ -304,16 +300,16 @@ ELSE
 !
     DO JK = 1,IKU
 !
-      IF (ZZM(JI,JJ,JK) >= ZZHATM(IKU)) THEN     ! copy out when  
-        PTHVREF(JI,JJ,JK) =  XTHVREFZ(IKU)       ! ZZM(IKU)= ZZHATM(IKU)  
+      IF (ZZM(JI,JJ,JK) >= PZHATM(IKU)) THEN     ! copy out when
+        PTHVREF(JI,JJ,JK) =  XTHVREFZ(IKU)       ! ZZM(IKU)= PZHATM(IKU)
         PRHODREF(JI,JJ,JK) =  XRHODREFZ(IKU)     ! (in case zs=0.)
 ! 
       ELSE              ! search levels on the mass grid without orography
-        IF (ZZM(JI,JJ,JK) < ZZHATM(2)) THEN
+        IF (ZZM(JI,JJ,JK) < PZHATM(2)) THEN
            IKS=3
         ELSE
           SEARCH : DO JKS = 3,IKU
-            IF((ZZM(JI,JJ,JK) >= ZZHATM(JKS-1)).AND.(ZZM(JI,JJ,JK) < ZZHATM(JKS))) &
+            IF((ZZM(JI,JJ,JK) >= PZHATM(JKS-1)).AND.(ZZM(JI,JJ,JK) < PZHATM(JKS))) &
             THEN          ! interpolation with the values on the grid without
                           ! orography
               IKS=JKS
@@ -321,7 +317,7 @@ ELSE
             END IF
           END DO SEARCH
         END IF
-        ZDZ1SDZ = (ZZM(JI,JJ,JK)-ZZHATM(IKS-1)) / (ZZHATM(IKS)-ZZHATM(IKS-1))
+        ZDZ1SDZ = (ZZM(JI,JJ,JK)-PZHATM(IKS-1)) / (PZHATM(IKS)-PZHATM(IKS-1))
         ZDZ2SDZ = 1. - ZDZ1SDZ
         PTHVREF(JI,JJ,JK) = ( ZDZ1SDZ* XTHVREFZ(IKS) )  &
                           + (ZDZ2SDZ* XTHVREFZ(IKS-1) )      
diff --git a/src/MNH/set_refz.f90 b/src/MNH/set_refz.f90
index f6e82cd85f4e20a2e30627462ed593af03d6d0b8..7822c7e60b1cc345e89cfa995891195401abd656 100644
--- a/src/MNH/set_refz.f90
+++ b/src/MNH/set_refz.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -193,8 +193,7 @@ XTHVREFZ(:)=-999.
 !ocl scalar
 !!!!!!!!!!!!!!!!  FUJI  compiler directive !!!!!!!!!!
 DO JK=IKB,IKU-1
-    XTHVREFZ(JK)= ZSECT(0.5*(XZHAT(JK)+XZHAT(JK+1)),             &
-                        ZZMASS(:,:,IKB:IKE+1),PTHV(:,:,IKB:IKE+1))
+    XTHVREFZ(JK) = ZSECT( XZHATM(JK), ZZMASS(:,:,IKB:IKE+1), PTHV(:,:,IKB:IKE+1) )
 END DO
     XTHVREFZ(IKU)=XTHVREFZ(IKU-1)                                 &
                  +(XTHVREFZ(IKU-1)-XTHVREFZ(IKU-2))               &
@@ -217,8 +216,7 @@ XTHVREFZ(1)=XTHVREFZ(2)
 ZRREFZ(:)=-999.
 !ocl scalar
 DO JK=IKB,IKU-1
-    ZRREFZ(JK)= ZSECT(0.5*(XZHAT(JK)+XZHAT(JK+1)),              &
-                      ZZMASS(:,:,IKB:IKE+1),PRV(:,:,IKB:IKE+1))
+    ZRREFZ(JK) = ZSECT( XZHATM(JK), ZZMASS(:,:,IKB:IKE+1), PRV(:,:,IKB:IKE+1) )
 END DO
     ZRREFZ(IKU)=ZRREFZ(IKU-1)                                   &
                +(ZRREFZ(IKU-1)-ZRREFZ(IKU-2))                   &
@@ -233,7 +231,7 @@ IF (ZRREFZ(IMINLEVEL)==0) THEN
 ELSE
   ZCOEFB=-(LOG(ZRREFZ(IMINLEVEL+1))-LOG(ZRREFZ(IMINLEVEL))) &
        /(0.5*(XZHAT(IMINLEVEL+2)-XZHAT(IMINLEVEL)))
-  ZCOEFA=ZRREFZ(IMINLEVEL)*EXP(ZCOEFB*0.5*(XZHAT(IMINLEVEL+1)+XZHAT(IMINLEVEL)))
+  ZCOEFA=ZRREFZ(IMINLEVEL)*EXP(ZCOEFB*XZHATM(IMINLEVEL))
   WHERE (ZRREFZ==-999.)
     ZRREFZ(:)=ZCOEFA*EXP(-ZCOEFB*0.5*(XZHAT(:)+EOSHIFT(XZHAT(:),1)))
   END WHERE
diff --git a/src/MNH/set_relfrc.f90 b/src/MNH/set_relfrc.f90
index d53c5ae3696c64ed0033d848869479f61bab9f88..857e92ede1b7d77173aff94e05a03c809f730a6e 100644
--- a/src/MNH/set_relfrc.f90
+++ b/src/MNH/set_relfrc.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2013-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2013-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.
@@ -126,7 +126,6 @@ CHARACTER(LEN=48)   :: CFNAM_MEANVAR_REL
 CHARACTER(LEN=28) :: CFNAM_REL
 !
 REAL, DIMENSION(:), ALLOCATABLE::     ZHEIGHTMFR,ZHEIGHTFR,ZTHVUFR,ZTHVUF
-REAL, DIMENSION(:), ALLOCATABLE::     ZZHATM
 REAL, DIMENSION(:), ALLOCATABLE::     ZPRESS_REL
 REAL, DIMENSION(:), ALLOCATABLE::     ZTHDFR,ZTHVFR,ZRVFR
 !
@@ -192,9 +191,6 @@ ALLOCATE(ZHEIGHTFR(NPRESSLEV_REL))
 ! allocations pour le module moddb_advfrc
 !  Adv forcing
 ! 
-! For reading in PRE_IDEA1.nam
-ALLOCATE(ZZHATM(IKU))
-!
 ! relaxation profile
 ALLOCATE(ZRVREL1D(IIU,IKU,NRELFRC))
 ALLOCATE(ZVREL1D(IIU,IKU,NRELFRC))
@@ -259,15 +255,11 @@ DO JKT = 1,NRELFRC
 
   END IF
   !
-  ZZHATM(1:IKU-1) = 0.5*(XZHAT(2:IKU)+XZHAT(1:IKU-1))
-  ZZHATM(IKU) = 2.*XZHAT(IKU)-ZZHATM(IKU-1)
-!
-
   DO JK = 1,IKU
-    IF (ZZHATM(JK) <= ZHEIGHTFR(1)) THEN
+    IF (XZHATM(JK) <= ZHEIGHTFR(1)) THEN
 !
  ! extrapolation below the first level
-      ZDZSDH  = (ZZHATM(JK)-ZHEIGHTFR(1)) / (ZHEIGHTFR(2)-ZHEIGHTFR(1))
+      ZDZSDH  = (XZHATM(JK)-ZHEIGHTFR(1)) / (ZHEIGHTFR(2)-ZHEIGHTFR(1))
       !
       ZRVREL1D(IIB:IIE,JK,JKT) = ZRVREL(IIB:IIE,1,JKT) + & 
       (ZRVREL(IIB:IIE,2,JKT) - ZRVREL(IIB:IIE,1,JKT)) * ZDZSDH
@@ -276,10 +268,10 @@ DO JKT = 1,NRELFRC
       ZTHREL1D(IIB:IIE,JK,JKT) = ZtHREL(IIB:IIE,1,JKT) + & 
       (ZTHREL(IIB:IIE,2,JKT) - ZTHREL(IIB:IIE,1,JKT)) * ZDZSDH
 
-    ELSE IF (ZZHATM(JK) > ZHEIGHTFR(NPRESSLEV_REL) ) THEN
+    ELSE IF (XZHATM(JK) > ZHEIGHTFR(NPRESSLEV_REL) ) THEN
 !
  ! extrapolation above the last level
-      ZDZSDH  = (ZZHATM(JK) - ZHEIGHTFR(NPRESSLEV_REL)) /      &
+      ZDZSDH  = (XZHATM(JK) - ZHEIGHTFR(NPRESSLEV_REL)) /      &
                 (ZHEIGHTFR(NPRESSLEV_REL) - ZHEIGHTFR(NPRESSLEV_REL-1))
       !
       ZRVREL1D(IIB:IIE,JK,JKT)  = ZRVREL(IIB:IIE,NPRESSLEV_REL,JKT) +  & 
@@ -294,9 +286,9 @@ DO JKT = 1,NRELFRC
  ! interpolation between first and last levels
 !
       DO JKLEV = 1,NPRESSLEV_REL-1
-        IF ( (ZZHATM(JK) > ZHEIGHTFR(JKLEV)).AND. &
-             (ZZHATM(JK) <= ZHEIGHTFR(JKLEV+1))   ) THEN
-          ZDZ1SDH = (ZZHATM(JK) - ZHEIGHTFR(JKLEV)) /  &
+        IF ( (XZHATM(JK) > ZHEIGHTFR(JKLEV)).AND. &
+             (XZHATM(JK) <= ZHEIGHTFR(JKLEV+1))   ) THEN
+          ZDZ1SDH = (XZHATM(JK) - ZHEIGHTFR(JKLEV)) /  &
                     (ZHEIGHTFR(JKLEV+1)-ZHEIGHTFR(JKLEV))
           ZDZ2SDH = 1.- ZDZ1SDH
           ZRVREL1D(IIB:IIE,JK,JKT)  = ZRVREL(IIB:IIE,JKLEV,JKT)*ZDZ2SDH & 
@@ -376,8 +368,6 @@ DEALLOCATE(ZTHVUFR)
 DEALLOCATE(ZRVFR)
 
 DEALLOCATE(ZHEIGHTFR)
-! pour lecture dans PREIDEA
-DEALLOCATE(ZZHATM)
 DEALLOCATE(ZPRESS_REL)
 
 DEALLOCATE(ZRVREL1D)
diff --git a/src/MNH/setlb_lg.f90 b/src/MNH/setlb_lg.f90
index 89bb3d3ecd0a556b7194a93f99b29c8540740641..bbdf91a69606292a347431757efb7364bc078968 100644
--- a/src/MNH/setlb_lg.f90
+++ b/src/MNH/setlb_lg.f90
@@ -1,6 +1,6 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2001-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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
 !     ###################
@@ -50,13 +50,13 @@
 !              ------------
 !
 !
-USE MODD_TIME
-USE MODD_TIME_n
-USE MODD_LBC_n
-USE MODD_LSFIELD_n
 USE MODD_GRID_n
+USE MODD_LBC_n
 USE MODD_LG
+USE MODD_LSFIELD_n
 USE MODD_NSV,         ONLY : NSV_LGBEG,NSV_LGEND
+USE MODD_TIME
+USE MODD_TIME_n
 !
 USE MODE_DATETIME
 USE MODE_ll
@@ -80,19 +80,12 @@ IKU=SIZE(XZZ,3)
 CALL DATETIME_DISTANCE(TDTEXP,TDTCUR,ZTEMP_DIST)
 !
 IF ( CLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
-  DO JK=1,IKU
-    DO JJ=1,IJU
-      XLBXSVM(1,JJ,JK,NSV_LGBEG)=0.5*(XXHAT(1)+XXHAT(2))-ZTEMP_DIST*XLGSPEED
-    END DO
-  END DO
-!
-  DO JK=1,IKU
-    DO JJ=1,IJU-1
-      XLBXSVM(1,JJ,JK,NSV_LGBEG+1)=0.5*(XYHAT(JJ)+XYHAT(JJ+1))
-    END DO
-    XLBXSVM(1,IJU,JK,NSV_LGBEG+1)=1.5*XYHAT(IJU)-0.5*XYHAT(IJU-1)
+  XLBXSVM(1,1:IJU,1:IKU,NSV_LGBEG)  = XXHATM(1) - ZTEMP_DIST * XLGSPEED
+
+  DO JK = 1, IKU
+    XLBXSVM(1,1:IJU,JK,NSV_LGBEG+1) = XYHATM(1:IJU)
   END DO
-  !
+
   DO JJ=1,IJU
     DO JK=1,IKU-1
       XLBXSVM(1,JJ,JK,NSV_LGEND)=0.5*(XZZ(1,JJ,JK)+XZZ(1,JJ,JK+1))
@@ -100,69 +93,44 @@ IF ( CLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
     XLBXSVM(1,JJ,IKU,NSV_LGEND)=1.5*XZZ(1,JJ,IKU)  -0.5*XZZ(1,JJ,IKU-1)
   END DO
 END IF
-  !
+
 IF ( CLBCX(1) /= "CYCL" .AND. LEAST_ll()) THEN
-  DO JK=1,IKU
-    DO JJ=1,IJU
-      XLBXSVM(SIZE(XLBXSVM,1),JJ,JK,NSV_LGBEG)=0.5*(XXHAT(IIU-1)+XXHAT(IIU))+ZTEMP_DIST*XLGSPEED
-    END DO
-  END DO
-!
-  DO JK=1,IKU
-    DO JJ=1,IJU-1
-      XLBXSVM(SIZE(XLBXSVM,1),JJ,JK,NSV_LGBEG+1)=0.5*(XYHAT(JJ)+XYHAT(JJ+1))
-    END DO
-    XLBXSVM(SIZE(XLBXSVM,1),IJU,JK,NSV_LGBEG+1)=1.5*XYHAT(IJU)-0.5*XYHAT(IJU-1)
+  XLBXSVM(SIZE(XLBXSVM,1),1:IJU,1:IKU,NSV_LGBEG)  = XXHATM(IIU-1) + ZTEMP_DIST * XLGSPEED
+
+  DO JK = 1, IKU
+    XLBXSVM(SIZE(XLBXSVM,1),1:IJU,JK,NSV_LGBEG+1) = XYHATM(1:IJU)
   END DO
-!
+
   DO JJ=1,IJU
     DO JK=1,IKU-1
       XLBXSVM(SIZE(XLBXSVM,1),JJ,JK,NSV_LGEND)=0.5*(XZZ(IIU,JJ,JK)+XZZ(IIU,JJ,JK+1))
     END DO
     XLBXSVM(SIZE(XLBXSVM,1),JJ,IKU,NSV_LGEND)=1.5*XZZ(IIU,JJ,IKU)-0.5*XZZ(IIU,JJ,IKU-1)
   END DO
-  !
-ENDIF
-  !
+END IF
+
 IF (SIZE(XLBYSVM,1) .NE. 0 .AND. CLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
-!
-  DO JK=1,IKU
-    DO JI=1,IIU-1
-      XLBYSVM(JI,1,JK,NSV_LGBEG)=0.5*(XXHAT(JI)+XXHAT(JI+1))
-    END DO
-    XLBYSVM(IIU,1,JK,NSV_LGBEG)=1.5*XXHAT(IIU)-0.5*XXHAT(IIU-1)
-  END DO
-      !
-  DO JK=1,IKU
-    DO JI=1,IIU
-      XLBYSVM(JI,1,JK,NSV_LGBEG+1)=0.5*(XYHAT(1)+XYHAT(2))-ZTEMP_DIST*XLGSPEED
-    END DO
+  DO JK = 1, IKU
+    XLBYSVM(1:IIU,1,JK,NSV_LGBEG)    = XXHATM(1:IIU)
   END DO
-      !
+
+  XLBYSVM(1:IIU,1,1:IKU,NSV_LGBEG+1) = XYHATM(1) - ZTEMP_DIST * XLGSPEED
+
   DO JI=1,IIU
     DO JK=1,IKU-1
       XLBYSVM(JI,1,JK,NSV_LGEND)=0.5*(XZZ(JI,1,JK)+XZZ(JI,1,JK+1))
     END DO
     XLBYSVM(JI,1,IKU,NSV_LGEND)=1.5*XZZ(JI,1,IKU)  -0.5*XZZ(JI,1,IKU-1)
   END DO
+END IF
 
-ENDIF
-    !
 IF (SIZE(XLBYSVM,1) .NE. 0 .AND. CLBCY(1) /= "CYCL" .AND. LNORTH_ll()) THEN
-!
-  DO JK=1,IKU
-    DO JI=1,IIU-1
-      XLBYSVM(JI,SIZE(XLBYSVM,2),JK,NSV_LGBEG)=0.5*(XXHAT(JI)+XXHAT(JI+1))
-    END DO
-    XLBYSVM(IIU,SIZE(XLBYSVM,2),JK,NSV_LGBEG)=1.5*XXHAT(IIU)-0.5*XXHAT(IIU-1)
-  END DO
-!
-  DO JK=1,IKU
-    DO JI=1,IIU
-      XLBYSVM(JI,SIZE(XLBYSVM,2),JK,NSV_LGBEG+1)=0.5*(XYHAT(IJU-1)+XYHAT(IJU))+ZTEMP_DIST*XLGSPEED
-    END DO
+  DO JK = 1, IKU
+    XLBYSVM(1:IIU,SIZE(XLBYSVM,2),JK,NSV_LGBEG)    = XXHATM(1:IIU)
   END DO
-!
+
+  XLBYSVM(1:IIU,SIZE(XLBYSVM,2),1:IKU,NSV_LGBEG+1) = XYHATM(IJU-1) + ZTEMP_DIST * XLGSPEED
+
   DO JI=1,IIU
     DO JK=1,IKU-1
       XLBYSVM(JI,SIZE(XLBYSVM,2),JK,NSV_LGEND)=0.5*(XZZ(JI,IJU,JK)+XZZ(JI,IJU,JK+1))
diff --git a/src/MNH/shallow_mf.f90 b/src/MNH/shallow_mf.f90
index 2ae315ad50a14bfbde7b9d63aa515abfd3646b0a..cac05f86aef9dfeb494f1feddbc2ba1a607cd4b9 100644
--- a/src/MNH/shallow_mf.f90
+++ b/src/MNH/shallow_mf.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -192,7 +192,6 @@ USE MODI_COMPUTE_FRAC_ICE
 USE MODI_SHUMAN_MF
 !
 USE MODI_COMPUTE_BL89_ML
-USE MODD_GRID_n, ONLY : XDXHAT, XDYHAT
 USE MODD_REF_n, ONLY : XTHVREF
 USE MODE_MSG
 !
diff --git a/src/MNH/spawn_grid2.f90 b/src/MNH/spawn_grid2.f90
index 4ba0d58a36220aa6703065ab8e356a3a0494bfaa..7131b7d5aab21f2ec4ee519f90497a91bce75c25 100644
--- a/src/MNH/spawn_grid2.f90
+++ b/src/MNH/spawn_grid2.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-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.
@@ -9,11 +9,11 @@ MODULE MODI_SPAWN_GRID2
 !
 INTERFACE
 !
-     SUBROUTINE SPAWN_GRID2 (KXOR,KYOR,KXEND,KYEND,KDXRATIO,KDYRATIO,         &
-                             PLONOR,PLATOR,PXHAT,PYHAT,PZHAT,PZTOP,           &
-                             OSLEVE,PLEN1,PLEN2,                              &
-                             PZS,PZSMT,PZS_LS,PZSMT_LS,                       &
-                             TPDTMOD,TPDTCUR                                  )
+     SUBROUTINE SPAWN_GRID2( KXOR, KYOR, KXEND, KYEND, KDXRATIO, KDYRATIO,                &
+                             PLONOR, PLATOR, PXHAT, PYHAT, PZHAT, PXHATM, PYHATM, PZHATM, &
+                             PZTOP, OSLEVE, PLEN1, PLEN2,                                 &
+                             PZS, PZSMT, PZS_LS, PZSMT_LS,                                &
+                             TPDTMOD, TPDTCUR                                             )
 !
 USE MODD_TIME
 !
@@ -23,22 +23,23 @@ INTEGER,   INTENT(IN)  :: KYOR,KYEND ! of the model 2 domain, relative to model
 INTEGER,   INTENT(IN)  :: KDXRATIO   !  x and y-direction Resolution ratio
 INTEGER,   INTENT(IN)  :: KDYRATIO   ! between model 2 and model 1
 !
-REAL,                 INTENT(INOUT) :: PLATOR            ! Latitude of the origine point
-REAL,                 INTENT(INOUT) :: PLONOR            ! Longitude of the origine point
-REAL, DIMENSION(:),   INTENT(INOUT) :: PXHAT,PYHAT,PZHAT ! positions x,y,z in the
+REAL,                 INTENT(OUT) :: PLATOR            ! Latitude of the origine point
+REAL,                 INTENT(OUT) :: PLONOR            ! Longitude of the origine point
+REAL, DIMENSION(:),   INTENT(OUT) :: PXHAT,PYHAT,PZHAT ! positions x,y,z in the
                                      ! conformal plane or on the cartesian plane
+REAL, DIMENSION(:),   INTENT(OUT) :: PXHATM, PYHATM, PZHATM ! positions x,y in the
+                                     ! conformal plane or on the cartesian plane at mass points
 REAL,                 INTENT(OUT)   :: PZTOP             ! model top (m)
 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(INOUT) :: PZS               ! orography
-REAL, DIMENSION(:,:), INTENT(INOUT) :: PZSMT             ! smooth orography
+REAL, DIMENSION(:,:), INTENT(OUT)   :: PZS               ! orography
+REAL, DIMENSION(:,:), INTENT(OUT)   :: PZSMT             ! smooth orography
 REAL, DIMENSION(:,:), INTENT(OUT)   :: PZS_LS            ! interpolated orography
 REAL, DIMENSION(:,:), INTENT(OUT)   :: PZSMT_LS          ! interpolated smooth orography
 !
-!
-TYPE (DATE_TIME),     INTENT(INOUT) :: TPDTMOD  ! Date and Time of MODel beginning
-TYPE (DATE_TIME),     INTENT(INOUT) :: TPDTCUR  ! CURent date and time
+TYPE (DATE_TIME),     INTENT(OUT) :: TPDTMOD  ! Date and Time of MODel beginning
+TYPE (DATE_TIME),     INTENT(OUT) :: TPDTCUR  ! CURent date and time
 !
 END SUBROUTINE SPAWN_GRID2
 !
@@ -47,13 +48,13 @@ END INTERFACE
 END MODULE MODI_SPAWN_GRID2
 !
 !
-!     #########################################################################
-     SUBROUTINE SPAWN_GRID2 (KXOR,KYOR,KXEND,KYEND,KDXRATIO,KDYRATIO,         &
-                             PLONOR,PLATOR,PXHAT,PYHAT,PZHAT,PZTOP,           &
-                             OSLEVE,PLEN1,PLEN2,                              &
-                             PZS,PZSMT,PZS_LS,PZSMT_LS,                       &
-                             TPDTMOD,TPDTCUR                                  )
-!     #########################################################################
+!    ######################################################################################
+     SUBROUTINE SPAWN_GRID2( KXOR, KYOR, KXEND, KYEND, KDXRATIO, KDYRATIO,                &
+                             PLONOR, PLATOR, PXHAT, PYHAT, PZHAT, PXHATM, PYHATM, PZHATM, &
+                             PZTOP, OSLEVE, PLEN1, PLEN2,                                 &
+                             PZS, PZSMT, PZS_LS, PZSMT_LS,                                &
+                             TPDTMOD, TPDTCUR                                             )
+!    ######################################################################################
 !
 !!****  *SPAWN_GRID2 * - subroutine to define spatial and temporal grid.
 !!
@@ -184,10 +185,12 @@ INTEGER,   INTENT(IN)  :: KYOR,KYEND ! of the model 2 domain, relative to model
 INTEGER,   INTENT(IN)  :: KDXRATIO   !  x and y-direction Resolution ratio
 INTEGER,   INTENT(IN)  :: KDYRATIO   ! between model 2 and model 1
 !
-REAL,                 INTENT(INOUT) :: PLATOR            ! Latitude of the origine point
-REAL,                 INTENT(INOUT) :: PLONOR            ! Longitude of the origine point
-REAL, DIMENSION(:),   INTENT(INOUT) :: PXHAT,PYHAT,PZHAT ! positions x,y,z in the
+REAL,                 INTENT(OUT) :: PLATOR            ! Latitude of the origine point
+REAL,                 INTENT(OUT) :: PLONOR            ! Longitude of the origine point
+REAL, DIMENSION(:),   INTENT(OUT) :: PXHAT,PYHAT,PZHAT ! positions x,y,z in the
                                      ! conformal plane or on the cartesian plane
+REAL, DIMENSION(:),   INTENT(OUT) :: PXHATM, PYHATM, PZHATM ! positions x,y in the
+                                     ! conformal plane or on the cartesian plane at mass points
 REAL,                 INTENT(OUT)   :: PZTOP             ! model top (m)
 LOGICAL,              INTENT(OUT)   :: OSLEVE            ! flag for SLEVE coordinate
 REAL,                 INTENT(OUT)   :: PLEN1             ! Decay scale for smooth topography
@@ -197,9 +200,8 @@ REAL, DIMENSION(:,:), INTENT(OUT)   :: PZSMT             ! smooth orography
 REAL, DIMENSION(:,:), INTENT(OUT)   :: PZS_LS            ! interpolated orography
 REAL, DIMENSION(:,:), INTENT(OUT)   :: PZSMT_LS          ! interpolated smooth orography
 !
-!
-TYPE (DATE_TIME),     INTENT(INOUT) :: TPDTMOD  ! Date and Time of MODel beginning
-TYPE (DATE_TIME),     INTENT(INOUT) :: TPDTCUR  ! CURent date and time
+TYPE (DATE_TIME),     INTENT(OUT) :: TPDTMOD  ! Date and Time of MODel beginning
+TYPE (DATE_TIME),     INTENT(OUT) :: TPDTCUR  ! CURent date and time
 !
 !*       0.2    Declarations of local variables for print on FM file
 !
@@ -304,7 +306,8 @@ END IF
 !              --------------------------------------
 !
 PZTOP    = XZTOP1
-PZHAT(:) = XZHAT1(:) 
+PZHAT(:)  = XZHAT1(:)
+PZHATM(:) = XZHATM1(:)
 OSLEVE   = LSLEVE1
 PLEN1    = XLEN11
 PLEN2    = XLEN21
@@ -391,6 +394,11 @@ 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
 !
@@ -449,6 +457,11 @@ PLEN2    = XLEN21
   DEALLOCATE(ZYHAT_2D_F)
   DEALLOCATE(ZYHAT_EXTENDED_C)
   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)
+
 !!$=======
 !!$  IXSIZE1=SIZE(XXHAT1)
 !!$  ALLOCATE(ZXHAT_EXTENDED(IXSIZE1+1))
diff --git a/src/MNH/spawn_model2.f90 b/src/MNH/spawn_model2.f90
index 7d286b6eefcd54cd57468413da56fb828b8fdbba..8837733b60b27b1b23673ddffbd5e506afbe8ce6 100644
--- a/src/MNH/spawn_model2.f90
+++ b/src/MNH/spawn_model2.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-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.
@@ -735,6 +735,7 @@ END IF
 !*       4.4   Grid variables (module MODD_GRID2 and MODD_METRICS2):
 !
 ALLOCATE(XXHAT(IIU),XYHAT(IJU),XZHAT(IKU))
+ALLOCATE(XXHATM(IIU),XYHATM(IJU),XZHATM(IKU))
 ALLOCATE(XZTOP)
 ALLOCATE(XMAP(IIU,IJU))
 ALLOCATE(XLAT(IIU,IJU))
@@ -1058,9 +1059,10 @@ ELSE
         NYEND_TMP = NYEND
 ENDIF
 XZS=0.
-CALL SPAWN_GRID2 (NXOR,NYOR,NXEND,NYEND,NDXRATIO,NDYRATIO,                    &
-                  XLONORI,XLATORI,XXHAT,XYHAT,XZHAT,XZTOP,LSLEVE,XLEN1,XLEN2, &
-                  XZS,XZSMT,ZZS_LS,ZZSMT_LS,TDTMOD,TDTCUR                     )
+CALL SPAWN_GRID2( NXOR, NYOR, NXEND, NYEND, NDXRATIO, NDYRATIO,                  &
+                  XLONORI, XLATORI, XXHAT, XYHAT, XZHAT, XXHATM, XYHATM, XZHATM, &
+                  XZTOP, LSLEVE, XLEN1, XLEN2,                                   &
+                  XZS, XZSMT, ZZS_LS, ZZSMT_LS, TDTMOD, TDTCUR                   )
 !
 CALL MPPDB_CHECK2D(ZZS_LS,"SPAWN_MOD2:ZZS_LS",PRECISION)
 CALL MPPDB_CHECK2D(ZZSMT_LS,"SPAWN_MOD2:ZZSMT_LS",PRECISION)
@@ -1079,10 +1081,12 @@ IF (LCARTESIAN) THEN
   CALL SM_GRIDCART(XXHAT,XYHAT,XZHAT,ZZS_LS,LSLEVE,XLEN1,XLEN2,ZZSMT_LS,XDXHAT,XDYHAT,ZZZ_LS,ZJ)
   CALL SM_GRIDCART(XXHAT,XYHAT,XZHAT,XZS   ,LSLEVE,XLEN1,XLEN2,XZSMT   ,XDXHAT,XDYHAT,XZZ   ,ZJ)
 ELSE
-  CALL SM_GRIDPROJ(XXHAT,XYHAT,XZHAT,ZZS_LS,LSLEVE,XLEN1,XLEN2,ZZSMT_LS,&
-                   XLATORI,XLONORI,XMAP,XLAT,XLON,XDXHAT,XDYHAT,ZZZ_LS,ZJ)
-  CALL SM_GRIDPROJ(XXHAT,XYHAT,XZHAT,XZS   ,LSLEVE,XLEN1,XLEN2,XZSMT   ,&
-                   XLATORI,XLONORI,XMAP,XLAT,XLON,XDXHAT,XDYHAT,XZZ   ,ZJ)
+  CALL SM_GRIDPROJ( XXHAT, XYHAT, XZHAT, XXHATM, XYHATM, ZZS_LS,      &
+                    LSLEVE, XLEN1, XLEN2, ZZSMT_LS, XLATORI, XLONORI, &
+                    XMAP, XLAT, XLON, XDXHAT, XDYHAT, ZZZ_LS, ZJ      )
+  CALL SM_GRIDPROJ( XXHAT, XYHAT, XZHAT, XXHATM, XYHATM, XZS,         &
+                    LSLEVE, XLEN1, XLEN2, XZSMT,    XLATORI, XLONORI, &
+                    XMAP, XLAT, XLON, XDXHAT, XDYHAT, XZZ,    ZJ      )
 END IF
 !
 !*       5.4  Compute the metric coefficients
@@ -1108,10 +1112,10 @@ CALL MPPDB_CHECK3D(XDZY,"spawnmod2-aftrupdate_metrics:XDZY",PRECISION)
 !
 !*       5.5    3D Reference state variables :
 !
-CALL SET_REF(0,TFILE_DUMMY,                        &
-             XZZ,XZHAT,ZJ,XDXX,XDYY,CLBCX,CLBCY,   &
-             XREFMASS,XMASS_O_PHI0,XLINMASS,       &
-             XRHODREF,XTHVREF,XRVREF,XEXNREF,XRHODJ)
+CALL SET_REF( 0, TFILE_DUMMY,                            &
+              XZZ, XZHATM, ZJ, XDXX, XDYY, CLBCX, CLBCY, &
+              XREFMASS, XMASS_O_PHI0, XLINMASS,          &
+              XRHODREF, XTHVREF, XRVREF, XEXNREF, XRHODJ )
 !
 CALL SECOND_MNH(ZTIME2)
 !
diff --git a/src/MNH/spawning.f90 b/src/MNH/spawning.f90
index 15c7d98b76f599f4f2b5329cc09333df5c484bfb..f0ef283104e3c88829e862df850bdd5a11a39c2e 100644
--- a/src/MNH/spawning.f90
+++ b/src/MNH/spawning.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-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.
@@ -280,6 +280,7 @@ USE MODD_PRECIP_n
 XXHAT1 => XXHAT
 XYHAT1 => XYHAT
 XZHAT1 => XZHAT
+XZHATM1 => XZHATM
 XZTOP1 => XZTOP
 XZS1 => XZS
 XZSMT1 => XZSMT
diff --git a/src/MNH/surf_rad_modif.f90 b/src/MNH/surf_rad_modif.f90
index a21ce0dc1b934549ae77da4bf3e3ebd704a21394..1cb4f085de74345d4701e2227a760187a1fe14e4 100644
--- a/src/MNH/surf_rad_modif.f90
+++ b/src/MNH/surf_rad_modif.f90
@@ -1,26 +1,23 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 param 2006/05/18 13:07:25
-!-----------------------------------------------------------------
 !    ##########################
      MODULE MODI_SURF_RAD_MODIF 
 !    ##########################
 !
 INTERFACE 
 !
-      SUBROUTINE SURF_RAD_MODIF ( PMAP, PXHAT, PYHAT,            &
-                  PCOSZEN, PSINZEN, PAZIMSOL,PZS,PZS_XY,         &
-                  PDIRFLASWD, PDIRSRFSWD                         )
+      SUBROUTINE SURF_RAD_MODIF ( PMAP, PDXHAT, PDYHAT, PXHATM, PYHATM,  &
+                                  PCOSZEN, PSINZEN, PAZIMSOL,PZS,PZS_XY, &
+                                  PDIRFLASWD, PDIRSRFSWD                 )
 !
 REAL, DIMENSION(:,:),     INTENT(IN) :: PMAP       ! map factor
-REAL, DIMENSION(:),       INTENT(IN) :: PXHAT      ! X coordinate
-REAL, DIMENSION(:),       INTENT(IN) :: PYHAT      ! Y coordinate
+REAL, DIMENSION(:),       INTENT(IN) :: PDXHAT     ! horizontal stretching in x
+REAL, DIMENSION(:),       INTENT(IN) :: PDYHAT     ! horizontal stretching in y
+REAL, DIMENSION(:),       INTENT(IN) :: PXHATM     ! X coordinates at mass points
+REAL, DIMENSION(:),       INTENT(IN) :: PYHATM     ! Y coordinates at mass points
 REAL, DIMENSION(:,:),     INTENT(IN) :: PCOSZEN    ! COS(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PSINZEN    ! SIN(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PAZIMSOL   ! azimuthal solar angle
@@ -36,11 +33,11 @@ END INTERFACE
 !
 END MODULE MODI_SURF_RAD_MODIF
 !
-!     ###################################################################
-      SUBROUTINE SURF_RAD_MODIF ( PMAP, PXHAT, PYHAT,            &
-                  PCOSZEN, PSINZEN, PAZIMSOL,PZS,PZS_XY,         &
-                  PDIRFLASWD, PDIRSRFSWD                         )
-!     ###################################################################
+!     ####################################################################
+      SUBROUTINE SURF_RAD_MODIF ( PMAP, PDXHAT, PDYHAT, PXHATM, PYHATM,  &
+                                  PCOSZEN, PSINZEN, PAZIMSOL,PZS,PZS_XY, &
+                                  PDIRFLASWD, PDIRSRFSWD                 )
+!     ####################################################################
 !
 !!****  * SURF_RAD_MODIF * - computes the modifications to the downwards
 !!                           radiative fluxes at the surface, due to
@@ -113,8 +110,10 @@ IMPLICIT NONE
 !
 !
 REAL, DIMENSION(:,:),     INTENT(IN) :: PMAP       ! map factor
-REAL, DIMENSION(:),       INTENT(IN) :: PXHAT      ! X coordinate
-REAL, DIMENSION(:),       INTENT(IN) :: PYHAT      ! Y coordinate
+REAL, DIMENSION(:),       INTENT(IN) :: PDXHAT     ! horizontal stretching in x
+REAL, DIMENSION(:),       INTENT(IN) :: PDYHAT     ! horizontal stretching in y
+REAL, DIMENSION(:),       INTENT(IN) :: PXHATM     ! X coordinates at mass points
+REAL, DIMENSION(:),       INTENT(IN) :: PYHATM     ! Y coordinates at mass points
 REAL, DIMENSION(:,:),     INTENT(IN) :: PCOSZEN    ! COS(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PSINZEN    ! SIN(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PAZIMSOL   ! azimuthal solar angle
@@ -158,14 +157,14 @@ ISWB = SIZE(PDIRFLASWD,3)
 !-------------------------------------------------------------------------------
 !
 DO JSWB = 1, ISWB
-  CALL SURF_SOLAR_SUM     (PXHAT, PYHAT, ZMAP, PDIRFLASWD(:,:,JSWB), ZENERGY1(JSWB) )
+  CALL SURF_SOLAR_SUM     (PDXHAT, PDYHAT, ZMAP, PDIRFLASWD(:,:,JSWB), ZENERGY1(JSWB) )
 END DO
 !
 !
 !*       2.    Slope direction direct SW effects
 !              ---------------------------------
 !
-CALL SURF_SOLAR_SLOPES  (ZMAP, PXHAT, PYHAT,                 &
+CALL SURF_SOLAR_SLOPES  (ZMAP, PDXHAT, PDYHAT,               &
                          PCOSZEN, PSINZEN, PAZIMSOL,         &
                          PZS, PZS_XY, PDIRFLASWD, ZDIRSWDT   )
 
@@ -173,7 +172,7 @@ CALL SURF_SOLAR_SLOPES  (ZMAP, PXHAT, PYHAT,                 &
 !*       3.    RESOLVED shadows for direct solar radiation
 !              -------------------------------------------
 !
-CALL SURF_SOLAR_SHADOWS (ZMAP, PXHAT, PYHAT,           &
+CALL SURF_SOLAR_SHADOWS (ZMAP, PXHATM, PYHATM,         &
                          PCOSZEN, PSINZEN, PAZIMSOL,   &
                          PZS, PZS_XY, ZDIRSWDT, ZDIRSWD)
 !
@@ -182,11 +181,11 @@ CALL SURF_SOLAR_SHADOWS (ZMAP, PXHAT, PYHAT,           &
 !              -------------------
 !
 DO JSWB = 1, ISWB
-  CALL SURF_SOLAR_SUM(PXHAT, PYHAT, ZMAP,  &
+  CALL SURF_SOLAR_SUM(PDXHAT, PDYHAT, ZMAP,  &
                       ZDIRSWD(:,:,JSWB),   &
                       ZENERGY2(JSWB)       )
   !
-  CALL SURF_SOLAR_SUM(PXHAT, PYHAT, ZMAP,                             &
+  CALL SURF_SOLAR_SUM(PDXHAT, PDYHAT, ZMAP,                             &
                       MAX(ZDIRSWD(:,:,JSWB)-PDIRFLASWD(:,:,JSWB),0.), &
                       ZENERGYP(JSWB)                                  )
   !
diff --git a/src/MNH/surf_solar_shadows.f90 b/src/MNH/surf_solar_shadows.f90
index b751989594ce9b2c0975f956fe7d0402802d7114..0a3eaabb77732cf8f574f43c7c9d3fbf8e8c3314 100644
--- a/src/MNH/surf_solar_shadows.f90
+++ b/src/MNH/surf_solar_shadows.f90
@@ -1,27 +1,22 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 param 2006/05/18 13:07:25
-!-----------------------------------------------------------------
 !    ##############################
      MODULE MODI_SURF_SOLAR_SHADOWS 
 !    ##############################
 !
 INTERFACE 
 !
-      SUBROUTINE SURF_SOLAR_SHADOWS ( PMAP, PXHAT, PYHAT,                     &
+      SUBROUTINE SURF_SOLAR_SHADOWS ( PMAP, PXHATM, PYHATM,                   &
                   PCOSZEN, PSINZEN, PAZIMSOL,                                 &
                   PZS, PZS_XY, PDIRSWDT, PDIRSRFSWD                           )
 !
 !
 REAL, DIMENSION(:,:),     INTENT(IN) :: PMAP         ! map factor
-REAL, DIMENSION(:),       INTENT(IN) :: PXHAT        ! X coordinate
-REAL, DIMENSION(:),       INTENT(IN) :: PYHAT        ! Y coordinate
+REAL, DIMENSION(:),       INTENT(IN) :: PXHATM       ! X coordinate at mass points
+REAL, DIMENSION(:),       INTENT(IN) :: PYHATM       ! Y coordinate at mass points
 REAL, DIMENSION(:,:),     INTENT(IN) :: PCOSZEN ! COS(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PSINZEN ! SIN(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PAZIMSOL! azimuthal solar angle
@@ -39,7 +34,7 @@ END INTERFACE
 !
 END MODULE MODI_SURF_SOLAR_SHADOWS
 !     #########################################################################
-      SUBROUTINE SURF_SOLAR_SHADOWS ( PMAP, PXHAT, PYHAT,                     &
+      SUBROUTINE SURF_SOLAR_SHADOWS ( PMAP, PXHATM, PYHATM,                   &
                   PCOSZEN, PSINZEN, PAZIMSOL,                                 &
                   PZS, PZS_XY, PDIRSWDT, PDIRSRFSWD                           )
 !     #########################################################################
@@ -95,8 +90,8 @@ IMPLICIT NONE
 !
 !
 REAL, DIMENSION(:,:),     INTENT(IN) :: PMAP         ! map factor
-REAL, DIMENSION(:),       INTENT(IN) :: PXHAT        ! X coordinate
-REAL, DIMENSION(:),       INTENT(IN) :: PYHAT        ! Y coordinate
+REAL, DIMENSION(:),       INTENT(IN) :: PXHATM       ! X coordinate at mass points
+REAL, DIMENSION(:),       INTENT(IN) :: PYHATM       ! Y coordinate at mass points
 REAL, DIMENSION(:,:),     INTENT(IN) :: PCOSZEN ! COS(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PSINZEN ! SIN(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PAZIMSOL! azimuthal solar angle
@@ -231,20 +226,20 @@ DO JJ=IJB,IJE
 !
       SELECT CASE (JT)
        CASE (1)
-         ZX=(5.*PXHAT(JI)+PXHAT(JI+1))/6.
-         ZY=0.5*(PYHAT(JJ)+PYHAT(JJ+1))
+         ZX=PXHATM(JI)/6.
+         ZY=PYHATM(JJ)
          ZZ=(PZS(JI,JJ)+PZS_XY(JI,JJ)+PZS_XY(JI,JJ+1))/3.
        CASE (2)
-         ZX=0.5*(PXHAT(JI)+PXHAT(JI+1))
-         ZY=(5.*PYHAT(JJ+1)+PYHAT(JJ))/6.
+         ZX=PXHATM(JI)
+         ZY=PYHATM(JJ)/6.
          ZZ=(PZS(JI,JJ)+PZS_XY(JI,JJ+1)+PZS_XY(JI+1,JJ+1))/3.
        CASE (3)
-         ZX=(5.*PXHAT(JI+1)+PXHAT(JI))/6.
-         ZY=0.5*(PYHAT(JJ)+PYHAT(JJ+1))
+         ZX=PXHATM(JI)/6.
+         ZY=PYHATM(JJ)
          ZZ=(PZS(JI,JJ)+PZS_XY(JI+1,JJ)+PZS_XY(JI+1,JJ+1))/3.
        CASE (4)
-         ZX=0.5*(PXHAT(JI)+PXHAT(JI+1))
-         ZY=(5.*PYHAT(JJ)+PYHAT(JJ+1))/6.
+         ZX=PXHATM(JI)
+         ZY=PYHATM(JJ)/6.
          ZZ=(PZS(JI,JJ)+PZS_XY(JI,JJ)+PZS_XY(JI+1,JJ))/3.
       END SELECT
 !
diff --git a/src/MNH/surf_solar_slopes.f90 b/src/MNH/surf_solar_slopes.f90
index e7ea4ef2284690164c6c4c0679c8621f245ec5dd..0ffcf81e32e18b4cceb8cc16350224eed428d22d 100644
--- a/src/MNH/surf_solar_slopes.f90
+++ b/src/MNH/surf_solar_slopes.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2002-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-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.
@@ -9,14 +9,14 @@
 !
 INTERFACE 
 !
-      SUBROUTINE SURF_SOLAR_SLOPES ( PMAP, PXHAT, PYHAT,                      &
+      SUBROUTINE SURF_SOLAR_SLOPES ( PMAP, PDXHAT, PDYHAT,                    &
                   PCOSZEN, PSINZEN, PAZIMSOL,                                 &
                   PZS, PZS_XY, PDIRSRFSWD, PDIRSWDT                           )
 !
 !
 REAL, DIMENSION(:,:),     INTENT(IN) :: PMAP         ! map factor
-REAL, DIMENSION(:),       INTENT(IN) :: PXHAT        ! X coordinate
-REAL, DIMENSION(:),       INTENT(IN) :: PYHAT        ! Y coordinate
+REAL, DIMENSION(:),       INTENT(IN) :: PDXHAT       ! horizontal stretching in x
+REAL, DIMENSION(:),       INTENT(IN) :: PDYHAT       ! horizontal stretching in y
 REAL, DIMENSION(:,:),     INTENT(IN) :: PCOSZEN ! COS(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PSINZEN ! SIN(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PAZIMSOL! azimuthal solar angle
@@ -34,7 +34,7 @@ END INTERFACE
 !
 END MODULE MODI_SURF_SOLAR_SLOPES
 !     #########################################################################
-      SUBROUTINE SURF_SOLAR_SLOPES ( PMAP, PXHAT, PYHAT,                      &
+      SUBROUTINE SURF_SOLAR_SLOPES ( PMAP, PDXHAT, PDYHAT,                    &
                   PCOSZEN, PSINZEN, PAZIMSOL,                                 &
                   PZS, PZS_XY, PDIRSRFSWD, PDIRSWDT                           )
 !     #########################################################################
@@ -86,8 +86,8 @@ IMPLICIT NONE
 !
 !
 REAL, DIMENSION(:,:),     INTENT(IN) :: PMAP         ! map factor
-REAL, DIMENSION(:),       INTENT(IN) :: PXHAT        ! X coordinate
-REAL, DIMENSION(:),       INTENT(IN) :: PYHAT        ! Y coordinate
+REAL, DIMENSION(:),       INTENT(IN) :: PDXHAT       ! horizontal stretching in x
+REAL, DIMENSION(:),       INTENT(IN) :: PDYHAT       ! horizontal stretching in y
 REAL, DIMENSION(:,:),     INTENT(IN) :: PCOSZEN ! COS(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PSINZEN ! SIN(zenithal solar angle)
 REAL, DIMENSION(:,:),     INTENT(IN) :: PAZIMSOL! azimuthal solar angle
@@ -144,27 +144,27 @@ DO JT=1,4
         CASE (1)
           ZDZSDX=(    2.* PZS   (JI,JJ)                   &
                    - (PZS_XY(JI,JJ)+PZS_XY(JI,JJ+1)) )    &
-                 / (PXHAT(JI+1)-PXHAT(JI))  * PMAP(JI,JJ)
+                 / PDXHAT(JI)  * PMAP(JI,JJ)
           ZDZSDY=(  PZS_XY(JI,JJ+1) - PZS_XY(JI,JJ) )     &
-                 / (PYHAT(JJ+1)-PYHAT(JJ))  * PMAP(JI,JJ)
+                 / PDYHAT(JJ)  * PMAP(JI,JJ)
         CASE (2)
            ZDZSDX=(  PZS_XY(JI+1,JJ+1) -PZS_XY(JI,JJ+1))  &
-                 / (PXHAT(JI+1)-PXHAT(JI))  * PMAP(JI,JJ)
+                 / PDXHAT(JI)  * PMAP(JI,JJ)
            ZDZSDY=(  (PZS_XY(JI+1,JJ+1)+PZS_XY(JI,JJ+1))  &
                      - 2.* PZS (JI,JJ) )                  &
-                 / (PYHAT(JJ+1)-PYHAT(JJ))  * PMAP(JI,JJ)
+                 / PDYHAT(JJ)  * PMAP(JI,JJ)
         CASE (3)
           ZDZSDX=(  (PZS_XY(JI+1,JJ)+PZS_XY(JI+1,JJ+1))   &
                    - 2.* PZS(JI,JJ)                    )  &
-                 / (PXHAT(JI+1)-PXHAT(JI))  * PMAP(JI,JJ)
+                 / PDXHAT(JI)  * PMAP(JI,JJ)
           ZDZSDY=(  PZS_XY(JI+1,JJ+1) - PZS_XY(JI+1,JJ) ) &
-                 / (PYHAT(JJ+1)-PYHAT(JJ))  * PMAP(JI,JJ)
+                 / PDYHAT(JJ)  * PMAP(JI,JJ)
         CASE (4)
            ZDZSDX=(  PZS_XY(JI+1,JJ) - PZS_XY(JI,JJ) )    &
-                 / (PXHAT(JI+1)-PXHAT(JI))  * PMAP(JI,JJ)
+                 / PDXHAT(JI)  * PMAP(JI,JJ)
            ZDZSDY=(  2.* PZS(JI,JJ)                       &
                    - (PZS_XY(JI+1,JJ)+PZS_XY(JI,JJ)) )    &
-                 / (PYHAT(JJ+1)-PYHAT(JJ))  * PMAP(JI,JJ)
+                 / PDYHAT(JJ)  * PMAP(JI,JJ)
       END SELECT
 !
 !* slope angles
diff --git a/src/MNH/surf_solar_sum.f90 b/src/MNH/surf_solar_sum.f90
index 3742792d5b0f261dfc6af9f88d7f71ee10c3de0a..11295f83a009955f5752b047d49f531b47609f10 100644
--- a/src/MNH/surf_solar_sum.f90
+++ b/src/MNH/surf_solar_sum.f90
@@ -1,24 +1,19 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2002-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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
 !-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source$ $Revision$
-! MASDEV4_7 param 2006/05/18 13:07:25
-!-----------------------------------------------------------------
 !    ##########################
      MODULE MODI_SURF_SOLAR_SUM 
 !    ##########################
 !
 INTERFACE 
 !
-      SUBROUTINE SURF_SOLAR_SUM ( PXHAT, PYHAT, PMAP, PDIRSWD, PENERGY )
+      SUBROUTINE SURF_SOLAR_SUM ( PDXHAT, PDYHAT, PMAP, PDIRSWD, PENERGY )
 !
 !
-REAL, DIMENSION(:),   INTENT(IN) :: PXHAT    ! X coordinate
-REAL, DIMENSION(:),   INTENT(IN) :: PYHAT    ! Y coordinate
+REAL, DIMENSION(:),   INTENT(IN) :: PDXHAT   ! horizontal stretching in x
+REAL, DIMENSION(:),   INTENT(IN) :: PDYHAT   ! horizontal stretching in y
 REAL, DIMENSION(:,:), INTENT(IN) :: PMAP     ! map factor
 REAL, DIMENSION(:,:), INTENT(IN) :: PDIRSWD  ! direct SW flux on hor. surf.
 REAL,                 INTENT(OUT):: PENERGY  ! energy received by the surface
@@ -30,9 +25,9 @@ END INTERFACE
 !
 END MODULE MODI_SURF_SOLAR_SUM
 !
-!     ##################################################################
-      SUBROUTINE SURF_SOLAR_SUM ( PXHAT, PYHAT, PMAP, PDIRSWD, PENERGY )
-!     ##################################################################
+!     ####################################################################
+      SUBROUTINE SURF_SOLAR_SUM ( PDXHAT, PDYHAT, PMAP, PDIRSWD, PENERGY )
+!     ####################################################################
 !
 !!****  * SURF_SOLAR_SUM * - computes the sum of energy received by
 !!                           the surface from direct solar radiation
@@ -78,8 +73,8 @@ IMPLICIT NONE
 !*       0.1   DECLARATIONS OF DUMMY ARGUMENTS :
 !
 !
-REAL, DIMENSION(:),   INTENT(IN) :: PXHAT    ! X coordinate
-REAL, DIMENSION(:),   INTENT(IN) :: PYHAT    ! Y coordinate
+REAL, DIMENSION(:),   INTENT(IN) :: PDXHAT   ! horizontal stretching in x
+REAL, DIMENSION(:),   INTENT(IN) :: PDYHAT   ! horizontal stretching in y
 REAL, DIMENSION(:,:), INTENT(IN) :: PMAP     ! map factor
 REAL, DIMENSION(:,:), INTENT(IN) :: PDIRSWD  ! direct SW flux on hor. surf.
 REAL,                 INTENT(OUT):: PENERGY  ! energy received by the surface
@@ -109,9 +104,7 @@ ALLOCATE(ZENERGY_2D(IIB:IIE,IJB:IJE))
 !
 DO JJ=IJB,IJE
   DO JI=IIB,IIE
-    ZENERGY_2D(JI,JJ) = PDIRSWD(JI,JJ)*(PXHAT(JI+1)-PXHAT(JI)) &
-                                      *(PYHAT(JJ+1)-PYHAT(JJ)) &
-                                      /PMAP(JI,JJ)**2
+    ZENERGY_2D(JI,JJ) = PDIRSWD(JI,JJ) * PDXHAT(JI) * PDYHAT(JJ) / PMAP(JI,JJ)**2
   END DO
 END DO
 !
diff --git a/src/MNH/turb_ver_thermo_flux.f90 b/src/MNH/turb_ver_thermo_flux.f90
index aa53f08222efb391a810c39dba484166dfeb9c61..596d78c009c37cdd8f5e6918e56dbca66f61ba2b 100644
--- a/src/MNH/turb_ver_thermo_flux.f90
+++ b/src/MNH/turb_ver_thermo_flux.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -339,7 +339,7 @@ END MODULE MODI_TURB_VER_THERMO_FLUX
 USE MODD_CST
 USE MODD_CTURB
 use modd_field,          only: tfieldmetadata, TYPEREAL
-USE MODD_GRID_n,         ONLY: XZS, XXHAT, XYHAT
+USE MODD_GRID_n,         ONLY: XZS, XXHAT, XXHATM, XYHAT, XYHATM
 USE MODD_IO,             ONLY: TFILEDATA
 USE MODD_METRICS_n,      ONLY: XDXX, XDYY, XDZX, XDZY, XDZZ
 USE MODD_PARAMETERS
@@ -545,10 +545,9 @@ IF (LOCEAN .AND. LDEEPOC) THEN
   CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
   DO JJ = IJB,IJE
     DO JI = IIB,IIE
-      ZDIST(JI,JJ) = SQRT(                         &
-      (( (XXHAT(JI)+XXHAT(JI+1))*0.5 - XCENTX_OC ) / XRADX_OC)**2 + &
-      (( (XYHAT(JJ)+XYHAT(JJ+1))*0.5 - XCENTY_OC ) / XRADY_OC)**2   &
-                                )
+      ZDIST(JI,JJ) = SQRT(                          &
+      ( ( XXHATM(JI) - XCENTX_OC ) / XRADX_OC )**2 + &
+      ( ( XYHATM(JJ) - XCENTY_OC ) / XRADY_OC )**2 )
     END DO
   END DO
   DO JJ=IJB,IJE
diff --git a/src/MNH/ver_interp_to_mixed_grid.f90 b/src/MNH/ver_interp_to_mixed_grid.f90
index 94b161a5d4989fd24bcc8dc2305436acfc304d03..1c1580cc9588bf11b59feda0d5f22d5fc4b2d63f 100644
--- a/src/MNH/ver_interp_to_mixed_grid.f90
+++ b/src/MNH/ver_interp_to_mixed_grid.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -274,7 +274,7 @@ ILU=SIZE(PZMASS_LS,3)
 !*       1.1   Grid definition
 !              ---------------
 !
-IF (MINVAL(PZMASS_LS (:,:,ILU))<0.5*(XZHAT(IKE)+XZHAT(IKE+1))) THEN
+IF (MINVAL(PZMASS_LS (:,:,ILU))<XZHATM(IKE)) THEN
   WRITE(ILUOUT0,*)
   WRITE(ILUOUT0,*) '+-----------------------------------------------------+'
   WRITE(ILUOUT0,*) '| MESONH highest mass level above highest input level |'
diff --git a/src/MNH/ver_thermo.f90 b/src/MNH/ver_thermo.f90
index a4e8ee654b3f55e81978ab3442703e3816afd7ce..186f230a42bb33b3b038d73bb0d58f9b8347c152 100644
--- a/src/MNH/ver_thermo.f90
+++ b/src/MNH/ver_thermo.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
@@ -332,9 +332,9 @@ ALLOCATE(XRVREF(IIU,IJU,IKU))
 ALLOCATE(XEXNREF(IIU,IJU,IKU))
 ALLOCATE(XRHODJ(IIU,IJU,IKU))
 XRVREF(:,:,:) = 0.
-CALL SET_REF(0,TFILE_DUMMY,XZZ,XZHAT,PJ,PDXX,PDYY,CLBCX,CLBCY,       &
-             XREFMASS,XMASS_O_PHI0,XLINMASS,XRHODREF,XTHVREF,XRVREF, &
-             XEXNREF,XRHODJ)
+CALL SET_REF( 0, TFILE_DUMMY, XZZ, XZHATM, PJ, PDXX, PDYY, CLBCX, CLBCY,   &
+              XREFMASS, XMASS_O_PHI0, XLINMASS, XRHODREF, XTHVREF, XRVREF, &
+              XEXNREF, XRHODJ                                              )
 
 CALL MPPDB_CHECK3D(XRHODREF,"VERTHERMO::XRHODREF",PRECISION)
 CALL MPPDB_CHECK3D(XTHVREF,"VERTHERMO::XTHVREF",PRECISION)
diff --git a/src/MNH/write_lfifm1_for_diag.f90 b/src/MNH/write_lfifm1_for_diag.f90
index 0bb1d984b2ca0d5c8a88ff9c8f53b2255dafbb69..d452fecafad20d6692b2651ff661dc2f9ccafe52 100644
--- a/src/MNH/write_lfifm1_for_diag.f90
+++ b/src/MNH/write_lfifm1_for_diag.f90
@@ -178,7 +178,7 @@ use modd_field,             only: tfieldmetadata, tfieldlist, TYPEINT, TYPEREAL
 USE MODD_FIELD_n,           ONLY: XCIT, XCLDFR, XPABSM, XPABST, XRT, XSIGS, XSRCT, XSVT, XTHT, XTKET, XUT, XVT, XWT, XZWS
 USE MODD_FRC,               ONLY: NFRC, XGXTHFRC, XGYTHFRC, XPGROUNDFRC, XRVFRC, XTENDRVFRC, XTENDTHFRC, XTHFRC, XUFRC, XVFRC, XWFRC
 USE MODD_GRID,              ONLY: XBETA, XLAT0, XLATORI, XLON0, XLONORI, XRPK
-USE MODD_GRID_n,            only: LSLEVE, XLAT, XLEN1, XLEN2, XLON, XZS, XXHAT, XYHAT, XZHAT, XZSMT, XZTOP, XZZ
+USE MODD_GRID_n,            only: LSLEVE, XLAT, XLEN1, XLEN2, XLON, XZS, XXHAT, XXHATM, XYHAT, XYHATM, XZHAT, XZSMT, XZTOP, XZZ
 USE MODD_IO,                ONLY: TFILEDATA
 USE MODD_LSFIELD_n,         ONLY: XLSRVM, XLSTHM, XLSUM, XLSVM, XLSWM
 USE MODD_LUNIT,             ONLY: TLUOUT0
@@ -1217,10 +1217,7 @@ IF (LTRAJ) THEN
   ! X coordinate
   DO JK=1,IKU
     DO JJ=1,IJU
-      DO JI=1,IIU-1
-       ZWORK31(JI,JJ,JK)=0.5*(XXHAT(JI)+XXHAT(JI+1))
-      END DO
-      ZWORK31(IIU,JJ,JK)=2.*ZWORK31(IIU-1,JJ,JK) - ZWORK31(IIU-2,JJ,JK)
+      ZWORK31(:,JJ,JK) = 1E-3*XXHATM(:)
     END DO
   END DO
 
@@ -1236,15 +1233,12 @@ IF (LTRAJ) THEN
     NDIMS      = 3,                    &
     LTIMEDEP   = .TRUE.                )
 
-  CALL IO_Field_write(TPFILE,TZFIELD,ZWORK31*1e-3)
+  CALL IO_Field_write(TPFILE,TZFIELD,ZWORK31)
 
   ! Y coordinate
   DO JK=1,IKU
     DO JI=1,IIU
-      DO JJ=1,IJU-1
-        ZWORK31(JI,JJ,JK)=0.5*(XYHAT(JJ)+XYHAT(JJ+1))
-      END DO
-      ZWORK31(JI,IJU,JK)=2.*ZWORK31(JI,IJU-1,JK) - ZWORK31(JI,IJU-2,JK)
+      ZWORK31(JI,:,JK) = 1E-3 * XYHATM(:)
     END DO
   END DO
 
@@ -1252,7 +1246,7 @@ IF (LTRAJ) THEN
   TZFIELD%CLONGNAME  = 'Y'
   TZFIELD%CCOMMENT   = 'X_Y_Z_Y coordinate'
 
-  CALL IO_Field_write(TPFILE,TZFIELD,ZWORK31*1e-3)
+  CALL IO_Field_write(TPFILE,TZFIELD,ZWORK31)
 END IF
 !
 ! Passive polluant scalar variables
diff --git a/src/MNH/write_surf_mnh.f90 b/src/MNH/write_surf_mnh.f90
index 216fefe1fd080a021b041b61c98a30e2683279ce..ef001e66192d2e198d90f5a7c05b6c27899f4653 100644
--- a/src/MNH/write_surf_mnh.f90
+++ b/src/MNH/write_surf_mnh.f90
@@ -286,7 +286,7 @@ END SUBROUTINE WRITE_SURFX0_MNH
 !
 USE MODD_CONF_n,        ONLY: CSTORAGE_TYPE
 use modd_field,         only: tfieldmetadata, tfieldlist, TYPEREAL
-USE MODD_GRID_n,        ONLY: XXHAT, XYHAT
+USE MODD_GRID_n,        ONLY: XXHAT, XXHATM, XYHAT, XYHATM
 USE MODD_IO,            ONLY: TFILE_SURFEX
 USE MODD_IO_SURF_MNH,   ONLY :NMASK, CMASK,                          &
                               NIU, NJU, NIB, NJB, NIE, NJE,          &
@@ -447,7 +447,12 @@ IF (      (CSTORAGE_TYPE=='PG' .OR. CSTORAGE_TYPE=='SU')  &
     IF (.NOT. (ASSOCIATED(XXHAT))) THEN
       !Store XXHAT if not yet done (necessary for PREP_PGD program when writing netCDF files)
       ALLOCATE(XXHAT(IIU-2*NHALO))
+      ALLOCATE(XXHATM(IIU-2*NHALO))
       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 )
     END IF
   END IF
   DEALLOCATE(ZW1D)
@@ -478,7 +483,12 @@ ELSE IF ( (CSTORAGE_TYPE=='PG' .OR. CSTORAGE_TYPE=='SU')  &
     IF (.NOT. (ASSOCIATED(XYHAT))) THEN
       !Store XYHAT if not yet done (necessary for PREP_PGD program when writing netCDF files)
       ALLOCATE(XYHAT(IJU-2*NHALO))
+      ALLOCATE(XYHATM(IJU-2*NHALO))
       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 )
     END IF
   END IF
   DEALLOCATE(ZW1D)