From c9bd25529fb0ea4d60be143d7870208d5099c623 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Thu, 11 Jan 2018 10:54:11 +0100
Subject: [PATCH] Philippe 11/01/2018: IO: CF+COMODO conventions: create
coordinate variables (preliminary version)
---
src/LIB/SURCOUCHE/src/modd_netcdf.f90 | 16 ++-
src/LIB/SURCOUCHE/src/mode_fm.f90 | 2 +
src/LIB/SURCOUCHE/src/mode_io.f90 | 10 +-
src/LIB/SURCOUCHE/src/mode_netcdf.f90 | 172 ++++++++++++++++++++++----
4 files changed, 164 insertions(+), 36 deletions(-)
diff --git a/src/LIB/SURCOUCHE/src/modd_netcdf.f90 b/src/LIB/SURCOUCHE/src/modd_netcdf.f90
index 58fa97adf..0eb3ac71a 100644
--- a/src/LIB/SURCOUCHE/src/modd_netcdf.f90
+++ b/src/LIB/SURCOUCHE/src/modd_netcdf.f90
@@ -4,11 +4,17 @@ IMPLICIT NONE
INTEGER,PARAMETER :: IDCDF_KIND = SELECTED_INT_KIND(8)
TYPE IOCDF
- TYPE(DIMCDF), POINTER :: DIMX
- TYPE(DIMCDF), POINTER :: DIMY
- TYPE(DIMCDF), POINTER :: DIMZ
- TYPE(DIMCDF), POINTER :: DIMSTR
- TYPE(DIMCDF), POINTER :: DIMLIST
+ TYPE(DIMCDF), POINTER :: DIM_NI => NULL()
+ TYPE(DIMCDF), POINTER :: DIM_NJ => NULL()
+ TYPE(DIMCDF), POINTER :: DIM_LEVEL => NULL()
+ TYPE(DIMCDF), POINTER :: DIM_NI_U => NULL()
+ TYPE(DIMCDF), POINTER :: DIM_NJ_U => NULL()
+ TYPE(DIMCDF), POINTER :: DIM_NI_V => NULL()
+ TYPE(DIMCDF), POINTER :: DIM_NJ_V => NULL()
+ TYPE(DIMCDF), POINTER :: DIM_LEVEL_W => NULL()
+ TYPE(DIMCDF), POINTER :: DIMTIME => NULL()
+ TYPE(DIMCDF), POINTER :: DIMSTR => NULL()
+ TYPE(DIMCDF), POINTER :: DIMLIST => NULL()
END TYPE IOCDF
TYPE DIMCDF
diff --git a/src/LIB/SURCOUCHE/src/mode_fm.f90 b/src/LIB/SURCOUCHE/src/mode_fm.f90
index ff77129ff..8cf959a71 100644
--- a/src/LIB/SURCOUCHE/src/mode_fm.f90
+++ b/src/LIB/SURCOUCHE/src/mode_fm.f90
@@ -343,6 +343,8 @@ IF (TPFILE%LMASTER) THEN
PRINT *, 'Error in opening (FMOPEN_ll/NF90_CREATE) ', TRIM(YFILEM)//'.nc', ' : ', NF90_STRERROR(INCERR)
STOP
END IF
+ CALL IO_SET_KNOWNDIMS_NC4(TPFILE)
+ CALL IO_WRITE_COORDVAR_NC4(TPFILE)
END IF
END IF
#endif
diff --git a/src/LIB/SURCOUCHE/src/mode_io.f90 b/src/LIB/SURCOUCHE/src/mode_io.f90
index 8cfa6aefb..a8df4bfb7 100644
--- a/src/LIB/SURCOUCHE/src/mode_io.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io.f90
@@ -646,8 +646,8 @@ CONTAINS
CALL PRINT_MSG(NVERB_DEBUG,'IO','OPEN_ll','NF90_OPEN(IO_ZSPLIT) for '//TRIM(TPFILE%CNAME)//CFILE//'.nc')
IOSCDF = NF90_OPEN(TRIM(TPFILE%CNAME)//CFILE//".nc", NF90_NOWRITE, TZSPLITFILE%NNCID)
IF (IOSCDF /= NF90_NOERR) THEN
- PRINT *, 'Error in opening (NF90_OPEN) ', TRIM(TPFILE%CNAME)//CFILE//'.nc', ' : ', NF90_STRERROR(IOSCDF)
- STOP
+ CALL PRINT_MSG(NVERB_FATAL,'IO','OPEN_ll','NF90_OPEN for '//TRIM(TPFILE%CNAME)//CFILE//'.nc: '// &
+ NF90_STRERROR(IOSCDF))
ELSE
IOS = 0
END IF
@@ -661,11 +661,13 @@ CONTAINS
IOSCDF = NF90_CREATE(TRIM(TPFILE%CNAME)//CFILE//".nc", &
&IOR(NF90_CLOBBER,NF90_NETCDF4), TZSPLITFILE%NNCID)
IF (IOSCDF /= NF90_NOERR) THEN
- PRINT *, 'Error in opening (NF90_CREATE) ', TRIM(TPFILE%CNAME)//CFILE//'.nc', ' : ', NF90_STRERROR(IOSCDF)
- STOP
+ CALL PRINT_MSG(NVERB_FATAL,'IO','OPEN_ll','NF90_CREATE for '//TRIM(TPFILE%CNAME)//CFILE//'.nc: '// &
+ NF90_STRERROR(IOSCDF))
ELSE
IOS = 0
END IF
+ CALL IO_SET_KNOWNDIMS_NC4(TPFILE)
+ CALL IO_WRITE_COORDVAR_NC4(TPFILE)
END IF
END IF
#endif
diff --git a/src/LIB/SURCOUCHE/src/mode_netcdf.f90 b/src/LIB/SURCOUCHE/src/mode_netcdf.f90
index b7f2835b3..d6568ce5d 100644
--- a/src/LIB/SURCOUCHE/src/mode_netcdf.f90
+++ b/src/LIB/SURCOUCHE/src/mode_netcdf.f90
@@ -49,7 +49,8 @@ END INTERFACE IO_READ_FIELD_NC4
! Public from module netcdf
PUBLIC NF90_CLOSE,NF90_OPEN,NF90_CREATE,NF90_NOWRITE,NF90_CLOBBER,NF90_NETCDF4,NF90_NOERR,NF90_STRERROR
! Public from this module :
-PUBLIC NEWIOCDF,CLEANIOCDF,IO_WRITE_FIELD_NC4,IO_READ_FIELD_NC4,IO_WRITE_HEADER_NC4
+PUBLIC NEWIOCDF,CLEANIOCDF,IO_SET_KNOWNDIMS_NC4,IO_WRITE_COORDVAR_NC4, &
+ IO_WRITE_FIELD_NC4,IO_READ_FIELD_NC4,IO_WRITE_HEADER_NC4
CONTAINS
@@ -74,12 +75,6 @@ IF (IRESP > 0) THEN
STOP
END IF
-NULLIFY(TZIOCDF%DIMX)
-NULLIFY(TZIOCDF%DIMY)
-NULLIFY(TZIOCDF%DIMZ)
-NULLIFY(TZIOCDF%DIMSTR)
-NULLIFY(TZIOCDF%DIMLIST)
-
NEWIOCDF=>TZIOCDF
END FUNCTION NEWIOCDF
@@ -113,6 +108,124 @@ END SUBROUTINE CLEANLIST
END SUBROUTINE CLEANIOCDF
+
+SUBROUTINE IO_SET_KNOWNDIMS_NC4(TPFILE)
+
+USE MODD_DIM_ll, ONLY: NIMAX_ll, NJMAX_ll, NKMAX_ll
+USE MODD_PARAMETERS_ll, ONLY: JPHEXT, JPVEXT
+
+TYPE(TFILEDATA),INTENT(IN) :: TPFILE
+
+INTEGER :: IIU_ll, IJU_ll, IKU
+TYPE(IOCDF), POINTER :: PIOCDF
+
+CALL PRINT_MSG(NVERB_DEBUG,'IO','SET_KNOWN_NC_DIMS','called')
+
+PIOCDF => TPFILE%TNCDIMS
+
+IIU_ll = NIMAX_ll + 2*JPHEXT
+IJU_ll = NJMAX_ll + 2*JPHEXT
+IKU = NKMAX_ll + 2*JPVEXT
+
+IF (.NOT. ASSOCIATED(PIOCDF%DIM_NI)) PIOCDF%DIM_NI => GETDIMCDF(TPFILE, IIU_ll, 'ni')
+IF (.NOT. ASSOCIATED(PIOCDF%DIM_NJ)) PIOCDF%DIM_NJ => GETDIMCDF(TPFILE, IJU_ll, 'nj')
+IF (.NOT. ASSOCIATED(PIOCDF%DIM_NI_U)) PIOCDF%DIM_NI_U => GETDIMCDF(TPFILE, IIU_ll, 'ni_u')
+IF (.NOT. ASSOCIATED(PIOCDF%DIM_NJ_U)) PIOCDF%DIM_NJ_U => GETDIMCDF(TPFILE, IJU_ll, 'nj_u')
+IF (.NOT. ASSOCIATED(PIOCDF%DIM_NI_V)) PIOCDF%DIM_NI_V => GETDIMCDF(TPFILE, IIU_ll, 'ni_v')
+IF (.NOT. ASSOCIATED(PIOCDF%DIM_NJ_V)) PIOCDF%DIM_NJ_V => GETDIMCDF(TPFILE, IJU_ll, 'nj_v')
+IF (.NOT. ASSOCIATED(PIOCDF%DIM_LEVEL)) PIOCDF%DIM_LEVEL => GETDIMCDF(TPFILE, IKU , 'level')
+IF (.NOT. ASSOCIATED(PIOCDF%DIM_LEVEL_W)) PIOCDF%DIM_LEVEL_W => GETDIMCDF(TPFILE, IKU , 'level_w')
+
+IF (.NOT. ASSOCIATED(PIOCDF%DIMTIME)) PIOCDF%DIMTIME => GETDIMCDF(TPFILE, NF90_UNLIMITED, 'time')
+
+END SUBROUTINE IO_SET_KNOWNDIMS_NC4
+
+
+SUBROUTINE IO_WRITE_COORDVAR_NC4(TPFILE)
+USE MODD_PARAMETERS, ONLY: JPHEXT
+
+TYPE(TFILEDATA),INTENT(IN) :: TPFILE
+
+INTEGER(KIND=IDCDF_KIND) :: INCID
+TYPE(IOCDF), POINTER :: PIOCDF
+
+CALL PRINT_MSG(NVERB_DEBUG,'IO','WRITE_NC_COORDS_VAR','called')
+
+PIOCDF => TPFILE%TNCDIMS
+
+! Get the Netcdf file ID
+INCID = TPFILE%NNCID
+
+IF (TPFILE%LMASTER) THEN
+ CALL WRITE_COORD(PIOCDF%DIM_NI,'x-dimension of the grid','x_grid_index','X',0.,JPHEXT)
+ CALL WRITE_COORD(PIOCDF%DIM_NJ,'y-dimension of the grid','y_grid_index','Y',0.,JPHEXT)
+ CALL WRITE_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_COORD(PIOCDF%DIM_NJ_U,'y-dimension of the grid at u location','y_grid_index_at_u_location','Y', 0., JPHEXT)
+ CALL WRITE_COORD(PIOCDF%DIM_NI_V,'x-dimension of the grid at v location','x_grid_index_at_v_location','X', 0., JPHEXT)
+ CALL WRITE_COORD(PIOCDF%DIM_NJ_V,'y-dimension of the grid at v location','y_grid_index_at_v_location','Y',-0.5,JPHEXT)
+ENDIF
+
+CONTAINS
+SUBROUTINE WRITE_COORD(TDIM,HLONGNAME,HSTDNAME,HAXIS,PSHIFT,KBOUND)
+ 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
+
+ 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
+
+ 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
+ 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))
+ 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
+ STATUS = NF90_PUT_VAR(INCID, IVARID, ZTAB)
+ IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'WRITE_NC_COORDS_VAR[NF90_PUT_VAR] '//TRIM(YVARNAME))
+
+ DEALLOCATE(ZTAB)
+
+END SUBROUTINE WRITE_COORD
+
+END SUBROUTINE IO_WRITE_COORDVAR_NC4
+
+
SUBROUTINE HANDLE_ERR(STATUS,LINE,TEXT,KRESP)
INTEGER(KIND=IDCDF_KIND),INTENT(IN) :: STATUS
INTEGER, INTENT(IN) :: LINE
@@ -267,11 +380,20 @@ CHARACTER(LEN=7) :: YSUFFIX
CHARACTER(LEN=8) :: YDIMNAME
CHARACTER(LEN=20) :: YLEN
INTEGER(KIND=IDCDF_KIND) :: STATUS
+LOGICAL :: GCHKLEN !Check if KLEN is valid
TYPE(IOCDF), POINTER :: PIOCDF
PIOCDF => TPFILE%TNCDIMS
-IF (KLEN < 1) THEN
+GCHKLEN = .TRUE.
+!Do not check KLEN if 'time' (because NF90_UNLIMITED = 0)
+IF (PRESENT(HDIMNAME)) THEN
+ IF (TRIM(HDIMNAME)=='time') THEN
+ GCHKLEN = .FALSE.
+ END IF
+END IF
+
+IF (GCHKLEN .AND. KLEN < 1) THEN
WRITE(YLEN,*) KLEN
CALL PRINT_MSG(NVERB_FATAL,'IO','GETDIMCDF','KLEN='//TRIM(YLEN))
END IF
@@ -367,16 +489,14 @@ DO II=1, SIZE(KSHAPE)
IF (II == 1) THEN
IF (HDIR == 'XX' .OR. HDIR == 'XY') THEN
- IF (.NOT. ASSOCIATED(PIOCDF%DIMX)) PIOCDF%DIMX => GETDIMCDF(TPFILE, KSHAPE(II), 'X')
- IF (KSHAPE(II) == PIOCDF%DIMX%LEN) THEN
- PTDIM => PIOCDF%DIMX
+ IF (KSHAPE(II) == PIOCDF%DIM_NI%LEN) THEN
+ PTDIM => PIOCDF%DIM_NI
ELSE
PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
END IF
ELSE IF (HDIR == 'YY') THEN
- IF (.NOT. ASSOCIATED(PIOCDF%DIMY)) PIOCDF%DIMY => GETDIMCDF(TPFILE, KSHAPE(II), 'Y')
- IF (KSHAPE(II) == PIOCDF%DIMY%LEN) THEN
- PTDIM => PIOCDF%DIMY
+ IF (KSHAPE(II) == PIOCDF%DIM_NJ%LEN) THEN
+ PTDIM => PIOCDF%DIM_NJ
ELSE
PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
END IF
@@ -386,9 +506,8 @@ DO II=1, SIZE(KSHAPE)
END IF
ELSE IF (II == 2) THEN
IF (HDIR == 'XY') THEN
- IF (.NOT. ASSOCIATED(PIOCDF%DIMY)) PIOCDF%DIMY => GETDIMCDF(TPFILE, KSHAPE(II), 'Y')
- IF (KSHAPE(II) == PIOCDF%DIMY%LEN) THEN
- PTDIM => PIOCDF%DIMY
+ IF (KSHAPE(II) == PIOCDF%DIM_NJ%LEN) THEN
+ PTDIM => PIOCDF%DIM_NJ
ELSE
PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
END IF
@@ -397,9 +516,8 @@ DO II=1, SIZE(KSHAPE)
END IF
ELSE IF (II == 3) THEN
IF (HDIR == 'XY') THEN
- IF (.NOT. ASSOCIATED(PIOCDF%DIMZ)) PIOCDF%DIMZ => GETDIMCDF(TPFILE, KSHAPE(II), 'Z')
- IF (KSHAPE(II) == PIOCDF%DIMZ%LEN) THEN
- PTDIM => PIOCDF%DIMZ
+ IF (KSHAPE(II) == PIOCDF%DIM_LEVEL%LEN) THEN
+ PTDIM => PIOCDF%DIM_LEVEL
ELSE
PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
END IF
@@ -889,21 +1007,21 @@ STATUS = NF90_PUT_VAR(INCID, IVARID, KFIELD)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__,'IO_WRITE_FIELD_NC4_N0[NF90_PUT_VAR] '//TRIM(TPFIELD%CMNHNAME),IRESP)
!
-! Use IMAX, JMAX, KMAX to define DIMX, DIMY, DIMZ
+! Use IMAX, JMAX, KMAX to define DIM_NI, DIM_NJ, DIM_LEVEL
! /!\ Can only work if IMAX, JMAX or KMAX are written before any array
!
#if 0
-IF (YVARNAME == 'IMAX' .AND. .NOT. ASSOCIATED(TPFILE%TNCDIMS%DIMX)) TPFILE%TNCDIMS%DIMX=>GETDIMCDF(TPFILE%TNCDIMS,KFIELD+2*JPHEXT,'X')
-IF (YVARNAME == 'JMAX' .AND. .NOT. ASSOCIATED(TPFILE%TNCDIMS%DIMY)) THEN
+IF (YVARNAME == 'IMAX' .AND. .NOT. ASSOCIATED(TPFILE%TNCDIMS%DIM_NI)) TPFILE%TNCDIMS%DIM_NI=>GETDIMCDF(TPFILE%TNCDIMS,KFIELD+2*JPHEXT,'X')
+IF (YVARNAME == 'JMAX' .AND. .NOT. ASSOCIATED(TPFILE%TNCDIMS%DIM_NJ)) THEN
IF (LPACK .AND. L2D) THEN
- TPFILE%TNCDIMS%DIMY=>GETDIMCDF(TPFILE, 1,'Y')
+ TPFILE%TNCDIMS%DIM_NJ=>GETDIMCDF(TPFILE, 1,'Y')
ELSE
- TPFILE%TNCDIMS%DIMY=>GETDIMCDF(TPFILE, KFIELD+2*JPHEXT, 'Y')
+ TPFILE%TNCDIMS%DIM_NJ=>GETDIMCDF(TPFILE, KFIELD+2*JPHEXT, 'Y')
END IF
END IF
#endif
-IF (YVARNAME == 'KMAX' .AND. .NOT. ASSOCIATED(TPFILE%TNCDIMS%DIMZ)) &
- TPFILE%TNCDIMS%DIMZ=>GETDIMCDF(TPFILE,INT(KFIELD+2*JPVEXT,KIND=IDCDF_KIND),'Z')
+IF (YVARNAME == 'KMAX' .AND. .NOT. ASSOCIATED(TPFILE%TNCDIMS%DIM_LEVEL)) &
+ TPFILE%TNCDIMS%DIM_LEVEL=>GETDIMCDF(TPFILE,INT(KFIELD+2*JPVEXT,KIND=IDCDF_KIND),'Z')
KRESP = IRESP
END SUBROUTINE IO_WRITE_FIELD_NC4_N0
--
GitLab