!
      module mkcdf_module
      use hist,only: history
      implicit none
      include 'netcdf.inc'
!
      integer,parameter :: mxflds = 100
      type netcdf_file
        character(len=120) :: filename
        integer :: ncid
        integer :: id_lon,id_lat,id_lev,id_nslice,id_zplev,id_hlev,
     |    id_hlevoh
        integer :: idv_lon,idv_lat,idv_lev,idv_slice,idv_zplev,idv_hlev,
     |    idv_hlevoh
        integer :: id_mtime,id_time,id_hvol_len,id_dtide,id_sdtide
        integer :: idv_mtime,idv_time,idv_ut,idv_step,idv_histvols,
     |    idv_yyddd,idv_year,idv_day,idv_gpi,idv_ncep,idv_gswmdi,
     |    idv_gswmsdi,idv_gswmnmidi,idv_gswmnmisdi,idv_f107a,idv_f107d,
     |    idv_hp,idv_cp,idv_by,idv_dtide,idv_sdtide
        integer :: nflds
        integer :: fids(mxflds) ! field id's
      end type netcdf_file
      type(netcdf_file) :: nc
!
      contains
!-----------------------------------------------------------------------
      subroutine nccreate(ncfilename,ncid)
!
! Args:
      character(len=*),intent(in) :: ncfilename
      integer,intent(out) :: ncid
!
! Local:
      integer :: istat,ier
      character(len=120) :: char120
!
! Create a netcdf file (file id ncid is returned by nf_create):
!
      istat = nf_create(ncfilename,NF_CLOBBER,ncid)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_create for netcdf ',
     |    'file ',a)") trim(ncfilename)
        call handle_ncerr(istat,char120)
        stop 'nf_create'
      else
        write(6,"('Created netcdf file ',a,' ncid=',i3)") 
     |    trim(ncfilename),ncid
      endif
      end subroutine nccreate
!-----------------------------------------------------------------------
      subroutine ncdefine(ncid,h,ncf)
      use input,only: modelhts
      use proc,only: diffs,spval
!
! Args:
      integer,intent(in) :: ncid ! opened netcdf file unit
      type(history),intent(in) :: h
      type(netcdf_file),intent(out) :: ncf
!
! Local:
      character(len=80) :: char80,units
      character(len=120) :: char120
      character(len=140) :: char140
      character(len=480) :: char480
      integer :: istat,idum
      integer :: ids1(1),ids2(2)
      integer,parameter :: ndtide=2, nsdtide=10
      character(len=8) ::
     |   rundate8,           ! current local date of run
     |   runtime8            ! current local time of run
      character(len=24) :: rundate
      character(len=16) ::
     |  host,              ! host machine
     |  system,            ! operating system of host (from pre-proc macros)
     |  logname            ! user login name
!
! Define time-dependent (non-dim) variables:
      ids1(1) = ncf%id_time	! for 1d time-unlimited vars
!
! Define model time var (3,unlimited):
! (ncl prints annoying error messages if this dimension is "mtime",
!  because it thinks coord vars can have only 1 dim)
!
      istat = nf_def_dim(ncid,"mtimedim",3,ncf%id_mtime)
      ids2(1) = ncf%id_mtime
      ids2(2) = ncf%id_time
      istat = nf_def_var(ncid,"mtime",NF_INT,2,ids2,ncf%idv_mtime)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining mtime variable')
      write(char80,"('model times (day, hour, minute)')")
      istat = nf_put_att_text(ncid,ncf%idv_mtime,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) then
        write(char80,"('ncdefine: Error defining long_name attribute ',
     |    'of mtime variable')")
        call handle_ncerr(istat,char80)
      endif
      istat = nf_put_att_text(ncid,ncf%idv_mtime,"units",17,
     |  'day, hour, minute')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining units of mtime variable')
