From 883a2d870417d6ecbe1baf1f9a47c6376b06f9c6 Mon Sep 17 00:00:00 2001
From: Edward Hartnett <Edward.Hartnett@noaa.gov>
Date: Thu, 25 Feb 2021 10:22:14 -0700
Subject: [PATCH 01/10] more doxygen for fvcom_tools.fd

---
 docs/Doxyfile.in                        |  2 +-
 sorc/fvcom_tools.fd/module_ncio.f90     |  5 +++--
 sorc/fvcom_tools.fd/module_nwp.f90      | 14 ++++++++------
 sorc/fvcom_tools.fd/module_nwp_base.f90 | 13 ++++++++-----
 4 files changed, 20 insertions(+), 14 deletions(-)

diff --git a/docs/Doxyfile.in b/docs/Doxyfile.in
index d596f9b16..512469f14 100644
--- a/docs/Doxyfile.in
+++ b/docs/Doxyfile.in
@@ -780,7 +780,7 @@ RECURSIVE              = YES
 # Note that relative paths are relative to the directory from which doxygen is
 # run.
 
-EXCLUDE                =
+EXCLUDE                = @abs_top_srcdir@/sorc/fre-nctools.fd  @abs_top_srcdir@/sorc/nst_tf_chg.fd
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded
diff --git a/sorc/fvcom_tools.fd/module_ncio.f90 b/sorc/fvcom_tools.fd/module_ncio.f90
index 1a790aaa9..d3d1093c9 100644
--- a/sorc/fvcom_tools.fd/module_ncio.f90
+++ b/sorc/fvcom_tools.fd/module_ncio.f90
@@ -1,8 +1,9 @@
 !> @file
 !! @brief Functions to read and write netcdf files.
-!!
 !! @author Ming Hu @date 2017-11-01
-!!
+
+!> @brief Functions to read and write netcdf files.
+!! @author Ming Hu @date 2017-11-01
 module module_ncio
 
   use netcdf
diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90
index e8a3d6897..69c586896 100644
--- a/sorc/fvcom_tools.fd/module_nwp.f90
+++ b/sorc/fvcom_tools.fd/module_nwp.f90
@@ -1,11 +1,13 @@
 !> @file
+!! @brief Defines FV3LAM and FVCOM forecast data structure.
+!! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
+
+!> This module defines FV3LAM and FVCOM forecast data structure and
+!! the method to read and write observations from and to those data
+!! structures. It is used by ingest_FVCOM.f90.
 !!
-!! This module defines FV3LAM and FVCOM forecast data structure and the method to
-!! read and write observations from and to those data structures.  It is used by
-!! ingest_FVCOM.f90.
-!!
-!! This script is strongly based upon Eric James' (ESRL/GSL) work with HRRR/WRF
-!!  to get FVCOM data into the model.
+!! This script is strongly based upon Eric James' (ESRL/GSL) work with
+!! HRRR/WRF to get FVCOM data into the model.
 !!
 !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
 !!
diff --git a/sorc/fvcom_tools.fd/module_nwp_base.f90 b/sorc/fvcom_tools.fd/module_nwp_base.f90
index 0f57bab64..f9526e928 100644
--- a/sorc/fvcom_tools.fd/module_nwp_base.f90
+++ b/sorc/fvcom_tools.fd/module_nwp_base.f90
@@ -1,10 +1,13 @@
 !> @file
+!! @brief Defines nwp observation data structure. 
+!! @author David Wright, University of Michigan and GLERL @date 17 Aug 2020
+
+!> This module defines nwp observation data structure and the method
+!! to read and write observations from and to those data
+!! structures. It is used by ingest_FVCOM.f90.
 !!
-!! This module defines nwp observation data structure and the method to
-!! read and write observations from and to those data structures. It is used by
-!! ingest_FVCOM.f90.
-!!
-!! This script is strongly based upon Eric James' (ESRL/GSL) work with HRRR/WRF.
+!! This script is strongly based upon Eric James' (ESRL/GSL) work with
+!! HRRR/WRF.
 !!
 !! @author David Wright, University of Michigan and GLERL @date 17 Aug 2020
 !!

From f8d2f0dd9902d78f1315d6b27191d851ba5415a1 Mon Sep 17 00:00:00 2001
From: Edward Hartnett <Edward.Hartnett@noaa.gov>
Date: Thu, 25 Feb 2021 10:37:33 -0700
Subject: [PATCH 02/10] more doxygen for fvcom

---
 sorc/fvcom_tools.fd/kinds.f90           | 10 +++---
 sorc/fvcom_tools.fd/module_ncio.f90     |  2 +-
 sorc/fvcom_tools.fd/module_nwp.f90      | 42 +++++++++++++++++++------
 sorc/fvcom_tools.fd/module_nwp_base.f90 | 24 +++++++++-----
 4 files changed, 57 insertions(+), 21 deletions(-)

diff --git a/sorc/fvcom_tools.fd/kinds.f90 b/sorc/fvcom_tools.fd/kinds.f90
index 37734c31c..915a0d94b 100644
--- a/sorc/fvcom_tools.fd/kinds.f90
+++ b/sorc/fvcom_tools.fd/kinds.f90
@@ -1,9 +1,11 @@
 !> @file
-!! @brief Module to hold specification kinds for variable declaration.
+!! @brief Hold specification kinds for variable declaration.
+!! @author David Wright, University of Michigan
+
+!> Module to hold specification kinds for variable declaration.
 !!
-!! This module is based on (copied from) Paul vanDelst's 
-!! type_kinds module found in the community radiative transfer
-!! model
+!! This module is based on (copied from) Paul vanDelst's type_kinds
+!! module found in the community radiative transfer model
 !!
 !! @author David Wright, University of Michigan
 module kinds
diff --git a/sorc/fvcom_tools.fd/module_ncio.f90 b/sorc/fvcom_tools.fd/module_ncio.f90
index d3d1093c9..940c060fd 100644
--- a/sorc/fvcom_tools.fd/module_ncio.f90
+++ b/sorc/fvcom_tools.fd/module_ncio.f90
@@ -2,7 +2,7 @@
 !! @brief Functions to read and write netcdf files.
 !! @author Ming Hu @date 2017-11-01
 
-!> @brief Functions to read and write netcdf files.
+!> Functions to read and write netcdf files.
 !! @author Ming Hu @date 2017-11-01
 module module_ncio
 
diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90
index 69c586896..21d46aa9a 100644
--- a/sorc/fvcom_tools.fd/module_nwp.f90
+++ b/sorc/fvcom_tools.fd/module_nwp.f90
@@ -57,11 +57,15 @@ module module_nwp
 
    contains
 
+     !> This subroutine defines the number of variables and their
+     !! names for each NWP data type. The indices of the variables are
+     !! also defined for later reference.
+     !!
+     !! @param this
+     !! @param[in] itype either ' FVCOM' or 'FV3LAM'.
+     !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
       subroutine initial_nwp(this,itype)
 
-!        This subroutine defines the number of variables and their names for
-!        each NWP data type.  The indices of the variables are
-!        also defined for later reference.
 
          class(fcst_nwp) :: this
 
@@ -141,11 +145,13 @@ subroutine initial_nwp(this,itype)
 
       end subroutine initial_nwp
 
+      !> This subroutine lists the setup for NWP data that was done by
+      !! the initial_nwp subroutine.
+      !!
+      !! @param this 
+      !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
       subroutine list_initial_nwp(this)
 
-!        This subroutine lists the setup for NWP data that was done by
-!        the initial_nwp subroutine.
-
          class(fcst_nwp) :: this
 
          integer :: k
@@ -165,11 +171,25 @@ subroutine list_initial_nwp(this)
 
       end subroutine list_initial_nwp
 
+      !> This subroutine initializes arrays to receive the NWP data,
+      !! and opens the file and gets the data.
+      !!
+      !! @param this
+      !! @param[in] filename
+      !! @param[in] itype
+      !! @param[inout] numlon
+      !! @param[inout] numlat
+      !! @param[inout] numtimes
+      !! @param[in] time_to_get
+      !! @param[inout] mask
+      !! @param[inout] sst
+      !! @param[inout] ice
+      !! @param[inout] sfcT
+      !! @param[inout] iceT
+      !! 
+      !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
       subroutine read_nwp(this,filename,itype,numlon,numlat,numtimes,time_to_get,mask,sst,ice,sfcT,iceT)
 
-!        This subroutine initializes arrays to receive the NWP data,
-!        and opens the file and gets the data.
-
          class(fcst_nwp) :: this
 
          character(len=5), intent(in) :: itype
@@ -242,6 +262,10 @@ subroutine read_nwp(this,filename,itype,numlon,numlat,numtimes,time_to_get,mask,
 
       end subroutine read_nwp
 
+      !> Finish.
+      !!
+      !! @param this 
+      !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
       subroutine finish_nwp(this)
 
          class(fcst_nwp) :: this
diff --git a/sorc/fvcom_tools.fd/module_nwp_base.f90 b/sorc/fvcom_tools.fd/module_nwp_base.f90
index f9526e928..87275c3da 100644
--- a/sorc/fvcom_tools.fd/module_nwp_base.f90
+++ b/sorc/fvcom_tools.fd/module_nwp_base.f90
@@ -49,9 +49,12 @@ module module_nwp_base
 
    contains
 
+     !> This subroutine lists the contents of a base nwp observation.
+     !!
+     !! @param this the base nwp obervation
+     !! @author David Wright, University of Michigan and GLERL @date 17 Aug 2020     
       subroutine list_obsbase(this)
 
-!     This subroutine lists the contents of a base nwp observation
 
          class(nwpbase) :: this
 
@@ -76,12 +79,16 @@ subroutine list_obsbase(this)
 
       end subroutine list_obsbase
 
+      !> This subroutine allocates memory for base nwp observation
+      !! variables.
+      !!
+      !! @param this the base nwp obervation
+      !! @param[in] numvar number of variables in this ob type
+      !! @param[in] itquality does this observation include quality
+      !! information?
+      !! @author David Wright, University of Michigan and GLERL @date 17 Aug 2020
       subroutine alloc_obsbase(this,numvar,ifquality)
 
-!     This subroutine allocates memory for base nwp observation variables
-!     Input variables:
-!        numvar : number of variables in this ob type
-!        itquality: does this observation include quality information?
 
          class(nwpbase) :: this
 
@@ -106,10 +113,13 @@ subroutine alloc_obsbase(this,numvar,ifquality)
 
       end subroutine alloc_obsbase
 
+      !> This subroutine releases memory associated with nwp
+      !! observations.
+      !!
+      !! @param this the base nwp obervation
+      !! @author David Wright, University of Michigan and GLERL @date 17 Aug 2020
       subroutine destroy_obsbase(this)
 
-!     This subroutine releases memory associated with nwp observations
-
          class(nwpbase) :: this
 
          this%numvar = 0

From ce3c1b9155f96775cdf42b6fe7b6ebcf3c2c803d Mon Sep 17 00:00:00 2001
From: Edward Hartnett <Edward.Hartnett@noaa.gov>
Date: Thu, 25 Feb 2021 10:41:20 -0700
Subject: [PATCH 03/10] more doxygen for fvcom

---
 sorc/fvcom_tools.fd/module_nwp.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90
index 21d46aa9a..43da62a90 100644
--- a/sorc/fvcom_tools.fd/module_nwp.f90
+++ b/sorc/fvcom_tools.fd/module_nwp.f90
@@ -52,7 +52,7 @@ module module_nwp
          procedure :: finish => finish_nwp
    end type fcst_nwp
 
-   type(ncio) :: ncdata
+   type(ncio) :: ncdata !< Wrapper object for netCDF data file.
 !   type(map_util) :: map
 
    contains

From c195be53820702fc6fcf9af7ea9a2a9ab749e984 Mon Sep 17 00:00:00 2001
From: Edward Hartnett <Edward.Hartnett@noaa.gov>
Date: Thu, 25 Feb 2021 10:58:14 -0700
Subject: [PATCH 04/10] more doxygen for fvcom

---
 sorc/fvcom_tools.fd/module_nwp.f90 | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90
index 43da62a90..9bf62c9c8 100644
--- a/sorc/fvcom_tools.fd/module_nwp.f90
+++ b/sorc/fvcom_tools.fd/module_nwp.f90
@@ -9,7 +9,8 @@
 !! This script is strongly based upon Eric James' (ESRL/GSL) work with
 !! HRRR/WRF to get FVCOM data into the model.
 !!
-!! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
+!! @author David Wright, University of Michigan and GLERL,
+!! @date 17 Aug 2020
 !!
 module module_nwp
 
@@ -63,7 +64,7 @@ module module_nwp
      !!
      !! @param this
      !! @param[in] itype either ' FVCOM' or 'FV3LAM'.
-     !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
+     !! @author David Wright, University of Michigan and GLERL
       subroutine initial_nwp(this,itype)
 
 
@@ -149,7 +150,7 @@ end subroutine initial_nwp
       !! the initial_nwp subroutine.
       !!
       !! @param this 
-      !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
+      !! @author David Wright, University of Michigan and GLERL
       subroutine list_initial_nwp(this)
 
          class(fcst_nwp) :: this
@@ -187,7 +188,7 @@ end subroutine list_initial_nwp
       !! @param[inout] sfcT
       !! @param[inout] iceT
       !! 
-      !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
+      !! @author David Wright, University of Michigan and GLERL
       subroutine read_nwp(this,filename,itype,numlon,numlat,numtimes,time_to_get,mask,sst,ice,sfcT,iceT)
 
          class(fcst_nwp) :: this
@@ -265,7 +266,7 @@ end subroutine read_nwp
       !> Finish.
       !!
       !! @param this 
-      !! @author David Wright, University of Michigan and GLERL, @date 17 Aug 2020
+      !! @author David Wright, University of Michigan and GLERL
       subroutine finish_nwp(this)
 
          class(fcst_nwp) :: this

From dfd7fe3b13f38f4fbe5a648dcdb2702b3803fda7 Mon Sep 17 00:00:00 2001
From: David Wright <dmwright526@gmail.com>
Date: Fri, 26 Feb 2021 15:05:09 -0500
Subject: [PATCH 05/10] Update module_ncio.f90

Additional information and clarification added to description of each subroutine
---
 sorc/fvcom_tools.fd/module_ncio.f90 | 274 ++++++++++++++--------------
 1 file changed, 137 insertions(+), 137 deletions(-)

diff --git a/sorc/fvcom_tools.fd/module_ncio.f90 b/sorc/fvcom_tools.fd/module_ncio.f90
index 940c060fd..bc975d0dd 100644
--- a/sorc/fvcom_tools.fd/module_ncio.f90
+++ b/sorc/fvcom_tools.fd/module_ncio.f90
@@ -260,12 +260,12 @@ subroutine get_dim_nc(this,dimname,dimvalue)
 !
 end subroutine get_dim_nc
 
-!> Read in one field.
+!> Replace 1D character type variable
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_char_1d(this,varname,nd1,field)
 
@@ -295,13 +295,13 @@ subroutine replace_var_nc_char_1d(this,varname,nd1,field)
 !
 end subroutine replace_var_nc_char_1d
 
-!> Read in one field.
+!> Replace 2D character type variable
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_char_2d(this,varname,nd1,nd2,field)
   use netcdf
@@ -342,14 +342,14 @@ subroutine replace_var_nc_char_2d(this,varname,nd1,nd2,field)
 !
 end subroutine replace_var_nc_char_2d
 
-!> Read in one field.
+!> Replace 3D character type variable
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] nd3
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] nd3 length of third dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_char_3d(this,varname,nd1,nd2,nd3,field)
   use netcdf
