From 50dc6c02b616bb253b0ceae0762fda1844e2ea70 Mon Sep 17 00:00:00 2001 From: Ron Date: Tue, 16 Apr 2024 17:29:20 -0700 Subject: [PATCH] Refactored dependent codes into a single code. Modified makefile template accordingly. Updated build scripts. --- LICENSE | 2 +- ...pu_mpi+multithread_gcc10+11_ubuntu20.04.sh | 59 - ...ld_cpu_mpi+multithread_gcc_ubuntu20.04.sh} | 0 ...cpu_mpi+multithread_nvidia_ubuntu20.04.sh} | 0 ...build_cpu_mpi-only_gcc10+11_ubuntu20.04.sh | 59 - ... => build_cpu_mpi-only_gcc_ubuntu20.04.sh} | 0 ... build_cpu_mpi-only_nvidia_ubuntu20.04.sh} | 0 ....04.sh => build_gpu_nvidia_ubuntu20.04.sh} | 2 +- src/Makefile.template | 20 +- src/number_types.f | 22 - src/pot3d.F | 23 +- src/psi_io.f90 | 1376 ++++++++++ src/stdpar/pot3d.F | 23 +- src/zm_parse.f | 1946 ------------- src/zm_parse_modules.f | 164 -- src/zm_sds.f | 2440 ----------------- src/zm_sds_modules.f | 150 - 17 files changed, 1424 insertions(+), 4862 deletions(-) delete mode 100755 build_examples/build_cpu_mpi+multithread_gcc10+11_ubuntu20.04.sh rename build_examples/{build_cpu_mpi+multithread_gcc9_ubuntu20.04.sh => build_cpu_mpi+multithread_gcc_ubuntu20.04.sh} (100%) rename build_examples/{build_cpu_mpi+multithread_nv22.3_ubuntu20.04.sh => build_cpu_mpi+multithread_nvidia_ubuntu20.04.sh} (100%) delete mode 100755 build_examples/build_cpu_mpi-only_gcc10+11_ubuntu20.04.sh rename build_examples/{build_cpu_mpi-only_gcc9_ubuntu20.04.sh => build_cpu_mpi-only_gcc_ubuntu20.04.sh} (100%) rename build_examples/{build_cpu_mpi-only_nv22.3_ubuntu20.04.sh => build_cpu_mpi-only_nvidia_ubuntu20.04.sh} (100%) rename build_examples/{build_gpu_nv22.3_ubuntu20.04.sh => build_gpu_nvidia_ubuntu20.04.sh} (96%) delete mode 100644 src/number_types.f create mode 100644 src/psi_io.f90 delete mode 100644 src/zm_parse.f delete mode 100644 src/zm_parse_modules.f delete mode 100644 src/zm_sds.f delete mode 100644 src/zm_sds_modules.f diff --git a/LICENSE b/LICENSE index 261eeb9..e30e94e 100644 --- a/LICENSE +++ b/LICENSE @@ -186,7 +186,7 @@ same "printed page" as the copyright notice for easier identification within third-party archives. - Copyright [yyyy] [name of copyright owner] + Copyright 2024 Predictive Science Inc. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. diff --git a/build_examples/build_cpu_mpi+multithread_gcc10+11_ubuntu20.04.sh b/build_examples/build_cpu_mpi+multithread_gcc10+11_ubuntu20.04.sh deleted file mode 100755 index be6b09c..0000000 --- a/build_examples/build_cpu_mpi+multithread_gcc10+11_ubuntu20.04.sh +++ /dev/null @@ -1,59 +0,0 @@ -#!/bin/bash -################################################################### -# This build assumes that you have an "mpif90" in your PATH -# set up to use your chosen MPI library and compiler. -################################################################## -################################################################# -# Please set the location of HDF5 include/library files and -# the linker flags to match your installed version. -# -# Note! The HDF5 library needs to have been compiled with -# the same compiler being used here and is loaded in the run-time -# environment (e.g. LD_LIBRARY_PATH). -################################################################# - -# Location of local hdf5 installed with same compiler being used for POT3D: -HDF5_INCLUDE_DIR="/usr/include/hdf5/serial" -HDF5_LIB_DIR="/usr/lib/x86_64-linux-gnu" -# Fortran HDF5 library flags (these can be version dependent): -HDF5_LIB_FLAGS="-lhdf5_serial_fortran -lhdf5_serialhl_fortran -lhdf5_serial -lhdf5_serial_hl" - -########################################################################### -# Please set the compile flags based on your compiler and hardware setup. -########################################################################### - -FFLAGS="-O3 -march=native -ftree-parallelize-loops=${OMP_NUM_THREADS} -fallow-argument-mismatch" - -########################################################################### -# If using NV HPC SDK for GPUs, with CUDA version >= 11.3, you can set -# POT3D_CUSPARSE to "1" to link the cuSparse library, allowing you to set -# 'ifprec=2' in 'pot3d.dat' to yield ~2x speed improvement! -# Warning! Using ifprec=2 takes much more GPU memory than ifprec=1. -# You must also set CCFLAGS to use OpenACC in the C code. -########################################################################### - -POT3D_CUSPARSE=0 -CCFLAGS="-O3" - -########################################################################### -########################################################################### -########################################################################### - -POT3D_HOME=$PWD - -cd ${POT3D_HOME}/src -echo "Making copy of Makefile..." -cp Makefile.template Makefile -echo "Modifying Makefile to chosen flags..." -sed -i "s##${FFLAGS}#g" Makefile -sed -i "s##${CCFLAGS}#g" Makefile -sed -i "s##${POT3D_CUSPARSE}#g" Makefile -sed -i "s##${HDF5_INCLUDE_DIR}#g" Makefile -sed -i "s##${HDF5_LIB_DIR}#g" Makefile -sed -i "s##${HDF5_LIB_FLAGS}#g" Makefile -echo "Building POT3D...." -make 1>build.log 2>build.err -echo "Copying POT3D executable from SRC to BIN..." -cp ${POT3D_HOME}/src/pot3d ${POT3D_HOME}/bin/pot3d -echo "Done!" - diff --git a/build_examples/build_cpu_mpi+multithread_gcc9_ubuntu20.04.sh b/build_examples/build_cpu_mpi+multithread_gcc_ubuntu20.04.sh similarity index 100% rename from build_examples/build_cpu_mpi+multithread_gcc9_ubuntu20.04.sh rename to build_examples/build_cpu_mpi+multithread_gcc_ubuntu20.04.sh diff --git a/build_examples/build_cpu_mpi+multithread_nv22.3_ubuntu20.04.sh b/build_examples/build_cpu_mpi+multithread_nvidia_ubuntu20.04.sh similarity index 100% rename from build_examples/build_cpu_mpi+multithread_nv22.3_ubuntu20.04.sh rename to build_examples/build_cpu_mpi+multithread_nvidia_ubuntu20.04.sh diff --git a/build_examples/build_cpu_mpi-only_gcc10+11_ubuntu20.04.sh b/build_examples/build_cpu_mpi-only_gcc10+11_ubuntu20.04.sh deleted file mode 100755 index 96e370b..0000000 --- a/build_examples/build_cpu_mpi-only_gcc10+11_ubuntu20.04.sh +++ /dev/null @@ -1,59 +0,0 @@ -#!/bin/bash -################################################################### -# This build assumes that you have an "mpif90" in your PATH -# set up to use your chosen MPI library and compiler. -################################################################## -################################################################# -# Please set the location of HDF5 include/library files and -# the linker flags to match your installed version. -# -# Note! The HDF5 library needs to have been compiled with -# the same compiler being used here and is loaded in the run-time -# environment (e.g. LD_LIBRARY_PATH). -################################################################# - -# Location of local hdf5 installed with same compiler being used for POT3D: -HDF5_INCLUDE_DIR="/usr/include/hdf5/serial" -HDF5_LIB_DIR="/usr/lib/x86_64-linux-gnu" -# Fortran HDF5 library flags (these can be version dependent): -HDF5_LIB_FLAGS="-lhdf5_serial_fortran -lhdf5_serialhl_fortran -lhdf5_serial -lhdf5_serial_hl" - -########################################################################### -# Please set the compile flags based on your compiler and hardware setup. -########################################################################### - -FFLAGS="-O3 -march=native -fallow-argument-mismatch" - -########################################################################### -# If using NV HPC SDK for GPUs, with CUDA version >= 11.3, you can set -# POT3D_CUSPARSE to "1" to link the cuSparse library, allowing you to set -# 'ifprec=2' in 'pot3d.dat' to yield ~2x speed improvement! -# Warning! Using ifprec=2 takes much more GPU memory than ifprec=1. -# You must also set CCFLAGS to use OpenACC in the C code. -########################################################################### - -POT3D_CUSPARSE=0 -CCFLAGS="-O3" - -########################################################################### -########################################################################### -########################################################################### - -POT3D_HOME=$PWD - -cd ${POT3D_HOME}/src -echo "Making copy of Makefile..." -cp Makefile.template Makefile -echo "Modifying Makefile to chosen flags..." -sed -i "s##${FFLAGS}#g" Makefile -sed -i "s##${CCFLAGS}#g" Makefile -sed -i "s##${POT3D_CUSPARSE}#g" Makefile -sed -i "s##${HDF5_INCLUDE_DIR}#g" Makefile -sed -i "s##${HDF5_LIB_DIR}#g" Makefile -sed -i "s##${HDF5_LIB_FLAGS}#g" Makefile -echo "Building POT3D...." -make 1>build.log 2>build.err -echo "Copying POT3D executable from SRC to BIN..." -cp ${POT3D_HOME}/src/pot3d ${POT3D_HOME}/bin/pot3d -echo "Done!" - diff --git a/build_examples/build_cpu_mpi-only_gcc9_ubuntu20.04.sh b/build_examples/build_cpu_mpi-only_gcc_ubuntu20.04.sh similarity index 100% rename from build_examples/build_cpu_mpi-only_gcc9_ubuntu20.04.sh rename to build_examples/build_cpu_mpi-only_gcc_ubuntu20.04.sh diff --git a/build_examples/build_cpu_mpi-only_nv22.3_ubuntu20.04.sh b/build_examples/build_cpu_mpi-only_nvidia_ubuntu20.04.sh similarity index 100% rename from build_examples/build_cpu_mpi-only_nv22.3_ubuntu20.04.sh rename to build_examples/build_cpu_mpi-only_nvidia_ubuntu20.04.sh diff --git a/build_examples/build_gpu_nv22.3_ubuntu20.04.sh b/build_examples/build_gpu_nvidia_ubuntu20.04.sh similarity index 96% rename from build_examples/build_gpu_nv22.3_ubuntu20.04.sh rename to build_examples/build_gpu_nvidia_ubuntu20.04.sh index 6deac8d..bb637e2 100755 --- a/build_examples/build_gpu_nv22.3_ubuntu20.04.sh +++ b/build_examples/build_gpu_nvidia_ubuntu20.04.sh @@ -23,7 +23,7 @@ HDF5_LIB_FLAGS="-lhdf5_fortran -lhdf5hl_fortran -lhdf5 -lhdf5_hl" # Please set the compile flags based on your compiler and hardware setup. ########################################################################### -FFLAGS="-O3 -march=native -stdpar=gpu -acc=gpu -gpu=nomanaged -Minfo=accel" +FFLAGS="-O3 -march=native -stdpar=gpu -acc=gpu -gpu=nomanaged,nounified -Minfo=accel" ########################################################################### # If using NV HPC SDK for GPUs, with CUDA version >= 11.3, you can set diff --git a/src/Makefile.template b/src/Makefile.template index 639e71b..b5f2a0b 100644 --- a/src/Makefile.template +++ b/src/Makefile.template @@ -14,11 +14,7 @@ endif FFLAGS = $(PARADD) -I -OBJS0 = number_types.o \ - zm_parse_modules.o \ - zm_parse.o \ - zm_sds_modules.o \ - zm_sds.o +OBJS0 = psi_io.o ifeq ($(POT3D_CUSPARSE),1) OBJS = $(OBJS0) lusol_cusparse.o pot3d_cpp.o @@ -39,19 +35,7 @@ clean: pot3d_cpp.f: pot3d.F $(FC) -E -cpp $(IF_DEF) > pot3d_cpp.f $< -number_types.o: number_types.f - $(FC) -c $(FFLAGS) $< - -zm_parse_modules.o: zm_parse_modules.f - $(FC) -c $(FFLAGS) $< - -zm_parse.o: zm_parse.f zm_parse_modules.f number_types.f - $(FC) -c $(FFLAGS) $< - -zm_sds_modules.o: zm_sds_modules.f - $(FC) -c $(FFLAGS) $< - -zm_sds.o: zm_sds.f zm_sds_modules.f number_types.f +psi_io.o: psi_io.f90 $(FC) -c $(FFLAGS) $< lusol_cusparse.o: lusol_cusparse.c diff --git a/src/number_types.f b/src/number_types.f deleted file mode 100644 index 02746e6..0000000 --- a/src/number_types.f +++ /dev/null @@ -1,22 +0,0 @@ -c####################################################################### - module number_types -c -c----------------------------------------------------------------------- -c ****** Basic number types. -c ****** This module is used to set the default precision for REALs. -c----------------------------------------------------------------------- -c - use iso_fortran_env -c -c----------------------------------------------------------------------- -c - implicit none -c - integer, parameter :: KIND_REAL_4=REAL32 - integer, parameter :: KIND_REAL_8=REAL64 - integer, parameter :: KIND_REAL_16=max(REAL128,REAL64) -c - integer, parameter :: r_typ=KIND_REAL_8 -c - end module -c####################################################################### diff --git a/src/pot3d.F b/src/pot3d.F index 3860c4e..68ef7a5 100644 --- a/src/pot3d.F +++ b/src/pot3d.F @@ -56,6 +56,27 @@ module ident character(*), parameter :: update='10/24/2022' c end module +c####################################################################### + module number_types +c +c----------------------------------------------------------------------- +c ****** Basic number types. +c ****** This module is used to set the default precision for REALs. +c----------------------------------------------------------------------- +c + use iso_fortran_env +c +c----------------------------------------------------------------------- +c + implicit none +c + integer, parameter :: KIND_REAL_4=REAL32 + integer, parameter :: KIND_REAL_8=REAL64 + integer, parameter :: KIND_REAL_16=max(REAL128,REAL64) +c + integer, parameter :: r_typ=KIND_REAL_8 +c + end module c####################################################################### module number_types_pc c @@ -6414,8 +6435,8 @@ subroutine readbr (fname,br0_g,ierr) use number_types use global_dims use global_mesh - use rdhdf_2d_interface use vars + use rdhdf_2d_interface c c----------------------------------------------------------------------- c diff --git a/src/psi_io.f90 b/src/psi_io.f90 new file mode 100644 index 0000000..21c8eb0 --- /dev/null +++ b/src/psi_io.f90 @@ -0,0 +1,1376 @@ +!####################################################################### +! ****** PSI I/O: PSI I/O Tools. +! +! Authors: Predictive Science Inc. +! +! Predictive Science Inc. +! www.predsci.com +! San Diego, California, USA 92121 +!####################################################################### +! Copyright 2024 Predictive Science Inc. +! +! Licensed under the Apache License, Version 2.0 (the "License"); +! you may not use this file except in compliance with the License. +! You may obtain a copy of the License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software +! distributed under the License is distributed on an "AS IS" BASIS, +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or +! implied. +! See the License for the specific language governing permissions and +! limitations under the License. +!####################################################################### +module rp1d_def +! +!----------------------------------------------------------------------- +! ****** Define a structure to hold a REAL 1D pointer. +!----------------------------------------------------------------------- +! + use iso_fortran_env +! + implicit none +! + type :: rp1d + real(REAL64), dimension(:), pointer, contiguous :: f + end type +! +end module +!#1###################################################################### +module sds_def +! +!----------------------------------------------------------------------- +! ****** Definition of the IO data structure. +!----------------------------------------------------------------------- +! + use iso_fortran_env + use rp1d_def +! + implicit none +! + integer, parameter, private :: mxdim = 3 +! + type :: sds + integer :: ndim + integer, dimension(mxdim) :: dims + logical :: scale + logical :: hdf32 + type(rp1d), dimension(mxdim) :: scales + real(REAL64), dimension(:,:,:), pointer, contiguous :: f + end type +! +end module +!####################################################################### +module rdhdf_1d_interface + interface + subroutine rdhdf_1d (fname,scale,nx,f,x,ierr) + use iso_fortran_env + implicit none + character(*) :: fname + logical :: scale + integer :: nx + real(REAL64), dimension(:), pointer :: f + real(REAL64), dimension(:), pointer :: x + integer :: ierr + intent(in) :: fname + intent(out) :: scale,nx,ierr + end subroutine + end interface +end module +!####################################################################### +module rdhdf_2d_interface + interface + subroutine rdhdf_2d (fname,scale,nx,ny,f,x,y,ierr) + use iso_fortran_env + implicit none + character(*) :: fname + logical :: scale + integer :: nx,ny + real(REAL64), dimension(:,:), pointer :: f + real(REAL64), dimension(:), pointer :: x,y + integer :: ierr + intent(in) :: fname + intent(out) :: scale,nx,ny,ierr + end subroutine + end interface +end module +!####################################################################### +module rdhdf_3d_interface + interface + subroutine rdhdf_3d (fname,scale,nx,ny,nz,f,x,y,z,ierr) + use iso_fortran_env + implicit none + character(*) :: fname + logical :: scale + integer :: nx,ny,nz + real(REAL64), dimension(:,:,:), pointer :: f + real(REAL64), dimension(:), pointer :: x,y,z + integer :: ierr + intent(in) :: fname + intent(out) :: scale,nx,ny,nz,ierr + end subroutine + end interface +end module +!####################################################################### +subroutine ffopen (iun,fname,mode,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Open file FNAME and link it to unit IUN. +! +! ****** If there is an error, this routine returns IERR.ne.0. +! +!----------------------------------------------------------------------- +! +! ****** When MODE='r', the file must exist. +! ****** When MODE='w', the file is created. +! ****** When MODE='rw', the file must exist, but can be overwritten. +! ****** When MODE='a', the file is created if it does not exist, +! ****** otherwise, it is appended. +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + integer :: iun + character(*) :: fname + character(*) :: mode + integer :: ierr + logical :: ex +! +!----------------------------------------------------------------------- +! + ierr=0 +! + if (mode.eq.'r') then + open (iun,file=fname,form="FORMATTED",status='old',err=900) + else if (mode.eq.'rw') then + open (iun,file=fname,form="FORMATTED",status='replace',err=900) + else if (mode.eq.'w') then + open (iun,file=fname,form="FORMATTED",status='new',err=900) + elseif (mode.eq.'a') then + inquire(file=fname, exist=ex) + if (ex) then + open (iun,file=fname,form="FORMATTED",position='append',err=900) + else + open (iun,file=fname,form="FORMATTED",status='new',err=900) + end if + else + write (*,*) + write (*,*) '### ERROR in FFOPEN:' + write (*,*) '### Invalid MODE requested.' + write (*,*) 'MODE = ',mode + write (*,*) 'File name: ',trim(fname) + ierr=2 + return + end if +! + return +! + 900 continue +! + write (*,*) + write (*,*) '### ERROR in FFOPEN:' + write (*,*) '### Error while opening the requested file.' + write (*,*) 'File name: ',trim(fname) + write (*,*) 'MODE = ',mode + ierr=1 +! +end subroutine +!####################################################################### +subroutine rdhdf (fname,s,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Read a 1D, 2D, or 3D scientific data set from an HDF file. +! ****** This routine uses the new SD API instead of the +! ****** outdated DFSD API. +! +!----------------------------------------------------------------------- +! +! ****** This routine allocates the required memory and returns +! ****** pointers to the data and scale arrays. +! +!----------------------------------------------------------------------- +! +! ****** Input arguments: +! +! FNAME : [character(*)] +! HDF data file name to read from. +! +! ****** Output arguments: +! +! S : [structure of type SDS] +! A structure that holds the field, its +! dimensions, and the scales, with the +! components described below. +! +! IERR : [integer] +! IERR=0 is returned if the data set was read +! successfully. Otherwise, IERR is set to a +! nonzero value. +! +! ****** Components of structure S: +! +! NDIM : [integer] +! Number of dimensions found in the data set. +! +! DIMS : [integer, dimension(3)] +! Number of points in the data set dimensions. +! For a 1D data set, DIMS(2)=DIMS(3)=1. +! For a 2D data set, DIMS(3)=1. +! +! SCALE : [logical] +! Flag to indicate the presence of scales (axes) +! in the data set. SCALE=.false. means that scales +! were not found; SCALE=.true. means that scales +! were found. +! +! HDF32 : [logical] +! Flag to indicate the precision of the data set +! read in. HDF32=.true. means that the data is +! 32-bit; HDF32=.false. means that the data is +! 64-bit. +! +! SCALES : [structure of type RP1D, dimension(3)] +! This array holds the pointers to the scales +! when SCALE=.true., and is undefined otherwise. +! +! F : [real, pointer to a rank-3 array] +! This array holds the data set values. +! +! ****** The storage for the arrays pointed to by F, and the +! ****** scales (if present) in structure SCALES, is allocated by +! ****** this routine. +! +!----------------------------------------------------------------------- +! + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + character, dimension(256) :: sds_name, dim_name + type(sds) :: s + integer :: ierr,i + intent(in) :: fname + intent(out) :: s,ierr +! +!----------------------------------------------------------------------- +! + ierr=0 +! +!----------------------------------------------------------------------- +! +! ****** Read hdf5 file if fname ends in '.h5'. +! + i=index(fname,'.h'); + if (fname(i+1:i+2).eq.'h5') then + call rdh5 (fname,s,ierr) + return + else + print*,"HDF4 has been disabled" + ierr=-1 + end if +! +end subroutine +!####################################################################### +subroutine rdh5 (fname,s,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Read a 1D, 2D, or 3D scientific data set from an HDF5 file. +! ****** The HDF5 file is currently assumed to contain only one +! ****** dataset (1D,2D,or 3D), with or without scales, in group "/", +! ****** and has no other data members. +! +!----------------------------------------------------------------------- +! +! ****** This routine allocates the required memory and returns +! ****** pointers to the data and scale arrays. +! +!----------------------------------------------------------------------- +! +! ****** Input arguments: +! +! FNAME : [character(*)] +! File name to read from. +! +! ****** Output arguments: +! +! S : [structure of type DS] +! A structure that holds the field, its +! dimensions, and the scales, with the +! components described below. +! +! IERR : [integer] +! IERR=0 is returned if the data set was read +! successfully. Otherwise, IERR is set to a +! nonzero value. +! +! ****** Components of structure S: +! +! NDIM : [integer] +! Number of dimensions found in the data set. +! +! DIMS : [integer, dimension(3)] +! Number of points in the data set dimensions. +! For a 1D data set, DIMS(2)=DIMS(3)=1. +! For a 2D data set, DIMS(3)=1. +! +! SCALE : [logical] +! Flag to indicate the presence of scales (axes) +! in the data set. SCALE=.false. means that scales +! were not found; SCALE=.true. means that scales +! were found. +! +! HDF32 : [logical] +! Flag to indicate the precision of the data set +! read in. HDF32=.true. means that the data is +! 32-bit; HDF32=.false. means that the data is +! 64-bit. +! +! SCALES : [structure of type RP1D, dimension(3)] +! This array holds the pointers to the scales +! when SCALE=.true., and is undefined otherwise. +! +! F : [real, pointer to a rank-3 array] +! This array holds the data set values. +! +! ****** The storage for the arrays pointed to by F, and the +! ****** scales (if present) in structure SCALES, is allocated by +! ****** this routine. +! +!----------------------------------------------------------------------- +! + use iso_fortran_env + use sds_def + use hdf5 + use h5ds +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + type(sds) :: s + character(*) :: fname +! +!----------------------------------------------------------------------- +! + integer :: ierr +! +!----------------------------------------------------------------------- +! + integer :: i,j,k,n_members,obj_type +! + integer(HID_T) :: file_id ! File identifier + integer(HID_T) :: dset_id ! Dataset identifier + integer(HID_T) :: dspace_id ! Dataspace identifier + integer(HID_T) :: dim_id ! Dimension identifiers + integer(HID_T) :: datatype_id ! Datatype identifiers +! + integer(SIZE_T) :: prec +! + integer(HSIZE_T),dimension(:), allocatable :: s_dims,maxpts + integer(HSIZE_T),dimension(1) :: s_dims_i +! + real(REAL32), dimension(:,:,:), allocatable :: f4 + real(REAL32), dimension(:), allocatable :: f4dim + real(REAL64), dimension(:,:,:), allocatable :: f8 + real(REAL64), dimension(:), allocatable :: f8dim +! + character(512) :: obj_name + character(4), parameter :: cname='RDH5' +! + logical :: is_scale +! +!----------------------------------------------------------------------- +! +! ****** Initialize dimension count and arrays. +! + s%ndim = 0 + s%dims(:) = 1 +! +! ****** Initialize hdf5 interface. +! + call h5open_f (ierr) +! +! ****** Open hdf5 file. +! + call h5Fopen_f (trim(fname),H5F_ACC_RDONLY_F,file_id,ierr) +! +! ****** Get information about the hdf5 file. +! + call h5Gn_members_f (file_id,"/",n_members,ierr) +! +! ****** Make sure there is (at maximum) one 3D dataset with scales. +! + if (n_members.eq.0.or.n_members.gt.4) then + write (*,*) + write (*,*) '### ERROR in ',cname,':' + write (*,*) '### Input file contains too few/many datasets.' + write (*,*) 'File name: ',trim(fname) + return + endif +! +! ****** Assume the Dataset is in index 0 and get its name. +! + call h5Gget_obj_info_idx_f (file_id,"/",0,obj_name,obj_type,ierr) +! +! ****** Open Dataset. +! + call h5Dopen_f (file_id,trim(obj_name),dset_id,ierr) +! +! ****** Make sure the Dataset is not a scale. +! + call h5DSis_scale_f(dset_id,is_scale,ierr) + if (is_scale) then + write (*,*) + write (*,*) '### ERROR in ',cname,':' + write (*,*) '### Input file Dataset at index 0 is a scale.' + write (*,*) 'File name: ',trim(fname) + return + endif +! +! ****** Get dimensions (need s_dims array for format requirements). +! + call h5Dget_space_f (dset_id,dspace_id,ierr) + call h5Sget_simple_extent_ndims_f (dspace_id,s%ndim,ierr) +! + allocate(s_dims(s%ndim)) +! + allocate(maxpts(s%ndim)) + call h5Sget_simple_extent_dims_f (dspace_id,s_dims,maxpts,ierr) + deallocate(maxpts) +! + do j=1,s%ndim + s%dims(j) = INT(s_dims(j)) + end do +! +! ****** Get the floating-point precision of the data and set flag. +! + call h5Dget_type_f (dset_id,datatype_id,ierr) + call h5Tget_precision_f (datatype_id,prec,ierr) +! + if (prec.eq.32) then + s%hdf32=.true. + elseif (prec.eq.64) then + s%hdf32=.false. + end if +! +! ****** Allocate the memory for the Dataset array in s. +! + allocate (s%f(s%dims(1),s%dims(2),s%dims(3))) +! +! ****** Need to read the file in its own datatype, and then convert +! ****** to datatype of s%f. +! + if (s%hdf32) then + allocate (f4(s%dims(1),s%dims(2),s%dims(3))) + call h5Dread_f (dset_id,datatype_id,f4,s_dims,ierr) + do k=1,s%dims(3) + do j=1,s%dims(2) + do i=1,s%dims(1) + s%f(i,j,k) = REAL(f4(i,j,k),REAL64) + enddo + enddo + enddo + deallocate (f4) + else + allocate (f8(s%dims(1),s%dims(2),s%dims(3))) + call h5Dread_f (dset_id,datatype_id,f8,s_dims,ierr) + do k=1,s%dims(3) + do j=1,s%dims(2) + do i=1,s%dims(1) + s%f(i,j,k) = REAL(f8(i,j,k),REAL64) + enddo + enddo + enddo + deallocate (f8) + end if +! + deallocate(s_dims) +! + if (ierr.ne.0) then + write (*,*) + write (*,*) '### ERROR in RDH5:' + write (*,*) '### Error while reading the dataset.' + write (*,*) 'File name: ',trim(fname) + write (*,*) '[Error return (from H5DREAD_F) = ',ierr,']' + ierr=4 + return + end if +! +! ****** Close the hdf5 type descriptor. +! + call h5Tclose_f (datatype_id,ierr) +! +! ****** Check if there might be scales present, if so, read them. +! + if (n_members.gt.1) then +! +! ***** First check that the number of scale datasets match the # dim. +! + if (n_members-1.ne.s%ndim) then + write (*,*) + write (*,*) '### ERROR in RDH5:' + write (*,*) '### # scales does not match # dims.' + write (*,*) 'File name: ',trim(fname) + return + end if +! + s%scale=.true. +! +! ****** Loop through scales, make sure each is a scale, and read them. +! + do i=1,n_members-1 +! +! ****** Get the name of scale dataset. +! + call h5Gget_obj_info_idx_f (file_id,"/",i, & + obj_name,obj_type,ierr) +! +! ****** Open scale dataset. +! + call h5Dopen_f (file_id,trim(obj_name),dim_id,ierr) +! +! ****** Make sure the scale is a scale. +! + call h5DSis_scale_f (dim_id,is_scale,ierr) + if (.not.is_scale) then + write (*,*) + write (*,*) '### ERROR in RDH5:' + write (*,*) '### Scale is not a scale.' + write (*,*) 'File name: ',trim(fname) + return + end if +! +! ****** Get dimension of scale. +! + s_dims_i = s%dims(i) +! +! ****** Allocate scale. +! + allocate (s%scales(i)%f(s_dims_i(1))) +! +! ****** Get the floating-point precision of the scale. +! + call h5Dget_type_f (dim_id,datatype_id,ierr) + call h5Tget_precision_f (datatype_id,prec,ierr) +! +! ****** Read in the scale data. +! + if (prec.eq.32) then + allocate (f4dim(s_dims_i(1))) + call h5Dread_f (dim_id,datatype_id,f4dim,s_dims_i,ierr) + do j=1,s%dims(i) + s%scales(i)%f(j) = REAL(f4dim(j),REAL64) + end do + deallocate (f4dim) + elseif (prec.eq.64) then + allocate (f8dim(s_dims_i(1))) + call h5Dread_f (dim_id,datatype_id,f8dim,s_dims_i,ierr) + do j=1,s%dims(i) + s%scales(i)%f(j) = REAL(f8dim(j),REAL64) + end do + deallocate (f8dim) + end if +! +! ****** Close the scale dataset. +! + call h5Dclose_f (dim_id,ierr) +! + enddo +! +! ****** Allocate dummy scales (of length 1) for empty dimensions. +! + do i=s%ndim+1,3 + allocate (s%scales(i)%f(1)) + enddo + else +! +! ****** If scales are not present, allocate dummy +! ****** scales (of length 1) so that the pointers to the scales +! ****** are valid. +! + s%scale = .false. +! + allocate (s%scales(1)%f(1)) + allocate (s%scales(2)%f(1)) + allocate (s%scales(3)%f(1)) + end if +! +! ****** Close the dataset. +! + call h5Dclose_f (dset_id,ierr) +! +! ****** Close the file. +! + call h5Fclose_f (file_id,ierr) +! +! ****** Close FORTRAN interface. +! + call h5close_f (ierr) +! +end subroutine +!####################################################################### +subroutine wrhdf (fname,s,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Write a 1D, 2D, or 3D scientific data set to an HDF file. +! +!----------------------------------------------------------------------- +! +! ****** Input arguments: +! +! FNAME : [character(*)] +! HDF data file name to write to. +! +! S : [structure of type SDS] +! A structure that holds the field, its +! dimensions, and the scales, with the +! components described below. +! +! ****** Output arguments: +! +! IERR : [integer] +! IERR=0 is returned if the data set was written +! successfully. Otherwise, IERR is set to a +! nonzero value. +! +! ****** Components of structure S: +! +! NDIM : [integer] +! Number of dimensions in the data set. +! +! DIMS : [integer, dimension(3)] +! Number of points in the data set dimensions. +! Only DIMS(1 .. NDIM) are referenced. +! +! SCALE : [logical] +! Flag to indicate the presence of scales (axes) +! in the data set. SCALE=.false. means that scales +! are not being supplied; SCALE=.true. means that +! scales are being supplied. +! +! HDF32 : [logical] +! Flag to specify the precision of the data to +! be written to the file. Set HDF32=.true. to +! write 32-bit data, and HDF32=.false. to write +! 64-bit data. +! +! SCALES : [structure of type RP1D, dimension(3)] +! This array holds the pointers to the scales +! when SCALE=.true., and is not referenced +! otherwise. +! +! F : [real, pointer to a rank-3 array] +! This array holds the data set values. +! +!----------------------------------------------------------------------- +! + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + type(sds) :: s + integer :: ierr + intent(in) :: fname,s + intent(out) :: ierr +! +!----------------------------------------------------------------------- +! + integer :: iret,i,n +! +!----------------------------------------------------------------------- +! + ierr=0 +! +! ****** Check the number of dimensions. +! + if (s%ndim.le.0.or.s%ndim.gt.3) then + write (*,*) + write (*,*) '### ERROR in WRHDF:' + write (*,*) '### Could not write the SDS data.' + write (*,*) 'Invalid number of dimensions.' + write (*,*) 'Number of dimensions = ',s%ndim + write (*,*) 'File name: ',trim(fname) + ierr=1 + return + end if +! +! ****** Write hdf5 file if fname ends in '.h5'. +! + i=index(fname,'.h') + if (fname(i+1:i+2).eq.'h5') then + call wrh5 (fname,s,ierr) + else + print*,"HDF4 has been disabled." + ierr=-1 + end if +! +end subroutine +!####################################################################### +subroutine wrh5 (fname,s,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Write a 1D, 2D, or 3D scientific data set to an HDF5 file. +! +!----------------------------------------------------------------------- +! +! ****** Input arguments: +! +! FNAME : [character(*)] +! File name to write to. +! +! S : [structure of type DS] +! A structure that holds the field, its +! dimensions, and the scales, with the +! components described below. +! +! ****** Output arguments: +! +! IERR : [integer] +! IERR=0 is returned if the data set was written +! successfully. Otherwise, IERR is set to a +! nonzero value. +! +!----------------------------------------------------------------------- +! + use iso_fortran_env + use sds_def + use hdf5 + use h5ds +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + type(sds) :: s + integer :: ierr + intent(in) :: fname,s + intent(out) :: ierr +! +!----------------------------------------------------------------------- +! + character(8) :: dimname + integer :: i,j,k + integer(HID_T) :: file_id ! File identifier + integer(HID_T) :: dset_id ! Dataset identifier + integer(HID_T) :: dspace_id,dspacedim_id ! Dataspace identifiers + integer(HID_T) :: dim_id ! Dimension identifiers + integer(HSIZE_T),dimension(3) :: s_dims + integer(HSIZE_T),dimension(1) :: s_dims_i +! + real(REAL32), dimension(:,:,:), allocatable :: f4 + real(REAL32), dimension(:), allocatable :: f4dim + real(REAL64), dimension(:,:,:), allocatable :: f8 + real(REAL64), dimension(:), allocatable :: f8dim +! +!----------------------------------------------------------------------- +! +! ****** HDF5 calls are picky about the integer format for the dims +! ****** so the s%dims need to be converted to HSIZE_T integers. +! +! ****** Also, sometimes calls to wrhdf() for 1D and 2D datasets +! ****** do not have the unused dims(i) set to 1 (unset). +! ****** To avoid needing a function call to implicitly reshape +! ****** f(n), set the dims here. +! + do i=1,3 + if (i.le.s%ndim) then + s_dims(i) = INT(s%dims(i),HSIZE_T) + else + s_dims(i) = 1 + endif + end do +! +! ****** Initialize hdf5 interface. +! + call h5open_f (ierr) +! +! ****** Create the file. +! + call h5Fcreate_f (trim(fname),H5F_ACC_TRUNC_F,file_id,ierr) +! +! ****** Create the dataspace. +! + call h5Screate_simple_f (s%ndim,s_dims,dspace_id,ierr) +! +! ****** Create and write the dataset (convert s%f to proper type). +! + if (s%hdf32) then + allocate (f4(s_dims(1),s_dims(2),s_dims(3))) + do k=1,s%dims(3) + do j=1,s%dims(2) + do i=1,s%dims(1) + f4(i,j,k) = REAL(s%f(i,j,k),REAL32) + enddo + enddo + enddo + call h5Dcreate_f (file_id,'Data',H5T_NATIVE_REAL, & + dspace_id,dset_id,ierr) + call h5Dwrite_f (dset_id,H5T_NATIVE_REAL,f4,s_dims,ierr) + deallocate (f4) + else + allocate (f8(s_dims(1),s_dims(2),s_dims(3))) + do k=1,s%dims(3) + do j=1,s%dims(2) + do i=1,s%dims(1) + f8(i,j,k) = REAL(s%f(i,j,k),REAL64) + enddo + enddo + enddo + call h5Dcreate_f (file_id,'Data',H5T_NATIVE_DOUBLE, & + dspace_id,dset_id,ierr) + call h5Dwrite_f (dset_id,H5T_NATIVE_DOUBLE,f8,s_dims,ierr) + deallocate (f8) + endif + if (ierr.ne.0) then + write (*,*) + write (*,*) '### ERROR in WRH5:' + write (*,*) '### Could not write the dataset.' + write (*,*) 'File name: ',trim(fname) + write (*,*) '[Error return (from h5Dwrite_f) = ',ierr,']' + ierr=4 + return + end if +! +! ****** Check for scales. If present, add them to the hdf5 dataset. +! + if (s%scale) then + do i=1,s%ndim + if (i.eq.1) then + dimname='dim1' + elseif (i.eq.2) then + dimname='dim2' + elseif (i.eq.3) then + dimname='dim3' + endif + s_dims_i = s_dims(i) + call h5Screate_simple_f(1,s_dims_i,dspacedim_id,ierr) + if (s%hdf32) then + allocate (f4dim(s_dims_i(1))) + do j=1,s%dims(i) + f4dim(j) = REAL(s%scales(i)%f(j),REAL32) + end do + call h5Dcreate_f (file_id,dimname,H5T_NATIVE_REAL, & + dspacedim_id,dim_id,ierr) + call h5Dwrite_f (dim_id,H5T_NATIVE_REAL, & + f4dim,s_dims_i,ierr) + deallocate (f4dim) + else + allocate (f8dim(s_dims_i(1))) + do j=1,s%dims(i) + f8dim(j) = REAL(s%scales(i)%f(j),REAL64) + end do + call h5Dcreate_f (file_id,dimname,H5T_NATIVE_DOUBLE, & + dspacedim_id,dim_id,ierr) + call h5Dwrite_f (dim_id,H5T_NATIVE_DOUBLE, & + f8dim,s_dims_i,ierr) + deallocate (f8dim) + endif + call h5DSset_scale_f (dim_id,ierr,dimname) + call h5DSattach_scale_f (dset_id,dim_id,i,ierr) + call h5DSset_label_f(dset_id, i, dimname, ierr) + call h5Dclose_f (dim_id,ierr) + call h5Sclose_f (dspacedim_id,ierr) + end do + if (ierr.ne.0) then + write (*,*) + write (*,*) '### ERROR in WRH5:' + write (*,*) '### Could not write the scales.' + write (*,*) 'File name: ',trim(fname) + ierr = 5 + return + endif + endif +! +! ****** Close the dataset. +! + call h5Dclose_f (dset_id,ierr) +! +! ****** Close the dataspace. +! + call h5Sclose_f (dspace_id,ierr) +! +! ****** Close the file. +! + call h5Fclose_f (file_id,ierr) +! +! ****** Close the hdf5 interface. +! + call h5close_f (ierr) +! +end subroutine +!####################################################################### +subroutine rdhdf_1d (fname,scale,nx,f,x,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Read a 1D scientific data set from an HDF file. +! +!----------------------------------------------------------------------- +! +! ****** This routine calls routine RDHDF to read the file. +! +!----------------------------------------------------------------------- +! + use iso_fortran_env + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + logical :: scale + integer :: nx,i + real(REAL64), dimension(:), pointer :: f + real(REAL64), dimension(:), pointer :: x + integer :: ierr + intent(in) :: fname + intent(out) :: scale,nx,ierr +! +!----------------------------------------------------------------------- +! +! ****** Declaration for the SDS structure. +! + type(sds) :: s +! +!----------------------------------------------------------------------- +! +! ****** Read the data set. +! + call rdhdf (fname,s,ierr) +! + if (ierr.ne.0) return +! +! ****** Check that this is a 1D data set. +! + if (s%ndim.ne.1) then + write (*,*) + write (*,*) '### ERROR in RDHDF_1D:' + write (*,*) '### The HDF file does not contain a 1D data set.' + write (*,*) 'File name: ',trim(fname) + ierr=3 + return + end if +! +! ****** Set the output arguments. +! + nx=s%dims(1) + scale=s%scale + x=>s%scales(1)%f +! + allocate (f(nx)) + do i=1,nx + f(i) = REAL(s%f(i,1,1),REAL64) + enddo + deallocate (s%f) +! +end subroutine +!####################################################################### +subroutine rdhdf_2d (fname,scale,nx,ny,f,x,y,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Read a 2D scientific data set from an HDF file. +! +!----------------------------------------------------------------------- +! +! ****** This routine calls routine RDHDF to read the file. +! +!----------------------------------------------------------------------- +! + use iso_fortran_env + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + logical :: scale + integer :: nx,ny,i + real(REAL64), dimension(:,:), pointer :: f + real(REAL64), dimension(:), pointer :: x,y + integer :: ierr + intent(in) :: fname + intent(out) :: scale,nx,ny,ierr +! +!----------------------------------------------------------------------- +! +! ****** Declaration for the SDS structure. +! + type(sds) :: s +! +!----------------------------------------------------------------------- +! +! ****** Read the data set. +! + call rdhdf (fname,s,ierr) +! + if (ierr.ne.0) return +! +! ****** Check that this is a 2D data set. +! + if (s%ndim.ne.2) then + write (*,*) + write (*,*) '### ERROR in RDHDF_2D:' + write (*,*) '### The HDF file does not contain a 2D data set.' + write (*,*) 'File name: ',trim(fname) + ierr=3 + return + end if +! +! ****** Set the output arguments. +! +! + nx=s%dims(1) + ny=s%dims(2) + scale=s%scale + x=>s%scales(1)%f + y=>s%scales(2)%f +! + allocate (f(nx,ny)) + f(:,:)=s%f(:,:,1) + deallocate (s%f) +! +end subroutine +!####################################################################### +subroutine rdhdf_3d (fname,scale,nx,ny,nz,f,x,y,z,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Read a 3D scientific data set from an HDF file. +! +!----------------------------------------------------------------------- +! +! ****** This routine calls routine RDHDF to read the file. +! +!----------------------------------------------------------------------- +! + use iso_fortran_env + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + logical :: scale + integer :: nx,ny,nz + real(REAL64), dimension(:,:,:), pointer :: f + real(REAL64), dimension(:), pointer :: x,y,z + integer :: ierr + intent(in) :: fname + intent(out) :: scale,nx,ny,nz,ierr +! +!----------------------------------------------------------------------- +! +! ****** Declaration for the SDS structure. +! + type(sds) :: s +! +!----------------------------------------------------------------------- +! +! ****** Read the data set. +! + call rdhdf (fname,s,ierr) +! + if (ierr.ne.0) return +! +! ****** Check that this is a 3D data set. +! + if (s%ndim.ne.3) then + write (*,*) + write (*,*) '### ERROR in RDHDF_3D:' + write (*,*) '### The HDF file does not contain a 3D data set.' + write (*,*) 'File name: ',trim(fname) + ierr=3 + return + end if +! +! ****** Set the output arguments. +! + nx=s%dims(1) + ny=s%dims(2) + nz=s%dims(3) + scale=s%scale + x=>s%scales(1)%f + y=>s%scales(2)%f + z=>s%scales(3)%f + f=>s%f +! +end subroutine +!####################################################################### +subroutine wrhdf_1d (fname,scale,nx,f,x,hdf32,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Write a 1D scientific data set to an HDF file. +! +!----------------------------------------------------------------------- +! +! ****** This routine calls routine WRHDF to write the file. +! +!----------------------------------------------------------------------- +! + use iso_fortran_env + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + logical :: scale + integer :: nx + real(REAL64), dimension(nx,1,1), target :: f + real(REAL64), dimension(nx), target :: x + logical :: hdf32 + integer :: ierr + intent(in) :: fname,scale,nx,f,x,hdf32 + intent(out) :: ierr +! +!----------------------------------------------------------------------- +! +! ****** Declaration for the SDS structure. +! + type(sds) :: s +! +!----------------------------------------------------------------------- +! +! ****** Set the structure components. +! + s%ndim=1 + s%dims(1)=nx + s%dims(2)=1 + s%dims(3)=1 + s%scale=scale + s%hdf32=hdf32 + if (scale) then + s%scales(1)%f=>x + else + nullify (s%scales(1)%f) + end if + nullify (s%scales(2)%f) + nullify (s%scales(3)%f) + s%f=>f +! +! ****** Write the data set. +! + call wrhdf (fname,s,ierr) +! + if (ierr.ne.0) then + write (*,*) + write (*,*) '### ERROR in WRHDF_1D:' + write (*,*) '### Could not write the 1D data set.' + write (*,*) 'File name: ',trim(fname) + return + end if +! +end subroutine +!####################################################################### +subroutine wrhdf_2d (fname,scale,nx,ny,f,x,y,hdf32,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Write a 2D scientific data set to an HDF file. +! +!----------------------------------------------------------------------- +! +! ****** This routine calls routine WRHDF to write the file. +! +!----------------------------------------------------------------------- +! + use iso_fortran_env + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + logical :: scale + integer :: nx,ny + real(REAL64), dimension(nx,ny,1), target :: f + real(REAL64), dimension(nx), target :: x + real(REAL64), dimension(ny), target :: y + logical :: hdf32 + integer :: ierr + intent(in) :: fname,scale,nx,ny,f,x,y,hdf32 + intent(out) :: ierr +! +!----------------------------------------------------------------------- +! +! ****** Declaration for the SDS structure. +! + type(sds) :: s +! +!----------------------------------------------------------------------- +! +! ****** Set the structure components. +! + s%ndim=2 + s%dims(1)=nx + s%dims(2)=ny + s%dims(3)=1 + s%scale=scale + s%hdf32=hdf32 + if (scale) then + s%scales(1)%f=>x + s%scales(2)%f=>y + else + nullify (s%scales(1)%f) + nullify (s%scales(2)%f) + end if + nullify (s%scales(3)%f) + s%f=>f +! +! ****** Write the data set. +! + call wrhdf (fname,s,ierr) +! + if (ierr.ne.0) then + write (*,*) + write (*,*) '### ERROR in WRHDF_2D:' + write (*,*) '### Could not write the 2D data set.' + write (*,*) 'File name: ',trim(fname) + return + end if +! +end subroutine +!####################################################################### +subroutine wrhdf_3d (fname,scale,nx,ny,nz,f,x,y,z,hdf32,ierr) +! +!----------------------------------------------------------------------- +! +! ****** Write a 3D scientific data set to an HDF file. +! +!----------------------------------------------------------------------- +! +! ****** This routine calls routine WRHDF to write the file. +! +!----------------------------------------------------------------------- +! + use iso_fortran_env + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + character(*) :: fname + logical :: scale + integer :: nx,ny,nz + real(REAL64), dimension(nx,ny,nz), target :: f + real(REAL64), dimension(nx), target :: x + real(REAL64), dimension(ny), target :: y + real(REAL64), dimension(nz), target :: z + logical :: hdf32 + integer :: ierr + intent(in) :: fname,scale,nx,ny,nz,f,x,y,z,hdf32 + intent(out) :: ierr +! +!----------------------------------------------------------------------- +! +! ****** Declaration for the SDS structure. +! + type(sds) :: s +! +!----------------------------------------------------------------------- +! +! ****** Set the structure components. +! + s%ndim=3 + s%dims(1)=nx + s%dims(2)=ny + s%dims(3)=nz + s%scale=scale + s%hdf32=hdf32 + if (scale) then + s%scales(1)%f=>x + s%scales(2)%f=>y + s%scales(3)%f=>z + else + nullify (s%scales(1)%f) + nullify (s%scales(2)%f) + nullify (s%scales(3)%f) + end if + s%f=>f +! +! ****** Write the data set. +! + call wrhdf (fname,s,ierr) +! + if (ierr.ne.0) then + write (*,*) + write (*,*) '### ERROR in WRHDF_3D:' + write (*,*) '### Could not write the 3D data set.' + write (*,*) 'File name: ',trim(fname) + return + end if +! +end subroutine +!####################################################################### +subroutine deallocate_sds (s) +! +!----------------------------------------------------------------------- +! +! ****** Deallocate the memory used by the SDS in structure S. +! +!----------------------------------------------------------------------- +! + use sds_def +! +!----------------------------------------------------------------------- +! + implicit none +! +!----------------------------------------------------------------------- +! + type(sds) :: s +! +!----------------------------------------------------------------------- +! + if (associated(s%f)) deallocate (s%f) +! + if (associated(s%scales(1)%f)) deallocate (s%scales(1)%f) + if (associated(s%scales(2)%f)) deallocate (s%scales(2)%f) + if (associated(s%scales(3)%f)) deallocate (s%scales(3)%f) +! +end subroutine +! diff --git a/src/stdpar/pot3d.F b/src/stdpar/pot3d.F index a01c3e3..7c839a0 100644 --- a/src/stdpar/pot3d.F +++ b/src/stdpar/pot3d.F @@ -56,6 +56,27 @@ module ident character(*), parameter :: update='11/01/2023' c end module +c####################################################################### + module number_types +c +c----------------------------------------------------------------------- +c ****** Basic number types. +c ****** This module is used to set the default precision for REALs. +c----------------------------------------------------------------------- +c + use iso_fortran_env +c +c----------------------------------------------------------------------- +c + implicit none +c + integer, parameter :: KIND_REAL_4=REAL32 + integer, parameter :: KIND_REAL_8=REAL64 + integer, parameter :: KIND_REAL_16=max(REAL128,REAL64) +c + integer, parameter :: r_typ=KIND_REAL_8 +c + end module c####################################################################### module number_types_pc c @@ -6358,8 +6379,8 @@ subroutine readbr (fname,br0_g,ierr) use number_types use global_dims use global_mesh - use rdhdf_2d_interface use vars + use rdhdf_2d_interface c c----------------------------------------------------------------------- c diff --git a/src/zm_parse.f b/src/zm_parse.f deleted file mode 100644 index 0b120ef..0000000 --- a/src/zm_parse.f +++ /dev/null @@ -1,1946 +0,0 @@ -c -c----------------------------------------------------------------------- -c -c ****** Source to build the parsing library. -c ****** These routines are used by Zoran Mikic's tools. -c -c ********************************************************************** -c -c Copyright 2018 Predictive Science Inc. -c -c Licensed under the Apache License, Version 2.0 (the "License"); -c you may not use this file except in compliance with the License. -c You may obtain a copy of the License at -c -c http://www.apache.org/licenses/LICENSE-2.0 -c -c Unless required by applicable law or agreed to in writing, software -c distributed under the License is distributed on an "AS IS" BASIS, -c WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or -c implied. -c See the License for the specific language governing permissions and -c limitations under the License. -c -c ********************************************************************** -c -c -c----------------------------------------------------------------------- -c -c 07/29/2003, ZM, Version 1.00: -c -c - Original version of the parsing library. -c This library was put together to facilitate the -c development of ZM's tools. -c It includes routines to parse the command line. -c The code was cleaned up to use standard FORTRAN90. -c -c 01/17/2005, ZM, Version 1.01: -c -c - Added the function NARGS_SPECIFIED to return the -c number of arguments specified on the command line. -c -c 10/28/2005, ZM, Version 1.02: -c -c - Added the functions LOAD_LIST_OF_REALS and -c LOAD_LIST_OF_INTS that can be used to parse -c arguments that contain lists of real or integer values. -c - Changed the length of temporary arguments to equal -c 512 characters to accomodate long file names. -c -c 10/31/2006, ZM, Version 1.03: -c -c - Removed the EXTERNAL declarations for GETARG and IARGC -c since these are now intrinsic routines in the -c Intel 9.1 compiler. -c -c 03/10/2008, ZM, Version 1.04: -c -c - Added the LCASE and UCASE functions to convert strings -c to lowercase and uppercase. -c -c 04/27/2018, ZM, Version 1.05: -c -c - Added the ability to specify group-sets of keywords. -c This extends the flexibility of the syntax allowed -c for keywords. For example, it allows a syntax in which -c only one of a set of keywords is allowed to be specified -c at a time. -c -c 01/08/2019, RC, Version 1.06: -c -c - Added "append" option to ffopen. -c -c----------------------------------------------------------------------- -c -c####################################################################### - module parselib_ident -c - character(*), parameter :: cname='PARSELIB' - character(*), parameter :: cvers='1.06' - character(*), parameter :: cdate='01/08/2019' -c - end module -c####################################################################### - module parse_args -c - use string_def -c - implicit none -c -c----------------------------------------------------------------------- -c ****** Argument descriptor and storage for command-line arguments. -c----------------------------------------------------------------------- -c -c ****** Structure to hold a group-set. -c - type :: group_set - integer :: group_type - integer :: n_members - integer, dimension(:), allocatable :: members - logical :: required - logical :: only_one - end type -c -c ****** Structure to hold an argument. -c - type :: arg_descriptor - integer :: group - logical :: set - logical :: required - type(string) :: keyword - type(string) :: name - type(string) :: value - end type -c -c ****** Maximum number of arguments. -c - integer, parameter :: mxarg=100 -c -c ****** Number of arguments defined. -c - integer :: nargs -c -c ****** Argument descriptor. -c - type(arg_descriptor), dimension(mxarg) :: args -c -c ****** Number of arguments specified. -c - integer :: nargs_spec -c -c ****** Maximum number of group-sets. -c - integer, parameter :: mxgroup_sets=10 -c -c ****** Number of group-sets defined. -c - integer :: ngroup_sets=0 -c -c ****** Group-set descriptor. -c - type(group_set), dimension(mxgroup_sets) :: group_sets -c - end module -c####################################################################### - subroutine ffopen (iun,fname,mode,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Open file FNAME and link it to unit IUN. -c -c ****** If there is an error, this routine returns IERR.ne.0. -c -c----------------------------------------------------------------------- -c -c ****** When MODE='r', the file must exist. -c ****** When MODE='w', the file is created. -c ****** When MODE='rw', the file must exist, but can be overwritten. -c ****** When MODE='a', the file is created if it does not exist, -c ****** otherwise, it is appended. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: iun - character(*) :: fname - character(*) :: mode - integer :: ierr - logical :: ex -c -c----------------------------------------------------------------------- -c - ierr=0 -c - if (mode.eq.'r') then - open (iun,file=fname,form="FORMATTED",status='old',err=900) - else if (mode.eq.'rw') then - open (iun,file=fname,form="FORMATTED",status='replace',err=900) - else if (mode.eq.'w') then - open (iun,file=fname,form="FORMATTED",status='new',err=900) - elseif (mode.eq.'a') then - inquire(file=fname, exist=ex) - if (ex) then - open (iun,file=fname,form="FORMATTED", - & position='append',err=900) - else - open (iun,file=fname,form="FORMATTED",status='new',err=900) - end if - else - write (*,*) - write (*,*) '### ERROR in FFOPEN:' - write (*,*) '### Invalid MODE requested.' - write (*,*) 'MODE = ',mode - write (*,*) 'File name: ',trim(fname) - ierr=2 - return - end if -c - return -c - 900 continue -c - write (*,*) - write (*,*) '### ERROR in FFOPEN:' - write (*,*) '### Error while opening the requested file.' - write (*,*) 'File name: ',trim(fname) - write (*,*) 'MODE = ',mode - ierr=1 -c - return - end -c####################################################################### - function lcase (s) -c -c----------------------------------------------------------------------- -c -c ****** Convert the string S into lowercase letters and return it as -c ****** the function result. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*), intent(in) :: s - character(len(s)) :: lcase -c -c----------------------------------------------------------------------- -c - integer :: i,ic -c -c----------------------------------------------------------------------- -c - lcase=' ' -c - do i=1,len_trim(s) - ic=iachar(s(i:i)) - if (ic.ge.65.and.ic.le.90) then - ic=ic+32 - end if - lcase(i:i)=achar(ic) - end do -c - return - end -c####################################################################### - function ucase (s) -c -c----------------------------------------------------------------------- -c -c ****** Convert the string S into uppercase letters and return it as -c ****** the function result. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*), intent(in) :: s - character(len(s)) :: ucase -c -c----------------------------------------------------------------------- -c - integer :: i,ic -c -c----------------------------------------------------------------------- -c - ucase=' ' -c - do i=1,len_trim(s) - ic=iachar(s(i:i)) - if (ic.ge.97.and.ic.le.122) then - ic=ic-32 - end if - ucase(i:i)=achar(ic) - end do -c - return - end -c####################################################################### - subroutine parse (errmsg,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Parse the command line. -c -c----------------------------------------------------------------------- -c -c ****** The syntax for the keyword/argument items can be defined -c ****** by using routine DEFARG. -c -c ****** On return, IERR=0 indicates that the command line was -c ****** parsed successfully. -c -c ****** IERR=1 indicates that no arguments were present. This -c ****** is usually used to print the usage line. -c -c ****** IERR=2 indicates that a syntax error occured. -c -c ****** IERR=3 indicates that one or more required arguments -c ****** was not supplied. -c -c ****** When IERR=2 or IERR=3, an error message is put into -c ****** character string ERRMSG. -c -c----------------------------------------------------------------------- -c - use syntax - use parse_args - use get_str_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: errmsg - integer :: ierr -c -c----------------------------------------------------------------------- -c -c ****** Command line arguments. -c - integer :: iargc -c -c----------------------------------------------------------------------- -c -c ****** Carriage return. -c - character(1), parameter :: cr=achar(10) -c -c----------------------------------------------------------------------- -c - integer, external :: matchkw - integer, external :: nextarg -c -c----------------------------------------------------------------------- -c - character(512) :: arg - integer :: na,ia,ia0,iarg,ls,i,j,n -c -c----------------------------------------------------------------------- -c -c ****** Initialization. -c - ierr=0 - nargs_spec=0 - errmsg=' ' -c -c ****** Get the number of command line arguments. -c - na=iargc() - if (na.eq.0) then - ierr=1 - return - end if -c - ia=1 - 200 continue -c - ia0=ia -c -c ****** Process arguments with syntax: -c - if (na-ia+1.ge.2) then - call getarg (ia,arg) - iarg=matchkw(GROUP_KA,trim(arg)) - if (iarg.gt.0) then - if (.not.args(iarg)%set) then - ia=ia+1 - call getarg (ia,arg) - call delete_str (args(iarg)%value) - call put_str (trim(arg),args(iarg)%value) - args(iarg)%set=.true. - ia=ia+1 - nargs_spec=nargs_spec+1 - go to 300 - end if - end if - end if -c -c ****** Process arguments with syntax: -c - if (na-ia+1.ge.3) then - call getarg (ia,arg) - iarg=matchkw(GROUP_KAA,trim(arg)) - if (iarg.gt.0) then - if (.not.args(iarg)%set) then - ia=ia+1 - call getarg (ia,arg) - ls=len_trim(arg) - ls=ls+1 - arg(ls:ls)=' ' - ia=ia+1 - call getarg (ia,arg(ls+1:)) - call delete_str (args(iarg)%value) - call put_str (trim(arg),args(iarg)%value) - args(iarg)%set=.true. - ia=ia+1 - nargs_spec=nargs_spec+1 - go to 300 - end if - end if - end if -c -c ****** Process arguments with syntax: -c - if (na-ia+1.ge.1) then - call getarg (ia,arg) - iarg=matchkw(GROUP_K,trim(arg)) - if (iarg.gt.0) then - if (.not.args(iarg)%set) then - call delete_str (args(iarg)%value) - call put_str (' ',args(iarg)%value) - args(iarg)%set=.true. - ia=ia+1 - nargs_spec=nargs_spec+1 - go to 300 - end if - end if - end if -c -c ****** Process arguments with syntax: -c - if (na-ia+1.ge.1) then - iarg=nextarg(GROUP_A) - if (iarg.gt.0) then - call getarg (ia,arg) - call delete_str (args(iarg)%value) - call put_str (trim(arg),args(iarg)%value) - args(iarg)%set=.true. - ia=ia+1 - nargs_spec=nargs_spec+1 - go to 300 - end if - end if -c - 300 continue -c -c ****** Check that an argument was found. -c - if (ia.eq.ia0) then - ierr=2 - errmsg='### Syntax error.' - return - end if -c -c ****** Keep processing arguments until done. -c - if (na-ia+1.gt.0) go to 200 -c -c ****** Check that the required arguments were supplied. -c - do i=1,nargs - if (args(i)%required.and..not.args(i)%set) then - ierr=3 - errmsg='### A required argument was not supplied.' - return - end if - enddo -c -c ****** Check that the group-set keywords were specified -c ****** properly. -c -c ****** Check the "only one" and "required" group-set keywords. -c - do i=1,ngroup_sets - n=0 - do j=1,group_sets(i)%n_members - iarg=group_sets(i)%members(j) - if (args(iarg)%set) n=n+1 - enddo - if (n.gt.1.and.group_sets(i)%only_one) then - ierr=2 - errmsg='### Syntax error.' - errmsg=trim(errmsg)//cr//cr - errmsg=trim(errmsg)//' ### You can only specify'// - & ' one of the following keywords:' - do j=1,group_sets(i)%n_members - iarg=group_sets(i)%members(j) - errmsg=trim(errmsg)//cr//' '// - & get_str(args(iarg)%keyword) - enddo - return - else if (n.eq.0.and.group_sets(i)%required) then - ierr=2 - errmsg='### Syntax error.' - errmsg=trim(errmsg)//cr//cr - errmsg=trim(errmsg)//' ### You must specify'// - & ' one of the following keywords:' - do j=1,group_sets(i)%n_members - iarg=group_sets(i)%members(j) - errmsg=trim(errmsg)//cr//' '// - & get_str(args(iarg)%keyword) - enddo - return - end if - enddo -c - return - end -c####################################################################### - subroutine get_usage_line (usage) -c -c----------------------------------------------------------------------- -c -c ****** Construct the usage line in paragraph USAGE. -c -c ****** Use routine PRINT_PAR to write the usage line. -c -c----------------------------------------------------------------------- -c - use syntax - use parse_args - use paragraph_def - use new_par_interface - use add_line_interface - use get_str_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(paragraph), pointer :: usage -c -c----------------------------------------------------------------------- -c -c ****** Right margin for printing the usage line. -c - integer, parameter :: rmargin=78 -c -c----------------------------------------------------------------------- -c - character(512) :: line - integer :: iarg,n0 - type(paragraph), pointer :: current_par -c - logical, dimension(:), allocatable :: done - integer :: i,j,ix - logical :: belongs_to_gs - character(512) :: str -c -c----------------------------------------------------------------------- -c -c ****** Construct the usage line in USAGE. -c - call new_par (usage) - current_par=>usage -c -c ****** Start with the command name (as invoked). -c - call getarg (0,line) -c -c ****** Allocate storage. -c - allocate (done(nargs)) - done=.false. -c - iarg=1 -c -c ****** Add the arguments one by one. -c - do while (iarg.le.nargs) -c - if (done(iarg)) then - iarg=iarg+1 - cycle - end if -c -c ****** Add the syntax for the next argument to LINE. -c - n0=len_trim(line) -c -c ****** Check if this argument belongs to a "only one" group-set. -c ****** If so, collect all members of the group-set. Otherwise, just -c ****** process this argument. -c - belongs_to_gs=.false. - do i=1,ngroup_sets - if (group_sets(i)%only_one) then - do j=1,group_sets(i)%n_members - if (group_sets(i)%members(j).eq.iarg) then - ix=i - belongs_to_gs=.true. - exit - end if - enddo - end if - if (belongs_to_gs) exit - enddo -c - if (belongs_to_gs) then -c -c ****** This argument belongs to a "only one" group-set: -c ****** process all its members. -c - str='' - do j=1,group_sets(ix)%n_members - i=group_sets(ix)%members(j) - done(i)=.true. - str=trim(str)//get_str(args(i)%keyword) - if (j.lt.group_sets(ix)%n_members) then - str=trim(str)//'|' - end if - enddo - if (group_sets(ix)%required) then - line=trim(line)//' '//trim(str) - else - line=trim(line)//' ['//trim(str)//']' - end if -c - else -c -c ****** This argument does not belong to a "only one" group-set: -c ****** process just this single argument. -c - if (args(iarg)%required) then - line=trim(line)//' '//get_str(args(iarg)%keyword) - else - line=trim(line)//' ['//get_str(args(iarg)%keyword) - end if - line=trim(line)//' '//get_str(args(iarg)%name) - if (.not.args(iarg)%required) then - line=trim(line)//']' - end if - done(iarg)=.true. -c - end if -c -c ****** Check if the addition of the argument causes the line -c ****** to wrap; if it does, break the line prior to the -c ****** argument text. -c - if (len_trim(line).gt.rmargin) then - call add_line (line(1:n0),current_par) - line=' '//line(n0+1:) - end if -c -c ****** If the line is still too long, force a break at RMARGIN -c ****** until the line is shorter than RMARGIN. -c - do while (len_trim(line).gt.rmargin) - call add_line (line(1:rmargin),current_par) - line=' '//line(rmargin+1:) - enddo -c -c ****** Process the next argument. -c - iarg=iarg+1 -c - enddo -c -c ****** Add the last line to the paragraph. -c - if (line.ne.' ') call add_line (trim(line),current_par) -c -c ****** Deallocate storage. -c - deallocate (done) -c - return - end -c####################################################################### - subroutine defarg (group,keyword,default,name) -c -c----------------------------------------------------------------------- -c -c ****** Define the syntax for a command line argument item. -c -c----------------------------------------------------------------------- -c -c ****** GROUP is the syntax group index; -c ****** KEYWORD is the keyword; -c ****** DEFAULT is the default value of the argument; -c ****** NAME is the name of the argument (for use in -c ****** constructing the usage line). -c -c----------------------------------------------------------------------- -c - use syntax - use parse_args -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: group - character(*) :: keyword - character(*) :: default - character(*) :: name -c -c----------------------------------------------------------------------- -c - integer :: i,n,ix - character(512), dimension(:), allocatable :: names -c -c----------------------------------------------------------------------- -c - integer, external :: number_of_names - integer, external :: match_keyword -c -c----------------------------------------------------------------------- -c -c ****** Check that the group index is valid. -c - if (group.lt.0.or.group.gt.ngroups) then - write (*,*) - write (*,*) '### ERROR in DEFARG:' - write (*,*) '### An invalid group index was specified.' - write (*,*) 'Group index = ',group - write (*,*) 'Keyword: ',trim(keyword) - if (name.ne.' ') write (*,*) 'Name: ',trim(name) - write (*,*) - write (*,*) '### This indicates a programming error'// - & ' in the syntax definition and use.' - call exit (1) - end if -c -c ****** Check for a null keyword. -c - if (keyword.eq.' ') then - write (*,*) - write (*,*) '### ERROR in DEFARG:' - write (*,*) '### The keyword is null.' - write (*,*) 'Group index = ',group - if (name.ne.' ') write (*,*) 'Name: ',trim(name) - write (*,*) - write (*,*) '### This indicates a programming error'// - & ' in the syntax definition and use.' - call exit (1) - end if -c -c ****** Check if this defines a group-set. -c - if (group.eq.GROUP_K_ONE_ONLY.or. - & group.eq.GROUP_K_ONE_OR_NONE) then -c -c ****** Increment the group-set counter. -c - if (ngroup_sets.ge.mxgroup_sets) then - write (*,*) - write (*,*) '### ERROR in DEFARG:' - write (*,*) '### Exceeded the number of allowed group-sets.' - write (*,*) 'Maximum number of group-sets = ',mxgroup_sets - write (*,*) 'Group index = ',group - write (*,*) 'Keyword: ',trim(keyword) - write (*,*) - write (*,*) '### This indicates a programming error'// - & ' in the syntax definition and use.' - call exit (1) - end if -c - ngroup_sets=ngroup_sets+1 -c -c ****** Set the group type. -c - group_sets(ngroup_sets)%group_type=group -c -c ****** Set the appropriate group-set flags. -c - if (group.eq.GROUP_K_ONE_ONLY) then - group_sets(ngroup_sets)%required=.true. - group_sets(ngroup_sets)%only_one=.true. - else if (group.eq.GROUP_K_ONE_OR_NONE) then - group_sets(ngroup_sets)%required=.false. - group_sets(ngroup_sets)%only_one=.true. - else - group_sets(ngroup_sets)%required=.false. - group_sets(ngroup_sets)%only_one=.false. - end if -c -c ****** Mark the members in the group-set. -c - n=number_of_names(keyword) - allocate (names(n)) - allocate (group_sets(ngroup_sets)%members(n)) - group_sets(ngroup_sets)%n_members=n - call load_list_of_names (keyword,n,names) - do i=1,n - ix=match_keyword(names(i)) - if (ix.eq.0) then - write (*,*) - write (*,*) '### ERROR in DEFARG:' - write (*,*) '### Error while defining a group set.' - write (*,*) 'Could not match a group-set keyword to'// - & ' an existing keyword:' - write (*,*) 'Group-set keyword: ',trim(keyword) - write (*,*) 'Keyword not matched: ',trim(names(i)) - write (*,*) - write (*,*) '### This indicates a programming error'// - & ' in the syntax definition and use.' - call exit (1) - end if - group_sets(ngroup_sets)%members(i)=ix - enddo - deallocate (names) -c - return -c - end if -c -c ****** Increment the argument counter. -c - if (nargs.ge.mxarg) then - write (*,*) - write (*,*) '### ERROR in DEFARG:' - write (*,*) '### Exceeded the number of allowed arguments.' - write (*,*) 'Maximum number of arguments = ',mxarg - write (*,*) 'Group index = ',group - write (*,*) 'Keyword: ',trim(keyword) - if (name.ne.' ') write (*,*) 'Name: ',trim(name) - write (*,*) - write (*,*) '### This indicates a programming error'// - & ' in the syntax definition and use.' - call exit (1) - end if -c - nargs=nargs+1 -c -c ****** Store the group index and keyword. -c -c ****** For group GROUP_A (single arguments), the name of the -c ****** argument is passed as the "keyword". -c - args(nargs)%group=group - call put_str (trim(keyword),args(nargs)%keyword) -c -c ****** Initialize the flag that indicates whether an argument -c ****** has been set. -c - args(nargs)%set=.false. -c -c ****** If a default argument was supplied, the argument -c ****** does not have to be set. Use DEFAULT=' ' to -c ****** indicate that an argument is required. -c -c ****** If a default argument has been supplied, store it in -c ****** ARGS(nargs)%VALUE. If there is no default, -c ****** set ARGS(nargs)%VALUE to an empty string. -c -c ****** Since group GROUP_K doesn't have an argument, -c ****** DEFAULT is ignored for this group. -c - if (group.eq.GROUP_K) then - args(nargs)%required=.false. - call put_str (' ',args(nargs)%value) - else - if (default.eq.' ') then - args(nargs)%required=.true. - call put_str (' ',args(nargs)%value) - else - args(nargs)%required=.false. - call put_str (trim(default),args(nargs)%value) - end if - end if -c -c ****** Store the argument name. For groups GROUP_K (keywords) -c ****** and GROUP_A (single arguments), there is no argument name, -c ****** so NAME is ignored. -c - if (group.eq.GROUP_K.or.group.eq.GROUP_A) then - call put_str (' ',args(nargs)%name) - else - call put_str (trim(name),args(nargs)%name) - end if -c - return - end -c####################################################################### - subroutine fetcharg (keyword,set,arg) -c -c----------------------------------------------------------------------- -c -c ****** Fetch the value of the argument corresponding to -c ****** keyword KEYWORD. -c -c----------------------------------------------------------------------- -c -c ****** If KEYWORD is a keyword-type argument (GROUP_K), return -c ****** its setting through variable SET. The variable ARG should -c ****** be ignored for this type of keyword. -c -c ****** For keywords with arguments (GROUP_A, GROUP_KA, and -c ****** GROUP_KAA), return the value of the arguments in ARG, -c ****** and return SET=.true. if they were set via the command line; -c ****** otherwise, return SET=.false.. -c -c----------------------------------------------------------------------- -c - use parse_args - use get_str_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: keyword - logical :: set - character(*) :: arg -c -c----------------------------------------------------------------------- -c - integer :: i -c -c----------------------------------------------------------------------- -c - do i=nargs,1,-1 - if (keyword.eq.get_str(args(i)%keyword)) go to 100 - enddo -c - write (*,*) - write (*,*) '### ERROR in FETCHARG:' - write (*,*) '### The requested keyword could not be matched.' - write (*,*) 'Keyword = ',trim(keyword) - write (*,*) - write (*,*) '### This indicates a programming error'// - & ' in the syntax definition and use.' - call exit (1) -c - 100 continue -c - set=args(i)%set - arg=get_str(args(i)%value) -c - return - end -c####################################################################### - function matchkw (group,keyword) -c -c----------------------------------------------------------------------- -c -c ****** Match keyword KEYWORD against the list of keywords in -c ****** group GROUP. -c -c ****** If found, set the function value to the corresponding -c ****** argument number. Otherwise, return MATCHKW=0. -c -c----------------------------------------------------------------------- -c - use parse_args - use get_str_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: group - character(*) :: keyword - integer :: matchkw -c -c----------------------------------------------------------------------- -c - integer :: i -c -c----------------------------------------------------------------------- -c - matchkw=0 -c - do i=nargs,1,-1 - if (group.eq.args(i)%group) then - if (keyword.eq.get_str(args(i)%keyword)) then - matchkw=i - return - end if - end if - enddo -c - return - end -c####################################################################### - function nextarg (group) -c -c----------------------------------------------------------------------- -c -c ****** Find the position of the next argument in group GROUP -c ****** that has not been set. -c -c----------------------------------------------------------------------- -c -c ****** If an empty slot is found, set the function value -c ****** to the corresponding argument number. -c -c ****** Otherwise, return NXTARG=0. -c -c----------------------------------------------------------------------- -c - use parse_args -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: group - integer :: nextarg -c -c----------------------------------------------------------------------- -c - integer :: i -c -c----------------------------------------------------------------------- -c - nextarg=0 -c - do i=1,nargs - if (group.eq.args(i)%group) then - if (.not.args(i)%set) then - nextarg=i - return - end if - end if - enddo -c - return - end -c####################################################################### - subroutine nargs_specified (n) -c -c----------------------------------------------------------------------- -c -c ****** Return the number of arguments specified on the command -c ****** line. -c -c----------------------------------------------------------------------- -c - use parse_args -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: n -c -c----------------------------------------------------------------------- -c - n=nargs_spec -c - return - end -c####################################################################### - function match_keyword (keyword) -c -c----------------------------------------------------------------------- -c -c ****** Match KEYWORD to the current list of defined keywords. -c -c ****** A successful match returns the index of the first matched -c ****** entry in TABLE; an unsuccessful match returns 0. -c -c----------------------------------------------------------------------- -c - use parse_args - use get_str_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*), intent(in) :: keyword - integer :: match_keyword -c -c----------------------------------------------------------------------- -c - integer :: i -c -c----------------------------------------------------------------------- -c - match_keyword=0 -c - do i=1,nargs - if (keyword.eq.get_str(args(i)%keyword)) then - match_keyword=i - return - end if - end do -c - return - end -c####################################################################### - subroutine new_par (par) -c -c----------------------------------------------------------------------- -c -c ****** Initialize paragraph PAR. -c -c----------------------------------------------------------------------- -c - use paragraph_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(paragraph), pointer :: par -c -c----------------------------------------------------------------------- -c - allocate (par) - nullify (par%line%c) - nullify (par%next) -c - return - end -c####################################################################### - subroutine delete_par (par) -c -c----------------------------------------------------------------------- -c -c ****** Delete paragraph PAR and deallocate its storage and that -c ****** of its linked lists. -c -c----------------------------------------------------------------------- -c - use paragraph_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(paragraph), pointer :: par -c -c----------------------------------------------------------------------- -c - type(paragraph), pointer :: current_par,previous_par -c -c----------------------------------------------------------------------- -c - current_par=>par -c - do -c -c ****** Deallocate the line buffer. -c - call delete_str (current_par%line) -c -c ****** Set the pointer to the next line (if it has been defined). -c - if (.not.associated(current_par%next)) exit - previous_par=>current_par - current_par=>current_par%next - deallocate (previous_par) -c - enddo -c - deallocate (current_par) -c - return - end -c####################################################################### - subroutine add_line (line,par) -c -c----------------------------------------------------------------------- -c -c ****** Add LINE to paragraph PAR. -c -c ****** On exit from this routine, PAR points to a new line, -c ****** and can be used to store the next line of text. -c -c----------------------------------------------------------------------- -c - use paragraph_def - use new_par_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: line - type(paragraph), pointer :: par -c -c----------------------------------------------------------------------- -c -c ****** Store LINE into the string buffer for the current line. -c - call put_str (line,par%line) -c -c ****** Allocate a pointer to the next line. -c - call new_par (par%next) -c -c ****** Set PAR to point to the next line. -c - par=>par%next -c - return - end -c####################################################################### - subroutine print_par (par) -c -c----------------------------------------------------------------------- -c -c ****** Print all lines of paragraph PAR to STDOUT. -c -c----------------------------------------------------------------------- -c - use paragraph_def - use get_str_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(paragraph), pointer :: par -c -c----------------------------------------------------------------------- -c - type(paragraph), pointer :: current_par -c -c----------------------------------------------------------------------- -c - current_par=>par -c - do -c -c ****** Print the line if it has been defined. -c - if (associated(current_par%line%c)) then - write (*,*) trim(get_str(current_par%line)) - end if -c -c ****** Set the pointer to the next line (if it has been defined). -c - if (.not.associated(current_par%next)) exit - current_par=>current_par%next -c - enddo -c - return - end -c####################################################################### - subroutine put_str (cval,str) -c -c----------------------------------------------------------------------- -c -c ****** Store character variable CVAL into string STR. -c ****** This routine allocates storage for the string. -c -c----------------------------------------------------------------------- -c - use string_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: cval - type(string) :: str -c -c----------------------------------------------------------------------- -c - integer :: l,i -c -c----------------------------------------------------------------------- -c - l=len(cval) -c - allocate (str%c(l)) -c - do i=1,l - str%c(i)=cval(i:i) - enddo -c - return - end -c####################################################################### - function get_str (str) -c -c----------------------------------------------------------------------- -c -c ****** Return the value of string STR as the function value -c ****** (as an assumed-length character variable). -c -c----------------------------------------------------------------------- -c - use string_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(string) :: str - character(size(str%c)) :: get_str -c -c----------------------------------------------------------------------- -c - integer :: i -c -c----------------------------------------------------------------------- -c - do i=1,size(str%c) - get_str(i:i)=str%c(i) - enddo -c - return - end -c####################################################################### - subroutine delete_str (str) -c -c----------------------------------------------------------------------- -c -c ****** Delete the storage for string STR. -c -c----------------------------------------------------------------------- -c - use string_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(string) :: str -c -c----------------------------------------------------------------------- -c - if (associated(str%c)) then - deallocate (str%c) - end if - nullify (str%c) -c - return - end -c####################################################################### - function intval (avalue,name) -c -c----------------------------------------------------------------------- -c -c ****** Get the value of the integer in character variable AVALUE. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: avalue - character(*) :: name - integer :: intval -c -c----------------------------------------------------------------------- -c - logical, external :: ifint -c -c----------------------------------------------------------------------- -c - integer :: ivalue -c -c----------------------------------------------------------------------- -c - if (.not.ifint(trim(avalue),ivalue)) then - write (*,*) - write (*,*) '### ERROR in INTVAL:' - write (*,*) '### Could not interpret an integer '// - & 'while setting: ',trim(name) - write (*,*) 'Invalid format: ',trim(avalue) - call exit (1) - end if -c - intval=ivalue -c - return - end -c####################################################################### - function fpval (avalue,name) -c -c----------------------------------------------------------------------- -c -c ****** Get the value of the floating point number in character -c ****** variable AVALUE. -c -c----------------------------------------------------------------------- -c - use number_types -c -c----------------------------------------------------------------------- -c - character(*) :: avalue - character(*) :: name - real(r_typ) :: fpval -c -c----------------------------------------------------------------------- -c - logical, external :: iffp -c -c----------------------------------------------------------------------- -c - real(r_typ) :: value -c -c----------------------------------------------------------------------- -c - if (.not.iffp(trim(avalue),value)) then - write (*,*) - write (*,*) '### ERROR in FPVAL:' - write (*,*) '### Could not interpret a floating point '// - & 'number while setting: ',trim(name) - write (*,*) 'Invalid format: ',trim(avalue) - call exit (1) - end if -c - fpval=value -c - return - end -c####################################################################### - function iffp (alpha,value) -c -c----------------------------------------------------------------------- -c -c ****** Determine if ALPHA represents a floating point number; -c ****** if so, return its value in VALUE. -c -c----------------------------------------------------------------------- -c -c ****** Set IFFP=.TRUE. if ALPHA contains an alphanumeric -c ****** string with the following format: -c -c ALPHA = '[A][B...B][.][B...B][e[A]B[B...B]]', -c -c ****** where A represents a + or - sign, and B represents a digit -c ****** between 0 and 9, inclusive. -c ****** The exponent may be denoted by a lower or upper case e. -c ****** The mantissa must have at least one digit, and the -c ****** the exponent, if present, must have between 1 and 3 digits. -c ****** Otherwise, set IFFP=.FALSE. -c -c----------------------------------------------------------------------- -c - use number_types -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: alpha - real(r_typ) :: value - logical :: iffp -c -c----------------------------------------------------------------------- -c - integer :: nmant,nexp,k,i,ke - logical :: ifpoint,ifexp - character(7) :: fmt -c -c----------------------------------------------------------------------- -c - iffp=.false. - ifpoint=.false. - ifexp=.false. - nmant=0 - nexp=0 -c - do k=1,len_trim(alpha) - i=iachar(alpha(k:k)) -c -c ****** Check for a sign in the first position. -c - if (k.eq.1.and.(i.eq.43.or.i.eq.45)) cycle -c -c ****** Check for a digit. -c - if (i.ge.48.and.i.le.57) then -c -c ****** Count digits in mantissa and exponent. -c - if (ifexp) then - nexp=nexp+1 - else - nmant=nmant+1 - end if - cycle -c - end if -c -c ****** Check for a decimal point. -c - if (.not.ifpoint.and.i.eq.46) then -c -c ****** Check that we are in the mantissa. -c - if (.not.ifexp) then - ifpoint=.true. - cycle - end if -c - end if -c -c ****** Check for an exponent. -c - if (.not.ifexp.and.(i.eq.101.or.i.eq.69)) then - ifexp=.true. - ke=k - cycle - end if -c -c ****** Check for an exponent sign. -c - if (ifexp.and.k.eq.ke+1.and.(i.eq.43.or.i.eq.45)) cycle -c -c ****** Failed check: fall through here. -c - iffp=.false. -c - return -c - enddo -c -c ****** Final check of validity: check number of digits in -c ****** the mantissa and exponent. -c - if (nmant.ge.1) iffp=.true. - if (ifexp.and.(nexp.lt.1.or.nexp.gt.3)) iffp=.false. -c -c ****** Obtain its numeric value. -c - fmt='(f .0)' - write (fmt(3:4),'(i2.2)') len_trim(alpha) -c - if (iffp) read (alpha,fmt) value -c - return - end -c####################################################################### - function ifint (alpha,ivalue) -c -c----------------------------------------------------------------------- -c -c ****** If ALPHA represents an integer, return IFINT=.true., and -c ****** put its value into IVALUE. -c -c ****** Otherwise, return IFINT=.false.. -c -c----------------------------------------------------------------------- -c -c ****** A valid integer has the format: -c -c ALPHA = '[A]B[B...B]', -c -c ****** where A represents a + or - sign, and B represents a digit -c ****** between 0 and 9, inclusive. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: alpha - integer :: ivalue - logical :: ifint -c -c----------------------------------------------------------------------- -c - integer :: k,i - character(5) :: fmt -c -c----------------------------------------------------------------------- -c - ifint=.false. -c - do k=1,len_trim(alpha) -c - i=iachar(alpha(k:k)) -c -c ****** Check for a sign in the first position. -c - if (k.eq.1.and.(i.eq.43.or.i.eq.45)) cycle -c -c ****** Check for a digit. -c - if (i.ge.48.and.i.le.57) then - ifint=.true. - cycle - end if -c -c ****** Failed check: fall through here. -c - ifint=.false. -c - return -c - enddo -c -c ****** Obtain its numeric value. -c - fmt='(i )' - write (fmt(3:4),'(i2.2)') len_trim(alpha) -c - if (ifint) read (alpha,fmt) ivalue -c - return - end -c####################################################################### - subroutine load_list_of_reals (s,label,n,f) -c -c----------------------------------------------------------------------- -c -c ****** Read N real values from character string S into -c ****** array F(N). The values in S may be either space or -c ****** comma separated. -c -c----------------------------------------------------------------------- -c - use number_types -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: s - character(*) :: label - integer :: n - real(r_typ), dimension(n) :: f - intent(in) :: s,label,n - intent(out) :: f -c -c----------------------------------------------------------------------- -c - integer :: i,i0,i1 - character :: delimiter - character(512) :: list -c -c----------------------------------------------------------------------- -c - real(r_typ), external :: fpval -c -c----------------------------------------------------------------------- -c -c ****** Make a local copy of the string (removing leading spaces). -c - list=adjustl(s) -c -c ****** If any commas are present, use a comma as the delimiter. -c ****** Otherwise, one or more spaces is used as a delimiter. -c ****** In this case, compress multiple spaces into a single space. -c - if (index(list,',').ne.0) then - delimiter=',' - else - delimiter=' ' - call delete_repeated_char (list,' ') - end if -c -c ****** Read the list of N numbers sequentially into F. -c - i0=1 - do i=1,n-1 - i1=scan(list(i0:),delimiter)+i0-2 - f(i)=fpval(adjustl(list(i0:i1)),label) - i0=i1+2 - enddo - f(n)=fpval(adjustl(list(i0:)),label) -c - return - end -c####################################################################### - subroutine load_list_of_ints (s,label,n,j) -c -c----------------------------------------------------------------------- -c -c ****** Read N integer values from character string S into -c ****** array J(N). The values in S may be either space or -c ****** comma separated. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: s - character(*) :: label - integer :: n - integer, dimension(n) :: j - intent(in) :: s,label,n - intent(out) :: j -c -c----------------------------------------------------------------------- -c - integer :: i,i0,i1 - character :: delimiter - character(512) :: list -c -c----------------------------------------------------------------------- -c - integer, external :: intval -c -c----------------------------------------------------------------------- -c -c ****** Make a local copy of the string (removing leading spaces). -c - list=adjustl(s) -c -c ****** If any commas are present, use a comma as the delimiter. -c ****** Otherwise, one or more spaces is used as a delimiter. -c ****** In this case, compress multiple spaces into a single space. -c - if (index(list,',').ne.0) then - delimiter=',' - else - delimiter=' ' - call delete_repeated_char (list,' ') - end if -c -c ****** Read the list of N numbers sequentially into J. -c - i0=1 - do i=1,n-1 - i1=scan(list(i0:),delimiter)+i0-2 - j(i)=intval(adjustl(list(i0:i1)),label) - i0=i1+2 - enddo - j(n)=intval(adjustl(list(i0:)),label) -c - return - end -c####################################################################### - function number_of_names (s) -c -c----------------------------------------------------------------------- -c -c ****** Return the number of names in string S. The values -c ****** can be separated by either spaces or commas. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: s - integer :: number_of_names - intent(in) :: s -c -c----------------------------------------------------------------------- -c - integer :: l,n,i0,i1 - character :: delimiter - character(len(s)) :: list -c -c----------------------------------------------------------------------- -c -c ****** Make a local copy of the string (removing leading spaces). -c - list=adjustl(s) -c -c ****** If any commas are present, use a comma as the delimiter. -c ****** Otherwise, one or more spaces is used as a delimiter. -c ****** In this case, compress multiple spaces into a single space. -c - if (index(list,',').ne.0) then - delimiter=',' - else - delimiter=' ' - call delete_repeated_char (list,' ') - end if -c -c ****** Find the number of names N. -c - l=len_trim(list) -c - n=0 -c - i0=1 - do - i1=scan(list(i0:l),delimiter) - if (i1.eq.0) then - if (.not.(delimiter.eq.' '.and.l.lt.i0)) then - n=n+1 - end if - exit - end if - i1=i1+i0-2 - i0=i1+2 - n=n+1 - enddo -c - number_of_names=n -c - return - end -c####################################################################### - subroutine load_list_of_names (s,n,names) -c -c----------------------------------------------------------------------- -c -c ****** Read up to N names from character string S into -c ****** array NAMES(N). The values can be separated by either -c ****** spaces or commas. -c -c----------------------------------------------------------------------- -c -c ****** This routine is designed to be used in conjunction with -c ****** routine NUMBER_OF_NAMES. A typical use would be: -c -c character(32) :: s -c character(32), dimension(:), allocatable :: names -c integer :: n -c -c n=number_of_names(s) -c allocate (names(n)) -c call load_list_of_names (s,n,names) -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: s - integer :: n - character(*), dimension(n) :: names - intent(in) :: s,n - intent(out) :: names -c -c----------------------------------------------------------------------- -c - integer :: l,i,i0,i1 - character :: delimiter - character(len(s)) :: list -c -c----------------------------------------------------------------------- -c -c ****** Make a local copy of the string (removing leading spaces). -c - list=adjustl(s) -c -c ****** If any commas are present, use a comma as the delimiter. -c ****** Otherwise, one or more spaces is used as a delimiter. -c ****** In this case, compress multiple spaces into a single space. -c - if (index(list,',').ne.0) then - delimiter=',' - else - delimiter=' ' - call delete_repeated_char (list,' ') - end if -c -c ****** Read the names. -c - l=len_trim(list) -c - i=0 -c - i0=1 - do - i1=scan(list(i0:l),delimiter) - if (i1.eq.0) then - if (.not.(delimiter.eq.' '.and.l.lt.i0)) then - i=i+1 - if (i.gt.n) exit - names(i)=adjustl(list(i0:l)) - end if - exit - end if - i=i+1 - if (i.gt.n) exit - i1=i1+i0-2 - names(i)=adjustl(list(i0:i1)) - i0=i1+2 - enddo -c - return - end -c####################################################################### - subroutine delete_repeated_char (s,c) -c -c----------------------------------------------------------------------- -c -c ****** Transform repeated adjoining occurrences of character C -c ****** in string S into single occurrences of C. -c -c ****** The string S is overwritten by the modified string. -c -c ****** Trailing blanks in S are ignored. -c -c----------------------------------------------------------------------- -c -c ****** For example, suppose this routine is called with C='d' and -c ****** S='abcdddeefdhdd'. On return, S will have the value -c ****** 'abcdeefdhd'. -c -c----------------------------------------------------------------------- -c -c ****** This routine uses the FORTRAN90 intrinsic SCAN. -c -c----------------------------------------------------------------------- -c - character(*) :: s - character :: c - intent(in) :: c - intent(inout) :: s -c -c----------------------------------------------------------------------- -c - integer :: i,i0 -c -c----------------------------------------------------------------------- -c - i0=1 - do - i=scan(trim(s(i0:)),c) - if (i.eq.0) exit - i0=i0+i - do - if (s(i0:i0).ne.c) exit - s(i0:)=s(i0+1:) - enddo - enddo -c - return - end diff --git a/src/zm_parse_modules.f b/src/zm_parse_modules.f deleted file mode 100644 index 575cbac..0000000 --- a/src/zm_parse_modules.f +++ /dev/null @@ -1,164 +0,0 @@ -c -c ********************************************************************** -c -c Copyright 2018 Predictive Science Inc. -c -c Licensed under the Apache License, Version 2.0 (the "License"); -c you may not use this file except in compliance with the License. -c You may obtain a copy of the License at -c -c http://www.apache.org/licenses/LICENSE-2.0 -c -c Unless required by applicable law or agreed to in writing, software -c distributed under the License is distributed on an "AS IS" BASIS, -c WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or -c implied. -c See the License for the specific language governing permissions and -c limitations under the License. -c -c ********************************************************************** -c -c####################################################################### - module syntax -c -c----------------------------------------------------------------------- -c ****** Group definitions for parsing command-line arguments. -c----------------------------------------------------------------------- -c -c GROUP 1: -c GROUP 2: -c GROUP 3: -c GROUP 4: -c GROUP 5: GROUP SET (only one must be specified) -c GROUP 6: GROUP SET (only one or none must be specified) -c - integer, parameter :: ngroups=6 -c - integer, parameter :: GROUP_K =1 - integer, parameter :: GROUP_A =2 - integer, parameter :: GROUP_KA =3 - integer, parameter :: GROUP_KAA =4 - integer, parameter :: GROUP_K_ONE_ONLY =5 - integer, parameter :: GROUP_K_ONE_OR_NONE=6 -c - end module -c####################################################################### - module string_def -c -c----------------------------------------------------------------------- -c ****** Define a structure to hold a string. -c----------------------------------------------------------------------- -c - implicit none -c - type :: string - character, dimension(:), pointer :: c - end type -c - end module -c####################################################################### - module paragraph_def -c - use string_def -c - implicit none -c -c----------------------------------------------------------------------- -c ****** Define a structure for a linked list of lines -c ****** (i.e., a paragraph). -c----------------------------------------------------------------------- -c - type :: paragraph - type(string) :: line - type(paragraph), pointer :: next - end type -c -c----------------------------------------------------------------------- -c ****** Define a structure to hold a list of paragraphs. -c----------------------------------------------------------------------- -c - type :: parlist - type(paragraph), pointer :: par - end type -c - end module -c####################################################################### - module lcase_interface - interface - function lcase (s) - character(*) :: s - character(len(s)) :: lcase - end function - end interface - end module -c####################################################################### - module ucase_interface - interface - function ucase (s) - character(*) :: s - character(len(s)) :: ucase - end function - end interface - end module -c####################################################################### - module new_par_interface - interface - subroutine new_par (par) - use paragraph_def - implicit none - type(paragraph), pointer :: par - end subroutine - end interface - end module -c####################################################################### - module delete_par_interface - interface - subroutine delete_par (par) - use paragraph_def - implicit none - type(paragraph), pointer :: par - end subroutine - end interface - end module -c####################################################################### - module add_line_interface - interface - subroutine add_line (line,par) - use paragraph_def - implicit none - character(*) :: line - type(paragraph), pointer :: par - end subroutine - end interface - end module -c####################################################################### - module print_par_interface - interface - subroutine print_par (par) - use paragraph_def - implicit none - type(paragraph), pointer :: par - end subroutine - end interface - end module -c####################################################################### - module get_str_interface - interface - function get_str (str) - use string_def - implicit none - type(string) :: str - character(size(str%c)) :: get_str - end function - end interface - end module -c####################################################################### - module get_usage_line_interface - interface - subroutine get_usage_line (usage) - use paragraph_def - implicit none - type(paragraph), pointer :: usage - end subroutine - end interface - end module diff --git a/src/zm_sds.f b/src/zm_sds.f deleted file mode 100644 index 5e9ab2f..0000000 --- a/src/zm_sds.f +++ /dev/null @@ -1,2440 +0,0 @@ -c -c----------------------------------------------------------------------- -c -c ****** Source to build the SDS library. -c ****** These routines are used by Zoran Mikic's tools. -c -c----------------------------------------------------------------------- -c -c ********************************************************************** -c -c Copyright 2018 Predictive Science Inc. -c -c Licensed under the Apache License, Version 2.0 (the "License"); -c you may not use this file except in compliance with the License. -c You may obtain a copy of the License at -c -c http://www.apache.org/licenses/LICENSE-2.0 -c -c Unless required by applicable law or agreed to in writing, software -c distributed under the License is distributed on an "AS IS" BASIS, -c WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or -c implied. -c See the License for the specific language governing permissions and -c limitations under the License. -c -c ********************************************************************** -c -c 07/29/2003, ZM, Version 1.00: -c -c - Original version of the SDS library. -c This library was put together to facilitate the -c development of ZM's tools. -c It includes the new read/write routines for -c scientific data sets (both text and HDF format). -c The code was cleaned up to use standard FORTRAN90. -c -c 02/20/2004, ZM, Version 1.01: -c -c - Added the ability to specify the format in writing -c floating point numbers in routine WRFP. This is used -c in writing SDS text files using routine WRTXT. -c For 32-bit data, WRTXT specifies 7 digits of precision, -c whereas for 64-bit data, 15 digits of precision are -c used. -c -c 04/02/2005, ZM, Version 1.02: -c -c - Added a call to DFSDclear() in routine WRHDF. This -c initializes the SDS interface to the default state -c for each new file. This is needed to prevent settings -c from previous files from interfering with each other. -c -c 06/16/2006, ZM, Version 1.03: -c -c - Fixed some pointer allocation and manipulation issues in -c the HDF read and write routines to make them obey -c the FORTRAN 90 standard. They were misbehaving on -c the IFORT compiler. -c -c 02/24/2009, ZM, Version 1.04: -c -c - Made a small change to the way an SDS is deallocated. -c - Added a routine to initialize an SDS. This is useful -c when deallocating SDS structures. -c -c 08/30/2016, RC, Version 2.00: -c -c - Added ability to read and write hdf5 files. -c Now, rdhdf and wrhdf will read or write an hdf5 file -c if the given fname ends in ".h5". -c The library now needs to be linked to the hdf5 libraries. -c - Modified rdhdf to be compatible with hdf4 files made using -c the SD API instead of the DFSD API. -c -c 09/05/2016, RC, Version 2.01: -c -c - Fixed problem with hdf5 writes when using 32-bit data. -c -c 05/22/2017, RC, Version 2.02: -c -c - Fixed problem with 1D and 2D hdf5 writes when -c the s%dims() are not set to 1 for the unused dimensions. -c -c 05/03/2019, RC, Version 2.03: -c -c - Bug fix for HDF5 reads. -c -c 05/03/2019, RC, Version 2.0.4: -c -c - Bug fix for HDF5 reads (can read mixed prec scales now). -c - Updated rdh5 and wrh5 to use proper casting syntax. -c -c----------------------------------------------------------------------- -c -c####################################################################### - module sdslib_ident -c - character(*), parameter :: cname='SDSLIB' - character(*), parameter :: cvers='2.0.4' - character(*), parameter :: cdate='03/08/2022' -c - end module -c####################################################################### - module assign_ptr_1d_interface - interface - subroutine assign_ptr_1d (from,to) - use number_types - implicit none - real(r_typ), dimension(:), target :: from - real(r_typ), dimension(:), pointer :: to - end subroutine - end interface - end module -c####################################################################### - module assign_ptr_3d_interface - interface - subroutine assign_ptr_3d (from,to) - use number_types - implicit none - real(r_typ), dimension(:,:,:), target :: from - real(r_typ), dimension(:,:,:), pointer :: to - end subroutine - end interface - end module -c####################################################################### - subroutine assign_ptr_1d (from,to) -c -c----------------------------------------------------------------------- -c - use number_types -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - real(r_typ), dimension(:), target :: from - real(r_typ), dimension(:), pointer :: to -c -c----------------------------------------------------------------------- -c - to=>from -c - return - end subroutine -c####################################################################### - subroutine assign_ptr_3d (from,to) -c -c----------------------------------------------------------------------- -c - use number_types -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - real(r_typ), dimension(:,:,:), target :: from - real(r_typ), dimension(:,:,:), pointer :: to -c -c----------------------------------------------------------------------- -c - to=>from -c - return - end subroutine -c####################################################################### - subroutine init_sds_pointer_status (s) -c -c----------------------------------------------------------------------- -c -c ****** Disassociate all the pointers in the SDS in structure S. -c -c----------------------------------------------------------------------- -c -c ****** This is useful when subsequently querying the association -c ****** status of these pointers (e.g., when deallocating storage). -c -c----------------------------------------------------------------------- -c - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c - nullify (s%f) -c - nullify (s%scales(1)%f) - nullify (s%scales(2)%f) - nullify (s%scales(3)%f) -c - return - end subroutine -c####################################################################### - subroutine deallocate_sds (s) -c -c----------------------------------------------------------------------- -c -c ****** Deallocate the memory used by the SDS in structure S. -c -c----------------------------------------------------------------------- -c - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c - if (associated(s%f)) deallocate (s%f) -c - if (associated(s%scales(1)%f)) deallocate (s%scales(1)%f) - if (associated(s%scales(2)%f)) deallocate (s%scales(2)%f) - if (associated(s%scales(3)%f)) deallocate (s%scales(3)%f) -c - return - end subroutine -c####################################################################### - subroutine rdhdf_1d (fname,scale,nx,f,x,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 1D scientific data set from an HDF file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine RDHDF to read the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx - real(r_typ), dimension(:), pointer :: f - real(r_typ), dimension(:), pointer :: x - integer :: ierr - intent(in) :: fname - intent(out) :: scale,nx,ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Read the data set. -c - call rdhdf (fname,s,ierr) -c - if (ierr.ne.0) return -c -c ****** Check that this is a 1D data set. -c - if (s%ndim.ne.1) then - write (*,*) - write (*,*) '### ERROR in RDHDF_1D:' - write (*,*) '### The HDF file does not contain a 1D data set.' - write (*,*) 'File name: ',trim(fname) - ierr=3 - return - end if -c -c ****** Set the output arguments. -c - nx=s%dims(1) - scale=s%scale - x=>s%scales(1)%f -c - allocate (f(nx)) - f=s%f(:,1,1) - deallocate (s%f) -c - return - end -c####################################################################### - subroutine rdhdf_2d (fname,scale,nx,ny,f,x,y,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 2D scientific data set from an HDF file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine RDHDF to read the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx,ny - real(r_typ), dimension(:,:), pointer :: f - real(r_typ), dimension(:), pointer :: x,y - integer :: ierr - intent(in) :: fname - intent(out) :: scale,nx,ny,ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Read the data set. -c - call rdhdf (fname,s,ierr) -c - if (ierr.ne.0) return -c -c ****** Check that this is a 2D data set. -c - if (s%ndim.ne.2) then - write (*,*) - write (*,*) '### ERROR in RDHDF_2D:' - write (*,*) '### The HDF file does not contain a 2D data set.' - write (*,*) 'File name: ',trim(fname) - ierr=3 - return - end if -c -c ****** Set the output arguments. -c - nx=s%dims(1) - ny=s%dims(2) - scale=s%scale - x=>s%scales(1)%f - y=>s%scales(2)%f -c - allocate (f(nx,ny)) - f(:,:)=s%f(:,:,1) - deallocate (s%f) -c - return - end -c####################################################################### - subroutine rdhdf_3d (fname,scale,nx,ny,nz,f,x,y,z,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 3D scientific data set from an HDF file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine RDHDF to read the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx,ny,nz - real(r_typ), dimension(:,:,:), pointer :: f - real(r_typ), dimension(:), pointer :: x,y,z - integer :: ierr - intent(in) :: fname - intent(out) :: scale,nx,ny,nz,ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Read the data set. -c - call rdhdf (fname,s,ierr) -c - if (ierr.ne.0) return -c -c ****** Check that this is a 3D data set. -c - if (s%ndim.ne.3) then - write (*,*) - write (*,*) '### ERROR in RDHDF_3D:' - write (*,*) '### The HDF file does not contain a 3D data set.' - write (*,*) 'File name: ',trim(fname) - ierr=3 - return - end if -c -c ****** Set the output arguments. -c - nx=s%dims(1) - ny=s%dims(2) - nz=s%dims(3) - scale=s%scale - x=>s%scales(1)%f - y=>s%scales(2)%f - z=>s%scales(3)%f - f=>s%f -c - return - end -c####################################################################### - subroutine wrhdf_1d (fname,scale,nx,f,x,hdf32,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 1D scientific data set to an HDF file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine WRHDF to write the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def - use assign_ptr_1d_interface - use assign_ptr_3d_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx - real(r_typ), dimension(nx,1,1) :: f - real(r_typ), dimension(nx) :: x - logical :: hdf32 - integer :: ierr - intent(in) :: fname,scale,nx,f,x,hdf32 - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Set the structure components. -c - s%ndim=1 - s%dims(1)=nx - s%dims(2)=1 - s%dims(3)=1 - s%scale=scale - s%hdf32=hdf32 - if (scale) then - call assign_ptr_1d (x,s%scales(1)%f) - else - nullify (s%scales(1)%f) - end if - nullify (s%scales(2)%f) - nullify (s%scales(3)%f) - call assign_ptr_3d (f,s%f) -c -c ****** Write the data set. -c - call wrhdf (fname,s,ierr) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRHDF_1D:' - write (*,*) '### Could not write the 1D data set.' - write (*,*) 'File name: ',trim(fname) - return - end if -c - return - end -c####################################################################### - subroutine wrhdf_2d (fname,scale,nx,ny,f,x,y,hdf32,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 2D scientific data set to an HDF file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine WRHDF to write the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def - use assign_ptr_1d_interface - use assign_ptr_3d_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx,ny - real(r_typ), dimension(nx,ny,1) :: f - real(r_typ), dimension(nx) :: x - real(r_typ), dimension(ny) :: y - logical :: hdf32 - integer :: ierr - intent(in) :: fname,scale,nx,ny,f,x,y,hdf32 - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Set the structure components. -c - s%ndim=2 - s%dims(1)=nx - s%dims(2)=ny - s%dims(3)=1 - s%scale=scale - s%hdf32=hdf32 - if (scale) then - call assign_ptr_1d (x,s%scales(1)%f) - call assign_ptr_1d (y,s%scales(2)%f) - else - nullify (s%scales(1)%f) - nullify (s%scales(2)%f) - end if - nullify (s%scales(3)%f) - call assign_ptr_3d (f,s%f) -c -c ****** Write the data set. -c - call wrhdf (fname,s,ierr) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRHDF_2D:' - write (*,*) '### Could not write the 2D data set.' - write (*,*) 'File name: ',trim(fname) - return - end if -c - return - end -c####################################################################### - subroutine wrhdf_3d (fname,scale,nx,ny,nz,f,x,y,z,hdf32,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 3D scientific data set to an HDF file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine WRHDF to write the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def - use assign_ptr_1d_interface - use assign_ptr_3d_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx,ny,nz - real(r_typ), dimension(nx,ny,nz) :: f - real(r_typ), dimension(nx) :: x - real(r_typ), dimension(ny) :: y - real(r_typ), dimension(nz) :: z - logical :: hdf32 - integer :: ierr - intent(in) :: fname,scale,nx,ny,nz,f,x,y,z,hdf32 - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Set the structure components. -c - s%ndim=3 - s%dims(1)=nx - s%dims(2)=ny - s%dims(3)=nz - s%scale=scale - s%hdf32=hdf32 - if (scale) then - call assign_ptr_1d (x,s%scales(1)%f) - call assign_ptr_1d (y,s%scales(2)%f) - call assign_ptr_1d (z,s%scales(3)%f) - else - nullify (s%scales(1)%f) - nullify (s%scales(2)%f) - nullify (s%scales(3)%f) - end if - call assign_ptr_3d (f,s%f) -c -c ****** Write the data set. -c - call wrhdf (fname,s,ierr) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRHDF_3D:' - write (*,*) '### Could not write the 3D data set.' - write (*,*) 'File name: ',trim(fname) - return - end if -c - return - end -c####################################################################### - subroutine rdsds (fmt,fname,s,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a scientific data set from file FNAME into -c ****** SDS structure S. -c -c ****** Use routine RDTXT or RDHDF, depending on the format -c ****** specified by FMT. -c -c----------------------------------------------------------------------- -c - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fmt - character(*) :: fname - type(sds) :: s - integer :: ierr - intent(in) :: fmt,fname - intent(out) :: s,ierr -c -c----------------------------------------------------------------------- -c - if (fmt.eq.'text') then - call rdtxt (fname,s,ierr) - else if (fmt.eq.'hdf'.or.fmt.eq.'h5') then - call rdhdf (fname,s,ierr) - else - write (*,*) - write (*,*) '### ERROR in RDSDS:' - write (*,*) '### Invalid file format specified.' - write (*,*) 'File name: ',trim(fname) - write (*,*) 'Format = ',fmt - ierr=5 - return - end if -c - return - end -c####################################################################### - subroutine wrsds (fmt,fname,s,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a scientific data set from SDS structure S to -c ****** file FNAME. -c -c ****** Use routine WRTXT or WRHDF, depending on the format -c ****** specified by FMT. -c -c----------------------------------------------------------------------- -c - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fmt - character(*) :: fname - type(sds) :: s - integer :: ierr - intent(in) :: fmt,fname,s - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c - if (fmt.eq.'text') then - call wrtxt (fname,s,ierr) - else if (fmt.eq.'hdf'.or.fmt.eq.'h5') then - call wrhdf (fname,s,ierr) - else - write (*,*) - write (*,*) '### ERROR in WRSDS:' - write (*,*) '### Invalid file format specified.' - write (*,*) 'File name: ',trim(fname) - write (*,*) 'Format = ',fmt - ierr=5 - return - end if -c - return - end -c####################################################################### - subroutine rdhdf (fname,s,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 1D, 2D, or 3D scientific data set from an HDF file. -c ****** This routine uses the new SD API instead of the -c ****** outdated DFSD API. -c -c----------------------------------------------------------------------- -c -c ****** This routine allocates the required memory and returns -c ****** pointers to the data and scale arrays. -c -c----------------------------------------------------------------------- -c -c ****** Input arguments: -c -c FNAME : [character(*)] -c HDF data file name to read from. -c -c ****** Output arguments: -c -c S : [structure of type SDS] -c A structure that holds the field, its -c dimensions, and the scales, with the -c components described below. -c -c IERR : [integer] -c IERR=0 is returned if the data set was read -c successfully. Otherwise, IERR is set to a -c nonzero value. -c -c ****** Components of structure S: -c -c NDIM : [integer] -c Number of dimensions found in the data set. -c -c DIMS : [integer, dimension(3)] -c Number of points in the data set dimensions. -c For a 1D data set, DIMS(2)=DIMS(3)=1. -c For a 2D data set, DIMS(3)=1. -c -c SCALE : [logical] -c Flag to indicate the presence of scales (axes) -c in the data set. SCALE=.false. means that scales -c were not found; SCALE=.true. means that scales -c were found. -c -c HDF32 : [logical] -c Flag to indicate the precision of the data set -c read in. HDF32=.true. means that the data is -c 32-bit; HDF32=.false. means that the data is -c 64-bit. -c -c SCALES : [structure of type RP1D, dimension(3)] -c This array holds the pointers to the scales -c when SCALE=.true., and is undefined otherwise. -c -c F : [real, pointer to a rank-3 array] -c This array holds the data set values. -c -c ****** The storage for the arrays pointed to by F, and the -c ****** scales (if present) in structure SCALES, is allocated by -c ****** this routine. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - character, dimension(256) :: sds_name, dim_name - type(sds) :: s - integer :: ierr,i - intent(in) :: fname - intent(out) :: s,ierr -c -c----------------------------------------------------------------------- -c - ierr=0 -c -c----------------------------------------------------------------------- -c -c ****** Read hdf5 file if fname ends in '.h5'. -c - i=index(fname,'.h'); - if (fname(i+1:i+2).eq.'h5') then - call rdh5 (fname,s,ierr) - return - else - print*,"HDF4 has been disabled" - ierr=-1 - end if -c - return - end -c####################################################################### - subroutine rdh5 (fname,s,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 1D, 2D, or 3D scientific data set from an HDF5 file. -c ****** The HDF5 file is currently assumed to contain only one -c ****** dataset (1D,2d,3D), with or without scales, in group "/", -c ****** and has no other data members. -c -c----------------------------------------------------------------------- -c -c ****** This routine allocates the required memory and returns -c ****** pointers to the data and scale arrays. -c -c----------------------------------------------------------------------- -c -c ****** Input arguments: -c -c FNAME : [character(*)] -c HDF5 data file name to read from. -c -c ****** Output arguments: -c -c S : [structure of type SDS] -c A structure that holds the field, its -c dimensions, and the scales, with the -c components described below. -c -c IERR : [integer] -c IERR=0 is returned if the data set was read -c successfully. Otherwise, IERR is set to a -c nonzero value. -c -c ****** Components of structure S: -c -c NDIM : [integer] -c Number of dimensions found in the data set. -c -c DIMS : [integer, dimension(3)] -c Number of points in the data set dimensions. -c For a 1D data set, DIMS(2)=DIMS(3)=1. -c For a 2D data set, DIMS(3)=1. -c -c SCALE : [logical] -c Flag to indicate the presence of scales (axes) -c in the data set. SCALE=.false. means that scales -c were not found; SCALE=.true. means that scales -c were found. -c -c HDF32 : [logical] -c Flag to indicate the precision of the data set -c read in. HDF32=.true. means that the data is -c 32-bit; HDF32=.false. means that the data is -c 64-bit. -c -c SCALES : [structure of type RP1D, dimension(3)] -c This array holds the pointers to the scales -c when SCALE=.true., and is undefined otherwise. -c -c F : [real, pointer to a rank-3 array] -c This array holds the data set values. -c -c ****** The storage for the arrays pointed to by F, and the -c ****** scales (if present) in structure SCALES, is allocated by -c ****** this routine. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def - use hdf5 - use h5ds -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - type(sds) :: s - character(*) :: fname -c -c----------------------------------------------------------------------- -c - integer :: ierr -c -c----------------------------------------------------------------------- -c - integer :: i,j,k,obj_type,n_members -c - integer(HID_T) :: file_id ! File identifier - integer(HID_T) :: dset_id ! Dataset identifier - integer(HID_T) :: dspace_id ! Dataspace identifier - integer(HID_T) :: dim_id ! Dimension identifiers - integer(HID_T) :: datatype_id ! Datatype identifiers -c - integer(SIZE_T) :: prec -c - integer(HSIZE_T),dimension(:), allocatable :: s_dims,maxpts - integer(HSIZE_T),dimension(1) :: s_dims_i -c - real(KIND_REAL_4), dimension(:,:,:), allocatable :: f4 - real(KIND_REAL_4), dimension(:), allocatable :: f4dim - real(KIND_REAL_8), dimension(:,:,:), allocatable :: f8 - real(KIND_REAL_8), dimension(:), allocatable :: f8dim -c - character(512) :: obj_name - character(4), parameter :: cname='RDH5' -c - logical :: is_scale -c -c----------------------------------------------------------------------- -c -c ****** Initialize dimension count and arrays. -c - s%ndim=0 - s%dims(:)=1 -c -c ****** Initialize hdf5 interface. -c - call h5open_f (ierr) -c -c ****** Open hdf5 file. -c - call h5Fopen_f (trim(fname),H5F_ACC_RDONLY_F,file_id,ierr) -c -c ****** Get information about the hdf5 file. -c - call h5Gn_members_f (file_id,"/",n_members,ierr) -c -c ****** Make sure there is (at maximum) one 3D dataset with scales. -c - if (n_members.eq.0.or.n_members.gt.4) then - write (*,*) - write (*,*) '### ERROR in ',cname,':' - write (*,*) '### Input file contains too few/many datasets.' - write (*,*) 'File name: ',trim(fname) - return - endif -c -c ****** Assume the Dataset is in index 0 and get its name. -c - call h5Gget_obj_info_idx_f (file_id,"/",0,obj_name,obj_type,ierr) -c -c ****** Open Dataset. -c - call h5Dopen_f (file_id,trim(obj_name),dset_id,ierr) -c -c ****** Make sure the Dataset is not a scale. -c - call h5DSis_scale_f(dset_id,is_scale,ierr) - if (is_scale) then - write (*,*) - write (*,*) '### ERROR in ',cname,':' - write (*,*) '### Input file Dataset at index 0 is a scale.' - write (*,*) 'File name: ',trim(fname) - return - endif -c -c ****** Get dimensions (need s_dims array for format requirements). -c - call h5Dget_space_f (dset_id,dspace_id,ierr) - call h5Sget_simple_extent_ndims_f (dspace_id,s%ndim,ierr) -c - allocate(s_dims(s%ndim)) -c - allocate(maxpts(s%ndim)) - call h5Sget_simple_extent_dims_f (dspace_id,s_dims,maxpts,ierr) - deallocate(maxpts) -c - do j=1,s%ndim - s%dims(j) = INT(s_dims(j)) - end do -c -c ****** Get the floating-point precision of the data and set flag. -c - call h5Dget_type_f (dset_id,datatype_id,ierr) - call h5Tget_precision_f (datatype_id,prec,ierr) -c - if (prec.eq.32) then - s%hdf32=.true. - elseif (prec.eq.64) then - s%hdf32=.false. - end if -c -c ****** Allocate the memory for the Dataset array in s. -c - allocate (s%f(s%dims(1),s%dims(2),s%dims(3))) -c -c ****** Need to read the file in its own datatype, and then convert -c ****** to datatype of s%f. -c - if (s%hdf32) then - allocate (f4(s%dims(1),s%dims(2),s%dims(3))) - call h5Dread_f (dset_id,datatype_id,f4,s_dims,ierr) - do k=1,s%dims(3) - do j=1,s%dims(2) - do i=1,s%dims(1) - s%f(i,j,k) = REAL(f4(i,j,k),r_typ) - enddo - enddo - enddo - deallocate (f4) - else - allocate (f8(s%dims(1),s%dims(2),s%dims(3))) - call h5Dread_f (dset_id,datatype_id,f8,s_dims,ierr) - s%f(:,:,:)=f8(:,:,:) - deallocate (f8) - end if -c - deallocate(s_dims) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in RDH5:' - write (*,*) '### Error while reading the dataset.' - write (*,*) 'File name: ',trim(fname) - write (*,*) '[Error return (from H5DREAD_F) = ',ierr,']' - ierr=4 - return - end if -c -c ****** Close the hdf5 type descriptor. -c - call h5Tclose_f (datatype_id,ierr) -c -c ****** Check if there might be scales present, if so, read them. -c - if (n_members.gt.1) then -c -c ***** First check that the number of scale datasets match the # dim. -c - if (n_members-1.ne.s%ndim) then - write (*,*) - write (*,*) '### ERROR in RDH5:' - write (*,*) '### # scales does not match # dims.' - write (*,*) 'File name: ',trim(fname) - return - end if -c - s%scale=.true. -c -c ****** Loop through scales, make sure each is a scale, and read them. -c - do i=1,n_members-1 -c -c ****** Get the name of scale dataset. -c - call h5Gget_obj_info_idx_f (file_id,"/",i, - & obj_name,obj_type,ierr) -c -c ****** Open scale dataset. -c - call h5Dopen_f (file_id,trim(obj_name),dim_id,ierr) -c -c ****** Make sure the scale is a scale. -c - call h5DSis_scale_f (dim_id,is_scale,ierr) - if (.not.is_scale) then - write (*,*) - write (*,*) '### ERROR in RDH5:' - write (*,*) '### Scale is not a scale.' - write (*,*) 'File name: ',trim(fname) - return - end if -c -c ****** Get dimension of scale. -c - s_dims_i=s%dims(i) -c -c ****** Allocate scale. -c - allocate (s%scales(i)%f(s_dims_i(1))) -c -c ****** Get the floating-point precision of the scale. -c - call h5Dget_type_f (dim_id,datatype_id,ierr) - call h5Tget_precision_f (datatype_id,prec,ierr) -c -c ****** Read in the scale data. -c - if (prec.eq.32) then - allocate (f4dim(s_dims_i(1))) - call h5Dread_f (dim_id,datatype_id,f4dim,s_dims_i,ierr) - do j=1,s%dims(i) - s%scales(i)%f(j) = REAL(f4dim(j),r_typ) - end do - deallocate (f4dim) - elseif (prec.eq.64) then - allocate (f8dim(s_dims_i(1))) - call h5Dread_f (dim_id,datatype_id,f8dim,s_dims_i,ierr) - s%scales(i)%f(:)=f8dim(:) - deallocate (f8dim) - end if -c -c ****** Close the scale dataset. -c - call h5Dclose_f (dim_id,ierr) -c - enddo -c -c ****** Allocate dummy scales (of length 1) for empty dimensions. -c - do i=s%ndim+1,3 - allocate (s%scales(i)%f(1)) - enddo - else -c -c ****** If scales are not present, allocate dummy -c ****** scales (of length 1) so that the pointers to the scales -c ****** are valid. -c - s%scale=.false. -c - allocate (s%scales(1)%f(1)) - allocate (s%scales(2)%f(1)) - allocate (s%scales(3)%f(1)) - end if -c -c ****** Close the dataset. -c - call h5Dclose_f (dset_id,ierr) -c -c ****** Close the file. -c - call h5Fclose_f (file_id,ierr) -c -c ****** Close FORTRAN interface. -c - call h5close_f (ierr) -c - return - end subroutine -c####################################################################### - subroutine wrhdf (fname,s,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 1D, 2D, or 3D scientific data set to an HDF file. -c -c----------------------------------------------------------------------- -c -c ****** Input arguments: -c -c FNAME : [character(*)] -c HDF data file name to write to. -c -c S : [structure of type SDS] -c A structure that holds the field, its -c dimensions, and the scales, with the -c components described below. -c -c ****** Output arguments: -c -c IERR : [integer] -c IERR=0 is returned if the data set was written -c successfully. Otherwise, IERR is set to a -c nonzero value. -c -c ****** Components of structure S: -c -c NDIM : [integer] -c Number of dimensions in the data set. -c -c DIMS : [integer, dimension(3)] -c Number of points in the data set dimensions. -c Only DIMS(1 .. NDIM) are referenced. -c -c SCALE : [logical] -c Flag to indicate the presence of scales (axes) -c in the data set. SCALE=.false. means that scales -c are not being supplied; SCALE=.true. means that -c scales are being supplied. -c -c HDF32 : [logical] -c Flag to specify the precision of the data to -c be written to the file. Set HDF32=.true. to -c write 32-bit data, and HDF32=.false. to write -c 64-bit data. -c -c SCALES : [structure of type RP1D, dimension(3)] -c This array holds the pointers to the scales -c when SCALE=.true., and is not referenced -c otherwise. -c -c F : [real, pointer to a rank-3 array] -c This array holds the data set values. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - type(sds) :: s - integer :: ierr - intent(in) :: fname,s - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c - integer :: iret,i,n -c -c----------------------------------------------------------------------- -c - ierr=0 -c -c ****** Check the number of dimensions. -c - if (s%ndim.le.0.or.s%ndim.gt.3) then - write (*,*) - write (*,*) '### ERROR in WRHDF:' - write (*,*) '### Could not write the SDS data.' - write (*,*) 'Invalid number of dimensions.' - write (*,*) 'Number of dimensions = ',s%ndim - write (*,*) 'File name: ',trim(fname) - ierr=1 - return - end if -c -c ****** Write hdf5 file if fname ends in '.h5'. -c - i=index(fname,'.h') - if (fname(i+1:i+2).eq.'h5') then - call wrh5 (fname,s,ierr) - else - print*,"HDF4 disabled." - ierr=-1 - end if -c - return - end -c####################################################################### - subroutine wrh5 (fname,s,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 1D, 2D, or 3D scientific data set to an HDF5 file. -c -c----------------------------------------------------------------------- -c -c ****** Input arguments: -c -c FNAME : [character(*)] -c HDF data file name to write to. -c -c S : [structure of type SDS] -c A structure that holds the field, its -c dimensions, and the scales, with the -c components described below. -c -c ****** Output arguments: -c -c IERR : [integer] -c IERR=0 is returned if the data set was written -c successfully. Otherwise, IERR is set to a -c nonzero value. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def - use hdf5 - use h5ds -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - type(sds) :: s - integer :: ierr - intent(in) :: fname,s - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c - character(8) :: dimname - integer :: i,j,k - integer(HID_T) :: file_id ! File identifier - integer(HID_T) :: dset_id ! Dataset identifier - integer(HID_T) :: dspace_id,dspacedim_id ! Dataspace identifiers - integer(HID_T) :: dim_id ! Dimension identifiers - integer(HSIZE_T),dimension(3) :: s_dims - integer(HSIZE_T),dimension(1) :: s_dims_i -c - real(KIND_REAL_4), dimension(:,:,:), allocatable :: f4 - real(KIND_REAL_4), dimension(:), allocatable :: f4dim - real(KIND_REAL_8), dimension(:,:,:), allocatable :: f8 - real(KIND_REAL_8), dimension(:), allocatable :: f8dim -c -c----------------------------------------------------------------------- -c -c ****** HDF5 calls are picky about the integer format for the dims -c ****** so the s%dims need to be converted to HSIZE_T integers. -c -c ****** Also, sometimes calls to wrhdf() for 1D and 2D datasets -c ****** do not have the unused dims(i) set to 1 (unset). -c ****** To avoid needing a function call to implicitly reshape -c ****** f(n), set the dims here. -c - do i=1,3 - if (i.le.s%ndim) then - s_dims(i)=INT(s%dims(i),HSIZE_T) - else - s_dims(i)=1 - endif - end do -c -c ****** Initialize hdf5 interface. -c - call h5open_f (ierr) -c -c ****** Create the file. -c - call h5Fcreate_f (trim(fname),H5F_ACC_TRUNC_F,file_id,ierr) -c -c ****** Create the dataspace. -c - call h5Screate_simple_f (s%ndim,s_dims,dspace_id,ierr) -c -c ****** Create and write the dataset (convert s%f to proper type). -c - if (s%hdf32) then - allocate (f4(s_dims(1),s_dims(2),s_dims(3))) - do k=1,s%dims(3) - do j=1,s%dims(2) - do i=1,s%dims(1) - f4(i,j,k) = REAL(s%f(i,j,k),KIND_REAL_4) - enddo - enddo - enddo - call h5Dcreate_f (file_id,'Data',H5T_NATIVE_REAL, - & dspace_id,dset_id,ierr) - call h5Dwrite_f (dset_id,H5T_NATIVE_REAL,f4,s_dims,ierr) - deallocate (f4) - else - allocate (f8(s_dims(1),s_dims(2),s_dims(3))) - f8(:,:,:)=s%f(:,:,:) - call h5Dcreate_f (file_id,'Data',H5T_NATIVE_DOUBLE, - & dspace_id,dset_id,ierr) - call h5Dwrite_f (dset_id,H5T_NATIVE_DOUBLE,f8,s_dims,ierr) - deallocate (f8) - endif - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRH5:' - write (*,*) '### Could not write the dataset.' - write (*,*) 'File name: ',trim(fname) - write (*,*) '[Error return (from h5Dwrite_f) = ',ierr,']' - ierr=4 - return - end if -c -c ****** Check for scales. If present, add them to the hdf5 dataset. -c - if (s%scale) then - do i=1,s%ndim - if (i.eq.1) then - dimname='dim1' - elseif (i.eq.2) then - dimname='dim2' - elseif (i.eq.3) then - dimname='dim3' - endif - s_dims_i=s_dims(i) - call h5Screate_simple_f(1,s_dims_i,dspacedim_id,ierr) - if (s%hdf32) then - allocate (f4dim(s_dims_i(1))) - do j=1,s%dims(i) - f4dim(j) = REAL(s%scales(i)%f(j),KIND_REAL_4) - end do - call h5Dcreate_f (file_id,dimname,H5T_NATIVE_REAL, - & dspacedim_id,dim_id,ierr) - call h5Dwrite_f (dim_id,H5T_NATIVE_REAL, - & f4dim,s_dims_i,ierr) - deallocate (f4dim) - else - allocate (f8dim(s_dims_i(1))) - f8dim(:)=s%scales(i)%f(:) - call h5Dcreate_f (file_id,dimname,H5T_NATIVE_DOUBLE, - & dspacedim_id,dim_id,ierr) - call h5Dwrite_f (dim_id,H5T_NATIVE_DOUBLE, - & f8dim,s_dims_i,ierr) - deallocate (f8dim) - endif - call h5DSset_scale_f (dim_id,ierr,dimname) - call h5DSattach_scale_f (dset_id,dim_id,i,ierr) - call h5DSset_label_f(dset_id, i, dimname, ierr) - call h5Dclose_f (dim_id,ierr) - call h5Sclose_f (dspacedim_id,ierr) - end do - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRH5:' - write (*,*) '### Could not write the scales.' - write (*,*) 'File name: ',trim(fname) - ierr=5 - return - endif - endif -c -c ****** Close the dataset. -c - call h5Dclose_f (dset_id,ierr) -c -c ****** Close the dataspace. -c - call h5Sclose_f (dspace_id,ierr) -c -c ****** Close the file. -c - call h5Fclose_f (file_id,ierr) -c -c ****** Close the hdf5 interface. -c - call h5close_f (ierr) -c - end subroutine -c####################################################################### - subroutine rdtxt_1d (fname,scale,nx,f,x,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 1D scientific data set from a text file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine RDTXT to read the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx - real(r_typ), dimension(:), pointer :: f - real(r_typ), dimension(:), pointer :: x - integer :: ierr - intent(in) :: fname - intent(out) :: scale,nx,ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Read the data set. -c - call rdtxt (fname,s,ierr) -c - if (ierr.ne.0) return -c -c ****** Check that this is a 1D data set. -c - if (s%ndim.ne.1) then - write (*,*) - write (*,*) '### ERROR in RDTXT_1D:' - write (*,*) '### The test file does not contain a 1D data set.' - write (*,*) 'File name: ',trim(fname) - ierr=3 - return - end if -c -c ****** Set the output arguments. -c - nx=s%dims(1) - scale=s%scale - x=>s%scales(1)%f -c - allocate (f(nx)) - f=s%f(:,1,1) - deallocate (s%f) -c - return - end -c####################################################################### - subroutine rdtxt_2d (fname,scale,nx,ny,f,x,y,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 2D scientific data set from a text file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine RDTXT to read the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx,ny - real(r_typ), dimension(:,:), pointer :: f - real(r_typ), dimension(:), pointer :: x,y - integer :: ierr - intent(in) :: fname - intent(out) :: scale,nx,ny,ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Read the data set. -c - call rdtxt (fname,s,ierr) -c - if (ierr.ne.0) return -c -c ****** Check that this is a 2D data set. -c - if (s%ndim.ne.2) then - write (*,*) - write (*,*) '### ERROR in RDTXT_2D:' - write (*,*) '### The text file does not contain a 2D data set.' - write (*,*) 'File name: ',trim(fname) - ierr=3 - return - end if -c -c ****** Set the output arguments. -c - nx=s%dims(1) - ny=s%dims(2) - scale=s%scale - x=>s%scales(1)%f - y=>s%scales(2)%f -c - allocate (f(nx,ny)) - f=s%f(:,:,1) - deallocate (s%f) -c - return - end -c####################################################################### - subroutine rdtxt_3d (fname,scale,nx,ny,nz,f,x,y,z,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 3D scientific data set from a text file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine RDTXT to read the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx,ny,nz - real(r_typ), dimension(:,:,:), pointer :: f - real(r_typ), dimension(:), pointer :: x,y,z - integer :: ierr - intent(in) :: fname - intent(out) :: scale,nx,ny,nz,ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Read the data set. -c - call rdtxt (fname,s,ierr) -c - if (ierr.ne.0) return -c -c ****** Check that this is a 3D data set. -c - if (s%ndim.ne.3) then - write (*,*) - write (*,*) '### ERROR in RDTXT_3D:' - write (*,*) '### The text file does not contain a 3D data set.' - write (*,*) 'File name: ',trim(fname) - ierr=3 - return - end if -c -c ****** Set the output arguments. -c - nx=s%dims(1) - ny=s%dims(2) - nz=s%dims(3) - scale=s%scale - x=>s%scales(1)%f - y=>s%scales(2)%f - z=>s%scales(3)%f - f=>s%f -c - return - end -c####################################################################### - subroutine wrtxt_1d (fname,scale,nx,f,x,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 1D scientific data set to a text file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine WRTXT to write the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def - use assign_ptr_1d_interface - use assign_ptr_3d_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx - real(r_typ), dimension(nx,1,1) :: f - real(r_typ), dimension(nx) :: x - integer :: ierr - intent(in) :: fname,scale,nx,f,x - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Set the structure components. -c - s%ndim=1 - s%dims(1)=nx - s%dims(2)=1 - s%dims(3)=1 - s%scale=scale - if (scale) then - call assign_ptr_1d (x,s%scales(1)%f) - else - nullify (s%scales(1)%f) - end if - nullify (s%scales(2)%f) - nullify (s%scales(3)%f) - call assign_ptr_3d (f,s%f) -c -c ****** Write the data set. -c - call wrtxt (fname,s,ierr) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRTXT_1D:' - write (*,*) '### Could not write the 1D data set.' - write (*,*) 'File name: ',trim(fname) - return - end if -c - return - end -c####################################################################### - subroutine wrtxt_2d (fname,scale,nx,ny,f,x,y,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 2D scientific data set to a text file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine WRTXT to write the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def - use assign_ptr_1d_interface - use assign_ptr_3d_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx,ny - real(r_typ), dimension(nx,ny,1) :: f - real(r_typ), dimension(nx) :: x - real(r_typ), dimension(ny) :: y - integer :: ierr - intent(in) :: fname,scale,nx,ny,f,x,y - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Set the structure components. -c - s%ndim=2 - s%dims(1)=nx - s%dims(2)=ny - s%dims(3)=1 - s%scale=scale - if (scale) then - call assign_ptr_1d (x,s%scales(1)%f) - call assign_ptr_1d (y,s%scales(2)%f) - else - nullify (s%scales(1)%f) - nullify (s%scales(2)%f) - end if - nullify (s%scales(3)%f) - call assign_ptr_3d (f,s%f) -c -c ****** Write the data set. -c - call wrtxt (fname,s,ierr) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRTXT_2D:' - write (*,*) '### Could not write the 2D data set.' - write (*,*) 'File name: ',trim(fname) - return - end if -c - return - end -c####################################################################### - subroutine wrtxt_3d (fname,scale,nx,ny,nz,f,x,y,z,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 3D scientific data set to a text file. -c -c----------------------------------------------------------------------- -c -c ****** This routine calls routine WRTXT to write the file. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def - use assign_ptr_1d_interface - use assign_ptr_3d_interface -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - logical :: scale - integer :: nx,ny,nz - real(r_typ), dimension(nx,ny,nz) :: f - real(r_typ), dimension(nx) :: x - real(r_typ), dimension(ny) :: y - real(r_typ), dimension(nz) :: z - integer :: ierr - intent(in) :: fname,scale,nx,ny,nz,f,x,y,z - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c -c ****** Declaration for the SDS structure. -c - type(sds) :: s -c -c----------------------------------------------------------------------- -c -c ****** Set the structure components. -c - s%ndim=3 - s%dims(1)=nx - s%dims(2)=ny - s%dims(3)=nz - s%scale=scale - if (scale) then - call assign_ptr_1d (x,s%scales(1)%f) - call assign_ptr_1d (y,s%scales(2)%f) - call assign_ptr_1d (z,s%scales(3)%f) - else - nullify (s%scales(1)%f) - nullify (s%scales(2)%f) - nullify (s%scales(3)%f) - end if - call assign_ptr_3d (f,s%f) -c -c ****** Write the data set. -c - call wrtxt (fname,s,ierr) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRTXT_3D:' - write (*,*) '### Could not write the 3D data set.' - write (*,*) 'File name: ',trim(fname) - return - end if -c - return - end -c####################################################################### - subroutine rdtxt (fname,s,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read a 1D, 2D, or 3D scientific data set from a text file. -c -c----------------------------------------------------------------------- -c -c ****** This routine allocates the required memory and returns -c ****** pointers to the data and scale arrays. -c -c----------------------------------------------------------------------- -c -c ****** Input arguments: -c -c FNAME : [character(*)] -c Text data file name to read from. -c -c ****** Output arguments: -c -c S : [structure of type SDS] -c A structure that holds the field, its -c dimensions, and the scales, with the -c components described below. -c -c IERR : [integer] -c IERR=0 is returned if the data set was read -c successfully. Otherwise, IERR is set to a -c nonzero value. -c -c ****** Components of structure S: -c -c NDIM : [integer] -c Number of dimensions found in the data set. -c -c DIMS : [integer, dimension(3)] -c Number of points in the data set dimensions. -c For a 1D data set, DIMS(2)=DIMS(3)=1. -c For a 2D data set, DIMS(3)=1. -c -c SCALE : [logical] -c Flag to indicate the presence of scales (axes) -c in the data set. SCALE=.false. means that scales -c were not found; SCALE=.true. means that scales -c were found. -c -c HDF32 : [logical] -c This flag is is not relevant to text files; -c it is used for HDF data. -c It is arbitrarily set to HDF32=.false.. -c -c SCALES : [structure of type RP1D, dimension(3)] -c This array holds the pointers to the scales -c when SCALE=.true., and is undefined otherwise. -c -c F : [real, pointer to a rank-3 array] -c This array holds the data set values. -c -c ****** The storage for the arrays pointed to by F, and the -c ****** scales (if present) in structure SCALES, is allocated by -c ****** this routine. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - type(sds) :: s - integer :: ierr - intent(in) :: fname - intent(out) :: s,ierr -c -c----------------------------------------------------------------------- -c - integer, parameter :: mxdim=3 - integer :: ifscale,i,n - integer :: int_tmp(1) -c -c----------------------------------------------------------------------- -c - ierr=0 -c -c ****** Open the file for reading. -c - call ffopen (1,fname,'r',ierr) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in RDTXT:' - write (*,*) '### Could not open the text file.' - write (*,*) 'File name: ',trim(fname) - ierr=1 - return - end if -c -c ****** Get the number of dimensions. -c - call rdint (1,1,int_tmp(1),ierr) - s%ndim=int_tmp(1) -c - if (ierr.ne.0) go to 910 -c - if (s%ndim.lt.1.or.s%ndim.gt.mxdim) then - write (*,*) - write (*,*) '### ERROR in RDTXT:' - write (*,*) '### Invalid number of dimensions in file.' - write (*,*) 'File name: ',trim(fname) - write (*,*) 'Number of dimensions = ',s%ndim - write (*,*) 'Maximum number of dimensions = ',mxdim - ierr=2 - return - end if -c -c ****** Read the dimensions. -c - s%dims(:)=1 -c - do i=1,s%ndim - call rdint (1,1,s%dims(i),ierr) - if (ierr.ne.0) go to 910 - if (s%dims(i).le.0) go to 920 - enddo -c -c ****** Check if the scales are present. -c - call rdint (1,1,int_tmp(1),ierr) - ifscale=int_tmp(1) - if (ierr.ne.0) go to 910 -c - s%scale=ifscale.ne.0 -c -c ****** Allocate memory and read the scales (if present). -c -c ****** If scales are not present, allocate dummy scales -c ****** (of length 1) so that the pointers to the scales -c ****** are valid. -c - if (s%scale) then - do i=1,s%ndim - allocate (s%scales(i)%f(s%dims(i))) - call rdfp (1,s%dims(i),s%scales(i)%f,ierr) - if (ierr.ne.0) go to 910 - enddo - do i=s%ndim+1,3 - allocate (s%scales(i)%f(1)) - enddo - else - allocate (s%scales(1)%f(1)) - allocate (s%scales(2)%f(1)) - allocate (s%scales(3)%f(1)) - end if -c -c ****** Allocate memory for the array. -c - allocate (s%f(s%dims(1),s%dims(2),s%dims(3))) -c -c ****** Read the data array. -c - n=product(s%dims(1:s%ndim)) - call rdfp (1,n,s%f,ierr) - if (ierr.ne.0) go to 910 -c - s%hdf32=.false. -c - close (1) -c - return -c - 910 continue -c - write (*,*) - write (*,*) '### ERROR in RDTXT:' - write (*,*) '### Error while reading text data.' - write (*,*) 'File name: ',trim(fname) - ierr=3 - return -c - 920 continue -c - write (*,*) - write (*,*) '### ERROR in RDTXT:' - write (*,*) '### Invalid value for dimension.' - write (*,*) 'Dimension number = ',i - write (*,*) 'Dimension value = ',s%dims(i) - ierr=4 -c - return - end -c####################################################################### - subroutine rdint (iun,n,i,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read N words of INTEGER data into array I from unit IUN -c ****** using a free format text read. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: iun - integer :: n - integer, dimension(n) :: i - integer :: ierr - intent(in) :: iun,n - intent(out) :: i,ierr -c -c----------------------------------------------------------------------- -c - ierr=0 -c - read (iun,*,err=100,end=100) i -c - return -c - 100 continue -c -c ****** Error in reading the data. -c - ierr=1 -c - return - end -c####################################################################### - subroutine rdfp (iun,n,f,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Read N words of REAL data into array F from unit IUN -c ****** using a free format text read. -c -c----------------------------------------------------------------------- -c - use number_types -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: iun - integer :: n - real(r_typ), dimension(n) :: f - integer :: ierr - intent(in) :: iun,n - intent(out) :: f,ierr -c -c----------------------------------------------------------------------- -c - ierr=0 -c - read (iun,*,err=100,end=100) f -c - return -c - 100 continue -c -c ****** Error in reading the data. -c - ierr=1 -c - return - end -c####################################################################### - subroutine wrtxt (fname,s,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write a 1D, 2D, or 3D scientific data set to a text file. -c -c----------------------------------------------------------------------- -c -c ****** Input arguments: -c -c FNAME : [character(*)] -c Text data file name to write to. -c -c S : [structure of type SDS] -c A structure that holds the field, its -c dimensions, and the scales, with the -c components described below. -c -c ****** Output arguments: -c -c IERR : [integer] -c IERR=0 is returned if the data set was written -c successfully. Otherwise, IERR is set to a -c nonzero value. -c -c ****** Components of structure S: -c -c NDIM : [integer] -c Number of dimensions in the data set. -c -c DIMS : [integer, dimension(3)] -c Number of points in the data set dimensions. -c Only DIMS(1 .. NDIM) are referenced. -c -c SCALE : [logical] -c Flag to indicate the presence of scales (axes) -c in the data set. SCALE=.false. means that scales -c are not being supplied; SCALE=.true. means that -c scales are being supplied. -c -c HDF32 : [logical] -c Flag that indicates the precision of the data. -c This flag is used to determine the format for data -c written to the text file. When HDF32=.TRUE., the -c data is assumed to originate from a 32-bit HDF data -c file, and is written with 7 digits to the text file. -c Otherwise, the data is assumed to originate from a -c 64-bit HDF data file, and is written with 14 digits -c to the text file. -c -c SCALES : [structure of type RP1D, dimension(3)] -c This array holds the pointers to the scales -c when SCALE=.true., and is not referenced -c otherwise. -c -c F : [real, pointer to a rank-3 array] -c This array holds the data set values. -c -c----------------------------------------------------------------------- -c - use number_types - use sds_def -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - character(*) :: fname - type(sds) :: s - integer :: ierr - integer :: int_tmp(1) - intent(in) :: fname,s - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c -c ****** Declarations for temporary variables. -c - integer :: i - integer :: n - character(32) :: fmt -c -c----------------------------------------------------------------------- -c - ierr=0 -c -c ****** Open the file for writing. -c - call ffopen (1,fname,'rw',ierr) -c - if (ierr.ne.0) then - write (*,*) - write (*,*) '### ERROR in WRTXT:' - write (*,*) '### Could not open the text file for writing.' - write (*,*) 'File name: ',trim(fname) - ierr=1 - return - end if -c -c ****** Check the number of dimensions. -c - if (s%ndim.le.0.or.s%ndim.gt.3) then - write (*,*) - write (*,*) '### ERROR in WRTXT:' - write (*,*) '### Could not write the SDS data.' - write (*,*) 'Invalid number of dimensions.' - write (*,*) 'NDIM = ',s%ndim - write (*,*) 'File name: ',trim(fname) - ierr=1 - return - end if -c -c ****** Construct the format string for writing floating point -c ****** numbers to the output file. -c - if (s%hdf32) then - fmt='(5(1x,1pe13.6))' - else - fmt='(3(1x,1pe21.14))' - end if -c -c ****** Write the number of dimensions. -c - int_tmp(1)=s%ndim - call wrint (1,1,int_tmp(1),ierr) - if (ierr.ne.0) go to 900 -c -c ****** Write the dimensions. -c - do i=1,s%ndim - call wrint (1,1,s%dims(i),ierr) - if (ierr.ne.0) go to 900 - enddo -c -c ****** Write the scales. -c - if (s%scale) then - int_tmp(1)=1 - call wrint (1,1,int_tmp(1),ierr) - if (ierr.ne.0) go to 900 - do i=1,s%ndim - call wrfp (1,s%dims(i),s%scales(i)%f,fmt,ierr) - if (ierr.ne.0) go to 900 - enddo - else - int_tmp(1)=0 - call wrint (1,1,int_tmp(1),ierr) - if (ierr.ne.0) go to 900 - end if -c -c ****** Write the array. -c - n=product(s%dims(1:s%ndim)) - call wrfp (1,n,s%f,fmt,ierr) - if (ierr.ne.0) go to 900 -c - close (1) -c - return -c - 900 continue -c - write (*,*) - write (*,*) '### ERROR in WRTXT:' - write (*,*) '### Error in writing data to the text file.' - write (*,*) 'File name: ',trim(fname) - ierr=2 -c - return - end -c####################################################################### - subroutine wrint (iun,n,i,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write N words of INTEGER data from array I to the file -c ****** connected to unit IUN using a free format text write. -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: iun - integer :: n - integer, dimension(n) :: i - integer :: ierr - intent(in) :: iun,n,i - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c - ierr=0 -c - write (iun,*,err=100) i -c - return -c - 100 continue -c -c ****** Error in writing the data. -c - ierr=1 -c - return - end -c####################################################################### - subroutine wrfp (iun,n,f,fmt,ierr) -c -c----------------------------------------------------------------------- -c -c ****** Write N words of REAL data from array F to the file -c ****** connected to unit IUN. -c -c ****** FMT specifies the format string to use. -c -c----------------------------------------------------------------------- -c - use number_types -c -c----------------------------------------------------------------------- -c - implicit none -c -c----------------------------------------------------------------------- -c - integer :: iun - integer :: n - real(r_typ), dimension(n) :: f - character(*) :: fmt - integer :: ierr - intent(in) :: iun,n,f,fmt - intent(out) :: ierr -c -c----------------------------------------------------------------------- -c - ierr=0 -c - write (iun,fmt=fmt,err=100) f -c - return -c - 100 continue -c -c ****** Error in writing the data. -c - ierr=1 -c - return - end diff --git a/src/zm_sds_modules.f b/src/zm_sds_modules.f deleted file mode 100644 index 0eb6258..0000000 --- a/src/zm_sds_modules.f +++ /dev/null @@ -1,150 +0,0 @@ -c -c ********************************************************************** -c -c Copyright 2018 Predictive Science Inc. -c -c Licensed under the Apache License, Version 2.0 (the "License"); -c you may not use this file except in compliance with the License. -c You may obtain a copy of the License at -c -c http://www.apache.org/licenses/LICENSE-2.0 -c -c Unless required by applicable law or agreed to in writing, software -c distributed under the License is distributed on an "AS IS" BASIS, -c WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or -c implied. -c See the License for the specific language governing permissions and -c limitations under the License. -c -c ********************************************************************** -c -c####################################################################### - module rp1d_def -c -c----------------------------------------------------------------------- -c ****** Define a structure to hold a REAL 1D pointer. -c----------------------------------------------------------------------- -c - use number_types -c - implicit none -c - type :: rp1d - real(r_typ), dimension(:), pointer :: f - end type -c - end module -c####################################################################### - module sds_def -c -c----------------------------------------------------------------------- -c ****** Definition of the SDS data structure. -c----------------------------------------------------------------------- -c - use number_types - use rp1d_def -c - implicit none -c - integer, parameter, private :: mxdim=3 -c - type :: sds - integer :: ndim - integer, dimension(mxdim) :: dims - logical :: scale - logical :: hdf32 - type(rp1d), dimension(mxdim) :: scales - real(r_typ), dimension(:,:,:), pointer :: f - end type -c - end module -c####################################################################### - module rdhdf_1d_interface - interface - subroutine rdhdf_1d (fname,scale,nx,f,x,ierr) - use number_types - implicit none - character(*) :: fname - logical :: scale - integer :: nx - real(r_typ), dimension(:), pointer :: f - real(r_typ), dimension(:), pointer :: x - integer :: ierr - end subroutine - end interface - end module -c####################################################################### - module rdhdf_2d_interface - interface - subroutine rdhdf_2d (fname,scale,nx,ny,f,x,y,ierr) - use number_types - implicit none - character(*) :: fname - logical :: scale - integer :: nx,ny - real(r_typ), dimension(:,:), pointer :: f - real(r_typ), dimension(:), pointer :: x,y - integer :: ierr - end subroutine - end interface - end module -c####################################################################### - module rdhdf_3d_interface - interface - subroutine rdhdf_3d (fname,scale,nx,ny,nz,f,x,y,z,ierr) - use number_types - implicit none - character(*) :: fname - logical :: scale - integer :: nx,ny,nz - real(r_typ), dimension(:,:,:), pointer :: f - real(r_typ), dimension(:), pointer :: x,y,z - integer :: ierr - end subroutine - end interface - end module -c####################################################################### - module rdtxt_1d_interface - interface - subroutine rdtxt_1d (fname,scale,nx,f,x,ierr) - use number_types - implicit none - character(*) :: fname - logical :: scale - integer :: nx - real(r_typ), dimension(:), pointer :: f - real(r_typ), dimension(:), pointer :: x - integer :: ierr - end subroutine - end interface - end module -c####################################################################### - module rdtxt_2d_interface - interface - subroutine rdtxt_2d (fname,scale,nx,ny,f,x,y,ierr) - use number_types - implicit none - character(*) :: fname - logical :: scale - integer :: nx,ny - real(r_typ), dimension(:,:), pointer :: f - real(r_typ), dimension(:), pointer :: x,y - integer :: ierr - end subroutine - end interface - end module -c####################################################################### - module rdtxt_3d_interface - interface - subroutine rdtxt_3d (fname,scale,nx,ny,nz,f,x,y,z,ierr) - use number_types - implicit none - character(*) :: fname - logical :: scale - integer :: nx,ny,nz - real(r_typ), dimension(:,:,:), pointer :: f - real(r_typ), dimension(:), pointer :: x,y,z - integer :: ierr - end subroutine - end interface - end module