!
! Define ut (decimal hours):
      ids1(1) = ncf%id_time
      istat = nf_def_var(ncid,"ut",NF_FLOAT,1,ids1,ncf%idv_ut)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining ut variable')
      write(char80,
     |  "('universal time (from model time hour and minute)')")
      istat = nf_put_att_text(ncid,ncf%idv_ut,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name attribute of ut variable')
      istat = nf_put_att_text(ncid,ncf%idv_ut,"units",5,'hours')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining units of ut var')
!
! Define time step (seconds):
      istat = nf_def_var(ncid,"timestep",NF_INT,0,idum,ncf%idv_step)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining timestep variable')
      istat = nf_put_att_text(ncid,ncf%idv_step,"long_name",8,
     |  'timestep')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name of timestep variable')
      istat = nf_put_att_text(ncid,ncf%idv_step,"units",7,'seconds')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining units of timestep variable')
!
! Define history volumes var (unlimited):
! (must define a character position dimension for char var)
! 2/16/06 btf: because h%mssvol may be blank, use disk file h%histfile
!
!     istat = nf_def_dim(ncid,"hvol_len",len(h%mssvol),ncf%id_hvol_len)
      istat=nf_def_dim(ncid,"hvol_len",len(h%histfile),ncf%id_hvol_len)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining ncf%id_hvol_len dimension')
      ids2(1) = ncf%id_hvol_len	! character position dimension
      ids2(2) = ncf%id_time
      istat = nf_def_var(ncid,"histvols",NF_CHAR,2,ids2,
     |  ncf%idv_histvols)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining histvols variable')
!     write(char80,"('mss paths to history volumes')")
      write(char80,"('history file')")
      istat = nf_put_att_text(ncid,ncf%idv_histvols,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) then
        write(char80,"('ncdefine: Error defining long_name attribute ',
     |    'of histvols variable')")
        call handle_ncerr(istat,char80)
      endif
!
! Define yyddd var (unlimited):
      istat = nf_def_var(ncid,"yyddd",NF_INT,1,ids1,ncf%idv_yyddd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining yyddd variable')
      write(char80,"('year and day of year (yyddd)')")
      istat = nf_put_att_text(ncid,ncf%idv_yyddd,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) then
        write(char80,"('ncdefine: Error defining long_name attribute ',
     |    'of yyddd variable')")
        call handle_ncerr(istat,char80)
      endif
!
! Define year var (unlimited):
      istat = nf_def_var(ncid,"year",NF_INT,1,ids1,ncf%idv_year)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining year variable')
      write(char80,"('calendar year')")
      istat = nf_put_att_text(ncid,ncf%idv_year,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name attribute of year variable')
!
! Define day var (unlimited):
      istat = nf_def_var(ncid,"day",NF_INT,1,ids1,ncf%idv_day)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefin: Error defining day variable')
      write(char80,"('calendar day')")
      istat = nf_put_att_text(ncid,ncf%idv_day,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefin: Error defining long_name attribute of day variable')
!
! Define gpi var (unlimited):
      istat = nf_def_var(ncid,"gpi",NF_INT,1,ids1,ncf%idv_gpi)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining gpi variable')
      write(char120,
     |  "('GPI flag (1 if using geophysical indices database for ',
     |  'f107d, f107a, hpower, and ctpoten)')")
      istat = nf_put_att_text(ncid,ncf%idv_gpi,"long_name",
     |  len_trim(char120),trim(char120))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name of gpi variable')
!
! Define ncep var (unlimited):
      istat = nf_def_var(ncid,"ncep",NF_INT,1,ids1,ncf%idv_ncep)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining ncep variable')
      write(char80,
     |  "('NCEP/NMC flag (1 if using ncep boundaries for t and z)')")
      istat = nf_put_att_text(ncid,ncf%idv_ncep,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name of ncep variable')
!
! Define gswmdi var (unlimited):
      istat = nf_def_var(ncid,"gswmdi",NF_INT,1,ids1,ncf%idv_gswmdi)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining gswmdi variable')
      write(char120,"('GSWMDI flag (1 if using GSWM diurnal tides as',
     |  ' boundary for Z,TN,UN,VN)')")
      istat = nf_put_att_text(ncid,ncf%idv_gswmdi,"long_name",
     |  len_trim(char120),trim(char120))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name of gswmdi variable')
!
! Define gswmsdi var (unlimited):
      istat = nf_def_var(ncid,"gswmsdi",NF_INT,1,ids1,ncf%idv_gswmsdi)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining gswmsdi variable')
      write(char120,"('GSWMSDI flag (1 if using GSWM semidiurnal tides',
     |   ' as boundary for Z,TN,UN,VN)')")
      istat = nf_put_att_text(ncid,ncf%idv_gswmsdi,"long_name",
     |  len_trim(char120),trim(char120))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name of gswmsdi variable')
!
! Define gswmnmidi var (unlimited):
      istat = nf_def_var(ncid,"gswmnmidi",NF_INT,1,ids1,
     |  ncf%idv_gswmnmidi)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining gswmnmidi variable')
      write(char120,"('GSWMNMIDI flag (1 if using GSWM nonmigrating',
     |  ' diurnal tides as boundary for Z,TN,UN,VN)')")
      istat = nf_put_att_text(ncid,ncf%idv_gswmnmidi,"long_name",
     |  len_trim(char120),trim(char120))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name of gswmnmidi variable')
!
! Define gswmnmisdi var (unlimited):
      istat = nf_def_var(ncid,"gswmnmisdi",NF_INT,1,ids1,
     |  ncf%idv_gswmnmisdi)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining gswmnmisdi variable')
      write(char120,"('GSWMNMISDI flag (1 if using GSWM nonmigrating',
     |  ' semidiurnal tides as boundary for Z,TN,UN,VN)')")
      istat = nf_put_att_text(ncid,ncf%idv_gswmnmisdi,"long_name",
     |  len_trim(char120),trim(char120))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name of gswmnmisdi variable')
!
! Define f107d var (unlimited):
      istat = nf_def_var(ncid,"f107d",NF_FLOAT,1,ids1,ncf%idv_f107d)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining f107d variable')
      write(char80,"('Daily f10.7 cm solar flux')")
      istat = nf_put_att_text(ncid,ncf%idv_f107d,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     | 'ncdefine: Error defining long_name attribute of f107d variable')
      units = ' '
      write(units,"('W/m^2Hz*1.e-22')")
      istat = nf_put_att_text(ncid,ncf%idv_f107d,"units",
     |  len_trim(units),units)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining units attribute of f107d variable')
!
! Define f107a var (unlimited):
      istat = nf_def_var(ncid,"f107a",NF_FLOAT,1,ids1,ncf%idv_f107a)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining f107a variable')
      write(char80,"('80-day average f10.7 cm solar flux')")
      istat = nf_put_att_text(ncid,ncf%idv_f107a,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     | 'ncdefine: Error defining long_name attribute of f107a variable')
      units = ' '
      write(units,"('W/m^2Hz*1.e-22')")
      istat = nf_put_att_text(ncid,ncf%idv_f107a,"units",
     |  len_trim(units),units)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining units attribute of f107a variable')
!
! Define hp var (unlimited):
      istat = nf_def_var(ncid,"hpower",NF_FLOAT,1,ids1,ncf%idv_hp)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefin: Error defining hp variable')
      write(char80,"('Hemispheric power')")
      istat = nf_put_att_text(ncid,ncf%idv_hp,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefin: Error defining long_name attribute of hp variable')
      istat = nf_put_att_text(ncid,ncf%idv_hp,"units",2,"Gw")
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefin: Error defining units attribute of hpower variable')
!
! Define cp var (unlimited):
      istat = nf_def_var(ncid,"ctpoten",NF_FLOAT,1,ids1,ncf%idv_cp)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefin: Error defining ctpoten variable')
      write(char80,"('Cross-tail potential (Volts)')")
      istat = nf_put_att_text(ncid,ncf%idv_cp,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |'ncdefin: Error defining long_name attribute of ctpoten variable')
      istat = nf_put_att_text(ncid,ncf%idv_cp,"units",5,"Volts")
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining units attribute of ctpoten variable')
!
! Define byimf var (unlimited):
      istat = nf_def_var(ncid,"byimf",NF_FLOAT,1,ids1,ncf%idv_by)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining byimf variable')
      write(char80,"('BY component of IMF')")
      istat = nf_put_att_text(ncid,ncf%idv_by,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     | 'ncdefine: Error defining long_name attribute of byimf variable')
!
! Dimension of diurnal and semidiurnal tide vars:
      istat = nf_def_dim(ncid,"dtidedim",ndtide,ncf%id_dtide)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining dtide dimension')
      istat = nf_def_dim(ncid,"sdtidedim",nsdtide,ncf%id_sdtide)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining sdtide dimension')
!
! Diurnal tide:
      ids2(1) = ncf%id_dtide
      ids2(2) = ncf%id_time
      istat = nf_def_var(ncid,"dtide",NF_FLOAT,2,ids2,ncf%idv_dtide)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining variable dtide')
      write(char80,"('amplitude and phase of diurnal tide')")
      istat = nf_put_att_text(ncid,ncf%idv_dtide,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     | 'ncdefine: Error defining long_name attribute of variable dtide')
!
! Semi-diurnal tide amp/phase:
      ids2(1) = ncf%id_sdtide
      ids2(2) = ncf%id_time
      istat = nf_def_var(ncid,"sdtide",NF_FLOAT,2,ids2,ncf%idv_sdtide)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining variable sdtide')
      write(char80,"('amplitude and phase of semi-diurnal tide')")
      istat = nf_put_att_text(ncid,ncf%idv_sdtide,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |'ncdefine: Error defining long_name attribute of variable sdtide')
!
! Define global attributes:
!
! Run date and time:
      call datetime(rundate8,runtime8)
      rundate = rundate8//' '//runtime8
      istat = nf_put_att_text(ncid,NF_GLOBAL,"rundate",
     |  len_trim(rundate),trim(rundate))
      if (istat /= NF_NOERR) then
        write(char120,"('ncdefine: Error return from nf_put_att_text ',
     |    'for rundate global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Logname:
      call getenv('LOGNAME',logname)
      if (len_trim(logname)==0) then
        write(6,"('>>> WARNING mkcdf: Cannot get LOGNAME environment ',
     |    'variable.',/)")
        logname = 'unknown'
      endif
      istat = nf_put_att_text(ncid,NF_GLOBAL,"logname",
     |  len_trim(logname),trim(logname))
      if (istat /= NF_NOERR) then
        write(char120,"('ncdefine: Error return from nf_put_att_text ',
     |    'for logname global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Model name:
      istat = nf_put_att_text(ncid,NF_GLOBAL,"model_name",
     |  len_trim(h%model_name),trim(h%model_name)) 
      if (istat /= NF_NOERR) then
        write(char120,"('ncdefine: Error return from nf_put_att_text ',
     |    'for model_name global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Model version:
      istat = nf_put_att_text(ncid,NF_GLOBAL,"model_version",
     |  len_trim(h%version),trim(h%version)) 
      if (istat /= NF_NOERR) then
        write(char120,"('ncdefine: Error return from nf_put_att_text ',
     |    'for version global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! History type (primary or secondary):
      if (h%issech) then
        istat = nf_put_att_text(ncid,NF_GLOBAL,"history_type",
     |    9,'secondary')
      else
        istat = nf_put_att_text(ncid,NF_GLOBAL,"history_type",
     |    7,'primary')
      endif
!
! Source is tgcmproc:
      istat = nf_put_att_text(ncid,NF_GLOBAL,"tgcmproc",
     |  3,'yes')
      write(char140,"('All fields have been post-processed according',
     |  ' to tgcmproc input parameters.')")
      istat = nf_put_att_text(ncid,NF_GLOBAL,"processed",
     |  len_trim(char140),trim(char140))
!
! Info about Z:
      if (modelhts<=0) then
        write(char140,"('Field Z was calculated with varying ',
     |    'gravity.')")
        istat = nf_put_att_text(ncid,NF_GLOBAL,"tgcmproc_modelhts",
     |    1,'0')
      else
        write(char140,"('Field Z was taken from the model history.')")
        istat = nf_put_att_text(ncid,NF_GLOBAL,"tgcmproc_modelhts",
     |    1,'1')
      endif
      istat = nf_put_att_text(ncid,NF_GLOBAL,"modelhts_explain",
     |  len_trim(char140),trim(char140))
!
! Difference fields?
      if (diffs) then
        istat = nf_put_att_text(ncid,NF_GLOBAL,"difference_fields",
     |    3,'yes') 
      else
        istat = nf_put_att_text(ncid,NF_GLOBAL,"difference_fields",
     |    2,'no') 
      endif
!
! Missing value:
      istat = nf_put_att_real(ncid,NF_GLOBAL,"missing_value",NF_FLOAT,
     |  1,spval)
      if (istat /= NF_NOERR) then
        write(char120,"('ncdefine: Error return from nf_put_att_real ',
     |    'for missing_value global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Explain zp_interp field attribute:
!  
      write(char480,"('If field attr zp_interp is yes then',
     |  ' field was interpolated from midpoint to interface zp levels.',
     |  ' If zp_interp is no then interpolation was not done ',
     |  ' (field was unknown or was already at interface zp on model',
     |  ' hist)')")
      istat = nf_put_att_text(ncid,NF_GLOBAL,"zp_interp_explain",
     |  len_trim(char480),trim(char480))
      if (istat /= NF_NOERR) then
        write(char120,"('ncdefine: Error return from nf_put_att_text ',
     |    'for zp_interp_explain global attribute')")
        call handle_ncerr(istat,char120)
      endif
      end subroutine ncdefine
!-----------------------------------------------------------------------
      subroutine ncvalues(ncid,h,nc,itime)
!
! Give values to history variables on open netcdf file ncid.
! These variables were earlier defined by ncdefine, and the
!   id's were put in nc.
!
! Args:
      integer,intent(in) :: ncid,itime
      type(history),intent(in) :: h
      type(netcdf_file),intent(in) :: nc
!
! Local:
      integer :: istat,mins,start_2d(2),count_2d(2),ipos
      character(len=120) :: char120
      character(len=80) :: char80
      real :: uthrs
      integer,parameter :: ndtide=2, nsdtide=10
      real :: var_sdtide(nsdtide),var_dtide(ndtide)
!
!     write(6,"('Enter ncvalues: itime=',i3,' h%mtime=',3i4)")
!    |  itime,h%mtime
!
! Give value to time var (total model minutes) for current model time:
      mins = h%mtime(1)*60*24+h%mtime(2)*60+h%mtime(3)
      istat = nf_put_var1_int(ncid,nc%idv_time,itime,mins)
      if (istat /= NF_NOERR) then

        write(6,"('Error defining time: ncid=',i3,' nc%idv_time=',
     |    i3)") ncid,nc%idv_time

        write(char120,"('Error from nf_put_var1_int defining ',
     |    'time at itime=',i3,': mins=',i10)") itime,mins
        call handle_ncerr(istat,char120)
      endif
!
! Give values to model time var:
      start_2d(1) = 1
      start_2d(2) = itime
      count_2d(1) = 3
      count_2d(2) = 1
      istat = nf_put_vara_int(ncid,nc%idv_mtime,start_2d,count_2d,
     |  h%mtime)
      if (istat /= NF_NOERR) then
        write(char120,"('Error defining mtime at itime=',i2,
     |    ' h%mtime=',3i4,' istat=',i5)") itime,h%mtime,istat
        call handle_ncerr(istat,char120)
      endif
!
! Give value to ut:
      uthrs = float(h%mtime(2))+float(h%mtime(3))/60.
      istat = nf_put_var1_real(ncid,nc%idv_ut,itime,uthrs)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'ut at mtime=',3i4,' uthrs=',f9.3)") h%mtime(:),uthrs
        call handle_ncerr(istat,char120)
      endif
!
! Give value to timestep var:
      istat = nf_put_var_int(ncid,nc%idv_step,h%step)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var_int defining ',
     |    'timestep: h%step=',i4)") h%step
        call handle_ncerr(istat,char120)
      endif
!
! Give value to histvols var for current model time:
! (note special treatment of character variable)
      start_2d(1) = 1
      start_2d(2) = itime
! 2/16/06 btf: because h%mssvol may be blank, use disk file h%histfile
!     count_2d(1) = len_trim(h%mssvol)+1 ! extra for null terminate
      count_2d(1) = len_trim(h%histfile)+1 ! extra for null terminate
      count_2d(2) = 1
!     char80 = h%mssvol
!     ipos = len_trim(h%mssvol)+1
      char80 = h%histfile
      ipos = len_trim(h%histfile)+1
      char80(ipos:ipos) = char(0)       ! null terminate
      istat = nf_put_vara_text(ncid,nc%idv_histvols,start_2d,count_2d,
     |  char80,80)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_vara_text defining ',
!    |    'histvol at itime=',i3,': h%mssvol=',a)") itime,trim(h%mssvol)
     |    'histvol at itime=',i3,': h%histfile=',a)") 
     |      itime,trim(h%histfile)
        call handle_ncerr(istat,char120)
      endif
!
! Give value to yyddd var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_yyddd,itime,h%iyd)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'yyddd at itime=',i3,': h%iyd=',i8)") itime,h%iyd
        call handle_ncerr(istat,char120)
      endif
!
! Give value to year var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_year,itime,h%year)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'year at itime=',i3,': h%year=',i8)") itime,h%year
        call handle_ncerr(istat,char120)
      endif
!
! Give value to day var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_day,itime,h%day)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'day at itime=',i3,': h%day=',i8)") itime,h%day
        call handle_ncerr(istat,char120)
      endif
!
! Give value to gpi var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_gpi,itime,h%gpi)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'gpi at itime=',i3,': h%gpi=',i8)") itime,h%gpi
        call handle_ncerr(istat,char120)
      endif
