#include <misc.h>
#include <params.h>
#if ( defined SCAM )
#include <max.h>
#endif

module history 52,16
!----------------------------------------------------------------------- 
! 
! Purpose: History module.  Contains data and functions for writing history files.
!
! Public functions/subroutines:
!   addfld, add_default
!   intht
!   write_restart_history
!   read_restart_history
!   outfld
!   wshist
!   scm_intht
!   scm_histfield_ini
!   scm_addfield
!   initialize_iop_history
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
! $Id: history.F90,v 1.26.2.48.4.1 2004/05/20 18:36:01 mvr Exp $
!-----------------------------------------------------------------------
   use shr_kind_mod, only: r8 => shr_kind_r8, r4 => shr_kind_r4
   use ppgrid,    only: pcols
   use constituents, only: pcnst, pnats, cnst_name, cnst_longname, &
                           dcconnam, sflxnam, cnst_get_ind
   use filenames, only: mss_wpass, mss_irt, interpret_filename_spec, get_archivedir
   use abortutils, only: endrun
#if ( defined STAGGERED )
   use pmgrid,    only: masterproc, beglat, endlat, plat, plon, plev, plevp, dyngrid_set, splon, beglev, endlev, endlevp
#else
   use pmgrid,    only: masterproc, beglat, endlat, plat, plon, plev, plevp, dyngrid_set
#endif
   use icarus_scops, only: ntau, npres, prlim, taulim
#if ( defined SCAM )
  use scamMod, only :initlonidx,initlatidx,divq3d,divt3d
#endif

   implicit none

PRIVATE

   include 'netcdf.inc'

   integer, parameter :: pflds = 1000           ! max number of fields 
   integer, parameter :: ptapes = 6             ! max number of tapes
   integer, parameter :: max_chars = 128        ! max chars for char variables

   real(r8), parameter :: fillvalue = 1.e36     ! fill value for reduced grid and others

   type field_info
      character*8 :: name                       ! field name
      character*(max_chars) :: long_name        ! long name
      character*(max_chars) :: units            ! units

      logical :: flag_xyfill                    ! non-applicable xy points flagged with fillvalue
      logical :: flag_isccplev                  ! levels dimension is ISCCP not CAM

      integer :: coldimin                       ! column dimension of model array
      integer :: numlev                         ! vertical dimension (.nc file and internal arr)
      integer :: begver                         ! on-node vert start index
      integer :: endver                         ! on-node vert end index
      integer :: begdim3                        ! on-node chunk or lat start index
      integer :: enddim3                        ! on-node chunk or lat end index
      integer :: decomp_type                    ! type of decomposition (physics or dynamics)
      integer, pointer :: colperdim3(:)         ! number of valid elements per chunk or lat
   end type field_info
!
! master_entry: elements of an entry in the master field list
!
   type master_entry
      type (field_info)     :: field            ! field information
      character*1           :: avgflag(ptapes)  ! averaging flag
      character*(max_chars) :: time_op(ptapes)  ! time operator (e.g. max, min, avg)
      logical               :: act_sometape     ! Field is active on some tape
      logical               :: actflag(ptapes)  ! Per tape active/inactive flag
      integer               :: htapeindx(ptapes)! This field's index on particular history tape
   end type master_entry

   type (master_entry) :: masterlist(pflds)     ! master field list
!
! hbuffer_2d, hbuffer_3d: 2-D and 3-D history buffer pointers.
!     Select either r4 or r8 kind buffer depending on hbuf_prec.
!
   type hbuffer_2d
      real(r8), pointer :: buf8(:,:)            ! 2-D history buffer for r8
      real(r4), pointer :: buf4(:,:)            ! 2-D history buffer for r4
   end type hbuffer_2d

   type hbuffer_3d
      real(r8), pointer :: buf8(:,:,:)          ! 3-D history buffer for r8
      real(r4), pointer :: buf4(:,:,:)          ! 3-D history buffer for r4
   end type hbuffer_3d
!
! arrays served as targets for history pointers
!
   integer,  target :: nothing_int(1,1)         ! 2-D integer target
   real(r8), target :: nothing_r8(1,1,1)        ! 3-D r8 target
   real(r4), target :: nothing_r4(1,1,1)        ! 3-D r4 target

   type column_info
      character*(max_chars) :: lat_name ! latitude name for this column or columns
      character*(max_chars) :: lon_name ! latitude name for this column or columns
      integer :: num_lats            ! number of lats in a group of contiguous columns
      integer :: num_lons            ! number of lons in a group of contiguous columns
      integer :: columnlat(2)       ! beginning and ending latitude (range) dimensioned by groups
      integer :: columnlon(2)       ! beginning and ending longitude (range) dimensioned by groups
   end type column_info
!
! hentry: elements of an entry in the list of active fields on a single history file
!
   type hentry
      type (field_info)     :: field            ! field information
      character*1           :: avgflag          ! averaging flag
      character*(max_chars) :: time_op          ! time operator (e.g. max, min, avg)
      character*(max_chars),pointer :: field_column_name(:) ! names of column groups written to tape

      integer :: hbuf_prec                      ! history buffer precision
      integer :: hwrt_prec                      ! history output precision

      type (hbuffer_3d)   :: hbuf               ! history buffer
      integer, pointer :: nacs(:,:)             ! accumulation counter
   end type hentry
!
! active_entry: vehicle for producing a ragged array
!
   type active_entry
      type (hentry) :: hlist(pflds)             ! array of history tape entries
      type (column_info),pointer :: column(:)             ! array of history tape entries
   end type active_entry

   type (active_entry) :: tape(ptapes)          ! history tapes
!
! dim_index_2d, dim_index_3d: 2-D & 3-D dimension index lower & upper bounds
!
   type dim_index_2d                   ! 2-D dimension index
      integer :: beg1, end1            ! lower & upper bounds of 1st dimension
      integer :: beg2, end2            ! lower & upper bounds of 2nd dimension
   end type dim_index_2d

   type dim_index_3d                   ! 3-D dimension index
      integer :: beg1, end1            ! lower & upper bounds of 1st dimension
      integer :: beg2, end2            ! lower & upper bounds of 2nd dimension
      integer :: beg3, end3            ! lower & upper bounds of 3rd dimension
   end type dim_index_3d

   integer :: ndm(12)                           ! number of days in each month (jan-dec)
   save ndm
   data ndm/31,28,31,30,31,30,31,31,30,31,30,31/

   integer :: nfmaster = 0             ! number of fields in master field list
   integer :: nflds(ptapes)            ! number of fields per tape

! per tape sampling frequency (0=monthly avg)

   integer :: i                        ! index for nhtfrq initialization
   integer :: nhtfrq(ptapes) = (/0, (-24, i=2,ptapes)/)  ! history write frequency (0 = monthly)
   integer :: mfilt(ptapes) = 30       ! number of time samples per tape
   integer :: nfils(ptapes)            ! Array of no. of files on current h-file
   integer :: ngroup(ptapes)           ! Array of no. of contiguous columns on current h-file
   integer :: mtapes = 0               ! index of max history file requested 
   integer :: nexcl(ptapes)            ! Actual number of excluded fields
   integer :: nincl(ptapes)            ! Actual number of included primary file fields
   integer :: nhstpr(ptapes) = 8       ! history buffer precision (8 or 4 bytes)
   integer :: ndens(ptapes) = 2        ! packing density (nf_float vs nf_double)
   integer :: ncprec(ptapes) = -999    ! netcdf packing parameter based on ndens
   real(r8) :: beg_time(ptapes)        ! time at beginning of an averaging interval
!
! Netcdf ids
!
   integer :: nfid(ptapes)             ! file id
   integer :: varid(pflds,ptapes)      ! variable ids
   integer :: mdtid(ptapes)            ! var id for timestep
   integer :: ndbaseid(ptapes)         ! var id for base day
   integer :: nsbaseid(ptapes)         ! var id for base seconds of base day
   integer :: nbdateid(ptapes)         ! var id for base date
   integer :: nbsecid(ptapes)          ! var id for base seconds of base date
   integer :: ndcurid(ptapes)          ! var id for current day
   integer :: nscurid(ptapes)          ! var id for current seconds of current day
   integer :: dateid(ptapes)           ! var id for current date
   integer :: co2vmrid(ptapes)         ! var id for co2 volume mixing ratio
   integer :: datesecid(ptapes)        ! var id for curent seconds of current date
#if ( defined BFB_CAM_SCAM_IOP )
   integer :: bdateid(ptapes)         ! var id for base date
   integer :: tsecid(ptapes)        ! var id for curent seconds of current date
#endif
   integer :: nstephid(ptapes)         ! var id for current timestep
   integer :: timeid(ptapes)           ! var id for time
   integer :: tbndid(ptapes)           ! var id for time_bnds
   integer :: gwid(ptapes)             ! var id for gaussian weights
   integer :: date_writtenid(ptapes)   ! var id for date time sample written
   integer :: time_writtenid(ptapes)   ! var id for time time sample written
   integer :: nlonid(ptapes)           ! var id for number of longitudes
   integer :: wnummaxid(ptapes)        ! var id for cutoff fourier wavenumber (reduced grid)
   integer :: nscurf(ptapes)           ! First "current" second of day for each h-file
   integer :: ncsecf(ptapes)           ! First "current" second of date for each h-file

   logical :: rgnht(ptapes) = .false.  ! flag array indicating regeneration volumes
   logical :: hstwr(ptapes) = .false.  ! Flag for history writes
   logical :: empty_htapes  = .false.  ! Namelist flag indicates no default history fields
   logical :: htapes_defined = .false. ! flag indicates history contents have been defined

   integer, parameter :: nlen = 256    ! Length of strings
   character(len=nlen) :: hrestpath(ptapes) = (/(' ',i=1,ptapes)/) ! Full history restart pathnames
   character(len=nlen) :: nfpath(ptapes) = (/(' ',i=1,ptapes)/) ! Array of first pathnames, for header
   character(len=nlen) :: cpath(ptapes)                   ! Array of current pathnames
   character(len=nlen) :: nhfil(ptapes)                   ! Array of current file names
   character(len=1)  :: avgflag_pertape(ptapes) = (/(' ',i=1,ptapes)/) ! per tape averaging flag
   character(len=8)  :: logname             ! user name
   character(len=16) :: host                ! host name
   character(len=nlen) :: ctitle = ' '      ! Case title
   character(len=8)  :: inithist = 'YEARLY' ! If set to '6-HOURLY, 'DAILY', 'MONTHLY' or
                                            ! 'YEARLY' then write IC file 
   character(len=10) :: fincl(pflds,ptapes) ! List of fields to add to primary h-file
   character(len=max_chars) :: fincllonlat(pflds,ptapes) ! List of fields to add to primary h-file
   character(len=8)  :: fexcl(pflds,ptapes) ! List of fields to rm from primary h-file
   character(len=10) :: fhstpr(pflds,ptapes) ! List of fields to change default hbuf size
   character(len=10) :: fwrtpr(pflds,ptapes) ! List of fields to change default history output prec
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Hashing.
!
!  Accelerate outfld processing by using a hash function of the field name
!  to index masterlist and determine whehter the particular field is to
!  be written to any history tape.
!
!
!  Note: the outfld hashing logic will fail if any of the following are true:
!
!         1) The lower bound on the dimension of 'masterlist' is less than 1.
!
!         2) 'outfld' is called with field names that are not defined on
!            masterlist.  This applies to both initial/branch and restart
!            runs.
!
!         3) An inconsistency between a field's tape active flag
!            'masterlist(ff)%actflag(t)' and active fields read from
!            restart files.
!
!         4) Invoking function 'gen_hash_key' before the primary and secondary
!            hash tables have been created (routine bld_outfld_hash_tbls).
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!  User definable constants for hash and overflow tables.
!  Define size of primary hash table (specified as 2**size).
!
   integer, parameter :: tbl_hash_pri_sz_lg2 = 12
!
!  Define size of overflow hash table % of primary hash table.
!
   integer, parameter :: tbl_hash_oflow_percent = 5
!
!  Do *not* modify the parameters below.
!
   integer, parameter :: tbl_hash_pri_sz = 2**tbl_hash_pri_sz_lg2
   integer, parameter :: tbl_hash_oflow_sz = tbl_hash_pri_sz * (tbl_hash_oflow_percent/100.0) 
!
!  The primary and overlow tables are organized to mimimize space (read:
!  try to maximimze cache line usage).
!
!  gen_hash_key(fieldname) will return an index on the interval
!  [0 ... tbl_hash_pri_sz-1].
!
!
!  Primary:
!  gen_hash_key(fieldname)-------+     +----------+
!                                |     |   -ii    | 1 ------>tbl_hash_oflow(ii)
!                                |     +----------+
!                                +-->  |    ff    | 2 ------>masterlist(ff)
!                                      +----------+
!                                      |          | ...
!                                      +----------+
!                                      |          | tbl_hash_pri_sz
!                                      +----------+
!
!  Overlow (if tbl_hash_pri() < 0):
!  tbl_hash_pri(gen_hash_key(fieldname))
!                         |
!                         |            +----------+
!                         |            |     1    | 1  (one entry on O.F. chain)
!                         |            +----------+
!                         |            |    ff_m  | 2
!                         |            +----------+
!                         +--------->  |     3    | 3  (three entries on chain)
!                                      +----------+
!                                      |    ff_x  | 4
!                                      +----------+
!                                      |    ff_y  | 5
!                                      +----------+
!                                      |    ff_z  | 6
!                                      +----------+
!                                      |          | ...
!                                      +----------+
!                                      |          | tbl_hash_oflow_sz
!                                      +----------+
!
!
   integer, dimension(0:tbl_hash_pri_sz-1) :: tbl_hash_pri ! Primary hash table
   integer, dimension(tbl_hash_oflow_sz) :: tbl_hash_oflow ! Overlow hash table
!
!  Constants used in hashing function gen_hash_key.
!  Note: if the constants in table 'tbl_gen_hash_key' below are modified,
!        changes are required to routine 'gen_hash_key' because of specific
!        logic in the routine that optimizes character strings of length 8.
!

   integer, parameter :: gen_hash_key_offset = z'000053db'

   integer, dimension(16), save :: tbl_gen_hash_key = &
   (/61,59,53,47,43,41,37,31,29,23,17,13,11,7,3,1/)

!
! Equivalence to please namelist on a wide variety of platforms
! NOTE: It is *ASSUMED* that ptapes is 6
!
   character*10 fincl1(pflds)
   character*10 fincl2(pflds)
   character*10 fincl3(pflds)
   character*10 fincl4(pflds)
   character*10 fincl5(pflds)
   character*10 fincl6(pflds)

   equivalence (fincl1,fincl(1,1))
   equivalence (fincl2,fincl(1,2))
   equivalence (fincl3,fincl(1,3))
   equivalence (fincl4,fincl(1,4))
   equivalence (fincl5,fincl(1,5))
   equivalence (fincl6,fincl(1,6))

   character(len=max_chars) fincl1lonlat(pflds)
   character(len=max_chars) fincl2lonlat(pflds)
   character(len=max_chars) fincl3lonlat(pflds)
   character(len=max_chars) fincl4lonlat(pflds)
   character(len=max_chars) fincl5lonlat(pflds)
   character(len=max_chars) fincl6lonlat(pflds)

   equivalence (fincl1lonlat,fincllonlat(1,1))
   equivalence (fincl2lonlat,fincllonlat(1,2))
   equivalence (fincl3lonlat,fincllonlat(1,3))
   equivalence (fincl4lonlat,fincllonlat(1,4))
   equivalence (fincl5lonlat,fincllonlat(1,5))
   equivalence (fincl6lonlat,fincllonlat(1,6))

   character*8 fexcl1(pflds)
   character*8 fexcl2(pflds)
   character*8 fexcl3(pflds)
   character*8 fexcl4(pflds)
   character*8 fexcl5(pflds)
   character*8 fexcl6(pflds)

   equivalence (fexcl1,fexcl(1,1))
   equivalence (fexcl2,fexcl(1,2))
   equivalence (fexcl3,fexcl(1,3))
   equivalence (fexcl4,fexcl(1,4))
   equivalence (fexcl5,fexcl(1,5))
   equivalence (fexcl6,fexcl(1,6))

   character*10 fhstpr1(pflds)
   character*10 fhstpr2(pflds)
   character*10 fhstpr3(pflds)
   character*10 fhstpr4(pflds)
   character*10 fhstpr5(pflds)
   character*10 fhstpr6(pflds)

   equivalence (fhstpr1,fhstpr(1,1))
   equivalence (fhstpr2,fhstpr(1,2))
   equivalence (fhstpr3,fhstpr(1,3))
   equivalence (fhstpr4,fhstpr(1,4))
   equivalence (fhstpr5,fhstpr(1,5))
   equivalence (fhstpr6,fhstpr(1,6))

   character*10 fwrtpr1(pflds)
   character*10 fwrtpr2(pflds)
   character*10 fwrtpr3(pflds)
   character*10 fwrtpr4(pflds)
   character*10 fwrtpr5(pflds)
   character*10 fwrtpr6(pflds)

   equivalence (fwrtpr1,fwrtpr(1,1))
   equivalence (fwrtpr2,fwrtpr(1,2))
   equivalence (fwrtpr3,fwrtpr(1,3))
   equivalence (fwrtpr4,fwrtpr(1,4))
   equivalence (fwrtpr5,fwrtpr(1,5))
   equivalence (fwrtpr6,fwrtpr(1,6))
!
! Overloading assignment operator
!

   interface assignment (=)
      module procedure hbuf_assigned_to_hbuf
      module procedure hbuf_assigned_to_real8
   end interface
!
! Generic procedures
!

   interface allocate_hbuf 9
      module procedure allocate_hbuf2d
      module procedure allocate_hbuf3d
   end interface


   interface deallocate_hbuf 7
      module procedure deallocate_hbuf2d
      module procedure deallocate_hbuf3d
   end interface


   interface nullify_hbuf 4
      module procedure nullify_hbuf2d
      module procedure nullify_hbuf3d
   end interface
!
! Public entities
!

!
! Filename specifiers for history, initial files and restart history files
! (%c = caseid, $y = year, $m = month, $d = day, $s = seconds in day, %t = tape number)
!
   character(len=256) :: ifilename_spec = '%c.cam2.i.%y-%m-%d-%s.nc'  ! Initial files
   character(len=256) :: rhfilename_spec = '%c.cam2.rh%t.%y-%m-%d-%s' ! history restart
   character(len=256), public :: hfilename_spec(ptapes) = (/ (' ', i=1, ptapes) /) ! filename specifyer

! Needed by anyone calling addfld

   integer, parameter, public :: phys_decomp = 1     ! flag indicates physics decomposition
   integer, parameter, public :: dyn_decomp  = 2     ! flag indicates dynamics decomposition

! To allow parameterizations to initialize arrays to the fillvalue
! THIS NEEDS TO BE FIXED.  No parameterization should be allowed access to fillvalue

   public :: fillvalue

! Needed by cam

   public :: bldfld

! Needed by initext

   public :: nhtfrq, mfilt, inithist, ctitle

! Needed by read_namelist

   public :: fincl, fincl1, fincl2, fincl3, fincl4, fincl5, fincl6
   public :: fincllonlat, fincl1lonlat, fincl2lonlat, fincl3lonlat, fincl4lonlat, fincl5lonlat, fincl6lonlat
   public :: fexcl, fexcl1, fexcl2, fexcl3, fexcl4, fexcl5, fexcl6
   public :: fhstpr, fhstpr1, fhstpr2, fhstpr3, fhstpr4, fhstpr5, fhstpr6
   public :: fwrtpr, fwrtpr1, fwrtpr2, fwrtpr3, fwrtpr4, fwrtpr5, fwrtpr6
   public :: pflds, ptapes, empty_htapes, nhstpr, ndens
   public :: avgflag_pertape

! Needed by stepon

   public :: hstwr
   public :: nfils

! Functions

#if ( defined SCAM )
   public :: scm_intht
   public :: scm_histfield_ini
   public :: scm_addfield
#endif
#if ( defined BFB_CAM_SCAM_IOP || defined SCAM )
   public :: initialize_iop_history
#endif
   public :: write_restart_history     ! Write restart history data
   public :: read_restart_history      ! Read restart history data
   public :: wshist                    ! Write files out
   public :: write_inithist            ! Write the initial file
   public :: outfld                    ! Output a field
   public :: intht                     ! Initialization
   public :: wrapup                    ! Archive history files at end of run
   public :: addfld                    ! Add a field to history file
   public :: add_default               ! Add the default fields
   public :: get_hfilepath             ! Return history filename
   public :: get_mtapes                ! Return the number of tapes being used
   public :: get_hist_restart_filepath ! Return the full filepath to the history restart file

CONTAINS


   subroutine intht () 2,12
!----------------------------------------------------------------------- 
! 
! Purpose: Initialize history file handler for initial or continuation run.
!          For example, on an initial run, this routine initializes "mtapes"
!          history files.  On a restart or regeneration  run, this routine 
!          only initializes history files declared beyond what existed on the 
!          previous run.  Files which already existed on the previous run have 
!          already been initialized (i.e. named and opened) in routine RESTRT.
! 
! Method: Loop over tapes and fields per tape setting appropriate variables and
!         calling appropriate routines
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
      use ioFileMod
      use shr_sys_mod, only: shr_sys_getenv
      use time_manager, only: get_curr_date, get_curr_time
#if ( defined SPMD )
      use mpishorthand
#endif
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
!
! Local workspace
!
      integer :: t, f              ! tape, field indices
      integer :: coldimin          ! column dimension of model array
      integer :: begver            ! on-node vert start index
      integer :: endver            ! on-node vert end index
      integer :: begdim3           ! on-node chunk or lat start index
      integer :: enddim3           ! on-node chunk or lat end index
      type (dim_index_3d) :: dimind  ! 3-D dimension index
      integer :: day, sec          ! day and seconds from base date
      integer :: i                 ! index
      integer :: rcode             ! shr_sys_getenv return code
!
! Get users logname and machine hostname
!
      if ( masterproc )then
         logname = ' '
         call shr_sys_getenv ('LOGNAME',logname,rcode)
         if (rcode /= 0) then
            call endrun ('INTHT: Cannot find LOGNAME environment variable')
         end if
         host = ' '
         call shr_sys_getenv ('HOST',host,rcode)
      end if
!
! Set default history history contents
!
      call h_default ()
!
! Override averaging flag for all fields on a particular tape if namelist input so specifies
!
      do t=1,ptapes
         if (avgflag_pertape(t) /= ' ') then
            call h_override (t)
         end if
      end do
!
! Define field list information for all history files.  
! Update mtapes to reflect *current* number of history files (note, 
! restart and regen runs can have additional auxiliary history files
! declared).
!
      call fldlst ()
!
! Loop over max. no. of history files permitted  
!
      call get_curr_time(day, sec)  ! elapased time since reference date
      do t=1,mtapes
         nfils(t) = 0            ! no. of time samples in hist. file no. t

! Time a beginning of current averaging interval.

         beg_time(t) = day + sec/86400._r8

      end do
!
! Check that the number of history files declared does not exceed
! the maximum allowed.
!
      if (mtapes > ptapes) then
         write(6,*) 'INTHT: Too many history files declared, max=',ptapes
         write(6,*)'To increase, change parameter ptapes.'
         call endrun
      end if