@@ -394,12 +394,12 @@ subroutine replace_var_nc_char_3d(this,varname,nd1,nd2,nd3,field)
 !
 end subroutine replace_var_nc_char_3d
 
-!> Read in one field.
+!> Replace character type variable
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength length of array
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_char(this,varname,ilength,field)
   use netcdf
@@ -501,12 +501,12 @@ subroutine replace_var_nc_char(this,varname,ilength,field)
 end subroutine replace_var_nc_char
 !--- replace_var_nc_char
 
-!> Replace real.
+!> Replace 1D real type variable
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_real_1d(this,varname,nd1,field)
   use netcdf
@@ -535,13 +535,13 @@ subroutine replace_var_nc_real_1d(this,varname,nd1,field)
 !
 end subroutine replace_var_nc_real_1d
 
-!> Replace real.
+!> Replace 2D real type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_real_2d(this,varname,nd1,nd2,field)
   use netcdf
@@ -582,14 +582,14 @@ subroutine replace_var_nc_real_2d(this,varname,nd1,nd2,field)
 !
 end subroutine replace_var_nc_real_2d
 
-!> Replace real.
+!> Replace 3D real type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] nd3
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] nd3 length of third dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_real_3d(this,varname,nd1,nd2,nd3,field)
 !
@@ -639,12 +639,12 @@ subroutine replace_var_nc_real_3d(this,varname,nd1,nd2,nd3,field)
 !
 end subroutine replace_var_nc_real_3d
 
-!> Read in one field.
+!> Replace real type variable
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength length of array
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_real(this,varname,ilength,field)
   use netcdf
@@ -745,12 +745,12 @@ subroutine replace_var_nc_real(this,varname,ilength,field)
 !
 end subroutine replace_var_nc_real
 
-!> Replace double.
+!> Replace 1D double type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_double_1d(this,varname,nd1,field)
 !
@@ -782,13 +782,13 @@ subroutine replace_var_nc_double_1d(this,varname,nd1,field)
 !
 end subroutine replace_var_nc_double_1d
 
-!> Replace double.
+!> Replace 2D double type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_double_2d(this,varname,nd1,nd2,field)
 !
@@ -832,14 +832,14 @@ subroutine replace_var_nc_double_2d(this,varname,nd1,nd2,field)
 !
 end subroutine replace_var_nc_double_2d
 
-!> Replace double.
+!> Replace 3D double type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] nd3
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] nd3 length of third dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_double_3d(this,varname,nd1,nd2,nd3,field)
 !
@@ -890,12 +890,12 @@ subroutine replace_var_nc_double_3d(this,varname,nd1,nd2,nd3,field)
 end subroutine replace_var_nc_double_3d
 !
 
-!> Read in one field.
+!> Replace double type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength size of array
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_double(this,varname,ilength,field)
   use netcdf
@@ -996,12 +996,12 @@ subroutine replace_var_nc_double(this,varname,ilength,field)
 !
 end subroutine replace_var_nc_double
 
-!> Read in one field.
+!> Replace 1D integer type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 lenth of first dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_int_1d(this,varname,nd1,field)
   use netcdf
@@ -1030,13 +1030,13 @@ subroutine replace_var_nc_int_1d(this,varname,nd1,field)
 !
 end subroutine replace_var_nc_int_1d
 
-!> Read in one field.
+!> Replace 2D integer type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_int_2d(this,varname,nd1,nd2,field)
 !
@@ -1080,14 +1080,14 @@ subroutine replace_var_nc_int_2d(this,varname,nd1,nd2,field)
 !
 end subroutine replace_var_nc_int_2d
 
-!> Read in one field.
+!> Replace 3D integer type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] nd3
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] nd3 length of third dimension
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_int_3d(this,varname,nd1,nd2,nd3,field)
   use netcdf
@@ -1134,12 +1134,12 @@ subroutine replace_var_nc_int_3d(this,varname,nd1,nd2,nd3,field)
 !
 end subroutine replace_var_nc_int_3d
 
-!> Read in one field.
+!> Replace integer type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength size of array
+!! @param[in] field replacement field
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine replace_var_nc_int(this,varname,ilength,field)
   use netcdf
@@ -1240,12 +1240,12 @@ subroutine replace_var_nc_int(this,varname,ilength,field)
 !
 end subroutine replace_var_nc_int
 
-!> Read in one field.
+!> Read in 1D double type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 lenth of first dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_double_1d(this,varname,nd1,field)
   use netcdf
@@ -1277,13 +1277,13 @@ subroutine get_var_nc_double_1d(this,varname,nd1,field)
 !
 end subroutine get_var_nc_double_1d
 
-!> Read in one field.
+!> Read in 2D double type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_double_2d(this,varname,nd1,nd2,field)
   use netcdf
@@ -1328,14 +1328,14 @@ subroutine get_var_nc_double_2d(this,varname,nd1,nd2,field)
 !
 end subroutine get_var_nc_double_2d
 
-!> Read in one field.
+!> Read in 3D double type field.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] nd3
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] nd3 length of third dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_double_3d(this,varname,nd1,nd2,nd3,field)
   use netcdf
@@ -1386,12 +1386,12 @@ subroutine get_var_nc_double_3d(this,varname,nd1,nd2,nd3,field)
 !
 end subroutine get_var_nc_double_3d
 
-!> Read in one field.
+!> Read in double type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength size of array
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_double(this,varname,ilength,field)
   use netcdf
@@ -1492,12 +1492,12 @@ subroutine get_var_nc_double(this,varname,ilength,field)
 !
 end subroutine get_var_nc_double
 
-!> Read in one field.
+!> Read in 1D real type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_real_1d(this,varname,nd1,field)
   use netcdf
@@ -1529,13 +1529,13 @@ subroutine get_var_nc_real_1d(this,varname,nd1,field)
 !
 end subroutine get_var_nc_real_1d
 
-!> Read in one field.
+!> Read in 2D real type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_real_2d(this,varname,nd1,nd2,field)
 !
@@ -1583,14 +1583,14 @@ subroutine get_var_nc_real_2d(this,varname,nd1,nd2,field)
 !
 end subroutine get_var_nc_real_2d
 
-!> Read in one field.
+!> Read in 3D real type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] nd3
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] nd3 length of third dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_real_3d(this,varname,nd1,nd2,nd3,field)
   use netcdf
@@ -1641,12 +1641,12 @@ subroutine get_var_nc_real_3d(this,varname,nd1,nd2,nd3,field)
 !
 end subroutine get_var_nc_real_3d
 
-!> Read in one field.
+!> Read in real type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength size of array
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_real(this,varname,ilength,field)
 !
@@ -1750,12 +1750,12 @@ subroutine get_var_nc_real(this,varname,ilength,field)
 !
 end subroutine get_var_nc_real
 
-!> Read in one field.
+!> Read in 1D integer variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_int_1d(this,varname,nd1,field)
   use netcdf
@@ -1787,13 +1787,13 @@ subroutine get_var_nc_int_1d(this,varname,nd1,field)
 !
 end subroutine get_var_nc_int_1d
 
-!> Read in one field.
+!> Read in 2D integer type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_int_2d(this,varname,nd1,nd2,field)
 !
@@ -1841,14 +1841,14 @@ subroutine get_var_nc_int_2d(this,varname,nd1,nd2,field)
 !
 end subroutine get_var_nc_int_2d
 
-!> Read in one field.
+!> Read in 3D integer type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] nd3
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] nd3 length of third dimension
+!! @param[in] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_int_3d(this,varname,nd1,nd2,nd3,field)
   use netcdf
@@ -1899,12 +1899,12 @@ subroutine get_var_nc_int_3d(this,varname,nd1,nd2,nd3,field)
 !
 end subroutine get_var_nc_int_3d
 
-!> Read in one field.
+!> Read in integer type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength size of array
+!! @param[in] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_int(this,varname,ilength,field)
 !
@@ -2008,12 +2008,12 @@ subroutine get_var_nc_int(this,varname,ilength,field)
 !
 end subroutine get_var_nc_int
 
-!> Read in one field.
+!> Read in 1D short type variable
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_short_1d(this,varname,nd1,field)
   use netcdf
@@ -2045,13 +2045,13 @@ subroutine get_var_nc_short_1d(this,varname,nd1,field)
 !
 end subroutine get_var_nc_short_1d
 
-!> Read in one field.
+!> Read in 2D short type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_short_2d(this,varname,nd1,nd2,field)
 !
@@ -2099,12 +2099,12 @@ subroutine get_var_nc_short_2d(this,varname,nd1,nd2,field)
 !
 end subroutine get_var_nc_short_2d
 !
-!> Read in one field.
+!> Read in short type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength size of array
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_short(this,varname,ilength,field)
   use netcdf
@@ -2205,12 +2205,12 @@ subroutine get_var_nc_short(this,varname,ilength,field)
 !
 end subroutine get_var_nc_short
 
-!> Read in one field.
+!> Read in 1D character type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_char_1d(this,varname,nd1,field)
   use netcdf
@@ -2242,13 +2242,13 @@ subroutine get_var_nc_char_1d(this,varname,nd1,field)
 !
 end subroutine get_var_nc_char_1d
 
-!> Read in one field.
+!> Read in 2D character type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_char_2d(this,varname,nd1,nd2,field)
   use netcdf
@@ -2293,14 +2293,14 @@ subroutine get_var_nc_char_2d(this,varname,nd1,nd2,field)
 !
 end subroutine get_var_nc_char_2d
 
-!> Read in one field.
+!> Read in 3D character type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] nd1
-!! @param[in] nd2
-!! @param[in] nd3
-!! @param[in] field
+!! @param[in] nd1 length of first dimension
+!! @param[in] nd2 length of second dimension
+!! @param[in] nd3 length of third dimension
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_char_3d(this,varname,nd1,nd2,nd3,field)
   use netcdf
@@ -2349,12 +2349,12 @@ subroutine get_var_nc_char_3d(this,varname,nd1,nd2,nd3,field)
 !
 end subroutine get_var_nc_char_3d
 !
-!> Read in one field.
+!> Read in character type variable.
 !!
 !! @param[in] this instance of an ncio class
 !! @param[in] varname name of the variable
-!! @param[in] ilength
-!! @param[in] field
+!! @param[in] ilength size of array
+!! @param[out] field output variable
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine get_var_nc_char(this,varname,ilength,field)
 !
@@ -2475,13 +2475,13 @@ subroutine handle_err(this,status)
   end if
 end subroutine handle_err
 
-!> Convert theta T to T.
+!> Convert theta T (Kelvin) to T (deg C).
 !!
 !! @param[in] this instance of an ncio class
-!! @param[in] nx
-!! @param[in] ny
-!! @param[in] ps
-!! @param[in] t2
+!! @param[in] nx number of grid points in x-dir
+!! @param[in] ny number of grid points in y-dir
+!! @param[in] ps Pressure (Pa)
+!! @param[inout] t2 Pot. Temperature (Kelvin)
 !! @author Ming Hu org: GSD/AMB @date 2017-11-01
 subroutine convert_theta2t_2dgrid(this,nx,ny,ps,t2)
          implicit none
@@ -2514,8 +2514,8 @@ end subroutine convert_theta2t_2dgrid
 !! @param[in] dname1 1st dimension name
 !! @param[in] dname2 2nd dimension name
 !! @param[in] dname3 3rd dimension name
-!! @param lname long name output for netcdf variable
-!! @param units units to use in netcdf variable
+!! @param[in] lname long name output for netcdf variable
+!! @param[in] units units to use in netcdf variable
 !!
 !! @author David.M.Wright org: UM/GLERL @date 2020-09-01
 subroutine add_new_var_3d(this,varname,dname1,dname2,dname3,lname,units)

From e1777fe22cff9366d19aaead5fe578ceefa38cf2 Mon Sep 17 00:00:00 2001
From: David Wright <dmwright526@gmail.com>
Date: Fri, 26 Feb 2021 15:10:29 -0500
Subject: [PATCH 06/10] Update module_nwp.f90

