diff --git a/src/LIB/SURCOUCHE/src/mode_netcdf.f90 b/src/LIB/SURCOUCHE/src/mode_netcdf.f90
index be013b70df9c178744e959b70a97fd1554d6c801..3740c2a282f1492e7aaf773439038625e28c2d3b 100644
--- a/src/LIB/SURCOUCHE/src/mode_netcdf.f90
+++ b/src/LIB/SURCOUCHE/src/mode_netcdf.f90
@@ -144,14 +144,18 @@ END SUBROUTINE IO_SET_KNOWNDIMS_NC4
SUBROUTINE IO_WRITE_COORDVAR_NC4(TPFILE)
-USE MODD_PARAMETERS, ONLY: JPHEXT
+USE MODD_CONF, ONLY: CPROGRAM
+USE MODD_GRID_n, ONLY: LSLEVE, XXHAT, XYHAT, XZHAT
+USE MODD_PARAMETERS, ONLY: JPHEXT, JPVEXT
TYPE(TFILEDATA),INTENT(IN) :: TPFILE
-INTEGER(KIND=IDCDF_KIND) :: INCID
-TYPE(IOCDF), POINTER :: PIOCDF
+INTEGER :: IIU, IJU, IKU
+INTEGER(KIND=IDCDF_KIND) :: INCID
+REAL,DIMENSION(:),ALLOCATABLE :: ZXHATM, ZYHATM,ZZHATM !Coordinates at mass points in the transformed space
+TYPE(IOCDF), POINTER :: PIOCDF
-CALL PRINT_MSG(NVERB_DEBUG,'IO','WRITE_NC_COORDS_VAR','called')
+CALL PRINT_MSG(NVERB_DEBUG,'IO','IO_WRITE_COORDVAR_NC4','called')
PIOCDF => TPFILE%TNCDIMS
@@ -159,22 +163,127 @@ PIOCDF => TPFILE%TNCDIMS
INCID = TPFILE%NNCID
IF (TPFILE%LMASTER) THEN
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NI,'x-dimension of the grid','x_grid_index','X',0.,JPHEXT)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ,'y-dimension of the grid','y_grid_index','Y',0.,JPHEXT)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_U,'x-dimension of the grid at u location','x_grid_index_at_u_location','X',-0.5,JPHEXT)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_U,'y-dimension of the grid at u location','y_grid_index_at_u_location','Y', 0., JPHEXT)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_V,'x-dimension of the grid at v location','x_grid_index_at_v_location','X', 0., JPHEXT)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_V,'y-dimension of the grid at v location','y_grid_index_at_v_location','Y',-0.5,JPHEXT)
+ IF (TRIM(CPROGRAM)/='PGD') THEN
+ IIU = SIZE(XXHAT)
+ IJU = SIZE(XYHAT)
+ ALLOCATE(ZXHATM(IIU),ZYHATM(IJU))
+ 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)
+ !
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NI,'x-dimension of the grid','projection_x_coordinate','X',0.,JPHEXT,ZXHATM)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ,'y-dimension of the grid','projection_y_coordinate','Y',0.,JPHEXT,ZYHATM)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_U,'x-dimension of the grid at u location','projection_x_coordinate_at_u_location',&
+ 'X',-0.5,JPHEXT,XXHAT)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_U,'y-dimension of the grid at u location','projection_y_coordinate_at_u_location',&
+ 'Y', 0., JPHEXT,ZYHATM)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_V,'x-dimension of the grid at v location','projection_x_coordinate_at_v_location',&
+ 'X', 0., JPHEXT,ZXHATM)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_V,'y-dimension of the grid at v location','projection_y_coordinate_at_v_location',&
+ 'Y',-0.5,JPHEXT,XYHAT)
+ !
+ DEALLOCATE(ZXHATM,ZYHATM)
+ ELSE
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NI,'x-dimension of the grid','x_grid_index','X',0.,JPHEXT)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ,'y-dimension of the grid','y_grid_index','Y',0.,JPHEXT)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_U,'x-dimension of the grid at u location','x_grid_index_at_u_location','X',-0.5,JPHEXT)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_U,'y-dimension of the grid at u location','y_grid_index_at_u_location','Y', 0., JPHEXT)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_V,'x-dimension of the grid at v location','x_grid_index_at_v_location','X', 0., JPHEXT)
+ CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_V,'y-dimension of the grid at v location','y_grid_index_at_v_location','Y',-0.5,JPHEXT)
+ ENDIF
+ !
+ IF (TRIM(CPROGRAM)/='PGD' .AND. TRIM(CPROGRAM)/='NESPGD') THEN
+ IKU = SIZE(XZHAT)
+ ALLOCATE(ZZHATM(IKU))
+ ZZHATM(1:IKU-1) = 0.5 * (XZHAT(2:IKU)+XZHAT(1:IKU-1))
+ ZZHATM(IKU) = 2.* XZHAT(IKU) - ZZHATM(IKU-1)
+ !
+ CALL WRITE_VER_COORD(PIOCDF%DIM_LEVEL, 'position z in the transformed space', '', 0., JPVEXT,ZZHATM)
+ !
+ CALL WRITE_VER_COORD(PIOCDF%DIM_LEVEL_W,'position z in the transformed space at w location','',-0.5,JPVEXT,XZHAT)
+ !
+ DEALLOCATE(ZZHATM)
+ END IF
+ !
ENDIF
+
CONTAINS
-SUBROUTINE WRITE_HOR_COORD(TDIM,HLONGNAME,HSTDNAME,HAXIS,PSHIFT,KBOUND)
+SUBROUTINE WRITE_HOR_COORD(TDIM,HLONGNAME,HSTDNAME,HAXIS,PSHIFT,KBOUND,PCOORDS)
+ TYPE(DIMCDF), POINTER, INTENT(IN) :: TDIM
+ CHARACTER(LEN=*), INTENT(IN) :: HLONGNAME
+ CHARACTER(LEN=*), INTENT(IN) :: HSTDNAME
+ CHARACTER(LEN=*), INTENT(IN) :: HAXIS
+ REAL, INTENT(IN) :: PSHIFT
+ INTEGER, INTENT(IN) :: KBOUND
+ REAL,DIMENSION(:),OPTIONAL,INTENT(IN) :: PCOORDS
+
+ CHARACTER(LEN=64) :: YRANGE
+ CHARACTER(LEN=:),ALLOCATABLE :: YVARNAME
+ INTEGER :: IRESP
+ INTEGER :: ISIZE
+ INTEGER :: JI
+ INTEGER(KIND=IDCDF_KIND) :: IVARID
+ INTEGER(KIND=IDCDF_KIND) :: IVDIM
+ INTEGER(KIND=IDCDF_KIND) :: STATUS
+ REAL,DIMENSION(:),ALLOCATABLE :: ZTAB
+
+ ISIZE = TDIM%LEN
+ YVARNAME = TRIM(TDIM%NAME)
+ IVDIM = TDIM%ID
+
+ IF (.NOT.PRESENT(PCOORDS)) THEN
+ ALLOCATE(ZTAB(ISIZE))
+ DO JI=1,ISIZE
+ ZTAB(JI) = REAL(JI,KIND=KIND(ZTAB(1)))+PSHIFT
+ END DO
+ END IF
+
+ STATUS = NF90_INQ_VARID(INCID, YVARNAME, IVARID)
+ IF (STATUS /= NF90_NOERR) THEN
+ ! Define the coordinate variable
+ STATUS = NF90_DEF_VAR(INCID, YVARNAME, NF90_DOUBLE, IVDIM, IVARID)
+ IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_DEF_VAR]')
+ ELSE
+ CALL PRINT_MSG(NVERB_ERROR,'IO','WRITE_NC_COORDS_VAR',TRIM(YVARNAME)//' already defined')
+ END IF
+
+ ! Write metadata
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'long_name',HLONGNAME)
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'standard_name',HSTDNAME)
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ IF (PRESENT(PCOORDS)) THEN
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'units','m')
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ END IF
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'axis',HAXIS)
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'c_grid_axis_shift',PSHIFT)
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ WRITE(YRANGE,'( I0,":",I0 )') 1+KBOUND,ISIZE-KBOUND
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'c_grid_dynamic_range',TRIM(YRANGE))
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+
+ ! Write the data
+ IF (PRESENT(PCOORDS)) THEN
+ STATUS = NF90_PUT_VAR(INCID, IVARID, PCOORDS)
+ ELSE
+ STATUS = NF90_PUT_VAR(INCID, IVARID, ZTAB)
+ DEALLOCATE(ZTAB)
+ END IF
+ IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+
+END SUBROUTINE WRITE_HOR_COORD
+
+SUBROUTINE WRITE_VER_COORD(TDIM,HLONGNAME,HSTDNAME,PSHIFT,KBOUND,PCOORDS)
TYPE(DIMCDF), POINTER, INTENT(IN) :: TDIM
CHARACTER(LEN=*), INTENT(IN) :: HLONGNAME
CHARACTER(LEN=*), INTENT(IN) :: HSTDNAME
- CHARACTER(LEN=*), INTENT(IN) :: HAXIS
REAL, INTENT(IN) :: PSHIFT
INTEGER, INTENT(IN) :: KBOUND
+ REAL,DIMENSION(:), INTENT(IN) :: PCOORDS
CHARACTER(LEN=64) :: YRANGE
CHARACTER(LEN=:),ALLOCATABLE :: YVARNAME
@@ -184,17 +293,11 @@ SUBROUTINE WRITE_HOR_COORD(TDIM,HLONGNAME,HSTDNAME,HAXIS,PSHIFT,KBOUND)
INTEGER(KIND=IDCDF_KIND) :: IVARID
INTEGER(KIND=IDCDF_KIND) :: IVDIM
INTEGER(KIND=IDCDF_KIND) :: STATUS
- REAL,DIMENSION(:),ALLOCATABLE :: ZTAB
ISIZE = TDIM%LEN
YVARNAME = TRIM(TDIM%NAME)
IVDIM = TDIM%ID
- ALLOCATE(ZTAB(ISIZE))
- DO JI=1,ISIZE
- ZTAB(JI) = REAL(JI,KIND=KIND(ZTAB(1)))+PSHIFT
- END DO
-
STATUS = NF90_INQ_VARID(INCID, YVARNAME, IVARID)
IF (STATUS /= NF90_NOERR) THEN
! Define the coordinate variable
@@ -209,21 +312,30 @@ SUBROUTINE WRITE_HOR_COORD(TDIM,HLONGNAME,HSTDNAME,HAXIS,PSHIFT,KBOUND)
IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
STATUS = NF90_PUT_ATT(INCID, IVARID, 'standard_name',HSTDNAME)
IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
- STATUS = NF90_PUT_ATT(INCID, IVARID, 'axis',HAXIS)
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'units','m')
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'axis','Z')
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'positive','up')
IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
STATUS = NF90_PUT_ATT(INCID, IVARID, 'c_grid_axis_shift',PSHIFT)
IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
WRITE(YRANGE,'( I0,":",I0 )') 1+KBOUND,ISIZE-KBOUND
STATUS = NF90_PUT_ATT(INCID, IVARID, 'c_grid_dynamic_range',TRIM(YRANGE))
IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ !Remark: ZS and ZTOP in the formula are the same for mass point or flux point
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'formula_terms','s: '//TRIM(YVARNAME)//' height: ZTOP orog: ZS')
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'formula_definition','z(n,k,j,i)=s(k)*(height-orog(j,i))/height+orog(j,i)')
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+ STATUS = NF90_PUT_ATT(INCID, IVARID, 'formula_terms','s: '//TRIM(YVARNAME)//' height: ZTOP orog: ZS')
+ IF (STATUS /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
! Write the data
- STATUS = NF90_PUT_VAR(INCID, IVARID, ZTAB)
+ STATUS = NF90_PUT_VAR(INCID, IVARID, PCOORDS)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
- DEALLOCATE(ZTAB)
-
-END SUBROUTINE WRITE_HOR_COORD
+END SUBROUTINE WRITE_VER_COORD
END SUBROUTINE IO_WRITE_COORDVAR_NC4