!
! Give value to ncep var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_ncep,itime,h%ncep)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'ncep at itime=',i3,': h%ncep=',i8)") itime,h%ncep
        call handle_ncerr(istat,char120)
      endif
!
! Give value to gswmdi var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_gswmdi,itime,h%gswmdi)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'gswmdi at itime=',i3,': h%gswmdi=',i8)") itime,h%gswmdi
        call handle_ncerr(istat,char120)
      endif
!
! Give value to gswmsdi var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_gswmsdi,itime,h%gswmsdi)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'gswmsdi at itime=',i3,': h%gswmsdi=',i8)") itime,h%gswmsdi
        call handle_ncerr(istat,char120)
      endif
!
! Give value to gswmnmidi var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_gswmnmidi,itime,h%gswmnmidi)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'gswmnmidi at itime=',i3,': h%gswmnmidi=',i8)") 
     |    itime,h%gswmnmidi
        call handle_ncerr(istat,char120)
      endif
!
! Give value to gswmnmisdi var for current model time:
      istat = nf_put_var1_int(ncid,nc%idv_gswmnmisdi,itime,h%gswmnmisdi)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'gswmnmisdi at itime=',i3,': h%gswmnmisdi=',i8)") 
     |    itime,h%gswmnmisdi
        call handle_ncerr(istat,char120)
      endif
!
! Give value to f107d var for current model time:
      istat = nf_put_var1_real(ncid,nc%idv_f107d,itime,h%f107d)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_real defining ',
     |    'f107d at itime=',i3,': h%f107d=',f8.2)") itime,h%f107d
        call handle_ncerr(istat,char120)
      endif
!
! Give value to f107a var for current model time:
      istat = nf_put_var1_real(ncid,nc%idv_f107a,itime,h%f107a)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_real defining ',
     |    'f107a at itime=',i3,': h%f107a=',f8.2)") itime,h%f107a
        call handle_ncerr(istat,char120)
      endif
!
! Give value to hpower var for current model time:
      istat = nf_put_var1_real(ncid,nc%idv_hp,itime,h%hpower)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_real defining ',
     |    'hp at itime=',i3,': h%hp=',f8.2)") itime,h%hpower
        call handle_ncerr(istat,char120)
      endif
!
! Give value to ctpoten var for current model time:
      istat = nf_put_var1_real(ncid,nc%idv_cp,itime,h%ctpoten)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_real defining ',
     |    'cp at itime=',i3,': h%cp=',f8.2)") itime,h%ctpoten
        call handle_ncerr(istat,char120)
      endif
!
! Give value to byimf var for current model time:
      istat = nf_put_var1_real(ncid,nc%idv_by,itime,h%byimf)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_real defining ',
     |    'byimf at itime=',i3,': h%byimf=',f8.2)") itime,h%byimf
        call handle_ncerr(istat,char120)
      endif
