Maintenance is scheduled between 16:00 and 23:59 CEST (14:00 and 21:59 UTC) on Thursday 2021-10-28. The system may be unavailable at any time during this timeframe. Please plan accordingly.

Commit 8624a75e authored by Don Morton's avatar Don Morton
Browse files

Enhancements to FPv9.3.2

Documented in Ticket #182 (as well as CTBTO ticket 357)
parent 4c0504c3
This diff is collapsed.
...@@ -10,37 +10,41 @@ GRIBAPI = /opt/grib-api ...@@ -10,37 +10,41 @@ GRIBAPI = /opt/grib-api
HDF5 = /opt/hdf5-1.8.16 HDF5 = /opt/hdf5-1.8.16
NETCDFF = /opt/netcdf-fortran-4.4.3 NETCDFF = /opt/netcdf-fortran-4.4.3
NETCDF = /opt/netcdf-c-4.4.0 NETCDF = /opt/netcdf-c-4.4.0
#GRIBAPI = /usr/local/grib-api
#HDF5 = /usr/local/hdf5-1.8.16
#NETCDFF = /usr/local/netcdf-fortran-4.4.3
#NETCDF = /usr/local/netcdf-c-4.4.0
BINARY = grib2nc4 BINARY = grib2nc4
OBJS = processmetfields.o CMP_BINARY = nc4cmp
FPMODOBJS = par_mod.o com_mod.o class_vtable_mod.o cmapf_mod.o conv_mod.o OBJS = processmetfields.o verttransform_grib2nc4_ecmwf.o verttransform_grib2nc4_gfs.o
FLXPRTOBJS = detectformat.o grib2check.o shift_field_0.o gridcheck.o \ FPMODOBJS = ${FLEXPART_SRC}/par_mod.o ${FLEXPART_SRC}/com_mod.o ${FLEXPART_SRC}/class_vtable_mod.o ${FLEXPART_SRC}/cmapf_mod.o ${FLEXPART_SRC}/conv_mod.o
readwind.o readwind_nests.o calcpar.o calcpar_nests.o \ FLXPRTOBJS = ${FLEXPART_SRC}/detectformat.o ${FLEXPART_SRC}/grib2check.o ${FLEXPART_SRC}/shift_field_0.o ${FLEXPART_SRC}/gridcheck.o \
shift_field.o pbl_profile.o scalev.o obukhov.o \ ${FLEXPART_SRC}/readwind.o ${FLEXPART_SRC}/readwind_nests.o ${FLEXPART_SRC}/calcpar.o ${FLEXPART_SRC}/calcpar_nests.o \
richardson.o ew.o getvdep.o calcpv.o obukhov_gfs.o \ ${FLEXPART_SRC}/shift_field.o ${FLEXPART_SRC}/pbl_profile.o ${FLEXPART_SRC}/scalev.o ${FLEXPART_SRC}/obukhov.o \
richardson_gfs.o getvdep_nests.o calcpv_nests.o psim.o psih.o \ ${FLEXPART_SRC}/richardson.o ${FLEXPART_SRC}/ew.o ${FLEXPART_SRC}/getvdep.o ${FLEXPART_SRC}/calcpv.o ${FLEXPART_SRC}/obukhov_gfs.o \
qvsat.o caldate.o getrb.o raerod.o getrc.o partdep.o \ ${FLEXPART_SRC}/richardson_gfs.o ${FLEXPART_SRC}/getvdep_nests.o ${FLEXPART_SRC}/calcpv_nests.o ${FLEXPART_SRC}/psim.o ${FLEXPART_SRC}/psih.o \
verttransform.o verttransform_nests.o readwind_gfs.o \ ${FLEXPART_SRC}/qvsat.o ${FLEXPART_SRC}/caldate.o ${FLEXPART_SRC}/getrb.o ${FLEXPART_SRC}/raerod.o ${FLEXPART_SRC}/getrc.o ${FLEXPART_SRC}/partdep.o \
calcpar_gfs.o verttransform_gfs.o gridcheck_gfs.o ${FLEXPART_SRC}/verttransform.o ${FLEXPART_SRC}/verttransform_nests.o ${FLEXPART_SRC}/readwind_gfs.o \
${FLEXPART_SRC}/calcpar_gfs.o ${FLEXPART_SRC}/verttransform_gfs.o ${FLEXPART_SRC}/gridcheck_gfs.o
VPATH = ${FLEXPART_SRC} VPATH = ${FLEXPART_SRC}
FFLAGS = -mcmodel=medium FFLAGS = -mcmodel=medium
INCLUDES_NETCDF = -I${NETCDFF}/include INCLUDES_NETCDF = -I${NETCDFF}/include
INCLUDES = -I${GRIBAPI}/include ${INCLUDES_NETCDF} INCLUDES = -I${GRIBAPI}/include ${INCLUDES_NETCDF} -I${FLEXPART_SRC}
### NetCDF link flags - use the first one for dynamic libs, the second ### NetCDF link flags - use the first one for dynamic libs, the second
### one for static libs ### one for static libs
#LDFLAGS_NETCDF = -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf LDFLAGS_NETCDF = -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf
LDFLAGS_NETCDF=-static -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf -lnetcdf -L${HDF5}/lib -lhdf5_fortran -lhdf5_hl -lhdf5hl_fortran -lhdf5 -ldl -lz #LDFLAGS_NETCDF=-static -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf -lnetcdf -L${HDF5}/lib -lhdf5_fortran -lhdf5_hl -lhdf5hl_fortran -lhdf5 -ldl -lz
LDFLAGS = -L${GRIBAPI}/lib -lgrib_api_f90 -lgrib_api ${LDFLAGS_NETCDF} -ljasper LDFLAGS = -L${GRIBAPI}/lib -lgrib_api_f90 -lgrib_api ${LDFLAGS_NETCDF} -ljasper -L${FLEXPART_SRC}
#------------ Creating the binary ------------------ #------------ Creating the binary ------------------
...@@ -51,6 +55,13 @@ ${BINARY} : ${BINARY}.o fp2nc4io_mod.o ${FLXPRTOBJS} ${FPMODOBJS} ${OBJS} ...@@ -51,6 +55,13 @@ ${BINARY} : ${BINARY}.o fp2nc4io_mod.o ${FLXPRTOBJS} ${FPMODOBJS} ${OBJS}
${BINARY}.o : ${BINARY}.F90 fp2nc4io_mod.mod ${FPMODOBJS} Makefile ${BINARY}.o : ${BINARY}.F90 fp2nc4io_mod.mod ${FPMODOBJS} Makefile
${FC} -c ${BINARY}.F90 ${FFLAGS} ${INCLUDES} ${FC} -c ${BINARY}.F90 ${FFLAGS} ${INCLUDES}
#----------- NC4 compare ------------------------
${CMP_BINARY} : ${CMP_BINARY}.o
${FC} -o ${CMP_BINARY} ${CMP_BINARY}.o ${LDFLAGS}
${CMP_BINARY}.o : ${CMP_BINARY}.F90 Makefile
${FC} -c ${CMP_BINARY}.F90 ${FFLAGS} ${INCLUDES}
fp2nc4io_mod.mod : ${FPMODOBJS} fp2nc4io_mod.mod : ${FPMODOBJS}
${FC} -c fp2nc4io_mod.F90 ${FFLAGS} ${INCLUDES} ${FC} -c fp2nc4io_mod.F90 ${FFLAGS} ${INCLUDES}
......
...@@ -6,7 +6,7 @@ Contact: Don Morton, don.morton@borealscicomp.com ...@@ -6,7 +6,7 @@ Contact: Don Morton, don.morton@borealscicomp.com
Delia Arnold, deliona.arnold@gmail.com Delia Arnold, deliona.arnold@gmail.com
M. Harustak M. Harustak
Last Update: 29 May 2016 Last Update: 05 July 2017
====================================================================== ======================================================================
1. Introduction 1. Introduction
...@@ -19,6 +19,10 @@ all the operations that FLEXPART normally does in order to prepare the met ...@@ -19,6 +19,10 @@ all the operations that FLEXPART normally does in order to prepare the met
data for computations. It then takes the data which has been stored in data for computations. It then takes the data which has been stored in
FLEXPART global arrays and selectively writes it to a NetCDF4 file. FLEXPART global arrays and selectively writes it to a NetCDF4 file.
Optional parameter allows user to request arrays calculated for specified
location instead of default algorithm selecting lowest left corner of meteorological
domain with a surface pressure above 100000Pa.
This distribution is integrated into the FLEXPART source code, in a grib2nc4 This distribution is integrated into the FLEXPART source code, in a grib2nc4
subdirectory, but it doesn't necessarily have to reside there. The subdirectory, but it doesn't necessarily have to reside there. The
FLEXPART_SRC variable in the Makefile allows for an alternate FP source code FLEXPART_SRC variable in the Makefile allows for an alternate FP source code
...@@ -41,6 +45,14 @@ write it to NC4 file, then read it back and compare with the preprocessed ...@@ -41,6 +45,14 @@ write it to NC4 file, then read it back and compare with the preprocessed
data still in the FLEXPART arrays. It does not test whether preprocessing data still in the FLEXPART arrays. It does not test whether preprocessing
is good, but only whether data was stored correctly in the NC4 files. is good, but only whether data was stored correctly in the NC4 files.
To be able to quickly compare variables in NetCDF files, simple utility
nc4cmp was introduced. It compares specified variable in two NetCDF files
(with dimension up to 3) and it prints average and maximal difference between
the values. It allows for specification of optional tolerance (in %).
Difference is calculated as abs((value1-value2)/value1)*100.
It is calculated for all values of variable multidimensional array.
====================================================================== ======================================================================
2. Quick-Start 2. Quick-Start
====================================================================== ======================================================================
...@@ -64,6 +76,8 @@ distribution, then there should be no need to change FLEXPART_SRC. ...@@ -64,6 +76,8 @@ distribution, then there should be no need to change FLEXPART_SRC.
Simply type "make" (noting that you may see some warnings related to Simply type "make" (noting that you may see some warnings related to
variables for nests) to produce grib2nc4. variables for nests) to produce grib2nc4.
NetCDF compare utility nc4cmp can be built using command "make nc4cmp".
Before trying to run, you might want to try "make test" to make sure that Before trying to run, you might want to try "make test" to make sure that
the simple test works as expected. the simple test works as expected.
...@@ -80,6 +94,11 @@ WARNING - there is currently nothing to protect you if you get these out ...@@ -80,6 +94,11 @@ WARNING - there is currently nothing to protect you if you get these out
of order. It is possible for you to overwrite your met file if you use it of order. It is possible for you to overwrite your met file if you use it
as the second argument. as the second argument.
Usage:
./grib2nc4 <inpath> <outpath> [lat=value lon=value] [optional varnames]
Sample usage: Sample usage:
./grib2nc4 GD15051200 def9_nc1p0.nc4 ./grib2nc4 GD15051200 def9_nc1p0.nc4
...@@ -92,6 +111,29 @@ them to the command. For example: ...@@ -92,6 +111,29 @@ them to the command. For example:
./grib2nc4 GD15051200 def9_nc1p0.nc4 q w ./grib2nc4 GD15051200 def9_nc1p0.nc4 q w
Calculation at user specified location can be requested by adding coord parameter
to the command line:
./grib2nc4 GD15051200 def9_nc1p0.nc4 lon=20 lat=45
To compare two NetCDF files, nc4cmp can be used.
Usage:
./nc4cmp <file1> <file2> <variable> [optional tolerance in %]
Sample usage:
./nc4cmp def9_nc1p0.nc4 def9_nc1p0_reference.nc4 U
Sample usage with tolerance 1%:
./nc4cmp def9_nc1p0.nc4 def9_nc1p0_reference.nc4 U 1.0
This will extract variable U from both NetCDF files, print and compare dimensions,
compare values and print average and maximal difference.
====================================================================== ======================================================================
3. Installation 3. Installation
====================================================================== ======================================================================
......
...@@ -156,6 +156,9 @@ CONTAINS ...@@ -156,6 +156,9 @@ CONTAINS
! Write the height field - variable 'height' is defined in com_mod ! Write the height field - variable 'height' is defined in com_mod
ncfunc_retval = nf90_def_var(ncid, 'height', NF90_DOUBLE, & ncfunc_retval = nf90_def_var(ncid, 'height', NF90_DOUBLE, &
& z_dimid, varid) & z_dimid, varid)
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","height of the FLEXPART model levels")
ncfunc_retval = nf90_put_att(ncid, varid, "units","m a.g.l")
ncfunc_retval = nf90_def_var_deflate(ncid, varid, & ncfunc_retval = nf90_def_var_deflate(ncid, varid, &
& shuffle=0, & & shuffle=0, &
...@@ -168,15 +171,31 @@ CONTAINS ...@@ -168,15 +171,31 @@ CONTAINS
! dx, dy, xlon0, xlat0 are all defined in com_mod ! dx, dy, xlon0, xlat0 are all defined in com_mod
ncfunc_retval = nf90_def_var(ncid, 'dx', NF90_DOUBLE, varid) ncfunc_retval = nf90_def_var(ncid, 'dx', NF90_DOUBLE, varid)
ncfunc_retval = nf90_put_var(ncid, varid, dx) ncfunc_retval = nf90_put_var(ncid, varid, dx)
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","grid distance in x direction")
ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
ncfunc_retval = nf90_def_var(ncid, 'dy', NF90_DOUBLE, varid) ncfunc_retval = nf90_def_var(ncid, 'dy', NF90_DOUBLE, varid)
ncfunc_retval = nf90_put_var(ncid, varid, dy) ncfunc_retval = nf90_put_var(ncid, varid, dy)
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","grid distance in y direction")
ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
ncfunc_retval = nf90_def_var(ncid, 'xlon0', NF90_DOUBLE, varid) ncfunc_retval = nf90_def_var(ncid, 'xlon0', NF90_DOUBLE, varid)
ncfunc_retval = nf90_put_var(ncid, varid, xlon0) ncfunc_retval = nf90_put_var(ncid, varid, xlon0)
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","longitude of the lowest left corner")
ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
ncfunc_retval = nf90_def_var(ncid, 'ylat0', NF90_DOUBLE, varid) ncfunc_retval = nf90_def_var(ncid, 'ylat0', NF90_DOUBLE, varid)
ncfunc_retval = nf90_put_var(ncid, varid, ylat0) ncfunc_retval = nf90_put_var(ncid, varid, ylat0)
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","latitude of the lowest left corner")
ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
! All done, close the NetCDF file ! All done, close the NetCDF file
ncfunc_retval = nf90_close(ncid) ncfunc_retval = nf90_close(ncid)
...@@ -295,18 +314,42 @@ CONTAINS ...@@ -295,18 +314,42 @@ CONTAINS
IF (nc_varname == 'U') THEN IF (nc_varname == 'U') THEN
ncfunc_retval = nf90_put_var(ncid, varid, & ncfunc_retval = nf90_put_var(ncid, varid, &
& uu(0:nx-1, 0:ny-1, 1:nz, 1)) & uu(0:nx-1, 0:ny-1, 1:nz, 1))
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","U component of wind in the X[horizontal] direction")
ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
ELSEIF (nc_varname == 'V') THEN ELSEIF (nc_varname == 'V') THEN
ncfunc_retval = nf90_put_var(ncid, varid, & ncfunc_retval = nf90_put_var(ncid, varid, &
& vv(0:nx-1, 0:ny-1, 1:nz, 1)) & vv(0:nx-1, 0:ny-1, 1:nz, 1))
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","V component of wind in the Y[horizontal] direction")
ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
ELSEIF (nc_varname == 'T') THEN ELSEIF (nc_varname == 'T') THEN
ncfunc_retval = nf90_put_var(ncid, varid, & ncfunc_retval = nf90_put_var(ncid, varid, &
& tt(0:nx-1, 0:ny-1, 1:nz, 1)) & tt(0:nx-1, 0:ny-1, 1:nz, 1))
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","temperature")
ncfunc_retval = nf90_put_att(ncid, varid, "units","k")
ELSEIF (nc_varname == 'W') THEN ELSEIF (nc_varname == 'W') THEN
ncfunc_retval = nf90_put_var(ncid, varid, & ncfunc_retval = nf90_put_var(ncid, varid, &
& ww(0:nx-1, 0:ny-1, 1:nz, 1)) & ww(0:nx-1, 0:ny-1, 1:nz, 1))
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","wind component in the Z[vertical] direction")
ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
ELSEIF (nc_varname == 'Q') THEN ELSEIF (nc_varname == 'Q') THEN
ncfunc_retval = nf90_put_var(ncid, varid, & ncfunc_retval = nf90_put_var(ncid, varid, &
& qv(0:nx-1, 0:ny-1, 1:nz, 1)) & qv(0:nx-1, 0:ny-1, 1:nz, 1))
! attributes
ncfunc_retval = nf90_put_att(ncid, varid, "description","specific humidity")
ncfunc_retval = nf90_put_att(ncid, varid, "units"," ")
ELSE ELSE
PRINT *, PRINT *,
PRINT *, 'fp2nc4io:private_dump_3d_field() bad var: ', nc_varname PRINT *, 'fp2nc4io:private_dump_3d_field() bad var: ', nc_varname
......
...@@ -10,6 +10,11 @@ PROGRAM grib2nc4 ...@@ -10,6 +10,11 @@ PROGRAM grib2nc4
! * ! *
! May 2016 * ! May 2016 *
!************************************************************************* !*************************************************************************
! M. Harustak *
! -) modification to generate the output in single precission *
! -) possibility to add a lat lon selection to obtain the met variables*
! in the vertical levels defined in that location *
!*************************************************************************
USE par_mod USE par_mod
USE com_mod USE com_mod
...@@ -20,20 +25,21 @@ PROGRAM grib2nc4 ...@@ -20,20 +25,21 @@ PROGRAM grib2nc4
IMPLICIT NONE IMPLICIT NONE
LOGICAL :: metfile_exists LOGICAL :: metfile_exists, coordinates_provided, lat_provided, lon_provided
INTEGER :: i, j, k INTEGER :: i, j, k
INTEGER :: num_optional_vars, num_vars INTEGER :: num_optional_vars, num_vars
INTEGER, PARAMETER :: DEFLATE_LEVEL = 2 ! Compression level (0-9) INTEGER, PARAMETER :: DEFLATE_LEVEL = 2 ! Compression level (0-9)
CHARACTER(LEN=512) :: met_filepath, netcdf4_filepath CHARACTER(LEN=512) :: met_filepath, netcdf4_filepath, param_str, coord_name_str, coord_val_str
CHARACTER, DIMENSION(:), ALLOCATABLE :: vars_list ! list of variables CHARACTER, DIMENSION(:), ALLOCATABLE :: vars_list ! list of variables
INTEGER :: coordX, coordY, stat
REAL :: coord_lat, coord_lon
INTEGER :: metdata_format = UNKNOWN_METDATA ! From FP par_mod INTEGER :: metdata_format = UNKNOWN_METDATA ! From FP par_mod
!-------------------------------------------------------- !--------------------------------------------------------
! Read in mandatory arguments ! Read in mandatory arguments
IF (IARGC() < 2) THEN IF (IARGC() < 2) THEN
PRINT *, 'Usage: grib2netcdf4 <inpath> <outpath> [optional varnames]' PRINT *, 'Usage: grib2netcdf4 <inpath> <outpath> [lon=X lat=Y] [optional varnames]'
STOP STOP
ELSE ELSE
CALL GETARG(1, met_filepath) CALL GETARG(1, met_filepath)
...@@ -47,23 +53,47 @@ PROGRAM grib2nc4 ...@@ -47,23 +53,47 @@ PROGRAM grib2nc4
! First, get the number of optional args and allocate vars_list, ! First, get the number of optional args and allocate vars_list,
! and fill the first three elements ! and fill the first three elements
IF (IARGC() > 2) THEN coordinates_provided = .FALSE.
num_optional_vars = IARGC() - 2 lat_provided = .FALSE.
ELSE lon_provided = .FALSE.
num_optional_vars = 0 ALLOCATE( vars_list(IARGC()+3),stat=stat )
ENDIF
num_vars = num_optional_vars + 3
ALLOCATE( vars_list(num_vars) )
vars_list(1) = 'u' vars_list(1) = 'u'
vars_list(2) = 'v' vars_list(2) = 'v'
vars_list(3) = 't' vars_list(3) = 't'
! Read in optional variable arguments, starting at element 4 of vars_list num_vars = 3
IF (IARGC() > 2) THEN DO i=3,IARGC()
DO i=1,num_optional_vars CALL GETARG(i,param_str)
CALL GETARG( i+2, vars_list(i+3) ) param_str = TRIM(param_str)
ENDDO j = SCAN(param_str,"=")
if (j>1) then
coord_name_str=param_str(1:j-1)
coord_val_str=param_str(j+1:)
IF ( coord_name_str == "lat" .or. coord_name_str == "LAT" ) THEN
read(coord_val_str,*,iostat=stat) coord_lat
if ( stat == 0 ) then
lat_provided = .TRUE.
else
print *, "Incorrect coordinates: ", coord_val_str
stop
endif
ELSE IF ( coord_name_str == "lon" .or. coord_name_str == "LON" ) THEN
read(coord_val_str,*,iostat=stat) coord_lon
if ( stat == 0 ) then
lon_provided = .TRUE.
else
print *, "Incorrect coordinates: ", coord_val_str
stop
endif
ENDIF
else
num_vars = num_vars + 1
vars_list(num_vars) = param_str
endif
ENDDO
IF (lat_provided .AND. lon_provided) THEN
coordinates_provided = .TRUE.
ENDIF ENDIF
! Before proceeding, let's make sure the vars_list is good - otherwise, ! Before proceeding, let's make sure the vars_list is good - otherwise,
...@@ -125,8 +155,18 @@ PROGRAM grib2nc4 ...@@ -125,8 +155,18 @@ PROGRAM grib2nc4
! to do so ! to do so
!CALL gridcheck_nests !CALL gridcheck_nests
! If the coordinates are provided, then we need to obtain
! the corresponding grid indexes since this is what verttransform needs
if ( coordinates_provided ) then
coordX = (coord_lon-xlon0)/dx
coordY = (coord_lat-ylat0)/dy
print *, "Coordinates: "
print *, "lon: ", coord_lon, ", lat: ",coord_lat
print *, "x: ", coordX, ", y: ", coordY
endif
PRINT *, 'Calling processmetfields()...' PRINT *, 'Calling processmetfields()...'
call processmetfields( 1, metdata_format) call processmetfields( 1, metdata_format, coordinates_provided, coordX, coordY)
PRINT *, 'Calling fp2nc4io_dump()...' PRINT *, 'Calling fp2nc4io_dump()...'
call fp2nc4io_dump( netcdf4_filepath, num_vars, vars_list, DEFLATE_LEVEL) call fp2nc4io_dump( netcdf4_filepath, num_vars, vars_list, DEFLATE_LEVEL)
......
PROGRAM nc4cmp
!*************************************************************************
! This program compares variable in two NetCDF files *
! (with dimension up to 3) and it prints average and maximal difference *
! between the values. It allows for specification of optional tolerance *
! *
! Usage: *
! *
! ./nc4cmp <file1> <file2> <variable> [optional tolerance in %] *
! *
! M. Harustak *
! *
! March 2017 *
!*************************************************************************
USE netcdf
IMPLICIT NONE
LOGICAL :: file1_exists, file2_exists
CHARACTER(LEN=512) :: file1, file2, precision_str, variable
CHARACTER(LEN=NF90_MAX_NAME) :: returned_name
REAL :: tolerance
INTEGER :: ncid, varid, retval, xtype, numDims1, numDims2, i, j, k, stat
INTEGER, dimension(NF90_MAX_VAR_DIMS) :: varDimIDs1, varDimIDs2, dims1, dims2
double precision :: variable_array_0d_1, variable_array_0d_2
double precision, allocatable, dimension(:) :: variable_array_1d_1, variable_array_1d_2
double precision, allocatable, dimension(:,:) :: variable_array_2d_1, variable_array_2d_2
double precision, allocatable, dimension(:,:,:) :: variable_array_3d_1, variable_array_3d_2
double precision :: v1, v2
double precision :: sum_diff, max_diff
integer :: count_diff
logical :: different = .FALSE.
! Read in mandatory arguments
IF (IARGC() < 3) THEN
PRINT *, 'Usage: nc4cmp <file1> <file2> <variable> [optional tolerance in %]'
STOP
ELSE
CALL GETARG(1, file1)
CALL GETARG(2, file2)
CALL GETARG(3, variable)
ENDIF
IF (IARGC() > 3) THEN
CALL GETARG( 4, precision_str)
READ(precision_str, '(f10.5)') tolerance
ELSE
tolerance = 0
ENDIF
print *,"Tolerance ",tolerance,"%"
INQUIRE( FILE=file1, EXIST=file1_exists )
IF ( file1_exists ) THEN
PRINT *, 'File 1: ', TRIM(file1)
ELSE
PRINT *, 'Unable to find file: ', TRIM(file1)
STOP
ENDIF
INQUIRE( FILE=file2, EXIST=file2_exists )
IF ( file2_exists ) THEN
PRINT *, 'File 2: ', TRIM(file2)
ELSE
PRINT *, 'Unable to find file: ', TRIM(file2)
STOP
ENDIF
print *, "Opening ",TRIM(file1)," for reading"
retval = nf90_open(file1, NF90_NOWRITE, ncid)
retval = nf90_inq_varid(ncid, TRIM(variable) , varid)
retval = nf90_inquire_variable(ncid, varid, ndims = numDims1, dimids = varDimIDs1)
print *, "Reading variable ",TRIM(variable), ", no of dimensions: ",numDims1
do i=1,numDims1
retval = nf90_inquire_dimension(ncid, varDimIDs1(i),len=j)
dims1(i) = j
print *," - dimension ",i,": ", j
enddo
if ( numDims1 == 0) then
retval = nf90_get_var(ncid, varid, variable_array_0d_1)
else if ( numDims1 == 1) then
allocate(variable_array_1d_1(dims1(1)),stat=stat)
retval = nf90_get_var(ncid, varid, variable_array_1d_1)
else if ( numDims1 == 2) then
allocate(variable_array_2d_1(dims1(1),dims1(2)),stat=stat)
retval = nf90_get_var(ncid, varid, variable_array_2d_1)
else if ( numDims1 == 3) then
allocate(variable_array_3d_1(dims1(1),dims1(2),dims1(3)),stat=stat)
retval = nf90_get_var(ncid, varid, variable_array_3d_1)
endif
retval = nf90_close(ncid)
print *, "...done"
print *, "Opening ",TRIM(file2)," for reading"
retval = nf90_open(file2, NF90_NOWRITE, ncid)
retval = nf90_inq_varid(ncid, TRIM(variable) , varid)
retval = nf90_inquire_variable(ncid, varid, ndims = numDims2, dimids = varDimIDs2)
print *, "Reading variable ",TRIM(variable), ", no of dimensions: ",numDims2
do i=1,numDims2
retval = nf90_inquire_dimension(ncid, varDimIDs2(i),len=j)
dims2(i) = j
print *," - dimension ",i,": ", j
enddo
if ( numDims2 == 0) then
retval = nf90_get_var(ncid, varid, variable_array_0d_2)
else if ( numDims2 == 1) then
allocate(variable_array_1d_2(dims2(1)),stat=stat)
retval = nf90_get_var(ncid, varid, variable_array_1d_2)
else if ( numDims2 == 2) then
allocate(variable_array_2d_2(dims2(1),dims2(2)),stat=stat)
retval = nf90_get_var(ncid, varid, variable_array_2d_2)
else if ( numDims2 == 3) then
allocate(variable_array_3d_2(dims2(1),dims2(2),dims2(3)),stat=stat)
retval = nf90_get_var(ncid, varid, variable_array_3d_2)
endif
retval = nf90_close(ncid)
print *, "...done"
if (numDims1 /= numDims2) THEN
print *, "Number of dimensions differs..."
STOP
endif
do i=1,numDims1
if (dims1(i) /= dims2(i)) then