From 6806bd4399c18338924420529921482aa8bffa0b Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Thu, 18 Jan 2018 16:30:31 +0100
Subject: [PATCH] Philippe 18/01/2018: IO: modify FILLVDIMS to choose right
coordinate variables for fields
---
src/LIB/SURCOUCHE/src/modd_netcdf.f90 | 11 ++-
src/LIB/SURCOUCHE/src/mode_netcdf.f90 | 124 ++++++++++++++------------
2 files changed, 75 insertions(+), 60 deletions(-)
diff --git a/src/LIB/SURCOUCHE/src/modd_netcdf.f90 b/src/LIB/SURCOUCHE/src/modd_netcdf.f90
index 0eb3ac71a..688797691 100644
--- a/src/LIB/SURCOUCHE/src/modd_netcdf.f90
+++ b/src/LIB/SURCOUCHE/src/modd_netcdf.f90
@@ -18,10 +18,13 @@ TYPE IOCDF
END TYPE IOCDF
TYPE DIMCDF
- CHARACTER(LEN=8) :: NAME
- INTEGER(KIND=IDCDF_KIND) :: LEN
- INTEGER(KIND=IDCDF_KIND) :: ID
- TYPE(DIMCDF), POINTER :: NEXT
+ CHARACTER(LEN=8) :: NAME = ''
+ INTEGER(KIND=IDCDF_KIND) :: LEN = 0
+ INTEGER(KIND=IDCDF_KIND) :: ID = -1
+ TYPE(DIMCDF), POINTER :: NEXT => NULL()
END TYPE DIMCDF
+TYPE(DIMCDF),DIMENSION(3,0:4) :: NCOORDID !X,Y,Z coordinates for the 4 Arakawa points
+ !0 2nd-dimension is to treat NGRID=0 case without crash
+
END MODULE MODD_NETCDF
diff --git a/src/LIB/SURCOUCHE/src/mode_netcdf.f90 b/src/LIB/SURCOUCHE/src/mode_netcdf.f90
index 9609697ee..b11e4b0a2 100644
--- a/src/LIB/SURCOUCHE/src/mode_netcdf.f90
+++ b/src/LIB/SURCOUCHE/src/mode_netcdf.f90
@@ -140,6 +140,24 @@ IF (.NOT. ASSOCIATED(PIOCDF%DIM_LEVEL_W)) PIOCDF%DIM_LEVEL_W => GETDIMCDF(TPFILE
IF (.NOT. ASSOCIATED(PIOCDF%DIMTIME)) PIOCDF%DIMTIME => GETDIMCDF(TPFILE, NF90_UNLIMITED, 'time')
+! Store X,Y,Z coordinates for the 4 Arakawa points
+! Mass point
+NCOORDID(1,1) = PIOCDF%DIM_NI
+NCOORDID(2,1) = PIOCDF%DIM_NJ
+NCOORDID(3,1) = PIOCDF%DIM_LEVEL
+! u point
+NCOORDID(1,2) = PIOCDF%DIM_NI_U
+NCOORDID(2,2) = PIOCDF%DIM_NJ_U
+NCOORDID(3,2) = PIOCDF%DIM_LEVEL
+! v point
+NCOORDID(1,3) = PIOCDF%DIM_NI_V
+NCOORDID(2,3) = PIOCDF%DIM_NJ_V
+NCOORDID(3,3) = PIOCDF%DIM_LEVEL
+! w point
+NCOORDID(1,4) = PIOCDF%DIM_NI
+NCOORDID(2,4) = PIOCDF%DIM_NJ
+NCOORDID(3,4) = PIOCDF%DIM_LEVEL_W
+
END SUBROUTINE IO_SET_KNOWNDIMS_NC4
@@ -597,66 +615,60 @@ TYPE(TFILEDATA), INTENT(IN) :: TPFILE
TYPE(TFIELDDATA), INTENT(IN) :: TPFIELD
INTEGER(KIND=IDCDF_KIND),DIMENSION(:),INTENT(IN) :: KSHAPE
INTEGER(KIND=IDCDF_KIND),DIMENSION(:),INTENT(OUT) :: KVDIMS
-
-INTEGER :: II
+!
+INTEGER :: IGRID
+INTEGER :: JI
+CHARACTER(LEN=32) :: YINT
CHARACTER(LEN=2) :: YDIR
TYPE(DIMCDF), POINTER :: PTDIM
-TYPE(IOCDF), POINTER :: PIOCDF
-
-CALL PRINT_MSG(NVERB_DEBUG,'IO','FILLVDIMS','called')
-
+!
+CALL PRINT_MSG(NVERB_DEBUG,'IO','FILLVDIMS','called for '//TRIM(TPFIELD%CMNHNAME))
+!
IF (SIZE(KSHAPE) < 1) CALL PRINT_MSG(NVERB_FATAL,'IO','FILLVDIMS','empty KSHAPE')
-
-YDIR = TPFIELD%CDIR
-PIOCDF => TPFILE%TNCDIMS
-
-DO II=1, SIZE(KSHAPE)
-
- IF (II == 1) THEN
- IF (YDIR == 'XX' .OR. YDIR == 'XY') THEN
- IF (KSHAPE(II) == PIOCDF%DIM_NI%LEN) THEN
- PTDIM => PIOCDF%DIM_NI
- ELSE
- PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
- END IF
- ELSE IF (YDIR == 'YY') THEN
- IF (KSHAPE(II) == PIOCDF%DIM_NJ%LEN) THEN
- PTDIM => PIOCDF%DIM_NJ
- ELSE
- PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
- END IF
- ELSE
- PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
- KVDIMS(II) = PTDIM%ID
- END IF
- ELSE IF (II == 2) THEN
- IF (YDIR == 'XY') THEN
- IF (KSHAPE(II) == PIOCDF%DIM_NJ%LEN) THEN
- PTDIM => PIOCDF%DIM_NJ
- ELSE
- PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
- END IF
- ELSE
- PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
- END IF
- ELSE IF (II == 3) THEN
- IF (YDIR == 'XY') THEN
- IF (KSHAPE(II) == PIOCDF%DIM_LEVEL%LEN) THEN
- PTDIM => PIOCDF%DIM_LEVEL
- ELSE
- PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
- END IF
- ELSE
- PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
- END IF
- ELSE
- PTDIM => GETDIMCDF(TPFILE, KSHAPE(II))
- END IF
-
- KVDIMS(II) = PTDIM%ID
-
+!
+IGRID = TPFIELD%NGRID
+YDIR = TPFIELD%CDIR
+!
+IF(SIZE(KSHAPE)/=TPFIELD%NDIMS) THEN
+ WRITE(YINT,'( I0,"/",I0 )') SIZE(KSHAPE),TPFIELD%NDIMS
+ CALL PRINT_MSG(NVERB_FATAL,'IO','FILLVDIMS','SIZE(KSHAPE)/=TPFIELD%NDIMS ('//TRIM(YINT)//') for field '//TRIM(TPFIELD%CMNHNAME))
+END IF
+!
+IF(IGRID<0 .OR. IGRID>4) THEN
+ WRITE(YINT,'( I0 )') IGRID
+ CALL PRINT_MSG(NVERB_FATAL,'IO','FILLVDIMS','invalid NGRID ('//TRIM(YINT)//') for field '//TRIM(TPFIELD%CMNHNAME))
+END IF
+!
+IF(IGRID==0 .AND. YDIR/='--' .AND. YDIR/='' ) THEN
+ CALL PRINT_MSG(NVERB_FATAL,'IO','FILLVDIMS','invalid YDIR ('//TRIM(YDIR)//') with NGRID=0 for field '//TRIM(TPFIELD%CMNHNAME))
+END IF
+!
+DO JI=1,TPFIELD%NDIMS
+ IF (JI == 1) THEN
+ IF ( (YDIR == 'XX' .OR. YDIR == 'XY') .AND. KSHAPE(1)==NCOORDID(1,IGRID)%LEN) THEN
+ KVDIMS(1) = NCOORDID(1,IGRID)%ID
+ ELSE IF ( YDIR == 'YY' .AND. KSHAPE(1)==NCOORDID(1,IGRID)%LEN) THEN
+ KVDIMS(1) = NCOORDID(2,IGRID)%ID
+ ELSE
+ PTDIM => GETDIMCDF(TPFILE, KSHAPE(1)); KVDIMS(1) = PTDIM%ID
+ END IF
+ ELSE IF (JI == 2) THEN
+ IF ( YDIR == 'XY' .AND. KSHAPE(2)==NCOORDID(2,IGRID)%LEN) THEN
+ KVDIMS(2) = NCOORDID(2,IGRID)%ID
+ ELSE
+ PTDIM => GETDIMCDF(TPFILE, KSHAPE(2)); KVDIMS(2) = PTDIM%ID
+ END IF
+ ELSE IF (JI == 3) THEN
+ IF ( YDIR == 'XY' .AND. KSHAPE(3)==NCOORDID(3,IGRID)%LEN) THEN
+ KVDIMS(3) = NCOORDID(3,IGRID)%ID
+ ELSE
+ PTDIM => GETDIMCDF(TPFILE, KSHAPE(3)); KVDIMS(3) = PTDIM%ID
+ END IF
+ ELSE
+ PTDIM => GETDIMCDF(TPFILE, KSHAPE(JI)); KVDIMS(JI) = PTDIM%ID
+ END IF
END DO
-
+!
END SUBROUTINE FILLVDIMS
--
GitLab