!
! Diurnal tide (dtide(time,2):
      start_2d(1) = 1
      start_2d(2) = itime
      count_2d(1) = ndtide
      count_2d(2) = 1
      var_dtide(:) = h%dtide(:)
      istat = nf_put_vara_real(ncid,nc%idv_dtide,start_2d,count_2d,
     |  var_dtide)
!
! Semi-diurnal tide (sdtide(time,10):
      start_2d(1) = 1
      start_2d(2) = itime
      count_2d(1) = nsdtide
      count_2d(2) = 1
      var_sdtide(:) = h%sdtide(:)
      istat = nf_put_vara_real(ncid,nc%idv_sdtide,start_2d,count_2d,
     |  var_sdtide)
! 
      end subroutine
!-----------------------------------------------------------------------
      subroutine mkcdf(f,fcntr,nf,h,flnm_cdf,itime,ier)
!
! Make netcdf file containing fields for each history. Fields are
!  3d or 4d including time (lon,lat,[lev],time).
! This routine is called at each time iteration (iter index itime).
! When itime == 1 the netcdf file is named and created, and dimensions
!   and variables are defined (flnm_cdf is output as the file name).
! When itime > 1, but < ntimes, fields are appended to the open file.
! When itime == ntimes, fields are appended, last_history global 
!   attribute is added, and the file is closed.
! Netcdf files made w/ this routine from time-gcm histories (nlon=73,
!   nlat=36,nlev=45) are about .473 Mb per field per history. E.g., 
!   a 15-field 25-time file will be about 169Mb.
!
      use proc
      use hist,only: history,tgcmheader
      use fields,only: field,oh_kmbot
      use input,only: mtimes,modelhts,ncid_out,
!
! real cdf_htscale(3): user-provided bottom,top,delta for y-axis height column
! real cdf_zprange(2): user-provided bottom,top for zp y-axis pressure column
!   (cdf_zprange is ignored if cdf_htscale is set)
!
     |  cdf_htscale,
     |  cdf_zprange,
     |  cdf_zonalmean
      implicit none
!
! Args:
      integer,intent(in) :: nf,itime
      type(field),intent(in) :: f(nf),fcntr(nf)
      type(history),intent(in) :: h
      character(len=120),intent(out) :: flnm_cdf	! netcdf file name
      integer,intent(out) :: ier
!
! Locals that must be saved for calls after 1st time:
      integer,save :: ncid 	! netcdf file id
      real,allocatable,save ::
     |  hlev(:),    ! column of heights (km) from user-provided bot,top,delta
     |  zplev(:),   ! column of log pressure coords (standard model zp)
     |  hlev_oh(:)  ! column of OH heights (for OH-VIB or OH-BAND fields)  
      integer,save :: nzplev,nhlev,nhlev_oh,ishape(3)
      integer,parameter :: ndtide=2, nsdtide=10
!
! Locals that need not be saved:
      integer :: 
     |  istat,idum,i,k,j,kk,ii,ipos,iier,ixz,! status, loop vars, char pos
     |  ids1(1),ids2(2),ids3(3),ids4(4),! vectors of dim id's
     |  idv_flds(nf),			! field variable id's
     |  nflds,				! number of fields on netcdf file
     |  nnlev,				! temp for number of levels
     |  mins,				! total model time in minutes
     |  start_4d(4),count_4d(4),	! start loc and count for 4d vars
     |  start_3d(3),count_3d(3),	! start loc and count for 3d vars
     |  start_2d(2),count_2d(2)		! start loc and count for 2d vars
      real :: uthrs                     ! decimal ut (hours)
      real,allocatable :: fdata(:,:,:)	! temp f var
      real,allocatable :: fcdata(:,:,:) ! temp f var for control (if diffs)
      real,allocatable :: var(:,:,:)	! temp var
      real,allocatable :: varcntr(:,:,:)! temp var for control (if diffs)
      real,allocatable :: fdiffs(:,:,:)	! difference fields
      real,allocatable :: var2(:,:)	! temp var
      real,allocatable :: var2cntr(:,:)	! temp var for control (2d)
      real,allocatable :: fz(:,:,:)     ! Z pert field, saved at Zp pressures
      real,allocatable :: fzcntr(:,:,:) ! Z cntr field, saved at Zp pressures
      character(len=4096) :: field_names
      character(len=6)   :: chnlev
      character(len=8)   :: charmtime,charmtime1,fname
      character(len=80)  :: char80,units
      character(len=120) :: char120
      character(len=140) :: char140
      character(len=8) ::
     |   rundate8,           ! current local date of run
     |   runtime8            ! current local time of run
      character(len=24) :: rundate
      character(len=16) ::
     |  host,              ! host machine
     |  system,            ! operating system of host (from pre-proc macros)
     |  logname            ! user login name
      logical :: 
     |  mkhts,   ! if true, user has provided cdf_htscale(3)
     |  htindep  ! true if a field is height-independent
      integer :: loginterp,length
      real :: fmin,fmax
      real :: var_sdtide(nsdtide),var_dtide(ndtide)
!
! External:
      integer,external :: ixfindc
      real,external :: fmean
!
!     write(6,"('Enter mkcdf: itime=',i4)") itime
!
! Error init:
      ier = 0
!
! Set mkhts if user provided cdf_htscale (cdf_htscale was init to
! spval by input before namelist read):
      mkhts = .false.
      if (any(cdf_htscale/=spval)) mkhts = .true.
!
! If cdf_zprange was not set by user, set it to zpbot,zptop
! (if cdf_zprange was given by user, it was validated by sub setflev)
      if (all(cdf_zprange==spval)) then
        cdf_zprange(1) = h%zpb
        cdf_zprange(2) = h%zpt
      endif 
      nzplev = ifix((cdf_zprange(2)-cdf_zprange(1)) / 
     |  dlev + 1.0000001)
!
! If this is not the 1st time, just grow the variables along the
! unlimited time dimension:
!
      if (itime > 1) goto 100
!
! If this is the 1st time iteration, create file, define variables, 
!   dimensions, and attributes:
! bf 1/26/03: use a simple disk file name for netcdf file. Later it
!   is sent to remote sendcdf.
! 8/29/05 btf: flnm_cdf is now namelist read parameter.
!     flnm_cdf = 'tgcmproc.nc'
!     nc%filename = 'tgcmproc.nc'
      nc%filename = flnm_cdf
!
! Create new netcdf data set:
      call nccreate(nc%filename,nc%ncid) ! filename is input, ncid is output
      ncid = nc%ncid
!
! Define dimensions (field vars will be (time,lev,lat,lon) in
! CDL notation). (time was defined by ncdefine)
!
! Time is the unlimited dimension:
      istat = nf_def_dim(ncid,"time",NF_UNLIMITED,nc%id_time)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining time dimension')
!
! Time (total model time in minutes):
      ids1(1) = nc%id_time
      istat = nf_def_var(ncid,"time",NF_INT,1,ids1,nc%idv_time)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining time dimension variable')
      write(char80,"('total model time (day*24*60+hour*60+min)')")
      istat = nf_put_att_text(ncid,nc%idv_time,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining long_name of time dimension variable')
      istat = nf_put_att_text(ncid,nc%idv_time,"units",7,'minutes')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'ncdefine: Error defining units of time dimension variable')
!
! Longitude dimension:
      istat = nf_def_dim(ncid,"lon",nlon,nc%id_lon)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining longitude dimension')
!
! Longitudes:
      ids1(1) = nc%id_lon
      istat = nf_def_var(ncid,"lon",NF_FLOAT,1,ids1,nc%idv_lon)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining longitude dimension variable')
      write(char80,"('geographic longitude (-west, +east)')")
      istat = nf_put_att_text(ncid,nc%idv_lon,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining long_name of longitude dimension variable')
      istat = nf_put_att_text(ncid,nc%idv_lon,"units",12,'degrees_east')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining units of longitude dimension variable')
!
! Latitude dimension:
      istat = nf_def_dim(ncid,"lat",nlat,nc%id_lat)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining latitude dimension')
!
! Latitudes:
      ids1(1) = nc%id_lat
      istat = nf_def_var(ncid,"lat",NF_FLOAT,1,ids1,nc%idv_lat)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining latitude dimension variable')
      write(char80,"('geographic latitude (-south, +north)')")
      istat = nf_put_att_text(ncid,nc%idv_lat,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining long_name of latitude dimension variable')
      istat = nf_put_att_text(ncid,nc%idv_lat,"units",13,
     |  'degrees_north')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining units of latitude dimension variable')
!
! Vertical dimension (zp, user-defined height scale, or OH height scale)
      do i=1,nf
        select case (trim(f(i)%vtype))
        case('ZP')
!
! Fields are either in zp or in heights, not both (except for 
! vtype==HEIGHT, i.e., OH, which is always in height). For
! vtype==ZP, if mkhts is set (fields in height), then fields
! will be interpolated to heights later. In either case, fields
! will be dimensioned "lev", and coord var will be lev(lev).
! (vtype=HEIGHT fields will be dimensioned "hlev")
!
          if (.not.mkhts) then ! fields in zp (not interp to height)
            if (.not.allocated(zplev)) then
!
! Lev dimension var (will be zp):
              istat = nf_def_dim(ncid,"lev",nzplev,nc%id_zplev)
              if (istat /= NF_NOERR) call handle_ncerr(istat,
     |          'Error defining zplev (log pressure) dimension')
!
! Set up zp pressure array zplev(nzplev) (may have been user-limited
! via cdf_zprange):
              allocate(zplev(nzplev),stat=ier)
              if (ier /= 0) call allocerr(ier,"mkcdf allocating zplev")
              do k=1,nzplev
                zplev(k) = cdf_zprange(1)+(k-1)*dlev
              enddo
!             if (.not.mkhts) write(6,"('mkcdf: nzplev=',i3,' zplev=',/,
!    |          (6f9.2))") nzplev,zplev
!
! Define zplev(nzplev) coordinate variable and attributes:
              ids1(1) = nc%id_zplev
              istat = nf_def_var(ncid,"lev",NF_FLOAT,1,ids1,
     |          nc%idv_zplev)
              if (istat /= NF_NOERR) call handle_ncerr(istat,
     |         'Error defining pressure (zplev) dimension variable')
              write(char80,"('LOG PRESSURE')")
              istat = nf_put_att_text(ncid,nc%idv_zplev,"long_name",
     |          len_trim(char80),trim(char80))
              if (istat /= NF_NOERR) call handle_ncerr(istat,
     |          'Error defining long_name of zplev coord variable')
              istat = nf_put_att_text(ncid,nc%idv_zplev,"units",8,
     |          'ln(p0/p)')
              if (istat /= NF_NOERR) call handle_ncerr(istat,
     |          'Error defining units of zplev coord variable')
!
! All fields are at "interface zp" (i.e., fields have been converted from 
!   "midpoint" levels of the model history)
              istat = nf_put_att_text(ncid,nc%idv_zplev,"zplev",16,
     |          'interface levels')
            endif ! .not.allocated(zplev)
!
! Set up optional user-requested height scale if cdf_heights(3) 
! was set by user. Allocate and define hlev(nhlev) based on cdf_heights.
!
          else ! mkhts is set -- fields will be interp to height
            if (.not.allocated(hlev)) then
!
! User requested height scale:
              nhlev = ifix((cdf_htscale(2)-cdf_htscale(1)) / 
     |          cdf_htscale(3) + 1.0000001)
              allocate(hlev(nhlev),stat=ier)
              if (ier /= 0) call allocerr(ier,"mkcdf allocating hlev")
              do k=1,nhlev
                hlev(k) = cdf_htscale(1)+(k-1)*cdf_htscale(3)
              enddo
              write(6,"('User requested height scale ',f10.2,' to ',
     |          f10.2,' by ',f10.2,': nhlev=',i4,' hlev = ',/,
     |          (6e12.4))") cdf_htscale(:),nhlev,hlev(:)
!
! Lev dimension variable (will be height):
              istat = nf_def_dim(ncid,"lev",nhlev,nc%id_hlev)
              if (istat /= NF_NOERR) call handle_ncerr(istat,
     |          'Error defining hlev (height) dimension')
!
! Lev(lev) coord variable and attributes:
              ids1(1) = nc%id_hlev
              istat = nf_def_var(ncid,"lev",NF_FLOAT,1,ids1,nc%idv_hlev)
              if (istat /= NF_NOERR) call handle_ncerr(istat,
     |          'Error defining height (hlev) coord variable')
              write(char80,"('HEIGHT')")
              istat = nf_put_att_text(ncid,nc%idv_hlev,"long_name",
     |          len_trim(char80),trim(char80))
              if (istat /= NF_NOERR) call handle_ncerr(istat,
     |          'Error defining long_name of height coord variable')
              istat = nf_put_att_text(ncid,nc%idv_hlev,"units",2,'km')
              if (istat /= NF_NOERR) call handle_ncerr(istat,
     |          'Error defining units of height dimension variable')
            endif
          endif ! mkhts
!
! OH height scale:
        case('HEIGHT')
          if (.not.allocated(hlev_oh)) then
!
! hlev_oh dimension:
            nhlev_oh = f(i)%nlev
            istat = nf_def_dim(ncid,"hlev",nhlev_oh,nc%id_hlevoh)
            if (istat /= NF_NOERR) call handle_ncerr(istat,
     |        'Error defining hlev_oh (oh heights) dimension')
!
! hlev_oh coord variable and attributes:
            allocate(hlev_oh(nhlev_oh),stat=ier)
            if (ier /= 0) call allocerr(ier,"mkcdf allocating hlev_oh")
            do k=1,nhlev_oh
              hlev_oh(k) = oh_kmbot+float(k-1)*f(i)%dlev
            enddo
            ids1(1) = nc%id_hlevoh
            istat =nf_def_var(ncid,"hlev",NF_FLOAT,1,ids1,nc%idv_hlevoh)
            if (istat /= NF_NOERR) call handle_ncerr(istat,
     |       'Error defining OH height (hlev_oh) coord variable')
            write(char80,"('OH height scale')")
            istat = nf_put_att_text(ncid,nc%idv_hlevoh,"long_name",
     |        len_trim(char80),trim(char80))
            if (istat /= NF_NOERR) call handle_ncerr(istat,
     |        'Error defining long_name of hlev_oh coord variable')
            istat = nf_put_att_text(ncid,nc%idv_hlevoh,"units",2,'km')
            if (istat /= NF_NOERR) call handle_ncerr(istat,
     |        'Error defining units of hlev_oh coord variable')
          endif
        case ('HT-INDEP')
! height-independent field -- do nothing here (field will be allocated later)
        case default
          write(6,"('>>> warning mkcdf: unknown vtype=',a,' for field ',
     |      a)") f(i)%vtype,f(i)%fname8
        end select
      enddo ! i=1,nf
!
! Define history vars on new file:
!
      call ncdefine(nc%ncid,h,nc) ! ncid and h are input, components of nc are output
!
! Define 4d field vars (always include Z):
!
      ids4(1) = nc%id_lon
      ids4(2) = nc%id_lat
      ids4(4) = nc%id_time
      nnlev = 0
      nflds = 0
!
! Define fields and their attributes:
!
      floop1: do i=1,nf
        if ((.not.f(i)%requested.or..not.associated(f(i)%data)).and.
     +       f(i)%fname8 /= 'Z       ') cycle floop1
        nflds = nflds+1
        htindep = .false.
        if (trim(f(i)%vtype)=='ZP') then
          ids4(3) = nc%id_zplev
          if (mkhts) ids4(3) = nc%id_hlev
        elseif (trim(f(i)%vtype)=='HEIGHT') then
          ids4(3) = nc%id_hlevoh
        elseif (trim(f(i)%vtype)=='HT-INDEP') then
          ids4(3) = 0
          htindep = .true.
        endif
!
! Netcdf does not allow '+' in a var field name, so change to 'P'
!   if necessary (e.g., O2+ will be O2P on the netcdf file).
! Also change any fields with '/' in the name to '_' (e.g.,
!   N2/O will become N2_O). 
!
        fname = f(i)%fname8
        ipos = index(fname,'+')
        if (ipos > 0) fname(ipos:ipos) = 'P'
        ipos = index(fname,'/')
        if (ipos > 0) fname(ipos:ipos) = '_'
!
! Define the field var:
!
! If neither htindep nor zonalmean, field is the normal 4-d:
        if (.not.htindep.and.cdf_zonalmean <= 0.) then
          istat = nf_def_var(ncid,fname,NF_FLOAT,4,ids4,idv_flds(i))
          if (istat /= NF_NOERR) then
            write(char120,"('Error return from nf_def_var for ',
     |        '3d field var ',a,': i=',i2,' htindep=',l1,
     |        ' cdf_zonalmean=',f3.1)") fname,i,htindep,cdf_zonalmean
            call handle_ncerr(istat,char120)
          else
            write(6,"('Defined 4-d field ',a,'(lon,lat,lev,time) on ',
     |        'netcdf file ',a,' idv=',i3)") fname,trim(nc%filename),
     |        idv_flds(i)
          endif
!
! Height independent field without zonal means is 3-d (lon,lat,time):
        elseif (htindep.and.cdf_zonalmean <= 0.) then
          ids3(1) = nc%id_lon
          ids3(2) = nc%id_lat
          ids3(3) = nc%id_time
          istat = nf_def_var(ncid,fname,NF_FLOAT,3,ids3,idv_flds(i))
          if (istat /= NF_NOERR) then
            write(char120,"('Error return from nf_def_var for ',
     |        '3d field var ',a,': i=',i2,' htindep=',l1,
     |        ' cdf_zonalmean=',f3.1)") fname,i,htindep,cdf_zonalmean
            call handle_ncerr(istat,char120)
          else
            write(6,"('Defined 3-d field ',a,'(lon,lat,time)',
     |        ' on netcdf file ',a)") fname,trim(nc%filename)
          endif
!
! Non height-independent field with zonal means is 3-d (lat,vert,time):
        elseif (.not.htindep.and.cdf_zonalmean >= 1.) then
          ids3(1) = nc%id_lat
          ids3(2) = ids4(3) ! zp or ht, see above
          ids3(3) = nc%id_time
          istat = nf_def_var(ncid,fname,NF_FLOAT,3,ids3,idv_flds(i))
          if (istat /= NF_NOERR) then
            write(char120,"('Error return from nf_def_var for ',
     |        '3d field var ',a,': i=',i2,' htindep=',l1,
     |        ' cdf_zonalmean=',f3.1)") fname,i,htindep,cdf_zonalmean
            call handle_ncerr(istat,char120)
          else
            write(6,"('Defined 3-d field ',a,'(lat,lev,time)',
     |        ' on netcdf file ',a)") fname,trim(nc%filename)
          endif
!
! If field is htindep and we are making zonal means, field is 2-d (lat,time):
        else ! both htindep and cdf_zonalmean
          ids2(1) = nc%id_lat
          ids2(2) = nc%id_time
          istat = nf_def_var(ncid,fname,NF_FLOAT,2,ids2,idv_flds(i))
          if (istat /= NF_NOERR) then
            write(char120,"('Error return from nf_def_var for ',
     |        '2d field var ',a,': i=',i2,' htindep=',l1,
     |        ' cdf_zonalmean=',f3.1)") fname,i,htindep,cdf_zonalmean
            call handle_ncerr(istat,char120)
          else
            write(6,"('Defined 2-d field ',a,'(lat,time)',
     |        ' on netcdf file ',a)") fname,trim(nc%filename)
          endif
        endif
!
! Assign attributes to field variables:
!
! Long name and units. Add zonal mean to long name if cdf_zonalmean is set.
        if (cdf_zonalmean <= 0.) then
          istat = nf_put_att_text(ncid,idv_flds(i),"long_name",
     |      len_trim(f(i)%fname56),trim(f(i)%fname56))
        else
          write(char80,"('ZONAL MEAN ',a)") trim(f(i)%fname56)
          istat = nf_put_att_text(ncid,idv_flds(i),"long_name",
     |      len_trim(char80),trim(char80))
        endif
        istat = nf_put_att_text(ncid,idv_flds(i),"units",
     |    len_trim(f(i)%units),trim(f(i)%units))
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_att_text',
     |      ' of long_name and units for field ',a)") fname
          call handle_ncerr(istat,char120)
        endif
!
! Field "type" (DENSITY, VECTOR, SCALAR, etc):
        istat = nf_put_att_text(ncid,idv_flds(i),"type",
     |    len(f(i)%type),f(i)%type)
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_att_text',
     |      ' of type ',a,' for field ',a)") f(i)%type,fname
          call handle_ncerr(istat,char120)
        endif
!
! If field was calculated at midpoint levels (f(i)%zptype=='MIDPOINTS'), then
!   the processor interpolated the zp dimension to provide the field
!   at interface levels (zp_interp=='yes'). Conversely if f(i)%zptype='INTERFACES',
!   then zp_interp=='no' (only fields NE,W,Z,POTEN). Note all "unknown"
!   fields (secondary histories) are zp_interp=='no', i.e., the processor
!   did not interpolate in zp.
!
        if (f(i)%known) then ! known field
          if (trim(f(i)%zptype)=="MIDPOINTS") then
            istat = nf_put_att_text(ncid,idv_flds(i),"zp_interp",
     |        3,'yes')
          elseif (trim(f(i)%zptype)=="INTERFACES") then
            istat = nf_put_att_text(ncid,idv_flds(i),"zp_interp",
     |        2,'no')
          else
            write(6,"('>>> mkcdf: unknown f(i)%zptype=',a)")
     |        trim(f(i)%zptype)
            istat = nf_put_att_text(ncid,idv_flds(i),"zp_interp",
     |        7,'unknown')
          endif
        else                 ! unknown field 
          istat = nf_put_att_text(ncid,idv_flds(i),"zp_interp",
     |      2,'no')
        endif
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_att_text',
     |      ' of attribute zp_interp for field ',a,': known=',l1,
     |      ' zptype=',a )") fname,f(i)%known,f(i)%zptype
          call handle_ncerr(istat,char120)
        endif
!
! If difference fields, define type (percent or raw):
        if (diffs) then
!         istat = nf_put_att_text(ncid,idv_flds(i),"differences",
!    |      len(f(i)%difftype),f(i)%difftype)
          if (f(i)%difftype == 'RAW') then
            istat = nf_put_att_text(ncid,idv_flds(i),"differences",
     |        5,'DIFFS')
          else
            istat = nf_put_att_text(ncid,idv_flds(i),"differences",
     |        6,'%DIFFS')
          endif
          if (istat /= NF_NOERR) then
            write(char120,"('Error return from nf_put_att_text',
     |        ' of difftype ',a,' for field ',a)") 
     |        f(i)%difftype,fname
            call handle_ncerr(istat,char120)
          endif
        endif
      enddo floop1
!
! Number of fields:
      istat = nf_put_att_int(ncid,NF_GLOBAL,"number_of_fields",
     |  NF_INT,1,nf)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_int ',
     |    'for number_of_fields global attribute: nf=',i4)") 
     |    nf
        call handle_ncerr(istat,char120)
      else
        write(6,"('mkcdf: number_of_fields attribute=',i4)") nf
      endif
