From e192c8dda41f59a28444e39aeab9825707f7f2f3 Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Thu, 5 Apr 2018 14:03:43 +0200
Subject: [PATCH] Philippe 05/04/2018: IO: modify IO_SET_KNOWNDIMS_NC4,
IO_WRITE_COORDVAR_NC4 and IO_WRITE_FIELD_ATTR_NC4 to support files from all
programs
---
src/LIB/SURCOUCHE/src/mode_netcdf.f90 | 172 ++++++++++++++++----------
1 file changed, 104 insertions(+), 68 deletions(-)
diff --git a/src/LIB/SURCOUCHE/src/mode_netcdf.f90 b/src/LIB/SURCOUCHE/src/mode_netcdf.f90
index a9c24f4af..5491a790b 100644
--- a/src/LIB/SURCOUCHE/src/mode_netcdf.f90
+++ b/src/LIB/SURCOUCHE/src/mode_netcdf.f90
@@ -120,6 +120,8 @@ END SUBROUTINE CLEANIOCDF
SUBROUTINE IO_SET_KNOWNDIMS_NC4(TPFILE)
+USE MODD_CONF, ONLY: CPROGRAM
+USE MODD_CONF_n, ONLY: CSTORAGE_TYPE
USE MODD_DIM_n, ONLY: NIMAX_ll, NJMAX_ll, NKMAX
USE MODD_PARAMETERS_ll, ONLY: JPHEXT, JPVEXT
@@ -142,10 +144,19 @@ IF (.NOT. ASSOCIATED(PIOCDF%DIM_NI_U)) PIOCDF%DIM_NI_U => GETDIMCDF(TPFILE
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')
+IF (TRIM(CPROGRAM)/='PGD' .AND. TRIM(CPROGRAM)/='NESPGD' .AND. TRIM(CPROGRAM)/='ZOOMPG' &
+ .AND. .NOT.(TRIM(CPROGRAM)=='REAL' .AND. CSTORAGE_TYPE=='SU') ) THEN !condition to detect PREP_SURFEX
+ 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')
+ELSE
+ !PGD and SURFEX files for MesoNH have no vertical levels or time scale
+ !These dimensions are allocated to default values
+ !(they need to be allocated when looking for dimensions of variables)
+ IF (.NOT. ASSOCIATED(PIOCDF%DIM_LEVEL)) ALLOCATE(PIOCDF%DIM_LEVEL)
+ IF (.NOT. ASSOCIATED(PIOCDF%DIM_LEVEL_W)) ALLOCATE(PIOCDF%DIM_LEVEL_W)
+ IF (.NOT. ASSOCIATED(PIOCDF%DIMTIME)) ALLOCATE(PIOCDF%DIMTIME)
+END IF
! Store X,Y,Z coordinates for the Arakawa points
! Mass point
@@ -193,6 +204,7 @@ USE MODD_PARAMETERS, ONLY: JPHEXT, JPVEXT
USE MODE_FIELD, ONLY: TFIELDLIST,FIND_FIELD_ID_FROM_MNHNAME
USE MODE_GRIDPROJ
+USE MODE_NEST_ll, ONLY: GET_MODEL_NUMBER_ll, GO_TOMODEL_ll
TYPE(TFILEDATA), INTENT(IN) :: TPFILE
CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: HPROGRAM_ORIG !To emulate a file coming from this program
@@ -200,16 +212,26 @@ CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: HPROGRAM_ORIG !To emulate a file coming
CHARACTER(LEN=:),ALLOCATABLE :: YSTDNAMEPREFIX
CHARACTER(LEN=:),ALLOCATABLE :: YPROGRAM
INTEGER :: IIU, IJU, IKU
-INTEGER :: ID, IRESP
+INTEGER :: ID, IID, IRESP
+INTEGER :: IMI
INTEGER(KIND=IDCDF_KIND) :: INCID
+LOGICAL :: GCHANGEMODEL
+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 :: ZLAT, ZLON
TYPE(IOCDF), POINTER :: PIOCDF
CALL PRINT_MSG(NVERB_DEBUG,'IO','IO_WRITE_COORDVAR_NC4','called for '//TRIM(TPFILE%CNAME))
+ZXHAT => NULL()
+ZYHAT => NULL()
+ZZHAT => NULL()
+
PIOCDF => TPFILE%TNCDIMS
+GCHANGEMODEL = .FALSE.
+
IF (PRESENT(HPROGRAM_ORIG)) THEN
YPROGRAM = HPROGRAM_ORIG
ELSE
@@ -219,87 +241,103 @@ ENDIF
! Get the Netcdf file ID
INCID = TPFILE%NNCID
-IF ( TRIM(YPROGRAM)/='PGD' &
- .AND. .NOT.(TRIM(YPROGRAM)=='REAL' .AND. CSTORAGE_TYPE=='SU') ) THEN !condition to detect PREP_SURFEX
- IIU = SIZE(XXHAT)
- IJU = SIZE(XYHAT)
- 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*(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)
+IF (TPFILE%NMODEL>0) THEN
+ CALL FIND_FIELD_ID_FROM_MNHNAME('XHAT',IID,IRESP)
+ 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('ZHAT',IID,IRESP)
+ ZZHAT => TFIELDLIST(IID)%TFIELD_X1D(TPFILE%NMODEL)%DATA
+ CALL FIND_FIELD_ID_FROM_MNHNAME('SLEVE',IID,IRESP)
+ GSLEVE => TFIELDLIST(IID)%TFIELD_L0D(TPFILE%NMODEL)%DATA
!
- IF (LCARTESIAN) THEN
- YSTDNAMEPREFIX = 'plane'
- ELSE
- YSTDNAMEPREFIX = 'projection'
- ENDIF
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NI,'x-dimension of the grid',TRIM(YSTDNAMEPREFIX)//'_x_coordinate','X',0.,JPHEXT,JPHEXT,ZXHATM)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ,'y-dimension of the grid',TRIM(YSTDNAMEPREFIX)//'_y_coordinate','Y',0.,JPHEXT,JPHEXT,ZYHATM)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_U,'x-dimension of the grid at u location', &
- TRIM(YSTDNAMEPREFIX)//'_x_coordinate_at_u_location','X',-0.5,JPHEXT,0, XXHAT)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_U,'y-dimension of the grid at u location', &
- TRIM(YSTDNAMEPREFIX)//'_y_coordinate_at_u_location','Y', 0., JPHEXT,JPHEXT,ZYHATM)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_V,'x-dimension of the grid at v location', &
- TRIM(YSTDNAMEPREFIX)//'_x_coordinate_at_v_location','X', 0., JPHEXT,JPHEXT,ZXHATM)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_V,'y-dimension of the grid at v location', &
- TRIM(YSTDNAMEPREFIX)//'_y_coordinate_at_v_location','Y',-0.5,JPHEXT,0, XYHAT)
- !
- IF (.NOT.LCARTESIAN) THEN
- ALLOCATE(ZLAT(IIU,IJU),ZLON(IIU,IJU))
- !
- !Compute latitude/longitude for the Arakawa points
- !
- ! Mass point
- CALL WRITE_HOR_2DCOORD(ZXHATM,ZYHATM,'latitude', 'longitude')
- ! u point
- CALL WRITE_HOR_2DCOORD(XXHAT, ZYHATM,'latitude_u','longitude_u')
- ! v point
- CALL WRITE_HOR_2DCOORD(ZXHATM,XYHAT, 'latitude_v','longitude_v')
- ! xi vorticity point (=f point =uv point)
- CALL WRITE_HOR_2DCOORD(XXHAT, XYHAT, 'latitude_f','longitude_f')
- !
- DEALLOCATE(ZLAT,ZLON)
+ CALL GET_MODEL_NUMBER_ll(IMI)
+ IF (IMI/=TPFILE%NMODEL) THEN
+ !This is necessary to have correct domain sizes (used by GATHER_XXFIELD)
+ CALL GO_TOMODEL_ll(TPFILE%NMODEL,IRESP)
+ GCHANGEMODEL = .TRUE.
END IF
- !
- DEALLOCATE(ZXHATM,ZYHATM)
ELSE
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NI,'x-dimension of the grid','x_grid_index','X',0.,JPHEXT,JPHEXT)
- CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ,'y-dimension of the grid','y_grid_index','Y',0.,JPHEXT,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,0)
- 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,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,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,0)
+ ZXHAT => XXHAT
+ ZYHAT => XYHAT
+ ZZHAT => XZHAT
+ 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'
+ELSE
+ YSTDNAMEPREFIX = 'projection'
ENDIF
+CALL WRITE_HOR_COORD(PIOCDF%DIM_NI,'x-dimension of the grid',TRIM(YSTDNAMEPREFIX)//'_x_coordinate','X',0.,JPHEXT,JPHEXT,ZXHATM)
+CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ,'y-dimension of the grid',TRIM(YSTDNAMEPREFIX)//'_y_coordinate','Y',0.,JPHEXT,JPHEXT,ZYHATM)
+CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_U,'x-dimension of the grid at u location', &
+ TRIM(YSTDNAMEPREFIX)//'_x_coordinate_at_u_location','X',-0.5,JPHEXT,0, ZXHAT)
+CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_U,'y-dimension of the grid at u location', &
+ TRIM(YSTDNAMEPREFIX)//'_y_coordinate_at_u_location','Y', 0., JPHEXT,JPHEXT,ZYHATM)
+CALL WRITE_HOR_COORD(PIOCDF%DIM_NI_V,'x-dimension of the grid at v location', &
+ TRIM(YSTDNAMEPREFIX)//'_x_coordinate_at_v_location','X', 0., JPHEXT,JPHEXT,ZXHATM)
+CALL WRITE_HOR_COORD(PIOCDF%DIM_NJ_V,'y-dimension of the grid at v location', &
+ TRIM(YSTDNAMEPREFIX)//'_y_coordinate_at_v_location','Y',-0.5,JPHEXT,0, ZYHAT)
+!
+IF (.NOT.LCARTESIAN) THEN
+ ALLOCATE(ZLAT(IIU,IJU),ZLON(IIU,IJU))
+ !
+ !Compute latitude/longitude for the Arakawa points
+ !
+ ! Mass point
+ CALL WRITE_HOR_2DCOORD(ZXHATM,ZYHATM,'latitude', 'longitude')
+ ! u point
+ CALL WRITE_HOR_2DCOORD(ZXHAT, ZYHATM,'latitude_u','longitude_u')
+ ! v point
+ CALL WRITE_HOR_2DCOORD(ZXHATM,ZYHAT, 'latitude_v','longitude_v')
+ ! xi vorticity point (=f point =uv point)
+ CALL WRITE_HOR_2DCOORD(ZXHAT, ZYHAT, 'latitude_f','longitude_f')
+ !
+ DEALLOCATE(ZLAT,ZLON)
+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' &
+ 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(XZHAT)
+ IKU = SIZE(ZZHAT)
ALLOCATE(ZZHATM(IKU))
- ZZHATM(1:IKU-1) = 0.5 * (XZHAT(2:IKU)+XZHAT(1:IKU-1))
- ZZHATM(IKU) = 2.* XZHAT(IKU) - ZZHATM(IKU-1)
+ 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(PIOCDF%DIM_LEVEL, 'position z in the transformed space', '', &
'altitude', 0., JPVEXT,JPVEXT,ZZHATM)
!
CALL WRITE_VER_COORD(PIOCDF%DIM_LEVEL_W,'position z in the transformed space at w location','', &
- 'altitude_at_w_location',-0.5,JPVEXT,0, XZHAT)
+ 'altitude_at_w_location',-0.5,JPVEXT,0, ZZHAT)
!
DEALLOCATE(ZZHATM)
END IF
END IF
!
!Write time scale
-IF (TPFILE%LMASTER) THEN
- CALL WRITE_TIME_COORD(PIOCDF%DIMTIME)
+IF (TPFILE%LMASTER) THEN !Time scale is 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
+ CALL WRITE_TIME_COORD(PIOCDF%DIMTIME)
+ END IF
END IF
+IF (GCHANGEMODEL) CALL GO_TOMODEL_ll(IMI,IRESP)
+
CONTAINS
SUBROUTINE WRITE_HOR_COORD(TDIM,HLONGNAME,HSTDNAME,HAXIS,PSHIFT,KBOUNDLOW,KBOUNDHIGH,PCOORDS)
USE MODE_ALLOCBUFFER_ll, ONLY: ALLOCBUFFER_ll
@@ -494,7 +532,7 @@ SUBROUTINE WRITE_VER_COORD(TDIM,HLONGNAME,HSTDNAME,HCOMPNAME,PSHIFT,KBOUNDLOW,KB
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))
!
- IF (LSLEVE) THEN
+ IF (GSLEVE) THEN
!Remark: ZS, ZSMT 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 oro_ls: ZSMT oro: ZS len1: LEN1 len2: LEN2')
@@ -744,9 +782,7 @@ IF(PRESENT(HCALENDAR)) THEN
ENDIF
!
! Coordinates (CF convention)
-IF (.NOT.GISCOORD &
- .AND. TRIM(CPROGRAM)/='PGD' &
- .AND. .NOT.(TRIM(CPROGRAM)=='REAL' .AND. CSTORAGE_TYPE=='SU' ) ) THEN !condition to detect PREP_SURFEX
+IF (.NOT.GISCOORD) THEN
!0D: nothing to do
!1D: no direct correspondance with latitude(_x)/longitude(_x) 2D variables => nothing to do
IF (.NOT.LCARTESIAN .AND. TPFIELD%NDIMS>1 .AND. TPFIELD%NGRID/=0) THEN
--
GitLab