Additional descriptions added for documentation
---
 sorc/fvcom_tools.fd/module_nwp.f90 | 34 +++++++++++++++---------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90
index 9bf62c9c8..c3f7682ed 100644
--- a/sorc/fvcom_tools.fd/module_nwp.f90
+++ b/sorc/fvcom_tools.fd/module_nwp.f90
@@ -62,7 +62,7 @@ module module_nwp
      !! names for each NWP data type. The indices of the variables are
      !! also defined for later reference.
      !!
-     !! @param this
+     !! @param this fcst_nwp object
      !! @param[in] itype either ' FVCOM' or 'FV3LAM'.
      !! @author David Wright, University of Michigan and GLERL
       subroutine initial_nwp(this,itype)
@@ -149,7 +149,7 @@ end subroutine initial_nwp
       !> This subroutine lists the setup for NWP data that was done by
       !! the initial_nwp subroutine.
       !!
-      !! @param this 
+      !! @param this fcst_nwp object
       !! @author David Wright, University of Michigan and GLERL
       subroutine list_initial_nwp(this)
 
@@ -175,19 +175,19 @@ end subroutine list_initial_nwp
       !> This subroutine initializes arrays to receive the NWP data,
       !! and opens the file and gets the data.
       !!
-      !! @param this
-      !! @param[in] filename
-      !! @param[in] itype
-      !! @param[inout] numlon
-      !! @param[inout] numlat
-      !! @param[inout] numtimes
-      !! @param[in] time_to_get
-      !! @param[inout] mask
-      !! @param[inout] sst
-      !! @param[inout] ice
-      !! @param[inout] sfcT
-      !! @param[inout] iceT
-      !! 
+      !! @param this fcst_nwp ojbect
+      !! @param[in] filename netcdf file name
+      !! @param[in] itype either ' FVCOM' or 'FV3LAM'
+      !! @param[inout] numlon number of grid points in x-direction
+      !! @param[inout] numlat number of grid poinst in y-direction
+      !! @param[inout] numtimes length of time dimension
+      !! @param[in] time_to_get integer of time dimension to read in
+      !! @param[inout] mask Water points mask
+      !! @param[inout] sst Water surface temperature
+      !! @param[inout] ice Ice concentration (%)
+      !! @param[inout] sfcT Skin Temperature
+      !! @param[inout] iceT Ice Skin Temperature
+      !!
       !! @author David Wright, University of Michigan and GLERL
       subroutine read_nwp(this,filename,itype,numlon,numlat,numtimes,time_to_get,mask,sst,ice,sfcT,iceT)
 
@@ -263,9 +263,9 @@ subroutine read_nwp(this,filename,itype,numlon,numlat,numtimes,time_to_get,mask,
 
       end subroutine read_nwp
 
-      !> Finish.
+      !> Finish and deallocate.
       !!
-      !! @param this 
+      !! @param this fcst_nwp object
       !! @author David Wright, University of Michigan and GLERL
       subroutine finish_nwp(this)
 

From 8ba983fa7087a88635234c78ecca12dfbb241105 Mon Sep 17 00:00:00 2001
From: Edward Hartnett <Edward.Hartnett@noaa.gov>
Date: Tue, 2 Mar 2021 07:52:51 -0700
Subject: [PATCH 07/10] fixing fvcom_tools doxygen warnings

---
 docs/Doxyfile.in                        |    2 +-
 sorc/fvcom_tools.fd/module_ncio.f90     | 5015 +++++++++++------------
 sorc/fvcom_tools.fd/module_nwp.f90      |   51 +-
 sorc/fvcom_tools.fd/module_nwp_base.f90 |   32 +-
 4 files changed, 2547 insertions(+), 2553 deletions(-)

diff --git a/docs/Doxyfile.in b/docs/Doxyfile.in
index 512469f14..71e7b4008 100644
--- a/docs/Doxyfile.in
+++ b/docs/Doxyfile.in
@@ -1560,7 +1560,7 @@ EXTRA_SEARCH_MAPPINGS  =
 # If the GENERATE_LATEX tag is set to YES doxygen will generate LaTeX output.
 # The default value is: YES.
 
-GENERATE_LATEX         = YES
+GENERATE_LATEX         = NO
 
 # The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. If a
 # relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
diff --git a/sorc/fvcom_tools.fd/module_ncio.f90 b/sorc/fvcom_tools.fd/module_ncio.f90
index bc975d0dd..77b0fc685 100644
--- a/sorc/fvcom_tools.fd/module_ncio.f90
+++ b/sorc/fvcom_tools.fd/module_ncio.f90
@@ -10,2545 +10,2536 @@ module module_ncio
   implicit none
 
   public :: ncio
-! set default to private
+  ! set default to private
   private
-!
+  !
   type :: ncio
-    character(len=256) :: filename
-    integer :: ncid, status
-    integer :: debug_level
-
-    integer :: nDims
-    integer :: ends(4)
-    integer :: xtype
-    character(len=40) :: dimname(4)
-    contains
-      procedure :: open => open_nc
-      procedure :: close => close_nc
-
-! read in dimension from the nc file
-      procedure :: get_dim => get_dim_nc
-
-! read in attribute from the nc file
-      generic   :: get_att => get_att_nc_int,get_att_nc_real,get_att_nc_string
-      procedure :: get_att_nc_int
-      procedure :: get_att_nc_real
-      procedure :: get_att_nc_string
-
-! read in a 1d, 2d, 3d, or 4d field from the nc file
-      generic   :: get_var => get_var_nc_double_1d, get_var_nc_double_2d, &
-                              get_var_nc_double_3d,                       &
-                              get_var_nc_real_1d,get_var_nc_real_2d,      &
-                              get_var_nc_real_3d,                         &
-                              get_var_nc_short_1d,get_var_nc_short_2d,    &
-                              get_var_nc_int_1d,get_var_nc_int_2d,        &
-                              get_var_nc_int_3d,                          &
-                              get_var_nc_char_1d,get_var_nc_char_2d,      &
-                              get_var_nc_char_3d
-      procedure :: get_var_nc_short
-      procedure :: get_var_nc_short_1d
-      procedure :: get_var_nc_short_2d
-      procedure :: get_var_nc_int
-      procedure :: get_var_nc_int_1d
-      procedure :: get_var_nc_int_2d
-      procedure :: get_var_nc_int_3d
-      procedure :: get_var_nc_real
-      procedure :: get_var_nc_real_1d
-      procedure :: get_var_nc_real_2d
-      procedure :: get_var_nc_real_3d
-      procedure :: get_var_nc_double
-      procedure :: get_var_nc_double_1d
-      procedure :: get_var_nc_double_2d
-      procedure :: get_var_nc_double_3d
-      procedure :: get_var_nc_char
-      procedure :: get_var_nc_char_1d
-      procedure :: get_var_nc_char_2d
-      procedure :: get_var_nc_char_3d
-
-! replace 1d, 2d, 3d, or 4d field from the nc file
-      generic   :: replace_var => replace_var_nc_double_1d, replace_var_nc_double_2d, &
-                                  replace_var_nc_double_3d,                           &
-                                  replace_var_nc_real_1d,replace_var_nc_real_2d,      &
-                                  replace_var_nc_real_3d,                             &
-                                  replace_var_nc_int_1d,replace_var_nc_int_2d,        &
-                                  replace_var_nc_int_3d,                              &
-                                  replace_var_nc_char_1d,replace_var_nc_char_2d,      &
-                                  replace_var_nc_char_3d
-      procedure :: replace_var_nc_int
-      procedure :: replace_var_nc_int_1d
-      procedure :: replace_var_nc_int_2d
-      procedure :: replace_var_nc_int_3d
-      procedure :: replace_var_nc_real
-      procedure :: replace_var_nc_real_1d
-      procedure :: replace_var_nc_real_2d
-      procedure :: replace_var_nc_real_3d
-      procedure :: replace_var_nc_double
-      procedure :: replace_var_nc_double_1d
-      procedure :: replace_var_nc_double_2d
-      procedure :: replace_var_nc_double_3d
-      procedure :: replace_var_nc_char
-      procedure :: replace_var_nc_char_1d
-      procedure :: replace_var_nc_char_2d
-      procedure :: replace_var_nc_char_3d
-
-      procedure :: handle_err
-
-      procedure :: convert_theta2t_2dgrid
-!Add a new 3d variable to output file (David.M.Wright)
-      procedure :: add_new_var => add_new_var_3d 
+     character(len=256) :: filename !< Name of data file.
+     integer :: ncid !< File ID.
+     integer :: status !< Return code.
+     integer :: debug_level !< Debug level.
+
+     integer :: nDims !< Number of dims.
+     integer :: ends(4) !< Counts of dims.
+     integer :: xtype !< Type of data.
+     character(len=40) :: dimname(4) !< Name of dims.
+   contains
+     procedure :: open => open_nc !< Open netCDF file. @return
+     procedure :: close => close_nc !< Close netCDF file. @return
+     procedure :: get_dim => get_dim_nc !< read in dimension from the nc file @return
+     generic   :: get_att => get_att_nc_int,get_att_nc_real,get_att_nc_string !< Get attribute. @return
+     procedure :: get_att_nc_int !< Get attribute. @return
+     procedure :: get_att_nc_real !< Get attribute. @return
+     procedure :: get_att_nc_string !< Get attribute. @return
+     generic   :: get_var => get_var_nc_double_1d, get_var_nc_double_2d, & 
+          get_var_nc_double_3d,                       &
+          get_var_nc_real_1d,get_var_nc_real_2d,      &
+          get_var_nc_real_3d,                         &
+          get_var_nc_short_1d,get_var_nc_short_2d,    &
+          get_var_nc_int_1d,get_var_nc_int_2d,        &
+          get_var_nc_int_3d,                          &
+          get_var_nc_char_1d,get_var_nc_char_2d,      &
+          get_var_nc_char_3d      !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_short !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_short_1d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_short_2d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_int !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_int_1d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_int_2d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_int_3d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_real !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_real_1d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_real_2d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_real_3d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_double !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_double_1d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_double_2d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_double_3d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_char !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_char_1d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_char_2d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: get_var_nc_char_3d !< Read in a 1d, 2d, 3d, or 4d field from the nc file. @return
+     generic   :: replace_var => replace_var_nc_double_1d, replace_var_nc_double_2d, &
+          replace_var_nc_double_3d,                           &
+          replace_var_nc_real_1d,replace_var_nc_real_2d,      &
+          replace_var_nc_real_3d,                             &
+          replace_var_nc_int_1d,replace_var_nc_int_2d,        &
+          replace_var_nc_int_3d,                              &
+          replace_var_nc_char_1d,replace_var_nc_char_2d,      &
+          replace_var_nc_char_3d file. !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_int !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_int_1d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_int_2d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_int_3d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_real !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_real_1d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_real_2d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_real_3d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_double !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_double_1d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_double_2d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_double_3d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+     procedure :: replace_var_nc_char !< Replace character type variable. @return
+     procedure :: replace_var_nc_char_1d !< Replace character type variable. @return
+     procedure :: replace_var_nc_char_2d !< Replace character type variable. @return
+     procedure :: replace_var_nc_char_3d !< Replace 3D character type variable. @return
+     procedure :: handle_err !< Handle netCDF errors. @return
+     procedure :: convert_theta2t_2dgrid !< Convert theta T (Kelvin) to T (deg C). @return
+     procedure :: add_new_var => add_new_var_3d      !< Add a new 3d variable to output file. @return
   end type ncio
 
 contains
 
-!> Open a netcdf file, set initial debug level.
-!!
-!! @param this instance of an ncio class
-!! @param filename the file to open
-!! @param action "r" for read, "w" for write
-!! @param debug_level set to non-zero for some verbose output
-!! @author Ming Hu @date 2017-11-01
-subroutine open_nc(this,filename,action,debug_level)
+  !> Open a netcdf file, set initial debug level.
+  !!
+  !! @param this instance of an ncio class
+  !! @param filename the file to open
+  !! @param action "r" for read, "w" for write
+  !! @param debug_level set to non-zero for some verbose output
+  !! @author Ming Hu @date 2017-11-01
+  subroutine open_nc(this,filename,action,debug_level)
+
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: filename
+    character(len=1),intent(in) :: action
+    integer,intent(in),optional :: debug_level
 
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: filename
-  character(len=1),intent(in) :: action
-  integer,intent(in),optional :: debug_level
-
-  integer :: ncid, status
-
-  this%debug_level=20
-  if(present(debug_level)) this%debug_level=debug_level
-
-  this%filename=trim(filename)
-! open existing netCDF dataset
-  if(action=="r" .or. action=="R") then
-     status = nf90_open(path = trim(filename), mode = nf90_nowrite, ncid = ncid)
-  elseif(action=="w" .or. action=="W") then
-     status = nf90_open(path = trim(filename), mode = nf90_write, ncid = ncid)
-  else
-     write(*,*) 'unknow action :', action
-     stop 123
-  endif
-  if (status /= nf90_noerr) call this%handle_err(status)
-  this%ncid=ncid
-
-  if(this%debug_level>0) then
-     write(*,*) '>>> open file: ',trim(this%filename)
-  endif
-
-end subroutine open_nc
-
-!> Close a netcdf file.
-!!
-!! @param this instance of an ncio class
-!! @author Ming Hu org: GSD/AMB @date 2017-04-10
-subroutine close_nc(this)
+    integer :: ncid, status
 
-  implicit none
-!
-  class(ncio) :: this
-
-  integer :: ncid, status
-
-  ncid=this%ncid
-!
-! close netCDF dataset
-  status = nf90_close(ncid)
-  if (status /= nf90_noerr) call this%handle_err(status)
-
-end subroutine close_nc
-
-!> Get attribute in wrf netcdf file.
-!!
-!! @param this instance of an ncio class
-!! @param attname name of the attribute to get
-!! @param rval return value
-!! @author Ming Hu org: GSD/AMB @date 2017-10-04
-subroutine get_att_nc_real(this,attname,rval)
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: attname
-  real, intent(out) :: rval
-
-  integer :: ncid, status
-
-! open existing netCDF dataset
-  ncid=this%ncid
-
-!  get date from exisiting NC file
-  status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), rval)
-  if (status /= nf90_noerr) call this%handle_err(status)
-!
-end subroutine get_att_nc_real
-
-!> Get integer attribute in wrf netcdf file.
-!!
-!! @param this instance of an ncio class
-!! @param attname name of the attribute to get
-!! @param ival value of attribute.
-!! @author Ming Hu org: GSD/AMB @date 2017-10-04
-subroutine get_att_nc_int(this,attname,ival)
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: attname
-  integer, intent(out) :: ival
-
-  integer :: ncid, status
-
-! open existing netCDF dataset
-  ncid=this%ncid
-
-!  get date from exisiting NC file
-  status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), ival)
-  if (status /= nf90_noerr) call this%handle_err(status)
-!
-end subroutine get_att_nc_int
-
-!> Get string attribute in wrf netcdf file.
-!!
-!! @param this instance of an ncio class
-!! @param attname name of the attribute to get
-!! @param string value of attribute.
-!! @author Ming Hu org: GSD/AMB @date 2017-10-04
-subroutine get_att_nc_string(this,attname,string)
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: attname
-  character(len=*), intent(out) :: string
-
-  integer :: ncid, status
-
-! open existing netCDF dataset
-  ncid=this%ncid
-
-!  get date from exisiting NC file
-  status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), string)
-  if (status /= nf90_noerr) call this%handle_err(status)
-!
-end subroutine get_att_nc_string
-
-
-!> Get dimensions in netcdf file.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] dimname name of the dimension
-!! @param[out] dimvalue length of the dimension
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_dim_nc(this,dimname,dimvalue)
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*), intent(in) :: dimname
-  integer,intent(out) :: dimvalue
-
-  integer :: ncid, status
-  integer :: DimId
-
-! open existing netCDF dataset
-  ncid=this%ncid
-
-!  get dimension from exisiting NC file
-  status = nf90_inq_dimid(ncid,trim(dimname), DimId)
-  if (status /= nf90_noerr) call this%handle_err(status)
-  status = nf90_Inquire_Dimension(ncid, DimId, len = dimvalue)
-  if (status /= nf90_noerr) call this%handle_err(status)
-!
-end subroutine get_dim_nc
-
-!> Replace 1D character type variable
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_char_1d(this,varname,nd1,field)
+    this%debug_level=20
+    if(present(debug_level)) this%debug_level=debug_level
+
+    this%filename=trim(filename)
+    ! open existing netCDF dataset
+    if(action=="r" .or. action=="R") then
+       status = nf90_open(path = trim(filename), mode = nf90_nowrite, ncid = ncid)
+    elseif(action=="w" .or. action=="W") then
+       status = nf90_open(path = trim(filename), mode = nf90_write, ncid = ncid)
+    else
+       write(*,*) 'unknow action :', action
+       stop 123
+    endif
+    if (status /= nf90_noerr) call this%handle_err(status)
+    this%ncid=ncid
 
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  character, intent(in) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='replace_var_nc_char_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) (field(i),i=1,min(nd1,10))
-  endif
-
-  call this%replace_var_nc_char(varname,ilength,field)
-!
-end subroutine replace_var_nc_char_1d
-
-!> Replace 2D character type variable
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_char_2d(this,varname,nd1,nd2,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  character, intent(in) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  character,allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='replace_var_nc_char_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  do j=1,nd2
-     istart=(j-1)*nd1+1
-     iend=(j-1)*nd1+nd1
-     temp(istart:iend)=field(:,j)
-  enddo
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) field(1,1)
-  endif
-!
-  call this%replace_var_nc_char(varname,ilength,temp)
-
-  deallocate(temp)
-!
-end subroutine replace_var_nc_char_2d
-
-!> Replace 3D character type variable
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] nd3 length of third dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_char_3d(this,varname,nd1,nd2,nd3,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
-  character, intent(in) :: field(nd1,nd2,nd3) !  values of the field read in
-  integer :: ilength
-!
-  character,allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='replace_var_nc_char_3d'
-!
-  integer :: i,j,k
-  integer :: length2d
-  integer :: istart,iend
-!
-!
-  length2d=nd1*nd2
-  ilength=length2d*nd3
-  allocate(temp(ilength))
-
-  do k=1,nd3
-     do j=1,nd2
-        istart=(k-1)*length2d+(j-1)*nd1+1
-        iend  =(k-1)*length2d+(j-1)*nd1+nd1
-        temp(istart:iend)=field(:,j,k)
-     enddo
-  enddo
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) field(1,1,1)
-  endif
-
-  call this%replace_var_nc_char(varname,ilength,temp)
-
-  deallocate(temp)
-!
-end subroutine replace_var_nc_char_3d
-
-!> Replace character type variable
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength length of array
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_char(this,varname,ilength,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  character, intent(in) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='replace_var_nc_char'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_CHAR) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+    if(this%debug_level>0) then
+       write(*,*) '>>> open file: ',trim(this%filename)
+    endif
+
+  end subroutine open_nc
+
+  !> Close a netcdf file.
+  !!
+  !! @param this instance of an ncio class
+  !! @author Ming Hu org: GSD/AMB @date 2017-04-10
+  subroutine close_nc(this)
+
+    implicit none
+    !
+    class(ncio) :: this
+
+    integer :: ncid, status
+
+    ncid=this%ncid
+    !
+    ! close netCDF dataset
+    status = nf90_close(ncid)
     if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+