!
! Field names:
      field_names = " " ! len=4096
      do i=1,nf
        if ((f(i)%requested.and.associated(f(i)%data)).or.
     |       f(i)%fname8 == 'Z       ') then
          fname = f(i)%fname8
          ipos = index(fname,'+')
          if (ipos > 0) fname(ipos:ipos) = 'P'
          ipos = index(fname,'/')
          if (ipos > 0) fname(ipos:ipos) = '_'
          ii = ii+1
          field_names = trim(field_names)//' '//trim(fname)
        endif
      enddo
      istat = nf_put_att_text(ncid,NF_GLOBAL,"field_names",
     |  len_trim(field_names),trim(field_names))
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     |    'for field_names')")
        call handle_ncerr(istat,char120)
      endif
!
! Number of times:
      istat = nf_put_att_int(ncid,NF_GLOBAL,"number_of_times",
     |  NF_INT,1,ntimes)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_int ',
     |    'for ntimes global attribute: ntimes=',i4)") ntimes
        call handle_ncerr(istat,char120)
      endif
!
! First history (mssvol and mtime):
! (note last history is added when itime==ntimes below)
! 2/16/06 btf: because h%mssvol may be blank, use disk file h%histfile
!
      write(char120,"(a,', ',i3.3,':',i2.2,':',i2.2)") 
     |  trim(h%histfile),h%mtime
      istat = nf_put_att_text(ncid,NF_GLOBAL,"first_history",
     |  len_trim(char120),trim(char120))
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     |    'for first_history global attribute: ',a)") trim(char80)
        call handle_ncerr(istat,char120)
      endif
