!-----------------------------------------------------------------------
! Copyright 2011 NOAA Geophysical Fluid Dynamics Lab, Princeton, NJ
! This program is distributed under the terms of the GNU General Public
! License. See the file COPYING contained in this directory
!
!  This program aggregates vertical data by pressure levels
!
!-----------------------------------------------------------------------
 
program run_pressure_interp

use netcdf

use    plev_constants_mod, only: GRAV, RDGAS, RVGAS

use pinterp_utilities_mod, only: copy_axis, copy_variable,  &
                                 set_verbose_level,         &
                                 define_new_variable,       &
                                 error_handler
use pressure_interp_mod


implicit none

!-----------------------------------------------------------------------

integer, parameter :: MAX_PLEVS = 200
integer, parameter :: MAX_FIELD_NAMES = 300
integer, parameter :: MAX_AXIS_NAMES  =  50
integer, parameter :: header_buffer_val = 16384
integer            ::  blksz = 65536  ! blksz must be writable for nf__create

 character(len=256)           ::  in_file_name, out_file_name
!character(len=NF90_MAX_NAME) ::  in_file_name, out_file_name
  real, dimension(MAX_PLEVS) :: pout
character(len=NF90_MAX_NAME) ::    ps_name = 'ps',     &
                                zsurf_name = 'zsurf',  &
                                 temp_name = 'temp',   &
                                sphum_name = 'sphum',  &
                                   bk_name = 'bk',     &
                                   pk_name = 'pk',     &
                                  res_name = 'res',    &
                                 hght_name = 'hght',   &
                                  slp_name = 'slp'

 character(len=NF90_MAX_NAME) :: pres_name_out = 'level'
!character(len=NF90_MAX_NAME) :: time_name_out = 'time'   ! not used

 character(len=NF90_MAX_NAME), dimension(MAX_FIELD_NAMES) :: field_names

type field_type
   character(len=NF90_MAX_NAME) :: name
   integer                      :: varid_in, varid_out, ndim, xtype
   integer                      :: dimids(NF90_MAX_VAR_DIMS), &
                                   dimlen(NF90_MAX_VAR_DIMS)
   integer                      :: y_axis, z_axis
   real                         :: pack(2), missval_in, unpkd_missval_in, missval_out, unpkd_missval_out
   logical                      :: static, skip, do_avg, do_pack, do_miss_in, do_miss_out
end type

type axis_type
   character(len=NF90_MAX_NAME) :: name_in, name_out
   integer                      :: varid_in, varid_out
   integer                      :: dimlen, start(1)
end type

type(field_type) :: Fields(MAX_FIELD_NAMES)
type(axis_type)  :: Axes(MAX_AXIS_NAMES)

logical :: mask_extrap = .false.
logical :: unpack = .false.  ! not implemented
integer :: verbose     = 1

logical :: do_all_3d_fields = .false.
logical :: do_all_fields = .false.
logical :: use_default_missing_value = .false.

logical :: allow_zero_topog = .false.
logical :: allow_zero_sphum = .false.

integer :: time_beg=1, time_end=0, time_inc=1

!#######################################################################
!         -------- NAMELIST -------
!
!  in_file_name = input file name/path
!               must contain atmospheric dynamics fields (required)
!
!  out_file_name = output file name/path
!
!  field_names = additional fields (up to 40) for output
!                (dynamics files are searched for these names)
!
!         pout = output pressure levels (max 200), in pascals
!
!--- specify these names for the following output fields ---
!         sphum_name = specific humidity
!--------------------------------------------------------
!
!  mask_extrap = set points below the bottom level to the missing value
!
!  use_default_missing_value = use the NetCDF default for missing/fill value
!                              for all fields
!
!   verbose   = integer verbose flag
!
!-----------------------------------------------------------------------
!
!  *** NOTES ***
!
!    1)  currently only ONE input atmospheric file and ONE output file
!        are allowed
!
!    2)  the namelist file must be passed to the program by
!        redirection: e.g., a.out < namelist
!
!#######################################################################

namelist /names/ temp_name, sphum_name, hght_name, slp_name, &
                 ps_name, zsurf_name, bk_name, pk_name 

namelist /input/ in_file_name, out_file_name,  field_names,  pout, &
                 do_all_fields,  do_all_3d_fields,                 &
                 mask_extrap,                          &
                 time_beg, time_end, time_inc,                     &
                 allow_zero_topog, allow_zero_sphum,               &
                 use_default_missing_value, verbose

!----- axis names -----

!character(len=NF90_MAX_NAME) :: xmass_in = 'lon',  &
!                                ymass_in = 'lat',  &
!character(len=NF90_MAX_NAME) :: pfull_in = 'pfull', &
!                                phalf_in = 'phalf', &
!                                time_in  = 'time'


character(len=NF90_MAX_NAME) :: xmass_in, ymass_in
character(len=NF90_MAX_NAME) :: pfull_in, time_in

character(len=NF90_MAX_NAME) :: xmass_out, ymass_out
character(len=NF90_MAX_NAME) :: edge_name
character(len=NF90_MAX_NAME) :: attname
character(len=NF90_MAX_NAME) :: file_title

integer :: dimids_mass(4)
!-----------------------------------------------------------------------

integer :: ncid_in, ncid_out
type (pres_interp_type) :: Interp(2)


character(len=NF90_MAX_NAME) :: version
character(len=NF90_MAX_NAME) :: time_avg_info
character(len=NF90_MAX_NAME) :: cell_methods
character(len=NF90_MAX_NAME) :: name_time_bounds, name_time_average

logical    :: do_time_average, do_time_bounds
character(len=NF90_MAX_NAME) :: avg_name = 'average'
character(len=8) :: avg_subnames(3) = (/ '_T1', '_T2', '_DT' /)

logical    :: temp_present, sphum_present
logical    :: exists

logical :: do_temp  = .false.
logical :: do_sphum = .false.
logical :: do_zsurf = .false.
logical :: do_hght  = .false.
logical :: do_slp   = .false.

logical :: need_zsurf, need_temp, need_sphum

integer :: varid_res,      varid_pout
integer :: varid_bk,       varid_pk,     &
           varid_temp_in,  varid_sphum_in, &
           varid_zsurf_in, varid_ps_in

integer :: varid_hght_out, varid_slp_out

integer :: varid_time_in,        varid_time_out
integer :: varid_time_bounds_in, varid_time_bounds_out
integer :: varid_tavg_in(3),     varid_tavg_out(3)

integer :: index_ps, index_zsurf, index_temp, index_sphum

integer :: istat, ivar, itime, otime, i, n, kv, attnum, numatts, recdim,  &
           nlon, nlat, nlev, ndim, dimids(NF90_MAX_VAR_DIMS), dimid, naxes, dimid_tbnds
integer :: nlev_out, start(4), kount(4), ntimes, num_field_names
integer :: itime_beg, itime_end, itime_inc
integer :: isg, ieg, jsg, jeg, is, ie, js, je
real    :: mval_hght, mval_slp
real(8) :: time, time_bnds(2)
logical :: done, do_phalf
integer :: idim, nlen
integer :: unit, ierr, io, in_format, cmode
real, allocatable :: axisdata(:)

real, allocatable, dimension(:)     :: bk, pk
real, allocatable, dimension(:,:)   :: res
real, allocatable, dimension(:,:,:) :: zsurf, ps

real, parameter :: tlapse = 6.5e-3
real, parameter :: tref   = 288.
real, parameter :: ginv = 1./GRAV
real, parameter :: gorg = GRAV / (RDGAS*tlapse)
real, parameter :: mrgog = -1./gorg
real, parameter :: d608 = (RVGAS-RDGAS)/RDGAS
real, parameter :: pref = 101325.
real, parameter :: rgog = RDGAS*tlapse/GRAV

!#######################################################################

!--------- initialize namelist variable and then read namelist ---------

      do i=1,len(out_file_name);  out_file_name(i:i) = ' '; enddo
      do i=1,len( in_file_name);   in_file_name(i:i) = ' '; enddo
      pout = 0.0

      do ivar = 1, MAX_FIELD_NAMES
         do i=1,NF90_MAX_NAME; field_names(ivar)(i:i) = ' '; enddo
      enddo

      inquire (file='plev.input.nml',exist=exists)
      if (exists) then
         open (10, file='plev.input.nml')
         read (10, nml=input, end=1, err=1)
      1  close (10)
         !write (*,input)
      endif

         do i = 1, MAX_PLEVS
            if (pout(i) < 1.e-6) exit
            nlev_out = i
         enddo
 