!
! Initialize history variables
!
      do t=1,mtapes
         do f=1,nflds(t)
            coldimin = tape(t)%hlist(f)%field%coldimin
            begver   = tape(t)%hlist(f)%field%begver
            endver   = tape(t)%hlist(f)%field%endver
            begdim3  = tape(t)%hlist(f)%field%begdim3
            enddim3  = tape(t)%hlist(f)%field%enddim3

            dimind = dim_index_3d (1,coldimin,begver,endver,begdim3,enddim3)
            call allocate_hbuf (tape(t)%hlist(f)%hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
            tape(t)%hlist(f)%hbuf = 0._r8

            allocate (tape(t)%hlist(f)%nacs(coldimin,begdim3:enddim3))
            do i = begdim3, enddim3
               tape(t)%hlist(f)%nacs(:coldimin,i) = 0
            end do
         end do
      end do
      
      return
   end subroutine intht

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


   subroutine write_restart_history (nrg, luhrest) 1,32

      use binary_io
      use ioFileMod
#ifdef SPMD
# if ( defined STAGGERED )
      use spmd_dyn, only: npes, compute_gsfactors, comm_y
      use pmgrid, only: myid_z, strip3dxzy, strip3dxzyp, strip2d
      use parutilitiesmodule, only: commglobal, pargatherreal, pargatherint, pargatherreal4
# else
      use spmd_dyn, only: npes, compute_gsfactors
# endif
      use mpishorthand
#endif
      use phys_grid, only: gather_chunk_to_field_int
      use dycore, only: dycore_is

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

!
! Arguments
!
      integer, intent(in) :: nrg        ! unit number
      integer, intent(in) :: luhrest    ! unit number
!
! Local workspace
!
      integer t,f                       ! Tape, field indices
      integer numlev                    ! number of vertical levels (dimension and loop)
      integer ioerr                     ! write error status
      integer coldimin                  ! column dimension of model array
      character(len=256) :: fname       ! History restart filename

#ifdef SPMD
      integer :: numsend                ! number of items to be sent
      integer :: numrecv(0:npes-1)      ! number of items to be received
      integer :: displs(0:npes-1)       ! displacement array
      integer :: numowned               ! number of items owned by this MPI task
      integer :: mpireal                ! MPI real data type
#endif

      type (hbuffer_3d) :: hbuf         ! full-field history buffer written by master
      integer, pointer :: fullnacs(:,:) ! full-field accumulation counter written by master
      type (dim_index_3d) :: dimind     ! 3-D dimension index
!
!-----------------------------------------------------------------------
! Write the history restart data if necessary
!-----------------------------------------------------------------------
!
      where (hstwr(:))
         rgnht(:) = .false.
      elsewhere
         rgnht(:) = .true.
      end where

      if (masterproc) then
         write (nrg) rgnht, mtapes, varid, fincl, fexcl
      end if
!
! If a history tape is not currently being disposed then write a history buffer restart file.
! History restart info: some f90 compilers (e.g. SGI) complain about I/O of derived types which 
! have pointer components, so explicitly write each one.
!
      do t=1,mtapes
         if (masterproc) then
            write (nrg, iostat=ioerr) nhtfrq(t), nflds(t), nfils(t), mfilt(t), &
                                      nfpath(t), cpath(t), nhfil(t),  &
                                      nhstpr(t), ndens(t), ncprec(t), beg_time(t)
            if (ioerr /= 0 ) then
               write (6,*) 'WRITE ioerror ',ioerr,' on i/o unit = ',nrg
               call endrun ('WRITE_RESTART_HISTORY')
            end if

            write (nrg,iostat=ioerr) ngroup(t)
            if (ngroup(t) .gt. 0) then 
               do f=1,ngroup(t)
                  write (nrg,iostat=ioerr) tape(t)%column(f)
               end do
            else
               write (nrg,iostat=ioerr) tape(t)%column(1)
            end if
            do f=1,nflds(t)
               write (nrg,iostat=ioerr) tape(t)%hlist(f)%field%name,        &
                                        tape(t)%hlist(f)%field%long_name,   &
                                        tape(t)%hlist(f)%field%units,       &
                                        tape(t)%hlist(f)%field%numlev,      &
                                        tape(t)%hlist(f)%field%decomp_type, &
                                        tape(t)%hlist(f)%avgflag,           &
                                        tape(t)%hlist(f)%time_op,           &
                                        tape(t)%hlist(f)%hbuf_prec,         &
                                        tape(t)%hlist(f)%hwrt_prec,         &
                                        tape(t)%hlist(f)%field%flag_xyfill, &
                                        tape(t)%hlist(f)%field%flag_isccplev 
               if (ioerr /= 0 ) then
                  write(6,*)'WRITE_RESTART_HISTORY: End or error condition writing ', &
                       'history restart field ',f,' from tape ',t
                  call endrun
               end if
            end do

         end if

         if (rgnht(t)) then
            if (masterproc) then
               fname = interpret_filename_spec( rhfilename_spec, number=(t-1) )
               hrestpath(t) = trim(get_archivedir('rest'))//fname
               write (nrg,iostat=ioerr) hrestpath(t)
               if (ioerr /= 0 ) then
                  call endrun ('WRITE_RESTART_HISTORY: End or error condition writing history restart filename')
               end if
               call opnfil (fname, luhrest, 'u')
            endif

            do f=1,nflds(t)
               coldimin =  tape(t)%hlist(f)%field%coldimin
               numlev   =  tape(t)%hlist(f)%field%numlev
#ifdef SPMD
               if (masterproc) then
                  dimind = dim_index_3d (1,plon,1,numlev,1,plat)
                  call allocate_hbuf (hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
                  allocate (fullnacs(plon,plat))
               else
                  call assoc_hbuf_with_nothing (hbuf,tape(t)%hlist(f)%hbuf_prec)
                  fullnacs => nothing_int
               end if

               if (tape(t)%hlist(f)%hbuf_prec == 8) then
                  mpireal = mpir8
               else
                  mpireal = mpir4
               end if

               select case (tape(t)%hlist(f)%field%decomp_type)
               case (phys_decomp)
                  call gather_chunk_to_field_hbuf (1, numlev, 1, plon, tape(t)%hlist(f)%hbuf, hbuf)
                  call gather_chunk_to_field_int (1, 1, 1, plon, tape(t)%hlist(f)%nacs, fullnacs)
               case (dyn_decomp)
                  if ( dycore_is('LR') )then
# if ( defined STAGGERED )
! NEW LR CODING
                     if (tape(t)%hlist(f)%hbuf_prec == 8) then
                        select case (numlev)
                        case (1)
                           if (myid_z .eq. 0) call pargatherreal(comm_y, 0,  &
                                  tape(t)%hlist(f)%hbuf%buf8, strip2d, hbuf%buf8)
                        case (plev)
                           call pargatherreal(commglobal, 0, tape(t)%hlist(f)%hbuf%buf8, &
                                       strip3dxzy, hbuf%buf8)
                        case (plevp)
                           call pargatherreal(commglobal, 0, tape(t)%hlist(f)%hbuf%buf8, &
                                       strip3dxzyp, hbuf%buf8)
                        case default
                           write(6,*)'WRITE_RESTART_HISTORY: invalid number of levels=', numlev
                           call endrun ()
                        end select
                        if (myid_z .eq. 0) call pargatherint(comm_y, 0, &
                           tape(t)%hlist(f)%nacs, strip2d, fullnacs)
                     else
                        select case (numlev)
                        case (1)
                           if (myid_z .eq. 0) call pargatherreal4(comm_y, 0,  &
                                  tape(t)%hlist(f)%hbuf%buf4, strip2d, hbuf%buf4)
                        case (plev)
                           call pargatherreal4(commglobal, 0, tape(t)%hlist(f)%hbuf%buf4, &
                                       strip3dxzy, hbuf%buf4)
                        case (plevp)
                           call pargatherreal4(commglobal, 0, tape(t)%hlist(f)%hbuf%buf4, &
                                       strip3dxzyp, hbuf%buf4)
                        case default
                           write(6,*)'WRITE_RESTART_HISTORY: invalid number of levels=', numlev
                           call endrun ()
                        end select
                        if (myid_z .eq. 0) call pargatherint(comm_y, 0, &
                           tape(t)%hlist(f)%nacs, strip2d, fullnacs)
                     endif
# endif
                  else
                     numowned = coldimin*numlev
                     call compute_gsfactors (numowned, numsend, numrecv, displs)
                     call mpigatherv_hbuf (tape(t)%hlist(f)%hbuf, numsend, mpireal, hbuf, numrecv, &
                                           displs, mpireal, 0, mpicom)

                     numowned = coldimin
                     call compute_gsfactors (numowned, numsend, numrecv, displs)
                     call mpigatherv (tape(t)%hlist(f)%nacs, numsend, mpiint, fullnacs, numrecv, &
                                      displs, mpiint, 0, mpicom)
                  endif
               end select

               if (masterproc) then
                  call write_hbuf (hbuf, luhrest, ioerr)
                  call deallocate_hbuf (hbuf)
                  write (luhrest,iostat=ioerr) fullnacs
                  deallocate (fullnacs)
               else
                  call nullify_hbuf (hbuf)
                  nullify (fullnacs)
               end if
#else
               select case (tape(t)%hlist(f)%field%decomp_type)
               case (phys_decomp)
                  dimind = dim_index_3d (1,plon,1,numlev,1,plat)
                  call allocate_hbuf (hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
                  call gather_chunk_to_field_hbuf (1, numlev, 1, plon, tape(t)%hlist(f)%hbuf, hbuf)
                  call write_hbuf (hbuf, luhrest, ioerr)
                  call deallocate_hbuf (hbuf)
                  allocate (fullnacs(plon,plat))
                  call gather_chunk_to_field_int (1, 1, 1, plon, tape(t)%hlist(f)%nacs, fullnacs)
                  write (luhrest,iostat=ioerr) fullnacs
                  deallocate (fullnacs)
               case (dyn_decomp)
                  call write_hbuf (tape(t)%hlist(f)%hbuf, luhrest, ioerr)
                  write (luhrest,iostat=ioerr) tape(t)%hlist(f)%nacs
               end select
#endif
            end do
            if (masterproc) then
               close (unit=luhrest)
               call putfil (fname, hrestpath(t), mss_wpass, mss_irt, .true.)
            end if
         end if
      end do

      return
   end subroutine write_restart_history

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


   subroutine read_restart_history (nrg, luhrest) 1,67

      use ppgrid,    only: begchunk, endchunk
      use phys_grid, only: get_ncols_p, scatter_field_to_chunk_int
      use rgrid,     only: nlon
      use dycore,    only: dycore_is
      use binary_io
      use ioFileMod
#ifdef SPMD
# if ( defined STAGGERED )
      use pmgrid, only: strip3dxzy, strip3dxzyp, strip2d
      use restart_dynamics, only: lrreadin, lrreadini, lrreadin4
# endif
      use mpishorthand
#endif
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
!
! Arguments
!
      integer, intent(in) :: nrg            ! unit number
      integer, intent(in) :: luhrest        ! unit number
!
! Local workspace
!
      integer t, f                     ! tape, field indices
      integer c                        ! chunk or lat index
      integer lenc                     ! length of useful character data
      integer numlev                   ! number of vertical levels (dimension and loop)
      integer begver                   ! on-node vert start index
      integer endver                   ! on-node vert end index
      integer ioerr                    ! error code from read()
      integer coldimin                 ! column dimension of model array
      integer begdim3                  ! on-node chunk or lat start index
      integer enddim3                  ! on-node chunk or lat end index
      integer ncol                     ! number of active columns per chunk
      integer lenarr                   ! global size of array to be read

      character(len=80)  :: locfn       ! Local filename

      integer, pointer :: nacs(:,:)    ! accumulation counter

      type (hbuffer_3d) :: hbuf             ! history buffer
      integer,  pointer :: fullnacs(:,:)    ! accumulation buffer
      type (dim_index_3d) :: dimind    ! 3-D dimension index

      if (masterproc) then

         read (nrg,iostat=ioerr) rgnht, mtapes, varid, fincl, fexcl
         if (ioerr /= 0 ) then
            write (6,*) 'READ ioerror ',ioerr,' on i/o unit = ',nrg
            call endrun ('READ_RESTART_HISTORY')
         end if

         do t=1,mtapes
            read (nrg,iostat=ioerr) nhtfrq(t), nflds(t), nfils(t), mfilt(t), &
                                    nfpath(t), cpath(t), nhfil(t), &
                                    nhstpr(t), ndens(t), ncprec(t), beg_time(t)
            if (ioerr /= 0) then
               write (6,*) 'READ ioerror ',ioerr,' on i/o unit = ',nrg
               call endrun ('READ_RESTART_HISTORY')
            end if

            read (nrg,iostat=ioerr) ngroup(t)
            if (ioerr /= 0) then
               write (6,*) 'READ ioerror for group number',ioerr,' on i/o unit = ',nrg
               call endrun ('READ_RESTART_HISTORY')
            end if

            if (ngroup(t) .gt. 0) then
               allocate(tape(t)%column(ngroup(t)))
               do f=1,ngroup(t)
                  read (nrg,iostat=ioerr) tape(t)%column(f)
               end do
               if (ioerr /= 0) then
                  write (6,*) 'READ ioerror for column data',ioerr,' on i/o unit = ',nrg
                  call endrun ('READ_RESTART_HISTORY')
               end if
            else
               allocate(tape(t)%column(1))
               read (nrg,iostat=ioerr) tape(t)%column(1)
            end if

            do f=1,nflds(t)
               read (nrg,iostat=ioerr) tape(t)%hlist(f)%field%name,        &
                                       tape(t)%hlist(f)%field%long_name,   &
                                       tape(t)%hlist(f)%field%units,       &
                                       tape(t)%hlist(f)%field%numlev,      &
                                       tape(t)%hlist(f)%field%decomp_type, &
                                       tape(t)%hlist(f)%avgflag,           &
                                       tape(t)%hlist(f)%time_op,           &
                                       tape(t)%hlist(f)%hbuf_prec,         &
                                       tape(t)%hlist(f)%hwrt_prec,         &
                                       tape(t)%hlist(f)%field%flag_xyfill, &
                                       tape(t)%hlist(f)%field%flag_isccplev 

               if (ioerr /= 0) then
                  write(6,*)'End or error condition reading history restart field ',f,' from tape ',t
                  write(6,*)'ioerr=',ioerr
                  call endrun ('READ_RESTART_HISTORY')
               end if
            end do
            if (rgnht(t)) then
               read (nrg,iostat=ioerr) hrestpath(t)
               if (ioerr /= 0) then
                  write (6,*) 'READ ioerror on read of filename ',ioerr,' on i/o unit = ',nrg
                  call endrun ('READ_RESTART_HISTORY')
               end if
            end if
         end do

      end if

#if ( defined SPMD )
      call mpibcast (rgnht  ,ptapes   ,mpilog ,0,mpicom)

      call mpibcast (mtapes ,1        ,mpiint ,0,mpicom)
      call mpibcast (nhtfrq ,ptapes   ,mpiint ,0,mpicom)
      call mpibcast (nflds  ,ptapes   ,mpiint ,0,mpicom)
      call mpibcast (nfils  ,ptapes   ,mpiint ,0,mpicom)
      call mpibcast (mfilt  ,ptapes   ,mpiint ,0,mpicom)

      do t=1,mtapes
         do f=1,nflds(t)
            call mpibcast (tape(t)%hlist(f)%field%numlev,       1,        mpiint, 0,mpicom)
            call mpibcast (tape(t)%hlist(f)%field%decomp_type,  1,        mpiint, 0,mpicom)
            call mpibcast (tape(t)%hlist(f)%field%name,         8,        mpichar,0,mpicom)
            call mpibcast (tape(t)%hlist(f)%field%units,        max_chars,mpichar,0,mpicom)
            call mpibcast (tape(t)%hlist(f)%avgflag,            1,        mpichar,0,mpicom)
            call mpibcast (tape(t)%hlist(f)%hbuf_prec,          1,        mpiint, 0,mpicom)
            call mpibcast (tape(t)%hlist(f)%hwrt_prec,          1,        mpiint, 0,mpicom)
            call mpibcast (tape(t)%hlist(f)%field%flag_xyfill,  1,        mpilog, 0,mpicom)
            call mpibcast (tape(t)%hlist(f)%field%flag_isccplev,1,        mpilog, 0,mpicom)
         end do
      end do
#endif
!
! Allocate space for history buffers and initialize
!
      do t=1,mtapes
         do f=1,nflds(t)
            numlev   = tape(t)%hlist(f)%field%numlev
            select case (tape(t)%hlist(f)%field%decomp_type)
            case (phys_decomp)
               tape(t)%hlist(f)%field%begdim3 = begchunk
               tape(t)%hlist(f)%field%enddim3 = endchunk
               tape(t)%hlist(f)%field%begver = 1
               tape(t)%hlist(f)%field%endver = numlev
               allocate (tape(t)%hlist(f)%field%colperdim3(begchunk:endchunk))
               do c=begchunk,endchunk
                  ncol = get_ncols_p(c)
                  tape(t)%hlist(f)%field%colperdim3(c) = ncol
               end do
               tape(t)%hlist(f)%field%coldimin = pcols
            case (dyn_decomp)
               tape(t)%hlist(f)%field%begdim3 = beglat
               tape(t)%hlist(f)%field%enddim3 = endlat
               if ( dycore_is('LR') )then
# if ( defined STAGGERED )
                  select case (numlev)
                  case (1)
                     tape(t)%hlist(f)%field%begver = 1
                     tape(t)%hlist(f)%field%endver = 1
                  case (plev)
                     tape(t)%hlist(f)%field%begver = beglev
                     tape(t)%hlist(f)%field%endver = endlev
                  case (plevp)
                     tape(t)%hlist(f)%field%begver = beglev
                     tape(t)%hlist(f)%field%endver = endlevp
                  case default
                     write(6,*)'READ_RESTART_HISTORY: invalid number of levels=', numlev
                     call endrun ()
                  end select
# endif
               else
                  tape(t)%hlist(f)%field%begver = 1
                  tape(t)%hlist(f)%field%endver = numlev
               endif
               allocate (tape(t)%hlist(f)%field%colperdim3(beglat:endlat))
               do c=beglat,endlat
                  tape(t)%hlist(f)%field%colperdim3(c) = nlon(c)
               end do
               tape(t)%hlist(f)%field%coldimin = plon
            case default
               write(6,*)'READ_RESTART_HISTORY: bad decomp_type=',tape(t)%hlist(f)%field%decomp_type
               call endrun ()
            end select

            coldimin = tape(t)%hlist(f)%field%coldimin
            begdim3  = tape(t)%hlist(f)%field%begdim3
            enddim3  = tape(t)%hlist(f)%field%enddim3
            begver   = tape(t)%hlist(f)%field%begver
            endver   = tape(t)%hlist(f)%field%endver

            dimind = dim_index_3d (1,coldimin,begver,endver,begdim3,enddim3)
            call allocate_hbuf (tape(t)%hlist(f)%hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
            tape(t)%hlist(f)%hbuf = 0._r8
            allocate (tape(t)%hlist(f)%nacs(coldimin,begdim3:enddim3))
            tape(t)%hlist(f)%nacs(:coldimin,begdim3:enddim3) = 0
         end do
      end do
!
! must allocate space and initialize field column names (this field is dimensioned by number of column groups)
!
      do t=1,mtapes
         do f=1,nflds(t)
            if (ngroup(t) .gt. 0) then
               allocate(tape(t)%hlist(f)%field_column_name(ngroup(t)))
               do i=1,ngroup(t)
                  tape(t)%hlist(f)%field_column_name(i) = trim(tape(t)%hlist(f)%field%name) // "_" // &
                       trim(tape(t)%column(i)%lon_name) // "_" // trim(tape(t)%column(i)%lat_name)
               end do
            else
               allocate(tape(t)%hlist(f)%field_column_name(1))
               tape(t)%hlist(f)%field_column_name(1) = ' '
            end if
         end do
      end do
!
!-----------------------------------------------------------------------
! Read history restart files
!-----------------------------------------------------------------------
!
! Loop over the total number of history files declared and
! read the pathname for any history restart files
! that are present (if any). Test to see if the run is a restart run
! AND if any history buffer regen files exist (rgnht=.T.). Note, rgnht 
! is preset to false, reset to true in routine WSDS if hbuf restart files
! are written and saved in the master restart file. Each history buffer
! restart file is then obtained.
! Note: some f90 compilers (e.g. SGI) complain about I/O of 
! derived types which have pointer components, so explicitly read each one.
! 
      do t=1,mtapes
         if (rgnht(t)) then
!
! Open history restart file
!
            if (masterproc) then
               call getfil (hrestpath(t), locfn)
               call opnfil (locfn, luhrest, 'u')
            end if
!
! Read history restart file
!
            do f=1,nflds(t)
               coldimin   =  tape(t)%hlist(f)%field%coldimin
               begdim3    =  tape(t)%hlist(f)%field%begdim3
               enddim3    =  tape(t)%hlist(f)%field%enddim3
               numlev     =  tape(t)%hlist(f)%field%numlev
               nacs       => tape(t)%hlist(f)%nacs(:coldimin,begdim3:enddim3)
#ifdef SPMD
               if (masterproc) then
                  dimind = dim_index_3d (1,plon,1,numlev,1,plat)
                  call allocate_hbuf (hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
                  allocate (fullnacs(plon,plat))
               else
                  call assoc_hbuf_with_nothing (hbuf,tape(t)%hlist(f)%hbuf_prec)
                  fullnacs => nothing_int
               end if

               select case (tape(t)%hlist(f)%field%decomp_type)
               case (phys_decomp)
                  if (masterproc) then
                     call read_hbuf (hbuf,luhrest,ioerr)
                     read (luhrest) fullnacs
                  end if
                  call scatter_field_to_chunk_hbuf (1, numlev, 1, plon, hbuf, tape(t)%hlist(f)%hbuf)
                  call scatter_field_to_chunk_int (1, 1, 1, plon, fullnacs, tape(t)%hlist(f)%nacs)
               case (dyn_decomp)
                  if ( dycore_is('LR') )then
# if ( defined STAGGERED )
! NEW LR CODING
                     if (tape(t)%hlist(f)%hbuf_prec == 8) then
                        lenarr = plon*numlev*plat
                        select case (numlev)
                        case (1)
                           call lrreadin(luhrest, strip2d, tape(t)%hlist(f)%hbuf%buf8,  &
                                         lenarr, 2)
                        case (plev)
                           call lrreadin(luhrest, strip3dxzy, tape(t)%hlist(f)%hbuf%buf8,  &
                                         lenarr, 3)
                        case (plevp)
                           call lrreadin(luhrest, strip3dxzyp, tape(t)%hlist(f)%hbuf%buf8,  &
                                         lenarr, 3)
                        case default
                           write(6,*)'READ_RESTART_HISTORY: invalid number of levels=', numlev
                           call endrun ()
                        end select
                        lenarr = plon*plat
                        call lrreadini(luhrest, strip2d, nacs,  &
                                      lenarr, 2)
                     else
                        lenarr = plon*numlev*plat
                        select case (numlev)
                        case (1)
                           call lrreadin4(luhrest, strip2d, tape(t)%hlist(f)%hbuf%buf4,  &
                                         lenarr, 2)
                        case (plev)
                           call lrreadin4(luhrest, strip3dxzy, tape(t)%hlist(f)%hbuf%buf4,  &
                                         lenarr, 3)
                        case (plevp)
                           call lrreadin4(luhrest, strip3dxzyp, tape(t)%hlist(f)%hbuf%buf4,  &
                                         lenarr, 3)
                        case default
                           write(6,*)'READ_RESTART_HISTORY: invalid number of levels=', numlev
                           call endrun ()
                        end select
                        lenarr = plon*plat
                        call lrreadini(luhrest, strip2d, nacs,  &
                                      lenarr, 2)
                     endif
# endif
                  else
                     call readin_hbuf(luhrest, tape(t)%hlist(f)%hbuf, coldimin*numlev)
                     call readin_int (luhrest, nacs, coldimin)
                  endif
               case default
                  write(6,*)'READ_RESTART_HISTORY: bad decomp_type=',tape(t)%hlist(f)%field%decomp_type
                  call endrun ()
               end select

               if (masterproc) then
                  call deallocate_hbuf (hbuf)
                  deallocate (fullnacs)
               else
                  call nullify_hbuf (hbuf)
                  nullify (fullnacs)
               end if
#else
               select case (tape(t)%hlist(f)%field%decomp_type)
               case (phys_decomp)
                  dimind = dim_index_3d (1,plon,1,numlev,1,plat)
                  call allocate_hbuf (hbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
                  call read_hbuf (hbuf,luhrest,ioerr)
                  call scatter_field_to_chunk_hbuf (1, numlev, 1, plon, hbuf, tape(t)%hlist(f)%hbuf)
                  call deallocate_hbuf (hbuf)
                  allocate (fullnacs(plon,plat))
                  read (luhrest) fullnacs
                  call scatter_field_to_chunk_int (1, 1, 1, plon, fullnacs, tape(t)%hlist(f)%nacs)
                  deallocate (fullnacs)
               case (dyn_decomp)
                  call readin_hbuf(luhrest, tape(t)%hlist(f)%hbuf, coldimin*numlev)
                  call readin_int (luhrest, nacs, coldimin)
               case default
                  write(6,*)'READ_RESTART_HISTORY: bad decomp_type=',tape(t)%hlist(f)%field%decomp_type
                  call endrun ()
               end select
#endif
            end do
!          
! Done reading this history restart file
!
            if (masterproc) close (luhrest)

         end if  ! rgnht(t)
      end do     ! end of do mtapes loop
!
! If the history files are partially complete (contain less than
! mfilt(t) time samples, then get the files and open them.
!
      do t=1,mtapes
         if (masterproc .and. nfils(t) > 0) then
            call getfil (cpath(t), locfn)
            call wrap_open (locfn, NF_WRITE, nfid(t))
            call h_inquire (t)
         end if
!
! If the history file is full, close the current unit
!
         if (nfils(t) >= mfilt(t)) then
            if (masterproc) then
               write(6,*)'READ_RESTART_HISTORY: nf_close(',nfid(t),')=',nhfil(t)
               call wrap_close (nfid(t))
            end if
            nfils(t) = 0
         end if
      end do
!
! set flag indicating h-tape contents are now defined (needed by addfld)
!
      htapes_defined = .true.      

      return
   end subroutine read_restart_history

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


   character(len=nlen) function get_hfilepath( tape )
!----------------------------------------------------------------------- 
! 
! Purpose: Return full filepath of history file for given tape number
! This allows public read access to the filenames without making
! the filenames public data.
!
!----------------------------------------------------------------------- 
  integer, intent(in) :: tape  ! Tape number

  get_hfilepath = cpath( tape )
  end function get_hfilepath

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


   character(len=nlen) function get_hist_restart_filepath( tape )
!----------------------------------------------------------------------- 
! 
! Purpose: Return full filepath of restart file for given tape number
! This allows public read access to the filenames without making
! the filenames public data.
!
!----------------------------------------------------------------------- 
  integer, intent(in) :: tape  ! Tape number

  get_hist_restart_filepath = hrestpath( tape )
  end function get_hist_restart_filepath

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


  integer function get_mtapes( ) 1
!----------------------------------------------------------------------- 
! 
! Purpose: Return the number of tapes being used.
! This allows public read access to the number of tapes without making
! mtapes public data.
!
!----------------------------------------------------------------------- 
  get_mtapes = mtapes
  end function get_mtapes


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


   subroutine fldlst () 1,22
!----------------------------------------------------------------------- 
! 
! Purpose: Define the contents of each history file based on namelist input for initial or branch
! run, and restart data if a restart run.
!          
! Method: Use arrays fincl and fexcl to modify default history tape contents.
!         Then sort the result alphanumerically for later use by OUTFLD to
!         allow an n log n search time.
! 
!-----------------------------------------------------------------------
#include <comctl.h>
!---------------------------Local variables-----------------------------
!
      integer t, f                   ! tape, field indices
      integer ff                     ! index into include, exclude and fprec list
      character(len=8) :: name       ! field name portion of fincl (i.e. no avgflag separator)
      character(len=8) :: mastername ! name from masterlist field
      character(len=max_chars) :: tmpcolumn_name,latlonname,latlonnamep1 ! tmp char fields
      character(len=1) :: avgflag    ! averaging flag
      character(len=1) :: prec_acc   ! history buffer precision flag
      character(len=1) :: prec_wrt   ! history buffer write precision flag

      type (hentry) :: tmp           ! temporary used for swapping
      type (column_info) :: tmpcolumn ! temporary used for swapping
!
! First ensure contents of fincl, fexcl, fhstpr and fwrtpr are all valid names
!
      do t=1,ptapes
         f = 1
         do while (f < pflds .and. fincl(f,t) /= ' ')
            name = getname (fincl(f,t))
            do ff=1,nfmaster
               mastername = masterlist(ff)%field%name
               if (name == mastername) exit
            end do
            if (name /= mastername) then
               write(6,*)'FLDLST: ', trim(name), ' in fincl(', f, ') not found'
               call endrun
            end if
            f = f + 1
         end do

         f = 1
         do while (f < pflds .and. fexcl(f,t) /= ' ')
            do ff=1,nfmaster
               mastername = masterlist(ff)%field%name
               if (fexcl(f,t) == mastername) exit
            end do
            if (fexcl(f,t) /= mastername) then
               write(6,*)'FLDLST: ', fexcl(f,t), ' in fexcl(', f, ') not found'
               call endrun
            end if
            f = f + 1
         end do

         f = 1
         do while (f < pflds .and. fhstpr(f,t) /= ' ')
            name = getname (fhstpr(f,t))
            do ff=1,nfmaster
               mastername = masterlist(ff)%field%name
               if (name == mastername) exit
            end do
            if (name /= mastername) then
               write(6,*)'FLDLST: ', trim(name), ' in fhstpr(', f, ') not found'
               call endrun
            end if
            do ff=1,f-1                 ! If duplicate entry is found, stop
               if (trim(name) == trim(getname(fhstpr(ff,t)))) then
                  write(6,*)'FLDLST: Duplicate field ', name, ' in fhstpr'
                  call endrun
               end if
            end do
            f = f + 1
         end do

         f = 1
         do while (f < pflds .and. fwrtpr(f,t) /= ' ')
            name = getname (fwrtpr(f,t))
            do ff=1,nfmaster
               mastername = masterlist(ff)%field%name
               if (name == mastername) exit
            end do
            if (name /= mastername) then
               write(6,*)'FLDLST: ', trim(name), ' in fwrtpr(', f, ') not found'
               call endrun
            end if
            do ff=1,f-1                 ! If duplicate entry is found, stop
               if (trim(name) == trim(getname(fwrtpr(ff,t)))) then
                  write(6,*)'FLDLST: Duplicate field ', name, ' in fwrtpr'
                  call endrun
               end if
            end do
            f = f + 1
         end do
      end do
!
! If kind values r8 and r4 are identical, set accumulation precision to 8 bytes
!
      if (r4 == r8 .and. any(nhstpr == 4)) then
         nhstpr(:) = 8
         if (masterproc) then
            write(6,*) 'FLDLST: Set nhstpr to 8 because kind values r8 and r4 are identical'
         end if
      end if

      nflds(:) = 0

      do t=1,ptapes
!
! Add the field to the tape if specified via namelist (FINCL[1-ptapes]), or if
! it is on by default and was not excluded via namelist (FEXCL[1-ptapes]).
! Also set history buffer accumulation and output precision values according
! to the values specified via namelist (FHSTPR[1-ptapes] and FWRTPR[1-ptapes])
! or, if not on the list, to the default values given by ndens(t) and
! nhstpr(t), respectively.
!
         do f=1,nfmaster
            mastername = masterlist(f)%field%name
            call list_index (fhstpr(1,t), mastername, ff)
            if (ff > 0) then
               prec_acc = getflag(fhstpr(ff,t))
            else
               prec_acc = ' '
            end if

            call list_index (fwrtpr(1,t), mastername, ff)
            if (ff > 0) then
               prec_wrt = getflag(fwrtpr(ff,t))
            else
               prec_wrt = ' '
            end if

            call list_index (fincl(1,t), mastername, ff)
            if (ff > 0) then
               avgflag = getflag (fincl(ff,t))
               call inifld (t, f, avgflag, prec_acc, prec_wrt)
            else if (.not. empty_htapes) then
               call list_index (fexcl(1,t), mastername, ff)
               if (ff == 0 .and. masterlist(f)%actflag(t)) then
                  call inifld (t, f, ' ', prec_acc, prec_wrt)
               end if
            end if
         end do
!
! If column output is specified make sure there are some fields defined
! for that tape
!
         if (nflds(t) .eq. 0 .and. fincllonlat(1,t) .ne. ' ') then
            write(6,*) 'FLDLST: Column output is specified for tape ',t,' but no fields defined for that tape.'
            call endrun()
         end if


!
! Specification of tape contents now complete.  Sort each list of active 
! entries for efficiency in OUTFLD.  Simple bubble sort.
!
         do f=nflds(t)-1,1,-1
            do ff=1,f

               if (tape(t)%hlist(ff)%field%name > tape(t)%hlist(ff+1)%field%name) then
            
                  tmp = tape(t)%hlist(ff)
                  tape(t)%hlist(ff  ) = tape(t)%hlist(ff+1)
                  tape(t)%hlist(ff+1) = tmp

               else if (tape(t)%hlist(ff  )%field%name == tape(t)%hlist(ff+1)%field%name) then
                  
                  write(6,*)'FLDLST: Duplicate field ', tape(t)%hlist(ff  )%field%name
                  write(6,*)'t,ff,name=',t,ff,tape(t)%hlist(ff  )%field%name
                  call endrun
            
               end if

            end do
         end do
!
! Bubble sort columns check for duplicates and rebuild field_column_names in newly sorted order
!
         if (ngroup(t) .gt. 0) then
            do i=ngroup(t)-1,1,-1
               do ff=1,i
                  latlonname=trim(tape(t)%column(ff)%lon_name) // "_" // trim(tape(t)%column(ff)%lat_name)
                  latlonnamep1=trim(tape(t)%column(ff+1)%lon_name) // "_" // trim(tape(t)%column(ff+1)%lat_name)
                  if (trim(latlonname) > trim(latlonnamep1)) then
                     tmpcolumn = tape(t)%column(ff)
                     tape(t)%column(ff) = tape(t)%column(ff+1)
                     tape(t)%column(ff+1) = tmpcolumn
                  else if (trim(latlonname) == trim(latlonnamep1)) then
                     write(6,*)'FLDLST: Duplicate column entry for tape ',t,'duplicate column name is ',trim(latlonname)
                     call endrun
                  end if
               end do
            end do
            do f=1,nflds(t)
               do i=1,ngroup(t)
                  if (ngroup(t) .gt. 0) then
                     tape(t)%hlist(f)%field_column_name(i) = trim(tape(t)%hlist(f)%field%name) // "_" // &
                          trim(tape(t)%column(i)%lon_name) // "_" // trim(tape(t)%column(i)%lat_name)
                  else 
                     tape(t)%hlist(f)%field_column_name(1) = ' '
                  end if
               end do
            end do
         end if

         if (masterproc) then
            if (nflds(t) > 0) then
               write(6,*)'FLDLST: Included fields tape ',t,'=',nflds(t)
            end if
   
            do f=1,nflds(t)
               if (ngroup(t) .gt. 0) then
                  if (f.eq.1) write(6,*)'Fields on this tape will be output as column data (FIELD_LON_LAT)'
                  do i=1,ngroup(t)
                     ff=(f-1)*ngroup(t)+i
                     write(6,*) ff,' ',trim(tape(t)%hlist(f)%field_column_name(i)), tape(t)%hlist(f)%field%numlev, &
                       ' ',tape(t)%hlist(f)%avgflag
                  end do
               else
                  write(6,*) f,' ',tape(t)%hlist(f)%field%name, tape(t)%hlist(f)%field%numlev, &
                            ' ',tape(t)%hlist(f)%avgflag
               end if
            end do
         end if
      end do
!
! Determine total number of active history tapes
!
      mtapes = 0
      do t=ptapes,1,-1
         if (nflds(t) > 0) then
            mtapes = t
            exit
         end if
      end do
!
! Ensure there are no "holes" in tape specification, i.e. empty tapes. Enabling
! holes should not be difficult if necessary.
!
      do t=1,mtapes
         if (nflds(t)  ==  0) then
            write(6,*)'FLDLST: Tape ',t,' is empty'
            call endrun
         end if
      end do
      
!
! Packing density, ndens: With netcdf, only 1 (nf_double) and 2 (nf_float)
! are allowed
! Accumulation precision, nhstpr, must be either 8 (real*8) or 4 (real*4)
!
      do t=1,mtapes
         if (ndens(t) == 1) then
            ncprec(t) = nf_double
         else if (ndens(t) == 2) then
            ncprec(t) = nf_float
         else
            call endrun ('FLDLST: ndens must be 1 or 2')
         end if

         if (nhstpr(t) /= 8 .and. nhstpr(t) /= 4) then
            call endrun ('FLDLST: nhstpr must be 8 or 4')
         end if
      end do

      if (masterproc) then
         if (nhtfrq(1) == 0) then
            write(6,*)'History File 1 write frequency MONTHLY'
         else
            write(6,*)'History File 1 write frequency ',nhtfrq(1)
         end if
         
         do t=2,mtapes
            write(6,*)'History File ',t,' write frequency ',nhtfrq(t)
         end do

         do t=1,mtapes
            write(6,*)'Accumulation precision history file ', t, '=', nhstpr(t)
            write(6,*)'Packing density history file ', t, '=', ndens(t)
            write(6,*)'Number of time samples per file (MFILT) for history file ',t,' is ',mfilt(t)
         end do
      end if
!
! set flag indicating h-tape contents are now defined (needed by addfld)
!
      htapes_defined = .true.      
!
!  Now that masterlist is defined, construct primary and secondary hashing
!  tables.
!
      call bld_outfld_hash_tbls()
      call bld_htapefld_indices()

      return
   end subroutine fldlst

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


   subroutine inifld (t, f, avgflag, prec_acc, prec_wrt) 2,9
!----------------------------------------------------------------------- 
! 
! Purpose: Add a field to the active list for a history tape
! 
! Method: Copy the data from the master field list to the active list for the tape
!         Also: define mapping arrays from (col,chunk) -> (lon,lat)
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
     use commap, only: latdeg,londeg 
!
! Arguments
!
      integer, intent(in) :: t            ! history tape index
      integer, intent(in) :: f            ! field index from master field list

      character*1, intent(in) :: avgflag  ! averaging flag
      character*1, intent(in) :: prec_acc ! history buffer precision flag
      character*1, intent(in) :: prec_wrt ! history output precision flag
!
! Local workspace
!
      integer :: n                        ! field index on defined tape
      integer :: ff            ! column index
      integer :: lonind(pflds,2) !   beginning and ending longitude range
      integer :: latind(pflds,2) !   beginning and ending latitude range
      integer :: group !   number of group
      character(len=max_chars) :: lonlatname(pflds), lonname(pflds), latname(pflds) ! variable name
!
! Ensure that it is not to late to add a field to the history tape
!
      if (htapes_defined) then
         call endrun ('INIFLD: Attempt to add field '//masterlist(f)%field%name//' after history files set')
      end if

      nflds(t) = nflds(t) + 1
      n = nflds(t)
!
! Copy field info.
!
      tape(t)%hlist(n)%field = masterlist(f)%field
!
! Set history buffer size and its output data type flags. Set them to
! the default values given by, respective, nhstpr(t) and ndens(t)
! if the input flags prec_acc and prec_wrt are blank; otherwise set to
! the specified values.
!
#ifdef CRAY
      tape(t)%hlist(n)%hbuf_prec = 8
      tape(t)%hlist(n)%hwrt_prec = 8

      if (prec_acc /= ' ' .and. prec_acc /= '8') then
         if (prec_acc == '4') then
            if (masterproc) then
               write(6,*) 'INIFLD: Requested change in history buffer size for ', &
                          tape(t)%hlist(n)%field%name, ' ignored'
            end if
         else
            call endrun ('INIFLD: unknown prec_acc='//prec_acc)
         end if
      end if

      if (prec_wrt /= ' ' .and. prec_wrt /= '8') then
         if (prec_wrt == '4') then
            if (masterproc) then
               write(6,*) 'INIFLD: Requested change in history output size for ', &
                          tape(t)%hlist(n)%field%name, ' ignored'
            end if
         else
            call endrun ('INIFLD: unknown prec_wrt='//prec_wrt)
         end if
      end if
#else
      select case (prec_acc)
      case (' ')
         tape(t)%hlist(n)%hbuf_prec = nhstpr(t)
      case ('4')
         if (r4 /= r8) then
            tape(t)%hlist(n)%hbuf_prec = 4
            if (masterproc) then
               write(6,*) 'INIFLD: History buffer for ', tape(t)%hlist(n)%field%name, &
                          ' is real*4'
            end if
         else       ! if kind values r4 and r8 are identical, ignore the request
            tape(t)%hlist(n)%hbuf_prec = 8
            if (masterproc) then
               write(6,*) 'INIFLD: Requested change in history output size for ', &
                           tape(t)%hlist(n)%field%name, ' ignored'
               write(6,*) '        because kind values r8 and r4 are identical'
            end if
         end if
      case ('8')
         tape(t)%hlist(n)%hbuf_prec = 8
         if (masterproc) then
            write(6,*) 'INIFLD: History buffer for ', tape(t)%hlist(n)%field%name, &
                       ' is real*8'
         end if
      case default
         call endrun ('INIFLD: unknown prec_acc='//prec_acc)
      end select

      select case (prec_wrt)
      case (' ')
         if (ndens(t) == 1) then
            tape(t)%hlist(n)%hwrt_prec = 8
         else
            tape(t)%hlist(n)%hwrt_prec = 4
         end if
      case ('4')
         tape(t)%hlist(n)%hwrt_prec = 4
         if (masterproc) then
            write(6,*) 'INIFLD: Output data type for ', tape(t)%hlist(n)%field%name, &
                       ' is real*4'
         end if
      case ('8')
         tape(t)%hlist(n)%hwrt_prec = 8
         if (masterproc) then
            write(6,*) 'INIFLD: Output data type for ', tape(t)%hlist(n)%field%name, &
                       ' is real*8'
         end if
      case default
         call endrun ('INIFLD: unknown prec_wrt='//prec_wrt)
      end select
#endif
!
! Override the default averaging (masterlist) averaging flag if non-blank
!
      if (avgflag == ' ') then
         tape(t)%hlist(n)%avgflag = masterlist(f)%avgflag(t)
         tape(t)%hlist(n)%time_op = masterlist(f)%time_op(t)
      else
         tape(t)%hlist(n)%avgflag = avgflag
         select case (avgflag)
         case ('A')
            tape(t)%hlist(n)%time_op = 'mean'
         case ('I')
            tape(t)%hlist(n)%time_op = ' '
         case ('X')
            tape(t)%hlist(n)%time_op = 'maximum'
         case ('M')
            tape(t)%hlist(n)%time_op = 'minimum'
         case default
            call endrun ('INIFLD: unknown avgflag='//avgflag)
         end select
      end if
!
! Setup column information if this field will be written as group
! First verify the column information in the namelist
!
      ff = 1
      do while (fincllonlat(ff,t) /= ' ')
         lonlatname(ff) = trim(fincllonlat(ff,t))
         call getlatind(lonlatname(ff),latind(ff,1),latind(ff,2),latname(ff))
         call getlonind(lonlatname(ff),lonind(ff,1),lonind(ff,2),lonname(ff))
#ifdef HDEBUG
         write(6,*)'closest lon lat range for group ',lonlatname(ff),' is ',lonind(ff,:),latind(ff,:)
#endif
         ff = ff + 1
      end do

      group = ff-1

      ngroup(t)=group

      if (group .gt. 0) then 

         allocate (tape(t)%column(group))
         allocate (tape(t)%hlist(n)%field_column_name(group))

         do ff = 1, group
            tape(t)%column(ff)%lat_name=latname(ff)
            tape(t)%column(ff)%lon_name=lonname(ff)
            tape(t)%column(ff)%columnlon(:)=lonind(ff,:)
            tape(t)%column(ff)%columnlat(:)=latind(ff,:)
            tape(t)%column(ff)%num_lats = &
                 tape(t)%column(ff)%columnlat(2)-tape(t)%column(ff)%columnlat(1)+1
            tape(t)%column(ff)%num_lons = &
                 tape(t)%column(ff)%columnlon(2)-tape(t)%column(ff)%columnlon(1)+1
            tape(t)%hlist(n)%field_column_name(ff) = &
                 trim(tape(t)%hlist(n)%field%name) // "_" // trim(lonname(ff)) // "_" // trim(latname(ff))
        end do


      else

         allocate (tape(t)%column(1))
         allocate (tape(t)%hlist(n)%field_column_name(1))

         tape(t)%hlist(n)%field_column_name(1) = ' '
         tape(t)%column(1)%lat_name=' '
         tape(t)%column(1)%lon_name=' '
         tape(t)%column(1)%columnlon(:)=0
         tape(t)%column(1)%columnlat(:)=0
         tape(t)%column(1)%num_lats=0
         tape(t)%column(1)%num_lons=0

      end if

#ifdef HDEBUG
      write(6,*)'INIFLD: field ', tape(t)%hlist(n)%field%name, ' added as ', 'field number ', n,' on tape ', t
      write(6,*)'units=',tape(t)%hlist(n)%field%units
      write(6,*)'numlev=',tape(t)%hlist(n)%field%numlev
      write(6,*)'avgflag=',tape(t)%hlist(n)%avgflag
      write(6,*)'time_op=',tape(t)%hlist(n)%time_op
      write(6,*)'hbuf_prec=',tape(t)%hlist(n)%hbuf_prec
      write(6,*)'hwrt_prec=',tape(t)%hlist(n)%hwrt_prec
#endif

      return
   end subroutine inifld

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


   character(len=8) function getname (inname),1
!----------------------------------------------------------------------- 
! 
! Purpose: retrieve name portion of inname
!          
! Method:  If an averaging flag separater character is present (":") in inname, 
!          lop it off
! 
!-------------------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: inname
!
! Local workspace
!
      integer :: length
      integer :: i
   
      length = len (inname)
      
      if (length < 8 .or. length > 10) then
         write(6,*) 'GETNAME: bad length=',length
         call endrun
      end if
   
      getname = ' '
      do i=1,8
         if (inname(i:i) == ':') exit
         getname(i:i) = inname(i:i)
      end do
      
      return
   end function getname

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


   subroutine getlatind(inname,beglatind,endlatind,latname) 1,5
!----------------------------------------------------------------------- 
! 
! Purpose: retrieve closes model lat index given north south latitude string
!          
! Author: John Truesdale
! 
!-------------------------------------------------------------------------------
     use commap, only: latdeg
!
! Arguments
!
      character(len=*) inname,latname
      integer :: beglatind,endlatind
!
! Local workspace
!
      character(len=max_chars) str1,str2,tmpstr(4)
      integer :: i,j,marker,latind(2),tmplen(4)
      real(r8) :: latdegree(2)

!
!-------------------------------------------------------------------------------
!
      str1=' '
      str2=' '
      tmpstr(:)= ' '

!
! make sure _ separator is present
!      
      if (scan(inname,'_').eq.0) then
         write(6,*)'GETLATIND: Improperly formatted column string.  Missing underscore character (xxxE_yyyS) ', &
                   inname,scan(inname,'_')
         call endrun
      end if
!
! split inname string into lat lon strings
!      
      marker=scan(inname,'_')
      str1=inname(:marker-1)
      str2=inname(marker+1:)

!
! split ranges of lats( or lons) into seperate substrings. Substrings 1,2 will contain lats(lons) from portion of string before underscore
! if a single column and not a range is specified substrings 1 and 2 will be the same
! and substrings 3,4 will contain lats(lons) from portion of main string after underscore character
! if a single column and not a range is specified substrings 3 and 4 will be the same
!
      if (scan(str1,':') .ne. 0) then
         marker=scan(str1,':')
         tmpstr(1)=str1(:marker-1)
         tmpstr(2)=str1(marker+1:)
      else
         tmpstr(1)=str1
         tmpstr(2)=str1
      end if
      if (scan(str2,':') .ne. 0) then
         marker=scan(str2,':')
         tmpstr(3)=str2(:marker-1)
         tmpstr(4)=str2(marker+1:)
      else
         tmpstr(3)=str2
         tmpstr(4)=str2
      end if

! check format of substrings - (number followed by single character north/south east/west designation)

      do i = 1,4
         tmplen(i)=len_trim(tmpstr(i))
         if (verify(tmpstr(i),"0123456789.").ne.tmplen(i) .or. verify(tmpstr(i)(tmplen(i):tmplen(i)),"eEwWnNsS") .ne. 0) then
            write(6,*)'GETLATIND (2): Improperly formatted column string. ',&
                 inname,verify(tmpstr(i),"0123456789."),tmplen(i),&
                 verify(tmpstr(i)(tmplen(i):tmplen(i)),"eEwWnNsS")
            call endrun
         end if
      end do

! find latitude substrings and put into temporary work space tmpstr(1:2)

      if (verify(tmpstr(1)(tmplen(1):tmplen(1)),"nNsSs").eq.0  .and. verify(tmpstr(2)(tmplen(2):tmplen(2)),"nNsSs").eq.0 ) then

      else if (verify(tmpstr(3)(tmplen(3):tmplen(3)),"nNsSs").eq.0 .and. verify(tmpstr(3)(tmplen(3):tmplen(3)),"nNsSs").eq.0) then
         tmpstr(1) = tmpstr(3)
         tmplen(1)=tmplen(3)
         tmpstr(2) = tmpstr(4)
         tmplen(2)=tmplen(4)
      else
         call endrun ('GETLATIND (3): Improperly
      end if
!
! convert lat substrings to real and if in southern hemisphere make negative
!
      do i = 1,2
         read(tmpstr(i)(1:tmplen(i)-1),*) latdegree(i)
         if (verify(tmpstr(i)(tmplen(i):tmplen(i)),'sS').eq.0) then
            latdegree(i)=(-1.)*latdegree(i)
         end if
!
! Make sure specified latitudes is in bounds
!
         if (latdegree(i) .lt. -90. .or. latdegree(i) .gt. 90.) then
            write(6,*)'GETLATIND: latitude for column namelist is out of range (-90 .. 90) value=',latdegree(i)
            call endrun
         endif
!
! Find closest lat index for each substring
!
         if (latdegree(i).ge.latdeg(plat)) then
            latind(i)=plat
         else
            do j = 1, plat-1
               if ( abs(latdeg(j)-latdegree(i)) .lt.  &
                    abs(latdeg(j+1)-latdegree(i))) then
                  latind(i)=j
                  exit
               endif
            enddo
         end if
      end do
!
! output begining and ending latitude indicies.   If just a column is specified then beginning latitude index will be the same as the
! ending latitude index  

!      
      if (latind(1) .le. latind(2) ) then
         beglatind =latind(1)
         endlatind =latind(2)
      else 
         beglatind = latind(2)
         endlatind = latind(1)
      end if

      if (beglatind .eq. endlatind) then
         latname = trim(tmpstr(1))
      else
         if (latind(1) .le. latind(2) ) then
            latname = trim(tmpstr(1)) // "_to_" // trim(tmpstr(2))
         else 
            latname = trim(tmpstr(2)) // "_to_" // trim(tmpstr(1))
         end if
      end if
         return
   end subroutine getlatind
!#######################################################################


   subroutine getlonind (inname,beglonind,endlonind,lonname) 1,5
!----------------------------------------------------------------------- 
! 
! Purpose: retrieve closes model lat index given north south latitude string
!          
! Method:  If an averaging flag separater character is present (":") in inname, 
!          lop it off
! 
!-------------------------------------------------------------------------------
     use commap, only: londeg
!
! Arguments
!
      character(len=*) inname,lonname
      integer :: beglonind,endlonind
!
! Local workspace
!
      character(len=max_chars) str1,str2,tmpstr(4)
      integer :: i,j,marker,lonind(2),tmplen(4)
      real(r8) :: londegree(2)

!
!-------------------------------------------------------------------------------
!
      str1=' '
      str2=' '
      tmpstr(:)= ' '
!
! make sure _ separator is present
!      
      if (scan(inname,'_').eq.0) then
         write(6,*)'GETLONIND: Improperly formatted column string.  Missing underscore character (xxxE_yyyS) ',&
              inname,scan(inname,'_')
         call endrun
      end if
!
! split string in to lat lon strings
!      
      marker=scan(inname,'_')
      str1=inname(:marker-1)
      str2=inname(marker+1:)

!
! split ranges of lats( or lons) into seperate substrings. Substrings 1,2 will contain lats(lons) from portion of string before underscore
! if a single column and not a range is specified substrings 1 and 2 will be the same
! and substrings 3,4 will contain lats(lons) from portion of main string after underscore character
! if a single column and not a range is specified substrings 3 and 4 will be the same
!
      if (scan(str1,':') .ne. 0) then
         marker=scan(str1,':')
         tmpstr(1)=str1(:marker-1)
         tmpstr(2)=str1(marker+1:)
      else
         tmpstr(1)=str1
         tmpstr(2)=str1
      end if

      if (scan(str2,':') .ne. 0) then
         marker=scan(str2,':')
         tmpstr(3)=str2(:marker-1)
         tmpstr(4)=str2(marker+1:)
      else
         tmpstr(3)=str2
         tmpstr(4)=str2
      end if

! check format of substrings - (number followed by single character north/south east/west designation)

      do i = 1,4
         tmplen(i)=len_trim(tmpstr(i))
         if (verify(tmpstr(i),"0123456789.").ne.tmplen(i) .or. verify(tmpstr(i)(tmplen(i):tmplen(i)),"eEwWnNsS") .ne. 0) then
            write(6,*)'GETLONIND (2): Improperly formatted column string. ', &
                 inname,verify(tmpstr(i),"0123456789."),tmplen(i), &
                 verify(tmpstr(i)(tmplen(i):tmplen(i)),"eEwWnNsS")
            call endrun
         end if
      end do

! find longitude substrings and put into temporary work space tmpstr(1:2)

      if (verify(tmpstr(1)(tmplen(1):tmplen(1)),"eEwW").eq.0  .and. verify(tmpstr(2)(tmplen(2):tmplen(2)),"eEwW").eq.0 ) then

      else if (verify(tmpstr(3)(tmplen(3):tmplen(3)),"eEwW").eq.0 .and. verify(tmpstr(3)(tmplen(3):tmplen(3)),"eEwW").eq.0) then
         tmpstr(1) = tmpstr(3)
         tmplen(1)=tmplen(3)
         tmpstr(2) = tmpstr(4)
         tmplen(2)=tmplen(4)
      else
         call endrun ('GETLONIND (3): Improperly
      end if
!
! convert lon substrings to real and make sure its degrees east
!
      do i = 1,2
         read(tmpstr(i)(1:tmplen(i)-1),*) londegree(i)
         if (verify(tmpstr(i)(tmplen(i):tmplen(i)),'wW').eq.0) then
            londegree(i) = 360. - londegree(i)
         end if
!
! Make sure specified longitudes are in bounds
!
         if (londegree(i) .lt. 0 .or. londegree(i) .gt. 360.) then
            write(6,*)'GETLONIND: longitude for column namelist is out of range (0 .. 360) value=',londegree(i)
            call endrun
         endif

!
! Find closest lon index for each substring.  If just a column is specified then beginning longitude index will be the same as the
! ending longitude index  
!
         if (londegree(i).ge.londeg(plon,1)) then
            lonind(i)=plon
         else
            do j = 1, plon-1
               if ( abs(londeg(j,1)-londegree(i)) .le. abs(londeg(j+1,1)-londegree(i))) then
                  lonind(i)=j
                  exit
               endif
            end do
         end if
      end do
!
! output begining and ending longitude indicies.   If just a column is specified then beginning longitude index will be the same as the
! ending longitude index  

!      
      if (lonind(1) .le. lonind(2) ) then
         beglonind =lonind(1)
         endlonind =lonind(2)
      else 
         beglonind = lonind(2)
         endlonind = lonind(1)
      end if


      if (beglonind .eq. endlonind) then
         lonname = trim(tmpstr(1))
      else
         if (lonind(1) .le. lonind(2) ) then
            lonname = trim(tmpstr(1)) // "_to_" // trim(tmpstr(2))
         else 
            lonname = trim(tmpstr(2)) // "_to_" // trim(tmpstr(1))
         end if
      end if


      return
    end subroutine getlonind

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


   character(len=1) function getflag (inname) 2,1
!----------------------------------------------------------------------- 
! 
! Purpose: retrieve flag portion of inname
!          
! Method:  If an averaging flag separater character is present (":") in inname, 
!          return the character after it as the flag
! 
!-------------------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: inname   ! character string
!
! Local workspace
!
      integer :: length         ! length of inname
      integer :: i              ! loop index

      length = len (inname)

      if (length /= 10) then
         write(6,*) 'GETFLAG: bad length=',length
         call endrun
      end if

      getflag = ' '
      do i=1,9
         if (inname(i:i) == ':') then
            getflag = inname(i+1:i+1)
            exit
         end if
      end do

      return
   end function getflag

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


   subroutine list_index (list, name, index) 4
!
! Input arguments
!
      character(len=*), intent(in) :: list(pflds) ! input list of names, possibly ":" delimited
      character(len=8), intent(in) :: name        ! name to be searched for
!
! Output arguments
!
      integer, intent(out) :: index               ! index of "name" in "list"
!
! Local workspace
!
      character(len=8) :: listname    ! input name with ":" stripped off.
      integer f                       ! field index

      index = 0
      do f=1,pflds
!
! Only list items
!
         listname = getname (list(f))
         if (listname == ' ') exit
         if (listname == name) then
            index = f
            exit
         end if
      end do
      
      return
   end subroutine list_index

!#######################################################################
#if (defined SCAM)


   subroutine scm_intht (),21
   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid
   use buffer
   use prognostics
   use comsrf
!------------------------------Includes---------------------------------
#include <comadj.h>
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
#include <comfrc.h>
!-----------------------------------------------------------------------
!#include <comtrcnm.h>
!-----------------------------------------------------------------------
#include <comsol.h>
   integer lat
!-----------------------------------------------------------------------
!
! Initialize history file handler
!
! $Id: history.F90,v 1.26.2.48.4.1 2004/05/20 18:36:01 mvr Exp $
! $Author: mvr $
!
   call scm_histfield_ini()
!
! Initialize the outfld data structures
!      
   call init_c_outfld
!      
! Build the Master Field List
!
   call bldfld
!      
! Add fields from Master Field to the scm data structures
!
   call scm_addfield()

!     Call outfld() for all of the modifiable fields so that
!     the outfield buffer field averages will include initial values
!     and also to initialize the pointers to the variables in
!     the outfld buffer.      

   call outfld('PHIS',  phis, plond, lat)
   call outfld('PS',    ps(1,1,n3),   plond, lat)
   call outfld('Q',     q3(1,1,1,1,n3),   plond, lat)
   call outfld('T',     t3(1,1,1,n3),   plond, lat)
   call outfld('U',     u3(1,1,1,n3),   plond, lat)
   call outfld('V',     v3(1,1,1,n3),   plond, lat)
   call outfld('OMEGA', wfld, plond, lat)
   call outfld('DIVQ',  divq, plond, lat)
   call outfld('DIVT',  divt, plond, lat)
   call outfld('DIVQ3D',divq3d, plond, lat)
   call outfld('DIVT3D',divt3d, plond, lat)
   call outfld('DIVU',  divu, plond, lat)
   call outfld('DIVV',  divv, plond, lat)

   return
end subroutine scm_intht

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


   subroutine scm_histfield_ini() 1,16
!----------------------------------------------------------------------- 
! 
! Purpose: 
!
! add master list fields to scm
! 
! Method: Call a subroutine to add each field
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
      use ppgrid, only: pver, pverp
      use comsrf, only: plevmx, tsnam
      use constituents
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
! Local variables
!
      integer m,j        ! Indices
      real(r8) dummy
!
! Call addfld to add each field to the Master Field List.
!
      call addfld ('TDIFF   ','K      ',plev,    'A','difference from observed temp', phys_decomp)

      call addfld ('TOBS    ','K      ',plev,    'A','observed temp', phys_decomp)
      call addfld ('QDIFF   ','kg/kg   ',plev,    'A','difference from observed water',phys_decomp)

      call addfld ('QOBS    ','kg/kg   ',plev,    'A','observed water',phys_decomp)
      call addfld ('PRECOBS','mm/day',plev,    'A','Total (convective and large-scale) precipitation rate', phys_decomp)
      call addfld ('DIVQ    ','kg/kg/s ',plev,    'A','Q advection tendency (horizontal)', phys_decomp)
      call addfld ('DIVQ3D  ','kg/kg/s ',pver,    'A','Q advection tendency (horiz/vert combined)', dyn_decomp)
      call addfld ('DIVV  ','m/s2    ',plev,    'A','V advection tendency (horizontal)', phys_decomp)
      call addfld ('DIVU  ','m/s2    ',plev,    'A','U advection tendency (horizontal)', phys_decomp)
      call addfld ('DIVT   ','K/s     ',plev,    'A','T advection tendency (horizontal)', phys_decomp)
      call addfld ('DIVT3D ','K/s     ',pver,    'A','T advection tendency (horiz/vert combined)', dyn_decomp)

      call addfld ('SHFLXOBS','W/m2    ',1,    'A','Obs Surface sensible heat flux',phys_decomp)
      call addfld ('LHFLXOBS','W/m2    ',1,    'A','Obs Surface latent heat flux',phys_decomp)
!      do m=1,pcnst
!         call addfld (tottnam(m), 'kg/kg/s ',pver, 'A',trim(cnst_name(m))//' Tot tendency due to moist processes',phys_decomp)
!         call addfld (tendnam(m), 'kg/kg/s ',pver, 'A',trim(cnst_name(m))//' tendency due to moist processes',phys_decomp)
!      end do
   end subroutine scm_histfield_ini


   subroutine scm_addfield() 1,2
!----------------------------------------------------------------------- 
! 
! Purpose: 
!
! add master list fields to scm
! 
! Method: Call a subroutine to add each field
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
      use ppgrid, only: pver, pverp
      use comsrf, only: plevmx, tsnam
!-----------------------------------------------------------------------
!
#define    SHOW        1  
#define    DONTSHOW    0
#define    MODIFIABLE  1
#define    DIAGNOSTIC  0
#define    AVERAGE     1
#define    INSTANTANEOUS 0
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
! Local variables
!
      integer m,j,idx1        ! Indices
      real(r8) dummy
      character*256 modifiable_fields
! for modifiable_fields the single letter fields must come first, 
! the two letter fields, then three letter etc.
      do j = 1,nfmaster
         modifiable_fields='!Q!T!U!V!PS!DIVQ!DIVQ3D!DIVT!DIVT3D!DIVU!DIVV!OMEGA!CWAT!PHIS!'
         idx1=index(trim(modifiable_fields),trim(masterlist(j)%field%name))
	 if (idx1.gt.1) then
            if( modifiable_fields(idx1-1:idx1-1) == '!' .and. &
                 modifiable_fields(idx1+len(trim(masterlist(j)%field%name)):idx1+len(trim(masterlist(j)%field%name))) == '!') then
               call addfield (&
                    trim(masterlist(j)%field%name) , &
                    trim(masterlist(j)%field%long_name), &
                    trim(masterlist(j)%field%units), &
                    trim(masterlist(j)%field%units), &
                    1.0_r8, 0._r8, 0._r8, &
                    SHOW, MODIFIABLE, INSTANTANEOUS, &
                    masterlist(j)%field%numlev, &
                    dummy)
            else
               call addfield (&
                    trim(masterlist(j)%field%name) , &
                    trim(masterlist(j)%field%long_name), &
                    trim(masterlist(j)%field%units), &
                    trim(masterlist(j)%field%units), &
                    1.0_r8, 0._r8, 0._r8, &
                    SHOW, DIAGNOSTIC, INSTANTANEOUS, &
                    masterlist(j)%field%numlev, &
                    dummy)
            end if
         else
            call addfield (&
                 trim(masterlist(j)%field%name) , &
                 trim(masterlist(j)%field%long_name), &
                 trim(masterlist(j)%field%units), &
                 trim(masterlist(j)%field%units), &
                 1.0_r8, 0._r8, 0._r8, &
                 SHOW, DIAGNOSTIC, INSTANTANEOUS, &
                 masterlist(j)%field%numlev, &
                 dummy)
         end if
      end do
        
   end subroutine scm_addfield


   subroutine outfld (fname, field, idim, c) 317,7
!-----------------------------------------------------------------------
!
! Purpose: Accumulate (or take min, max, etc. as appropriate) input field
!          into its history buffer for appropriate tapes
!
! Method: Search for fname among fields on history tapes.  If found, do the
!         accumulation.  If not found, return silently.
!
! Author: CCM Core Group
!
!-----------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: fname ! Field name--should be 8 chars long

      integer, intent(in) :: idim           ! Longitude dimension of field array
      integer, intent(in) :: c              ! chunk (physics) or latitude (dynamics) index

      real(r8), intent(in) :: field(idim,*) ! Array containing field values
!
!-----------------------------------------------------------------------

      call c_outfld(fname, field, idim, c)
      return
   end subroutine outfld

!#######################################################################
#else


   subroutine outfld (fname, field, idim, c) 317,7
!----------------------------------------------------------------------- 
! 
! Purpose: Accumulate (or take min, max, etc. as appropriate) input field
!          into its history buffer for appropriate tapes
! 
! Method: Check 'masterlist' whether the requested field 'fname' is active
!         on one or more history tapes, and if so do the accumulation.
!         If not found, return silently.
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: fname ! Field name--should be 8 chars long

      integer, intent(in) :: idim           ! Longitude dimension of field array
      integer, intent(in) :: c              ! chunk (physics) or latitude (dynamics) index

      real(r8), intent(in) :: field(idim,*) ! Array containing field values
!
! Local variables
!
      integer :: t, f                ! tape, field indices
      integer coldimin               ! column dimension of model array
      integer :: fl, fu              ! upper, lower indices used in binary search thru sorted list
      integer :: begver              ! on-node vert start index
      integer :: endver              ! on-node vert end index
      integer :: endi                ! ending longitude index (reduced grid)

      character*8 :: fname8          ! 8-char equivalent of fname
      character*1 :: avgflag         ! averaging flag
      
      type (hbuffer_2d) :: hbuf      ! history buffer
      integer, pointer :: nacs(:)    ! accumulation counter
      type (dim_index_2d) :: dimind  ! 2-D dimension index
      logical :: flag_xyfill         ! non-applicable xy points flagged with fillvalue
      integer :: ff                  ! masterlist index pointer
!-----------------------------------------------------------------------
!      call t_startf ('outfld')
      fname8 = fname
      ff = get_masterlist_indx(fname8)
!
!  If ( ff < 0 ), the field is not defined on the masterlist. This check
!  is necessary because of coding errors calling outfld without first defining
!  the field on masterlist.
!
      if ( ff < 0 ) return
!
!  Next, check to see whether this field is active on one or more history
!  tapes.
!

      if ( .not. masterlist(ff)%act_sometape ) return

!
! Note, the field may be on any or all of the history files (primary
! and auxiliary).
!
!      write(6,*)'fname8=',fname8
      do 40 t=1,ptapes
         if ( .not. masterlist(ff)%actflag(t)) cycle
         f = masterlist(ff)%htapeindx(t)
!
! Update history buffer
!
         begver  = tape(t)%hlist(f)%field%begver
         endver  = tape(t)%hlist(f)%field%endver
         endi    = tape(t)%hlist(f)%field%colperdim3(c)
         avgflag = tape(t)%hlist(f)%avgflag
         flag_xyfill = tape(t)%hlist(f)%field%flag_xyfill
         coldimin= tape(t)%hlist(f)%field%coldimin
         nacs   => tape(t)%hlist(f)%nacs(:coldimin,c)
         call assoc_hbuf2d_with_hbuf3d (hbuf, tape(t)%hlist(f)%hbuf, c)
         
         dimind = dim_index_2d (1,endi,begver,endver)

         select case (avgflag)

         case ('I') ! Instantaneous

            call hbuf_accum_inst (hbuf%buf4, hbuf%buf8, field, nacs, dimind, idim, flag_xyfill)

         case ('A') ! Time average

            call hbuf_accum_add (hbuf%buf4, hbuf%buf8, field, nacs, dimind, idim, flag_xyfill)

         case ('X') ! Maximum over time

            call hbuf_accum_max (hbuf%buf4, hbuf%buf8, field, nacs, dimind, idim, flag_xyfill)

         case ('M') ! Minimum over time

            call hbuf_accum_min (hbuf%buf4, hbuf%buf8, field, nacs, dimind, idim, flag_xyfill)

         case default

            call endrun ('OUTFLD: invalid avgflag='//avgflag)

         end select
40    continue
!      call t_stopf ('outfld')
      return
   end subroutine outfld
#endif

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


   integer function strcmpf (name1, name2)
!----------------------------------------------------------------------- 
! 
! Purpose: Return the lexical difference between two strings
! 
! Method: Use ichar() intrinsic as we loop through the names
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      character(len=8), intent(in) :: name1, name2  ! strings to compare
      integer n                                     ! loop index
   
      do n=1,8
         strcmpf = ichar(name1(n:n)) - ichar(name2(n:n))
         if (strcmpf /= 0) exit
      end do

      return
   end function strcmpf

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


   subroutine h_inquire (t) 1,19
!----------------------------------------------------------------------- 
! 
! Purpose: Ensure that the proper variables are on a history file
! 
! Method: Issue the appropriate netcdf wrapper calls
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      integer, intent(in) :: t   ! tape index
!
! Local workspace
!
      integer f,ff               ! field index
      integer ret                ! return value from function call
      integer londim             ! longitude dimension id
      integer latdim             ! latitude dimension id
      integer levdim             ! level dimension id
      integer ilevdim            ! intfc dimension id
      integer tbnddim            ! time_bnds dimension id
      integer old_mode           ! returned from nf_set_fill
      integer marker             ! index for string marker
!
! Setup netcdf file - create the dimensions of lat,lon,time,level
!     
      ret = nf_set_fill (nfid(t), nf_nofill, old_mode)
!
! Dimension id's
!
      call wrap_inq_dimid (nfid(t), 'lat', latdim)
      call wrap_inq_dimid (nfid(t), 'lon', londim)
      call wrap_inq_dimid (nfid(t), 'lev', levdim)
      call wrap_inq_dimid (nfid(t), 'ilev', ilevdim)
      call wrap_inq_dimid (nfid(t), 'tbnd', tbnddim)
!
! Create variables for model timing and header information 
!
      call wrap_inq_varid (nfid(t),'ndcur   ',    ndcurid(t))
      call wrap_inq_varid (nfid(t),'nscur   ',    nscurid(t))
      call wrap_inq_varid (nfid(t),'date    ',    dateid(t))
      call wrap_inq_varid (nfid(t),'co2vmr  ',    co2vmrid(t))
      call wrap_inq_varid (nfid(t),'datesec ',    datesecid(t))
      call wrap_inq_varid (nfid(t),'nsteph  ',    nstephid(t))
      call wrap_inq_varid (nfid(t),'time    ',    timeid(t))
      call wrap_inq_varid (nfid(t),'time_bnds',   tbndid(t))
      call wrap_inq_varid (nfid(t),'date_written',date_writtenid(t))
      call wrap_inq_varid (nfid(t),'time_written',time_writtenid(t))
#if ( defined BFB_CAM_SCAM_IOP )
      call wrap_inq_varid (nfid(t),'tsec    ',tsecid(t))
      call wrap_inq_varid (nfid(t),'bdate   ',bdateid(t))
#endif
!
! Obtain variable name from ID which was read from restart file
!
      do f=1,nflds(t)
!
! If this field will be put out as columns then get column names for field
!
         if (ngroup(t) .gt. 0) then
            do i=1,ngroup(t)
               ff=(f-1)*ngroup(t)+i
               call wrap_inq_varname (nfid(t), varid(ff,t), tape(t)%hlist(f)%field_column_name(i))
               marker=scan(tape(t)%hlist(f)%field_column_name(i),'_')
               tape(t)%hlist(f)%field%name=tape(t)%hlist(f)%field_column_name(i)(:marker-1)
            end do
         else	
            call wrap_inq_varname (nfid(t), varid(f,t), tape(t)%hlist(f)%field%name)
         end if
      end do
!
      ret = nf_enddef(nfid(t))
      write(6,*)'H_INQUIRE: Successfully opened netcdf file '

      return
   end subroutine h_inquire

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


   subroutine h_default () 1,119
!----------------------------------------------------------------------- 
! 
! Purpose: Define default contents of history files
! 
! Method: Call add_default for each field.  Arguments are field name, tape index, and
!         modification to averaging flag (blank means no modification)
! 
!-----------------------------------------------------------------------
      use dycore, only: dycore_is

#include <comctl.h>

!
! Local workspace
!
      integer m      ! tracer index
!
! First tape (monthly output by default)
!
      call add_default ('PHIS    ', 1, ' ')
      call add_default ('PS      ', 1, ' ')
      call add_default ('PSDRY   ', 1, ' ')
      call add_default ('PDELDRY ', 1, ' ')
      call add_default ('T       ', 1, ' ')
      call add_default ('U       ', 1, ' ')
      call add_default ('V       ', 1, ' ')

#if ( defined STAGGERED )
      call add_default ('US      ', 1, ' ')
      call add_default ('VS      ', 1, ' ')
#endif

#if ( defined COUP_CSM )
      call add_default ('CPLRAINC', 1, ' ')
      call add_default ('CPLRAINL', 1, ' ')
      call add_default ('CPLSNOWC', 1, ' ')
      call add_default ('CPLSNOWL', 1, ' ')
      call add_default ('CPLPRCER', 1, ' ')
#endif

      do m=1,pcnst+pnats
         call add_default (cnst_name(m), 1, ' ')
         call add_default (sflxnam(m), 1, ' ')
      end do

      call add_default (dcconnam(1), 1, ' ')

      if ( .not. dycore_is('LR') )then
         call add_default ('DTH     ', 1, ' ')
      end if
      call add_default ('SNOWHICE', 1, ' ')
      call add_default ('SNOWHLND', 1, ' ')
      call add_default ('PRECL   ', 1, ' ')
      call add_default ('PRECC   ', 1, ' ')
      call add_default ('PRECSH  ', 1, ' ')
      call add_default ('PRECSL  ', 1, ' ')
      call add_default ('PRECSC  ', 1, ' ')
      call add_default ('SHFLX   ', 1, ' ')
      call add_default ('LHFLX   ', 1, ' ')
      call add_default ('QFLX    ', 1, ' ')
      call add_default ('PBLH    ', 1, ' ')
      call add_default ('TREFHT  ', 1, ' ')
#if ( defined COUP_CSM )
      call add_default ('QREFHT  ', 1, ' ')
#endif

      call add_default ('DTV     ', 1, ' ')
      call add_default ('FSNS    ', 1, ' ')
      call add_default ('FLNS    ', 1, ' ')
      call add_default ('FLNT    ', 1, ' ')
      call add_default ('FLUT    ', 1, ' ')
      call add_default ('FLUTC   ', 1, ' ')
      call add_default ('FSNT    ', 1, ' ')
      call add_default ('FSNTOA  ', 1, ' ')
      call add_default ('FSNTOAC ', 1, ' ')
      call add_default ('LWCF    ', 1, ' ')
      call add_default ('SWCF    ', 1, ' ')
      call add_default ('CLOUD   ', 1, ' ')
      call add_default ('CLDTOT  ', 1, ' ')
      call add_default ('CLDLOW  ', 1, ' ')
      call add_default ('CLDMED  ', 1, ' ')
      call add_default ('CLDHGH  ', 1, ' ')
      call add_default ('FLNTC   ', 1, ' ')
!      call add_default ('FLN200  ', 1, ' ')
!      call add_default ('FLN200C ', 1, ' ')
!      call add_default ('FSN200  ', 1, ' ')
!      call add_default ('FSN200C ', 1, ' ')
      call add_default ('FSNTC   ', 1, ' ')
      call add_default ('FLNSC   ', 1, ' ')
      call add_default ('FSNSC   ', 1, ' ')
      call add_default ('FSDSC   ', 1, ' ')
      call add_default ('OMEGA   ', 1, ' ')
      call add_default ('TAUX    ', 1, ' ')
      call add_default ('TAUY    ', 1, ' ')
      call add_default ('SRFRAD  ', 1, ' ')
      call add_default ('QRS     ', 1, ' ')
      call add_default ('QRL     ', 1, ' ')
      call add_default ('TS      ', 1, ' ')
      call add_default ('TSMN    ', 1, ' ')
      call add_default ('TSMX    ', 1, ' ')
      call add_default ('SOLIN   ', 1, ' ')
      call add_default ('DTCOND  ', 1, ' ')
      call add_default ('CMFDT   ', 1, ' ')
      call add_default ('CMFDQ   ', 1, ' ')
      call add_default ('CMFDQR  ', 1, ' ')
      call add_default ('QC      ', 1, ' ')
      call add_default ('CMFMC   ', 1, ' ')
      call add_default ('FSDS    ', 1, ' ')
      call add_default ('CONCLD  ', 1, ' ')
      call add_default ('RELHUM  ', 1, ' ')

!++scyc  add fields associated with sulfur cycle

      if ( indirect ) then
         call add_default ('MSO4    ', 1, ' ')
         call add_default ('LWC     ', 1, ' ')
         call add_default ('CLDFRQ  ', 1, ' ')
         call add_default ('WREL    ', 1, ' ')
         call add_default ('WLWC    ', 1, ' ')
      end if

      call add_default ('VT      ', 1, ' ')
      call add_default ('VU      ', 1, ' ')
      call add_default ('OMEGAT  ', 1, ' ')
      call add_default ('VQ      ', 1, ' ')
      call add_default ('VV      ', 1, ' ')
      call add_default ('UU      ', 1, ' ')
      call add_default ('Z3      ', 1, ' ')
      call add_default ('TMQ     ', 1, ' ')
      call add_default ('PSL     ', 1, ' ')

      call add_default ('ICEFRAC ', 1, ' ')
      call add_default ('OCNFRAC ', 1, ' ')
      call add_default ('LANDFRAC ', 1, ' ')

#if ( defined COUP_SOM )
      call add_default ('MELTB   ', 1, ' ')
      call add_default ('MELTT   ', 1, ' ')
      call add_default ('MELTL   ', 1, ' ')
      call add_default ('GROWB   ', 1, ' ')
      call add_default ('FRAZIL  ', 1, ' ')
      call add_default ('FLOOD   ', 1, ' ')
      call add_default ('FRZMLT  ', 1, ' ')
      call add_default ('QFLUX   ', 1, ' ')
      call add_default ('QFLUX_FT', 1, ' ')
      call add_default ('QFLUX_TH', 1, ' ')
      call add_default ('QFLUX_A2', 1, ' ')
      call add_default ('FOCN    ', 1, ' ')
      call add_default ('EICEIN  ', 1, ' ')
      call add_default ('EICEOUT ', 1, ' ')
      call add_default ('F_ICE   ', 1, ' ')
      call add_default ('F_OCN   ', 1, ' ')
      call add_default ('FRZMLTMX', 1, ' ')
      call add_default ('DELTAICE', 1, ' ')
      call add_default ('IMBAL   ', 1, ' ')
      call add_default ('NRGERROR', 1, ' ')
      call add_default ('MLDANN  ', 1, ' ')
      call add_default ('ONF     ', 1, ' ')
      call add_default ('OIE     ', 1, ' ')
      call add_default ('OIERATE ', 1, ' ')
      call add_default ('NRGICE  ', 1, ' ')
      call add_default ('IIERATE ', 1, ' ')
#endif
      call add_default ('FSNSOI', 1, ' ')
      call add_default ('FLNSOI', 1, ' ')
      call add_default ('LHFLXOI', 1, ' ')
      call add_default ('SHFLXOI', 1, ' ')

! These fields removed 10/30/2003 per CRB decision.
!      call add_default ('CLDST   ', 1, ' ')
!      call add_default ('CNVCLD  ', 1, ' ')
!      call add_default ('EVAPPCT ', 1, ' ')
!      call add_default ('MQ      ', 1, ' ')
!      call add_default ('O3VMR   ', 1, ' ')
!      call add_default ('OMEGAU  ', 1, ' ')
!      call add_default ('PRECTMX ', 1, ' ')
!      call add_default ('QPERT   ', 1, ' ')
!      call add_default ('SICTHK  ', 1, ' ')
!      call add_default ('TPERT   ', 1, ' ')
!      call add_default ('TREFMNAV', 1, ' ')
!      call add_default ('TREFMXAV', 1, ' ')
!      call add_default ('TSICE   ', 1, ' ')
!      call add_default ('TT      ', 1, ' ')
!      call add_default ('USTAR   ', 1, ' ')
!      call add_default ('VZ      ', 1, ' ')
!      call add_default ('WSPEED  ', 1, ' ')
!      call add_default ('ZMDQ    ', 1, ' ')
!      call add_default ('ZMDT    ', 1, ' ')
!      call add_default ('ZZ      ', 1, ' ')
!      call add_default ('FLNT    ', 2, ' ')
!      call add_default ('OMEGA500', 2, ' ')
!      call add_default ('OMEGA850', 2, ' ')
!      call add_default ('PRECT   ', 2, ' ')
!      call add_default ('PSL     ', 2, 'I')
!      call add_default ('QFLX    ', 2, ' ')
!      call add_default ('QREFHT  ', 2, ' ')
!      call add_default ('T300    ', 2, 'I')
!      call add_default ('T850    ', 2, 'I')
!      call add_default ('TMQ     ', 2, 'I')
!      call add_default ('TREFHT  ', 2, ' ')
!      call add_default ('TREFHTMN', 2, ' ')
!      call add_default ('TREFHTMX', 2, ' ')
!      call add_default ('TS      ', 2, ' ')
!      call add_default ('TSMN    ', 2, ' ')
!      call add_default ('TSMX    ', 2, ' ')
!      call add_default ('U200    ', 2, 'I')
!      call add_default ('U850    ', 2, 'I')
!      call add_default ('V200    ', 2, 'I')
!      call add_default ('V850    ', 2, 'I')
!      call add_default ('Z300    ', 2, 'I')
!      call add_default ('Z500    ', 2, 'I')
!      call add_default ('Z700    ', 2, 'I')

      return
   end subroutine h_default

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


   subroutine add_default (name, tindex, flag) 153,4
!----------------------------------------------------------------------- 
! 
! Purpose: Add a field to the default "on" list for a given history file
! 
! Method: 
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      character(len=*), intent(in) :: name  ! field name
      character(len=1), intent(in) :: flag  ! averaging flag

      integer, intent(in) :: tindex         ! history tape index
!
! Local workspace
!
      integer :: f            ! field index
      logical :: found        ! flag indicates field found in masterlist
!
! Check validity of input arguments
!
      if (tindex > ptapes) then
         write(6,*)'ADD_DEFAULT: tape index=', tindex, ' is too big'
         call endrun
      end if

      if (flag /= ' ' .and. flag /= 'A' .and. flag /= 'I' .and. &
          flag /= 'X' .and. flag /= 'M') then

         call endrun ('ADD_DEFAULT: unknown averaging flag='//flag)
      end if
!
! Look through master list for input field name.  When found, set active
! flag for that tape to true.  Also set averaging flag if told to use other
! than default.
!
      found = .false.
      do f=1,nfmaster
         if (trim(name) == trim(masterlist(f)%field%name)) then
            masterlist(f)%actflag(tindex) = .true.
            if (flag /= ' ') then
               masterlist(f)%avgflag(tindex) = flag
               select case (flag)
               case ('A')
                  masterlist(f)%time_op(tindex) = 'mean'
               case ('I')
                  masterlist(f)%time_op(tindex) = ' '
               case ('X')
                  masterlist(f)%time_op(tindex) = 'maximum'
               case ('M')
                  masterlist(f)%time_op(tindex) = 'minimum'
               case default
                  call endrun ('H_DEFAULT: unknown avgflag='//flag)
               end select
            end if
            found = .true.
            exit
         end if
      end do

      if (.not. found) then
         call endrun ('ADD_DEFAULT: field='//name//' not found')
      end if

      return
   end subroutine add_default

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


   subroutine h_override (t) 1,1
!----------------------------------------------------------------------- 
! 
! Purpose: Override default history tape contents for a specific tape
!
! Method: Copy the flag into the master field list
! 
!-----------------------------------------------------------------------
!
! Arguments
!
      integer, intent(in) :: t         ! history tape index
!
! Local workspace
!
      integer :: f                     ! index over fields
      character(len=1) :: avgflg       ! lcl equiv of avgflag_pertape(t) (to address xlf90 compiler bug)

      avgflg = avgflag_pertape(t)
      do f=1,nfmaster
         select case (avgflg)
         case ('A')
            masterlist(f)%avgflag(t) = avgflag_pertape(t)
            masterlist(f)%time_op(t) = 'mean'
         case ('I')
            masterlist(f)%avgflag(t) = avgflag_pertape(t)
            masterlist(f)%time_op(t) = ' '
         case ('X')
            masterlist(f)%avgflag(t) = avgflag_pertape(t)
            masterlist(f)%time_op(t) = 'maximum'
         case ('M')
            masterlist(f)%avgflag(t) = avgflag_pertape(t)
            masterlist(f)%time_op(t) = 'minimum'
         case default
            call endrun ('H_OVERRIDE: unknown avgflag='//avgflag_pertape(t))
         end select
      end do
   end subroutine h_override
         
!#######################################################################


   subroutine h_define (t) 1,176
!----------------------------------------------------------------------- 
! 
! Purpose: Define contents of history file t
! 
! Method: Issue the required netcdf wrapper calls to define the history file contents
! 
!-----------------------------------------------------------------------
      use pspect
      use rgrid
      use commap
      use time_manager, only: get_step_size, get_ref_date
      use filenames,    only: caseid

#include <comctl.h>
#include <comhyb.h>

!-----------------------------------------------------------------------
!
! Input arguments
!
      integer, intent(in) :: t   ! tape index
!
! Local workspace
!
      integer :: i, j            ! longitude, latitude indices
      integer :: k               ! ISCCP vertical index
      integer :: l               ! ISCCP optical depth index
      integer :: kl              ! ISCCP merged k and l indices
      integer :: f               ! field index
      integer :: ff              ! varid index for fields output by column
      integer :: numlev          ! number of vertical levels (dimension and loop)
      integer :: ncreal          ! netCDF real data type
      integer :: dtime           ! timestep size
      integer :: ndbase = 0      ! days component of base time
      integer :: nsbase = 0      ! seconds component of base time
      integer :: nbdate          ! base date in yyyymmdd format
      integer :: bdate           ! base date in yyyymmdd format
      integer :: nbsec           ! time of day component of base date [seconds]
      integer :: yr, mon, day    ! year, month, day components of a date

#ifdef STAGGERED
      integer :: slatdim         ! staggered latitude dimension
      integer :: slondim         ! staggered longitude dimension
      integer :: slatvar         ! variable id for staggered lat
      integer :: slonvar         ! variable id for staggered lon
      integer :: dimen4us(4)     ! dimension arrary for staggered U winds
      integer :: dimen4vs(4)     ! dimension arrary for staggered V winds
      integer :: wsid            ! Staggered latitude weight ID
      
      real(r8) slons( splon )    ! Staggered grid point array (lon)
      real(r8) slats( plat-1 )   ! Staggered grid point array (lat)
#endif

      real(r8) ailev(plevp)      ! interface level values
      real(r8) gausslat(plat)    ! gaussian latitudes
      real(r8) pie               ! 3.14159...
      
      character(len=max_chars) str ! character temporary 
!
! netcdf variables
!
      integer ret                ! function return value
      integer timdim             ! unlimited dimension id
      integer latvar             ! latitude variable id
      integer lonvar             ! longitude variable id
      integer glatvar(pflds)     ! column latitude variable id
      integer glonvar(pflds)     ! column longitude variable id
      integer rlonvar            ! reduced longitude variable id
      integer ps0var             ! variable id for PS0
      integer chardim            ! character dimension id
      
      real(r8) alon(plon)        ! longitude values (degrees)
      real(r8) alev(plev)        ! level values (pascals)
      real(r8) rlon(plon,plat)   ! reduced longitudes (degrees)
      real(r8) alat(plat)        ! latitude values (degrees)
      real(r8) prmid(npres)      ! pressure midpoints of ISCCP data
      real(r8) taumid(ntau)      ! optical depth midpoints of ISCCP data
      real(r8) prstau(npres*ntau)! prmid + taumid/1000

      integer dimenchar(2)       ! character dimension ids
      integer dimen1(1)          ! dimension ids (1d)
      integer dimen2(2)          ! dimension ids (2d)
      integer dimen2t(2)         ! dimension ids (2d time boundaries)
      integer dimen3(3)          ! dimension ids (3d)
      integer dimen4f(4)         ! dimension ids (4d at levels)
      integer dimen4g(4)         ! temp array holding dimension ids for groups of contiguous columns (4d)
      integer dimen4i(4)         ! dimension ids (4d at interfaces)
      integer dimen4n(4)         ! dimension ids (4d at isccp pressure levels)
      integer londim             ! longitude dimension id
      integer latdim             ! latitude dimension id
      integer grouplondim(pflds) ! longitude dimension id
      integer grouplatdim(pflds) ! latitude dimension id
      integer levdim             ! level dimension id
      integer ilevdim            ! interface dimension id
      integer isccp_prs_dim      ! dimension variable for ISCCP pressure levels
      integer isccp_tau_dim      ! dimension variable for ISCCP tau values
      integer isccp_prstau_dim   ! dimension variable for ISCCP pressure*tau levels
      integer tbnddim            ! time_bnds dimension id
      integer levvar             ! level variable id
      integer ilevvar            ! intfc variable id
      integer isccp_prs_var      ! ISCCP mean pressure variable id
      integer isccp_tau_var      ! ISCCP mean optical depth variable id
      integer isccp_prstau_var   ! ISCCP mixed variable id
      integer old_mode           ! returned mode from netcdf call
      integer hyaiid             ! hybrid A coef. intfc var id
      integer hybiid             ! hybrid B coef. intfc var id
      integer hyamid             ! hybrid A coef. level var id
      integer hybmid             ! hybrid B coef. level var id
      integer ntrmid             ! M truncation parameter var id
      integer ntrnid             ! N truncation parameter var id
      integer ntrkid             ! K truncation parameter var id
!
      write(6,*)'Opening netcdf history file ', trim(nhfil(t))
      call wrap_create (nhfil(t), nf_clobber, nfid(t))
!
! Setup netcdf file - create the dimensions of lat,lon,time,level
!     
      ret = nf_set_fill (nfid(t), nf_nofill, old_mode)
      ret = nf_def_dim (nfid(t), 'lat', plat, latdim)
      ret = nf_def_dim (nfid(t), 'lon', plon, londim)
      
#ifdef STAGGERED
      ret = nf_def_dim (nfid(t), 'slat', plat-1, slatdim)
      ret = nf_def_dim (nfid(t), 'slon', splon, slondim)
#endif

      ret = nf_def_dim (nfid(t), 'lev', plev, levdim)
      ret = nf_def_dim (nfid(t), 'ilev', plevp, ilevdim)
      ret = nf_def_dim (nfid(t), 'isccp_prs', npres, isccp_prs_dim)
      ret = nf_def_dim (nfid(t), 'isccp_tau', ntau, isccp_tau_dim)
      ret = nf_def_dim (nfid(t), 'isccp_prstau', npres*ntau, isccp_prstau_dim)
      ret = nf_def_dim (nfid(t), 'time', nf_unlimited, timdim)
      ret = nf_def_dim (nfid(t), 'tbnd', 2, tbnddim)
      ret = nf_def_dim (nfid(t), 'chars', 8, chardim)
!
! create dimensions for groups of contiguous columns
! variables with for single columns are created later in this routine
!
      if (ngroup(t).ne.0) then
         do i = 1, ngroup(t)
            ret = nf_def_dim (nfid(t), tape(t)%column(i)%lat_name, tape(t)%column(i)%num_lats, grouplatdim(i))
            call wrap_def_var (nfid(t),tape(t)%column(i)%lat_name,nf_double,1,grouplatdim(i),glatvar(i))
            call wrap_put_att_text (nfid(t), glatvar(i),'long_name','latitude')
            call wrap_put_att_text (nfid(t), glatvar(i),'units','degrees_north')
            ret = nf_def_dim (nfid(t), tape(t)%column(i)%lon_name, tape(t)%column(i)%num_lons, grouplondim(i))
            call wrap_def_var (nfid(t),tape(t)%column(i)%lon_name,nf_double,1,grouplondim(i),glonvar(i))
            call wrap_put_att_text (nfid(t), glonvar(i),'long_name','longitude')
            call wrap_put_att_text (nfid(t), glonvar(i),'units','degrees_east')
         end do
      end if
!
! setup dimension arrays for 1,2,3,4d variables 
!     
      dimen1(1) = timdim

      dimen2(1) = londim
      dimen2(2) = latdim
      
      dimen2t(1) = tbnddim
      dimen2t(2) = timdim
      
      dimen3(1) = londim
      dimen3(2) = latdim
      dimen3(3) = timdim
!            
      dimen4f(1) = londim
      dimen4f(4) = timdim

      dimen4i(1) = londim
      dimen4i(4) = timdim

      dimen4f(2) = latdim
      dimen4f(3) = levdim
         
      dimen4i(2) = latdim
      dimen4i(3) = ilevdim

      dimen4n(1) = londim
      dimen4n(2) = latdim
      dimen4n(3) = isccp_prstau_dim
      dimen4n(4) = timdim

#ifdef STAGGERED
      dimen4us(1) = londim
      dimen4us(2) = slatdim
      dimen4us(3) = levdim
      dimen4us(4) = timdim
      
      dimen4vs(1) = slondim
      dimen4vs(2) = latdim
      dimen4vs(3) = levdim
      dimen4vs(4) = timdim
#endif

! define variables to label the dimensions, use the same names as the
! dimensions

      call wrap_def_var (nfid(t),'P0',nf_double,0,0,ps0var)
      str = 'reference pressure'
      call wrap_put_att_text (nfid(t), ps0var, 'long_name', str)
      call wrap_put_att_text (nfid(t), ps0var, 'units', 'Pa')
      
      call wrap_def_var (nfid(t),'lat',nf_double,1,LATDIM,latvar)
      call wrap_put_att_text (nfid(t), latvar, 'long_name', 'latitude')
      call wrap_put_att_text (nfid(t), latvar, 'units', 'degrees_north')
      
      if (fullgrid) then

         call wrap_def_var (nfid(t),'lon',nf_double,1,LONDIM,lonvar)
         call wrap_put_att_text (nfid(t), lonvar,'long_name','longitude')
         call wrap_put_att_text (nfid(t), lonvar,'units','degrees_east')

      else

         call wrap_def_var (nfid(t),'rlon',nf_double,2,dimen2,rlonvar)
         call wrap_put_att_text (nfid(t), rlonvar, 'long_name', 'reduced_longitude')
         call wrap_put_att_text (nfid(t), rlonvar,'units','degrees_east')

      end if

! If a staggered grid is in use, output the lat and lon arrays variables.
! Note: the staggered grid is currently independent of the reduced grid,
! and the staggered grid is a full grid.

#ifdef STAGGERED
      call wrap_def_var (nfid(t), 'slat', nf_double,1,SLATDIM,slatvar)
      call wrap_put_att_text (nfid(t), slatvar, 'long_name', 'staggered latitude')
      call wrap_put_att_text (nfid(t), slatvar, 'units', 'degrees_north')

      call wrap_def_var (nfid(t),'slon',nf_double,1,SLONDIM,slonvar)
      call wrap_put_att_text (nfid(t), slonvar, 'long_name', 'staggered longitude')
      call wrap_put_att_text(nfid(t), slonvar, 'units', 'degrees_east')
      
      call wrap_def_var (nfid(t),'w_stag',nf_double,1,slatdim,wsid)
      str = 'staggered latitude weights'
      call wrap_put_att_text (nfid(t), wsid, 'long_name', str)
#endif

      call wrap_def_var (nfid(t),'lev',nf_double,1,LEVDIM,levvar)
      str = 'hybrid level at midpoints (1000*(A+B))'
      call wrap_put_att_text (nfid(t), levvar, 'long_name', str)
      str = 'level'
      call wrap_put_att_text (nfid(t), levvar, 'units', str)
      call wrap_put_att_text (nfid(t), levvar, 'positive', 'down')
      call wrap_put_att_text (nfid(t), levvar, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate')
      call wrap_put_att_text (nfid(t), levvar, 'formula_terms', 'a: hyam b: hybm p0: P0 ps: PS')

      call wrap_def_var (nfid(t),'ilev',nf_double,1,ILEVDIM,ilevvar)
      str = 'hybrid level at interfaces (1000*(A+B))'
      call wrap_put_att_text (nfid(t), ilevvar, 'long_name', str)
      str = 'level'
      call wrap_put_att_text (nfid(t), ilevvar, 'units', str)
      call wrap_put_att_text (nfid(t), ilevvar, 'positive', 'down')
      call wrap_put_att_text (nfid(t), ilevvar, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate')
      call wrap_put_att_text (nfid(t), ilevvar, 'formula_terms', 'a: hyai b: hybi p0: P0 ps: PS')

! ISCCP pressure, optical depth, and mixed dimension

      call wrap_def_var (nfid(t), 'isccp_prs', NF_DOUBLE, 1, isccp_prs_dim, isccp_prs_var)
      str = 'Mean ISCCP pressure'
      call wrap_put_att_text (nfid(t), isccp_prs_var, 'long_name', str)
      str = 'mb'
      call wrap_put_att_text (nfid(t), isccp_prs_var, 'units', str)
      call wrap_put_att_realx (nfid(t), isccp_prs_var, 'isccp_prs_bnds', NF_DOUBLE, npres+1, prlim)

      call wrap_def_var (nfid(t), 'isccp_tau', NF_DOUBLE, 1, isccp_tau_dim, isccp_tau_var)
      str = 'Mean ISCCP optical depth'
      call wrap_put_att_text (nfid(t), isccp_tau_var, 'long_name', str)
      str = 'unitless'
      call wrap_put_att_text (nfid(t), isccp_tau_var, 'units', str)
      call wrap_put_att_realx (nfid(t), isccp_tau_var, 'isccp_tau_bnds', NF_DOUBLE, ntau+1, taulim)

      call wrap_def_var (nfid(t), 'isccp_prstau', NF_DOUBLE, 1, isccp_prstau_dim, isccp_prstau_var)
      str = 'Mean pressure (mb).mean optical depth (unitless)/1000'
      call wrap_put_att_text (nfid(t), isccp_prstau_var, 'long_name', str)
      str = 'mixed'
      call wrap_put_att_text (nfid(t), isccp_prstau_var, 'units', str)

      call get_ref_date(yr, mon, day, nbsec)
      nbdate = yr*10000 + mon*100 + day
      call wrap_def_var (nfid(t),'time',nf_double,1,TIMDIM,timeid(t))
      call wrap_put_att_text (nfid(t), timeid(t), 'long_name', 'time')
      str = 'days since ' // date2yyyymmdd(nbdate) // ' ' // sec2hms(nbsec)
      call wrap_put_att_text (nfid(t), timeid(t), 'units', str)
      call wrap_put_att_text (nfid(t), timeid(t), 'calendar', 'noleap')
      call wrap_put_att_text (nfid(t), timeid(t), 'bounds', 'time_bnds')

      call wrap_def_var (nfid(t),'time_bnds',nf_double,2,dimen2t,tbndid(t))
      call wrap_put_att_text (nfid(t), tbndid(t), 'long_name', 'time interval endpoints')
!
! Character
!
      dimenchar(1) = chardim
      dimenchar(2) = timdim
      call wrap_def_var (nfid(t),'date_written',NF_CHAR,2,dimenchar, date_writtenid(t))
      call wrap_def_var (nfid(t),'time_written',NF_CHAR,2,dimenchar, time_writtenid(t))
!
! Integer Header
!
      call wrap_def_var (nfid(t),'ntrm',NF_INT,0,0,ntrmid)
      str = 'spectral truncation parameter M'
      call wrap_put_att_text (nfid(t), ntrmid, 'long_name', str)

      call wrap_def_var (nfid(t),'ntrn',NF_INT,0,0,ntrnid)
      str = 'spectral truncation parameter N'
      call wrap_put_att_text (nfid(t), ntrnid, 'long_name', str)

      call wrap_def_var (nfid(t),'ntrk',NF_INT,0,0,ntrkid)
      str = 'spectral truncation parameter K'
      call wrap_put_att_text (nfid(t), ntrkid, 'long_name', str)

      call wrap_def_var (nfid(t),'ndbase',NF_INT,0,0,ndbaseid(t))
      str = 'base day'
      call wrap_put_att_text (nfid(t), ndbaseid(t), 'long_name', str)

      call wrap_def_var (nfid(t),'nsbase',NF_INT,0,0,nsbaseid(t))
      str = 'seconds of base day'
      call wrap_put_att_text (nfid(t), nsbaseid(t), 'long_name', str)

      call wrap_def_var (nfid(t),'nbdate',NF_INT,0,0,nbdateid(t))
      str = 'base date (YYYYMMDD)'
      call wrap_put_att_text (nfid(t), nbdateid(t), 'long_name', str)

#if ( defined BFB_CAM_SCAM_IOP )
      call wrap_def_var (nfid(t),'bdate',NF_INT,0,0,bdateid(t))
      str = 'base date (YYYYMMDD)'
      call wrap_put_att_text (nfid(t), bdateid(t), 'long_name', str)
#endif
      call wrap_def_var (nfid(t),'nbsec',NF_INT,0,0,nbsecid(t))
      str = 'seconds of base date'
      call wrap_put_att_text (nfid(t), nbsecid(t), 'long_name', str)

      call wrap_def_var (nfid(t),'mdt',NF_INT,0,0,mdtid)
      call wrap_put_att_text (nfid(t), mdtid, 'long_name', 'timestep')
      call wrap_put_att_text (nfid(t), mdtid, 'units', 's')

      call wrap_def_var (nfid(t),'nlon',NF_INT,1,latdim,nlonid(t))
      str = 'number of longitudes'
      call wrap_put_att_text (nfid(t), nlonid(t), 'long_name', str)

      call wrap_def_var (nfid(t),'wnummax',NF_INT,1,latdim,wnummaxid(t))
      str = 'cutoff Fourier wavenumber'
      call wrap_put_att_text (nfid(t), wnummaxid(t), 'long_name', str)
!
! Floating point time-invariant
!
      call wrap_def_var (nfid(t),'hyai',NF_DOUBLE,1,ilevdim,hyaiid)
      str = 'hybrid A coefficient at layer interfaces'
      call wrap_put_att_text (nfid(t), hyaiid, 'long_name', str)

      call wrap_def_var (nfid(t),'hybi',NF_DOUBLE,1,ilevdim,hybiid)
      str = 'hybrid B coefficient at layer interfaces'
      call wrap_put_att_text (nfid(t), hybiid, 'long_name', str)

      call wrap_def_var (nfid(t),'hyam',NF_DOUBLE,1,levdim,hyamid)
      str = 'hybrid A coefficient at layer midpoints'
      call wrap_put_att_text (nfid(t), hyamid, 'long_name', str)

      call wrap_def_var (nfid(t),'hybm',NF_DOUBLE,1,levdim,hybmid)
      str = 'hybrid B coefficient at layer midpoints'
      call wrap_put_att_text (nfid(t), hybmid, 'long_name', str)

      call wrap_def_var (nfid(t),'gw',NF_DOUBLE,1,latdim,gwid)
      str = 'gauss weights'
      call wrap_put_att_text (nfid(t), gwid, 'long_name', str)
!     
! Character header information 
!
      str = 'CF-1.0'
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'Conventions', str)
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'source', 'CAM')
#if ( defined BFB_CAM_SCAM_IOP )
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'CAM_GENERATED_FORCING','create SCAM IOP dataset')
#endif
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'case',caseid)
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'title',ctitle)
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'logname',logname)
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'host', host)
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'Version', &
           '$Name: cam3_0_brnchT_release01 $')
      call wrap_put_att_text (nfid(T), NF_GLOBAL, 'revision_Id', &
           '$Id: history.F90,v 1.26.2.48.4.1 2004/05/20 18:36:01 mvr Exp $')
!
! Create variables for model timing and header information 
!
      call wrap_def_var (nfid(t),'ndcur   ',nf_int,1,dimen1,ndcurid(t))
      str = 'current day (from base day)'
      call wrap_put_att_text (nfid(t), ndcurid(t), 'long_name', str)

      call wrap_def_var (nfid(t),'nscur   ',nf_int,1,dimen1,nscurid(t))
      str = 'current seconds of current day'
      call wrap_put_att_text (nfid(t), nscurid(t), 'long_name', str)

      call wrap_def_var (nfid(t),'date    ',nf_int,1,dimen1,dateid(t))
      str = 'current date (YYYYMMDD)'
      call wrap_put_att_text (nfid(t), dateid(t), 'long_name', str)

      call wrap_def_var (nfid(t),'co2vmr  ',nf_double,1,dimen1,co2vmrid(t))
      str = 'co2 volume mixing ratio'
      call wrap_put_att_text (nfid(t), co2vmrid(t), 'long_name', str)

      call wrap_def_var (nfid(t),'datesec ',nf_int,1,dimen1, datesecid(t))
      str = 'current seconds of current date'
      call wrap_put_att_text (nfid(t), datesecid(t), 'long_name', str)

#if ( defined BFB_CAM_SCAM_IOP )
      call wrap_def_var (nfid(t),'tsec ',nf_int,1,dimen1, tsecid(t))
      str = 'current seconds of current date needed for scam'
      call wrap_put_att_text (nfid(t), tsecid(t), 'long_name', str)
#endif
      call wrap_def_var (nfid(t),'nsteph  ',nf_int,1,dimen1,nstephid(t))
      str = 'current timestep'
      call wrap_put_att_text (nfid(t), nstephid(t), 'long_name', str)
!
! Create variables and attributes for field list
!
      do f=1,nflds(t)
         numlev = tape(t)%hlist(f)%field%numlev
         if (tape(t)%hlist(f)%hwrt_prec == 8) then
            ncreal = nf_double
         else
            ncreal = nf_float
         end if
!
!  Create variables and atributes for fields written out as columns
!         
         if (ngroup(t).ne.0) then
            do i=1,ngroup(t)
               ff=(f-1)*ngroup(t)+i
               dimen4g(1)=grouplondim(i)               
               dimen4g(2)=grouplatdim(i)               
               dimen4g(4)=timdim

               if (numlev == npres*ntau .and. tape(t)%hlist(f)%field%flag_isccplev) then
                  dimen4g(3)=isccp_prstau_dim
                  call wrap_def_var(nfid(t), tape(t)%hlist(f)%field_column_name(i),ncreal,4,dimen4g,varid(ff,t))
               else if (numlev == 1) then
                  dimen4g(3)=timdim
                  call wrap_def_var(nfid(t),trim(tape(t)%hlist(f)%field_column_name(i)),ncreal,3,dimen4g,varid(ff,t))
                  
               else if (numlev == plev) then
                  dimen4g(3)=levdim
                  call wrap_def_var(nfid(t), trim(tape(t)%hlist(f)%field_column_name(i)), ncreal,4,dimen4g,varid(ff,t))
                  
               else if (numlev == plevp) then
                  dimen4g(3)=ilevdim
                  call wrap_def_var(nfid(t), trim(tape(t)%hlist(f)%field_column_name(i)), ncreal,4,dimen4g,varid(ff,t))
               else
                  write(6,*)'H_DEFINE: bad numlev=',numlev
                  call endrun
               end if

               if (.not. fullgrid .or. tape(t)%hlist(f)%field%flag_xyfill) then
                  call wrap_put_att_realx (nfid(t), varid(ff,t), '_FillValue', ncreal, 1, fillvalue)
               end if

!JR Add missing_value for nco operators

               if (tape(t)%hlist(f)%field%flag_xyfill) then
                  call wrap_put_att_realx (nfid(t), varid(ff,t), 'missing_value', ncreal, 1, fillvalue)
               end if
               
               str = tape(t)%hlist(f)%field%units
               if ( str(1:1) /= ' ' ) then
                  call wrap_put_att_text (nfid(t), varid(ff,t), 'units', str)
               end if
               
               str = tape(t)%hlist(f)%field%long_name
               call wrap_put_att_text (nfid(t), varid(ff,t), 'long_name', str)
!
! Assign field attributes defining valid levels and averaging info
!
               str = tape(t)%hlist(f)%time_op
               select case (str)
               case ('mean', 'maximum', 'minimum' )
                  call wrap_put_att_text (nfid(t), varid(ff,t),'cell_method', 'time: '//str)
               end select
            end do
!
!  else create variables and atributes for fields written out as a full model grid
!         
         else
         
            if (numlev == npres*ntau .and. tape(t)%hlist(f)%field%flag_isccplev) then

               call wrap_def_var(nfid(t), tape(t)%hlist(f)%field%name, ncreal,4,dimen4n,varid(f,t))

            else if (numlev == 1) then

               call wrap_def_var(nfid(t),tape(t)%hlist(f)%field%name,ncreal,3,dimen3,varid(f,t))

            else if (numlev == plev) then	

#ifdef STAGGERED
               select case (tape(t)%hlist(f)%field%name)
               case ('US')
                  call wrap_def_var(nfid(t),tape(t)%hlist(f)%field%name,ncreal,4,dimen4us,varid(f,t))
               case ('VS')
                  call wrap_def_var(nfid(t),tape(t)%hlist(f)%field%name,ncreal,4,dimen4vs,varid(f,t))
               case default
                  call wrap_def_var(nfid(t),tape(t)%hlist(f)%field%name,ncreal,4,dimen4f,varid(f,t))
               end select
#else
               call wrap_def_var(nfid(t), tape(t)%hlist(f)%field%name, ncreal,4,dimen4f,varid(f,t))
               
#endif
            else if (numlev == plevp) then	
               
               call wrap_def_var(nfid(t), tape(t)%hlist(f)%field%name, ncreal,4,dimen4i,varid(f,t))
            else
               write(6,*)'H_DEFINE: bad numlev=',numlev
               call endrun
            end if
!	
            if (.not. fullgrid .or. tape(t)%hlist(f)%field%flag_xyfill) then
               call wrap_put_att_realx (nfid(t), varid(f,t), '_FillValue', ncreal, 1, fillvalue)
            end if

!JR Add missing_value for nco operators

            if (tape(t)%hlist(f)%field%flag_xyfill) then
               call wrap_put_att_realx (nfid(t), varid(f,t), 'missing_value', ncreal, 1, fillvalue)
            end if

            str = tape(t)%hlist(f)%field%units
            if ( str(1:1) /= ' ' ) then
               call wrap_put_att_text (nfid(t), varid(f,t), 'units', str)
            end if

            str = tape(t)%hlist(f)%field%units
            if ( str(1:1) /= ' ' ) then
               call wrap_put_att_text (nfid(t), varid(f,t), 'units', str)
            end if
            
            str = tape(t)%hlist(f)%field%long_name
            call wrap_put_att_text (nfid(t), varid(f,t), 'long_name', str)
!
! Assign field attributes defining valid levels and averaging info
!
            str = tape(t)%hlist(f)%time_op
            select case (str)
            case ('mean', 'maximum', 'minimum' )
               call wrap_put_att_text (nfid(t), varid(f,t),'cell_method', 'time: '//str)
            end select
         endif
      end do
!
      ret = nf_enddef(nfid(t))
      write(6,*)'H_DEFINE: Successfully opened netcdf file '
!
! Write time-invariant portion of history header
!
      call wrap_put_var_realx (nfid(t), ps0var, ps0)

      pie = 4.*atan(1.)
      do j=1,plat
         gausslat(j) = (180./pie)*clat(j)
      end do
      call wrap_put_var_realx (nfid(t), latvar, gausslat)
!
! Coordinate var defined only in full grid: otherwise define 2-d var.
!
      if (fullgrid) then
         do i=1,plon
            alon(i) = (i-1) * 360.0 / plon
         end do
         call wrap_put_var_realx (nfid(t), lonvar, alon)

      else

         do j=1,plat
            rlon(:,j) = fillvalue
            do i=1,nlon(j)
               rlon(i,j) = (i-1) * 360.0 / nlon(j)
            end do
         end do
         call wrap_put_var_realx (nfid(t), rlonvar, rlon)
      end if

#ifdef STAGGERED

! Calculate latitudes for the staggered grid.

      do j = 1, plat-1

! Staggered latitudes fall in between the regular latitudes. Note that
! for Lin-Rood dynamics, gausslat latitudes are not gaussian!

         slats(j) = (180./pie) * clat_staggered(j)

      end do

! Sanity check.

      if (slats(plat-1) .gt. 90.0) then
         write(6,*) "H_DEFINE: WARNING - last staggered grid latitude"
         write(6,*) "was calculated to be: ", slats(plat-1)
         write(6,*) "Point has been redefined to be 90.0"
         slats(plat-1) = 90.0
      endif

      call wrap_put_var_realx (nfid(t), slatvar, slats)
      
      do i = 1, splon
         slons(i) = ((i-1) * 360.0 / splon) - 180.0/splon
      enddo

      call wrap_put_var_realx (nfid(t), slonvar, slons)
      call wrap_put_var_realx (nfid(t), wsid, w_staggered)
#endif
!
! write out coordinate variables for columns
!
!jt This needs to be fixed for reduced grid : note the bogus dimension 1 on londeg
!

      if (ngroup(t).ne.0) then
         do i = 1, ngroup(t)
            call wrap_put_var_realx (nfid(t), glonvar(i), londeg(tape(t)%column(i)%columnlon(1):tape(t)%column(i)%columnlon(2),1))
            call wrap_put_var_realx (nfid(t), glatvar(i), latdeg(tape(t)%column(i)%columnlat(1):tape(t)%column(i)%columnlat(2)))
         end do
      end if
!
! 0.01 converts Pascals to millibars
!
      alev(:plev) = 0.01*ps0*(hyam(:plev) + hybm(:plev))
      ailev(:plevp) = 0.01*ps0*(hyai(:plevp) + hybi(:plevp))

      do k=1,npres
         prmid(k) = 0.5*(prlim(k) + prlim(k+1))
      end do

      do k=1,ntau
         taumid(k) = 0.5*(taulim(k) + taulim(k+1))
      end do

!JR Kludgey way of combining pressure and optical depth into a single dimension:
!JR pressure in millibars will show up on the left side of the decimal point, 
!JR optical depth/1000 will show up on the right

      do k=1,npres
         do l=1,ntau
            kl = (k-1)*ntau + l
            prstau(kl) = prmid(k) + taumid(l)*0.001
         end do
      end do

      call wrap_put_var_realx (nfid(t), levvar, alev)
      call wrap_put_var_realx (nfid(t), ilevvar, ailev)
      call wrap_put_var_realx (nfid(t), hyaiid, hyai)
      call wrap_put_var_realx (nfid(t), hybiid, hybi)
      call wrap_put_var_realx (nfid(t), hyamid, hyam)
      call wrap_put_var_realx (nfid(t), hybmid, hybm)
      call wrap_put_var_realx (nfid(t), gwid, w)
      call wrap_put_var_realx (nfid(t), isccp_prs_var, prmid)
      call wrap_put_var_realx (nfid(t), isccp_tau_var, taumid)
      call wrap_put_var_realx (nfid(t), isccp_prstau_var, prstau)
      
      call wrap_put_var_int (nfid(t), ntrmid, ptrm)
      call wrap_put_var_int (nfid(t), ntrnid, ptrn)
      call wrap_put_var_int (nfid(t), ntrkid, ptrk)
      dtime = get_step_size()
      call wrap_put_var_int (nfid(t), mdtid, dtime)
!
! Model date info
!
      call wrap_put_var_int (nfid(t), ndbaseid(t), ndbase)
      call wrap_put_var_int (nfid(t), nsbaseid(t), nsbase)

      call wrap_put_var_int (nfid(t), nbdateid(t), nbdate)
#if ( defined BFB_CAM_SCAM_IOP )
      call wrap_put_var_int (nfid(t), bdateid(t), nbdate)
#endif
      call wrap_put_var_int (nfid(t), nbsecid(t), nbsec)
!
! Reduced grid info
!
      call wrap_put_var_int (nfid(t), nlonid(t), nlon)
      call wrap_put_var_int (nfid(t), wnummaxid(t), wnummax)
      
      return
   end subroutine h_define

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


character(len=10) function date2yyyymmdd (date),3

   implicit none
   
! Input arguments

   integer, intent(in) :: date

! Local workspace

   integer :: year    ! year of yyyy-mm-dd
   integer :: month   ! month of yyyy-mm-dd
   integer :: day     ! day of yyyy-mm-dd

   if (date < 0) then
      call endrun ('DATE2YYYYMMDD: negative date not allowed')
   end if

   year  = date / 10000
   month = (date - year*10000) / 100
   day   = date - year*10000 - month*100

   if (month < 1 .or. month > 12) then
      write(6,*)'DATE2YYYYMMDD: invalid month encountered:', month
      call endrun ()
   end if

   if (day < 1 .or. day > ndm(month)) then
      write(6,*)'DATE2YYYYMMDD: invalid day encountered:', day
      call endrun ()
   end if

   write(date2yyyymmdd,80) year, month, day
80 format(i4.4,'-',i2.2,'-',i2.2)
   return
end function date2yyyymmdd

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


character(len=8) function sec2hms (seconds),3

   implicit none
   
! Input arguments

   integer, intent(in) :: seconds

! Local workspace

   integer :: hours     ! hours of hh:mm:ss
   integer :: minutes   ! minutes of hh:mm:ss
   integer :: secs      ! seconds of hh:mm:ss

   if (seconds < 0 .or. seconds > 86400) then
      write(6,*)'SEC2HRS: bad input seconds:', seconds
      call endrun ()
   end if

   hours   = seconds / 3600
   minutes = (seconds - hours*3600) / 60
   secs    = (seconds - hours*3600 - minutes*60)

   if (minutes < 0 .or. minutes > 60) then
      write(6,*)'SEC2HRS: bad minutes = ',minutes
      call endrun ()
   end if

   if (secs < 0 .or. secs > 60) then
      write(6,*)'SEC2HRS: bad secs = ',secs
      call endrun ()
   end if

   write(sec2hms,80) hours, minutes, secs
80 format(i2.2,':',i2.2,':',i2.2)
   return
end function sec2hms

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


   subroutine h_normalize (f, t) 1,1
!----------------------------------------------------------------------- 
! 
! Purpose: Normalize fields on a history file by the number of accumulations
! 
! Method: Loop over fields on the tape.  Need averaging flag and number of
!         accumulations to perform normalization.
! 
!-----------------------------------------------------------------------
!
! Input arguments
!
      integer, intent(in) :: f       ! field index
      integer, intent(in) :: t       ! tape index
!
! Local workspace
!
      integer begver                 ! on-node vert start index
      integer endver                 ! on-node vert end index
      integer c                      ! chunk (or lat) index
      integer coldimin               ! column dimension of model array
      integer endi                   ! terminating column index
      integer begdim3                ! on-node chunk or lat start index
      integer enddim3                ! on-node chunk or lat end index  
      integer, pointer :: nacs(:)    ! accumulation counter
      integer :: ib, ie    ! beginning and ending indices of first dimension
      integer :: jb, je    ! beginning and ending indices of second dimension
      integer :: ieu, jeu  ! number of elements in each dimension
      integer k

      logical :: flag_xyfill         ! non-applicable xy points flagged with fillvalue
      character*1 avgflag            ! averaging flag

      type (hbuffer_2d) :: hbuf      ! Per field history buffer
      type (dim_index_2d) :: dimind  ! 2-D dimension index
      
      call t_startf ('h_normalize')

      begver  = tape(t)%hlist(f)%field%begver
      endver  = tape(t)%hlist(f)%field%endver
      avgflag = tape(t)%hlist(f)%avgflag
!
! normalize by number of accumulations for averaged case
!
      begdim3 = tape(t)%hlist(f)%field%begdim3
      enddim3 = tape(t)%hlist(f)%field%enddim3
      flag_xyfill = tape(t)%hlist(f)%field%flag_xyfill

      do c=begdim3,enddim3
         call assoc_hbuf2d_with_hbuf3d (hbuf, tape(t)%hlist(f)%hbuf, c)
         coldimin = tape(t)%hlist(f)%field%coldimin
         nacs => tape(t)%hlist(f)%nacs(:coldimin,c)
         endi =  tape(t)%hlist(f)%field%colperdim3(c)

         dimind = dim_index_2d (1,endi,begver,endver)
         ib = dimind%beg1
         ie = dimind%end1
         jb = dimind%beg2
         je = dimind%end2
         
         ieu = ie-ib+1
         jeu = je-jb+1

         if (flag_xyfill) then
            if (associated(hbuf%buf8)) then
               do k=1,jeu
                  where (nacs(1:ieu) == 0)
                     hbuf%buf8(1:ieu,k) = fillvalue
                  endwhere
               end do
            else if (associated(hbuf%buf4)) then
               do k=1,jeu
                  where (nacs(1:ieu) == 0)
                     hbuf%buf4(1:ieu,k) = fillvalue
                  endwhere
               end do
            end if
         end if
         
         if (avgflag == 'A') then
            if (associated(hbuf%buf8)) then
               do k=1,jeu
                  where (nacs(1:ieu) /= 0)
                     hbuf%buf8(1:ieu,k) = hbuf%buf8(1:ieu,k) / nacs(1:ieu)
                  endwhere
               end do
            else if (associated(hbuf%buf4)) then
               do k=1,jeu
                  where (nacs(1:ieu) /= 0)
                     hbuf%buf4(1:ieu,k) = hbuf%buf4(1:ieu,k) / nacs(1:ieu)
                  endwhere
               end do
            end if
         end if
      end do

      call t_stopf ('h_normalize')
      
      return
   end subroutine h_normalize

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


   subroutine h_zero (f, t) 1,1
!----------------------------------------------------------------------- 
! 
! Purpose: Zero out accumulation buffers for a tape
! 
! Method: Loop through fields on the tape
! 
!-----------------------------------------------------------------------

      integer, intent(in) :: f     ! field index
      integer, intent(in) :: t     ! tape index
!
! Local workspace
!
      integer begver               ! on-node vert start index
      integer endver               ! on-node vert end index
      integer c                    ! chunk index
      integer endi                 ! terminating column index
      integer begdim3              ! on-node chunk or lat start index
      integer enddim3              ! on-node chunk or lat end index  
      type (dim_index_3d) :: dimind ! 3-D dimension index
      
      call t_startf ('h_zero')

      begdim3 = tape(t)%hlist(f)%field%begdim3
      enddim3 = tape(t)%hlist(f)%field%enddim3
      begver  = tape(t)%hlist(f)%field%begver
      endver  = tape(t)%hlist(f)%field%endver

      do c=begdim3,enddim3
         endi   = tape(t)%hlist(f)%field%colperdim3(c)

         dimind = dim_index_3d (1,endi,begver,endver,c,c)
         call set_hbuf_section_to_val (tape(t)%hlist(f)%hbuf,dimind,0._r8)
         tape(t)%hlist(f)%nacs(:endi,c) = 0
      end do

      call t_stopf ('h_zero')

      return
   end subroutine h_zero

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


   subroutine dump_field (f, t) 1,24
!----------------------------------------------------------------------- 
! 
! Purpose: Write a variable to a history tape
! 
! Method: If SPMD, first gather the data to the master processor.  
!         Next, transpose the data to COORDS order (the default).
!         Finally, issue the netcdf call to write the variable
! 
!-----------------------------------------------------------------------
#if ( defined SPMD )
# if ( defined STAGGERED )
      use spmd_dyn, only: npes, compute_gsfactors, comm_y
      use pmgrid, only: myid_z, strip3dxzy, strip3dxzyp, strip2d
      use parutilitiesmodule, only: commglobal, pargatherreal, pargatherreal4
# else
      use spmd_dyn, only: npes, compute_gsfactors
# endif
      use mpishorthand
#endif
      use dycore, only: dycore_is
!-----------------------------------------------------------------------
!
! Input arguments
!
      integer, intent(in) :: t, f ! tape, field indices
!
! Local workspace
!
      integer :: numlev           ! number of vertical levels (dimension and loop)
      integer :: start(4)         ! array of starting field indices for the field in the netcdf file
      integer :: count(4)         ! array of count values for the field in the netcdf file
      integer :: columnlon1       ! tmp
      integer :: columnlon2       ! tmp
      integer :: columnlat1       ! tmp
      integer :: columnlat2       ! tmp

      type (hbuffer_3d) :: tmpxyzbuf ! temporary array to hold contiguous column data in COORDS orDer
      type (hbuffer_3d) :: xyzbuf ! temporary array to hold data in COORDS order
      type (hbuffer_3d) :: xzybuf ! temporary array holding full x, y, z dims

#ifdef SPMD
      integer :: coldimin         ! column dimension of model array
      integer :: numsend          ! number of items to be sent
      integer :: numrecv(0:npes-1)! number of items to be received
      integer :: displs(0:npes-1) ! displacement array
      integer :: numowned         ! number of items owned by this MPI task
      integer :: mpireal          ! MPI real data type
#endif
      integer :: ff,j               ! varid index for fields with groups of columns
      type (dim_index_3d) :: dimind  ! 3-D dimension index
      logical :: flag_xyfill         ! non-applicable xy points flagged with fillvalue
      character*1 avgflag            ! averaging flag

      avgflag     = tape(t)%hlist(f)%avgflag
      flag_xyfill = tape(t)%hlist(f)%field%flag_xyfill
      numlev      = tape(t)%hlist(f)%field%numlev
#ifdef HDEBUG
      write(6,*)'DUMP_FIELD: writing ',tape(t)%hlist(f)%field%name,masterproc
#endif
      start(:) = 0
      start(:) = 0

      start(1) = 1
      start(2) = 1
      
      count(1) = plon

      if (numlev == 1) then

         start(3) = nfils(t)
         count(2) = plat
         count(3) = 1
         
      else if (numlev > 1) then

         start(3) = 1
         start(4) = nfils(t)
         count(2) = plat
         count(3) = numlev
         count(4) = 1

      else

         write(6,*)'DUMP_FIELD: bad numlev=', numlev
         stop 999
         
      end if

#ifdef STAGGERED
      if (tape(t)%hlist(f)%field%name == 'US') then
         count(2) = plat - 1
      else if (tape(t)%hlist(f)%field%name == 'VS') then
         count(1) = splon
      endif
#endif

#ifdef HDEBUG
      write(6,*)'DUMP_FIELD: writing time indx ', nfils(t), ' field id', varid(f,t), &
                ' name ', tape(t)%hlist(f)%field%name,masterproc
      write(6,*)'start=',start
      write(6,*)'count=',count
      call print_memusage ('dump_field')
#endif
!
! Transpose to COORDS order
!
! Note: this may or may not work when outputting the US variable
! in Lin-Rood because it has one less latitude (plat-1).
! xyzbuf is the memory for field to be written.  Initialize entire array to fillvalue since 
! history module does not know how to map unused portions of the physics domain
!
      if (masterproc) then
         dimind = dim_index_3d (1,plon,1,numlev,1,plat)
         call allocate_hbuf (xzybuf,dimind,tape(t)%hlist(f)%hbuf_prec)
         dimind = dim_index_3d (1,plon,1,plat,1,numlev)
         ! fvitt - for US field want xyzbuf to start at 2nd latitude
         if (tape(t)%hlist(f)%field%name == 'US') dimind = &
           dim_index_3d (1,plon,2,plat,1,numlev)
         call allocate_hbuf (xyzbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
      else
         call assoc_hbuf_with_nothing (xzybuf,tape(t)%hlist(f)%hbuf_prec)
         call nullify_hbuf (xyzbuf)
      end if
!
! physics decomposition: convert from col,lev,chunk -> lon,lev,lat -> lon,lat,lev
! dynamics decomposition: convert from lon,lev,lat -> lon,lat,lev
!
      select case (tape(t)%hlist(f)%field%decomp_type)
      case (phys_decomp)
         call gather_chunk_to_field_hbuf (1, numlev, 1, plon, tape(t)%hlist(f)%hbuf, xzybuf)
      case (dyn_decomp)
#ifdef SPMD
         if (tape(t)%hlist(f)%hbuf_prec == 8) then
            mpireal = mpir8
         else
            mpireal = mpir4
         end if

         if ( dycore_is('LR') )then
# if ( defined STAGGERED )
! NEW LR CODING
            if (tape(t)%hlist(f)%hbuf_prec == 8) then
               select case (numlev)
               case (1)
                  if (myid_z .eq. 0) call pargatherreal(comm_y, 0,  &
                         tape(t)%hlist(f)%hbuf%buf8, strip2d, xzybuf%buf8)
               case (plev)
                  call pargatherreal(commglobal, 0, tape(t)%hlist(f)%hbuf%buf8, &
                         strip3dxzy, xzybuf%buf8)
               case (plevp)
                  call pargatherreal(commglobal, 0, tape(t)%hlist(f)%hbuf%buf8, &
                         strip3dxzyp, xzybuf%buf8)
               case default
                  write(6,*)'DUMP_FIELD: invalid number of levels=', numlev
                  call endrun ()
               end select
            else
               select case (numlev)
               case (1)
                  if (myid_z .eq. 0) call pargatherreal4(comm_y, 0,  &
                         tape(t)%hlist(f)%hbuf%buf4, strip2d, xzybuf%buf4)
               case (plev)
                  call pargatherreal4(commglobal, 0, tape(t)%hlist(f)%hbuf%buf4, &
                         strip3dxzy, xzybuf%buf4)
               case (plevp)
                  call pargatherreal4(commglobal, 0, tape(t)%hlist(f)%hbuf%buf4, &
                         strip3dxzyp, xzybuf%buf4)
               case default
                  write(6,*)'DUMP_FIELD: invalid number of levels=', numlev
                  call endrun ()
               end select
            endif
# endif
         else
            coldimin = tape(t)%hlist(f)%field%coldimin
            numowned = coldimin*numlev
            call compute_gsfactors (numowned, numsend, numrecv, displs)
            call mpigatherv_hbuf (tape(t)%hlist(f)%hbuf, numsend, mpireal, xzybuf, numrecv, &
                                  displs, mpireal, 0, mpicom)
         endif
#else
!
! Need to copy instead of pointer associate to keep phys_decomp and dyn_decomp consistent
!
         xzybuf = tape(t)%hlist(f)%hbuf       ! Overloaded assignment
#endif
      case default
         write(6,*) 'DUMP_FIELD: bad decomp_type=',tape(t)%hlist(f)%field%decomp_type
         call endrun
      end select

      if (masterproc) then
         dimind = dim_index_3d (1,plon,1,numlev,1,plat) ! xzybuf order
         if (tape(t)%hlist(f)%field%name == 'US') dimind = dim_index_3d (1,plon,1,numlev,2,plat)
         call xzy_to_xyz (xyzbuf, xzybuf, dimind)
!
! Max/min fields got initialized to -huge and huge, respectively.  If fillvalue was in play
! these might still exist and need to be changed to fillvalue before output
!
         if (flag_xyfill .and. (avgflag == 'X' .or. avgflag == 'M')) then
            call fill_unset (xyzbuf%buf4, xyzbuf%buf8, dimind)
         end if
!
!  Check whether we are writing out an entire field or partial field
!  If writing partial field determine whether it is a single column or group of contiguous columns
!

         if (ngroup(t) .gt. 0) then
            do i=1,ngroup(t)
               ff=(f-1)*ngroup(t)+i     ! this is the variable id for the columm (or groups of columns)
#ifdef HDEBUG
               write(6,*)'DUMP_FIELD: writing column time indx ', &
                    nfils(t), ' field id', varid(ff,t), &
                    ' name ', trim(tape(t)%hlist(f)%field_column_name(i)), &
                    'numlons=',tape(t)%column(i)%num_lons,'numlats=', &
                    tape(t)%column(i)%num_lats,' base name ', &
                    trim(tape(t)%hlist(f)%field%name),' columnlon=', &
                    tape(t)%column(i)%columnlon,' columnlat=', &
                    tape(t)%column(i)%columnlat
#endif
               dimind = dim_index_3d (1,tape(t)%column(i)%num_lons,1,tape(t)%column(i)%num_lats,1,numlev)
               ! fvitt - for US field want tmpxyzbuf to have one less latitude ??
               if (tape(t)%hlist(f)%field%name == 'US')                   &
                 dimind = dim_index_3d (1,tape(t)%column(i)%num_lons,1,   &
                                          tape(t)%column(i)%num_lats-1,1, &
                                          numlev)
               call allocate_hbuf (tmpxyzbuf,dimind,tape(t)%hlist(f)%hbuf_prec)
               columnlon1=tape(t)%column(i)%columnlon(1)
               columnlon2=tape(t)%column(i)%columnlon(2)
               columnlat1=tape(t)%column(i)%columnlat(1)
               columnlat2=tape(t)%column(i)%columnlat(2)
 
              if (tape(t)%hlist(f)%hbuf_prec == 8) then
                  tmpxyzbuf%buf8(1:tape(t)%column(i)%num_lons,1:tape(t)%column(i)%num_lats,:)= &
                       xyzbuf%buf8(columnlon1:columnlon2,columnlat1:columnlat2,:)
               else if(tape(t)%hlist(f)%hbuf_prec == 4) then
                  tmpxyzbuf%buf4(1:tape(t)%column(i)%num_lons,1:tape(t)%column(i)%num_lats,:)= &
                       xyzbuf%buf4(columnlon1:columnlon2,columnlat1:columnlat2,:)
               end if
               count(1)=tape(t)%column(i)%num_lons
               count(2)=tape(t)%column(i)%num_lats
#ifdef HDEBUG
               write(6,*)'DUMP_FIELD: writing column variable ', &
                    varid(ff,t),' start=',start,' count=',count
#endif
               call wrap_put_vara_hbuf (nfid(t), varid(ff,t), start, count, tmpxyzbuf)
               call deallocate_hbuf (tmpxyzbuf)
            end do
!
! else if writing full fields
!
         else 

#ifdef HDEBUG
            write(6,*)'DUMP_FIELD: writing full field for time indx ', nfils(t), ' field id', varid(f,t), &
                 ' name ', tape(t)%hlist(f)%field%name
#endif
            call wrap_put_vara_hbuf (nfid(t), varid(f,t), start, count, xyzbuf)
         end if

         call deallocate_hbuf (xzybuf)
         call deallocate_hbuf (xyzbuf)
      else
         call nullify_hbuf (xzybuf)
      end if

      return
   end subroutine dump_field

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


   subroutine wshist () 1,24
!----------------------------------------------------------------------- 
! 
! Purpose: Driver routine to write fields on history tape t
! 
! Method: For variables which do not need to be gathered (SPMD) just issue the netcdf call
!         For those that do need to be gathered, call "dump_field" to do the operation.
!         Finally, zero the history buffers for each field written.
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------

      use time_manager, only: get_nstep, get_curr_date, get_curr_time,get_step_size
      use ghg_surfvals, only: co2vmr
!
! Local workspace
!
      character(len=8) :: cdate  ! system date
      character(len=8) :: ctime  ! system time
      
      integer t, f               ! tape, field indices
      integer ret                ! return value from netcdf call
      integer start              ! starting index required by nf_put_vara
      integer count1             ! count values required by nf_put_vara
      integer startc(2)          ! start values required by nf_put_vara (character)
      integer countc(2)          ! count values required by nf_put_vara (character)
#ifdef HDEBUG
      integer begdim3
      integer enddim3
#endif
      
      integer :: yr, mon, day      ! year, month, and day components of a date
      integer :: nstep             ! current timestep number
      integer :: ncdate            ! current date in integer format [yyyymmdd]
      integer :: ncsec             ! current time of day [seconds]
      integer :: ndcur             ! day component of current time
      integer :: nscur             ! seconds component of current time
      real(r8) :: time             ! current time
      real(r8) :: tdata(2)         ! time interval boundaries
      character(len=nlen) :: fname ! Filename
      logical :: prev              ! Label file with previous date rather than current
#if ( defined BFB_CAM_SCAM_IOP )
      integer :: tsec             ! day component of current time
      integer :: dtime            ! seconds component of current time
#endif
!-----------------------------------------------------------------------

      nstep = get_nstep()
      call get_curr_date(yr, mon, day, ncsec)
      ncdate = yr*10000 + mon*100 + day
      call get_curr_time(ndcur, nscur)
!
! Write time-varying portion of history file header
!
      do t=1,mtapes
         if (nhtfrq(t) == 0) then
            hstwr(t) = nstep /= 0 .and. day == 1 .and. ncsec == 0
            prev = .true.
         else
            hstwr(t) = mod(nstep,nhtfrq(t)) == 0
            prev = .false.
         end if

         if (hstwr(t)) then
            if (masterproc) then
               write(6,100) nfils(t),t,yr,mon,day,ncsec
100            format('WSHIST: writing time sample ',i3,' to h-file ', &
                      i1,' DATE=',i4.4,'/',i2.2,'/',i2.2,' NCSEC=',i6)
               write(6,*)
!
! Starting a new volume => define the metadata
!
               if (nfils(t)==0) then
                  fname    = interpret_filename_spec( hfilename_spec(t), number=(t-1), &
                             prev=prev )
!
! Check that this new filename isn't the same as a previous or current filename
!
                  do f = 1, mtapes
                    if ( trim(fname) == trim(nhfil(f)) )then
                       write(6,*)'WSHIST: New filename same as old file = ', trim(fname)
                       write(6,*)'Is there an error in your filename specifiers?'
                       write(6,*)'hfilename_spec(', t, ') = ', hfilename_spec(t)
                       if ( t /= f )then
                          write(6,*)'hfilename_spec(', f, ') = ', hfilename_spec(f)
                       end if
                       call endrun
                    end if
                  end do
                  nhfil(t) = fname
                  write(6,*)'WSHIST: nhfil(',t,')=',trim(nhfil(t))
                  cpath(t) = trim(get_archivedir('hist')) // nhfil(t)
                  if ( len_trim(nfpath(t)) == 0 ) nfpath(t) = cpath(t)
                  call h_define (t)
               end if
            end if
            
            nfils(t) = nfils(t) + 1
            start = nfils(t)
            count1 = 1
      
            if (masterproc) then
               call wrap_put_vara_int (nfid(t), ndcurid(t),start, count1,ndcur)
               call wrap_put_vara_int (nfid(t), nscurid(t),start, count1,nscur)
               call wrap_put_vara_int (nfid(t), dateid(t),start, count1,ncdate)
               call wrap_put_vara_realx (nfid(t), co2vmrid(t),start, count1,co2vmr)
               call wrap_put_vara_int (nfid(t), datesecid(t),start,count1,ncsec)
#if ( defined BFB_CAM_SCAM_IOP )
               dtime = get_step_size()
               tsec=dtime*nstep
               call wrap_put_vara_int (nfid(t), tsecid(t),start,count1,tsec)
#endif
               call wrap_put_vara_int (nfid(t), nstephid(t),start, count1,nstep)
               time = ndcur + nscur/86400._r8
               call wrap_put_vara_realx (nfid(t), timeid(t), start, count1,time)

               startc(1) = 1
               startc(2) = nfils(t)
               countc(1) = 2
               countc(2) = 1
               tdata(1) = beg_time(t)
               tdata(2) = time
               call wrap_put_vara_realx (nfid(t), tbndid(t), startc, countc, tdata)
               beg_time(t) = time  ! update beginning time of next interval

               startc(1) = 1
               startc(2) = nfils(t)
               countc(1) = 8
               countc(2) = 1
               call datetime (cdate, ctime)
               call wrap_put_vara_text (nfid(t), date_writtenid(t), startc, countc, cdate)
               call wrap_put_vara_text (nfid(t), time_writtenid(t), startc, countc, ctime)
#ifdef HDEBUG
               ret = nf_sync (nfid(t))
#endif
            end if

!$OMP PARALLEL DO PRIVATE (F)

            do f=1,nflds(t)
               call h_normalize (f, t)  ! Normalized averaged fields and/or put fillvalue
            end do
!
! Write field to history tape.  Note that this is NOT threaded due to netcdf limitations
!
            call t_startf ('dump_field')
            do f=1,nflds(t)
#ifdef HDEBUG
               begdim3 = tape(t)%hlist(f)%field%begdim3
               enddim3 = tape(t)%hlist(f)%field%enddim3
               if (tape(t)%hlist(f)%hbuf_prec == 8) then
                  write(6,*)'WSHIST:',tape(t)%hlist(f)%field%name,'(1,1,1)=',tape(t)%hlist(f)%hbuf%buf8(1,1,begdim3),masterproc
               else
                  write(6,*)'WSHIST:',tape(t)%hlist(f)%field%name,'(1,1,1)=',tape(t)%hlist(f)%hbuf%buf4(1,1,begdim3),masterproc
               end if
#endif
               call dump_field (f, t)
            end do
            call t_stopf ('dump_field')
!
! Zero history buffers and accumulators now that the fields have been written.
!
! JR: recoded away from array syntax because pgf90 generated bad code. Found
! by compiling with -Mbounds.  Appeared to be only multitasking that caused
! problems.
!
!$OMP PARALLEL DO PRIVATE (F)

            do f=1,nflds(t)
               call h_zero (f, t)
            end do
      
         end if
      end do

   end subroutine wshist

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


   subroutine write_inithist () 1,226
!----------------------------------------------------------------------- 
! 
! Purpose: Write an initial dataset
! 
! Method: Make appropriate netcdf calls to put required fields on tape
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
      use pspect
      use comsrf
      use buffer
      use prognostics
      use rgrid
      use ioFileMod
      use commap
      use phys_buffer, only: pbuf, pbuf_next_tim_idx, pbuf_old_tim_idx, pbuf_get_fld_idx
#if ( defined SPMD )
      use mpishorthand
# if ( defined STAGGERED )
      use pmgrid, only: strip3dxzy, strip2d, myid_z
      use spmd_dyn, only: comm_y
      use parutilitiesmodule, only: commglobal, pargatherreal
# endif
#endif
      use dycore, only: dycore_is
      use phys_grid, only: gather_chunk_to_field
      use time_manager, only: get_step_size, get_nstep, get_ref_date, &
                              get_curr_date, get_curr_time

      use filenames,    only: caseid
#if ( ! defined COUP_SOM ) && ( ! defined COUP_CSM )
      use sst_data, only: sst
#endif

#include <comctl.h>
#include <comhyb.h>
!
! Local variables
!
      character*80 name              ! variable name
      character*80 filename          ! filename to write
      character*256 path             ! mass store path
      character str*(max_chars)      ! string to hold netcdf attribute
      character*8 cdate              ! date
      character*8 ctime              ! time
      
      integer :: dtime               ! timestep size
      integer :: nstep               ! current timestep number
      integer :: ndbase = 0          ! days component of base time
      integer :: nsbase = 0          ! seconds component of base time
      integer :: nbdate              ! base date in integer format [yyyymmdd]
      integer :: nbsec               ! time of day relative to base date [seconds]
      integer :: ncdate              ! current date in integer format [yyyymmdd]
      integer :: ncsec               ! time of day relative to current date [seconds]
      integer :: ndcur               ! days component of current time
      integer :: nscur               ! seconds component of current time
      integer :: yr, mon, day        ! year, month, day components of a date

      real(r8) alon(plon)            ! longitude values
      real(r8) rlon(plon,plat)       ! reduced longitude values
      real(r8) alev(plev)            ! vertical levels
      real(r8) ailev(plevp)          ! vertical interfaces
      real(r8) gausslat(plat)        ! gaussian latitudes
      real(r8) pie                   ! 3.14...
      real(r8) time                  ! current time value
      real(r8) qpert_tmp (pcols,begchunk:endchunk)

      real(r8), allocatable :: arr2d(:,:)   ! 2d array to output
      real(r8), allocatable :: arr3d(:,:,:) ! 3d array to output
      real(r8), allocatable :: arr2dm(:,:,:)! multicomponent 2d array to output
      real(r8), allocatable :: arr2dchnk(:,:)   ! 2d array to output
      real(r8), allocatable :: arr2dmchnk(:,:,:)! multicomponent 2d array to output
!
! pointer to physics buffer fields
!
      integer itim, ifld
      real(r8), pointer, dimension(:,:,:) :: arr3dchnk_ptr

      integer i, m, j                ! indices
      integer ret                    ! return from netcdf lib call

      integer id                     ! netcdf file id
      integer old_mode               ! returned from nf_set_fill
      integer londim                 ! longitude dimension id
      integer latdim                 ! latitude dimension id
      integer levdim                 ! level dimension id
      integer ilevdim                ! interface dimension id
      integer timdim                 ! time dimension id
      integer dimenchar(2)           ! character dimension ids
      integer dimen1                 ! 1d dimension id
      integer dimen2(2)              ! 2d dimension ids
      integer dimen3(3)              ! 3d dimension ids
      integer dimen4(4)              ! 4d dimension ids
      integer ps0var                 ! ps0 variable id
      integer chardim                ! character dimension id
      integer lonvar                 !----------------------------------------------------
      integer rlonvar                ! 
      integer latvar                 ! dimension variable ids
      integer levvar                 ! 
      integer ilevvar                !----------------------------------------------------

      integer timeid                 !------------------------------------------------------
      integer date_writtenid         ! 
      integer time_writtenid         ! 
      integer ntrmid                 ! 
      integer ntrnid                 ! 
      integer ntrkid                 ! 
      integer ndbaseid               ! 
      integer nsbaseid               ! 
      integer nbdateid               ! 
      integer nbsecid                ! 
      integer mdtid                  ! 
      integer hyaiid                 !  
      integer hybiid                 ! 
      integer hyamid                 ! 
      integer hybmid                 ! 
      integer tracid(pcnst+pnats)    ! netcdf variable ids for variables of the given name
      integer gwid                   ! with "id" tacked on to the end
      integer ndcurid                ! 
      integer nscurid                ! 
      integer dateid                 ! 
      integer datesecid              ! 
      integer nstephid               ! 
      integer uid                    ! 
      integer vid                    ! 
      integer tid                    ! 
      integer pblhtid                ! 
      integer tpertid                ! 
      integer qpertid                ! 
      integer cldid                  ! 
      integer qcwatid                ! 
      integer tcwatid                ! 
      integer lcwatid                ! 
      integer qid                    ! 
      integer psid                   ! 
      integer phisid                 ! 
      integer sghid                  ! 
      integer landmid                ! 
      integer tsid                   ! 
      integer tsiceid                ! 
      integer tsice_rad_id           ! 
      integer tbotid                 ! 
      integer tssubid(4)             ! 
      integer sicthkid               !
      integer snowhlndid             !
      integer snowhiceid             !
      integer landfracid             ! land-fraction id
      integer icefracid              ! ice-fraction id
      integer tsocnid
      integer nlonid, wnummaxid      !------------------------------------------------------

      integer start(1)               ! starting index array for nf_put_vara (1d vars)
      integer count(1)               ! count array for nf_put_vara (1d vars)
      data start/1/
      data count/1/

      integer start2d(3)             ! starting index array for nf_put_vara (2d vars)
      integer count2d(3)             ! count array for nf_put_vara (2d vars)
      data start2d/1,1,1/
      data count2d/plon,plat,1/

      integer start3d(4)             ! starting index array for nf_put_vara (3d vars)
      integer count3d(4)             ! count array for nf_put_vara (3d vars)
      data start3d/1,1,1,1/
      data count3d/plon,plev,plat,1/

      logical haslon                 ! flag indicates longitude dimension is present
      logical haslat                 ! flag indicates latitude dimension is present

      integer n, nn                  ! indices
      integer nvars                  ! number of variables returned from nf_inq_nvars
      integer ndims                  ! number of dimensions returned from nf_inq_var
      integer natts                  ! number of attributes
      integer xtype                  ! variable type flag (netcdf)
      integer dimids(nf_max_var_dims)! dimension ids

#ifdef STAGGERED
      real(r8) slons(splon)
      real(r8) slats(plat-1)

      integer slonvar                ! Staggered longitude variable
      integer slatvar                ! Staggered latitude variable
      integer usid                   ! Staggered latidude wind variable id
      integer vsid                   !     "     longitude "      "     "
      integer slatdim                ! Dimension of the staggered lat array
      integer slondim                !     "     "   "     "      lon   "
      integer dimen4us(4)            ! NetCDF array holding array dimensions
      integer dimen4vs(4)            !    "
      integer count3dus(4)           !    "
      integer count3dvs(4)           !    "
      integer wsid                   ! Staggered latitude weight ID

      data count3dvs/splon,plev,plat,1/

      real(r8), allocatable :: arr3dxyz_local(:,:,:) ! help array
      real(r8), allocatable :: arr3dxzy_local(:,:,:) ! help array
      real(r8), allocatable :: arr3dxyz(:,:,:) ! help array
      real(r8), allocatable :: arr3dxzy(:,:,:) ! help array
      real(r8), allocatable :: arr3dus(:,:,:)
      real(r8), allocatable :: arr2dxy_local(:,:) ! help array
#endif

      integer startc(2)              ! starting index array for nf_put_vara (character)
      integer countc(2)              ! count array for nf_put_vara (character)

      integer k                      ! level index
      integer endi                   ! ending longitude index
      integer jcen                   ! latitude index (extended array)
      integer numi                   ! number of longitude values
      integer numperlat              ! number of values per latitude band

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

#ifdef STAGGERED
!
! Set count3dus in executable code since "plat-1" cannot be computed in a data stmt
!
      count3dus(1) = plon
      count3dus(2) = plev
      count3dus(3) = plat - 1
      count3dus(4) = 1
#endif
      if (masterproc) then

         dtime = get_step_size()
         nstep = get_nstep()
         call get_ref_date(yr, mon, day, nbsec)
         nbdate = yr*10000 + mon*100 + day
         call get_curr_date(yr, mon, day, ncsec)
         if (dycore_is ('SLD') ) then
            call get_curr_date(yr, mon, day, ncsec,offset=dtime)
         endif
         ncdate = yr*10000 + mon*100 + day
         call get_curr_time(ndcur, nscur)

         write(6,*)'WRITE_INITHIST: ncdate=',ncdate
         filename = interpret_filename_spec( ifilename_spec )
         path = trim(get_archivedir( 'init' )) // trim(filename)
         
         ret = nf_create (filename, nf_clobber, id)
         if (ret /= nf_noerr) call handle_error (ret)
         
         ret = nf_set_fill (id, nf_nofill, old_mode)
         ret = nf_def_dim (id, 'lat', plat, latdim)
         ret = nf_def_dim (id, 'lon', plon, londim)
         
#ifdef STAGGERED
         ret = nf_def_dim (id, 'slat', plat-1, slatdim)
         ret = nf_def_dim (id, 'slon', splon, slondim)
#endif

         ret = nf_def_dim (id, 'lev', plev, levdim)
         ret = nf_def_dim (id, 'ilev', plevp, ilevdim)
         ret = nf_def_dim (id, 'time', nf_unlimited, timdim)
         ret = nf_def_dim (id, 'chars', 8, chardim)
!
! setup dimension arrays for 1,2,3,4d variables 
!     
         dimen1 = timdim

         dimen2(1) = londim
         dimen2(2) = latdim

         dimen3(1) = londim
         dimen3(2) = latdim
         dimen3(3) = timdim
!            
         dimen4(1) = londim
         dimen4(2) = levdim
         dimen4(3) = latdim
         dimen4(4) = timdim

#ifdef STAGGERED
         dimen4us(1) = londim
         dimen4us(2) = levdim
         dimen4us(3) = slatdim
         dimen4us(4) = timdim

         dimen4vs(1) = slondim
         dimen4vs(2) = levdim
         dimen4vs(3) = latdim
         dimen4vs(4) = timdim
#endif
!           
! define variables to label the dimensions, use the same names as the
! dimensions
!
         call wrap_def_var (id,'P0',nf_double,0,0,ps0var)
         str = 'reference pressure'
         call wrap_put_att_text (id, ps0var, 'long_name', str)
         call wrap_put_att_text (id, ps0var, 'units', 'Pa')

         call wrap_def_var (id,'lat',nf_double,1,latdim,latvar)
         call wrap_put_att_text (id, latvar, 'long_name', 'latitude')
         call wrap_put_att_text (id, latvar, 'units', 'degrees_north')

         if (fullgrid) then

            call wrap_def_var (id,'lon',nf_double,1,londim,lonvar)
            call wrap_put_att_text (id, lonvar, 'long_name', 'longitude')
            call wrap_put_att_text (id, lonvar, 'units', 'degrees_east')

         else

            call wrap_def_var (id,'rlon',nf_double,2,dimen2,rlonvar)
            call wrap_put_att_text (id, rlonvar, 'long_name', 'reduced_longitude')
            call wrap_put_att_text (id, rlonvar, 'units', 'degrees_east')

         end if

! If a staggered grid is in use, output the lat and lon arrays variables.
! Note: the staggered grid is currently independent of any reduced grid
! information.

#ifdef STAGGERED
         call wrap_def_var (id, 'slat', nf_double, 1,slatdim, slatvar)
         call wrap_put_att_text (id, slatvar, 'long_name','staggered latitude')
         call wrap_put_att_text(id, slatvar, 'units', 'degrees_north')

         call wrap_def_var(id, 'slon', nf_double, 1, slondim, slonvar)
         call wrap_put_att_text (id, slonvar, 'long_name', 'staggered longitude')
         call wrap_put_att_text (id, slonvar, 'units', 'degrees_east')

         call wrap_def_var (id,'w_stag',nf_double,1,slatdim,wsid)
         str = 'staggered latitude weights'
         call wrap_put_att_text (id, wsid, 'long_name', str)
#endif

         call wrap_def_var (id,'lev',nf_double,1,levdim,levvar)
         str = 'hybrid level at midpoints (1000*(A+B))'
         call wrap_put_att_text (id, levvar, 'long_name', str)
         str = 'level'
         call wrap_put_att_text (id, levvar, 'units', str)
         call wrap_put_att_text (id, levvar, 'positive', 'down')
         call wrap_put_att_text (id, levvar, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate')
         call wrap_put_att_text (id, levvar, 'formula_terms', 'a: hyam b: hybm p0: P0 ps: PS')

         call wrap_def_var (id,'ilev',nf_double,1,ilevdim,ilevvar)
         str = 'hybrid level at interfaces (1000*(A+B))'
         call wrap_put_att_text (id, ilevvar, 'long_name', str)
         str = 'level'
         call wrap_put_att_text (id, ilevvar, 'units', str)
         call wrap_put_att_text (id, ilevvar, 'positive', 'down')
         call wrap_put_att_text (id, ilevvar, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate')
         call wrap_put_att_text (id, ilevvar, 'formula_terms', 'a: hyai b: hybi p0: P0 ps: PS')

         call wrap_def_var (id,'time',nf_double,1,timdim,timeid)
         call wrap_put_att_text (id, timeid, 'long_name', 'time')
         str = 'days since ' // date2yyyymmdd(nbdate) // ' ' // sec2hms(nbsec)
         call wrap_put_att_text (id, timeid, 'units', str)
         call wrap_put_att_text (id, timeid, 'calendar', 'noleap')
!
! Character
!
         dimenchar(1) = chardim
         dimenchar(2) = timdim
         call wrap_def_var (id,'date_written',nf_char,2,dimenchar, date_writtenid)
         call wrap_def_var (id,'time_written',nf_char,2,dimenchar, time_writtenid)
!
! Integer Header
!
         call wrap_def_var (id,'ntrm',nf_int,0,0,ntrmid)
         str = 'spectral truncation parameter M'
         call wrap_put_att_text (id, ntrmid, 'long_name', str)

         call wrap_def_var (id,'ntrn',nf_int,0,0,ntrnid)
         str = 'spectral truncation parameter N'
         call wrap_put_att_text (id, ntrnid, 'long_name', str)

         call wrap_def_var (id,'ntrk',nf_int,0,0,ntrkid)
         str = 'spectral truncation parameter K'
         call wrap_put_att_text (id, ntrkid, 'long_name', str)

         call wrap_def_var (id,'ndbase',nf_int,0,0,ndbaseid)
         str = 'base day'
         call wrap_put_att_text (id, ndbaseid, 'long_name', str)

         call wrap_def_var (id,'nsbase',nf_int,0,0,nsbaseid)
         str = 'seconds of base day'
         call wrap_put_att_text (id, nsbaseid, 'long_name', str)

         call wrap_def_var (id,'nbdate',nf_int,0,0,nbdateid)
         str = 'base date (YYYYMMDD)'
         call wrap_put_att_text (id, nbdateid, 'long_name', str)

         call wrap_def_var (id,'nbsec',nf_int,0,0,nbsecid)
         str = 'seconds of base date'
         call wrap_put_att_text (id, nbsecid, 'long_name', str)

         call wrap_def_var (id,'mdt',nf_int,0,0,mdtid)
         call wrap_put_att_text (id, mdtid, 'long_name', 'timestep')
         call wrap_put_att_text (id, mdtid, 'units', 's')

         if (.not. fullgrid) then
            call wrap_def_var (id,'nlon',nf_int,1,latdim,nlonid)
            str = 'number of longitudes'
            call wrap_put_att_text (id, nlonid, 'long_name', str)

            call wrap_def_var (id,'wnummax',nf_int,1,latdim,wnummaxid)
            str = 'cutoff Fourier wavenumber'
            call wrap_put_att_text (id, wnummaxid, 'long_name', str)
         end if
!
! Floating point time-invariant
!
         call wrap_def_var (id,'hyai',nf_double,1,ilevdim,hyaiid)
         str = 'hybrid A coefficient at layer interfaces'
         call wrap_put_att_text (id, hyaiid, 'long_name', str)

         call wrap_def_var (id,'hybi',nf_double,1,ilevdim,hybiid)
         str = 'hybrid B coefficient at layer interfaces'
         call wrap_put_att_text (id, hybiid, 'long_name', str)

         call wrap_def_var (id,'hyam',nf_double,1,levdim,hyamid)
         str = 'hybrid A coefficient at layer midpoints'
         call wrap_put_att_text (id, hyamid, 'long_name', str)

         call wrap_def_var (id,'hybm',nf_double,1,levdim,hybmid)
         str = 'hybrid B coefficient at layer midpoints'
         call wrap_put_att_text (id, hybmid, 'long_name', str)

         call wrap_def_var (id,'gw',nf_double,1,latdim,gwid)
         str = 'gauss weights'
         call wrap_put_att_text (id, gwid, 'long_name', str)
!     
! Character header information 
!
         str = 'CF-1.0'
         call wrap_put_att_text (id, nf_global, 'Conventions', str)
         call wrap_put_att_text (id, nf_global, 'source', 'CAM')
         call wrap_put_att_text (id, nf_global, 'case',caseid)
         call wrap_put_att_text (id, nf_global, 'title',ctitle)
         call wrap_put_att_text (id, nf_global, 'logname',logname)
         call wrap_put_att_text (id, nf_global, 'host', host)
!
! Create variables for model timing and header information 
!
         call wrap_def_var (id,'ndcur   ',nf_int,1,dimen1,ndcurid)
         str = 'current day (from base day)'
         call wrap_put_att_text (id, ndcurid, 'long_name', str)

         call wrap_def_var (id,'nscur   ',nf_int,1,dimen1,nscurid)
         str = 'current seconds of current day'
         call wrap_put_att_text (id, nscurid, 'long_name', str)

         call wrap_def_var (id,'date    ',nf_int,1,dimen1,dateid)
         str = 'current date (YYYYMMDD)'
         call wrap_put_att_text (id, dateid, 'long_name', str)

         call wrap_def_var (id,'datesec ',nf_int,1,dimen1, datesecid)
         str = 'current seconds of current date'
         call wrap_put_att_text (id, datesecid, 'long_name', str)

         call wrap_def_var (id,'nsteph  ',nf_int,1,dimen1,nstephid)
         str = 'current timestep'
         call wrap_put_att_text (id, nstephid, 'long_name', str)
!
! Call wrapper routine to add attributes as well as define variable for 
! fields which are in master field list
!

#ifdef STAGGERED
         call addvar (id,'US      ',nf_double,4,dimen4us,usid)
         call addvar (id,'VS      ',nf_double,4,dimen4vs,vsid)
#else
         call addvar (id,'U       ',nf_double,4,dimen4,uid)
         call addvar (id,'V       ',nf_double,4,dimen4,vid)
#endif
         call addvar (id,'T       ',nf_double,4,dimen4,tid)
         call addvar (id,'Q       ',nf_double,4,dimen4,qid)
         call addvar (id,'PS      ',nf_double,3,dimen3,psid)
         call addvar (id,'PHIS    ',nf_double,3,dimen3,phisid)
         call addvar (id,'SGH     ',nf_double,3,dimen3,sghid)
         call addvar (id,'LANDM_COSLAT',nf_double,3,dimen3,landmid)
         call addvar (id,'PBLH    ',nf_double,3,dimen3,pblhtid)
         call addvar (id,'TPERT   ',nf_double,3,dimen3,tpertid)
         call addvar (id,'QPERT   ',nf_double,3,dimen3,qpertid)
         call addvar (id,'CLOUD   ',nf_double,4,dimen4,cldid)
         call addvar (id,'QCWAT   ',nf_double,4,dimen4,qcwatid)
         call addvar (id,'TCWAT   ',nf_double,4,dimen4,tcwatid)
         call addvar (id,'LCWAT   ',nf_double,4,dimen4,lcwatid)

#if ( ! defined COUP_CSM )
         call addvar (id,'TSICERAD',nf_double,3,dimen3,tsice_rad_id)
         call addvar (id,'TS      ',nf_double,3,dimen3,tsid)
         call addvar (id,'TSICE   ',nf_double,3,dimen3,tsiceid)
         call addvar (id,'TS1     ',nf_double,3,dimen3,tssubid(1))
         call addvar (id,'TS2     ',nf_double,3,dimen3,tssubid(2))
         call addvar (id,'TS3     ',nf_double,3,dimen3,tssubid(3))
         call addvar (id,'TS4     ',nf_double,3,dimen3,tssubid(4))
         call addvar (id,'SNOWHICE',nf_double,3,dimen3,snowhiceid)
         call addvar (id,'LANDFRAC',nf_double,3,dimen3,landfracid)
         call addvar (id,'TBOT'    ,nf_double,3,dimen3,tbotid)
         call addvar (id,'ICEFRAC ',nf_double,3,dimen3,icefracid)
         call addvar (id,'SICTHK  ',nf_double,3,dimen3,sicthkid)
         call addvar (id,'TSOCN   ',nf_double,3,dimen3,tsocnid)
#endif
!
! Advected and non-advected constituents.
!
         do m=2,pcnst+pnats
            call addvar (id,cnst_name(m),nf_double,4,dimen4,tracid(m))
         end do
!
! Add fillvalue attribute to arrays which have lat and lon dimensions
!
         if (nf_inq_nvars (id, nvars) /= nf_noerr) then
            call endrun ('WRITE_INITHIST: bad call to nf_inq_nvars')
         end if

         do n=1,nvars
            call wrap_inq_var (id, n, name, xtype, ndims, dimids, natts)

            haslon = .false.
            haslat = .false.

            do nn=1,ndims
               if (dimids(nn) == londim) then
                  haslon = .true.
               else if (dimids(nn) == latdim) then
                  haslat = .true.
               end if

#ifdef STAGGERED
               if (dimids(nn) == slondim) then
                  haslon = .true.
               else if (dimids(nn) == slatdim) then
                  haslat = .true.
               end if
#endif

            end do

            if (haslon .and. haslat .and. .not. fullgrid) then
               call wrap_put_att_realx (id, n, '_FillValue', nf_double, 1, fillvalue)
            end if
         end do

         ret = nf_enddef(id)
         write(6,*)'WRITE_INITHIST: Successfully opened netcdf file'
!
! Write time-invariant portion of history header
!
         call wrap_put_var_realx (id, ps0var, ps0)

         ret = nf_put_vara_int (id, ndcurid,start, count,ndcur)
         ret = nf_put_vara_int (id, nscurid,start, count,nscur)
         ret = nf_put_vara_int (id, dateid,start, count,ncdate)
         ret = nf_put_vara_int (id, datesecid,start, count,ncsec)
         ret = nf_put_vara_int (id, nstephid,start, count,nstep)
         time = ndcur + nscur/86400._r8
         call wrap_put_vara_realx (id, timeid, start, count, time)

         startc(1) = 1
         startc(2) = 1
         countc(1) = 8
         countc(2) = 1
         call datetime (cdate, ctime)
         ret = nf_put_vara_text (id, date_writtenid, startc, countc, cdate)
         ret = nf_put_vara_text (id, time_writtenid, startc, countc, ctime)

         pie = 4.*atan(1.)
         do j=1,plat
            gausslat(j) = (180./pie)*clat(j)
         end do
         call wrap_put_var_realx (id, latvar, gausslat)
!
! Coordinate var defined only in full grid: otherwise define 2-d var.
!
         if (fullgrid) then
            do i=1,plon
               alon(i) = (i-1) * 360.0 / plon
            end do
            call wrap_put_var_realx (id, lonvar, alon)

         else

            do j=1,plat
               rlon(:,j) = fillvalue
               do i=1,nlon(j)
                  rlon(i,j) = (i-1) * 360.0 / nlon(j)
               end do
            end do
            call wrap_put_var_realx (id, rlonvar, rlon)

            call wrap_put_var_int (id, nlonid, nlon)
            call wrap_put_var_int (id, wnummaxid, wnummax)
         end if

#ifdef STAGGERED
! Calculate the staggered grid latitudes.
! Staggered latitudes fall in between the regular latitudes. Note that
! for Lin-Rood, gausslat latitudes are not gaussian!

         do j = 1, plat-1
            slats(j) = (gausslat(j) + gausslat(j+1))/2.
         end do

! Sanity check.

         if (slats(plat-1) .gt. 90.0) then
            write(6,*)"WRITE_INITHIST: WARNING - last staggered grid latitude"
            write(6,*) "was calculated to be:", slats(plat-1)
            write(6,*) "Point has been redefined to be 90.0"
            slats(plat-1) = 90.0
         endif

         call wrap_put_var_realx (id, slatvar, slats)

         do i = 1, splon
            slons(i) = ((i-1) * 360.0 / splon) - 180.0/splon
         enddo

         call wrap_put_var_realx (id, slonvar, slons)
         call wrap_put_var_realx (id, wsid, w_staggered)
#endif

!
! 0.01 converts Pascals to millibars
!
         alev(:plev) = 0.01*ps0*(hyam(:plev) + hybm(:plev))
         ailev(:plevp) = 0.01*ps0*(hyai(:plevp) + hybi(:plevp))

         call wrap_put_var_realx (id, levvar, alev)
         call wrap_put_var_realx (id, ilevvar, ailev)
         call wrap_put_var_realx (id, hyaiid, hyai)
         call wrap_put_var_realx (id, hybiid, hybi)
         call wrap_put_var_realx (id, hyamid, hyam)
         call wrap_put_var_realx (id, hybmid, hybm)
         call wrap_put_var_realx (id, gwid, w)

         call wrap_put_var_int (id, ntrmid, ptrm)
         call wrap_put_var_int (id, ntrnid, ptrn)
         call wrap_put_var_int (id, ntrkid, ptrk)
         call wrap_put_var_int (id, mdtid, dtime)
!
! Model date info
!
         call wrap_put_var_int (id, ndbaseid, ndbase)
         call wrap_put_var_int (id, nsbaseid, nsbase)
         call wrap_put_var_int (id, nbdateid, nbdate)
         call wrap_put_var_int (id, nbsecid, nbsec)

      end if  ! end of if-masterproc
!
! Allocate space for arrays to be written to initial history file.  In SPMD mode,
! masterproc needs all plat latitudes for writing
!
! If the staggered arrays are a different size than the regular arrays,
! yer gonna need something like "begsplat" and "endsplat" to define the
! portions used for MPI output below. -gg

      if (masterproc) then
         allocate (arr2d(plon,plat))
         allocate (arr3d(plon,plev,plat))
         allocate (arr2dchnk(pcols,begchunk:endchunk))
         allocate (arr2dmchnk(pcols,plevmx,begchunk:endchunk))
         allocate (arr2dm(plon,plevmx,plat))
#ifdef STAGGERED
         allocate (arr3dxyz_local(plon,beglat:endlat,beglev:endlev))
         allocate (arr3dxzy_local(plon,beglev:endlev,beglat:endlat))
         allocate (arr3dxyz(plon,plat,plev))
         allocate (arr3dxzy(plon,plev,plat))
         allocate (arr3dus(plon,plev,plat-1))
         allocate (arr2dxy_local(plon,beglat:endlat))
#endif
      else
         allocate (arr2d(plon,beglat:endlat))
         allocate (arr3d(plon,plev,beglat:endlat))
         allocate (arr2dchnk(pcols,begchunk:endchunk))
         allocate (arr2dmchnk(pcols,plevmx,begchunk:endchunk))
         allocate (arr2dm(1,1,1))
#ifdef STAGGERED
         allocate (arr3dxyz_local(plon,beglat:endlat,beglev:endlev))
         allocate (arr3dxzy_local(plon,beglev:endlev,beglat:endlat))
         allocate (arr3dxyz(1,1,1))
         allocate (arr3dxzy(1,1,1))
         allocate (arr3dus(1,1,1))
         allocate (arr2dxy_local(plon,beglat:endlat))
#endif
      end if

!----------------------------------------------------------------------------
! output variables
! 3-d
!----------------------------------------------------------------------------

      numperlat = plnlv
      arr3d = fillvalue

#ifdef STAGGERED
!
! output u
!
!
! WS 2001.01.29 :  Fixed Glenn's "DANGER, WILL ROBINSON! DANGER!" problem
! was:  call gather_put (id, usid, arr3dus, num, start3d, count3dus)
! Temporary: requires better solution
!

#if defined( SPMD )
!
!  Copy U3S (ghosted) to arr3dxyz_local (hack)
!
      do k=beglev,endlev
         do j=beglat,endlat
            do i=1,plon
               arr3dxyz_local(i,j,k) = u3s(i,j,k)
            enddo
         enddo
      enddo
      call gather( arr3dxyz_local, strip3dxyz, arr3dxyz )
#else
      do k=1,plev
         do j=1,plat
            do i=1,plon
               arr3dxyz(i,j,k) = u3s(i,j,k)
            enddo
         enddo
      enddo
#endif

      if (masterproc) then
         do j = 1, plat-1         ! shift the indices (what a mess...)
            do k = 1, plev
               do i = 1, plon
                  arr3dus(i,k,j) = arr3dxyz(i,j+1,k)
               enddo
            enddo
         enddo
         call wrap_put_vara_realx (id, usid, start3d, count3dus, arr3dus)
      end if

!
! output v
!
      do j = beglat, endlat
         do k = beglev, endlev
            arr3dxzy_local(:,k,j) = v3s(:,j,k)
         enddo
      enddo
      
#if defined (SPMD)
      call pargatherreal(commglobal, 0, arr3dxzy_local, strip3dxzy, arr3dxzy)
#else
      arr3dxzy(:,:,:) = arr3dxzy_local(:,:,:)
#endif
      if (masterproc) then
         call wrap_put_vara_realx (id, vsid, start3d, count3dvs, arr3dxzy)
      end if

#else

      arr3d(:,:,:) = fillvalue
      arr2d(:,:) = fillvalue

!
! output u
!
      do j=beglat,endlat
         endi = i1+nlon(j)-1
         numi = nlon(j)
         jcen = j1 - 1 + j
         arr3d(:numi,:plev,j) = u3(i1:endi,:plev,jcen,n3m1)
      end do
      call gather_put (id, uid, arr3d, numperlat, start3d, count3d)

!
! output v
!
      do j=beglat,endlat
         endi = i1+nlon(j)-1
         numi = nlon(j)
         jcen = j1 - 1 + j
         arr3d(:numi,:plev,j) = v3(i1:endi,:plev,jcen,n3m1)
      end do
      call gather_put (id, vid, arr3d, numperlat, start3d, count3d)

#endif

!
! output t
!
#ifdef STAGGERED
      do j=beglat,endlat
         do k=beglev,endlev
            arr3dxzy_local(:plon,k,j) = t3(:plon,k,j)
         enddo
      end do

#if defined (SPMD)
      call pargatherreal(commglobal, 0, arr3dxzy_local, strip3dxzy, arr3dxzy)
#else
      arr3dxzy(:,:,:) = arr3dxzy_local(:,:,:)
#endif

      if (masterproc) then
         call wrap_put_vara_realx (id, tid, start3d, count3d, arr3dxzy)
      end if
#else

      do j=beglat,endlat
         endi = i1+nlon(j)-1
         numi = nlon(j)
         jcen = j1 - 1 + j
         arr3d(:numi,:plev,j) = t3(i1:endi,:plev,jcen,n3m1)
      end do
      call gather_put (id, tid, arr3d, numperlat, start3d, count3d)
#endif

!
! output q (first component)
!
#ifdef STAGGERED
      do j=beglat,endlat
         if ( dycore_is('LR') )then
            do k=beglev,endlev
               do i=1,plon
                  arr3dxzy_local(i,k,j) = q3(i,j,k,1)
               enddo
            enddo
         else
            arr3dxzy_local(:plon,:plev,j) = q3(:plon,:plev,1,j)
         end if
      end do

#if defined (SPMD)
      call pargatherreal(commglobal, 0, arr3dxzy_local, strip3dxzy, arr3dxzy)
#else
      arr3dxzy(:,:,:) = arr3dxzy_local(:,:,:)
#endif

      if (masterproc) then
         call wrap_put_vara_realx (id, qid, start3d, count3d, arr3dxzy)
      end if
#else

      do j=beglat,endlat
         endi = i1+nlon(j)-1
         numi = nlon(j)
         jcen = j1 - 1 + j
         arr3d(:numi,:plev,j) = q3(i1:endi,:plev,1,jcen,n3m1)
      end do
      call gather_put (id, qid, arr3d, numperlat, start3d, count3d)
#endif

!
! output q (other components)
!
#ifdef STAGGERED
      do m=2,pcnst+pnats
         do j=beglat,endlat
            if ( dycore_is('LR') )then
               do k=beglev,endlev
                  do i=1,plon
                     arr3dxzy_local(i,k,j) = q3(i,j,k,m)
                  enddo
               enddo
            else
               arr3dxzy_local(:plon,:plev,j) = q3(:plon,:plev,m,j)
            endif
         end do

#if defined (SPMD)
         call pargatherreal(commglobal, 0, arr3dxzy_local, strip3dxzy, arr3dxzy)
#else
         arr3dxzy(:,:,:) = arr3dxzy_local(:,:,:)
#endif

         if (masterproc) then
            call wrap_put_vara_realx (id, tracid(m), start3d, count3d, arr3dxzy)
         end if
      end do
#else
      do m=2,pcnst+pnats
         do j=beglat,endlat
            endi = i1+nlon(j)-1
            numi = nlon(j)
            jcen = j1 - 1 + j
            arr3d(:numi,:plev,j) = q3(i1:endi,:plev,m,jcen,n3m1)
         end do
         call gather_put (id, tracid(m), arr3d, numperlat, start3d, count3d)
      end do
#endif

!
! 2-d fields
!
      numperlat = plon
      arr2d = fillvalue

#ifdef STAGGERED
      do j=beglat,endlat
         arr2dxy_local(:plon,j) = ps(:plon,j)
      end do

#if defined (SPMD)
      if (myid_z .eq. 0) call pargatherreal(comm_y, 0,      &
                  arr2dxy_local, strip2d, arr2d)
#else
      arr2d(:,:) = arr2dxy_local(:,:)
#endif

      if (masterproc) then
         call wrap_put_vara_realx (id, psid, start2d, count2d, arr2d)
      end if
#else

      do j=beglat,endlat
         numi = nlon(j)
         arr2d(:numi,j) = ps(:numi,j,n3m1)
      end do
      call gather_put (id, psid, arr2d, numperlat, start2d, count2d)
#endif

      do j=beglat,endlat
         numi = nlon(j)
#ifdef STAGGERED
         arr2dxy_local(:numi,j) = phis(:numi,j)
#else
         arr2d(:numi,j) = phis(:numi,j)
#endif
      end do

#ifdef STAGGERED
#if defined (SPMD)
      if (myid_z .eq. 0) call pargatherreal(comm_y, 0,      &
                  arr2dxy_local, strip2d, arr2d)
#else
      arr2d(:,:) = arr2dxy_local(:,:)
#endif

      if (masterproc) then
         call wrap_put_vara_realx (id, phisid, start2d, count2d, arr2d)
      end if
#else

      call gather_put (id, phisid, arr2d, numperlat, start2d, count2d)
#endif

      call gather_chunk_to_field(1,1,1,plon,sgh,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, sghid, start2d, count2d, arr2d)
      end if

      call gather_chunk_to_field(1,1,1,plon,landm,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, landmid, start2d, count2d, arr2d)
      end if
!
! 2d physics arrays
!
      call gather_chunk_to_field(1,1,1,plon,pblht,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, pblhtid, start2d, count2d, arr2d)
      end if

      call gather_chunk_to_field(1,1,1,plon,tpert,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, tpertid, start2d, count2d, arr2d)
      end if
         
      qpert_tmp(:,:) = qpert(:,1,:)
      call gather_chunk_to_field(1,1,1,plon,qpert_tmp,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, qpertid, start2d, count2d, arr2d)
      end if
!
! 3d physics arrays
!
      itim = pbuf_old_tim_idx()
      itim = pbuf_next_tim_idx(itim)
 
      ifld = pbuf_get_fld_idx('CLD')
      arr3dchnk_ptr => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,itim)
      call gather_chunk_to_field(1,plev,1,plon,arr3dchnk_ptr(:,:,:),arr3d)
      if (masterproc) then
         call wrap_put_vara_realx (id, cldid  , start3d, count3d, arr3d)
      end if

      ifld = pbuf_get_fld_idx('QCWAT')
      arr3dchnk_ptr => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,itim)
      call gather_chunk_to_field(1,plev,1,plon,arr3dchnk_ptr(:,:,:),arr3d)
      if (masterproc) then
         call wrap_put_vara_realx (id, qcwatid, start3d, count3d, arr3d)
      end if

      ifld = pbuf_get_fld_idx('TCWAT')
      arr3dchnk_ptr => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,itim)
      call gather_chunk_to_field(1,plev,1,plon,arr3dchnk_ptr(:,:,:),arr3d)
      if (masterproc) then
         call wrap_put_vara_realx (id, tcwatid, start3d, count3d, arr3d)
      end if

      ifld = pbuf_get_fld_idx('LCWAT')
      arr3dchnk_ptr => pbuf(ifld)%fld_ptr(1,1:pcols,1:pver,begchunk:endchunk,itim)
      call gather_chunk_to_field(1,plev,1,plon,arr3dchnk_ptr(:,:,:),arr3d)
      if (masterproc) then
         call wrap_put_vara_realx (id, lcwatid, start3d, count3d, arr3d)
      end if

#if ( ! defined COUP_CSM )
      call gather_chunk_to_field(1,1,1,plon,tsice_rad,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, tsice_rad_id, start2d, count2d, arr2d)
      end if

      call gather_chunk_to_field(1,1,1,plon,tsice,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, tsiceid, start2d, count2d, arr2d)
      end if

      do j=begchunk,endchunk
         arr2dchnk(:,j) = srfflx_state2d(j)%ts(:)
      end do
      call gather_chunk_to_field(1,1,1,plon,arr2dchnk,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, tsid, start2d, count2d, arr2d)
      end if


      do k=1,4
         do j=begchunk,endchunk
            arr2dmchnk(:,k,j) = surface_state2d(j)%tssub(:,k)
         end do
      end do
      call gather_chunk_to_field(1,plevmx,1,plon,arr2dmchnk,arr2dm)
      if (masterproc) then
         do k=1,plevmx
            do j=1,plat
               numi = nlon(j)
               arr2d(:numi,j) = arr2dm(:numi,k,j)
            end do
            call wrap_put_vara_realx (id, tssubid(k), start2d, count2d, arr2d)
         end do
      end if

      call gather_chunk_to_field (1,1,1,plon,sicthk,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, sicthkid, start2d, count2d, arr2d)
      end if
      call gather_chunk_to_field(1,1,1,plon,snowhice,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, snowhiceid, start2d, count2d, arr2d)
      end if
      call gather_chunk_to_field(1,1,1,plon,landfrac,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, landfracid, start2d, count2d, arr2d)
      end if
      do j=begchunk,endchunk
         arr2dchnk(:,j) = surface_state2d(j)%tbot(:)
      end do
      call gather_chunk_to_field(1,1,1,plon,arr2dchnk,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, tbotid, start2d, count2d, arr2d)
      end if
      call gather_chunk_to_field(1,1,1,plon,icefrac,arr2d)
      if (masterproc) then
         call wrap_put_vara_realx (id, icefracid, start2d, count2d, arr2d)
      end if
!
!JR Write out tsocn only in full-physics mode.  Problem is array sst has not been
!JR allocated in adiabatic and ideal_phys modes.
!
      if (.not. adiabatic .and. .not. ideal_phys) then
#if ( defined COUP_SOM )
         call gather_chunk_to_field(1,1,1,plon,tsocn,arr2d)
#else
         call gather_chunk_to_field(1,1,1,plon,sst,arr2d)
#endif
         if (masterproc) then
            call wrap_put_vara_realx (id, tsocnid, start2d, count2d, arr2d)
         end if
      end if
#endif

      if (masterproc) then
         ret = nf_close (id)
         call putfil (filename, path, mss_wpass, mss_irt, .true.)
      end if

      deallocate (arr2d)
      deallocate (arr3d)
      deallocate (arr2dm)
#ifdef STAGGERED
      deallocate (arr3dxyz_local)
      deallocate (arr3dxzy_local)
      deallocate (arr3dxyz)
      deallocate (arr3dxzy)
      deallocate (arr3dus)
      deallocate (arr2dxy_local)
#endif      
      return
   end subroutine write_inithist

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


   subroutine gather_put (fileid, fieldid, arr, numperlat, start, count) 7,4
!----------------------------------------------------------------------- 
! 
! Purpose: gather an array on the master processor and issue the proper netcdf call
!          to write it out
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
#if ( defined SPMD )
      use mpishorthand
      use spmd_dyn, only: npes, compute_gsfactors
#endif
!
! Input Arguments
!
      integer, intent(in) :: fileid    ! netcdf file id
      integer, intent(in) :: fieldid   ! netcdf field id
      integer, intent(in) :: numperlat ! number of array values to gather per latitude
      integer, intent(in) :: start(*)  ! array of start values for netcdf call
      integer, intent(in) :: count(*)  ! array of count values for netcdf call

      real(r8), intent(inout) :: arr(*)! array to be gathered then written

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

#ifdef SPMD
      integer :: numrecv(0:npes-1)     ! number of items to be received
      integer :: displs(0:npes-1)      ! displacement array
      integer :: numsend               ! number of items to send
      
      call compute_gsfactors (numperlat, numsend, numrecv, displs)
      call mpigatherv (arr, numsend, mpir8, arr, numrecv, displs, mpir8, 0, mpicom)
#endif

      if (masterproc) then
         call wrap_put_vara_realx (fileid, fieldid, start, count, arr)
      end if
      
      return
   end subroutine gather_put

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


   subroutine addvar (ncid, name, xtype , ndims , dimids, vid) 31,3
!----------------------------------------------------------------------- 
! 
! Purpose: Issue the netcdf call to add a variable to the dataset
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
!
! Input arguments
!
      integer, intent(in) :: ncid       ! netcdf file id
      integer, intent(in) :: xtype      ! netcdf type flag
      integer, intent(in) :: ndims      ! number of dimensions
      integer, intent(in) :: dimids(:)  ! dimension ids
      integer, intent(in) :: vid        ! variable ids
      
      character(len=*), intent(in) :: name ! variable name
!
! Local workspace
!
      integer f                         ! field index

      call wrap_def_var (ncid, name, xtype , ndims , dimids, vid)
      
      do f=1,nfmaster
         if (masterlist(f)%field%name == name) then
            call wrap_put_att_text (ncid, vid, 'long_name', masterlist(f)%field%long_name)
            if (masterlist(f)%field%units(1:1) /= ' ') then
               call wrap_put_att_text (ncid, vid, 'units', masterlist(f)%field%units)
            end if
            return
         end if
      end do
!
! Field not found in masterlist: long_name and units attributes unknown so just
! return
!
      return
   end subroutine addvar


   subroutine bldfld () 4,229

!----------------------------------------------------------------------- 
! 
! Purpose: 
!
! Build Master Field List of all possible fields in a history file.  Each field has 
! associated with it a "long_name" netcdf attribute that describes what the field is, 
! and a "units" attribute.
! 
! Method: Call a subroutine to add each field
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
      use ppgrid, only: pver, pverp
      use comsrf, only: plevmx, tsnam
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
!
! Local variables
!
      integer m,j                   ! Indices
      character(len=20) :: string   ! Character string to use
!
! Call addfld to add each field to the Master Field List.
!
      call addfld ('PHIS    ','m2/s2   ',1,    'I','Surface geopotential',phys_decomp)
      call addfld ('PS      ','Pa      ',1,    'A','Surface pressure',phys_decomp)
      call addfld ('PSDRY   ','Pa      ',1,    'A','Surface pressure',phys_decomp)
      call addfld ('T       ','K       ',plev, 'A','Temperature',phys_decomp)
      call addfld ('U       ','m/s     ',plev, 'A','Zonal wind',phys_decomp)
      call addfld ('V       ','m/s     ',plev, 'A','Meridional wind',phys_decomp)
      call addfld ('ETADOT  ','1/s     ',plevp,'A','Vertical (eta) velocity',dyn_decomp)
      call addfld ('SGH     ','m       ',1,    'I','Standard deviation of orography',phys_decomp)
      call addfld ('US      ','m/s     ',plev, 'A','Zonal wind, staggered',dyn_decomp)
      call addfld ('VS      ','m/s     ',plev, 'A','Meridional wind, staggered',dyn_decomp)
!
! Constituent tracers
!
      do m=1,pcnst+pnats
         call addfld (cnst_name(m),'kg/kg   ',pver, 'A',cnst_longname(m),phys_decomp)
         call addfld (sflxnam(m),  'kg/m2/s ',1,    'A',trim(cnst_name(m))//' surface flux',phys_decomp)
         call addfld (dcconnam(m), 'kg/kg/s ',pver, 'A',trim(cnst_name(m))//' tendency due to moist processes',phys_decomp)
      end do
      call addfld ('PDELDRY ','Pa      ',plev, 'A','Dry pressure difference between levels',phys_decomp)

      call addfld ('DUH     ','K/s     ',plev, 'A','U horizontal diffusive heating',dyn_decomp)
      call addfld ('DVH     ','K/s     ',plev, 'A','V horizontal diffusive heating',dyn_decomp)
      call addfld ('DTH     ','K/s     ',plev, 'A','T horizontal diffusive heating',dyn_decomp)
      call addfld ('SNOWHICE','m       ',1,    'A','Water equivalent snow depth',phys_decomp)
      call addfld ('SNOWHLND','m       ',1,    'A','Water equivalent snow depth',phys_decomp)
      call addfld ('SICTHK  ','m       ',1,    'A','Sea ice thickness',phys_decomp)
      call addfld ('PRECL   ','m/s     ',1,    'A','Large-scale (stable) precipitation rate',phys_decomp)
      call addfld ('PRECC   ','m/s     ',1,    'A','Convective precipitation rate',phys_decomp)
      call addfld ('PRECSH  ','m/s     ',1,    'A','Shallow Convection precipitation rate',phys_decomp)
      write(string,'(F4.2)') precc_thresh*3600.*1000.0
      call addfld ('PRECCFRQ','fraction',1,    'A',&
      'Convective precipitation frequency (fraction of time where rate is > '//trim(string)//'mm/hr)' &
      ,phys_decomp)
      call addfld ('PRECCINT','mm/hr   ',1,    'A',&
      'Convective precipitation rate (less than '//trim(string)// &
      'mm/hr is set to zero -- to get intensity divide by PRECCFRQ)',phys_decomp)
      write(string,'(F4.2)') precl_thresh*3600.*1000.0
      call addfld ('PRECLFRQ','fraction',1,    'A',&
      'Large-scale (stable) precipitation frequency (fraction of time where rate is > '//trim(string)//'mm/hr)' &
      ,phys_decomp)
      call addfld ('PRECLINT','mm/hr   ',1,    'A',&
      'Large-scale (stable) precipitation rate (less than '//trim(string)// &
      'mm/hr is set to zero -- to get intensity divide by PRECLFRQ)',phys_decomp)
      call addfld ('PRECT   ','m/s     ',1,    'A','Total (convective and large-scale) precipitation rate',phys_decomp)
      call addfld ('EVAPPCT ','percent ',1,    'A','Percentage of Zhang-McFarlane precipitation going into evaporation',phys_decomp)
      call addfld ('PRECTMX ','m/s     ',1,    'X','Maximum (convective and large-scale) precipitation rate',phys_decomp)
      call addfld ('PRECSL  ','m/s     ',1,    'A','Large-scale (stable) snow rate (water equivalent)',phys_decomp)
      call addfld ('PRECSC  ','m/s     ',1,    'A','Convective snow rate (water equivalent)',phys_decomp)
      call addfld ('SHFLX   ','W/m2    ',1,    'A','Surface sensible heat flux',phys_decomp)
      call addfld ('LHFLX   ','W/m2    ',1,    'A','Surface latent heat flux',phys_decomp)
      call addfld ('QFLX    ','kg/m2/s ',1,    'A','Surface water flux',phys_decomp)
      call addfld ('PBLH    ','m       ',1,    'A','PBL height',phys_decomp)
      call addfld ('USTAR   ','m/s     ',1,    'A','Surface friction velocity',phys_decomp)
      call addfld ('TREFHT  ','K       ',1,    'A','Reference height temperature',phys_decomp)
#if ( defined COUP_CSM )
      call addfld ('QREFHT  ','kg/kg   ',1,    'A','Reference height humidity',phys_decomp)
#endif
      call addfld ('CGH     ','K/m     ',pverp,'A','Counter-gradient term for heat in PBL',phys_decomp)
      call addfld ('CGQ     ','1/m     ',pverp,'A','Counter-gradient term for moisture in PBL',phys_decomp)
      call addfld ('CGS     ','s/m2    ',pverp,'A','Counter-gradient coeff on surface kinematic fluxes',phys_decomp)
      call addfld ('TPERT   ','K       ',1,    'A','Perturbation temperature (eddies in PBL)',phys_decomp)
      call addfld ('QPERT   ','kg/kg   ',1,    'A','Perturbation specific humidity (eddies in PBL)',phys_decomp)
      call addfld ('KVH     ','m2/s    ',pverp,'A','Vertical diffusion diffusivities (heat/moisture)',phys_decomp)
      call addfld ('KVM     ','m2/s    ',pverp,'A','Vertical diffusion diffusivities (momentum)',phys_decomp)
      call addfld ('DUV     ','m/s2    ',pver, 'A','U vertical diffusion',phys_decomp)
      call addfld ('DVV     ','m/s2    ',pver, 'A','V vertical diffusion',phys_decomp)
      call addfld ('DTV     ','K/s     ',pver, 'A','T vertical diffusion',phys_decomp)
      call addfld ('DTVKE   ','K/s     ',pver, 'A','dT/dt vertical diffusion KE dissipation',phys_decomp)
      call addfld ('FSNS    ','W/m2    ',1,    'A','Net solar flux at surface',phys_decomp)
      call addfld ('FLNS    ','W/m2    ',1,    'A','Net longwave flux at surface',phys_decomp)
      call addfld ('FLNT    ','W/m2    ',1,    'A','Net longwave flux at top of model',phys_decomp)
      call addfld ('FLUT    ','W/m2    ',1,    'A','Upwelling longwave flux at top of model',phys_decomp)
      call addfld ('FLUTC   ','W/m2    ',1,    'A','Clearsky upwelling longwave flux at top of model',phys_decomp)
      call addfld ('FSNT    ','W/m2    ',1,    'A','Net solar flux at top of model',phys_decomp)
      call addfld ('FSNTOA  ','W/m2    ',1,    'A','Net solar flux at top of atmosphere',phys_decomp)
      call addfld ('FSNTOAC ','W/m2    ',1,    'A','Clearsky net solar flux at top of atmosphere',phys_decomp)
      call addfld ('LWCF    ','W/m2    ',1,    'A','Longwave cloud forcing',phys_decomp)
      call addfld ('SWCF    ','W/m2    ',1,    'A','Shortwave cloud forcing',phys_decomp)
      call addfld ('CLOUD   ','fraction',pver, 'A','Cloud fraction',phys_decomp)
      call addfld ('FLNTC   ','W/m2    ',1,    'A','Clearsky net longwave flux at top of model',phys_decomp)
      call addfld ('FLN200  ','W/m2    ',1,    'A','Net longwave flux at 200 mb',phys_decomp)
      call addfld ('FLN200C ','W/m2    ',1,    'A','Clearsky net longwave flux at 200 mb',phys_decomp)
      call addfld ('FSN200  ','W/m2    ',1,    'A','Net shortwave flux at 200 mb',phys_decomp)
      call addfld ('FSN200C ','W/m2    ',1,    'A','Clearsky net shortwave flux at 200 mb',phys_decomp)
      call addfld ('FSNTC   ','W/m2    ',1,    'A','Clearsky net solar flux at top of model',phys_decomp)
      call addfld ('FLNSC   ','W/m2    ',1,    'A','Clearsky net longwave flux at surface',phys_decomp)
      call addfld ('FSNSC   ','W/m2    ',1,    'A','Clearsky net solar flux at surface',phys_decomp)
      call addfld ('FSDSC   ','W/m2    ',1,    'A','Clearsky downwelling solar flux at surface',phys_decomp)
      call addfld ('OMEGA   ','Pa/s    ',pver, 'A','Vertical velocity (pressure)',phys_decomp)
      call addfld ('DQP     ','kg/kg/s ',pver, 'A','Specific humidity tendency due to precipitation',phys_decomp)
      call addfld ('TAUX    ','N/m2    ',1,    'A','Zonal surface stress',phys_decomp)
      call addfld ('TAUY    ','N/m2    ',1,    'A','Meridional surface stress',phys_decomp)
      call addfld ('SRFRAD  ','W/m2    ',1,    'A','Net radiative flux at surface',phys_decomp)
      call addfld ('QRS     ','K/s     ',pver, 'A','Solar heating rate',phys_decomp)
      call addfld ('QRL     ','K/s     ',pver, 'A','Longwave heating rate',phys_decomp)
      call addfld ('CLDTOT  ','fraction',1,    'A','Vertically-integrated total cloud',phys_decomp)
      call addfld ('CLDLOW  ','fraction',1,    'A','Vertically-integrated low cloud',phys_decomp)
      call addfld ('CLDMED  ','fraction',1,    'A','Vertically-integrated mid-level cloud',phys_decomp)
      call addfld ('CLDHGH  ','fraction',1,    'A','Vertically-integrated high cloud',phys_decomp)
      call addfld ('ENGYCORR','W/m2    ',plev, 'A','Energy correction for over-all conservation',dyn_decomp)
!
! Put TS1 on history file since it will be different than TS in a coupled
! run
!
      do m=1,plevmx
         call addfld (tsnam(m),'K       ',1,    'A',tsnam(m)//' subsoil temperature',phys_decomp)
      end do

      call addfld ('TS      ','K       ',1,    'A','Surface temperature (radiative)',phys_decomp)
      call addfld ('TSMN    ','K       ',1,    'M','Minimum surface temperature over output period',phys_decomp)
      call addfld ('TSMX    ','K       ',1,    'X','Maximum surface temperature over output period',phys_decomp)
      call addfld ('SOLIN   ','W/m2    ',1,    'A','Solar insolation',phys_decomp)
      call addfld ('FU      ','m/s     ',plev, 'I','Zonal wind forcing term',dyn_decomp)
      call addfld ('FV      ','m/s     ',plev, 'I','Meridional wind forcing term',dyn_decomp)
      call addfld ('UTEND   ','m/s2    ',plev, 'A','U tendency',dyn_decomp)
      call addfld ('VTEND   ','m/s2    ',plev, 'A','V tendency',dyn_decomp)
      call addfld ('TTEND   ','K/s     ',plev, 'A','T tendency',dyn_decomp)
      call addfld ('LPSTEN  ','Pa/s    ',1,    'A','Surface pressure tendency',dyn_decomp)
      call addfld ('DTCOND  ','K/s     ',pver, 'A','T tendency - moist processes',phys_decomp)
      call addfld ('CMFDT   ','K/s     ',pver, 'A','T tendency - Hack convection',phys_decomp)
      call addfld ('CMFDQ   ','kg/kg/s ',pver, 'A','Q tendency - Hack convection',phys_decomp)
      call addfld ('ZMDT    ','K/s     ',pver, 'A','T tendency - Zhang-McFarlane moist convection',phys_decomp)
      call addfld ('ZMDQ    ','kg/kg/s ',pver, 'A','Q tendency - Zhang-McFarlane moist convection',phys_decomp)
      call addfld ('CMFDQR  ','kg/kg/s ',pver, 'A','Q tendency - shallow convection rainout',phys_decomp)
      call addfld ('QC      ','kg/kg/s ',pver, 'A','Q tendency - shallow convection LW export',phys_decomp)
      call addfld ('CMFMC   ','kg/m2/s ',pverp,'A','Moist convection mass flux',phys_decomp)
      call addfld ('CMFSL   ','W/m2    ',pverp,'A','Moist convection liquid water static energy flux',phys_decomp)
      call addfld ('CMFLQ   ','W/m2    ',pverp,'A','Moist convection total water flux',phys_decomp)
      call addfld ('CNVCLD  ','fraction',1,    'A','Vertically integrated convective cloud amount',phys_decomp)
      call addfld ('FSDS    ','W/m2    ',1,    'A','Downwelling solar flux at surface',phys_decomp)
      call addfld ('FSNIRTOA','W/m2    ',1,    'A','Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', &
            phys_decomp)
      call addfld ('FSNRTOAC','W/m2    ',1,    'A','Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', &
            phys_decomp)
      call addfld ('FSNRTOAS','W/m2    ',1,    'A','Net near-infrared flux (>= 0.7 microns) at top of atmosphere',phys_decomp)
      call addfld ('O3VMR   ','m3/m3   ',pver, 'A','Ozone volume mixing ratio',phys_decomp)
      call addfld ('CLDST   ','fraction',pver, 'A','Stratus cloud fraction',phys_decomp)
      call addfld ('CONCLD  ','fraction',pver, 'A','Convective cloud cover',phys_decomp)
!++pcw                                  
      call addfld ('RELHUM  ','percent ',pver, 'A','Relative humidity',phys_decomp)
!--pcw

      call addfld ('SULFBIO ','kg/kg   ',pver, 'A','Biogenic sulfate mass mixing ratio' ,phys_decomp)
      call addfld ('SULFANT ','kg/kg   ',pver, 'A','Anthropogenic sulfate mass mixing ratio' ,phys_decomp)
      call addfld ('SULFMMR ','kg/kg   ',pver, 'A','Sulfate mass mixing ratio' ,phys_decomp)
!++tls                                  
      call addfld ('TAUVIS  ','unitless',1,    'A','Total column aerosol extinction, vis band [aerosol optical depth]',phys_decomp)
      call addfld ('MSO4    ','gram/cm3',pver, 'A','Mass concentration of SO4',phys_decomp)
      call addfld ('LWC     ','kg/m3   ',pver, 'A','Liquid Water Content',phys_decomp)
      call addfld ('CLDFRQ  ','fraction',pver, 'A','Frequency of occurance of clouds (CLOUD > 0.01)',phys_decomp)
      call addfld ('WREL    ','um      ',pver, 'A','Weighted effective radius (by CLDFRQ)',phys_decomp)
      call addfld ('WLWC    ','kg/m3   ',pver, 'A','Weighted Liquid Water Content, prognostic (by CLDFRQ)',phys_decomp)
      call addfld ('TREFHTMN','K       ',1,    'M','Minimum reference height temperature over output period',phys_decomp)
      call addfld ('TREFHTMX','K       ',1,    'X','Maximum reference height temperature over output period',phys_decomp)
      call addfld ('TREFMNAV','K       ',1,    'A','Average of TREFHT daily minimum',phys_decomp)
      call addfld ('TREFMXAV','K       ',1,    'A','Average of TREFHT daily maximum',phys_decomp)
      call addfld ('LANDFRAC','fraction',1,    'A','Fraction of sfc area covered by land',phys_decomp)
      call addfld ('ICEFRAC ','fraction',1,    'A','Fraction of sfc area covered by sea-ice',phys_decomp)
      call addfld ('OCNFRAC ','fraction',1,    'A','Fraction of sfc area covered by ocean',phys_decomp)
      call addfld ('ASDIR',   '1',       1,    'A','albedo: shortwave, direct', phys_decomp)
      call addfld ('ASDIF',   '1',       1,    'A','albedo: shortwave, diffuse', phys_decomp)
      call addfld ('ALDIR',   '1',       1,    'A','albedo: longwave, direct', phys_decomp)
      call addfld ('ALDIF',   '1',       1,    'A','albedo: longwave, diffuse', phys_decomp)
      call addfld ('SST',     'K',       1,    'A','sea surface temperature', phys_decomp)
      call addfld ('UBOT    ','m/s     ',1,    'A','Lowest model level zonal wind',phys_decomp)
      call addfld ('VBOT    ','m/s     ',1,    'A','Lowest model level meridional wind',phys_decomp)
      call addfld ('QBOT    ','kg/kg   ',1,    'A','Lowest model level water vapor mixing ratio',phys_decomp)
      call addfld ('TBOT    ','K       ',1,    'A','Lowest model level temperature', phys_decomp)
      call addfld ('PBOT    ','Pa      ',1,    'A','Lowest model level pressure', phys_decomp)
      call addfld ('ZBOT    ','m       ',1,    'A','Lowest model level height', phys_decomp)
#if ( defined COUP_CSM )
      call addfld ('CPLRAINC','kg/m2/s ',1,    'A','Convective rainfall sent to coupler' ,phys_decomp)
      call addfld ('CPLRAINL','kg/m2/s ',1,    'A','Large-scale rainfall sent to coupler' ,phys_decomp)
      call addfld ('CPLSNOWC','kg/m2/s ',1,    'A','Convective snowfall sent to coupler' ,phys_decomp)
      call addfld ('CPLSNOWL','kg/m2/s ',1,    'A','Large-scale snowfall sent to coupler' ,phys_decomp)
      call addfld ('CPLPRCER','kg/m2/s ',1,    'A','Error in precipitation state (rain or snow) sent to coupler' ,phys_decomp)
#endif
      call addfld ('PRECCav ','m/s     ',1,    'A','Average large-scale precipitation',phys_decomp)
      call addfld ('PRECLav ','m/s     ',1,    'A','Average convective precipitation',phys_decomp)
!
! Variables output on specific pressure surfaces
!
      call addfld ('Z700    ','m       ',1,    'A','Geopotential Z at 700 mbar pressure surface',phys_decomp)
      call addfld ('Z500    ','m       ',1,    'A','Geopotential Z at 500 mbar pressure surface',phys_decomp)
      call addfld ('Z300    ','m       ',1,    'A','Geopotential Z at 300 mbar pressure surface',phys_decomp)
      call addfld ('Z050    ','m       ',1,    'A','Geopotential Z at 50 mbar pressure surface',phys_decomp)
      call addfld ('Q850    ','kg/kg   ',1,    'A','Specific Humidity at 850 mbar pressure surface',phys_decomp)
      call addfld ('Q200    ','kg/kg   ',1,    'A','Specific Humidity at 700 mbar pressure surface',phys_decomp)
      call addfld ('T850    ','K       ',1,    'A','Temperature at 850 mbar pressure surface',phys_decomp)
      call addfld ('T300    ','K       ',1,    'A','Temperature at 300 mbar pressure surface',phys_decomp)
      call addfld ('U850    ','m/s     ',1,    'A','Zonal wind at 850 mbar pressure surface',phys_decomp)
      call addfld ('U200    ','m/s     ',1,    'A','Zonal wind at 200 mbar pressure surface',phys_decomp)
      call addfld ('V850    ','m/s     ',1,    'A','Meridional wind at 850 mbar pressure surface',phys_decomp)
      call addfld ('V200    ','m/s     ',1,    'A','Meridional wind at 200 mbar pressure surface',phys_decomp)
      call addfld ('OMEGA850','Pa/s    ',1,    'A','Vertical velocity at 850 mbar pressure surface',phys_decomp)
      call addfld ('OMEGA500','Pa/s    ',1,    'A','Vertical velocity at 500 mbar pressure surface',phys_decomp)
!
! Fields needed such that all variables passed to flux coupler
! also appear on history tape
!
      call addfld ('SOLL    ','W/m2    ',1,    'A','Solar downward near infrared direct  to surface',phys_decomp)
      call addfld ('SOLS    ','W/m2    ',1,    'A','Solar downward visible direct  to surface',phys_decomp)
      call addfld ('SOLLD   ','W/m2    ',1,    'A','Solar downward near infrared diffuse to surface',phys_decomp)
      call addfld ('SOLSD   ','W/m2    ',1,    'A','Solar downward visible diffuse to surface',phys_decomp)
!
! Heating rate needed for d(theta)/dt computation
!
      call addfld ('HR      ','K/s     ',pver, 'A','Heating rate needed for d(theta)/dt computation',phys_decomp)
!
! Fields which go out active monthly average history files
!
      call addfld ('VT      ','K m/s   ',plev, 'A','Meridional heat transport',phys_decomp)
      call addfld ('VU      ','m2/s2   ',plev, 'A','Meridional flux of zonal momentum' ,phys_decomp)
      call addfld ('VV      ','m2/s2   ',plev, 'A','Meridional velocity squared' ,phys_decomp)
      call addfld ('WSPEED  ','m/s     ',plev, 'X','Horizontal total wind speed' ,phys_decomp)
      call addfld ('UU      ','m2/s2   ',plev, 'A','Zonal velocity squared' ,phys_decomp)
      call addfld ('ZZ      ','m2      ',pver ,'A','Eddy height variance' ,phys_decomp)
      call addfld ('TT      ','K2      ',plev ,'A','Eddy temperature variance' ,phys_decomp)
      call addfld ('OMEGAT  ','K Pa/s  ',pver ,'A','Vertical heat flux' ,phys_decomp)
      call addfld ('OMEGAU  ','K Pa/s  ',pver ,'A','Vertical flux of zonal momentum' ,phys_decomp)
      call addfld ('VZ      ','m2/s    ',plev, 'A','Meridional transport of geopotential energy',phys_decomp)
      call addfld ('VQ      ','m/skg/kg',pver, 'A','Meridional water transport',phys_decomp)
      call addfld ('Z3      ','m       ',pver, 'A','Geopotential Height (above sea level)',phys_decomp)
      call addfld ('MQ      ','kg/m2   ',pver, 'A','Water vapor mass in layer',phys_decomp)
      call addfld ('TMQ     ','kg/m2   ',1,    'A','Total (vertically integrated) precipitatable water',phys_decomp)
      call addfld ('PSL     ','Pa      ',1,    'A','Sea level pressure',phys_decomp)
!
! Fields which are only needed for recording attributes in master field list
!
      call addfld ('LANDMCOS','unitless',1,    'I', &
                   'Land ocean transition mask: ocean (0), continent (1), transition (0-1)',phys_decomp)
#if ( defined COUP_SOM )
      call addfld ('MELTB   ','m/s   ',1,'A','Sea ice basal melt rate',phys_decomp)
      call addfld ('MELTT   ','m/s   ',1,'A','Sea ice top surface melt rate',phys_decomp)
      call addfld ('MELTL   ','m/s   ',1,'A','Sea ice lateral surface melt rate',phys_decomp)
      call addfld ('GROWB   ','m/s   ',1,'A','Sea ice basal growth rate',phys_decomp)
      call addfld ('FRAZIL  ','m/s   ',1,'A','Sea ice frazil growth rate',phys_decomp)
      call addfld ('FLOOD   ','units?',1,'A','Description of FLOOD goes here',phys_decomp)
      call addfld ('FRZMLT  ','W/m2  ',1,'A','Potential to grow/melt sea ice',phys_decomp)
      call addfld ('QFLUX   ','W/m2  ',1,'A','Ocean mixed layer heat flux',phys_decomp)
      call addfld ('QFLUX_FT','W/m2  ',1,'A','M.L. heat flux after ft adjustment',phys_decomp)
      call addfld ('QFLUX_TH','W/m2  ',1,'A','M.L. heat flux after thickness factor adjustment',phys_decomp)
      call addfld ('QFLUX_A2','W/m2  ',1,'A','Adjusted ocean mixed layer heat flux 2',phys_decomp)
      call addfld ('FOCN    ','W/m2  ',1,'A','Melting heat flux from ice model',phys_decomp)
      call addfld ('EICEIN  ','W/m2  ',1,'A','Ice internal energy init/dtime',phys_decomp)
      call addfld ('EICEOUT ','W/m2  ',1,'A','Ice internal energy final/dtime',phys_decomp)
      call addfld ('F_ICE   ','W/m2  ',1,'A','net ice sfc flx over water/ice',phys_decomp)
      call addfld ('F_OCN   ','W/m2  ',1,'A','ice melt flx',phys_decomp)
      call addfld ('FRZMLTMX','W/m2  ',1,'A','ice
      call addfld ('DELTAICE','W/m2  ',1,'A','change in ice area',phys_decomp)
      call addfld ('IMBAL   ','W/m2  ',1,'A','local ice imbalance',phys_decomp)
      call addfld ('NRGERROR','W/m2  ',1,'A','local ice imbalance',phys_decomp)
      call addfld ('MLDANN  ','m     ',1,'I','mixed layer depth',phys_decomp)
      call addfld ('ONF     ','W/m2  ',1,'A','net ocn sfc flx over water/ice',phys_decomp)
      call addfld ('OIE     ','J/m2  ',1,'A','ocean internal energy',phys_decomp)
      call addfld ('OIERATE ','W/m2  ',1,'A','ocean internal energy change rate',phys_decomp)
      call addfld ('NRGICE  ','J/m2  ',1,'A','ice internal energy',phys_decomp)
      call addfld ('IIERATE ','W/m2  ',1,'A','ice internal energy change rate',phys_decomp)
#endif
      call addfld ('TSICE   ','K     ',1,'A','Ice temperature',phys_decomp)
!
! Fields to calculate surface energy budget for each component model
!
      call addfld ('FSNSLND ','W/m2  ',1,'A','FSNS over land',phys_decomp)
      call addfld ('FLNSLND ','W/m2  ',1,'A','FLNS over land',phys_decomp)
      call addfld ('LHFLXLND','W/m2  ',1,'A','LHFLX over land',phys_decomp)
      call addfld ('SHFLXLND','W/m2  ',1,'A','SHFLX over land',phys_decomp)

      call addfld ('FSNSOCN ','W/m2  ',1,'A','FSNS over open ocn',phys_decomp)
      call addfld ('FLNSOCN ','W/m2  ',1,'A','FLNS over open ocn',phys_decomp)
      call addfld ('LHFLXOCN','W/m2  ',1,'A','LHFLX over open ocn',phys_decomp)
      call addfld ('SHFLXOCN','W/m2  ',1,'A','SHFLX over open ocn',phys_decomp)

      call addfld ('FSNSICE ','W/m2  ',1,'A','FSNS over sea ice',phys_decomp)
      call addfld ('FLNSICE ','W/m2  ',1,'A','FLNS over sea ice',phys_decomp)
      call addfld ('LHFLXICE','W/m2  ',1,'A','LHFLX over sea ice',phys_decomp)
      call addfld ('SHFLXICE','W/m2  ',1,'A','SHFLX over sea ice',phys_decomp)

! Fields to calculate ocean/ice forcing for som

      call addfld ('FSNSOI ','W/m2  ',1,'A','FSNS over open ocn and ice',phys_decomp)
      call addfld ('FLNSOI ','W/m2  ',1,'A','FLNS over open ocn and ice',phys_decomp)
      call addfld ('LHFLXOI','W/m2  ',1,'A','LHFLX over open ocn and ice',phys_decomp)
      call addfld ('SHFLXOI','W/m2  ',1,'A','SHFLX over open ocn and ice',phys_decomp)

!
! Add NSTEP to history tape for debugging
!
      call addfld ('NSTEP   ','timestep',1,'A','Model timestep',phys_decomp)
!
! Print master field list
!
      if (masterproc) then
         write(6,*)'BLDFLD:nfmaster=',nfmaster
         write(6,*)' ******* MASTER FIELD LIST *******'

         do j=1,nfmaster
            write(6,9000)j,masterlist(j)%field%name,masterlist(j)%field%units
9000        format (i5,1x,a8,1x,a8)
         end do
      end if
!
!  Now that masterlist is defined and we are performing a restart run
!  (htapes_defined == .true.), construct primary and secondary hashing tables.
!
      if ( htapes_defined ) then
         call bld_outfld_hash_tbls()
         call bld_htapefld_indices()
      end if

      return
   end subroutine bldfld

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


   subroutine addfld (fname, units, numlev, avgflag, long_name, & 377,12
                      decomp_type, flag_xyfill, flag_isccplev)
!----------------------------------------------------------------------- 
! 
! Purpose: Add a field to the master field list
! 
! Method: Put input arguments of field name, units, number of levels, averaging flag, and 
!         long name into a type entry in the global master field list (masterlist).
! 
! Author: CCM Core Group
! 
!-----------------------------------------------------------------------
      use ppgrid,    only: begchunk, endchunk
      use rgrid,     only: nlon
      use phys_grid, only: get_ncols_p, physgrid_set
      use dycore, only: dycore_is
!
! Arguments
!
      character(len=*), intent(in) :: fname      ! field name--should be 8 characters
      character(len=*), intent(in) :: units      ! units of fname--should be 8 chars
      character(len=1), intent(in) :: avgflag    ! averaging flag
      character(len=*), intent(in) :: long_name  ! long name of field
      
      integer, intent(in) :: numlev              ! number of vertical levels (dimension and loop)
      integer, intent(in) :: decomp_type         ! decomposition type

      logical, intent(in), optional :: flag_xyfill ! non-applicable xy points flagged with fillvalue
      logical, intent(in), optional :: flag_isccplev ! levels are ISCCP levels not vertical
!
! Local workspace
!
      integer :: n     ! loop index
      integer :: c     ! chunk (physics ) or latitude (dynamics) index
      integer :: ncol  ! number of columns per chunk
!
! Ensure that required grid info is available now
!
      select case (decomp_type)
      case (phys_decomp)
         if (.not. physgrid_set) then
            call endrun ('ADDFLD: Attempt to add field '//fname//' to masterlist before physics grid set')
         end if
      case (dyn_decomp)
         if (.not. dyngrid_set) then
            call endrun ('ADDFLD: Attempt to add field '//fname//' to masterlist before dynamics grid set')
         end if
      end select
!
! Ensure that new field is not all blanks
!
      if (fname == ' ') then
         write(6,*)'ADDFLD: blank field name not allowed'
      end if
!
! Ensure that new field doesn't already exist
!
      do n=1,nfmaster
         if (masterlist(n)%field%name == fname) then
            call endrun ('ADDFLD:'//fname//' already on list')
         end if
      end do
!
! Add field to Master Field List arrays fieldn and iflds
!
      nfmaster = nfmaster + 1  ! Increase number of fields on Master F.L.
!
! Check number of fields against max for master list.
!
      if (nfmaster > pflds) then
         write(6,*)'ADDFLD: Too many fields for primary history file -- pflds,nfmaster=', &
                   pflds,nfmaster
         stop 999
      end if
      masterlist(nfmaster)%field%name        = fname
      masterlist(nfmaster)%field%long_name   = long_name
      masterlist(nfmaster)%field%units       = units
      masterlist(nfmaster)%field%numlev      = numlev
      masterlist(nfmaster)%field%decomp_type = decomp_type
!
! Whether to apply xy fillvalue: default is false
!
      if (present(flag_xyfill)) then
         masterlist(nfmaster)%field%flag_xyfill = flag_xyfill
      else
         masterlist(nfmaster)%field%flag_xyfill = .false.
      end if
!
! Whether level dimension is ISCCP (currently 49) or CAM
!
      if (present(flag_isccplev)) then
         masterlist(nfmaster)%field%flag_isccplev = flag_isccplev
      else
         masterlist(nfmaster)%field%flag_isccplev = .false.
      end if

      if (masterlist(nfmaster)%field%flag_isccplev .and. numlev /= npres*ntau) then
         write(6,*)'ADDFLD: Number of ISCCP levels must be ',npres*ntau, ' got ', numlev
         call endrun ()
      end if
!
! Dimension history info based on decomposition type (dynamics or physics)
!
      select case (decomp_type)
      case (phys_decomp)
         masterlist(nfmaster)%field%begdim3  = begchunk
         masterlist(nfmaster)%field%enddim3  = endchunk
         masterlist(nfmaster)%field%begver   = 1
         masterlist(nfmaster)%field%endver   = numlev
         allocate (masterlist(nfmaster)%field%colperdim3(begchunk:endchunk))
         do c=begchunk,endchunk
            ncol = get_ncols_p(c)
            masterlist(nfmaster)%field%colperdim3(c) = ncol
         end do
         masterlist(nfmaster)%field%coldimin   = pcols
      case (dyn_decomp)
         masterlist(nfmaster)%field%begdim3  = beglat
         masterlist(nfmaster)%field%enddim3  = endlat

         if ( dycore_is('LR') )then
# if ( defined STAGGERED )
            select case (numlev)
            case (1)
               masterlist(nfmaster)%field%begver   = 1
               masterlist(nfmaster)%field%endver   = 1
            case (plev)
               masterlist(nfmaster)%field%begver   = beglev
               masterlist(nfmaster)%field%endver   = endlev
            case (plevp)
               masterlist(nfmaster)%field%begver   = beglev
               masterlist(nfmaster)%field%endver   = endlevp
            case default
               write(6,*)'ADDFLD: invalid number of levels=', numlev
               call endrun ()
            end select
# endif
         else
            masterlist(nfmaster)%field%begver   = 1
            masterlist(nfmaster)%field%endver   = numlev
         endif

         allocate (masterlist(nfmaster)%field%colperdim3(beglat:endlat))
         do c=beglat,endlat
            masterlist(nfmaster)%field%colperdim3(c) = nlon(c)
         end do
         masterlist(nfmaster)%field%coldimin = plon
      case default
         write(6,*)'ADDFLD: unknown decomp_type=', decomp_type
         call endrun ()
      end select
!
! These 2 fields are used only in master field list, not runtime field list
!
      masterlist(nfmaster)%avgflag(:) = avgflag
      masterlist(nfmaster)%actflag(:) = .false.

      select case (avgflag)
      case ('A')
         masterlist(nfmaster)%time_op(:) = 'mean'
      case ('I')
         masterlist(nfmaster)%time_op(:) = ' '
      case ('X')
         masterlist(nfmaster)%time_op(:) = 'maximum'
      case ('M')
         masterlist(nfmaster)%time_op(:) = 'minimum'
      case default
         call endrun ('ADDFLD: unknown avgflag='//avgflag)
      end select
      
      return
   end subroutine addfld

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


   subroutine wrapup () 1,10
!-----------------------------------------------------------------------
!
! Purpose: 
! Close and dispose of history files.
! 
! Method: 
! This routine will close and dispose to Mass Store any full hist. files
! or any hist. file that has data on it when restart files are being 
! written.
! If a partially full history file was disposed (for restart 
! purposes), then wrapup will open that unit back up and position 
! it for appending new data. 
!
! Original version: CCM2
!
!-----------------------------------------------------------------------
      use shr_kind_mod, only: r8 => shr_kind_r8
      use pspect
      use ioFileMod
      use time_manager, only: get_nstep, get_curr_date, get_curr_time

#include <comctl.h>
#include <comlun.h>

!
! Local workspace
!
      integer :: nstep          ! current timestep number
      integer :: ncdate         ! current date in integer format [yyyymmdd]
      integer :: ncsec          ! time of day relative to current date [seconds]
      integer :: ndcur          ! days component of current time
      integer :: nscur          ! seconds component of current time
      integer :: yr, mon, day   ! year, month, day components of a date

      logical lfill   (ptapes)  ! Is history file ready to dispose?
      logical hdispose(ptapes)  ! Primary file disposed
      logical lremov            ! Rm file after dispose? 
      logical lhdisp            ! true => history file is disposed
      logical lhfill            ! true => history file is full

      integer t                 ! History file number

      real(r8) tday             ! Model day number for printout
!
! Externals
!
      logical, external :: rstwr ! whether it is time to write restarts
!-----------------------------------------------------------------------

      nstep = get_nstep()
      call get_curr_date(yr, mon, day, ncsec)
      ncdate = yr*10000 + mon*100 + day
      call get_curr_time(ndcur, nscur)
!
!-----------------------------------------------------------------------
! Dispose history files.
!-----------------------------------------------------------------------
!
! Begin loop over mtapes (the no. of declared history files - primary
! and auxiliary).  This loop disposes a history file to Mass Store
! when appropriate.
!
      do t=1,mtapes

         hdispose(t) = .false.
         lfill(t) = .false.
!
! Find out if file is full
!
         if (hstwr(t) .and. nfils(t) >= mfilt(t)) then
            lfill(t) = .true.
         endif
!
! Dispose history file to mass store if 
!    1) file is filled or 
!    2) this is the end of run and file has data on it or
!    3) restarts are being put out and history file has data on it
!
         if (lfill(t) .or. (nlend .and. nfils(t) >= 1) .or. (rstwr() .and. nfils(t) >= 1)) then
!
! Dispose history file
!
            hdispose(t) = .true.
!
! Remove the disk history file unless this is the end of the run, or
! if the history file is not full.
!
            lremov = .true.
            if (nlend.or..not.lfill(t)) lremov = .false.
!
! Is this the 0 timestep data of a monthly run?
! If so, just close primary unit do not dispose. 
!
            if (masterproc) then
               write(6,*)'WRAPUP: nf_close(',nfid(t),')=',nhfil(t)
               call wrap_close (nfid(t))

               if (nhtfrq(t) /= 0 .or. nstep > 0) then

! Dispose history file.  

                  call putfil (nhfil(t), cpath(t), mss_wpass, mss_irt, lremov)
! 
! Print information concerning model output.
! Model day number = iteration number of history file data * delta-t / (seconds per day)
! 
                  tday = ndcur + nscur/86400._r8
                  if (t==1) then
                     write(6,*)'   Primary history file'
                  else
                     write(6,*)'   Auxiliary history file number ', t-1
                  end if
                  write(6,9003)nstep,nfils(t),tday
                  write(6,9004)
!                      
! Auxilary files may have been closed and saved off without being full. 
! We must reopen the files and position them for more data.
! Must position auxiliary files if not full
!              
                  if (.not.nlend .and. .not.lfill(t)) then
                     call wrap_open (nhfil(t), NF_WRITE, nfid(t))
                  end if
               endif                 ! if 0 timestep of montly run****
            end if                   ! if (masterproc)
         end if                      ! if time dispose history fiels***   
      end do                         ! do mtapes
!
! Reset number of files on each history tape
!
      do t=1,mtapes
         lhfill = hstwr(t) .and. nfils(t) >= mfilt(t)
         lhdisp = lhfill .or. (nlend .and. nfils(t) >= 1) .or. &
                              (rstwr() .and. nfils(t) >= 1)
         if (lhfill.and.lhdisp) then
            nfils(t) = 0
         endif
      end do
      
      return
9003  format('    Output at NSTEP     = ',i10,/, &
             '    Number of time samples on this file = ',i10,/, &
             '    Model Day           = ',f10.2)
9004  format('---------------------------------------')
   end subroutine wrapup

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


   subroutine allocate_hbuf2d (hbuf, dimind, hbuf_prec) 1
!-----------------------------------------------------------------------
!
! Purpose: Allocate memory for the 2-D history buffer of the right precision.
!          Nullify the other buffer.
!
!-----------------------------------------------------------------------
      type (hbuffer_2d) :: hbuf                   ! history buffer
      type (dim_index_2d) :: dimind               ! 2-D dimension index
      integer, intent(in) :: hbuf_prec            ! precision for history buffer

      if (hbuf_prec == 8) then
         nullify  (hbuf%buf4)
         allocate (hbuf%buf8(dimind%beg1:dimind%end1, dimind%beg2:dimind%end2))
      else
         allocate (hbuf%buf4(dimind%beg1:dimind%end1, dimind%beg2:dimind%end2))
         nullify  (hbuf%buf8)
      end if

   end subroutine allocate_hbuf2d

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


   subroutine allocate_hbuf3d (hbuf, dimind, hbuf_prec) 1
!-----------------------------------------------------------------------
!
! Purpose: Allocate memory for the 3-D history buffer of the right precision
!          Nullify the other buffer.
!
!-----------------------------------------------------------------------
      type (hbuffer_3d) :: hbuf                   ! history buffer
      type (dim_index_3d) :: dimind               ! 3-D dimension index
      integer, intent(in) :: hbuf_prec            ! precision for history buffer

      if (hbuf_prec == 8) then
         nullify  (hbuf%buf4)
         allocate (hbuf%buf8(dimind%beg1:dimind%end1, dimind%beg2:dimind%end2, &
                             dimind%beg3:dimind%end3))
      else
         allocate (hbuf%buf4(dimind%beg1:dimind%end1, dimind%beg2:dimind%end2, &
                             dimind%beg3:dimind%end3))
         nullify  (hbuf%buf8)
      end if

   end subroutine allocate_hbuf3d

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


   subroutine deallocate_hbuf2d (hbuf) 1
!-----------------------------------------------------------------------
!
! Purpose: Deallocate memory for 2-D history buffer
!
!-----------------------------------------------------------------------
      type (hbuffer_2d) :: hbuf                   ! history buffer

      if (associated(hbuf%buf8)) then
         deallocate (hbuf%buf8)
      else if (associated(hbuf%buf4)) then
         deallocate (hbuf%buf4)
      end if

   end subroutine deallocate_hbuf2d

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


   subroutine deallocate_hbuf3d (hbuf) 1
!-----------------------------------------------------------------------
!
! Purpose: Deallocate memory for 3-D history buffer
!
!-----------------------------------------------------------------------
      type (hbuffer_3d) :: hbuf                   ! history buffer

      if (associated(hbuf%buf8)) then
         deallocate (hbuf%buf8)
      else if (associated(hbuf%buf4)) then
         deallocate (hbuf%buf4)
      end if

   end subroutine deallocate_hbuf3d

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


   subroutine nullify_hbuf2d (hbuf) 1
!-----------------------------------------------------------------------
!
! Purpose: Nullify 2-D history buffer pointers
!
!-----------------------------------------------------------------------
      type (hbuffer_2d) :: hbuf                   ! history buffer

      nullify (hbuf%buf8, hbuf%buf4)

   end subroutine nullify_hbuf2d

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


   subroutine nullify_hbuf3d (hbuf) 1
!-----------------------------------------------------------------------
!
! Purpose: Nullify 3-D history buffer pointers
!
!-----------------------------------------------------------------------
      type (hbuffer_3d) :: hbuf                   ! history buffer

      nullify (hbuf%buf8, hbuf%buf4)

   end subroutine nullify_hbuf3d

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


   subroutine read_hbuf (hbuf, luhrest, ioerr) 2
!-----------------------------------------------------------------------
!
! Purpose: Read history information using the appropriate buffer
!
!-----------------------------------------------------------------------
      type (hbuffer_3d), intent(inout) :: hbuf    ! history buffer
      integer, intent(in ) :: luhrest
      integer, intent(out) :: ioerr

      if (associated(hbuf%buf8)) then
         read (luhrest,iostat=ioerr) hbuf%buf8
      else if (associated(hbuf%buf4)) then
         read (luhrest,iostat=ioerr) hbuf%buf4
      end if

   end subroutine read_hbuf

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


   subroutine write_hbuf (hbuf, luhrest, ioerr) 3
!-----------------------------------------------------------------------
!
! Purpose: Write (unformatted) history information to file unit luhrest
!
!-----------------------------------------------------------------------
      type (hbuffer_3d), intent(in) :: hbuf       ! history buffer
      integer, intent(in ) :: luhrest             ! file unit number
      integer, intent(out) :: ioerr               ! I/O error status

      if (associated(hbuf%buf8)) then
         write (luhrest,iostat=ioerr) hbuf%buf8
      else if (associated(hbuf%buf4)) then
         write (luhrest,iostat=ioerr) hbuf%buf4
      end if

   end subroutine write_hbuf

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


   subroutine hbuf_assigned_to_hbuf (hbuf1, hbuf2) 1
!-----------------------------------------------------------------------
!
! Purpose: Set hbuf1 to hbuf2 (copy the contents).
!
!-----------------------------------------------------------------------
      type (hbuffer_3d), intent(inout) :: hbuf1   ! history buffer
      type (hbuffer_3d), intent(in   ) :: hbuf2   ! history buffer

      if (associated(hbuf2%buf8)) then
         hbuf1%buf8 = hbuf2%buf8
      else if (associated(hbuf2%buf4)) then
         hbuf1%buf4 = hbuf2%buf4
      end if

   end subroutine hbuf_assigned_to_hbuf

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


   subroutine hbuf_assigned_to_real8 (hbuf, scalar) 1
!-----------------------------------------------------------------------
!
! Purpose: Set appropriate history buffer to the value, scalar
!
!-----------------------------------------------------------------------
      type (hbuffer_3d), intent(inout) :: hbuf    ! history buffer
      real(r8), intent(in) :: scalar              ! scalar real

      if (associated(hbuf%buf8)) then
         hbuf%buf8 = scalar
      else if (associated(hbuf%buf4)) then
         hbuf%buf4 = scalar
      end if

   end subroutine hbuf_assigned_to_real8

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


   subroutine set_hbuf_section_to_val (hbuf, dimind, scalar) 1
!-----------------------------------------------------------------------
!
! Purpose: Set section of appropriate history buffer to scalar
!
!-----------------------------------------------------------------------
      type (hbuffer_3d), intent(inout) :: hbuf    ! history buffer
      type (dim_index_3d), intent(in ) :: dimind  ! 3-D dimension index
      real(r8), intent(in) :: scalar              ! scalar real

      if (associated(hbuf%buf8)) then
         hbuf%buf8(dimind%beg1:dimind%end1, dimind%beg2:dimind%end2, &
                   dimind%beg3:dimind%end3) = scalar
      else if (associated(hbuf%buf4)) then
         hbuf%buf4(dimind%beg1:dimind%end1, dimind%beg2:dimind%end2, &
                   dimind%beg3:dimind%end3) = scalar
      end if

   end subroutine set_hbuf_section_to_val

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


   subroutine assoc_hbuf2d_with_hbuf3d (hbuf2d, hbuf3d, c) 2
!-----------------------------------------------------------------------
!
! Purpose: Associate 2-D hbuf2d with 3-D hbuf3d of column c.
!          Nullify the other buffer of hbuf2d.
!
!-----------------------------------------------------------------------
      type (hbuffer_2d), intent(out) :: hbuf2d    ! 2-D history buffer
      type (hbuffer_3d), intent(in ) :: hbuf3d    ! 3-D history buffer
      integer, intent(in) :: c                    ! chunk (or lat) index

      if (associated(hbuf3d%buf8)) then
         hbuf2d%buf8 => hbuf3d%buf8(:,:,c)
         nullify (hbuf2d%buf4)
      else if (associated(hbuf3d%buf4)) then
         hbuf2d%buf4 => hbuf3d%buf4(:,:,c)
         nullify (hbuf2d%buf8)
      end if

   end subroutine assoc_hbuf2d_with_hbuf3d

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


   subroutine assoc_hbuf_with_nothing (hbuf, hbuf_prec) 3
!-----------------------------------------------------------------------
!
! Purpose: Associate the appropriate 3-D hbuf pointer with nothing.
!          Nullify the other pointer.
!
!-----------------------------------------------------------------------
      type (hbuffer_3d), intent(out) :: hbuf      ! 3-D history buffer
      integer, intent(in) :: hbuf_prec            ! hbuffer precision

      if (hbuf_prec == 8) then
         hbuf%buf8 => nothing_r8
         nullify (hbuf%buf4)
      else if (hbuf_prec == 4) then
         hbuf%buf4 => nothing_r4
         nullify (hbuf%buf8)
      end if
      return
   end subroutine assoc_hbuf_with_nothing

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


   subroutine xzy_to_xyz (xyzbuf, xzybuf, dimind) 1,1
!-----------------------------------------------------------------------
!
! Purpose: Set the appropriate buffer of xyzbuf to the transpose of
!          xzybuf with respect to the last two indices.
!
!-----------------------------------------------------------------------
   use rgrid,     only: nlon

   type (hbuffer_3d), intent(inout) :: xyzbuf      ! data in COORDS order
   type (hbuffer_3d), intent(in   ) :: xzybuf      ! data in xzy order
   type (dim_index_3d), intent(in ) :: dimind      ! dim index (xzybuf order)
   integer j, k                                    ! loop indices
   integer endi                                    ! ending longitude index

   if (associated(xzybuf%buf8)) then

      do k=dimind%beg2,dimind%end2
         do j=dimind%beg3,dimind%end3
            endi = nlon(j)
            xyzbuf%buf8(dimind%beg1:endi,j,k) = xzybuf%buf8(dimind%beg1:endi,k,j)
            if (endi < dimind%end1) then
               xyzbuf%buf8(endi+1:,j,k) = fillvalue
            end if
         end do
      end do

   else if (associated(xzybuf%buf4)) then

      do k=dimind%beg2,dimind%end2
         do j=dimind%beg3,dimind%end3
            endi = nlon(j)
            xyzbuf%buf4(dimind%beg1:endi,j,k) = xzybuf%buf4(dimind%beg1:endi,k,j)
            if (endi < dimind%end1) then
               xyzbuf%buf4(endi+1:,j,k) = fillvalue
            end if
         end do
      end do

   end if
   return
   end subroutine xzy_to_xyz

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


   subroutine fill_unset (buf4, buf8, dimind) 1,1
!-----------------------------------------------------------------------
!
! Purpose: For max/min calcs which initialized hbuf to +-huge, set those values to fillvalue
!
!-----------------------------------------------------------------------
   use rgrid,     only: nlon

   real(r4), pointer :: buf4(:,:,:)
   real(r8), pointer :: buf8(:,:,:)
   type (dim_index_3d), intent(in ) :: dimind      ! dim index (xzybuf order)
   integer j, k                                    ! loop indices
   integer endi                                    ! ending longitude index

   if (associated(buf8)) then

      do k=dimind%beg2,dimind%end2
         do j=dimind%beg3,dimind%end3
            endi = nlon(j)
            do i=dimind%beg1,endi
               if (buf8(i,j,k) == -huge(buf8) .or. buf8(i,j,k) == +huge(buf8)) then
                  buf8(i,j,k) = fillvalue
               end if
            end do
         end do
      end do

   else if (associated(buf4)) then

      do k=dimind%beg2,dimind%end2
         do j=dimind%beg3,dimind%end3
            endi = nlon(j)
            do i=dimind%beg1,endi
               if (buf4(i,j,k) == -huge(buf4) .or. buf4(i,j,k) == +huge(buf4)) then
                  buf4(i,j,k) = fillvalue
               end if
            end do
         end do
      end do

   end if

   return
   end subroutine fill_unset

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


   subroutine hbuf_accum_inst (buf4, buf8, field, nacs, dimind, idim, flag_xyfill) 1,1
!-----------------------------------------------------------------------
!
! Purpose: Accumulate instantaneous values of field in 2-D hbuf.
!          Set accumulation counter to 1.
!
!-----------------------------------------------------------------------
      real(r4), pointer :: buf4(:,:)    ! 2-D history buffer
      real(r8), pointer :: buf8(:,:)    ! 2-D history buffer
      integer, pointer                 :: nacs(:) ! accumulation counter
      integer, intent(in)              :: idim    ! Longitude dimension of field array
      logical, intent(in)              :: flag_xyfill ! non-applicable xy points flagged with fillvalue
      real(r8),          intent(in ) :: field(idim,*)   ! real*8 array
      type (dim_index_2d), intent(in ) :: dimind  ! 2-D dimension index
!
! Local indices
!
      integer :: ib, ie    ! beginning and ending indices of first dimension
      integer :: jb, je    ! beginning and ending indices of second dimension
      integer :: ieu, jeu  ! number of elements in each dimension
      integer :: i, k      ! loop indices

      logical :: bad       ! flag indicates input field fillvalues not applied consistently
                           ! with vertical level

      ib = dimind%beg1
      ie = dimind%end1
      jb = dimind%beg2
      je = dimind%end2

      ieu = ie-ib+1
      jeu = je-jb+1

      if (associated(buf8)) then
         do k=1,jeu
            do i=1,ieu
               buf8(i,k) = field(i,k)
            end do
         end do
      else if (associated(buf4)) then
         do k=1,jeu
            do i=1,ieu
               buf4(i,k) = field(i,k)
            end do
         end do
      end if

      if (flag_xyfill) then
         do i=1,ieu
            if (field(i,1) == fillvalue) then
               nacs(i) = 0
            else
               nacs(i) = 1
            end if
            call check_accum (field, idim, ieu, jeu)
         end do
      else
         do i=1,ieu
            nacs(i) = 1
         end do
      end if

      return
   end subroutine hbuf_accum_inst

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


   subroutine check_accum (field, idim, ieu, jeu) 4,1
      integer, intent(in)  :: idim
      real(r8), intent(in) :: field(idim,*)   ! real*8 array
      integer, intent(in)  :: ieu,jeu         ! loop ranges

      logical :: bad
      integer :: i,k
!
! For multilevel fields ensure that all levels have fillvalue applied consistently
!
      bad = .false.
      do k=2,jeu
         do i=1,ieu
            if (field(i,1) == fillvalue .and. field(i,k) /= fillvalue .or. &
                field(i,1) /= fillvalue .and. field(i,k) == fillvalue) then
               bad = .true.
            end if
         end do
      end do

      if (bad) then
         call endrun ('CHECK_ACCUM: inconsistent level values')
      end if

      return
   end subroutine check_accum

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


   subroutine hbuf_accum_add (buf4, buf8, field, nacs, dimind, idim, flag_xyfill) 1,1
!-----------------------------------------------------------------------
!
! Purpose: Add the values of field to 2-D hbuf.
!          Increment accumulation counter by 1.
!
!-----------------------------------------------------------------------
      real(r4), pointer :: buf4(:,:)    ! 2-D history buffer
      real(r8), pointer :: buf8(:,:)    ! 2-D history buffer
      integer, pointer                 :: nacs(:) ! accumulation counter
      integer, intent(in) :: idim           ! Longitude dimension of field array
      logical, intent(in)              :: flag_xyfill ! non-applicable xy points flagged with fillvalue
      real(r8),          intent(in ) :: field(idim,*)   ! real*8 array
      type (dim_index_2d), intent(in ) :: dimind  ! 2-D dimension index
!
! Local indices
!
      integer :: ib, ie    ! beginning and ending indices of first dimension
      integer :: jb, je    ! beginning and ending indices of second dimension
      integer :: ieu, jeu  ! number of elements in each dimension
      integer :: i,k       ! indices

      ib = dimind%beg1
      ie = dimind%end1
      jb = dimind%beg2
      je = dimind%end2

      ieu = ie-ib+1
      jeu = je-jb+1

      if (flag_xyfill) then
         if (associated(buf8)) then
            do k=1,jeu
               do i=1,ieu
                  if (field(i,k) /= fillvalue) then
                     buf8(i,k) = buf8(i,k) + field(i,k)
                  end if
               end do
            end do
         else if (associated(buf4)) then
            do k=1,jeu
               do i=1,ieu
                  if (field(i,k) /= fillvalue) then
                     buf4(i,k) = buf4(i,k) + field(i,k)
                  end if
               end do
            end do
         end if
!
! Ensure input field has fillvalue defined invariant in the z-direction, then increment nacs
!
         call check_accum (field, idim, ieu, jeu)
         do i=1,ieu
            if (field(i,1) /= fillvalue) then
               nacs(i) = nacs(i) + 1
            end if
         end do
      else
         if (associated(buf8)) then
            do k=1,jeu
               do i=1,ieu
                  buf8(i,k) = buf8(i,k) + field(i,k)
               end do
            end do
         else if (associated(buf4)) then
            do k=1,jeu
               do i=1,ieu
                  buf4(i,k) = buf4(i,k) + field(i,k)
               end do
            end do
         end if
         do i=1,ieu
            nacs(i) = nacs(i) + 1
         end do
      end if

      return
   end subroutine hbuf_accum_add

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


   subroutine hbuf_accum_max (buf4, buf8, field, nacs, dimind, idim, flag_xyfill) 1,1
!-----------------------------------------------------------------------
!
! Purpose: Accumulate the maximum values of field in 2-D hbuf
!          Set accumulation counter to 1.
!
!-----------------------------------------------------------------------
      real(r4), pointer :: buf4(:,:)    ! 2-D history buffer
      real(r8), pointer :: buf8(:,:)    ! 2-D history buffer
      integer, pointer                 :: nacs(:) ! accumulation counter
      integer, intent(in) :: idim           ! Longitude dimension of field array
      logical, intent(in)              :: flag_xyfill ! non-applicable xy points flagged with fillvalue
      real(r8),          intent(in ) :: field(idim,*)   ! real*8 array
      type (dim_index_2d), intent(in ) :: dimind  ! 2-D dimension index
!
! Local indices
!
      integer :: ib, ie    ! beginning and ending indices of first dimension
      integer :: jb, je    ! beginning and ending indices of second dimension
      integer :: ieu, jeu  ! number of elements in each dimension
      integer :: i, k

      ib = dimind%beg1
      ie = dimind%end1
      jb = dimind%beg2
      je = dimind%end2

      ieu = ie-ib+1
      jeu = je-jb+1

      if (associated(buf8)) then

         if (flag_xyfill) then
            do k=1,jeu
               do i=1,ieu
                  if (nacs(i) == 0) then
                     buf8(i,k) = -huge (buf8)
                  end if
                  if (field(i,k) > buf8(i,k) .and. field(i,k) /= fillvalue) then
                     buf8(i,k) = field(i,k)
                  end if
               end do
            end do
         else
            do k=1,jeu
               do i=1,ieu
                  if (nacs(i) == 0) then
                     buf8(i,k) = -huge (buf8)
                  end if
                  if (field(i,k) > buf8(i,k)) then
                     buf8(i,k) = field(i,k)
                  end if
               end do
            end do
         end if

      else if (associated(buf4)) then

         if (flag_xyfill) then
            do k=1,jeu
               do i=1,ieu
                  if (nacs(i) == 0) then
                     buf4(i,k) = -huge (buf4)
                  end if
                  if (field(i,k) > buf4(i,k) .and. field(i,k) /= fillvalue) then
                     buf4(i,k) = field(i,k)
                  end if
               end do
            end do
         else
            do k=1,jeu
               do i=1,ieu
                  if (nacs(i) == 0) then
                     buf4(i,k) = -huge (buf4)
                  end if
                  if (field(i,k) > buf4(i,k)) then
                     buf4(i,k) = field(i,k)
                  end if
               end do
            end do
         end if
      end if

      if (flag_xyfill) then
         call check_accum (field, idim, ieu, jeu)
         do i=1,ieu
            if (field(i,1) /= fillvalue) then
               nacs(i) = 1
            end if
         end do
      else
         do i=1,ieu
            nacs(i) = 1
         end do
      end if

      return
   end subroutine hbuf_accum_max

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


   subroutine hbuf_accum_min (buf4, buf8, field, nacs, dimind, idim, flag_xyfill) 1,1
!-----------------------------------------------------------------------
!
! Purpose: Accumulate the minimum values of field in 2-D hbuf
!          Set accumulation counter to 1.
!
!-----------------------------------------------------------------------
      real(r4), pointer :: buf4(:,:)    ! 2-D history buffer
      real(r8), pointer :: buf8(:,:)    ! 2-D history buffer
      integer, pointer                 :: nacs(:) ! accumulation counter
      integer, intent(in) :: idim           ! Longitude dimension of field array
      logical, intent(in)              :: flag_xyfill ! non-applicable xy points flagged with fillvalue
      real(r8),          intent(in ) :: field(idim,*)   ! real*8 array
      type (dim_index_2d), intent(in ) :: dimind  ! 2-D dimension index
!
! Local indices
!
      integer :: ib, ie    ! beginning and ending indices of first dimension
      integer :: jb, je    ! beginning and ending indices of second dimension
      integer :: ieu, jeu  ! number of elements in each dimension
      integer :: i, k

      ib = dimind%beg1
      ie = dimind%end1
      jb = dimind%beg2
      je = dimind%end2

      ieu = ie-ib+1
      jeu = je-jb+1

      if (associated(buf8)) then

         if (flag_xyfill) then
            do k=1,jeu
               do i=1,ieu
                  if (nacs(i) == 0) then
                     buf8(i,k) = +huge (buf8)
                  end if
                  if (field(i,k) < buf8(i,k) .and. field(i,k) /= fillvalue) then
                     buf8(i,k) = field(i,k)
                  end if
               end do
            end do
         else
            do k=1,jeu
               do i=1,ieu
                  if (nacs(i) == 0) then
                     buf8(i,k) = +huge (buf8)
                  end if
                  if (field(i,k) < buf8(i,k)) then
                     buf8(i,k) = field(i,k)
                  end if
               end do
            end do
         end if

      else if (associated(buf4)) then

         if (flag_xyfill) then
            do k=1,jeu
               do i=1,ieu
                  if (nacs(i) == 0) then
                     buf4(i,k) = +huge (buf4)
                  end if
                  if (field(i,k) < buf4(i,k) .and. field(i,k) /= fillvalue) then
                     buf4(i,k) = field(i,k)
                  end if
               end do
            end do
         else
            do k=1,jeu
               do i=1,ieu
                  if (nacs(i) == 0) then
                     buf4(i,k) = +huge (buf4)
                  end if
                  if (field(i,k) < buf4(i,k)) then
                     buf4(i,k) = field(i,k)
                  end if
               end do
            end do
         end if
      end if

      if (flag_xyfill) then
         call check_accum (field, idim, ieu, jeu)
         do i=1,ieu
            if (field(i,1) /= fillvalue) then
               nacs(i) = 1
            end if
         end do
      else
         do i=1,ieu
            nacs(i) = 1
         end do
      end if

      return
   end subroutine hbuf_accum_min

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


   subroutine readin_hbuf (iu, hbuf, numperlat) 2,3
!-----------------------------------------------------------------------
!
! Purpose: Read history restart binary file 
!
!-----------------------------------------------------------------------
      use binary_io
      integer, intent(in) :: iu                   ! Logical unit
      integer, intent(in) :: numperlat            ! Length of arr
      type (hbuffer_3d), intent(inout) :: hbuf    ! hbuf to be read

      if (associated(hbuf%buf8)) then
         call readin_r8 (iu, hbuf%buf8, numperlat)
      else if (associated(hbuf%buf4)) then
         call readin_r4 (iu, hbuf%buf4, numperlat)
      end if

   end subroutine readin_hbuf

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

#ifdef SPMD

   subroutine mpigatherv_hbuf (hbuf_send, numsend, mpireal1, hbuf_recv, & 2,2
                               numrecv, displs, mpireal2, root, comm)
!-----------------------------------------------------------------------
!
! Purpose: Collects history buffer from each thread on masterproc
!
!-----------------------------------------------------------------------
      type (hbuffer_3d), intent(in   ) :: hbuf_send  ! send buffer
      type (hbuffer_3d), intent(inout) :: hbuf_recv  ! receive buffer
      integer :: numsend                ! number of items to be sent
      integer :: mpireal1               ! MPI real data type for hbuf_send
      integer :: mpireal2               ! MPI real data type for hbuf_recv
      integer :: numrecv(*)             ! number of items to be received
      integer :: displs(*)              ! displacement array
      integer, intent(in) :: root
      integer, intent(in) :: comm
 
      if (associated(hbuf_send%buf8)) then
         call mpigatherv (hbuf_send%buf8, numsend, mpireal1, hbuf_recv%buf8, numrecv, &
                          displs, mpireal2, root, comm)
      else if (associated(hbuf_send%buf4)) then
         call mpigatherv (hbuf_send%buf4, numsend, mpireal1, hbuf_recv%buf4, numrecv, &
                          displs, mpireal2, root, comm)
      end if

   end subroutine mpigatherv_hbuf
#endif

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


   subroutine wrap_put_vara_hbuf (nfid, varid, start, count, hbuf) 2,2
!-----------------------------------------------------------------------
!
! Purpose: Output the given portion of the real array.
!
!-----------------------------------------------------------------------
      type (hbuffer_3d), intent(in) :: hbuf    ! buffer to write
      integer, intent(in):: nfid               ! file ID
      integer, intent(in):: varid              ! variable ID
      integer, intent(in):: start(*)           ! array of starting field indices
      integer, intent(in):: count(*)           ! array of count values

      if (associated(hbuf%buf8)) then
         call wrap_put_vara_realx (nfid, varid, start, count, hbuf%buf8)
      else if (associated(hbuf%buf4)) then
         call wrap_put_vara_real  (nfid, varid, start, count, hbuf%buf4)
      end if

   end subroutine wrap_put_vara_hbuf

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


   subroutine gather_chunk_to_field_hbuf (fdim,mdim,ldim, & 3,3
                                          nlond,localchunks,globalfield)
!-----------------------------------------------------------------------
!
! Purpose: Reconstruct longitude/latitude field
!          from decomposed chunk data structure
!
!-----------------------------------------------------------------------
      use phys_grid, only: gather_chunk_to_field, gather_chunk_to_field4

      type (hbuffer_3d), intent(in   ) :: localchunks
      type (hbuffer_3d), intent(inout) :: globalfield
      integer, intent(in) :: fdim      ! declared length of first dimension
      integer, intent(in) :: mdim      ! declared length of middle dimension
      integer, intent(in) :: ldim      ! declared length of last dimension
      integer, intent(in) :: nlond     ! declared number of longitudes

      if (associated(localchunks%buf8)) then
         call gather_chunk_to_field (fdim,mdim,ldim, &
                                 nlond,localchunks%buf8,globalfield%buf8)
      else if (associated(localchunks%buf4)) then
         call gather_chunk_to_field4(fdim,mdim,ldim, &
                                 nlond,localchunks%buf4,globalfield%buf4)
      end if

   end subroutine gather_chunk_to_field_hbuf

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


   subroutine scatter_field_to_chunk_hbuf (fdim,mdim,ldim, & 2,3
                                      nlond,globalfield,localchunks)
!-----------------------------------------------------------------------
!
! Purpose: Reconstruct longitude/latitude field
!          from decomposed chunk data structure
!
!-----------------------------------------------------------------------
      use phys_grid, only: scatter_field_to_chunk, scatter_field_to_chunk4

      type (hbuffer_3d), intent(inout) :: localchunks
      type (hbuffer_3d), intent(in   ) :: globalfield
      integer, intent(in) :: fdim      ! declared length of first dimension
      integer, intent(in) :: mdim      ! declared length of middle dimension
      integer, intent(in) :: ldim      ! declared length of last dimension
      integer, intent(in) :: nlond     ! declared number of longitudes

      if (associated(globalfield%buf8)) then
         call scatter_field_to_chunk (fdim,mdim,ldim, &
                                   nlond,globalfield%buf8,localchunks%buf8)
      else if (associated(globalfield%buf4)) then
         call scatter_field_to_chunk4(fdim,mdim,ldim, &
                                   nlond,globalfield%buf4,localchunks%buf4)
      end if

   end subroutine scatter_field_to_chunk_hbuf
#if ( defined BFB_CAM_SCAM_IOP || defined SCAM )

  subroutine initialize_iop_history() 2,24
!
! !DESCRIPTION: 
! !USES:
    use iop
    use ppgrid, only: pver, pverp
! !ARGUMENTS:
    implicit none
!
! !CALLED FROM:
!
! !REVISION HISTORY:
!
!EOP
!
! !LOCAL VARIABLES:
    integer m
!-----------------------------------------------------------------------

  call addfld ('q','kg/kg   ',plev, 'A','Q for scam',dyn_decomp)
  call addfld ('u','m/s     ',plev, 'A','U for scam',dyn_decomp)
  call addfld ('v','m/s     ',plev, 'A','V for scam',dyn_decomp)
  call addfld ('t','K       ',plev, 'A','Temperature for scam',dyn_decomp)
  call addfld ('Tg','K      ',1,    'A','Surface temperature (radiative) for scam',phys_decomp)
  call addfld ('Ps','Pa      ',1, 'A','Ps for scam',dyn_decomp)
  call addfld ('divq3d','kg/kg   ',plev, 'A','Dynamics Residual for Q',dyn_decomp)
  call addfld ('dcldliq','kg/kg   ',plev, 'A','Dynamics Residual for liq cld',dyn_decomp)
  call addfld ('dcldice','kg/kg   ',plev, 'A','Dynamics Residual for ice cld',dyn_decomp)
  call addfld ('divT3d','K       ',plev, 'A','Dynamics Residual for T',dyn_decomp)
  call addfld ('fixmas','percent',1, 'A','Mass fixer',dyn_decomp)
  call addfld ('beta','percent  ',1, 'A','Mass fixer',dyn_decomp)
  do m=1,pcnst
    call addfld (alphanam(m),'kg/kg   ',1, 'A',trim(cnst_name(m))//' alpha constituent fixer',dyn_decomp)
    call addfld (dqfxnam(m),'kg/kg   ',plev, 'A',trim(cnst_name(m))//' dqfx3 fixer',dyn_decomp)
  end do
  call addfld ('CLAT ','none      ',1,    'A','cos lat for bfb testing', dyn_decomp)
  call addfld ('shflx ','W/m2    ',1,    'A','Surface sensible heat flux for scam',phys_decomp)
  call addfld ('lhflx   ','W/m2    ',1,    'A','Surface latent heat flux for scam',phys_decomp)
  call addfld ('trefht  ','K       ',1,    'A','Reference height temperature',phys_decomp)
  call addfld ('Tsair  ','K       ',1,    'A','Reference height temperature for scam',phys_decomp)
  call addfld ('phis   ','m2/s2   ',1,    'I','Surface geopotential for scam',phys_decomp)
  call addfld ('Prec   ','m/s     ',1,    'A','Total (convective and large-scale) precipitation rate for scam',phys_decomp)
  call addfld ('omega   ','Pa/s    ',pver, 'A','Vertical velocity (pressure)',phys_decomp)

  end subroutine initialize_iop_history
#endif
!#######################################################################


   integer function gen_hash_key(string) 3
!-----------------------------------------------------------------------
!
! Purpose: Generate a hash key on the interval [0 .. tbl_hash_pri_sz-1]
!          given a character string.
!
! Algorithm is a variant of perl's internal hashing function.
!
!-----------------------------------------------------------------------

   implicit none
!
!  Arguments:
!
   character(len=*), intent(in) :: string
!
!  Local.
!
   integer :: hash
   integer :: i

   hash = gen_hash_key_offset

   if ( len(string) /= 8 ) then
!
!     Process arbitrary string length.
!
      do i = 1, len(string)
         hash = ieor(hash , (ichar(string(i:i)) * tbl_gen_hash_key(iand(i,16-1))))
      end do
   else
!
!     Special case string length = 8.
!
      hash = ieor(hash , ichar(string(1:1)) * 61)
      hash = ieor(hash , ichar(string(2:2)) * 59)
      hash = ieor(hash , ichar(string(3:3)) * 53)
      hash = ieor(hash , ichar(string(4:4)) * 47)
      hash = ieor(hash , ichar(string(5:5)) * 43)
      hash = ieor(hash , ichar(string(6:6)) * 41)
      hash = ieor(hash , ichar(string(7:7)) * 37)
      hash = ieor(hash , ichar(string(8:8)) * 31)
   end if

!  gen_hash_key = iand(hash, ishft(1, tbl_hash_pri_sz_lg2)-1)
   gen_hash_key = iand(hash, (2**tbl_hash_pri_sz_lg2)-1)

   return
   end function gen_hash_key

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


   integer function get_masterlist_indx(fldname) 2,3
!-----------------------------------------------------------------------
!
! Purpose: Return the the index of the field's name on the master file list.
!
!          If the field is not found on the masterlist, return -1.
!
!-----------------------------------------------------------------------

   implicit none
!
!  Arguments:
!
   character(len=*), intent(in) :: fldname
!
!  Local.
!
   integer :: hash_key
   integer :: ff
   integer :: ii
   integer :: io   ! Index of overflow chain in overlow table
   integer :: in   ! Number of entries on overflow chain

   hash_key = gen_hash_key(fldname)
   ff = tbl_hash_pri(hash_key)
   if ( ff < 0 ) then
      io = abs(ff)
      in = tbl_hash_oflow(io)
      do ii = 1, in
         ff = tbl_hash_oflow(io+ii)
         if ( masterlist(ff)%field%name == fldname ) exit
      end do
   end if

   if (ff == 0) then
      ! fldname generated a hash key that doesn't have an entry in tbl_hash_pri.
      ! This means that fldname isn't in the masterlist
      call endrun ('GET_MASTERLIST_INDX: attemping to output field '//fldname//' not on master list')
   end if

   if ( masterlist(ff)%field%name /= fldname ) then
      call endrun ('GET_MASTERLIST_INDX: error finding field '//fldname//' on master list')
   end if

   get_masterlist_indx = ff
   return
   end function get_masterlist_indx

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


   subroutine bld_outfld_hash_tbls() 2,3
!-----------------------------------------------------------------------
!
! Purpose: Build primary and overlow hash tables for outfld processing.
!
! Steps:
!  1) Foreach field on masterlist, find all collisions.
!  2) Given the number of collisions, verify overflow table has sufficient
!     space.
!  3) Build primary and overlow indices.
!
!-----------------------------------------------------------------------

   implicit none
!
!  Local.
!
   integer :: ff
   integer :: ii
   integer :: itemp
   integer :: ncollisions
   integer :: hash_key

!
!  1) Find all collisions.
!
   tbl_hash_pri = 0

   do ff = 1, nfmaster
      hash_key = gen_hash_key(masterlist(ff)%field%name)
      tbl_hash_pri(hash_key) = tbl_hash_pri(hash_key) + 1
   end do

!
!  2) Count number of collisions and define start of a individual
!     collision's chain in overflow table. A collision is defined to be any
!     location in tbl_hash_pri that has a value > 1.
!
   ncollisions = 0
   do ii = 0, tbl_hash_pri_sz-1
      if ( tbl_hash_pri(ii) > 1 ) then  ! Define start of chain in O.F. table
         itemp = tbl_hash_pri(ii)
         tbl_hash_pri(ii) = -(ncollisions + 1)
         ncollisions = ncollisions + itemp + 1
      end if
   end do

   if ( ncollisions > tbl_hash_oflow_sz ) then
      write(6,*) 'BLD_OUTFLD_HASH_TBLS: ncollisions > tbl_hash_oflow_sz', &
              ncollisions, tbl_hash_oflow_sz
      call endrun()
   end if

!
!  3) Build primary and overflow tables.
!     i - set collisions in tbl_hash_pri to point to their respective
!         chain in the overflow table.
!
   tbl_hash_oflow = 0

   do ff = 1, nfmaster
      hash_key = gen_hash_key(masterlist(ff)%field%name)
      if ( tbl_hash_pri(hash_key) < 0 ) then
         ii = abs(tbl_hash_pri(hash_key))
         tbl_hash_oflow(ii) = tbl_hash_oflow(ii) + 1
         tbl_hash_oflow(ii+tbl_hash_oflow(ii)) = ff
      else
         tbl_hash_pri(hash_key) = ff
      end if
   end do

!
!  Dump out primary and overflow hashing tables.
!
!   if ( masterproc ) then
!      do ii = 0, tbl_hash_pri_sz-1
!         if ( tbl_hash_pri(ii) /= 0 ) write(6,666) 'tbl_hash_pri', ii, tbl_hash_pri(ii)
!      end do
!
!      do ii = 1, tbl_hash_oflow_sz
!         if ( tbl_hash_oflow(ii) /= 0 ) write(6,666) 'tbl_hash_oflow', ii, tbl_hash_oflow(ii)
!      end do
!
!      itemp = 0
!      ii = 1
!      do 
!         if ( tbl_hash_oflow(ii) == 0 ) exit
!         itemp = itemp + 1
!         write(6,*) 'Overflow chain ', itemp, ' has ', tbl_hash_oflow(ii), ' entries:'
!         do ff = 1, tbl_hash_oflow(ii)  ! dump out colliding names on this chain
!            write(6,*) '     ', ff, ' = ', tbl_hash_oflow(ii+ff), &
!                       ' ', masterlist(tbl_hash_oflow(ii+ff))%field%name
!         end do
!         ii = ii + tbl_hash_oflow(ii) +1 !advance pointer to start of next chain
!      end do
!   end if

   return
666 format(1x, a, '(', i4, ')', 1x, i6)

   end subroutine bld_outfld_hash_tbls

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


   subroutine bld_htapefld_indices 2,3
!-----------------------------------------------------------------------
!
! Purpose: Set history tape field indicies in masterlist for each
!          field defined on every tape.
!
! Note: because of restart processing, the actflag field is cleared and
!       then set only for active output fields on the different history
!       tapes.
!
!-----------------------------------------------------------------------

   implicit none
!
!  Arguments:
!

!
!  Local.
!
   integer :: f
   integer :: ff
   integer :: t

!
!  Initialize htapeindx to an invalid value.
!
   do ff = 1, nfmaster
      masterlist(ff)%htapeindx = -1
      masterlist(ff)%act_sometape = .false.
      masterlist(ff)%actflag = .false.
   end do

   do t = 1, ptapes
      do f = 1, nflds(t)
         ff = get_masterlist_indx(tape(t)%hlist(f)%field%name)
         if ( ff < 0 ) then
            write(6,*) 'BLD_HTAPEFLD_INDICES: something wrong, field not found on masterlist'
            write(6,*) 'BLD_HTAPEFLD_INDICES: t, f, ff = ', t, f, ff
            write(6,*) 'BLD_HTAPEFLD_INDICES: tape%name = ', tape(t)%hlist(f)%field%name
            call endrun
         end if
         masterlist(ff)%act_sometape = .true.
         masterlist(ff)%actflag(t) = .true.
         masterlist(ff)%htapeindx(t) = f
         if ( tape(t)%hlist(f)%field%name /= masterlist(ff)%field%name ) then
            write(6,*) 'BLD_HTAPEFLD_INDICES: hlist/masterlist inconsistency'
            write(6,*) 'BLD_HTAPEFLD_INDICES: t, f, ff = ', t, f, ff
            write(6,*) 'BLD_HTAPEFLD_INDICES: tape%name = ', tape(t)%hlist(f)%field%name
            write(6,*) 'BLD_HTAPEFLD_INDICES: masterlist%name = ', masterlist(ff)%field%name
            call endrun
         end if
      end do
    end do

!   if ( masterproc ) then
!      do ff = 1, nfmaster
!         write(6,*) masterlist(ff)%field%name, masterlist(ff)%act_sometape, &
!                    masterlist(ff)%actflag, masterlist(ff)%htapeindx
!      end do
!   end if
   return
   end subroutine bld_htapefld_indices
end module history