!
! Exit define mode:
      istat = nf_enddef(ncid)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error return from nf_enddef')
!
! Give values to dimension variables:
      istat = nf_put_var_real(ncid,nc%idv_lon,gcmlon) 
      istat = nf_put_var_real(ncid,nc%idv_lat,gcmlat) 
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_var_real ',
     |    ' to assign values to lon,lat dimension vars')")
        call handle_ncerr(istat,char120)
      endif
!
! Give values to zp or height coord vars:
      if (allocated(zplev)) then
        istat = nf_put_var_real(ncid,nc%idv_zplev,zplev) 
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_var_real ',
     |      ' to assign values to zplev coord var')")
          call handle_ncerr(istat,char120)
        endif
      endif
      if (allocated(hlev)) then
        istat = nf_put_var_real(ncid,nc%idv_hlev,hlev) 
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_var_real ',
     |      ' to assign values to hlev coord var')")
          call handle_ncerr(istat,char120)
        endif
      endif
      if (allocated(hlev_oh)) then
        istat = nf_put_var_real(ncid,nc%idv_hlevoh,hlev_oh)
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_var_real ',
     |      ' to assign values to hlev_oh coord var')")
          call handle_ncerr(istat,char120)
        endif
      endif
!
! End itime==1
!
      if (ier /= 0) return
 100  continue
!
! Give values to history vars:
      call ncvalues(ncid,h,nc,itime)
!
! Save Z in pressure:
      ixz = ixfindc(f%fname8,nf,'Z       ')
      if (ixz <= 0) write(6,"('>>> WARNING mkcdf: cannot find ',
     |  'Z on history: ixz=',i3)") ixz
      if (allocated(fz)) deallocate(fz)
      allocate(fz(nlon,nlat,nzplev),stat=iier)
      if (iier /= 0) call allocerr(ier,
     |  "mkcdf allocate fz var with (nlon,nlat,nzplev)")
      fz(:,:,:) = f(ixz)%data(:,:,:)
      call fminmax(fz,nlon*nlat*nzplev,fmin,fmax)
!     write(6,"('mkcdf: ixz=',i3,' allocated fz: nlon=',i3,' nlat=',i3,
!    |  ' nzplev=',i3,' fz min,max=',2e12.4)") ixz,nlon,nlat,nzplev,
!    |  fmin,fmax
      if (diffs) then
        ixz = ixfindc(fcntr%fname8,nf,'Z       ')
        if (ixz <= 0) write(6,"('>>> WARNING mkcdf: cannot find ',
     |    'Z on cntr history: ixz=',i3)") ixz
        if (allocated(fzcntr)) deallocate(fzcntr)
        allocate(fzcntr(nlon,nlat,nzplev),stat=iier)
        if (iier /= 0) call allocerr(ier,
     |    "mkcdf allocate fzcntr var with (nlon,nlat,nzplev)")
        fzcntr(:,:,:) = fcntr(ixz)%data(:,:,:)
      endif
!
! If doing height interp, add Z at pressure levels.
! This means adding new dimension zplev, coord var zplev(zplev), 
!   and 3-d field ZP(lon,lat,zplev,time). This is called every time.
!
      if (mkhts) then
        call addzplev(ncid,fz,nlon,nlat,nzplev,cdf_zprange,dlev,itime)
      endif
!
! Give values to field variables for current model time:
!
      start_4d(1:3) = 1
      start_4d(4) = itime
      count_4d(1) = nlon
      count_4d(2) = nlat
      count_4d(4) = 1