!     ---- count the number of field names ----
         num_field_names = 0
      do ivar = 1, MAX_FIELD_NAMES
         if (field_names(ivar)(1:1) == ' ') exit
         num_field_names = ivar
      enddo

      call set_verbose_level (verbose)

  ! create version string (may replace with CVS $Id: run_pressure_interp.F90,v 20.0 2013/12/14 00:29:42 fms Exp $)
    version = 'pressure level interpolator, version 3.0'
    if (precision(ps) == precision(time)) then
        version = trim(version)//', precision=double'
    else
        version = trim(version)//', precision=float'
    endif

!-----------------------------------------------------------------------
!----------------------- loop through files ----------------------------
!-----------------------------------------------------------------------
!   --- if any file name has no length then terminate the program ---
      
     if (  in_file_name(1:1) == ' ' .or.   &
          out_file_name(1:1) == ' ' )  call error_handler  &
           ('must specifiy atmospheric and output file names')

!-----------------------------------------------------------------------
!-----------------  get input file information -------------------------

    !--- opening input files ---
     istat = NF90_OPEN (trim(in_file_name), NF90_NOWRITE, ncid_in, blksz)
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)

    !--- get the time axis identifier ---
     istat = NF90_INQUIRE ( ncid_in, nattributes=numatts, unlimiteddimid=recdim,formatnum=in_format )
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)

!    --------- set flags for special computed fields --------

     if (verbose > 1) print *, 'num_field_names=',num_field_names
     do i = 1, num_field_names
        if (verbose > 1) print *, 'field_name=',trim(field_names(i))
        if (field_names(i) ==  hght_name) then
            do_hght  = .true.; cycle; endif
        if (field_names(i) ==   slp_name) then
            do_slp   = .true.; cycle; endif
     enddo

!    ------- check to see if required fields are in atmos file -------

     call check_fields ( ncid_in )

!-----------------------------------------------------------------------
!--- get the size of the vertical axis ---
     istat = NF90_INQUIRE_VARIABLE (ncid_in, varid_bk, ndims=ndim, dimids=dimids)
     if (istat /= NF90_NOERR) call error_handler ('getting dimids for '//trim(bk_name), ncode=istat)
     istat = NF90_INQUIRE_DIMENSION (ncid_in, dimids(1), len=nlev)
     if (istat /= NF90_NOERR) call error_handler ('getting dim length for '//trim(bk_name), ncode=istat)
     nlev = nlev-1

    !--- get the names of the input horizontal mass axes ---
     istat = NF90_INQUIRE_VARIABLE (ncid_in, varid_ps_in, ndims=ndim, dimids=dimids)
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)

     if (ndim /= 3) call error_handler ('variable '//trim(ps_name)//&
                                    &' must have three dimensions: x,y,t')

    !--- size and name of input horizontal mass axes ---
     istat = NF90_INQUIRE_DIMENSION (ncid_in, dimids(1), name=xmass_in, len=nlon)
     if (istat /= NF90_NOERR) call error_handler ('getting name and length of first mass dimension', ncode=istat)
     istat = NF90_INQUIRE_DIMENSION (ncid_in, dimids(2), name=ymass_in, len=nlat)
     if (istat /= NF90_NOERR) call error_handler ('getting name and length of second mass dimension', ncode=istat)
     istat = NF90_INQUIRE_DIMENSION (ncid_in, dimids(3), name=time_in)
     if (istat /= NF90_NOERR) call error_handler ('getting name of last mass dimension (time)', ncode=istat)

!-----------------------------------------------------------------------
!--- set output fields using all input fields ---

     if (do_all_fields .or. do_all_3d_fields) then
         call set_all_field_names ( ncid_in, field_names, num_field_names )
     endif

!--- error: if there were no requested fields ---
     if (num_field_names == 0) call error_handler &
                               ('no field names have been selected')

!--- set flags for special non-computed fields ---

     if (verbose > 1) print *, 'num_field_names=',num_field_names
     do i = 1, num_field_names
        if (verbose > 1) print *, 'field_name=',trim(field_names(i))
        if (field_names(i) == sphum_name) then
            if (sphum_present) do_sphum = .true.; cycle; endif
        if (field_names(i) ==  temp_name) then
            if (temp_present)  do_temp  = .true.; cycle; endif
        if (field_names(i) == zsurf_name) then
            do_zsurf = .true.; cycle; endif
     enddo

!-----------------------------------------------------------------------

     isg = 1;   ieg = nlon; jsg = 1;   jeg = nlat
     is  = isg; ie  = ieg;  js  = jsg; je  = jeg

!    ---- allocate space ----

     allocate (bk(nlev+1), pk(nlev+1))
     allocate (res  (is:ie,js:je),    &
               zsurf(is:ie,js:je,1),  &
               ps   (is:ie,js:je,1)   )

!    ---- read vert coord variables and set up vert coord ----

     istat = NF90_GET_VAR (ncid_in, varid_bk, bk)
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)
     if (varid_pk > 0) then
        istat = NF90_GET_VAR (ncid_in, varid_pk, pk)
        if (istat /= NF90_NOERR) call error_handler (ncode=istat)
     else
        pk = 0.0
     endif

!-----------------------------------------------------------------------
!----------------- output file header info -----------------------------
     if (in_format==NF90_FORMAT_NETCDF4) then
        cmode = NF90_NETCDF4
     elseif (in_format==NF90_FORMAT_NETCDF4_CLASSIC) then
        cmode=IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL)
     elseif (in_format==NF90_FORMAT_64BIT) then
        cmode=IOR(NF90_CLOBBER,NF90_64BIT_OFFSET)
     elseif (in_format==NF90_FORMAT_CLASSIC) then
        cmode=IOR(NF90_CLOBBER,NF90_CLASSIC_MODEL)
     else
        call error_handler ('Unknown netCDF format')
     endif

     istat = NF90_CREATE (trim(out_file_name), cmode, ncid_out, chunksize=blksz)
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)

   ! copy global attributes
     do attnum = 1, numatts
        ! get name
        istat = NF90_INQ_ATTNAME (ncid_in, NF90_GLOBAL, attnum, attname)
        if (istat /= NF90_NOERR) call error_handler (ncode=istat)
        ! change global filename attribute
        if (trim(attname) == 'filename') then
            istat = NF90_PUT_ATT (ncid_out, NF90_GLOBAL, 'filename', trim(out_file_name))
            if (istat /= NF90_NOERR) call error_handler (ncode=istat)
        ! append to history attribute
        else if (trim(attname) == 'history') then
           !istat = NF90_INQUIRE_ATTRIBUTE (ncid_in, NF90_GLOBAL, trim(attname), len=nc)
           !istat = NF90_GET_ATT (ncid_in, NF90_GLOBAL, 'history', history)
        else    
            ! copy all others
            istat = NF90_COPY_ATT (ncid_in, NF90_GLOBAL, attname, ncid_out, NF90_GLOBAL)
            if (istat /= NF90_NOERR) call error_handler ('copying global attributes', &
                                                          ncode=istat)
        endif   
     enddo
     ! add comment attribute
     istat = NF90_PUT_ATT (ncid_out, NF90_GLOBAL, 'comment', trim(version))
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)

!-----------------------------------------------------------------------
!------------------ generate axes for output file ----------------------

    xmass_out = xmass_in;  ymass_out = ymass_in

    !---- copy mass axes ----
     naxes = 1
     call copy_axis_information (ncid_in, xmass_in, ncid_out, (/isg,ieg,is,ie/), Axes(naxes), &
                                 dimid_out=dimids_mass(1), name_out=xmass_out, name_edges=edge_name)
         ! copy edge axis if available
         if (edge_name(1:1) .ne. ' ') then
             naxes = naxes + 1
             call copy_axis_information (ncid_in, edge_name, ncid_out, (/isg,ieg+1,is,ie+1/), Axes(naxes))
         endif

     naxes = naxes + 1
     call copy_axis_information (ncid_in, ymass_in, ncid_out, (/jsg,jeg,js,je/), Axes(naxes), &
                                 dimid_out=dimids_mass(2), name_out=ymass_out, name_edges=edge_name)
         ! copy edge axis if available
         if (edge_name(1:1) .ne. ' ') then
             naxes = naxes + 1
             call copy_axis_information (ncid_in, edge_name, ncid_out, (/jsg,jeg+1,js,je+1/), Axes(naxes))
         endif