+  end subroutine close_nc
+
+  !> Get attribute in wrf netcdf file.
+  !!
+  !! @param this instance of an ncio class
+  !! @param attname name of the attribute to get
+  !! @param rval return value
+  !! @author Ming Hu org: GSD/AMB @date 2017-10-04
+  subroutine get_att_nc_real(this,attname,rval)
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: attname
+    real, intent(out) :: rval
+
+    integer :: ncid, status
+
+    ! open existing netCDF dataset
+    ncid=this%ncid
+
+    !  get date from exisiting NC file
+    status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), rval)
+    if (status /= nf90_noerr) call this%handle_err(status)
+    !
+  end subroutine get_att_nc_real
+
+  !> Get integer attribute in wrf netcdf file.
+  !!
+  !! @param this instance of an ncio class
+  !! @param attname name of the attribute to get
+  !! @param ival value of attribute.
+  !! @author Ming Hu org: GSD/AMB @date 2017-10-04
+  subroutine get_att_nc_int(this,attname,ival)
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: attname
+    integer, intent(out) :: ival
+
+    integer :: ncid, status
+
+    ! open existing netCDF dataset
+    ncid=this%ncid
+
+    !  get date from exisiting NC file
+    status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), ival)
+    if (status /= nf90_noerr) call this%handle_err(status)
+    !
+  end subroutine get_att_nc_int
+
+  !> Get string attribute in wrf netcdf file.
+  !!
+  !! @param this instance of an ncio class
+  !! @param attname name of the attribute to get
+  !! @param string value of attribute.
+  !! @author Ming Hu org: GSD/AMB @date 2017-10-04
+  subroutine get_att_nc_string(this,attname,string)
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: attname
+    character(len=*), intent(out) :: string
+
+    integer :: ncid, status
+
+    ! open existing netCDF dataset
+    ncid=this%ncid
+
+    !  get date from exisiting NC file
+    status = nf90_get_att(ncid, NF90_GLOBAL, trim(attname), string)
+    if (status /= nf90_noerr) call this%handle_err(status)
+    !
+  end subroutine get_att_nc_string
+
+
+  !> Get dimensions in netcdf file.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] dimname name of the dimension
+  !! @param[out] dimvalue length of the dimension
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_dim_nc(this,dimname,dimvalue)
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*), intent(in) :: dimname
+    integer,intent(out) :: dimvalue
+
+    integer :: ncid, status
+    integer :: DimId
+
+    ! open existing netCDF dataset
+    ncid=this%ncid
+
+    !  get dimension from exisiting NC file
+    status = nf90_inq_dimid(ncid,trim(dimname), DimId)
+    if (status /= nf90_noerr) call this%handle_err(status)
+    status = nf90_Inquire_Dimension(ncid, DimId, len = dimvalue)
+    if (status /= nf90_noerr) call this%handle_err(status)
+    !
+  end subroutine get_dim_nc
+
+  !> Replace 1D character type variable
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_char_1d(this,varname,nd1,field)
+
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    character, intent(in) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='replace_var_nc_char_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) (field(i),i=1,min(nd1,10))
+    endif
+
+    call this%replace_var_nc_char(varname,ilength,field)
+    !
+  end subroutine replace_var_nc_char_1d
+
+  !> Replace 2D character type variable
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_char_2d(this,varname,nd1,nd2,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    character, intent(in) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    character,allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='replace_var_nc_char_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    do j=1,nd2
+       istart=(j-1)*nd1+1
+       iend=(j-1)*nd1+nd1
+       temp(istart:iend)=field(:,j)
+    enddo
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) field(1,1)
+    endif
+    !
+    call this%replace_var_nc_char(varname,ilength,temp)
+
+    deallocate(temp)
+    !
+  end subroutine replace_var_nc_char_2d
+
+  !> Replace 3D character type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] nd3 length of third dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_char_3d(this,varname,nd1,nd2,nd3,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
+    character, intent(in) :: field(nd1,nd2,nd3) !  values of the field read in
+    integer :: ilength
+    !
+    character,allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='replace_var_nc_char_3d'
+    !
+    integer :: i,j,k
+    integer :: length2d
+    integer :: istart,iend
+    !
+    !
+    length2d=nd1*nd2
+    ilength=length2d*nd3
+    allocate(temp(ilength))
+
+    do k=1,nd3
+       do j=1,nd2
+          istart=(k-1)*length2d+(j-1)*nd1+1
+          iend  =(k-1)*length2d+(j-1)*nd1+nd1
+          temp(istart:iend)=field(:,j,k)
+       enddo
+    enddo
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) field(1,1,1)
+    endif
+
+    call this%replace_var_nc_char(varname,ilength,temp)
+
+    deallocate(temp)
+    !
+  end subroutine replace_var_nc_char_3d
+
+  !> Replace character type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength length of array
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_char(this,varname,ilength,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    character, intent(in) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='replace_var_nc_char'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_CHAR) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_put_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_put_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>replace variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(8x,a,I10)') 'data type : ',this%xtype
-     write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine replace_var_nc_char
-!--- replace_var_nc_char
-
-!> Replace 1D real type variable
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_real_1d(this,varname,nd1,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  real(4), intent(in) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='replace_var_nc_real_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) (field(i),i=1,min(nd1,10))
-  endif
-!
-  call this%replace_var_nc_real(varname,ilength,field)
-!
-end subroutine replace_var_nc_real_1d
-
-!> Replace 2D real type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_real_2d(this,varname,nd1,nd2,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  real(4), intent(in) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  real(4),allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='replace_var_nc_real_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  do j=1,nd2
-     istart=(j-1)*nd1+1
-     iend=(j-1)*nd1+nd1
-     temp(istart:iend)=field(:,j)
-  enddo
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
-  endif
-
-  call this%replace_var_nc_real(varname,ilength,temp)
-
-  deallocate(temp)
-!
-end subroutine replace_var_nc_real_2d
-
-!> Replace 3D real type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] nd3 length of third dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_real_3d(this,varname,nd1,nd2,nd3,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
-  real(4), intent(in) :: field(nd1,nd2,nd3) !  values of the field read in
-  integer :: ilength
-!
-  real(4),allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='replace_var_nc_real_3d'
-!
-  integer :: i,j,k
-  integer :: length2d
-  integer :: istart,iend
-!
-!
-  length2d=nd1*nd2
-  ilength=length2d*nd3
-  allocate(temp(ilength))
-
-  do k=1,nd3
-     do j=1,nd2
-        istart=(k-1)*length2d+(j-1)*nd1+1
-        iend  =(k-1)*length2d+(j-1)*nd1+nd1
-        temp(istart:iend)=field(:,j,k)
-     enddo
-  enddo
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     do k=1,nd3
-        write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
-     enddo
-  endif
-
-  call this%replace_var_nc_real(varname,ilength,temp)
-
-  deallocate(temp)
-!
-end subroutine replace_var_nc_real_3d
-
-!> Replace real type variable
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength length of array
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_real(this,varname,ilength,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  real(4), intent(in) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='replace_var_nc_real'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_FLOAT) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
-    if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>replace variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(8x,a,I10)') 'data type : ',this%xtype
+       write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine replace_var_nc_char
+  !--- replace_var_nc_char
+
+  !> Replace 1D real type variable
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_real_1d(this,varname,nd1,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    real(4), intent(in) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='replace_var_nc_real_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) (field(i),i=1,min(nd1,10))
+    endif
+    !
+    call this%replace_var_nc_real(varname,ilength,field)
+    !
+  end subroutine replace_var_nc_real_1d
+
+  !> Replace 2D real type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_real_2d(this,varname,nd1,nd2,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    real(4), intent(in) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    real(4),allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='replace_var_nc_real_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    do j=1,nd2
+       istart=(j-1)*nd1+1
+       iend=(j-1)*nd1+nd1
+       temp(istart:iend)=field(:,j)
+    enddo
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
+    endif
+
+    call this%replace_var_nc_real(varname,ilength,temp)
+
+    deallocate(temp)
+    !
+  end subroutine replace_var_nc_real_2d
+
+  !> Replace 3D real type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] nd3 length of third dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_real_3d(this,varname,nd1,nd2,nd3,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
+    real(4), intent(in) :: field(nd1,nd2,nd3) !  values of the field read in
+    integer :: ilength
+    !
+    real(4),allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='replace_var_nc_real_3d'
+    !
+    integer :: i,j,k
+    integer :: length2d
+    integer :: istart,iend
+    !
+    !
+    length2d=nd1*nd2
+    ilength=length2d*nd3
+    allocate(temp(ilength))
+
+    do k=1,nd3
+       do j=1,nd2
+          istart=(k-1)*length2d+(j-1)*nd1+1
+          iend  =(k-1)*length2d+(j-1)*nd1+nd1
+          temp(istart:iend)=field(:,j,k)
+       enddo
+    enddo
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       do k=1,nd3
+          write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
+       enddo
+    endif
+
+    call this%replace_var_nc_real(varname,ilength,temp)
+
+    deallocate(temp)
+    !
+  end subroutine replace_var_nc_real_3d
+
+  !> Replace real type variable
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength length of array
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_real(this,varname,ilength,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    real(4), intent(in) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='replace_var_nc_real'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_FLOAT) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_put_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_put_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>replace variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(8x,a,I10)') 'data type : ',this%xtype
-     write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine replace_var_nc_real
-
-!> Replace 1D double type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_double_1d(this,varname,nd1,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  real(8), intent(in) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='replace_var_nc_double_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) (field(i),i=1,min(nd1,10))
-  endif
-!
-  call this%replace_var_nc_double(varname,ilength,field)
-!
-end subroutine replace_var_nc_double_1d
-
-!> Replace 2D double type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_double_2d(this,varname,nd1,nd2,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  real(8), intent(in) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  real(8),allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='replace_var_nc_double_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  do j=1,nd2
-     istart=(j-1)*nd1+1
-     iend=(j-1)*nd1+nd1
-     temp(istart:iend)=field(:,j)
-  enddo
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
-  endif
-
-  call this%replace_var_nc_double(varname,ilength,temp)
-
-  deallocate(temp)
-!
-end subroutine replace_var_nc_double_2d
-
-!> Replace 3D double type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] nd3 length of third dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_double_3d(this,varname,nd1,nd2,nd3,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
-  real(8), intent(in) :: field(nd1,nd2,nd3) !  values of the field read in
-  integer :: ilength
-!
-  real(8),allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='replace_var_nc_double_3d'
-!
-  integer :: i,j,k
-  integer :: length2d
-  integer :: istart,iend
-!
-!
-  length2d=nd1*nd2
-  ilength=length2d*nd3
-  allocate(temp(ilength))
-
-  do k=1,nd3
-     do j=1,nd2
-        istart=(k-1)*length2d+(j-1)*nd1+1
-        iend  =(k-1)*length2d+(j-1)*nd1+nd1
-        temp(istart:iend)=field(:,j,k)
-     enddo
-  enddo
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     do k=1,nd3
-        write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
-     enddo
-  endif
-
-  call this%replace_var_nc_double(varname,ilength,temp)
-
-  deallocate(temp)
-!
-end subroutine replace_var_nc_double_3d
-!
-
-!> Replace double type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength size of array
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_double(this,varname,ilength,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  real(8), intent(in) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='replace_var_nc_double'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_DOUBLE) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
-    if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>replace variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(8x,a,I10)') 'data type : ',this%xtype
+       write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine replace_var_nc_real
+
+  !> Replace 1D double type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_double_1d(this,varname,nd1,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    real(8), intent(in) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='replace_var_nc_double_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) (field(i),i=1,min(nd1,10))
+    endif
+    !
+    call this%replace_var_nc_double(varname,ilength,field)
+    !
+  end subroutine replace_var_nc_double_1d
+
+  !> Replace 2D double type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_double_2d(this,varname,nd1,nd2,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    real(8), intent(in) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    real(8),allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='replace_var_nc_double_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    do j=1,nd2
+       istart=(j-1)*nd1+1
+       iend=(j-1)*nd1+nd1
+       temp(istart:iend)=field(:,j)
+    enddo
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
+    endif
+
+    call this%replace_var_nc_double(varname,ilength,temp)
+
+    deallocate(temp)
+    !
+  end subroutine replace_var_nc_double_2d
+
+  !> Replace 3D double type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] nd3 length of third dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_double_3d(this,varname,nd1,nd2,nd3,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
+    real(8), intent(in) :: field(nd1,nd2,nd3) !  values of the field read in
+    integer :: ilength
+    !
+    real(8),allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='replace_var_nc_double_3d'
+    !
+    integer :: i,j,k
+    integer :: length2d
+    integer :: istart,iend
+    !
+    !
+    length2d=nd1*nd2
+    ilength=length2d*nd3
+    allocate(temp(ilength))
+
+    do k=1,nd3
+       do j=1,nd2
+          istart=(k-1)*length2d+(j-1)*nd1+1
+          iend  =(k-1)*length2d+(j-1)*nd1+nd1
+          temp(istart:iend)=field(:,j,k)
+       enddo
+    enddo
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       do k=1,nd3
+          write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
+       enddo
+    endif
+
+    call this%replace_var_nc_double(varname,ilength,temp)
+
+    deallocate(temp)
+    !
+  end subroutine replace_var_nc_double_3d
+  !
+
+  !> Replace double type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength size of array
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_double(this,varname,ilength,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    real(8), intent(in) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='replace_var_nc_double'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_DOUBLE) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_put_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_put_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>replace variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(8x,a,I10)') 'data type : ',this%xtype
-     write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine replace_var_nc_double
-
-!> Replace 1D integer type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 lenth of first dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_int_1d(this,varname,nd1,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  integer, intent(in) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='get_var_nc_int_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) (field(i),i=1,min(nd1,10))
-  endif
-
-  call this%replace_var_nc_int(varname,ilength,field)
-!
-end subroutine replace_var_nc_int_1d
-
-!> Replace 2D integer type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_int_2d(this,varname,nd1,nd2,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  integer, intent(in) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  integer,allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='replace_var_nc_int_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  do j=1,nd2
-     istart=(j-1)*nd1+1
-     iend=(j-1)*nd1+nd1
-     temp(istart:iend)=field(:,j)
-  enddo
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
-  endif
-
-  call this%replace_var_nc_int(varname,ilength,temp)
-
-  deallocate(temp)
-!
-end subroutine replace_var_nc_int_2d
-
-!> Replace 3D integer type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] nd3 length of third dimension
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_int_3d(this,varname,nd1,nd2,nd3,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
-  integer, intent(in) :: field(nd1,nd2,nd3) !  values of the field read in
-  integer :: ilength
-!
-  integer,allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='replace_var_nc_int_3d'
-!
-  integer :: i,j,k
-  integer :: length2d
-  integer :: istart,iend
-!
-!
-  length2d=nd1*nd2
-  ilength=length2d*nd3
-  allocate(temp(ilength))
-
-  do k=1,nd3
-     do j=1,nd2
-        istart=(k-1)*length2d+(j-1)*nd1+1
-        iend  =(k-1)*length2d+(j-1)*nd1+nd1
-        temp(istart:iend)=field(:,j,k)
-     enddo
-  enddo
-!
-  if(this%debug_level>100) then
-     write(*,*) trim(thissubname),' show samples:'
-     do k=1,nd3
-        write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
-     enddo
-  endif
-
-  call this%replace_var_nc_int(varname,ilength,temp)
-
-  deallocate(temp)
-!
-end subroutine replace_var_nc_int_3d
-
-!> Replace integer type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength size of array
-!! @param[in] field replacement field
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine replace_var_nc_int(this,varname,ilength,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  integer, intent(in) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='replace_var_nc_int'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_INT) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
-    if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>replace variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(8x,a,I10)') 'data type : ',this%xtype
+       write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine replace_var_nc_double
+
+  !> Replace 1D integer type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 lenth of first dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_int_1d(this,varname,nd1,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    integer, intent(in) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='get_var_nc_int_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) (field(i),i=1,min(nd1,10))
+    endif
+
+    call this%replace_var_nc_int(varname,ilength,field)
+    !
+  end subroutine replace_var_nc_int_1d
+
+  !> Replace 2D integer type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_int_2d(this,varname,nd1,nd2,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    integer, intent(in) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    integer,allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='replace_var_nc_int_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    do j=1,nd2
+       istart=(j-1)*nd1+1
+       iend=(j-1)*nd1+nd1
+       temp(istart:iend)=field(:,j)
+    enddo
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
+    endif
+
+    call this%replace_var_nc_int(varname,ilength,temp)
+
+    deallocate(temp)
+    !
+  end subroutine replace_var_nc_int_2d
+
+  !> Replace 3D integer type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] nd3 length of third dimension
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_int_3d(this,varname,nd1,nd2,nd3,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
+    integer, intent(in) :: field(nd1,nd2,nd3) !  values of the field read in
+    integer :: ilength
+    !
+    integer,allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='replace_var_nc_int_3d'
+    !
+    integer :: i,j,k
+    integer :: length2d
+    integer :: istart,iend
+    !
+    !
+    length2d=nd1*nd2
+    ilength=length2d*nd3
+    allocate(temp(ilength))
+
+    do k=1,nd3
+       do j=1,nd2
+          istart=(k-1)*length2d+(j-1)*nd1+1
+          iend  =(k-1)*length2d+(j-1)*nd1+nd1
+          temp(istart:iend)=field(:,j,k)
+       enddo
+    enddo
+    !
+    if(this%debug_level>100) then
+       write(*,*) trim(thissubname),' show samples:'
+       do k=1,nd3
+          write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
+       enddo
+    endif
+
+    call this%replace_var_nc_int(varname,ilength,temp)
+
+    deallocate(temp)
+    !
+  end subroutine replace_var_nc_int_3d
+
+  !> Replace integer type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength size of array
+  !! @param[in] field replacement field
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine replace_var_nc_int(this,varname,ilength,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    integer, intent(in) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='replace_var_nc_int'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_INT) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_put_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_put_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>replace variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(8x,a,I10)') 'data type : ',this%xtype
-     write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine replace_var_nc_int
-
-!> Read in 1D double type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 lenth of first dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_double_1d(this,varname,nd1,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  real(8), intent(out) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='get_var_nc_double_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-  call this%get_var_nc_double(varname,ilength,field)
-!
-  if(nd1==this%ends(1)) then
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) (field(i),i=1,min(nd1,10))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-  endif
-!
-end subroutine get_var_nc_double_1d
-
-!> Read in 2D double type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_double_2d(this,varname,nd1,nd2,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  real(8), intent(out) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  real(8),allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='get_var_nc_double_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  call this%get_var_nc_double(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
-     do j=1,nd2
-        istart=(j-1)*nd1+1
-        iend=(j-1)*nd1+nd1
-        field(:,j)=temp(istart:iend)
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_double_2d
-
-!> Read in 3D double type field.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] nd3 length of third dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_double_3d(this,varname,nd1,nd2,nd3,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
-  real(8), intent(out) :: field(nd1,nd2,nd3) !  values of the field read in
-  integer :: ilength
-!
-  real(8),allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='get_var_nc_double_3d'
-!
-  integer :: i,j,k
-  integer :: length2d
-  integer :: istart,iend
-!
-!
-  length2d=nd1*nd2
-  ilength=length2d*nd3
-  allocate(temp(ilength))
-
-  call this%get_var_nc_double(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then
-     do k=1,nd3
-        do j=1,nd2
-           istart=(k-1)*length2d+(j-1)*nd1+1
-           iend  =(k-1)*length2d+(j-1)*nd1+nd1
-           field(:,j,k)=temp(istart:iend)
-        enddo
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        do k=1,nd3
-           write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
-        enddo
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_double_3d
-
-!> Read in double type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength size of array
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_double(this,varname,ilength,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  real(8), intent(out) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='get_var_nc_double'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_DOUBLE) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_DOUBLE,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
-    if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>replace variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(8x,a,I10)') 'data type : ',this%xtype
+       write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine replace_var_nc_int
+
+  !> Read in 1D double type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 lenth of first dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_double_1d(this,varname,nd1,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    real(8), intent(out) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='get_var_nc_double_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    call this%get_var_nc_double(varname,ilength,field)
+    !
+    if(nd1==this%ends(1)) then
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) (field(i),i=1,min(nd1,10))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+    endif
+    !
+  end subroutine get_var_nc_double_1d
+
+  !> Read in 2D double type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_double_2d(this,varname,nd1,nd2,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    real(8), intent(out) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    real(8),allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='get_var_nc_double_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    call this%get_var_nc_double(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
+       do j=1,nd2
+          istart=(j-1)*nd1+1
+          iend=(j-1)*nd1+nd1
+          field(:,j)=temp(istart:iend)
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_double_2d
+
+  !> Read in 3D double type field.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] nd3 length of third dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_double_3d(this,varname,nd1,nd2,nd3,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
+    real(8), intent(out) :: field(nd1,nd2,nd3) !  values of the field read in
+    integer :: ilength
+    !
+    real(8),allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='get_var_nc_double_3d'
+    !
+    integer :: i,j,k
+    integer :: length2d
+    integer :: istart,iend
+    !
+    !
+    length2d=nd1*nd2
+    ilength=length2d*nd3
+    allocate(temp(ilength))
+
+    call this%get_var_nc_double(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then
+       do k=1,nd3
+          do j=1,nd2
+             istart=(k-1)*length2d+(j-1)*nd1+1
+             iend  =(k-1)*length2d+(j-1)*nd1+nd1
+             field(:,j,k)=temp(istart:iend)
+          enddo
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          do k=1,nd3
+             write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
+          enddo
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_double_3d
+
+  !> Read in double type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength size of array
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_double(this,varname,ilength,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    real(8), intent(out) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='get_var_nc_double'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_DOUBLE) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_DOUBLE,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_get_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_get_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>read in variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(a,I10)') '        data type : ',this%xtype
-     write(*,'(a,I10)')'        dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(a,I5,I10,2x,a)') '        rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine get_var_nc_double
-
-!> Read in 1D real type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_real_1d(this,varname,nd1,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  real(4), intent(out) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='get_var_nc_real_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-  call this%get_var_nc_real(varname,ilength,field)
-!
-  if(nd1==this%ends(1)) then
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) (field(i),i=1,min(nd1,10))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-  endif
-!
-end subroutine get_var_nc_real_1d
-
-!> Read in 2D real type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_real_2d(this,varname,nd1,nd2,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  real(4), intent(out) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  real(4),allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='get_var_nc_real_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  call this%get_var_nc_real(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
-     do j=1,nd2
-        istart=(j-1)*nd1+1
-        iend=(j-1)*nd1+nd1
-        field(:,j)=temp(istart:iend)
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_real_2d
-
-!> Read in 3D real type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] nd3 length of third dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_real_3d(this,varname,nd1,nd2,nd3,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
-  real(4), intent(out) :: field(nd1,nd2,nd3) !  values of the field read in
-  integer :: ilength
-!
-  real(4),allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='get_var_nc_real_3d'
-!
-  integer :: i,j,k
-  integer :: length2d
-  integer :: istart,iend
-!
-!
-  length2d=nd1*nd2
-  ilength=length2d*nd3
-  allocate(temp(ilength))
-
-  call this%get_var_nc_real(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then
-     do k=1,nd3
-        do j=1,nd2
-           istart=(k-1)*length2d+(j-1)*nd1+1
-           iend  =(k-1)*length2d+(j-1)*nd1+nd1
-           field(:,j,k)=temp(istart:iend)
-        enddo
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        do k=1,nd3
-           write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
-        enddo
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_real_3d
-
-!> Read in real type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength size of array
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_real(this,varname,ilength,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  real(4), intent(out) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='get_var_nc_real'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_FLOAT) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_FLOAT,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
-    if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>read in variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(a,I10)') '        data type : ',this%xtype
+       write(*,'(a,I10)')'        dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(a,I5,I10,2x,a)') '        rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine get_var_nc_double
+
+  !> Read in 1D real type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_real_1d(this,varname,nd1,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    real(4), intent(out) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='get_var_nc_real_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    call this%get_var_nc_real(varname,ilength,field)
+    !
+    if(nd1==this%ends(1)) then
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) (field(i),i=1,min(nd1,10))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+    endif
+    !
+  end subroutine get_var_nc_real_1d
+
+  !> Read in 2D real type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_real_2d(this,varname,nd1,nd2,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    real(4), intent(out) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    real(4),allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='get_var_nc_real_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    call this%get_var_nc_real(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
+       do j=1,nd2
+          istart=(j-1)*nd1+1
+          iend=(j-1)*nd1+nd1
+          field(:,j)=temp(istart:iend)
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_real_2d
+
+  !> Read in 3D real type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] nd3 length of third dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_real_3d(this,varname,nd1,nd2,nd3,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
+    real(4), intent(out) :: field(nd1,nd2,nd3) !  values of the field read in
+    integer :: ilength
+    !
+    real(4),allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='get_var_nc_real_3d'
+    !
+    integer :: i,j,k
+    integer :: length2d
+    integer :: istart,iend
+    !
+    !
+    length2d=nd1*nd2
+    ilength=length2d*nd3
+    allocate(temp(ilength))
+
+    call this%get_var_nc_real(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then
+       do k=1,nd3
+          do j=1,nd2
+             istart=(k-1)*length2d+(j-1)*nd1+1
+             iend  =(k-1)*length2d+(j-1)*nd1+nd1
+             field(:,j,k)=temp(istart:iend)
+          enddo
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          do k=1,nd3
+             write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
+          enddo
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_real_3d
+
+  !> Read in real type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength size of array
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_real(this,varname,ilength,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    real(4), intent(out) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='get_var_nc_real'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_FLOAT) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_FLOAT,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_get_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_get_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>read in variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(8x,a,I10)') 'data type : ',this%xtype
-     write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine get_var_nc_real
-
-!> Read in 1D integer variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_int_1d(this,varname,nd1,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  integer, intent(out) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='get_var_nc_int_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-  call this%get_var_nc_int(varname,ilength,field)
-!
-  if(nd1==this%ends(1)) then
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) (field(i),i=1,min(nd1,10))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-  endif
-!
-end subroutine get_var_nc_int_1d
-
-!> Read in 2D integer type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_int_2d(this,varname,nd1,nd2,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  integer, intent(out) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  integer,allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='get_var_nc_int_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  call this%get_var_nc_int(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
-     do j=1,nd2
-        istart=(j-1)*nd1+1
-        iend=(j-1)*nd1+nd1
-        field(:,j)=temp(istart:iend)
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_int_2d
-
-!> Read in 3D integer type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] nd3 length of third dimension
-!! @param[in] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_int_3d(this,varname,nd1,nd2,nd3,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
-  integer, intent(out) :: field(nd1,nd2,nd3) !  values of the field read in
-  integer :: ilength
-!
-  integer,allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='get_var_nc_int_3d'
-!
-  integer :: i,j,k
-  integer :: length2d
-  integer :: istart,iend
-!
-!
-  length2d=nd1*nd2
-  ilength=length2d*nd3
-  allocate(temp(ilength))
-
-  call this%get_var_nc_int(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then
-     do k=1,nd3
-        do j=1,nd2
-           istart=(k-1)*length2d+(j-1)*nd1+1
-           iend  =(k-1)*length2d+(j-1)*nd1+nd1
-           field(:,j,k)=temp(istart:iend)
-        enddo
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        do k=1,nd3
-           write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
-        enddo
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_int_3d
-
-!> Read in integer type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength size of array
-!! @param[in] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_int(this,varname,ilength,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  integer, intent(out) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='get_var_nc_int'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_INT) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
-    if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>read in variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(8x,a,I10)') 'data type : ',this%xtype
+       write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine get_var_nc_real
+
+  !> Read in 1D integer variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_int_1d(this,varname,nd1,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    integer, intent(out) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='get_var_nc_int_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    call this%get_var_nc_int(varname,ilength,field)
+    !
+    if(nd1==this%ends(1)) then
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) (field(i),i=1,min(nd1,10))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+    endif
+    !
+  end subroutine get_var_nc_int_1d
+
+  !> Read in 2D integer type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_int_2d(this,varname,nd1,nd2,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    integer, intent(out) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    integer,allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='get_var_nc_int_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    call this%get_var_nc_int(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
+       do j=1,nd2
+          istart=(j-1)*nd1+1
+          iend=(j-1)*nd1+nd1
+          field(:,j)=temp(istart:iend)
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_int_2d
+
+  !> Read in 3D integer type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] nd3 length of third dimension
+  !! @param[in] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_int_3d(this,varname,nd1,nd2,nd3,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
+    integer, intent(out) :: field(nd1,nd2,nd3) !  values of the field read in
+    integer :: ilength
+    !
+    integer,allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='get_var_nc_int_3d'
+    !
+    integer :: i,j,k
+    integer :: length2d
+    integer :: istart,iend
+    !
+    !
+    length2d=nd1*nd2
+    ilength=length2d*nd3
+    allocate(temp(ilength))
+
+    call this%get_var_nc_int(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then
+       do k=1,nd3
+          do j=1,nd2
+             istart=(k-1)*length2d+(j-1)*nd1+1
+             iend  =(k-1)*length2d+(j-1)*nd1+nd1
+             field(:,j,k)=temp(istart:iend)
+          enddo
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          do k=1,nd3
+             write(*,*) 'k,max,min:',k,maxval(field(:,:,k)),minval(field(:,:,k))
+          enddo
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_int_3d
+
+  !> Read in integer type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength size of array
+  !! @param[in] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_int(this,varname,ilength,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    integer, intent(out) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='get_var_nc_int'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_INT) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_INT,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_get_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_get_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>read in variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(8x,a,I10)') 'data type : ',this%xtype
-     write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine get_var_nc_int
-
-!> Read in 1D short type variable
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_short_1d(this,varname,nd1,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  integer(2), intent(out) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='get_var_nc_short_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-  call this%get_var_nc_short(varname,ilength,field)
-!
-  if(nd1==this%ends(1)) then
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) (field(i),i=1,min(nd1,10))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-  endif
-!
-end subroutine get_var_nc_short_1d
-
-!> Read in 2D short type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_short_2d(this,varname,nd1,nd2,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  integer(2), intent(out) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  integer(2),allocatable :: temp(:)
-!
-  character*40,parameter :: thissubname='get_var_nc_short_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  call this%get_var_nc_short(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
-     do j=1,nd2
-        istart=(j-1)*nd1+1
-        iend=(j-1)*nd1+nd1
-        field(:,j)=temp(istart:iend)
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_short_2d
-!
-!> Read in short type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength size of array
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_short(this,varname,ilength,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  integer(2), intent(out) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='get_var_nc_short'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_SHORT) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_SHORT,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
-    if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>read in variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(8x,a,I10)') 'data type : ',this%xtype
+       write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine get_var_nc_int
+
+  !> Read in 1D short type variable
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_short_1d(this,varname,nd1,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    integer(2), intent(out) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='get_var_nc_short_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    call this%get_var_nc_short(varname,ilength,field)
+    !
+    if(nd1==this%ends(1)) then
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) (field(i),i=1,min(nd1,10))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+    endif
+    !
+  end subroutine get_var_nc_short_1d
+
+  !> Read in 2D short type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_short_2d(this,varname,nd1,nd2,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    integer(2), intent(out) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    integer(2),allocatable :: temp(:)
+    !
+    character*40,parameter :: thissubname='get_var_nc_short_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    call this%get_var_nc_short(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
+       do j=1,nd2
+          istart=(j-1)*nd1+1
+          iend=(j-1)*nd1+nd1
+          field(:,j)=temp(istart:iend)
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) 'max,min:',maxval(field(:,:)),minval(field(:,:))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_short_2d
+  !
+  !> Read in short type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength size of array
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_short(this,varname,ilength,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    integer(2), intent(out) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='get_var_nc_short'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_SHORT) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_SHORT,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_get_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_get_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>read in variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(8x,a,I10)') 'data type : ',this%xtype
-     write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine get_var_nc_short
-
-!> Read in 1D character type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_char_1d(this,varname,nd1,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1              !  size of array dval
-  character, intent(out) :: field(nd1)     !  values of the field read in
-  integer :: ilength
-!
-  character*40,parameter :: thissubname='get_var_nc_char_1d'
-!
-  integer :: i
-!
-!
-  ilength=nd1
-  call this%get_var_nc_char(varname,ilength,field)
-!
-  if(nd1==this%ends(1)) then
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) (field(i),i=1,min(nd1,10))
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-  endif
-!
-end subroutine get_var_nc_char_1d
-
-!> Read in 2D character type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_char_2d(this,varname,nd1,nd2,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2          !  size of array dval
-  character, intent(out) :: field(nd1,nd2) !  values of the field read in
-  integer :: ilength
-!
-  character,allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='get_var_nc_char_2d'
-!
-  integer :: i,j,k
-  integer :: istart,iend
-!
-!
-  ilength=nd1*nd2
-  allocate(temp(ilength))
-
-  call this%get_var_nc_char(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
-     do j=1,nd2
-        istart=(j-1)*nd1+1
-        iend=(j-1)*nd1+nd1
-        field(:,j)=temp(istart:iend)
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) field(1,1)
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_char_2d
-
-!> Read in 3D character type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] nd1 length of first dimension
-!! @param[in] nd2 length of second dimension
-!! @param[in] nd3 length of third dimension
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_char_3d(this,varname,nd1,nd2,nd3,field)
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
-  character, intent(out) :: field(nd1,nd2,nd3) !  values of the field read in
-  integer :: ilength
-!
-  character,allocatable :: temp(:) 
-!
-  character*40,parameter :: thissubname='get_var_nc_char_3d'
-!
-  integer :: i,j,k
-  integer :: length2d
-  integer :: istart,iend
-!
-!
-  length2d=nd1*nd2
-  ilength=length2d*nd3
-  allocate(temp(ilength))
-
-  call this%get_var_nc_char(varname,ilength,temp)
-
-  if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then
-     do k=1,nd3
-        do j=1,nd2
-           istart=(k-1)*length2d+(j-1)*nd1+1
-           iend  =(k-1)*length2d+(j-1)*nd1+nd1
-           field(:,j,k)=temp(istart:iend)
-        enddo
-     enddo
-!
-     if(this%debug_level>100) then
-        write(*,*) trim(thissubname),' show samples:'
-        write(*,*) field(1,1,1)
-     endif
-  else
-     write(*,*) trim(thissubname),' ERROR: dimension does not match.'
-     write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3)
-  endif
-  deallocate(temp)
-!
-end subroutine get_var_nc_char_3d
-!
-!> Read in character type variable.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] varname name of the variable
-!! @param[in] ilength size of array
-!! @param[out] field output variable
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine get_var_nc_char(this,varname,ilength,field)
-!
-! read in one field 
-!
-  use netcdf
-!
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname  ! name of the field to read
-  integer, intent(in) :: ilength          !  size of array dval
-  character, intent(out) :: field(ilength)   !  values of the field read in
-!
-  integer :: ncid
-! 
-  integer :: status
-  integer :: varid 
-  integer :: ends(4),start(4)
-
-  integer :: length4d,length3d,length2d
-  integer :: nDims,ndim
-  integer :: dimids(4)
-  integer :: xtype
-  character*40 :: dimname
-
-  character*40,parameter :: thissubname='get_var_nc_char'
-!
-  integer :: i,k
-!
-!
-  ncid=this%ncid
-  
-! get variable IDs
-  status = nf90_inq_varid(ncid, trim(varname), VarId)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-
-!  get dimensions
-  ends=1
-  start=1
-  this%ends=1
-
-  this%dimname="                           "
-! get variable type
-  status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  if(xtype==NF90_CHAR) then
-     this%xtype=xtype
-  else
-     write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_CHAR,' but read in ',xtype
-     stop 123
-  endif
-  
-! get dimension size
-  status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  this%ndims=nDims
-!
-  status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
-  if(status /= nf90_NoErr) call this%handle_err(status)
-  do i=1,nDims
-    dimname="       "
-    status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
-    if (status /= nf90_noerr) call this%handle_err(status)
-    ends(i)=ndim
-    this%ends(i)=ends(i)
-    this%dimname(i)=trim(dimname)
-    if(this%ends(i) < 1) then
-       write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>read in variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(8x,a,I10)') 'data type : ',this%xtype
+       write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine get_var_nc_short
+
+  !> Read in 1D character type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_char_1d(this,varname,nd1,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1              !  size of array dval
+    character, intent(out) :: field(nd1)     !  values of the field read in
+    integer :: ilength
+    !
+    character*40,parameter :: thissubname='get_var_nc_char_1d'
+    !
+    integer :: i
+    !
+    !
+    ilength=nd1
+    call this%get_var_nc_char(varname,ilength,field)
+    !
+    if(nd1==this%ends(1)) then
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) (field(i),i=1,min(nd1,10))
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+    endif
+    !
+  end subroutine get_var_nc_char_1d
+
+  !> Read in 2D character type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_char_2d(this,varname,nd1,nd2,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2          !  size of array dval
+    character, intent(out) :: field(nd1,nd2) !  values of the field read in
+    integer :: ilength
+    !
+    character,allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='get_var_nc_char_2d'
+    !
+    integer :: i,j,k
+    integer :: istart,iend
+    !
+    !
+    ilength=nd1*nd2
+    allocate(temp(ilength))
+
+    call this%get_var_nc_char(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2)) then
+       do j=1,nd2
+          istart=(j-1)*nd1+1
+          iend=(j-1)*nd1+nd1
+          field(:,j)=temp(istart:iend)
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) field(1,1)
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_char_2d
+
+  !> Read in 3D character type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] nd1 length of first dimension
+  !! @param[in] nd2 length of second dimension
+  !! @param[in] nd3 length of third dimension
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_char_3d(this,varname,nd1,nd2,nd3,field)
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: nd1,nd2,nd3      !  size of array dval
+    character, intent(out) :: field(nd1,nd2,nd3) !  values of the field read in
+    integer :: ilength
+    !
+    character,allocatable :: temp(:) 
+    !
+    character*40,parameter :: thissubname='get_var_nc_char_3d'
+    !
+    integer :: i,j,k
+    integer :: length2d
+    integer :: istart,iend
+    !
+    !
+    length2d=nd1*nd2
+    ilength=length2d*nd3
+    allocate(temp(ilength))
+
+    call this%get_var_nc_char(varname,ilength,temp)
+
+    if(nd1==this%ends(1) .and. nd2==this%ends(2) .and. nd3==this%ends(3)) then
+       do k=1,nd3
+          do j=1,nd2
+             istart=(k-1)*length2d+(j-1)*nd1+1
+             iend  =(k-1)*length2d+(j-1)*nd1+nd1
+             field(:,j,k)=temp(istart:iend)
+          enddo
+       enddo
+       !
+       if(this%debug_level>100) then
+          write(*,*) trim(thissubname),' show samples:'
+          write(*,*) field(1,1,1)
+       endif
+    else
+       write(*,*) trim(thissubname),' ERROR: dimension does not match.'
+       write(*,*) nd1,this%ends(1),nd2,this%ends(2),nd3,this%ends(3)
+    endif
+    deallocate(temp)
+    !
+  end subroutine get_var_nc_char_3d
+  !
+  !> Read in character type variable.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] varname name of the variable
+  !! @param[in] ilength size of array
+  !! @param[out] field output variable
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine get_var_nc_char(this,varname,ilength,field)
+    !
+    ! read in one field 
+    !
+    use netcdf
+    !
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname  ! name of the field to read
+    integer, intent(in) :: ilength          !  size of array dval
+    character, intent(out) :: field(ilength)   !  values of the field read in
+    !
+    integer :: ncid
+    ! 
+    integer :: status
+    integer :: varid 
+    integer :: ends(4),start(4)
+
+    integer :: length4d,length3d,length2d
+    integer :: nDims,ndim
+    integer :: dimids(4)
+    integer :: xtype
+    character*40 :: dimname
+
+    character*40,parameter :: thissubname='get_var_nc_char'
+    !
+    integer :: i,k
+    !
+    !
+    ncid=this%ncid
+
+    ! get variable IDs
+    status = nf90_inq_varid(ncid, trim(varname), VarId)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+
+    !  get dimensions
+    ends=1
+    start=1
+    this%ends=1
+
+    this%dimname="                           "
+    ! get variable type
+    status = nf90_inquire_variable(ncid, VarId, xtype = xtype)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    if(xtype==NF90_CHAR) then
+       this%xtype=xtype
+    else
+       write(*,*) trim(thissubname),' ERROR: wrong data type, expect ',NF90_CHAR,' but read in ',xtype
+       stop 123
+    endif
+
+    ! get dimension size
+    status = nf90_inquire_variable(ncid, VarId, ndims = nDims)
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    this%ndims=nDims
+    !
+    status = nf90_inquire_variable(ncid, VarId, dimids = dimids(1:nDims))
+    if(status /= nf90_NoErr) call this%handle_err(status)
+    do i=1,nDims
+       dimname="       "
+       status = nf90_inquire_dimension(ncid, dimids(i), dimname, len = ndim)
+       if (status /= nf90_noerr) call this%handle_err(status)
+       ends(i)=ndim
+       this%ends(i)=ends(i)
+       this%dimname(i)=trim(dimname)
+       if(this%ends(i) < 1) then
+          write(*,*) trim(thissubname),' Error, ends dimension should larger than 0 :', ends(i)
+          stop 1234
+       endif
+    enddo
+    length2d=ends(1)*ends(2)
+    length3d=length2d*ends(3)
+    length4d=length3d*ends(4)
+    if(ilength .ne. length4d) then
+       write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
+       stop 123
+    endif
+    !
+    if(nDims <=4 ) then
+       status = nf90_get_var(ncid, VarId, field, &
+            start = start(1:4) , &
+            count = ends(1:4))
+       if(status /= nf90_NoErr) call this%handle_err(status)
+    else
+       write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
        stop 1234
     endif
-  enddo
-  length2d=ends(1)*ends(2)
-  length3d=length2d*ends(3)
-  length4d=length3d*ends(4)
-  if(ilength .ne. length4d) then
-     write(*,*) trim(thissubname),'ERROR: ',ilength,' should equal to ',length4d
-     stop 123
-  endif
-!
-  if(nDims <=4 ) then
-     status = nf90_get_var(ncid, VarId, field, &
-                         start = start(1:4) , &
-                         count = ends(1:4))
-     if(status /= nf90_NoErr) call this%handle_err(status)
-  else
-     write(*,*) trim(thissubname),'Error: too many dimensions:',nDims
-     stop 1234
-  endif
-!
-  if(this%debug_level>0) then
-     write(*,'(a,a)') '>>>read in variable: ',trim(varname)
-  endif
-  if(this%debug_level>10) then
-     write(*,'(8x,a,I10)') 'data type : ',this%xtype
-     write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
-     do i=1,this%nDims
-       write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
-     enddo
-  endif
-!
-end subroutine get_var_nc_char
-
-!> Handle netCDF errors.
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] status return code from neCDF
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine handle_err(this,status)
-  use netcdf
-  implicit none
-  class(ncio) :: this
-!
-  integer, intent ( in) :: status
-  if(status /= nf90_noerr) then
-      print *, trim(nf90_strerror(status))
-      stop "Stopped"
-  end if
-end subroutine handle_err
-
-!> Convert theta T (Kelvin) to T (deg C).
-!!
-!! @param[in] this instance of an ncio class
-!! @param[in] nx number of grid points in x-dir
-!! @param[in] ny number of grid points in y-dir
-!! @param[in] ps Pressure (Pa)
-!! @param[inout] t2 Pot. Temperature (Kelvin)
-!! @author Ming Hu org: GSD/AMB @date 2017-11-01
-subroutine convert_theta2t_2dgrid(this,nx,ny,ps,t2)
-         implicit none
-         class(ncio) :: this
-
-         integer :: nx,ny
-         real, intent(in   ) :: ps(nx,ny)
-         real, intent(inout) :: t2(nx,ny)
-
-         integer :: i,j
-         real(8) :: rd,cp,rd_over_cp
-
-
-         rd     = 2.8705e+2_8
-         cp     = 1.0046e+3_8  !  specific heat of air @pressure (J/kg/K)
-         rd_over_cp = rd/cp
-
-         do j=1,ny
-            do i=1,nx
-               t2(i,j)=t2(i,j)*(ps(i,j)/1000.0)**rd_over_cp - 273.15
-            enddo
-         enddo
-
-end subroutine convert_theta2t_2dgrid
-
-!> Add a new variable to sfc_data.nc with dimensions (Time, yaxis_1,
-!! xaxis_1).
-!!
-!! @param[in] varname Name of variable to be created in netcdf file
-!! @param[in] dname1 1st dimension name
-!! @param[in] dname2 2nd dimension name
-!! @param[in] dname3 3rd dimension name
-!! @param[in] lname long name output for netcdf variable
-!! @param[in] units units to use in netcdf variable
-!!
-!! @author David.M.Wright org: UM/GLERL @date 2020-09-01
-subroutine add_new_var_3d(this,varname,dname1,dname2,dname3,lname,units)
-  implicit none
-!
-  class(ncio) :: this
-  character(len=*),intent(in) :: varname,dname1,dname2,dname3 &
-                                 ,lname,units
-  integer :: status, ncid, dim1id, dim2id, dim3id, varid
-
-  status = nf90_redef(this%ncid) !Enter Define Mode
-  if (status /= nf90_noerr) call this%handle_err(status)
-
-  status = nf90_inq_dimid(this%ncid, dname1, dim1id)
-  if (status /= nf90_noerr) call this%handle_err(status)
-  status = nf90_inq_dimid(this%ncid, dname2, dim2id)
-  if (status /= nf90_noerr) call this%handle_err(status)
-  status = nf90_inq_dimid(this%ncid, dname3, dim3id)
-  if (status /= nf90_noerr) call this%handle_err(status)
-
-  status = nf90_def_var(this%ncid, varname, nf90_double, &
-                         (/ dim1id, dim2id, dim3id /), varid)
-  if (status /= nf90_noerr) call this%handle_err(status)
-
-  status = nf90_put_att(this%ncid, varid, 'long_name', lname)
-  if (status /= nf90_noerr) call this%handle_err(status)
-  status = nf90_put_att(this%ncid, varid, 'units', units)
-  if (status /= nf90_noerr) call this%handle_err(status)
-
-  status = nf90_enddef(this%ncid) !Exit Define Mode and
-                                  ! return to Data Mode
-  if (status /= nf90_noerr) call this%handle_err(status)
-
-end subroutine add_new_var_3d
+    !
+    if(this%debug_level>0) then
+       write(*,'(a,a)') '>>>read in variable: ',trim(varname)
+    endif
+    if(this%debug_level>10) then
+       write(*,'(8x,a,I10)') 'data type : ',this%xtype
+       write(*,'(8x,a,I10)') 'dimension size: ',this%nDims
+       do i=1,this%nDims
+          write(*,'(8x,a,I5,I10,2x,a)') 'rank, ends, name=',i,this%ends(i),trim(this%dimname(i))
+       enddo
+    endif
+    !
+  end subroutine get_var_nc_char
+
+  !> Handle netCDF errors.
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] status return code from neCDF
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine handle_err(this,status)
+    use netcdf
+    implicit none
+    class(ncio) :: this
+    !
+    integer, intent ( in) :: status
+    if(status /= nf90_noerr) then
+       print *, trim(nf90_strerror(status))
+       stop "Stopped"
+    end if
+  end subroutine handle_err
+
+  !> Convert theta T (Kelvin) to T (deg C).
+  !!
+  !! @param[in] this instance of an ncio class
+  !! @param[in] nx number of grid points in x-dir
+  !! @param[in] ny number of grid points in y-dir
+  !! @param[in] ps Pressure (Pa)
+  !! @param[inout] t2 Pot. Temperature (Kelvin)
+  !! @author Ming Hu org: GSD/AMB @date 2017-11-01
+  subroutine convert_theta2t_2dgrid(this,nx,ny,ps,t2)
+    implicit none
+    class(ncio) :: this
+
+    integer :: nx,ny
+    real, intent(in   ) :: ps(nx,ny)
+    real, intent(inout) :: t2(nx,ny)
+
+    integer :: i,j
+    real(8) :: rd,cp,rd_over_cp
+
+
+    rd     = 2.8705e+2_8
+    cp     = 1.0046e+3_8  !  specific heat of air @pressure (J/kg/K)
+    rd_over_cp = rd/cp
+
+    do j=1,ny
+       do i=1,nx
+          t2(i,j)=t2(i,j)*(ps(i,j)/1000.0)**rd_over_cp - 273.15
+       enddo
+    enddo
+
+  end subroutine convert_theta2t_2dgrid
+
+  !> Add a new variable to sfc_data.nc with dimensions (Time, yaxis_1,
+  !! xaxis_1).
+  !!
+  !! @param this instance of an ncio class
+  !! @param[in] varname Name of variable to be created in netcdf file
+  !! @param[in] dname1 1st dimension name
+  !! @param[in] dname2 2nd dimension name
+  !! @param[in] dname3 3rd dimension name
+  !! @param[in] lname long name output for netcdf variable
+  !! @param[in] units units to use in netcdf variable
+  !!
+  !! @author David.M.Wright org: UM/GLERL @date 2020-09-01
+  subroutine add_new_var_3d(this,varname,dname1,dname2,dname3,lname,units)
+    implicit none
+    !
+    class(ncio) :: this
+    character(len=*),intent(in) :: varname,dname1,dname2,dname3 &
+         ,lname,units
+    integer :: status, ncid, dim1id, dim2id, dim3id, varid
+
+    status = nf90_redef(this%ncid) !Enter Define Mode
+    if (status /= nf90_noerr) call this%handle_err(status)
+
+    status = nf90_inq_dimid(this%ncid, dname1, dim1id)
+    if (status /= nf90_noerr) call this%handle_err(status)
+    status = nf90_inq_dimid(this%ncid, dname2, dim2id)
+    if (status /= nf90_noerr) call this%handle_err(status)
+    status = nf90_inq_dimid(this%ncid, dname3, dim3id)
+    if (status /= nf90_noerr) call this%handle_err(status)
+
+    status = nf90_def_var(this%ncid, varname, nf90_double, &
+         (/ dim1id, dim2id, dim3id /), varid)
+    if (status /= nf90_noerr) call this%handle_err(status)
+
+    status = nf90_put_att(this%ncid, varid, 'long_name', lname)
+    if (status /= nf90_noerr) call this%handle_err(status)
+    status = nf90_put_att(this%ncid, varid, 'units', units)
+    if (status /= nf90_noerr) call this%handle_err(status)
+
+    status = nf90_enddef(this%ncid) !Exit Define Mode and
+    ! return to Data Mode
+    if (status /= nf90_noerr) call this%handle_err(status)
+
+  end subroutine add_new_var_3d
 
 end module module_ncio
diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90
index c3f7682ed..be9b22e12 100644
--- a/sorc/fvcom_tools.fd/module_nwp.f90
+++ b/sorc/fvcom_tools.fd/module_nwp.f90
@@ -26,31 +26,38 @@ module module_nwp
 
    private
    type :: nwp_type
-      character(len=6) :: datatype
-      integer :: numvar, xlat, xlon, xtime
-      integer :: i_mask, i_sst, i_ice, i_sfcT, i_iceT
-      character(len=20), allocatable :: varnames(:)
-      character(len=20), allocatable :: latname
-      character(len=20), allocatable :: lonname
-      character(len=20), allocatable :: dimnameEW
-      character(len=20), allocatable :: dimnameNS
-      character(len=20), allocatable :: dimnameTIME
-
-      real(r_kind), allocatable :: nwp_mask(:,:,:)
-      real(r_kind), allocatable :: nwp_sst(:,:,:)
-      real(r_kind), allocatable :: nwp_ice(:,:,:)
-      real(r_kind), allocatable :: nwp_sfcT(:,:,:)
-      real(r_kind), allocatable :: nwp_iceT(:,:,:) 
+      character(len=6) :: datatype !< Data type.
+      integer :: numvar !< Number of variabls.
+      integer :: xlat !< Number of latitudes.
+      integer :: xlon !< Number of longitudes.
+      integer :: xtime !< Number of times.
+      integer :: i_mask !< Is var visible (always 1).
+      integer :: i_sst !< Index of sst var.
+      integer :: i_ice !< Index of ice var.
+      integer :: i_sfcT !< Index of sst temp var.
+      integer :: i_iceT !< Index of ice temp var.
+      character(len=20), allocatable :: varnames(:) !< Variable names.
+      character(len=20), allocatable :: latname !< Latitude name.
+      character(len=20), allocatable :: lonname !< Longitude name.
+      character(len=20), allocatable :: dimnameEW !< East/West dimension name.
+      character(len=20), allocatable :: dimnameNS !< North/South dimension name.
+      character(len=20), allocatable :: dimnameTIME !< Time dimension name.
+
+      real(r_kind), allocatable :: nwp_mask(:,:,:) !< ???
+      real(r_kind), allocatable :: nwp_sst(:,:,:) !< ???
+      real(r_kind), allocatable :: nwp_ice(:,:,:) !< ???
+      real(r_kind), allocatable :: nwp_sfcT(:,:,:) !< ???
+      real(r_kind), allocatable :: nwp_iceT(:,:,:)  !< ???
    end type nwp_type
 
    type, extends(nwp_type) :: fcst_nwp
-      type(nwpbase), pointer :: head => NULL()
-      type(nwpbase), pointer :: tail => NULL()
+      type(nwpbase), pointer :: head => NULL() !< ???
+      type(nwpbase), pointer :: tail => NULL() !< ???
       contains
-         procedure :: initial => initial_nwp
-         procedure :: list_initial => list_initial_nwp
-         procedure :: read_n => read_nwp
-         procedure :: finish => finish_nwp
+         procedure :: initial => initial_nwp !< Defines vars and names. @return
+         procedure :: list_initial => list_initial_nwp !< List the setup. @return
+         procedure :: read_n => read_nwp !< Initialize arrays, get data. @return
+         procedure :: finish => finish_nwp !< Finish and deallocate. @return
    end type fcst_nwp
 
    type(ncio) :: ncdata !< Wrapper object for netCDF data file.
@@ -66,8 +73,6 @@ module module_nwp
      !! @param[in] itype either ' FVCOM' or 'FV3LAM'.
      !! @author David Wright, University of Michigan and GLERL
       subroutine initial_nwp(this,itype)
-
-
          class(fcst_nwp) :: this
 
          character(len=6), intent(in) :: itype
diff --git a/sorc/fvcom_tools.fd/module_nwp_base.f90 b/sorc/fvcom_tools.fd/module_nwp_base.f90
index 87275c3da..382826458 100644
--- a/sorc/fvcom_tools.fd/module_nwp_base.f90
+++ b/sorc/fvcom_tools.fd/module_nwp_base.f90
@@ -25,26 +25,25 @@ module module_nwp_base
 !  Define a nwp observation type.
 
    type nwplocation
-      real(r_single) :: lon  ! stroke longitude
-      real(r_single) :: lat  ! stroke latitiude
+      real(r_single) :: lon  !< stroke longitude
+      real(r_single) :: lat  !< stroke latitiude
    end type nwplocation
 
 !  Define a nwp observation type to contain actual data.
 
    type, extends(nwplocation) :: nwpbase
 !  HOW DOES THIS POINTER THING WORK?
-      type(nwpbase), pointer :: next => NULL()
-      real(r_single) :: time                 ! observation time
-      integer :: numvar                      ! number of variables in this obs type
-!      real(r_single), allocatable :: obs(:)  ! observation value (# numvar)
-      real(r_kind), allocatable :: obs(:)
-      logical :: ifquality                   ! do these obs include quality info?
-!                                              GLM has flash_quality_flag
-      integer, allocatable :: quality(:)     ! if so, quality flags
+      type(nwpbase), pointer :: next => NULL() !< Pointer.
+      real(r_single) :: time                 !< observation time.
+      integer :: numvar                      !< number of variables in this obs type.
+!      real(r_single), allocatable :: obs(:)  !< observation value (# numvar).
+      real(r_kind), allocatable :: obs(:) !< Observations.
+      logical :: ifquality                   !< do these obs include quality info? GLM has flash_quality_flag.
+      integer, allocatable :: quality(:)     !< if so, quality flags.
       contains
-         procedure :: list => list_obsbase
-         procedure :: alloc => alloc_obsbase
-         procedure :: destroy => destroy_obsbase
+         procedure :: list => list_obsbase !< List contents of obs. @return
+         procedure :: alloc => alloc_obsbase !< Allocate memory for observations. @return
+         procedure :: destroy => destroy_obsbase  !< Release memory. @return
    end type nwpbase
 
    contains
@@ -52,10 +51,9 @@ module module_nwp_base
      !> This subroutine lists the contents of a base nwp observation.
      !!
      !! @param this the base nwp obervation
-     !! @author David Wright, University of Michigan and GLERL @date 17 Aug 2020     
+     !! @author David Wright, University of Michigan and GLERL
+     !! @date 17 Aug 2020
       subroutine list_obsbase(this)
-
-
          class(nwpbase) :: this
 
          integer :: i, numvar
@@ -84,7 +82,7 @@ end subroutine list_obsbase
       !!
       !! @param this the base nwp obervation
       !! @param[in] numvar number of variables in this ob type
-      !! @param[in] itquality does this observation include quality
+      !! @param[in] ifquality does this observation include quality
       !! information?
       !! @author David Wright, University of Michigan and GLERL @date 17 Aug 2020
       subroutine alloc_obsbase(this,numvar,ifquality)

From 83960192094a15e51f55c76ab71ee1a856bbbdd1 Mon Sep 17 00:00:00 2001
From: Edward Hartnett <Edward.Hartnett@noaa.gov>
Date: Tue, 2 Mar 2021 08:04:03 -0700
Subject: [PATCH 08/10] fixed typo

---
 sorc/fvcom_tools.fd/module_ncio.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sorc/fvcom_tools.fd/module_ncio.f90 b/sorc/fvcom_tools.fd/module_ncio.f90
index 77b0fc685..318d9f172 100644
--- a/sorc/fvcom_tools.fd/module_ncio.f90
+++ b/sorc/fvcom_tools.fd/module_ncio.f90
@@ -66,7 +66,7 @@ module module_ncio
           replace_var_nc_int_1d,replace_var_nc_int_2d,        &
           replace_var_nc_int_3d,                              &
           replace_var_nc_char_1d,replace_var_nc_char_2d,      &
-          replace_var_nc_char_3d file. !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
+          replace_var_nc_char_3d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
      procedure :: replace_var_nc_int !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
      procedure :: replace_var_nc_int_1d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return
      procedure :: replace_var_nc_int_2d !< Replace 1d, 2d, 3d, or 4d field from the nc file. @return

From ea40794e28bb244fa9b3031bb5caec197c47d0ca Mon Sep 17 00:00:00 2001
From: David Wright <dmwright526@gmail.com>
Date: Tue, 2 Mar 2021 13:43:41 -0500
Subject: [PATCH 09/10] Update module_nwp.f90 with comments

Identified arrays nwp_mask, nwp_sst, nwp_ice, nwp_sfcT, and nwp_iceT with more detailed comments.
---
 sorc/fvcom_tools.fd/module_nwp.f90 | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90
index be9b22e12..9addf8870 100644
--- a/sorc/fvcom_tools.fd/module_nwp.f90
+++ b/sorc/fvcom_tools.fd/module_nwp.f90
@@ -43,11 +43,11 @@ module module_nwp
       character(len=20), allocatable :: dimnameNS !< North/South dimension name.
       character(len=20), allocatable :: dimnameTIME !< Time dimension name.
 
-      real(r_kind), allocatable :: nwp_mask(:,:,:) !< ???
-      real(r_kind), allocatable :: nwp_sst(:,:,:) !< ???
-      real(r_kind), allocatable :: nwp_ice(:,:,:) !< ???
-      real(r_kind), allocatable :: nwp_sfcT(:,:,:) !< ???
-      real(r_kind), allocatable :: nwp_iceT(:,:,:)  !< ???
+      real(r_kind), allocatable :: nwp_mask(:,:,:) !< Land/water mask 3D array
+      real(r_kind), allocatable :: nwp_sst(:,:,:) !< SST 3D array
+      real(r_kind), allocatable :: nwp_ice(:,:,:) !< Over water ice concentration 3D array
+      real(r_kind), allocatable :: nwp_sfcT(:,:,:) !< Skin temperature 3D array
+      real(r_kind), allocatable :: nwp_iceT(:,:,:)  !< Ice skin temperature 3D array
    end type nwp_type
 
    type, extends(nwp_type) :: fcst_nwp

From 97428239989afdf5ba4e7193dcd49ca2d823d964 Mon Sep 17 00:00:00 2001
From: Edward Hartnett <Edward.Hartnett@noaa.gov>
Date: Tue, 2 Mar 2021 12:43:51 -0700
Subject: [PATCH 10/10] more doxygen changes

---
 sorc/fvcom_tools.fd/module_nwp.f90 | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/sorc/fvcom_tools.fd/module_nwp.f90 b/sorc/fvcom_tools.fd/module_nwp.f90
index 9addf8870..c3c925333 100644
--- a/sorc/fvcom_tools.fd/module_nwp.f90
+++ b/sorc/fvcom_tools.fd/module_nwp.f90
@@ -51,8 +51,12 @@ module module_nwp
    end type nwp_type
 
    type, extends(nwp_type) :: fcst_nwp
-      type(nwpbase), pointer :: head => NULL() !< ???
-      type(nwpbase), pointer :: tail => NULL() !< ???
+      ! The pointers are carryover from when I inherited the code from
+      ! GSL's work with HRRR for a similar use. I am not sure with
+      ! object based coding in Fortran if it needs to have parts
+      ! initialized to gain access to the procedures within it. - D. Wright.
+      type(nwpbase), pointer :: head => NULL() !< Pointer to head of list.
+      type(nwpbase), pointer :: tail => NULL() !< Pointer to tail of list.
       contains
          procedure :: initial => initial_nwp !< Defines vars and names. @return
          procedure :: list_initial => list_initial_nwp !< List the setup. @return