!
! Fields loop:
      floop2: do i=1,nf
        if ((.not.f(i)%requested.or..not.associated(f(i)%data)).and.
     |       f(i)%fname8 /= 'Z       ') cycle floop2
!
! Get var id from name (actually only necessary when itime > 1)
        fname = f(i)%fname8
        ipos = index(fname,'+')
        if (ipos > 0) fname(ipos:ipos) = 'P'
        ipos = index(fname,'/')
        if (ipos > 0) fname(ipos:ipos) = '_'
        istat = nf_inq_varid(ncid,fname,idv_flds(i))
        if (istat /= NF_NOERR) then
          write(char120,"('Error getting id of field var ',a)") fname
          call handle_ncerr(istat,char120)
        endif
!
! Netcdf doesn't like the pointer f(i)%data to be passed directly
!  to nf_put_vara_real, so transfer to tmp array var first. 
! Note alloc/dealloc of var must be inside loop because f(i)%nlev 
!   may vary.
!
        if (allocated(fdata))    deallocate(fdata)
        if (allocated(fcdata))   deallocate(fcdata)
        if (allocated(var))      deallocate(var)
        if (allocated(var2))     deallocate(var2)
        if (allocated(fdiffs))   deallocate(fdiffs)
        if (allocated(varcntr))  deallocate(varcntr)
        if (allocated(var2cntr)) deallocate(var2cntr)
!
! Allocate temp arrays to hold f and fcntr data:
        ishape = shape(f(i)%data(:,:,:)) ! lon,lat,lev
        allocate(fdata(ishape(1),ishape(2),ishape(3)),stat=iier)
        if (iier /= 0) then 
          write(6,"('>>> Error allocating fdata with ishape=',3i4)")
     |      ishape
          call allocerr(iier,"mkcdf allocate fdata with ishape")
        endif
        if (diffs) then ! assume shape of fcntr is same as f
          allocate(fcdata(ishape(1),ishape(2),ishape(3)),stat=iier)
          if (iier /= 0) then 
            write(6,"('>>> Error allocating fcdata with ishape=',3i4)")
     |        ishape
            call allocerr(ier,"mkcdf allocate fcdata with ishape")
          endif
        endif
!
! Transfer f(i)%data to fdata. If doing zonal means, store means
! in first index of longitude dimension (i==1).
        if (cdf_zonalmean <= 0.) then
          fdata(:,:,:) = f(i)%data(:,:,:)
          if (diffs) fcdata(:,:,:) = fcntr(i)%data(:,:,:)
        else ! find zonal means
          do k=1,ishape(3)
            do j=1,ishape(2)
              fdata(1,j,k) = fmean(f(i)%data(:,j,k),ishape(1),spval,1)
              if (diffs) fcdata(:,j,k) = fmean(fcntr(i)%data(:,j,k),
     |          ishape(1),spval,1)
            enddo
          enddo
        endif
!
! Allocate temp var according to vertical dimension:
        if (trim(f(i)%vtype)=='ZP') then
          if (.not.mkhts) then
            count_4d(3) = nzplev
            allocate(var(nlon,nlat,nzplev),stat=iier)
            if (iier /= 0) call allocerr(ier,
     |        "mkcdf allocate tmp var with (nlon,nlat,nzplev)")
!
! Take diffs:
            if (diffs) then
!
! 6/11/07 btf: allocate fdiffs at full npress, to be compatable with
! fdata,fcdata for mkdiffs. This fixed a seg fault under AIX blueice,
! however, it is still seg faulting on hao Linux oreo.
! 
              allocate(fdiffs(nlon,nlat,npress),stat=iier)
              if (iier /= 0) call allocerr(ier,
     |          "mkcdf allocate fdiffs with (nlon,nlat,npress)")
              length = size(fdata)
              call mkdiffs(fdata,fcdata,fdiffs,length,
     |          f(i)%difftype)
!
! Full or partial zp dimension (diffs):
              if (cdf_zprange(1)==h%zpb.and.cdf_zprange(2)==h%zpt) then
                var(:,:,:) = fdiffs(:,:,:)
              else
                kk = 0
                do k=1,npress
                  if (gcmlev(k)>=zplev(1).and.gcmlev(k)<=zplev(nzplev))
     |              then
                    kk = kk+1
                    var(:,:,kk) = fdiffs(:,:,k)
                  endif
                enddo
              endif
            else ! not diffs
! Full or partial zp dimension (no diffs):
              if (cdf_zprange(1)==h%zpb.and.cdf_zprange(2)==h%zpt) then
                var(:,:,:) = fdata(:,:,:)
              else
                kk = 0
                do k=1,npress
                  if (gcmlev(k)>=zplev(1).and.gcmlev(k)<=zplev(nzplev))
     |              then
                    kk = kk+1
                    var(:,:,kk) = fdata(:,:,k)
                  endif
                enddo
              endif
            endif ! diffs
!
! User-specified height scale (interpolate from zp to ht):
! If diffs, do interpolation first, then take diffs.
          else
            count_4d(3) = nhlev
            allocate(var(nlon,nlat,nhlev),stat=iier)
            if (iier /= 0) then
              write(6,"('mkcdf allocate error for var: nlon=',i3,
     |          ' nlat=',i3,' nhlev=',i3)") nlon,nlat,nhlev
              call allocerr(ier,
     |          "mkcdf allocate tmp var with (nlon,nlat,nhlev)")
            endif
!
! Do global height interpolation of field to user defined height scale:
!
            loginterp = 0
            if (trim(f(i)%type)=='DENSITY') loginterp = 1
            call glbhtin(fdata,fz,nlon,nlat,npress,
     |        var,hlev,nhlev,loginterp,spval,ier,1)
!
! If doing diffs, interpolate control field:
            if (diffs) then
              allocate(varcntr(nlon,nlat,nhlev),stat=iier)
              if (iier /= 0) then
                write(6,"('mkcdf allocate error for varcntr: nlon=',i3,
     |            ' nlat=',i3,' nhlev=',i3)") nlon,nlat,nhlev
                call allocerr(ier,
     |            "mkcdf allocate varcntr with (nlon,nlat,nhlev)")
              endif
              call glbhtin(fcdata,fzcntr,nlon,nlat,npress,
     |          varcntr,hlev,nhlev,loginterp,spval,ier,1)
!
! Take diffs of height interpolated fields (leaving diffs in var):
              length = size(var)
              call mkdiffs(var,varcntr,var,length,f(i)%difftype)
            endif ! diffs of ZP fields at heights
          endif ! zp or ht within ZP field type
!
! OH height scale:
        elseif (trim(f(i)%vtype)=='HEIGHT') then
          count_4d(3) = nhlev_oh
          allocate(var(nlon,nlat,nhlev_oh),stat=iier)
          if (iier /= 0) call allocerr(ier,
     |      "mkcdf allocate tmp var with (nlon,nlat,nhlev_oh)")
          var(:,:,:) = fdata(:,:,:)
          if (diffs) then
            allocate(varcntr(nlon,nlat,nhlev_oh),stat=iier)
            if (iier /= 0) call allocerr(ier,
     |          "mkcdf allocate varcntr for with (nlon,nlat,nhlev_oh)")
            varcntr(:,:,:) = fcdata(:,:,:)
            length = size(varcntr)
            call mkdiffs(var,varcntr,var,length,f(i)%difftype)
          endif
!
! Height-independent fields (will be 3d (lon,lat,time) on netcdf file):
        elseif (trim(f(i)%vtype)=='HT-INDEP') then
          start_3d(1:2) = 1
          start_3d(3) = itime
          count_3d(1) = nlon
          count_3d(2) = nlat
          count_3d(3) = 1
          allocate(var2(nlon,nlat),stat=iier)
          if (iier /= 0) then
            write(6,"('>>> Error allocating var2 with nlon=',i3,
     |        ' nlat=',i3)") nlon,nlat
            call allocerr(ier,
     |      "mkcdf allocate tmp var2 with (nlon,nlat)")
          endif
          var2(:,:) = fdata(:,:,1)
          if (diffs) then
            allocate(var2cntr(nlon,nlat),stat=iier)
            if (iier /= 0) call allocerr(ier,
     |          "mkcdf allocate var2cntr for with (nlon,nlat)")
            var2cntr(:,:) = fcdata(:,:,1)
            length = size(var2cntr)
            call mkdiffs(var2,var2cntr,var2,length,f(i)%difftype)
          endif
        endif