!    -------- vertical axis (convert to mb) -------

   ! define new output pressure dimension and variable
     istat = NF90_DEF_DIM (ncid_out, trim(pres_name_out), nlev_out, dimids_mass(3))
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)
     istat = NF90_DEF_VAR (ncid_out, trim(pres_name_out), NF90_REAL4, (/dimids_mass(3)/), varid_pout)
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)
   ! attributes for this axis
     istat = NF90_PUT_ATT (ncid_out, varid_pout, 'units', 'hPa')
     if (istat /= NF90_NOERR) call error_handler ('putting attributes for pout dimension', ncode=istat)
     istat = NF90_PUT_ATT (ncid_out, varid_pout, 'long_name', 'pressure')
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)
     istat = NF90_PUT_ATT (ncid_out, varid_pout, 'cartesian_axis', 'Z')
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)
     istat = NF90_PUT_ATT (ncid_out, varid_pout, 'axis', 'Z') ! for CF compliance
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)
     istat = NF90_PUT_ATT (ncid_out, varid_pout, 'positive', 'down')
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)

   ! copy the time axis
   ! get variable identifiers and length
     call copy_axis (ncid_in, time_in, ncid_out, dimids_mass(4), &
                     varid_time_in, varid_time_out, ntimes)

   ! set up time axis indexing
     itime_beg = time_beg
     itime_end = time_end
     itime_inc = time_inc
     if (itime_end == 0 .or. itime_end > ntimes) itime_end = ntimes
     if (itime_beg <= 0) itime_beg = 1
     if (itime_inc <= 0) itime_inc = 1
    !if (itime_beg < itime_end) itime_end = itime_beg  ! process at least one time level


   ! check for time averaging "bounds"
   ! make sure axis is copied (time axis already copied)
     do i=1,NF90_MAX_NAME;  name_time_bounds(i:i) = ' '; enddo
     do i=1,NF90_MAX_NAME;  cell_methods(i:i) = ' '    ; enddo
     istat = NF90_GET_ATT (ncid_in, varid_time_in, 'bounds', name_time_bounds)
     if (istat /= NF90_NOERR) then
         istat = NF90_GET_ATT (ncid_in, varid_time_in, 'climatology', name_time_bounds)
     endif
     if (istat == NF90_NOERR) then
         if (verbose > 1) print *, 'name_time_bounds=',trim(name_time_bounds)
         istat = NF90_INQ_VARID (ncid_in, trim(name_time_bounds), varid_time_bounds_in)
         if (istat /= NF90_NOERR) call error_handler ('getting varid for time bounds', ncode=istat)
         istat = NF90_INQUIRE_VARIABLE (ncid_in, varid_time_bounds_in, ndims=ndim, dimids=dimids)
         if (istat /= NF90_NOERR) call error_handler ('getting dimensions for time bounds', ncode=istat)
         if (ndim == 1 .or. ndim == 2) then
             istat = NF90_INQUIRE_DIMENSION (ncid_in, dimids(1), name=edge_name)
             if (istat == NF90_NOERR) then
                 naxes = naxes + 1
                 call copy_axis_information (ncid_in, edge_name, ncid_out, (/1,1,1,1/), &
                                             Axes(naxes), dimid_out=dimid_tbnds)
             endif
         else
             call error_handler ('number of dimensions not consistent for '//trim(name_time_bounds))
         endif
         ! get cell methods for new variables
         istat = NF90_GET_ATT (ncid_in, varid_ps_in, 'cell_methods', cell_methods)
         do_time_bounds = .true.
     else
         do_time_bounds = .false.
     endif

!-----------------------------------------------------------------------
!-------------- define new variables for output file -----------------

!  ---- geopotential height ----

if (do_hght) then
   mval_hght = -1000.
   if (use_default_missing_value) mval_hght = NF90_FILL_REAL4

   call define_new_variable ( ncid_out, 'hght', NF90_REAL4, dimids_mass, varid_hght_out, &
                              units='m', long_name='height', missing_value=mval_hght,    &
                              time_avg_info=time_avg_info, cell_methods=cell_methods)
   
endif

!  ---- sea level pressure ----

if (do_slp) then
   mval_slp = -1.
   if (use_default_missing_value) mval_slp = NF90_FILL_REAL4

   call define_new_variable ( ncid_out, 'slp', NF90_REAL4,                      &
                              (/dimids_mass(1),dimids_mass(2),dimids_mass(4)/), &
                              varid_slp_out, units='hPa', long_name='sea level pressure', &
                              missing_value=mval_slp, time_avg_info=time_avg_info,   &
                              cell_methods=cell_methods)
endif

!-----------------------------------------------------------------------
!---- copy variable info for other fields ------

    Fields(:)%varid_in = 0
    Fields(:)%y_axis = 0
    Fields(:)%z_axis = 0
    Fields(:)%static = .false.   ! only allowing 2d static fields
    do_phalf = .false.

    do i = 1, num_field_names

       call copy_variable_information ( field_names(i), ncid_in, ncid_out, Fields(i) )

       ! save index for special fields
       if (trim(field_names(i)) == trim(   ps_name)) index_ps    = i
       if (trim(field_names(i)) == trim(zsurf_name)) index_zsurf = i
       if (trim(field_names(i)) == trim( temp_name)) index_temp  = i
       if (trim(field_names(i)) == trim(sphum_name)) index_sphum = i

       ! check valid fields for half levels then go to end of loop
       if (.not.Fields(i)%skip) then
           if (Fields(i)%z_axis == 2) do_phalf = .true.
           cycle
       endif

       if (trim(field_names(i)) == trim( hght_name)) cycle
       if (trim(field_names(i)) == trim(  slp_name)) cycle
       print *, 'WARNING: unable to process field '//trim(field_names(i))
    enddo

!-----------------------------------------------------------------------
!---- input variable only are needed to process some field ----
!---- check to see if these variables are in the field_names list ----

      if (.not.field_name_present(ps_name,field_names,num_field_names)) then
          num_field_names = num_field_names + 1
          index_ps = num_field_names
          call get_input_variable_info ( ps_name, ncid_in, Fields(num_field_names) )
      endif

      if (need_zsurf) then
          if (.not.field_name_present(zsurf_name,field_names,num_field_names)) then
              num_field_names = num_field_names + 1
              index_zsurf = num_field_names
              call get_input_variable_info ( zsurf_name, ncid_in, Fields(num_field_names) )
          endif
      endif

      if (need_temp) then
          if (.not.field_name_present(temp_name,field_names,num_field_names)) then
              num_field_names = num_field_names + 1
              index_temp = num_field_names
              call get_input_variable_info ( temp_name, ncid_in, Fields(num_field_names) )
          endif
      endif

      if (need_sphum) then
          if (.not.field_name_present(sphum_name,field_names,num_field_names)) then
              num_field_names = num_field_names + 1
              index_sphum = num_field_names
              call get_input_variable_info ( sphum_name, ncid_in, Fields(num_field_names) )
          endif
      endif

!-----------------------------------------------------------------------
! ---- copy variable info for time averaging -----

      if (do_time_bounds) then
          call copy_variable ( ncid_in, trim(name_time_bounds), ncid_out,    &
                               (/dimid_tbnds,dimids_mass(4)/), varid_time_bounds_out )
      endif

      if (do_time_average) then
          do i = 1, 3
             name_time_average = trim(avg_name)//trim(avg_subnames(i))
             istat = NF90_INQ_VARID (ncid_in, trim(name_time_average), varid_tavg_in(i))
             if (istat /= NF90_NOERR) call error_handler (ncode=istat)
             call copy_variable ( ncid_in, trim(name_time_average), ncid_out,    &
                                  (/dimids_mass(4)/), varid_tavg_out(i) )
          enddo
       endif
          
!-----------------------------------------------------------------------
! ---- write header info for output file ----

      istat = NF90_ENDDEF (ncid_out,header_buffer_val,4,0,4)
      if (istat /= NF90_NOERR) call error_handler (ncode=istat)

!-----------------------------------------------------------------------
!------------ read and write static axis data -------------

       do idim = 1, naxes
          if (Axes(idim)%varid_in < 0) cycle
          nlen = Axes(idim)%dimlen
          allocate (axisdata(nlen))
          istat = NF90_GET_VAR (ncid_in, Axes(idim)%varid_in, axisdata, start=Axes(idim)%start)
          if (istat /= NF90_NOERR) call error_handler ('getting axis data for var='//trim(Axes(idim)%name_in), ncode=istat)
          istat = NF90_PUT_VAR (ncid_out, Axes(idim)%varid_out, axisdata)
          if (istat /= NF90_NOERR) call error_handler ('putting axis data for var='//trim(Axes(idim)%name_out), ncode=istat)
          deallocate (axisdata)
       enddo
        ! output pressure levels in millibars
          istat = NF90_PUT_VAR (ncid_out, varid_pout, pout(1:nlev_out)*0.01)
          if (istat /= NF90_NOERR) call error_handler ('putting axis data for var='//trim(pres_name_out), ncode=istat)

!-----------------------------------------------------------------------
!------------ first time: read constant fields -------------------------

       start(1:3) = (/ is, js, 1 /)
       istat = NF90_INQ_VARID (ncid_in, trim(res_name), varid_res)
       if ( istat == NF90_NOERR ) then
          istat = NF90_GET_VAR (ncid_in, varid_res, res, start=start(1:2))
          if (istat /= NF90_NOERR) call error_handler (ncode=istat)
       else
          res = 1.0
       endif
       if (need_zsurf .or. do_zsurf) then
          call read_variable (ncid_in, Fields(index_zsurf), zsurf, start=start(1:2))
       else
          zsurf = 0.0
       endif
       if (do_zsurf) then
          call write_variable (ncid_out, Fields(index_zsurf), zsurf)
       endif
!      --- eta coord fix ----
       if (count(res > 1.0001) > 0) zsurf = 0.0

!-----------------------------------------------------------------------
!      ---- read/write static fields ----

    do i = 1, num_field_names
       if ( Fields(i)%skip )         cycle
       if ( .not.Fields(i)%static )  cycle
       if (trim(Fields(i)%name) == trim(zsurf_name)) cycle
       call read_write_static (ncid_in, ncid_out, is, ie, js, je, Fields(i))
    enddo

!-----------------------------------------------------------------------
!-------------------- loop through time axis ---------------------------
!-----------------------------------------------------------------------

      otime = 0
      do itime = itime_beg, itime_end, itime_inc
          otime = otime+1
          if ( verbose > 0 ) then
              if (itime == otime) then
                 print '(a,i8)', 'itime =', itime
              else
                 print '(a,2i8)', 'itime, otime =', itime, otime
              endif
          endif

!---------------- read/write time level --------------

          istat =  NF90_GET_VAR (ncid_in, varid_time_in, time, start=(/itime/))
          if (istat /= NF90_NOERR) call error_handler ('getting time coord value', &
                                                       ncode=istat)

          istat =  NF90_PUT_VAR (ncid_out, varid_time_out, time, start=(/otime/))
          if (istat /= NF90_NOERR) call error_handler ('putting time coord value', &
                                                       ncode=istat)

!-----------------------------------------------------------------------
!-------------- interpolation setup ------------------------------------

     !---- surface pressure ----

          start(1:3) = (/ is, js, itime /)
          call read_variable (ncid_in, Fields(index_ps), ps, start=start(1:3))
          ps(is:ie,js:je,1) = ps(is:ie,js:je,1) * res

!-----------------------------------------------------------------------
!----------------- mass field interpolation ---------------------

!     ---- interpolation initialization for mass fields ------

     call interp_field_setup ( ncid_in, ncid_out, bk, pk,   &
                               is, js, itime, otime,        &
                               ps(is:ie,js:je,1),           &
                               zsurf(:,:,1), pout(1:nlev_out) )

!     ---- interpolation of mass fields ----

     do i = 1, num_field_names
        if ( Fields(i)%skip   )    cycle
        if ( Fields(i)%static )    cycle
        if ( Fields(i)%y_axis /= nlat ) cycle

        kv = max(1,Fields(i)%z_axis)
        if ( verbose > 1 ) then
           print *, 'Interp field = ',trim(Fields(i)%name)
        endif

!     ---- interpolation of mass fields ----
        call interp_field ( ncid_in, ncid_out, Interp(kv),  &
                            is, js, itime, otime, Fields(i) )

     enddo

!     ---- free up allocated space ----
                    call pres_interp_free (Interp(1))
      if (do_phalf) call pres_interp_free (Interp(2))

!-----------------------------------------------------------------------
!   ---- copy time average data ----

    if (do_time_bounds) then
        istat = NF90_GET_VAR (ncid_in,  varid_time_bounds_in,  time_bnds(:), start=(/1,itime/))
        if (istat /= NF90_NOERR) call error_handler ('getting time_bounds', ncode=istat)
        istat = NF90_PUT_VAR (ncid_out, varid_time_bounds_out, time_bnds(:), start=(/1,otime/))
        if (istat /= NF90_NOERR) call error_handler ('putting time_bounds', ncode=istat)
    endif

    if (do_time_average) then
        do i = 1, 3
           istat = NF90_GET_VAR (ncid_in,  varid_tavg_in(i),  time_bnds(1), start=(/itime/))
           if (istat /= NF90_NOERR) call error_handler ('getting tavg_info', ncode=istat)
           istat = NF90_PUT_VAR (ncid_out, varid_tavg_out(i), time_bnds(1), start=(/otime/))
           if (istat /= NF90_NOERR) call error_handler ('putting tavg_info', ncode=istat)
        enddo
    endif

!-----------------------------------------------------------------------

   enddo  ! end of time loop

!-----------------------------------------------------------------------

       istat = NF90_CLOSE (ncid_in)
       if (istat /= NF90_NOERR) call error_handler (ncode=istat)
       istat = NF90_CLOSE (ncid_out)
       if (istat /= NF90_NOERR) call error_handler (ncode=istat)

!-----------------------------------------------------------------------

contains

!#######################################################################

   subroutine interp_field_setup ( ncid_in, ncid_out, bk, pk, &
                                   is, js, itime, otime, ps,  &
                                   zsurf, pout )

   integer, intent(in) :: ncid_in, ncid_out
   real,    intent(in) :: bk(:), pk(:)
   integer, intent(in) :: is, js, itime, otime
   real   , intent(in) :: ps(:,:), zsurf(:,:), pout(:)

 real, dimension(size(ps,1),size(ps,2),size(bk(:))-1) :: &
            temp_in, sphum_in, mask_in, pfull, hght_in
 real, dimension(size(ps,1),size(ps,2),size(bk(:)))   :: phalf
 real, dimension(size(ps,1),size(ps,2),size(pout)) ::  &
                         temp_out, hght_out, sphum_out
 real   , dimension(size(ps,1),size(ps,2)) :: zbot, tbot, pbot, slp
 integer, dimension(size(ps,1),size(ps,2)) :: kbot
   character(len=NF90_MAX_NAME) :: name
   integer :: start(4), kount(4), ndim, ivar, nlev, kb, kr, i, j, k
   real    :: sig


      do k = 1, size(bk(:))
         phalf(:,:,k) = pk(k) + bk(k)*ps(:,:)
      enddo
      call compute_pres_full (pk, bk, phalf, pfull)

  !---- read in temperature data ----
   if (do_temp .or. need_temp) then
       start = (/ is, js, 1, itime /)
       call read_variable (ncid_in, Fields(index_temp), temp_in, start)
      !---- mask ----
       call set_up_input_mask (ncid_in, Fields(index_temp), temp_in, mask_in, kbot)
   endif

  !---- read in specific humidity data ----
   if (do_sphum .or. need_sphum) then
       if (sphum_present) then
           start = (/ is, js, 1, itime /)
           call read_variable (ncid_in, Fields(index_sphum), sphum_in, start)
          !---- mask ----
           if (.not.do_temp) &
           call set_up_input_mask (ncid_in, Fields(index_sphum), sphum_in, mask_in, kbot)
       else
           sphum_in = 0.0
       endif
   endif

  !---- interpolation setup for standard full pressure levels -----
   if (do_hght) then
       if ( verbose > 1 ) print *, 'Computing field = hght'
      !---- geopotential height (in meters) -----
       call compute_height ( zsurf, temp_in, sphum_in, pfull, phalf, hght_in, mask_in)
      !---- interpolate temp, sphum, hght ----
       if ( verbose > 1 ) then
           print *, 'Interp field = ',trim(Fields(index_temp)%name)
           print *, 'Interp field = ',trim(Fields(index_sphum)%name)
       endif
       Interp(1) = pres_interp_init ( pfull, pout, kbot,   &
                                      tin=temp_in,  tout=temp_out,  &
                                      zin=hght_in,  zout=hght_out,  &
                                      qin=sphum_in, qout=sphum_out, &
                                      use_extrap=.not.mask_extrap   )
   else if ((do_temp.and.do_sphum).and.(.not.do_hght)) then
       if ( verbose > 1 ) then
           print *, 'Interp field = ',trim(Fields(index_temp)%name)
           print *, 'Interp field = ',trim(Fields(index_sphum)%name)
       endif
      !---- interpolate temp, sphum ----
       Interp(1) = pres_interp_init ( pfull, pout, kbot,   &
                                      tin=temp_in,  tout=temp_out,  &
                                      qin=sphum_in, qout=sphum_out, &
                                      use_extrap=.not.mask_extrap   )
   else if ((do_temp).and.(.not.do_hght)) then
       if ( verbose > 1 ) print *, 'Interp field = ',trim(Fields(index_temp)%name)
      !---- interpolate temp ----
       Interp(1) = pres_interp_init ( pfull, pout, kbot,   &
                                      tin=temp_in,  tout=temp_out,  &
                                      use_extrap=.not.mask_extrap   )
   else if ((do_sphum).and.(.not.do_hght)) then
       if ( verbose > 1 ) print *, 'Interp field = ',trim(Fields(index_sphum)%name)
      !---- interpolate sphum ----
       Interp(1) = pres_interp_init ( pfull, pout, kbot,   &
                                      qin=sphum_in, qout=sphum_out, &
                                      use_extrap=.not.mask_extrap   )
   else
      !---- initialize only (need kbot) ----
       start = (/ is, js, 1, itime /)
       call read_and_set_up_mask (ncid_in, start, mask_in, kbot)
       Interp(1) = pres_interp_init ( pfull, pout, kbot,   &
                                      use_extrap=.not.mask_extrap   )
   endif
   
  !---- interpolation setup for standard full pressure levels -----
   if (do_phalf) then
       Interp(2) = pres_interp_init ( phalf, pout, kbot )
   endif

! ---- write out data ----

     !start = (/ is, js,  1, otime /)
      start = (/  1,  1,  1, otime /)

!     ---- temperature ----

   if (do_temp) then
      if ( verbose > 1 ) print *, 'Writing field = ',trim(Fields(index_temp)%name)
      if (Fields(index_temp)%do_miss_out) then
          call mask_data (Fields(index_temp)%unpkd_missval_out, Interp(1)%mask, temp_out)
      endif
      call write_variable (ncid_out, Fields(index_temp), temp_out, start)
   endif

!        ---- specific humidity ----

   if (do_sphum) then
      if ( verbose > 1 ) print *, 'Writing field = ',trim(Fields(index_sphum)%name)
      if (Fields(index_sphum)%do_miss_out) then
          call mask_data (Fields(index_sphum)%unpkd_missval_out, Interp(1)%mask, sphum_out)
      endif
      call write_variable (ncid_out, Fields(index_sphum), sphum_out, start)
   endif

!        ---- height ----

   if (do_hght) then
      if ( verbose > 1 ) print *, 'Writing field = hght'
      call mask_data (mval_hght, Interp(1)%mask, hght_out)
      istat = NF90_PUT_VAR (ncid_out, varid_hght_out, hght_out, start)
      if (istat /= NF90_NOERR) call error_handler ('putting hght data', ncode=istat)
   endif

!        ---- sea level pressure ----

   if (do_slp) then
     if ( verbose > 1 ) print *, 'Computing field = slp'
     do j = 1, size(ps,2)
     do i = 1, size(ps,1)
        kb = kbot(i,j)
        pbot(i,j) = phalf(i,j,kb+1)

        if ( abs(zsurf(i,j)) > 0.0001 ) then

!            ---- get ref level for temp (as in spectral model) ----
             do k = 1, kb
                sig = pfull(i,j,k)/pbot(i,j)
                if ( sig > 0.8 ) then
                     kr = k
                     exit
                endif
             enddo

            !zbot(i,j) = zhalf(i,j,kb+1)
!!!          tbot(i,j) = tref - tlapse*zbot(i,j)           ! original
             tbot(i,j) = temp_in(i,j,kr) * sig ** mrgog    ! spectral
             slp(i,j) = 0.01 * pbot(i,j) *  &
                       ( 1.0 + tlapse * zsurf(i,j) / tbot(i,j) ) ** gorg
                      !( 1.0 + tlapse * zbot(i,j) / tbot(i,j) ) ** gorg
        else
             slp(i,j) = pbot(i,j) * 0.01
        endif
     enddo
     enddo

     if ( verbose > 1 ) print *, 'Writing field = slp'
     start(1:3) = (/ is, js, otime /)
     istat = NF90_PUT_VAR (ncid_out, varid_slp_out, slp, start)
     if (istat /= NF90_NOERR) call error_handler ('putting slp data', ncode=istat)
   endif


   end subroutine interp_field_setup

!#######################################################################

   subroutine interp_field ( ncid_in, ncid_out, Interp, is, js, itime, otime, Field )

   integer               , intent(in) :: ncid_in, ncid_out
   type(pres_interp_type), intent(in) :: Interp
   integer               , intent(in) :: is, js, itime, otime
   type(field_type)      , intent(in) :: Field

   real, dimension(Interp%nx,Interp%ny,Interp%kx) :: data_in
   real, dimension(Interp%nx,Interp%ny,Interp%nz) :: data_out

   integer :: start_in(4), start_out(4), istat

!     ---- skip fields have been processed in the initialization routine
       if (trim(Field%name) == trim(temp_name))  return
       if (trim(Field%name) == trim(sphum_name)) return

!     ---- read in data ----

       start_in (1:3) = (/ is, js, 1 /)
      !start_out(1:3) = (/ is, js, 1 /)
       start_out(1:3) = (/  1,  1, 1 /)
       start_in (Field%ndim) = itime
       start_out(Field%ndim) = otime

!     ---- check horiz axis lengths ----

      !if ( Field%dimlen(1) /= Interp%nx .or. Field%dimlen(2) /= Interp%ny ) then
      !    print *, 'x,y=',Field%dimlen(1),Interp%nx,Field%dimlen(2),Interp%ny
      !    call error_handler (   &
      !      'horizontal coordinate resolution inconsistent, var = '//trim(Field%name))
      !endif


!  ------ 2-d or 3-d data allowed -----

 select case ( Field%ndim )

   case (4)
!     ---- check vert axis length ----
     if ( Field%dimlen(3) /= Interp%kx ) print *, trim(Field%name)
     if ( Field%dimlen(3) /= Interp%kx ) call error_handler (    &
                'vertical coordinate resolution inconsistent')
     call read_variable (ncid_in, Field, data_in, start_in(1:4))
!    ---- interpolation ----
     call pres_interp (Interp, data_in, data_out)
     if (Field%do_miss_out) &
     call mask_data (Field%unpkd_missval_out, Interp%mask, data_out)
!    ---- write out data ----
     call write_variable (ncid_out, Field, data_out, start_out(1:4))

   case (3)
     call read_variable (ncid_in, Field, data_in(:,:,1:1), start_in(1:3))
     call write_variable (ncid_out, Field, data_in(:,:,1:1), start_out(1:3))

   case default
     print *, 'ERROR: only 2-d or 3-d fields allowed'
     stop

 end select



 end subroutine interp_field

!#######################################################################

  subroutine mask_data (mval, mask, data)

   real            , intent(in)    :: mval
   logical         , intent(in)    :: mask(:,:,:)
   real            , intent(inout) :: data(:,:,:)

     if (.not.mask_extrap) return
     where (mask) data = mval

  end subroutine mask_data

!#######################################################################

 subroutine check_fields ( ncid )

   integer, intent(in) :: ncid

   integer :: id(6), ierr, iwarn, i
   integer :: avg(6), istat(6)
   character(len=NF90_MAX_NAME), dimension(6) :: names

!-----------------------------------------------------------------------

     names = (/ bk_name, pk_name, zsurf_name, ps_name, temp_name, sphum_name /)

     need_zsurf = .false.
     need_temp  = .false.
     need_sphum = .false.

     id = 0

   ! get variable index for all important fields
     do i = 1, 6
       istat(i) = NF90_INQ_VARID (ncid, trim(names(i)), id(i))
     enddo

!------ perform checks ------

       ierr  = 0
       iwarn = 0

    ! minimum: bk, ps (, pk)
      if (istat(1) /= NF90_NOERR) then
          print *, 'ERROR: required field does not exist: ',trim(names(1))
          ierr = 1
      endif
      if (istat(2) /= NF90_NOERR) then
          print *, 'WARNING: required field does not exist: ',trim(names(2))
          iwarn = 1
      endif
      if (istat(4) /= NF90_NOERR) then
          print *, 'ERROR: required field does not exist: ',trim(names(4))
          ierr = 1
      else
          ! check if time average to be done
          avg(4) = NF90_INQUIRE_ATTRIBUTE (ncid, id(4), 'time_avg_info')
      endif
    ! time averaging for all variable must be consistent
    ! initialize all variables to surf pres value
      avg(5:6) = avg(4)

    ! surf height & temperature
      temp_present = .false.;  if (istat(5) == NF90_NOERR ) temp_present = .true.
      if (do_hght .or. do_slp) then
          if (istat(3) /= NF90_NOERR) then
              if (allow_zero_topog) then
                  print *, 'NOTE: field does not exist: ',trim(names(3)),' = 0 assumed'
              else
                  print *, 'ERROR: required field does not exist: ',trim(names(3))
                  ierr = 1
                  do_hght = .false.
                  do_slp  = .false.
              endif
          else
              need_zsurf = .true.
          endif
          if (temp_present) then
              avg(5) = NF90_INQUIRE_ATTRIBUTE (ncid, id(5), 'time_avg_info')
              need_temp = .true.
          else
              print *, 'ERROR: required field does not exist: ',trim(names(5))
              ierr = 1
              do_hght = .false.
              do_slp  = .false.
          endif
      endif

    ! specific humidity
      sphum_present = .false.;  if (istat(6) == NF90_NOERR ) sphum_present = .true.
      if (do_hght) then
          if (sphum_present) then
              avg(6) = NF90_INQUIRE_ATTRIBUTE (ncid, id(6), 'time_avg_info')
              need_sphum = .true.
          else
              if (allow_zero_sphum) then
                  print *, 'NOTE: field does not exist: ',trim(names(6)),' = 0 assumed'
              else
                  print *, 'ERROR: required field does not exist: ',trim(names(6))
                  ierr = 1
                  do_hght = .false.
              endif
          endif
      
      endif


   ! set variable identifiers in the input file
      varid_bk       = id(1)
      varid_pk       = id(2)
      varid_zsurf_in = id(3)
      varid_ps_in    = id(4)
      varid_temp_in  = id(5)
      varid_sphum_in = id(6)

      if ( ierr  /= 0 ) call error_handler ('run_pressure_interp')
      if ( iwarn /= 0 ) print *, 'WARNING: run_pressure_interp'
             
!-----------------------------------------------------------------------
   ! check time averaging flag for consistency
     do i = 5, 6
        if (avg(4) /= avg(i)) call error_handler ('inconsistent time averages')
     enddo

    ! set averaging flag for program ----
      do_time_average = avg(4) == NF90_NOERR

    ! set metadata time average string
      do i=1,NF90_MAX_NAME;  time_avg_info(i:i) = ' '; enddo
      if (do_time_average) then
         time_avg_info = trim(avg_name)//trim(avg_subnames(1))//','// &
                         trim(avg_name)//trim(avg_subnames(2))//','// &
                         trim(avg_name)//trim(avg_subnames(3))
      endif
      
 end subroutine check_fields

!#######################################################################

 subroutine set_all_field_names ( ncid, field_names, num_field_names )

   integer,          intent(in)    :: ncid
   character(len=*), intent(inout) :: field_names(:)
   integer,          intent(inout) :: num_field_names

   character(len=NF90_MAX_NAME) :: name
   integer :: ivar, idim, rdim, nvar, ndim, i, istat
   integer :: dimids(NF90_MAX_VAR_DIMS), axlen(NF90_MAX_VAR_DIMS)
   logical :: static, threed

   if ( num_field_names >= MAX_FIELD_NAMES ) return

   istat = NF90_INQUIRE (ncid, nvariables=nvar, unlimiteddimid=rdim)
   if (istat /= NF90_NOERR) call error_handler (ncode=istat)

     do ivar = 1, nvar

        ! get number and names of dimensions
        istat = NF90_INQUIRE_VARIABLE (ncid, ivar, name=name, ndims=ndim, dimids=dimids)
        if (istat /= NF90_NOERR) call error_handler (ncode=istat)
        ! get dimension lengths
        do idim = 1, ndim
           istat = NF90_INQUIRE_DIMENSION (ncid, dimids(idim), len=axlen(idim))
           if (istat /= NF90_NOERR) call error_handler (ncode=istat)
        enddo

!     ----  has this name been specified already? ----
       if (field_name_present(name,field_names,num_field_names)) cycle

!     ---- skip surface pressure and topography fields ----
        if ( trim(name) == 'res') cycle    ! for eta coordinate

!     ---- must be 2-d, 3-d, or 4-d field ----
        static = dimids(ndim) .ne. rdim
        threed = .false.
!!!del  if (static) cycle   ! will not allow static for now (see below)

!     ---- skip 3-d fields if flag set ----

        select case (ndim)
          case (4)
            if (static) cycle   ! not possible
            threed = .true. 
          case (3)
           !if (.not.static .and. do_all_3d_fields) cycle
            if (static) threed = .true.
            if (.not.static .and. .not.do_all_fields) cycle
           !if (static) cycle           ! no 3d static fields allowed
          case (2)
            if (.not.static) cycle
            if (.not.do_all_fields) cycle
          case default
            cycle
        end select

!     ---- axes must have correct size ----
        if (axlen(1) /= nlon .or.  axlen(2) /= nlat)   cycle
        if (threed) then
        if (axlen(3) /= nlev .and. axlen(3) /= nlev+1) cycle
        endif

!     ---- output this field ----
        num_field_names = num_field_names + 1
        field_names(num_field_names) = name
        if ( num_field_names == MAX_FIELD_NAMES ) return

     enddo

 end subroutine set_all_field_names

!#######################################################################

 function field_name_present ( name, field_names, num_field_names ) 
 character(len=*), intent(in) :: name, field_names(:)
 integer,          intent(in) :: num_field_names
 logical :: field_name_present
 integer :: i

     field_name_present = .false.

    !----  has this name been specified already? ----
     do i = 1, num_field_names
        if (trim(name) == trim(field_names(i))) then
           field_name_present = .true.
           return
        endif
     enddo

 end function field_name_present

!#######################################################################

 function skip_field_name ( field_name )

   character(len=NF90_MAX_NAME), intent(in) :: field_name
   logical :: skip_field_name
   integer :: i

!  ------ skip fields that are written separately ------

       !skip_field_name = .false.

       !if (trim(field_name) == trim(sphum_name)) then
       !    skip_field_name = .true.; return; endif
       !if (trim(field_name) ==  trim(temp_name)) then
       !    skip_field_name = .true.; return; endif
       !if (trim(field_name) == trim(zsurf_name)) then
       !    skip_field_name = .true.; return; endif

        if (trim(field_name) ==  trim(hght_name)) then
            skip_field_name = .true.; return; endif
        if (trim(field_name) ==   trim(slp_name)) then
            skip_field_name = .true.; return; endif

 end function skip_field_name

!#######################################################################

 subroutine get_input_variable_info ( name, ncid_in, Field )

   character(len=NF90_MAX_NAME), intent(in)    :: name
   integer          ,       intent(in)    :: ncid_in
   type(field_type) ,       intent(inout) :: Field
   logical :: answer
   integer :: rdim, zcoord, idim
   integer :: i, istat, ipack, ioffs, imiss, ifill
   real    :: mval, fval

 ! reset skip flag
 ! skip all fields unless properties are acceptable
   Field%skip = .true.
   Field%varid_in = 0

 ! skip certain field names
   if (trim(name) == trim(hght_name))  return
   if (trim(name) == trim(slp_name))   return

 ! get time axis id
   istat = NF90_INQUIRE (ncid_in, unlimiteddimid=rdim)
   if (istat /= NF90_NOERR) call error_handler (ncode=istat)
 ! get variable id 
   istat = NF90_INQ_VARID (ncid_in, trim(name), Field%varid_in)

 ! proceed if this variable exists
   if (istat /= NF90_NOERR) return

       if (verbose > 2) print *, 'getting input variable info, name=',trim(name)

       ! get number and names of dimensions
       istat = NF90_INQUIRE_VARIABLE (ncid_in, Field%varid_in, name=Field%name, &
                        xtype=Field%xtype, ndims=Field%ndim, dimids=Field%dimids)
       if (istat /= NF90_NOERR) call error_handler (ncode=istat)
       ! get dimension lengths
       do idim = 1, Field%ndim
          istat = NF90_INQUIRE_DIMENSION (ncid_in, Field%dimids(idim), len=Field%dimlen(idim))
          if (istat /= NF90_NOERR) call error_handler (ncode=istat)
       enddo
       Field%static = Field%dimids(Field%ndim) .ne. rdim
       Field%z_axis = 0
       Field%y_axis = 0

     ! check time averaging
       istat = NF90_INQUIRE_ATTRIBUTE (ncid_in, Field%varid_in, 'time_avg_info')
       Field%do_avg = istat == NF90_NOERR

     ! get packing attributes of input data
       ipack = NF90_GET_ATT (ncid_in, Field%varid_in, 'scale_factor', Field%pack(1))
       ioffs = NF90_GET_ATT (ncid_in, Field%varid_in, 'add_offset',   Field%pack(2))
       if (ipack == NF90_NOERR .and. ioffs == NF90_NOERR) then
           Field%do_pack = .true.
       else
           Field%do_pack = .false.
           Field%pack = (/1.,0./)
       endif

     ! get missing value attribute or fill value attribute
     ! NOTE: conversion of missing value data type to default real type
       imiss = NF90_GET_ATT (ncid_in, Field%varid_in, 'missing_value', mval)
       ifill = NF90_GET_ATT (ncid_in, Field%varid_in, '_FillValue',    fval)
       if (imiss == NF90_NOERR) then
          Field%do_miss_in = .true.
          Field%missval_in = mval
       else if (ifill == NF90_NOERR) then
          Field%do_miss_in = .true.
          Field%missval_in = fval
       else
          Field%do_miss_in = .false.
          Field%missval_in = 0.
       endif

     ! unpack input missing value (for comparison with unpacked data)
       Field%unpkd_missval_in = Field%missval_in
       if (Field%do_miss_in .and. Field%do_pack) then
           Field%unpkd_missval_in = Field%unpkd_missval_in*Field%pack(1) + Field%pack(2)
       endif

       Field%skip = .false.

 end subroutine get_input_variable_info

!#######################################################################

 subroutine copy_variable_information ( name, ncid_in, ncid_out, Field )

   character(len=NF90_MAX_NAME), intent(in)    :: name
   integer          ,       intent(in)    :: ncid_in, ncid_out
   type(field_type) ,       intent(inout) :: Field
   logical :: answer
   integer :: rdim, zcoord, idim
   integer :: i, istat, ipack, ioffs, imiss, ifill
   integer :: dimids(4)
   real    :: mval, fval

   call get_input_variable_info ( name, ncid_in, Field )

   if (Field%skip) return
   Field%skip = .true.
   zcoord = 0

   ! set output missing value

     Field%do_miss_out = Field%do_miss_in
     Field%missval_out = Field%missval_in

     if (use_default_missing_value) then
         if (Field%xtype == NF90_REAL8) Field%missval_out = NF90_FILL_REAL8
         if (Field%xtype == NF90_REAL4) Field%missval_out = NF90_FILL_REAL4
         if (Field%xtype == NF90_INT)   Field%missval_out = NF90_FILL_INT
         if (Field%xtype == NF90_INT2)  Field%missval_out = NF90_FILL_INT2
     endif

     ! unpack output missing value (for comparison with unpacked data)
       Field%unpkd_missval_out = Field%missval_out
       if (Field%do_miss_out .and. Field%do_pack) then
           Field%unpkd_missval_out = Field%unpkd_missval_out*Field%pack(1) + Field%pack(2)
       endif

   ! do not process fields with one dimension
     if (Field%ndim == 1) return

   !---- check horizontal dimensions (again) ---
     if (Field%ndim >= 2) then
         if (Field%dimlen(1) /= nlon .or. Field%dimlen(2) /= nlat) then
               if (Field%ndim == 2 .and. Field%static) return
               if (Field%ndim >= 3)                    return
         endif
     endif

   ! set dimension identifiers for output file

     if (Field%dimlen(2) == nlat)   dimids = dimids_mass

   !---- copy input variable metadata to output file ----

     select case (Field%ndim)
       case (4)
            if (Field%dimlen(3) == nlev)   zcoord = 1
            if (Field%dimlen(3) == nlev+1) zcoord = 2
            if (zcoord == 0) return
            if (Field%static) return   ! 4d static fields not allowed
            call copy_variable ( ncid_in, trim(name), ncid_out, &
                                 dimids, Field%varid_out,       &
                                 missvalue=Field%missval_out  )
       case (3)
          if (Field%static) then
            if (Field%dimlen(3) == nlev)   zcoord = 1
            if (Field%dimlen(3) == nlev+1) zcoord = 2
            if (zcoord == 0) return
            return               ! no 3d static fields allowed
            call copy_variable ( ncid_in, trim(name), ncid_out, &
                                 dimids(1:3), Field%varid_out,  &
                                 missvalue=Field%missval_out )
          else
             !---- 2d non-static field ----
             if (Field%do_miss_in) then
                call copy_variable ( ncid_in, trim(name), ncid_out,      &
                      (/dimids(1),dimids(2),dimids(4)/), Field%varid_out,&
                      missvalue=Field%missval_out )
             else
                call copy_variable ( ncid_in, trim(name), ncid_out,      &
                      (/dimids(1),dimids(2),dimids(4)/), Field%varid_out )
             endif
          endif
       case(2)
          !---- 2d static field ----
          if (Field%static) then
             if (Field%do_miss_in) then
                call copy_variable ( ncid_in, trim(name), ncid_out, &
                                     dimids(1:2), Field%varid_out,  &
                                     missvalue=Field%missval_out    )
             else
                call copy_variable ( ncid_in, trim(name), ncid_out, &
                                     dimids(1:2), Field%varid_out   )
             endif
          endif
       case default
            return
     end select

     Field%skip = .false.
     Field%y_axis = Field%dimlen(2)
     Field%z_axis = zcoord

 end subroutine copy_variable_information

!#######################################################################

 subroutine copy_axis_information (ncid_in, name_in, ncid_out, decomp, &
                                   Axis, dimid_out, name_out, name_edges)

 integer,          intent(in)    :: ncid_in, ncid_out, decomp(4)
 character(len=*), intent(in)    :: name_in
 type(axis_type),  intent(inout) :: Axis
 integer,          intent(out), optional :: dimid_out
 character(len=*), intent(out), optional :: name_out, name_edges

 integer :: varid_in, varid_out, dimlen_out, dimlen_in, did_out

   if (verbose > 2) print *, 'COPYING, name_in=',trim(name_in)

   if (decomp(1) == decomp(3) .and. decomp(2) == decomp(4)) then
      ! global domain indices == compute domain indices
      call copy_axis ( ncid_in, name_in, ncid_out, did_out, varid_in, varid_out, &
                       dimlen_out, name_out=name_out, name_edges=name_edges)
   else
      dimlen_in = decomp(4)-decomp(3)+1
      call copy_axis ( ncid_in, name_in, ncid_out, did_out, varid_in, varid_out, dimlen_out, &
                       dimlen_in=dimlen_in, name_out=name_out, name_edges=name_edges)
   endif

   Axis%name_in   = name_in
   Axis%name_out  = name_in
   if (present(name_out)) Axis%name_out  = name_out
   Axis%varid_in  = varid_in 
   Axis%varid_out = varid_out
   Axis%dimlen    = dimlen_out
   Axis%start     = (/decomp(3)/)

   if (present(dimid_out)) dimid_out = did_out

 end subroutine copy_axis_information

!#######################################################################

 subroutine set_up_input_mask (ncid, Field, data, mask, kbot)

   integer,          intent(in)  :: ncid
   type(field_type), intent(in)  :: Field
   real,             intent(in)  :: data(:,:,:)
   real,             intent(out) :: mask(:,:,:)
   integer,          intent(out) :: kbot(:,:)

   integer :: k, nlev

!     ----- set up mask and bottom level index -----

      nlev = size(data,3)

      kbot = nlev
      mask = 1.0

      if (Field%do_miss_in) then
         where ( data == Field%unpkd_missval_in )
                 mask = 0.0
         elsewhere
                 mask = 1.0
         endwhere

         do k = nlev, 1, -1
            if (count(mask(:,:,k) < 0.5) == 0) exit
            where (mask(:,:,k) < 0.5) kbot(:,:) = k-1
         enddo
      endif

 end subroutine set_up_input_mask

!#######################################################################

 subroutine read_and_set_up_mask (ncid, start, mask, kbot)

   integer,           intent(in)  :: ncid, start(4)
   real,              intent(out) :: mask(:,:,:)
   integer,           intent(out) :: kbot(:,:)

   character(len=NF90_MAX_NAME) :: name
   real,  dimension(size(mask,1),size(mask,2),size(mask,3)) :: values
   integer :: ivar, nvar, istat
   logical :: done
   type(field_type) :: Field

     done = .false.

   ! get number of variables and time axis id
     istat = NF90_INQUIRE (ncid, nvariables=nvar)
     if (istat /= NF90_NOERR) call error_handler (ncode=istat)

     do ivar = 1, nvar
        istat = NF90_INQUIRE_VARIABLE ( ncid, ivar, name=name )
        if (istat /= NF90_NOERR) call error_handler (ncode=istat)
        call get_input_variable_info ( name, ncid, Field )

      ! find first 3d mass field
        if (Field%ndim <= 2) cycle
        if (Field%ndim <= 3 .and. .not.Field%static) cycle
        if (Field%ndim == 4 .and.      Field%static) cycle

      ! must have mass dimension lengths and full vertical levels
        if (Field%dimlen(1) /= nlon .or. Field%dimlen(2) /= nlat .or. Field%dimlen(3) /= nlev) cycle

      ! must have input missing value
        if (.not.Field%do_miss_in) cycle

        call read_variable     (ncid, Field, values, start )
        call set_up_input_mask (ncid, Field, values, mask, kbot)
        done = .true.
        exit

     enddo

     if (.not.done) then
        mask = 1.0
        kbot = nlev
     endif

 end subroutine read_and_set_up_mask

!#######################################################################

 subroutine read_write_static (ncid_in, ncid_out, is, ie, js, je, Field)
   integer         , intent(in) :: ncid_in, ncid_out, is, ie, js, je
   type(field_type), intent(in) :: Field

   real :: values(is:ie,js:je,1)

 ! read and write static fields that are on the mass grid (nlon x nlat)

 ! read data
   call read_variable ( ncid_in, Field, values, start=(/is, js/))

 ! write data
   call write_variable ( ncid_out, Field, values)

 end subroutine read_write_static

!#######################################################################

 subroutine read_variable ( ncid, Field, values, start )
 integer,           intent(in)  :: ncid
 type(field_type),  intent(in)  :: Field
 real,              intent(out) :: values(:,:,:)
 integer, optional, intent(in)  :: start(:)

 integer :: ndim, nd, begin(4)

 ! number of dimensions (less the time dimension)
   ndim = Field%ndim
   if (.not.Field%static) ndim = ndim-1

 ! start indices
   nd    = Field%ndim
   begin = 1
   if (present(start)) then
       nd = size(start)
       begin(1:nd) = start
   endif

   select case (ndim)
      case(1)
          istat = NF90_GET_VAR (ncid, Field%varid_in, values(:,1,1), start=begin(1:nd))
          if (istat /= NF90_NOERR) call error_handler ('getting 1d var='//trim(Field%name), ncode=istat)
      case(2)
          istat = NF90_GET_VAR (ncid, Field%varid_in, values(:,:,1), start=begin(1:nd))
          if (istat /= NF90_NOERR) call error_handler ('getting 2d var='//trim(Field%name), ncode=istat)
      case(3)
          istat = NF90_GET_VAR (ncid, Field%varid_in, values(:,:,:), start=begin(1:nd))
          if (istat /= NF90_NOERR) call error_handler ('getting 3d var='//trim(Field%name), ncode=istat)
      case default
          call error_handler ('invalid number of dimensions in read_variable, var='//trim(Field%name))
   end select

 ! unpack values
   if (Field%do_pack) then
       if (Field%do_miss_in) then
           where (values .ne. Field%missval_in)
                 values = values*Field%pack(1) + Field%pack(2)
           elsewhere
                 values = Field%unpkd_missval_in
           endwhere
       else
           values = values*Field%pack(1) + Field%pack(2)
       endif
   endif

 end subroutine read_variable

!#######################################################################

 subroutine write_variable ( ncid, Field, values, start )
 integer,           intent(in)    :: ncid
 type(field_type),  intent(in)    :: Field
 real,              intent(inout) :: values(:,:,:)
 integer, optional, intent(in)    :: start(:)

 integer :: ndim, nd, begin(4)

 ! pack values (overwrite input data)
   if (Field%do_pack) then
       if (Field%do_miss_out) then
           where (values .ne. Field%unpkd_missval_out)
               values = nint((values-Field%pack(2))/Field%pack(1))
           elsewhere
               values = Field%missval_out
           endwhere
       else
           values = nint((values-Field%pack(2))/Field%pack(1))
       endif
   endif

 ! number of dimensions (less the time dimension)
   ndim = Field%ndim
   if (.not.Field%static) ndim = ndim-1

 ! start indices
   nd    = Field%ndim
   begin = 1
   if (present(start)) then
       nd = size(start)
       begin(1:nd) = start
   endif

   select case (ndim)
      case(1)
          istat = NF90_PUT_VAR (ncid, Field%varid_out, values(:,1,1), start=begin(1:nd))
          if (istat /= NF90_NOERR) call error_handler ('putting 1d var='//trim(Field%name), ncode=istat)
      case(2)
          istat = NF90_PUT_VAR (ncid, Field%varid_out, values(:,:,1), start=begin(1:nd))
          if (istat /= NF90_NOERR) call error_handler ('putting 2d var='//trim(Field%name), ncode=istat)
      case(3)
          istat = NF90_PUT_VAR (ncid, Field%varid_out, values(:,:,:), start=begin(1:nd))
          if (istat /= NF90_NOERR) call error_handler ('putting 3d var='//trim(Field%name), ncode=istat)
      case default
          call error_handler ('invalid number of dimensions in write_variable, var='//trim(Field%name))
   end select

 end subroutine write_variable

!#######################################################################

 subroutine compute_pres_full (pk, bk, ph, pf)
 real, intent(in)  :: pk(:), bk(:), ph(:,:,:)
 real, intent(out) :: pf(:,:,:)
 real, dimension(size(pf,1),size(pf,2))            :: dp, lpf
 real, dimension(size(pf,1),size(pf,2),size(ph,3)) :: lph
 integer :: k, nlev

 nlev = size(pf,3)
 
 ! compute p*logp at half levels
   where (ph(:,:,1) > 0.)
      lph(:,:,1) = ph(:,:,1) * log(ph(:,:,1))
   elsewhere
      lph(:,:,1) = 0.0
   endwhere
   lph(:,:,2:nlev+1) = ph(:,:,2:nlev+1) * log(ph(:,:,2:nlev+1))

! compute pressure at full levels
  do k = 1, nlev
      dp(:,:) = ph(:,:,k+1)-ph(:,:,k)
      lpf(:,:) = (lph(:,:,k+1)-lph(:,:,k))/dp(:,:) - 1.0
      pf(:,:,k) = exp((lph(:,:,k+1)-lph(:,:,k))/dp(:,:) - 1.0)
  enddo

 end subroutine compute_pres_full

!#######################################################################

 subroutine compute_height ( zsurf, temp, sphum, pfull, phalf, zfull, mask)
 real, intent(in),  dimension(:,:)   :: zsurf
 real, intent(in),  dimension(:,:,:) :: temp, sphum, pfull, phalf
 real, intent(out), dimension(:,:,:) :: zfull
 real, intent(in),  optional         :: mask(:,:,:)
 real, dimension(size(temp,1),size(temp,2)) ::  wta, wtb, vt, zb, zt, &
                                                lpb, lpt, lpf
 integer :: k
 logical :: flag

! flag for zero ptop
flag = count(phalf(:,:,1)<=0.) > 0

! vertical integration (bottom to top)
 zb(:,:) = zsurf(:,:)*GRAV
 lpb(:,:) = log(phalf(:,:,size(zfull,3)+1))
 do k = size(zfull,3), 1, -1
    lpf = log(pfull(:,:,k))
    wtb = lpb - lpf
    ! check for ptop=0
    if (k == 1 .and. flag) then
       wta = wtb
    else
       lpt = log(phalf(:,:,k))
       wta = lpf - lpt
    endif

    vt = temp(:,:,k)*(1.0+d608*sphum(:,:,k))*RDGAS ! virtual temperature * Rgas
    zt = zb + vt*(wta+wtb)
    zfull(:,:,k) = (zb + vt*wtb)/GRAV
    if (present(mask)) then
       where (mask(:,:,k) < 0.01)
          zt = GRAV*tref*(1.0-(pfull(:,:,k)/pref)**rgog)/tlapse
          zfull(:,:,k) = 0.5*(zt+zb)/GRAV
       endwhere
    endif
    zb  = zt
    lpb = lpt
 enddo

 end subroutine compute_height

!#######################################################################

end program run_pressure_interp