!
! Put values into the 4d variables (lon,lat,lev,time):
! (if zonal means, 3d field is (lat,lev,time))
        if (trim(f(i)%vtype)/='HT-INDEP') then ! NOT height-indep
          if (cdf_zonalmean <= 0.) then
            istat = nf_put_vara_real(ncid,idv_flds(i),start_4d,count_4d,
     |        var)
            if (istat /= NF_NOERR) then
              write(char120,"('Error return from nf_put_vara_real',
     |          ' for field var ',a,' itime=',i2)") f(i)%fname8,itime
              call handle_ncerr(istat,char120)
            endif
          else ! zonal means 3d field (lat,lev,time)
            start_3d(1:2) = 1
            start_3d(3) = itime
            count_3d(1) = nlat
            count_3d(2) = count_4d(3) ! vertical
            count_3d(3) = 1
            istat = nf_put_vara_real(ncid,idv_flds(i),start_3d,count_3d,
     |        var(1,:,:))
            if (istat /= NF_NOERR) then
              write(char120,"('Error return from nf_put_vara_real',
     |          ' for zm field var ',a,' itime=',i2)") f(i)%fname8,itime
              call handle_ncerr(istat,char120)
            endif
          endif
!
! Height independent 3d field is (lon,lat,time)
! (if zonal means, 2d field is (lat,time))
        else
          if (cdf_zonalmean <= 0.) then
            istat = nf_put_vara_real(ncid,idv_flds(i),start_3d,count_3d,
     |        var2)
            if (istat /= NF_NOERR) then
              write(char120,"('Error return from nf_put_vara_real',
     |          ' for ht-indep field var2 ',a,' itime=',i2)") 
     |          f(i)%fname8,itime
              call handle_ncerr(istat,char120)
            endif
          else ! height independent zonal means 2-d field (lat,time)
            start_2d(1) = 1
            start_2d(2) = itime
            count_2d(1) = nlat
            count_2d(2) = 1
            istat = nf_put_vara_real(ncid,idv_flds(i),start_2d,count_2d,
     |        var2(1,:))
            if (istat /= NF_NOERR) then
              write(char120,"('Error return from nf_put_vara_real',
     |          ' for ht-indep zonal mean field var2 ',a,' itime=',i2)") 
     |          f(i)%fname8,itime
              call handle_ncerr(istat,char120)
            endif
          endif
        endif ! ht-independent field
      enddo floop2
!
! Release allocated memory:
      if (allocated(fdata))    deallocate(fdata)
      if (allocated(fcdata))   deallocate(fcdata)
      if (allocated(var))      deallocate(var)
      if (allocated(var2))     deallocate(var2)
      if (allocated(fdiffs))   deallocate(fdiffs)
      if (allocated(varcntr))  deallocate(varcntr)
      if (allocated(var2cntr)) deallocate(var2cntr)
      if (ier /= 0) return
      write(6,"('Wrote to netcdf file ',a,' at mtime=',
     |  i3,':',i2,':',i2,' itime=',i3)") trim(nc%filename),h%mtime,itime
!
! If last time iteration, add last_history attribute and close the file:
      if (itime == ntimes) then
!
! Add last history (mssvol and mtime) as global attribute:
! 2/16/06 btf: because h%mssvol may be blank, use disk file h%histfile
        istat = nf_redef(ncid)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |      'Error reentering define mode at last history')
        write(char80,"(a,', ',i3.3,':',i2.2,':',i2.2)") 
     |    trim(h%histfile),h%mtime
!    |    trim(h%mssvol),h%mtime
        istat = nf_put_att_text(ncid,NF_GLOBAL,"last_history",
     |    len_trim(char80),trim(char80))
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_att_text ',
     |      'for last_history global attribute: ',a)") trim(char80)
          call handle_ncerr(istat,char120)
        endif
        istat = nf_enddef(ncid)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'Error return from nf_enddef')
!
! Close netcdf file:
        istat = nf_close(ncid)
        if (istat /= NF_NOERR) then
          call handle_ncerr(istat,
     |      'Error return from nf_close to close netcdf file')
        else
!         write(6,"('Closed netcdf file ',a)") trim(nc%filename)
        endif
      endif
!
! Release local space:
      if (allocated(fz))     deallocate(fz)
      if (allocated(fzcntr)) deallocate(fzcntr)
!
! Debug allocate error:
!     write(6,"('allocated(hlev)    =',l1)") allocated(hlev)
!     write(6,"('allocated(zplev)   =',l1)") allocated(zplev)
!     write(6,"('allocated(hlev_oh) =',l1)") allocated(hlev_oh)
!     write(6,"('allocated(var)     =',l1)") allocated(var)
!     write(6,"('allocated(varcntr) =',l1)") allocated(varcntr)
!     write(6,"('allocated(fdiffs)  =',l1)") allocated(fdiffs)
!     write(6,"('allocated(var2)    =',l1)") allocated(var2)
!     write(6,"('allocated(var2cntr)=',l1)") allocated(var2cntr)
!     write(6,"('allocated(fz)      =',l1)") allocated(fz)
!     write(6,"('allocated(fzcntr)  =',l1)") allocated(fzcntr)
!
      return
      end subroutine mkcdf
!-------------------------------------------------------------------
      subroutine addzplev(ncid,fz,nlon,nlat,nzplev,zprange,dlev,
     |  itime)
      implicit none
!
! Args:
      integer,intent(in) :: ncid,nlon,nlat,nzplev,itime
      real,intent(in) :: fz(nlon,nlat,nzplev),zprange(2),dlev
!
! Local:
      integer :: k,istat,id_zplev,idv_zplev,ids1(1),ids4(4),ier,
     |  start_4d(4),count_4d(4)
      integer,save :: idv_zp
      real,allocatable :: zplev(:)
      character(len=80) :: char80
!
!     write(6,"('Enter addzplev: itime=',i3)") itime
!
! If first time, define dimension zplev, coord var zplev(zplev), 
! and 4-d var ZP(nlon,nlat,nzplev,time) on the netcdf file:
!
      if (itime==1) then 
!
! Put nc file in define mode:
        istat = nf_redef(ncid)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |      'addzplev: Error reentering define mode')
!
! Define zplev dimension on nc file (id_zplev is returned):
        istat = nf_def_dim(ncid,"zplev",nzplev,id_zplev)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining zplev dimension')
!
! Define zplev(zplev) coord array (idv_zplev is defined):
        ids1(1) = id_zplev
        istat = nf_def_var(ncid,"zplev",NF_FLOAT,1,ids1,idv_zplev)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining zplev coord variable')
!
! Long name of zplev coord var:
        istat = nf_put_att_text(ncid,idv_zplev,"long_name",
     |    12,'LOG PRESSURE')
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining zplev long name attribute')
!
! Units of zplev coord var:
        istat = nf_put_att_text(ncid,idv_zplev,"units",
     |    8,'LN(P0/P)')
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining zplev units attribute')
!
! Define ZP(lon,lat,zplev,time) (nc structure is module data)
        ids4(1) = nc%id_lon
        ids4(2) = nc%id_lat
        ids4(3) = id_zplev
        ids4(4) = nc%id_time
        istat = nf_def_var(ncid,'ZP',NF_FLOAT,4,ids4,idv_zp)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining 4-d variable ZP')
!
! Long name of ZP var:
        istat = nf_put_att_text(ncid,idv_zp,"long_name",
     |    19,'GEOPOTENTIAL HEIGHT')
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining ZP long name attribute')
!
! Units of ZP var:
        istat = nf_put_att_text(ncid,idv_zp,"units",
     |    2,'KM')
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining ZP units attribute')
!
! ZP is at interface levels:
        istat = nf_put_att_text(ncid,idv_zp,"zp_interp",
     |    2,'no')
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining ZP zp_interp attribute')
!
! ZP type SCALAR:
        istat = nf_put_att_text(ncid,idv_zp,"type",
     |    6,'SCALAR')
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error defining ZP type attribute')
!
! End define mode
        istat = nf_enddef(ncid)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error return from nf_enddef')
!
! Allocate and define local zplev array:
        allocate(zplev(nzplev),stat=ier) 
        if (ier /= 0) call allocerr(ier,"addzplev allocating zplev")
        do k=1,nzplev
          zplev(k) = zprange(1)+(k-1)*dlev
        enddo
!
! Give values to zplev(zplev) coord array (still 1st time only):
        istat = nf_put_var_real(ncid,idv_zplev,zplev)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'addzplev: Error from nf_put_var_real for zplev')
!
        write(6,"('addzplev: defined dimension zplev, coord var ',
     |    'zplev(zplev), and 4-d geopotential var ZP on netcdf file')")
      endif ! itime==1
!
! This part is called every time (dimension and vars have been defined):
!
! Give values to ZP variable:
      start_4d(1:3) = 1
      start_4d(4) = itime
      count_4d(1) = nlon
      count_4d(2) = nlat
      count_4d(3) = nzplev
      count_4d(4) = 1
      istat = nf_put_vara_real(ncid,idv_zp,start_4d,count_4d,fz)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'addzplev: Error from nf_put_vara_real for ZP (local fz)')

!     write(6,"('addzplev: itime=',i3,' fz 3d min,max=',2e12.4)")
!    |  itime,minval(fz),maxval(fz)

      end subroutine addzplev
!-------------------------------------------------------------------
      end module mkcdf_module
