!
      module netcdf_module
      use fields_module,only: nfprog
      implicit none
!
! Module to create, read and write netcdf history files:
!
      include "netcdf.inc"
!
      integer,private ::
     |  id_time,                   ! dimension id for time (unlimited)
     |  id_mtime,                  ! dimension id for model time (3)
     |  id_lon,id_lat,id_lev,      ! dimension ids for 3d model grid
     |  id_mlon,id_mlat,id_mlev,   ! dimension ids for 3d magnetic grid
     |  id_latlon,                 ! dimension id for coords (2)
     |  id_dtide,                  ! dimension id for dtide (2)
     |  id_sdtide,                 ! dimension id for sdtide (10)
     |  id_fake,                   ! fake dimension id (1)
     |  ids1(1),ids2(2),ids4(4),	! vectors of dim id's
! These next 4 should probably be made local
     |  start_2d(2),count_2d(2),	! start loc and count for 2d vars
     |  start_4d(4),count_4d(4),	! start loc and count for 4d vars
!
! Variable id's (idv_xxxx). Note lon, lat, lev, and time are
!   coordinate variables (lon(lon), lat(lat), lev(lev), time(time)).
! Variables starting at mtime are "primary" variables.
! idv_fprog(nfprog) are id's for prognostic fields.
!
     |  idv_time,idv_mtime,idv_ut,      ! variable id's
     |  idv_lon,idv_lat,idv_lev,
     |  idv_mlon,idv_mlat,idv_mlev,
     |  idv_step,idv_iter,idv_year,idv_day,idv_p0,
     |  idv_hpower,idv_ctpoten,idv_byimf,idv_colfac,
     |  idv_f107d,idv_f107a,idv_mag,idv_dtide,idv_sdtide,
     |  idv_nflds,idv_fprog(nfprog)
      integer,allocatable :: idv_sech(:)
!
      contains
!-----------------------------------------------------------------------
      subroutine nc_open(ncid,flnm,status,rdwr)
!
! Open or create a netcdf dataset flnm, returning netcdf file
!   id handle in ncid.
! If status=='NEW' or 'REPLACE', a new netcdf file is created,
!   otherwise (status=='OLD'), open existing dataset).
! If an existing file is opened and rdwr=='READ', the dataset is
!   opened read-only (NF_NOWRITE), otherwise (rdwr=='WRITE' or 
!   'APPEND') the dataset is opened read/write (NF_WRITE).
! If the open fails, print error msg and return ncid == 0.
!
! Args:
      integer,intent(out) :: ncid    ! netcdf file id
      character(len=*),intent(in) :: flnm,status,rdwr
!
! Local:
      integer :: istat
      character(len=120) :: char120
!
      ncid = 0
      if (status=='NEW'.or.status=='REPLACE') then
!
! Create new dataset (clobber old one if it pre-exists):
!
        istat = nf_create(flnm,NF_CLOBBER,ncid)
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_create for netcdf ',
     |      'file ',a,' status=',a)") trim(flnm),status
          call handle_ncerr(istat,char120)
          ncid = 0
        else
          write(6,"('Created netcdf file ',a,' ncid=',i8)") 
     |      trim(flnm),ncid
        endif
!
! Open pre-existing dataset:
! NF_WRITE opens with read-write, where write means append, change, etc
! NF_NOWRITE is read-only.
!
      elseif (status=='OLD') then
        if (rdwr=='WRITE'.or.rdwr=='APPEND') then
          istat = nf_open(flnm,NF_WRITE,ncid)
        elseif (rdwr=='READ') then
          istat = nf_open(flnm,NF_NOWRITE,ncid)
        else
          write(6,"('>>> nc_open: unrecognized rdwr=',a)") trim(rdwr)
          write(6,"('  Will open NF_WRITE (read/write).')")
          istat = nf_open(flnm,NF_WRITE,ncid)
        endif
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_open for netcdf ',
     |      'file ',a,' status=',a)") trim(flnm),status
          call handle_ncerr(istat,char120)
          ncid = 0
        else
          write(6,"('Opened existing netcdf file ',a,' ncid=',i8)") 
     |      trim(flnm),ncid
        endif
!
! Unknown status request:
      else
        write(6,"(/,'>>> nc_open: unrecognized status = ',a)")
     |    status
        ncid = 0
      endif
      end subroutine nc_open
!-------------------------------------------------------------------
      subroutine nc_close(ncid)
!
! Close a netcdf dataset file:
!
! Arg:
      integer,intent(in) :: ncid
!
! Local:
      integer :: istat
!
      istat = nf_close(ncid)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error return from nf_close')
      end subroutine nc_close
!-----------------------------------------------------------------------
      subroutine nc_rdhist(ncid,diskfile,mtime,itime,ier)
!
! Read history at model time mtime from netcdf file diskfile,
!   currently open under ncid. Define output history in global
!   history structure h (module hist_module in history.f).
!   Prognostic fields are transferred to global fg-array (wrfprog).
! If there is a problem (e.g., mtime not found), return ier > 0
!
      use hist_module,only: h,hist_initype,hist_print
      use fields_module,only: fprog
      use init_module,only: istep
      include "params.h"
      include "grid.h"
!
! Args:
      integer,intent(in) :: ncid,mtime(3)
      character(len=*) :: diskfile
      integer,intent(out) :: itime,ier
!
! Local:
      integer,parameter :: mxdims=20 ! assume <= mxdims dimensions
      integer :: istat,id_time,id_mtime,mtime_read(3),ntimes,
     |  i,ix,istart2(2),icount2(2),j
      integer :: ndims,nvars,ngatts,idunlim,itype,natts,len,
     |  iddims(mxdims),iscalar
      character(len=80) :: varname
      real :: scalar,var1(1),var22(2,2),var2(2),var10(10)
!
      write(6,"('Reading source history from diskfile ',a,':')")
     |  trim(diskfile)
      ier = 0
!
! Get id for mtime variable:
      istat = nf_inq_varid(ncid,"mtime",idv_mtime) 
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,'Error getting id of variable mtime')
        ier = 1
        return
      endif
!
! Get number of histories (times) on the file:
      istat = nf_inq_unlimdim(ncid,id_time) ! id of unlimited record var
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting unlimited dimension id')
      istat = nf_inq_dimlen(ncid,id_time,ntimes) ! length of time var
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error getting length of time record variable')
!
! Search for requested model time. When found, retain time index itime:
      istart2(1) = 1
      icount2(1) = 3
      icount2(2) = 1
      itime = 0
      do i=1,ntimes
        istart2(2) = i
        istat = nf_get_vara_int(ncid,idv_mtime,istart2,icount2,
     |    mtime_read)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'Error getting mtime values')
        write(6,"('nc_rdhist: seeking ',3i4,' found ',3i4,'  n=',i3)")
     |    mtime,mtime_read,i
        if (all(mtime_read==mtime)) then
          itime = i
          exit
        endif
      enddo
      if (itime==0) then
        write(6,"(/,'>>> nc_rdhist: mtime ',3i4,' NOT found on ',
     |    'diskfile ',a)") mtime,trim(diskfile)
        ier = 1
        return
      endif
!
! Define history structure from netcdf vars:
!
! Initialize history structure:
      call hist_initype(h,istep)  ! initialize history structure
!
! Define vars already found above:
      h%ihist = itime
      h%modeltime(1:3) = mtime
      h%modeltime(4) = 0
!
      istat = nf_get_att_text(ncid,NF_GLOBAL,"rundate",h%rundate)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"createdate",h%createdate)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"writedate",h%writedate)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"logname",h%logname)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"host",h%host)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"system",h%system)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"model_name",h%model_name)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"model_version",
     |  h%model_version)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"mss_path",h%mss_path)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"mss_source",h%mss_source)
      istat = nf_get_att_text(ncid,NF_GLOBAL,"history_type",h%type)
!
! Get number of dims, vars, atts, and id of unlimited dimension:
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
!
! Get dimensions and check against current model parameters:
! (i.e., current model run at resolution A should not read a 
!  source history from a model that was run with resolution B) 
!
      do i=1,ndims
        istat = nf_inq_dim(ncid,i,varname,len)
        if (trim(varname)=='lon') then
          if (len /= nlon) then
            write(6,"(/,'>>> nc_rdhist: length of lon read from',
     |        ' history does not match parameter nlon.')")
            write(6,"('  len=',i4,' nlon=',i4)") len,nlon
          endif
        endif
        if (trim(varname)=='lat') then
          if (len /= nlat) then
            write(6,"(/,'>>> nc_rdhist: length of lat read from',
     |        ' history does not match parameter nlat.')")
            write(6,"('  len=',i4,' nlat=',i4)") len,nlat
          endif
        endif
        if (trim(varname)=='lev') then
          if (len /= nlev) then
            write(6,"(/,'>>> nc_rdhist: length of lev read from',
     |        ' history does not match parameter nlev.')")
            write(6,"('  len=',i4,' nlev=',i4)") len,nlev
          endif
        endif
      enddo
!
! Get variables (should not need to get coord vars lat,lon,lev)
      do i=1,nvars
        istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts)
        select case(trim(varname))
          case('lon')
            istat = nc_get_var_real(ncid,i,glon)
          case('lat') 
            istat = nc_get_var_real(ncid,i,glat)
          case('lev')
            istat = nc_get_var_real(ncid,i,plev)
          case('mlon')
            istat = nc_get_var_real(ncid,i,gmlon)
          case('mlat') 
            istat = nc_get_var_real(ncid,i,gmlat)
          case('mlev')
            istat = nc_get_var_real(ncid,i,pmlev)
          case('time')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%time = iscalar
          case('mtime') ! already got it from hist search above
          case('year')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%year = iscalar
          case('day')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%calday = iscalar
          case('iter')
            istat = nf_get_var1_int(ncid,i,itime,iscalar)
            h%iter = iscalar
          case('ut')
            istat = nc_get_var1_real(ncid,i,itime,scalar)
            h%ut = scalar
          case('mag')
            istat = nc_get_var_real(ncid,i,var22)
            h%mag(:,:) = var22(:,:)
          case('dtide')
            istat = nc_get_var_real(ncid,i,var2)
            h%dtide(:) = var2(:)
          case('sdtide')
            istat = nc_get_var_real(ncid,i,var10)
            h%sdtide(:) = var10(:)
          case('f107d')
            istat = nc_get_var_real(ncid,i,var1)
            h%f107d = var1(1)
          case('f107a')
            istat = nc_get_var_real(ncid,i,var1)
            h%f107a = var1(1)
          case('hpower')
            istat = nc_get_var_real(ncid,i,var1)
            h%hpower = var1(1)
          case('ctpoten')
            istat = nc_get_var_real(ncid,i,var1)
            h%ctpoten = var1(1)
          case('byimf')
            istat = nc_get_var_real(ncid,i,var1)
            h%byimf = var1(1)
          case('colfac')
            istat = nc_get_var_real(ncid,i,var1)
            h%colfac = var1(1)
          case('p0')
            istat = nc_get_var_real(ncid,i,var1)
            h%p0 = var1(1)
          case('timestep')
            istat = nf_get_var_int(ncid,i,iscalar)
            h%step = iscalar
          case('nflds')
            istat = nf_get_var_int(ncid,i,iscalar)
            h%nflds = iscalar
          case default
            if (any(fprog%name==varname)) then
              call rdfprog(ncid,varname,itime)
            else
              write(6,"('Note nc_rdhist: unknown variable: ',a)") 
     |          trim(varname)
            endif
        end select
      enddo ! i=1,nvars
!
! Report to stdout:
      call hist_print(h,'READ')
!     
      end subroutine nc_rdhist
!-----------------------------------------------------------------------
      subroutine rdfprog(ncid,name,itime)
!
! Read prognostic field "name" from netcdf file on ncid into fg-array:
!
      use fields_module,only: fprog,fakeflds,findnx,len_name,nfprog
      include "params.h" ! for grid sizes and spval
      include "grid.h"
      include "fgcom.h"
      include "dynphi.h" ! for dynpot
!
! Args:
      character(len=*),intent(in) :: name
      integer,intent(in) :: ncid,itime
!
! Local:
      integer :: i,j,k,ix,istat,idvar,itype,ndims,iddims(4),natts,
     |  iscalar,idimsizes(4),nx,nxk
      character(len=16) :: rdname
      character(len=120) :: char120
      real :: f3d(nlon,nlat,nlev),fmin,fmax
      real :: fakevar(1,1,1)
      character(len=len_name) :: fprog_names(nfprog)
!
! External:
      integer,external :: strloc
!
! Get info about the field:
      istat = nf_inq_varid(ncid,name,idvar)
      istat = nf_inq_var(ncid,idvar,rdname,itype,ndims,iddims,natts)
      if (ndims /= 4) then
        write(6,"(/,'>>> WARNING rdfprog: bad ndims=',i3,
     |    ' (every prognostic should have 4 dimensions)')") ndims
      endif
!
! Get info about the dimensions:
      do i=1,ndims
        istat = nf_inq_dim(ncid,iddims(i),rdname,idimsizes(i))
      enddo 
!
! Check dimension sizes against current model grid resolution:
      if (.not.fakeflds) then
        if (idimsizes(1) /= nlon .or. idimsizes(2) /= nlat .or. 
     |      idimsizes(3) /= nlev) then
          write(6,"(/,'>>> WARNING rdfprog: bad dimension sizes',
     |      ' for prognostic ',a)") trim(rdname)
          write(6,"('  dim sizes=',4i4,' but should be ',
     |      'nlon,nlat,nlev=',4i4)") idimsizes,nlon,nlat,nlev
        endif
      endif
      if (fakeflds) then
        start_4d(1:3) = 1
        start_4d(4) = itime
        count_4d(1) = 1
        count_4d(2) = 1
        count_4d(3) = 1
        count_4d(4) = 1
        istat = nf_get_vara_double(ncid,idvar,start_4d,count_4d,fakevar)
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_get_vara_double',
     +      ' for fake field var ',a,' itime=',i2)") name,itime
          call handle_ncerr(istat,char120)
        endif
        return
      endif
!
! Get index to the field in fprog%name:
!     ix = strloc(fprog%name,size(fprog%name),name)

!     write(6,"('rdfprog: strloc fprog: loc(fprog%name)=',
!    |  i12,' len(fprog%name)=',i4,' name=',a,' loc(name)=',
!    |  i12,' len(name)=',i4)") loc(fprog%name),len(fprog%name),
!    |  trim(name),loc(name),len(name)

      fprog_names(:) = fprog(:)%name
      ix = strloc(fprog_names,nfprog,name)
!
!     write(6,"('Read prognostic field ',a,' dimensions = ',4i4)")
!    |  fprog(ix)%name,idimsizes
!
! Read field into local f3d:
      start_4d(1:3) = 1
      start_4d(4) = itime
      count_4d(1) = nlon
      count_4d(2) = nlat
      count_4d(3) = nlev
      count_4d(4) = 1
      istat = nf_get_vara_double(ncid,idvar,start_4d,count_4d,f3d)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_get_vara_double',
     +    ' for field var ',a,' itime=',i2)") name,itime
        call handle_ncerr(istat,char120)
      endif
!
! In old cray-blocked histories, NOPNM was not defined at the
! top level. Define it here with spval for compatability with
! old sources (s.a., read_oldsrc). This is not strictly necessary,
! but is convenient for debug and comparisons between old and new 
! histories with fminmax, etc.
!
      if (trim(name)=='OP_NM') then
        f3d(:,:,nlev) = spval
      endif
!
!     call fminmaxspv(f3d,nlon*nlat*nlev,fmin,fmax,spval)
!     write(6,"('rdfprog: Read field ',a,' itime=',i2,
!    |  ' 3d min,max=',2e12.4)") name(1:8),itime,fmin,fmax
!
! Load local f3d into global fg-array (lons 3-74) and add periodic 
! points:
      nx = findnx(name)
      if (nx < 0) write(6,"(/,'>>> rdfprog: bad return from findnx',
     |  ' for prognostic field ',a)") trim(name)
      do j=1,nlat
        nxk = nx
        do k=1,nlev
          nxk = nxk+1
          do i=1,nlon
            fg(i+2,nxk,j,ixtimep) = f3d(i,j,k) 
          enddo
! periodic points:
          fg(1,nxk,j,ixtimep) = fg(nlon+1,nxk,j,ixtimep)
          fg(2,nxk,j,ixtimep) = fg(nlon+2,nxk,j,ixtimep)
          fg(nlon+3,nxk,j,ixtimep) = fg(3,nxk,j,ixtimep)
          fg(nlon+4,nxk,j,ixtimep) = fg(4,nxk,j,ixtimep)
        enddo
      enddo
!
! Define dynpot:
      if (trim(name)=='POTEN') then
        do j=1,nlat
          nxk = nx
          do k = 1,nlev
            nxk = nxk+1
            do i = 1,nlon+1
              dynpot(i,j,k) = fg(i+2,nxk,j,ixtimep)
            enddo
          enddo
        enddo
      endif
!
      end subroutine rdfprog
!-----------------------------------------------------------------------
      subroutine nc_define(ncid)
!
! Define dimensions, variables and attributes in a new history
!   netcdf dataset (file has been opened, but no histories written yet).
! On input:
!   ncid = file id for open netcdf dataset (in define mode).
!   h = defined history structure (in hist_module)
! On output:
!   Variables and dimensions for the netcdf file are defined.
!   Values are given to the dimension variables, and the file
!     is taken out of define mode (ready to receive histories)
!
      use hist_module,only: h,nsource
      use fields_module,only: nfprog,fprog,fakeflds
      include "params.h"
      include "grid.h"
      include "consts.h" ! for mag grid
!
! Args:
      integer,intent(in) :: ncid
!
! Local:
      integer :: i,istat,idum,ivar1(1)
      character(len=120) :: char120
      character(len=80) :: char80
      real :: var1(1),rtdeg
      character(len=8) :: createdate,createtime
!
      rtdeg = 180./(4.*atan(1.))
!
! Define dimensions (field vars will be (time,lev,lat,lon) in
! CDL notation. Mag fields will be (time,mlev,mlat,mlon): 
!
      istat = nf_def_dim(ncid,"time",NF_UNLIMITED,id_time)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining time dimension')
!
! Geographic grid dimensions:
      istat = nf_def_dim(ncid,"lon",nlon,id_lon)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining longitude dimension')
      istat = nf_def_dim(ncid,"lat",nlat,id_lat)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'Error defining latitude dimension')
      istat = nf_def_dim(ncid,"lev",nlev,id_lev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining levels dimension')
!
! Magnetic grid dimensions:
      istat = nf_def_dim(ncid,"mlon",imaxmp,id_mlon)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining magnetic longitude dimension')
      istat = nf_def_dim(ncid,"mlat",jmaxm,id_mlat)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'Error defining magnetic latitude dimension')
      istat = nf_def_dim(ncid,"mlev",zkmxp+3,id_mlev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'Error defining magnetic levels dimension')
!
      istat = nf_def_dim(ncid,"mtime",3,id_mtime)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining mtime dimension')
      istat = nf_def_dim(ncid,"latlon",2,id_latlon)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining latlon dimension')
      istat = nf_def_dim(ncid,"dtidedim",2,id_dtide)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining dtide dimension')
      istat = nf_def_dim(ncid,"sdtidedim",10,id_sdtide)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining sdtide dimension')
      if (fakeflds) then ! fake dimension for testing
        istat = nf_def_dim(ncid,"fakedim",1,id_fake)
        if (istat /= NF_NOERR) call handle_ncerr(istat,
     +    'Error defining fakedim dimension')
      endif
!
! Define dimension (coordinate) variables and attributes:
!
! Geographic longitude (deg) (coordinate variable lon(lon)):
      ids1(1) = id_lon
      istat = nf_def_var(ncid,"lon",NF_DOUBLE,1,ids1,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,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,idv_lon,"units",7,'degrees')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of longitude dimension variable')
!
! Geographic latitude (deg) (coordinate variable lat(lat)):
      ids1(1) = id_lat
      istat = nf_def_var(ncid,"lat",NF_DOUBLE,1,ids1,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,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,idv_lat,"units",7,'degrees')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of latitude dimension variable')
!
! Vertical levels (log pressure) (coordinate variable lev(lev)):
      ids1(1) = id_lev
      istat = nf_def_var(ncid,"lev",NF_DOUBLE,1,ids1,idv_lev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining levels dimension variable')
      write(char80,"('log pressure levels')")
      istat = nf_put_att_text(ncid,idv_lev,"long_name",
     +  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of levels dimension variable')
      istat = nf_put_att_text(ncid,idv_lev,"units",8,'ln(p0/p)')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of levels dimension variable')
!
! Magnetic longitude (deg) (coordinate variable mlon(mlon)):
      ids1(1) = id_mlon
      istat = nf_def_var(ncid,"mlon",NF_DOUBLE,1,ids1,idv_mlon)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining magnetic longitude dimension variable')
      write(char80,"('magnetic longitude (-west, +east)')")
      istat = nf_put_att_text(ncid,idv_mlon,"long_name",
     +  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of mag longitude dimension variable')
      istat = nf_put_att_text(ncid,idv_lon,"units",7,'degrees')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of magnetic longitude dimension variable')
!
! Magnetic latitude (deg) (coordinate variable mlat(mlat)):
      ids1(1) = id_mlat
      istat = nf_def_var(ncid,"mlat",NF_DOUBLE,1,ids1,idv_mlat)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining magnetic latitude dimension variable')
      write(char80,"('magnetic latitude (-south, +north)')")
      istat = nf_put_att_text(ncid,idv_mlat,"long_name",
     +  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of mag latitude dimension variable')
      istat = nf_put_att_text(ncid,idv_mlat,"units",7,'degrees')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of magnetic latitude dimension variable')
!
! Magnetic levels (log pressure) (coordinate variable mlev(mlev)):
      ids1(1) = id_mlev
      istat = nf_def_var(ncid,"mlev",NF_DOUBLE,1,ids1,idv_mlev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining magnetic levels dimension variable')
      write(char80,"('log pressure levels (magnetic grid)')")
      istat = nf_put_att_text(ncid,idv_mlev,"long_name",
     +  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of mag levels dimension variable')
      istat = nf_put_att_text(ncid,idv_mlev,"units",8,'ln(p0/p)')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of magnetic levels dimension variable')
!
! ids1(1) is id of time dimension for the following defines:
      ids1(1) = id_time	! for 1d time-unlimited vars
!
! Time (total model time in minutes) (coordinate variable time(time)):
      istat = nf_def_var(ncid,"time",NF_INT,1,ids1,idv_time)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining time dimension variable')
      write(char80,"('total model time (day*24*60+hour*60+min)')")
      istat = nf_put_att_text(ncid,idv_time,"long_name",
     +  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of time dimension variable')
      istat = nf_put_att_text(ncid,idv_time,"units",7,'minutes')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of time dimension variable')
!
! Define "primary" (non-coordinate) time dependent variables:
! (these variables get a value for each history) 
!
! Model time variable mtime(3,time), where 3 is the mtime dimension 
! (id_mtime) (day,hour,minute), and time is the unlimited dimension 
! (id_time):
!
      ids2(1) = id_mtime
      ids2(2) = id_time
      istat = nf_def_var(ncid,"mtime",NF_INT,2,ids2,idv_mtime)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining mtime variable')
      write(char80,"('model times (day, hour, minute)')")
      istat = nf_put_att_text(ncid,idv_mtime,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name attribute of mtime variable')
      istat = nf_put_att_text(ncid,idv_mtime,"units",17,
     |  'day, hour, minute')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of mtime variable')  
!
! Calendar year and day (this will be advanced if model is advancing
! in calendar time, otherwise will be constant from input)
      istat = nf_def_var(ncid,"year",NF_INT,1,ids1,idv_year)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining calendar year variable')
      istat = nf_put_att_text(ncid,idv_year,"long_name",
     |  13,'calendar year')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of calendar year variable')
!
      istat = nf_def_var(ncid,"day",NF_INT,1,ids1,idv_day)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining calendar day variable')
      istat = nf_put_att_text(ncid,idv_day,"long_name",
     |  12,'calendar day')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of calendar day variable')
!
! Number of time steps from model time 0,0,0 to current model time
! (iter(time)):
      istat = nf_def_var(ncid,"iter",NF_INT,1,ids1,idv_iter)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining iter variable')
      istat = nf_put_att_text(ncid,idv_iter,"long_name",
     |  42,'number of time steps from model time 0,0,0')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of iter variable')
!
! Decimal ut(time) (ids1==id_time is unlimited time dimension):
      call defvar_time_dbl(ncid,"ut",ids1,idv_ut,
     |  "universal time (from model time hour and minute)",
     |  "hours")
!
! Time step (seconds) (integer scalar) (time-independent):
      istat = nf_def_var(ncid,"timestep",NF_INT,0,idum,idv_step)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining timestep variable')
      istat = nf_put_att_text(ncid,idv_step,"long_name",8,'timestep')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name of timestep variable')
      istat = nf_put_att_text(ncid,idv_step,"units",7,'seconds')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of timestep variable')
!
! f107 daily and f107 average solar fluxes:
      call defvar_time_dbl(ncid,"f107d",ids1,idv_f107d,
     |  "Daily 10.7 cm solar flux","W/m^2Hz*1.e-22")
      call defvar_time_dbl(ncid,"f107a",ids1,idv_f107a,
     |  "81-day average 10.7 cm solar flux","W/m^2Hz*1.e-22")
!
! Hemispheric power: 
      call defvar_time_dbl(ncid,"hpower",ids1,idv_hpower,
     |  "Hemispheric power","Gw")
!
! Cross-tail potential:
      call defvar_time_dbl(ncid,"ctpoten",ids1,idv_ctpoten,
     |  "Cross-tail potential","Volts")
!
! Byimf:
      call defvar_time_dbl(ncid,"byimf",ids1,idv_byimf,
     |  "BY component of IMF"," ")
!
! Time-independent "primary" variables:
!
! Magnetic pole coordinates (real(2,2)):
      ids2(1) = id_latlon
      ids2(2) = id_latlon
      istat = nf_def_var(ncid,"mag",NF_DOUBLE,2,ids2,idv_mag)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining variable mag')
      write(char80,"('lat,lon coordinates of S,N magnetic poles')")
      istat = nf_put_att_text(ncid,idv_mag,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name attribute of variable mag')
      istat = nf_put_att_text(ncid,idv_mtime,"units",7,
     |  'degrees')
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining units of variable mag')  
!
! Diurnal tide amp/phase:
      ids1(1) = id_dtide
      istat = nf_def_var(ncid,"dtide",NF_DOUBLE,1,ids1,idv_dtide)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining variable dtide')
      write(char80,"('amplitude and phase of diurnal tide')")
      istat = nf_put_att_text(ncid,idv_dtide,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name attribute of variable dtide')
!
! Semi-diurnal tide amps/phases:
      ids1(1) = id_sdtide
      istat = nf_def_var(ncid,"sdtide",NF_DOUBLE,1,ids1,idv_sdtide)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining variable sdtide')
      write(char80,"('amplitude and phase of semi-diurnal tide')")
      istat = nf_put_att_text(ncid,idv_sdtide,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name attribute of variable sdtide')
!
! Collision factor:
      istat = nf_def_var(ncid,"colfac",NF_DOUBLE,0,idum,idv_colfac)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining variable colfac')
      istat = nf_put_att_text(ncid,idv_colfac,"long_name",
     |  16,"Collision factor")
!
! Standard pressure:
      istat = nf_def_var(ncid,"p0",NF_DOUBLE,0,idum,idv_p0)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining variable p0')
      istat = nf_put_att_text(ncid,idv_p0,"long_name",
     |  17,"Standard pressure")
!
! Number of fields (will be nfprog for primary histories):
      istat = nf_def_var(ncid,"nflds",NF_INT,0,idum,idv_nflds)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining nflds variable')
      write(char80,"('number of 3-d model fields')")
      istat = nf_put_att_text(ncid,idv_nflds,"long_name",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'Error defining long_name attribute of variable nflds')
!
! All prognostic fields are on primary history:
      if (h%type(1:3)=='pri') then
        do i=1,nfprog
          call defvar_fprog(ncid,fprog(i),idv_fprog(i))
        enddo
!
! Define prognostics and diagnostic (geographic and magnetic) 
!   fields on secondary history:
      elseif (h%type(1:3)=='sec') then
        call defvar_sech(ncid,h%nfgeo,h%nfmag)
      else
        write(6,"('>>> nc_define: unknown h%type=',a)") h%type
      endif
!
! Define global file attributes:
!
! Run date and time:
      istat = nf_put_att_text(ncid,NF_GLOBAL,"rundate",
     +  len_trim(h%rundate),trim(h%rundate))
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     +    'for rundate global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Creation date:
      call datetime(createdate,createtime)
      h%createdate = createdate//' '//createtime
      istat = nf_put_att_text(ncid,NF_GLOBAL,"createdate",
     +  len_trim(h%createdate),trim(h%createdate))
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     +    'for createdate global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Most recent write date: 
! (This global attribute is blank now, but is updated by nc_wrhist,
!  when histories are written. h%writedate is updated by define_hist 
!  in output.f when histories are written)
!
      istat = nf_put_att_text(ncid,NF_GLOBAL,"writedate",
     +  len_trim(h%writedate),trim(h%writedate))
      if (istat /= NF_NOERR) then
        write(char120,"('nc_define: Error return from nf_put_att_text ',
     +    'for writedate global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! User login name:
      istat = nf_put_att_text(ncid,NF_GLOBAL,"logname",
     +  len_trim(h%logname),trim(h%logname))
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     +    'for logname global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Machine host name:
      istat = nf_put_att_text(ncid,NF_GLOBAL,"host",
     +  len_trim(h%host),trim(h%host)) 
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     +    'for host global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Machine operating system (from pre-proc macro)
      istat = nf_put_att_text(ncid,NF_GLOBAL,"system",
     +  len_trim(h%system),trim(h%system)) 
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     +    'for system global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Model name string:
      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,"('Error return from nf_put_att_text ',
     +    'for model name global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Model version string:
      istat = nf_put_att_text(ncid,NF_GLOBAL,"model_version",
     +  len_trim(h%model_version),trim(h%model_version)) 
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     +    'for model version global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Mss path to history file:
      istat = nf_put_att_text(ncid,NF_GLOBAL,"mss_path",
     |  len_trim(h%mss_path),trim(h%mss_path))
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     |    'for mss_path global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Primary or Secondary histories:
      if (h%type(1:3)=='pri') then
        istat = nf_put_att_text(ncid,NF_GLOBAL,"history_type",
     |    7,'primary')
      elseif (h%type(1:3)=='sec') then
        istat = nf_put_att_text(ncid,NF_GLOBAL,"history_type",
     |    9,'secondary')
      else
        write(6,"('>>> nc_define: unknown h%type=',a)") h%type
        istat = nf_put_att_text(ncid,NF_GLOBAL,"history_type",
     |    7,'unknown')
      endif
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     |    'for history_type global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Mss path to source history file:
      if (nsource > 0) then ! initial source file
        istat = nf_put_att_text(ncid,NF_GLOBAL,"mss_source",
     |    len_trim(h%mss_source)+10,trim(h%mss_source)//' (initial)')
      else ! continuation source file
        istat = nf_put_att_text(ncid,NF_GLOBAL,"mss_source",
     |    len_trim(h%mss_source)+15,
     |    trim(h%mss_source)//' (continuation)')
      endif
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     |    'for mss_source global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Missing value (nc_put_att_real is in util.F):
      var1(1) = spval
      istat = nc_put_att_real(ncid,NF_GLOBAL,"missing_value",1,var1)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nc_put_att_real ',
     |    'for missing_value global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Conversion from lev (ln(p0/p)) to pressure mb:
      write(char80,"('p0*e(-lev(k))*1.e-3')")
      istat = nf_put_att_text(ncid,NF_GLOBAL,"lev_to_mb",
     |  len_trim(char80),trim(char80))
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_att_text ',
     |    'for lev_to_mb global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Number of histories on the file:
      ivar1(1) = h%ihist ! this will be updated in nc_wrhist
      istat = nf_put_att_int(ncid,NF_GLOBAL,"nhist",NF_INT,1,ivar1)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nc_put_att_int ',
     |    'for nhist global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Delta minutes between histories on the file:
      ivar1(1) = h%delhmins
      istat = nf_put_att_int(ncid,NF_GLOBAL,"delhist_mins",NF_INT,
     |  1,ivar1)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nc_put_att_int ',
     |    'for delhist_mins global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Exit define mode:
      istat = nf_enddef(ncid)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'nc_define: Error return from nf_enddef')
!
! Define grid vars (this will be done earlier in the real model):
      do i=1,nlon
        glon(i) = glon1+(i-1)*dlon
      enddo
      do i=1,nlat
        glat(i) = glat1+(i-1)*dlat
      enddo
      do i=1,nlev
        plev(i) = plev1+(i-1)*dlev
      enddo
!
! Mag grid coordinates (consts.h) were set in cnstnt.f 
! (convert to degrees):
      do i=1,nmlon
        gmlon(i) = ylonm(i)*rtdeg
      enddo
      do i=1,nmlat
        gmlat(i) = ylatm(i)*rtdeg
      enddo
      do i=1,nlev
        pmlev(i+3) = plev(i)
      enddo
      do i=3,1,-1
        pmlev(i) = pmlev(i+1)-dlev
      enddo
!
! Give values to dimension variables:
!
! Geographic:
      istat = nc_put_var_real(ncid,idv_lon,glon)
      istat = nc_put_var_real(ncid,idv_lat,glat)
      istat = nc_put_var_real(ncid,idv_lev,plev)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nc_put_var_real ',
     +    ' to assign values to geographic dimension vars')")
        call handle_ncerr(istat,char120)
      endif
!
! Magnetic:
      istat = nc_put_var_real(ncid,idv_mlon,gmlon)
      istat = nc_put_var_real(ncid,idv_mlat,gmlat)
      istat = nc_put_var_real(ncid,idv_mlev,pmlev)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nc_put_var_real ',
     +    ' to assign values to magnetic dimension vars')")
        call handle_ncerr(istat,char120)
      endif
      end subroutine nc_define
!-------------------------------------------------------------------
      subroutine defvar_time_dbl(ncid,name,idtime,idvar,
     |  longname,units)

!     call defvar_time_dbl(ncid,"ut",ids1,idv_ut,
!    |  "universal time (from model time hour and minute)",
!    |  "hours")
!
! Define time-dependent variable "name" on open netcdf file ncid.
! Variable has a single dimension in time (e.g., ut(time)).
! Also define long_name and units attributes.
!
! Args:
      integer,intent(in)  :: ncid,idtime(1)
      integer,intent(out) :: idvar
      character(len=*),intent(in) :: name,longname,units
!
! Local:
      integer :: istat
      character(len=80) :: char80
!
! Define the variable:
      istat = nf_def_var(ncid,name,NF_DOUBLE,1,idtime,idvar)
      if (istat /= NF_NOERR) then
        write(char80,"('Error defining double variable ',a)") name
        call handle_ncerr(istat,trim(char80))
      endif
!
! Give long name attribute:
      if (len_trim(longname) > 0) then
        istat = nf_put_att_text(ncid,idvar,"long_name",
     |    len_trim(longname),longname)
        if (istat /= NF_NOERR) then
          write(char80,"('Error defining long_name of double variable ',
     |      a)") name
          call handle_ncerr(istat,trim(char80))
        endif
      endif
!
! Give units attribute:
      if (len_trim(units) > 0) then
        istat = nf_put_att_text(ncid,idvar,"units",len_trim(units),
     |    units)
        if (istat /= NF_NOERR) then
          write(char80,"('Error defining units of variable ',a)") name
          call handle_ncerr(istat,trim(char80))
        endif
      endif
      end subroutine defvar_time_dbl
!-------------------------------------------------------------------
      subroutine defvar_sech(ncid,nfgeo,nfmag)
!
! Define secondary history fields on currently open ncid
! (not data -- see wrfdiag for data write)
!
      use fields_module,only: fprog,fsech,fsechmag,nfprog,len_name
!
! Args:
      integer,intent(in) :: 
     |  ncid,  ! netcdf file id
     |  nfgeo, ! number of fields on geographic grid
     |  nfmag  ! number of fields on magnetic grid
!
! Local:
      integer :: i,ii,ier,istat,iprog
      character(len=80) :: char80
      character(len=len_name) :: fprog_names(nfprog)
!
! External:
       integer,external :: strloc
!
      if (allocated(idv_sech)) deallocate(idv_sech)
      allocate(idv_sech(nfgeo+nfmag),stat=ier)
      if (ier /= 0) then
        write(6,"(/,'>>> WARNING defvar_sech: ',
     |  'error return from allocate for idv_sech: nfgeo=',i3,
     |  ' nfmag=',i3,' ier=',i4)") nfgeo,nfmag,ier
      else
        write(6,"('defvar_sech allocated idv_sech: nfgeo=',i3,
     |    ' nfmag=',i3)") nfgeo,nfmag
      endif
      if (nfgeo > 0) then
        fprog_names(:) = fprog(:)%name
        do i=1,nfgeo

!     write(6,"('defvar_sech: i=',i3,' strloc fprog: loc(fprog%name)=',
!    |  i12,' len(fprog%name)=',i4,' fsech(i)%name=',a,
!    |  ' loc(fsech(i)%name)=',i12,' len(fsech(i)%name)=',i4)") 
!    |  i,loc(fprog%name),len(fprog%name),trim(fsech(i)%name),
!    |  loc(fsech(i)%name),len(fsech(i)%name)

          iprog = strloc(fprog_names,size(fprog),fsech(i)%name)
!
! Define prognostic secondary history field:
          if (iprog > 0) then
            call defvar_fprog(ncid,fprog(iprog),idv_sech(i))
!
! Define diagnostic secondary history field:
          else
            call defvar_fdiag(ncid,fsech(i),idv_sech(i))
          endif
        enddo
      endif
      if (nfmag > 0) then
        do i=1,nfmag
          call defvar_fdiag(ncid,fsechmag(i),idv_sech(nfgeo+i))
        enddo
      endif
      end subroutine defvar_sech
!-------------------------------------------------------------------
      subroutine defvar_fprog(ncid,f,idvar)
!
! Define a prognostic variable on current netcdf history file:
! (not data -- see wrfprog for data write)
!
      use fields_module,only: field,fakeflds
!
! Args:
      integer,intent(in)  :: ncid
      integer,intent(out) :: idvar
      type(field),intent(in) :: f
!
! Local:
      integer :: istat,idimids(4)
      character(len=80) :: char80
!
! Prognostics are dimensioned (lon,lat,lev,time):
      if (.not.fakeflds) then
        idimids(1) = id_lon
        idimids(2) = id_lat
        idimids(3) = id_lev
        idimids(4) = id_time
      else
        idimids(1:3) = id_fake
        idimids(4) = id_time
      endif
! This call returns the variable id, idvar.
      istat = nf_def_var(ncid,f%name,NF_DOUBLE,4,idimids,idvar)
      if (istat /= NF_NOERR) then
        write(char80,"('Error defining prognostic field ',a)") 
     |    f%name
        call handle_ncerr(istat,trim(char80))
      endif
!
! Add long_name attribute:
      istat = nf_put_att_text(ncid,idvar,"long_name",
     |  len_trim(f%long_name),trim(f%long_name))
      if (istat /= NF_NOERR) then
        write(char80,"('Error defining long_name of prognostic',
     |    ' variable ',a,': long_name = ',/,a)") trim(f%name),
     |    trim(f%long_name)
        call handle_ncerr(istat,trim(char80))
      endif
!
! Add units attribute:
      istat = nf_put_att_text(ncid,idvar,"units",len_trim(f%units),
     |  f%units)
      if (istat /= NF_NOERR) then
        write(char80,"('Error defining units of prognostic',
     |    ' variable ',a,': units = ',a)") trim(f%name),trim(f%units)
        call handle_ncerr(istat,trim(char80))
      endif
      end subroutine defvar_fprog
!-------------------------------------------------------------------
      subroutine defvar_fdiag(ncid,f,idvar)
!
! Define a diagnostic variable (geographic or magnetic) on the current 
! netcdf history file:
!
      use fields_module,only: field,fakeflds
!     use init_module,only: istep
!
! Args:
      integer,intent(in)  :: ncid
      integer,intent(out) :: idvar
      type(field),intent(inout) :: f
!
! Local:
      integer :: istat,idimids(4)
      character(len=80) :: char80
!
      if (.not.fakeflds) then
        if (.not.f%magnetic) then ! field is on geographic grid
          idimids(1) = id_lon
          idimids(2) = id_lat
          idimids(3) = id_lev
        else                      ! field is on magnetic grid
          idimids(1) = id_mlon
          idimids(2) = id_mlat
          idimids(3) = id_mlev
        endif
          idimids(4) = id_time
      else
        idimids(1:3) = id_fake
        idimids(4) = id_time
      endif
!
      istat = nf_def_var(ncid,f%name,NF_DOUBLE,4,idimids,idvar)
      if (istat /= NF_NOERR) then
        write(char80,"('Error defining diagnostic field ',a)") 
     |    f%name
        call handle_ncerr(istat,trim(char80))
!     else
!       write(6,"('defvar_fdiag: defined diagnostic field ',a)")
!    |    f%name
      endif
      end subroutine defvar_fdiag
!-------------------------------------------------------------------
      subroutine nc_wrhist(ncid,iprint)
!
! Write current history structure to open netcdf output file:
!
      use hist_module,only: h,hist_print
      use fields_module,only: nfprog,fprog,fsech,fsechmag,fakeflds,
     |  len_name,nfprog
!     include "fgcom.h"  ! temp for fgmnmx
!
! Args:
      integer,intent(in) :: ncid,iprint
!
! Local:
      integer :: mins,istat,ivar1(1),i,j,iprog
      character(len=120) :: char120
      integer :: ndims,nvars,ngatts,unlimdimid
      real :: var1(1),var22(2,2),var2(2),var10(10)
      character(len=len_name) :: fprog_names(nfprog)
!
! External:
      integer,external :: strloc
!
! Total model time in minutes:
      mins = h%modeltime(1)*60*24+h%modeltime(2)*60+h%modeltime(3)
      if (idv_time <= 0) istat = nf_inq_varid(ncid,"time",idv_time)
      istat = nf_put_var1_int(ncid,idv_time,h%ihist,mins)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     +    'time at h%ihist=',i3,': mins=',i10)") h%ihist,mins
        call handle_ncerr(istat,char120)
      endif
!
! Model time (day,hour,minute):
      start_2d(1) = 1
      start_2d(2) = h%ihist
      count_2d(1) = 3
      count_2d(2) = 1
      if (idv_mtime <= 0) istat = nf_inq_varid(ncid,"mtime",idv_mtime)
      istat = nf_put_vara_int(ncid,idv_mtime,start_2d,count_2d,
     |  h%modeltime(1:3))
      if (istat /= NF_NOERR) then
        write(char120,"('Error defining mtime at h%ihist=',i2,
     +    ' h%modeltime=',4i4,' istat=',i5)") h%ihist,h%modeltime,istat
        call handle_ncerr(istat,char120)
      endif
!
! Calendar year (initially from input, then updated if model is
! advancing calendar time):
      if (idv_year <= 0) istat = nf_inq_varid(ncid,"year",idv_year)
      istat = nf_put_var1_int(ncid,idv_year,h%ihist,h%year)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'calendar year at h%ihist=',i3,': h%year=',i5)") 
     |    h%ihist,h%year
        call handle_ncerr(istat,char120)
      endif
!
! Calendar day (initially from input, then updated if model is
! advancing calendar time):
      if (idv_day <= 0) istat = nf_inq_varid(ncid,"day",idv_day)
      istat = nf_put_var1_int(ncid,idv_day,h%ihist,h%calday)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     |    'calendar day at h%ihist=',i3,': h%calday=',i5)") 
     |    h%ihist,h%calday
        call handle_ncerr(istat,char120)
      endif
!
! Iteration number (number of steps from 0,0,0 to current mtime):
      if (idv_iter <= 0) istat = nf_inq_varid(ncid,"iter",idv_iter)
      istat = nf_put_var1_int(ncid,idv_iter,h%ihist,h%iter)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     +    'iter at h%ihist=',i3,': h%iter=',i8)") h%ihist,h%iter
        call handle_ncerr(istat,char120)
      endif
!
! Universal time (decimal hours):
      if (idv_ut <= 0) istat = nf_inq_varid(ncid,"ut",idv_ut)
      istat = nc_put_var1_real(ncid,idv_ut,h%ihist,h%ut,'ut')
!
! Magnetic pole coords:
      if (idv_mag <= 0) istat = nf_inq_varid(ncid,"mag",idv_mag)
      var22(:,:) = h%mag(:,:)
      istat = nc_put_var_real(ncid,idv_mag,var22)
!
! Diurnal tide:
      if (idv_dtide <= 0) istat = nf_inq_varid(ncid,"dtide",idv_dtide)
      var2(:) = h%dtide(:)
      istat = nc_put_var_real(ncid,idv_dtide,var2)
!
! Semi-diurnal tide:
      if (idv_sdtide <= 0) istat = nf_inq_varid(ncid,"sdtide",
     |  idv_sdtide)
      var10(:) = h%sdtide(:)
      istat = nc_put_var_real(ncid,idv_sdtide,var10)
!
! f107d and f107a (daily and average solar fluxes):
      if (idv_f107d <= 0) istat = nf_inq_varid(ncid,"f107d",
     |  idv_f107d)
      istat = nc_put_var1_real(ncid,idv_f107d,h%ihist,h%f107d,'f107d')
      if (idv_f107a <= 0) istat = nf_inq_varid(ncid,"f107a",
     |  idv_f107a)
      istat = nc_put_var1_real(ncid,idv_f107a,h%ihist,h%f107a,'f107a')
!
! Hemispheric power:
      if (idv_hpower <= 0) istat = nf_inq_varid(ncid,"hpower",
     |  idv_hpower)
      istat = nc_put_var1_real(ncid,idv_hpower,h%ihist,h%hpower,
     |  'hpower')
!
! Cross-tail potential:
      if (idv_ctpoten <= 0) istat = nf_inq_varid(ncid,"ctpoten",
     |  idv_ctpoten)
      istat = nc_put_var1_real(ncid,idv_ctpoten,h%ihist,h%ctpoten,
     |  'ctpoten')
!
! BY imf:
      if (idv_byimf <= 0) istat = nf_inq_varid(ncid,"byimf",idv_byimf)
      istat = nc_put_var1_real(ncid,idv_byimf,h%ihist,h%byimf,
     |  'byimf')
!
! Collision factor:
      if (idv_colfac <= 0) istat = nf_inq_varid(ncid,"colfac",
     |  idv_colfac)
      istat = nc_put_var1_real(ncid,idv_colfac,h%ihist,h%colfac,
     |  'colfac')
!
! Standard pressure:
      if (idv_p0 <= 0) istat = nf_inq_varid(ncid,"p0",idv_p0)
      istat = nc_put_var1_real(ncid,idv_p0,h%ihist,h%p0,'p0')
!
! Timestep (seconds):
      if (idv_step <= 0) istat = nf_inq_varid(ncid,"timestep",idv_step)
      istat = nf_put_var1_int(ncid,idv_step,h%ihist,h%step)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var1_int defining ',
     +    'ut at h%ihist=',i3,': h%step=',i4)") h%ihist,h%step
        call handle_ncerr(istat,char120)
      endif
!
! Update number of histories on the file (global attribute):
      ivar1(1) = h%ihist
      istat = nf_put_att_int(ncid,NF_GLOBAL,"nhist",NF_INT,1,ivar1)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nc_put_att_int ',
     +    'for nhist global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Update most recent write date (global attribute):
      istat = nf_put_att_text(ncid,NF_GLOBAL,"writedate",
     +  len_trim(h%writedate),trim(h%writedate))
      if (istat /= NF_NOERR) then
        write(char120,"('nc_wrhist: Error return from nf_put_att_text ',
     +    'for writedate global attribute')")
        call handle_ncerr(istat,char120)
      endif
!
! Number of fields (scalar):
      if (idv_nflds <= 0) istat = nf_inq_varid(ncid,"nflds",idv_nflds)
!     istat = nf_put_var1_int(ncid,idv_nflds,h%ihist,h%nflds)
      istat = nf_put_var_int(ncid,idv_nflds,h%nflds)
      if (istat /= NF_NOERR) then
        write(char120,"('Error from nf_put_var_int defining ',
     +    'nflds: h%nflds=',i4)") h%nflds
        call handle_ncerr(istat,char120)
      endif
!
! Add fields data to primary history:
! Currently (2/00), only prognostic fields are written to the
!   primary histories:
!
      if (h%type(1:3)=='pri') then
        do i=1,nfprog
          call wrfprog(ncid,fprog(i)%name,h%ihist,fakeflds)
        enddo
!
! Add fields data to secondary history:
! If a prognostic field, write field from fg to the history 
!   (sub wrfprog).
! If a diagnostic field, write field from fsech%data
!   (fsech%data was init to spval by sub set_fsech. If istep > 0,
!    then fsech%data should have been defined by addfsech.
!
      elseif (h%type(1:3)=='sec') then
        if (h%nfgeo > 0) then
          fprog_names(:) = fprog(:)%name
          do i=1,h%nfgeo ! geographic fields

!     write(6,"('ncwrhist: i=',i3,' strloc fprog: loc(fprog%name)=',
!    |  i12,' len(fprog%name)=',i4,' fsech(i)%name=',a,
!    |  ' loc(fsech(i)%name)=',i12,' len(fsech(i)%name)=',i4)") 
!    |  i,loc(fprog%name),len(fprog%name),trim(fsech(i)%name),
!    |  loc(fsech(i)%name),len(fsech(i)%name)

            iprog = strloc(fprog_names,size(fprog),fsech(i)%name)
            if (iprog > 0) then
              call wrfprog(ncid,fsech(i)%name,h%ihist,fakeflds)
            else
              call wrfdiag(ncid,fsech(i),h%ihist,fakeflds)
            endif
          enddo
        endif
        if (h%nfmag > 0) then
          do i=1,h%nfmag ! magnetic fields
            call wrfdiag(ncid,fsechmag(i),h%ihist,fakeflds)
          enddo
        endif
      else
        write(6,"(/,'>>> WARNING nc_wrhist: unknown h%type=',a)") 
     |    h%type
      endif
!
! Report to stdout:
      if (iprint > 0) call hist_print(h,'WRITE')
      end subroutine nc_wrhist
!-------------------------------------------------------------------
      subroutine wrfdiag(ncid,f,itime,fake)
!
! Write diagnostic field to current open history file:
!
! f%data should be allocated and defined.
! If the model has taken at least one timestep (istep > 0), then 
!   f%data should have been defined by user-called sub addfsech. 
! If istep==0, then f%data was init to spval by set_fsech at beginning 
!   of run.
!
      use fields_module,only: field,fprog,fakeflds
      use init_module,only: istep
      include "params.h" ! for spval and grid dims
      include "grid.h"   ! for grid dims
!
! Args:
      integer,intent(in) :: ncid,itime
      type(field),intent(in) :: f 
      logical,intent(in) :: fake
!
! Local:
      integer :: i,istat,idvar,itype,ndims,iddims(4),natts,idimsizes(4)
      character(len=120) :: char120
      character(len=80) :: char80
      character(len=16) :: rdname
      real :: f3d(nlon,nlat,nlev),f3dmag(nmlon,nmlat,nmlev)
      real :: fmin,fmax
      real :: fakevar(1,1,1)
!
! Fake means 3d fields are dimensioned (1,1,1) for testing:
! (see fakeflds in flds_mod.f)
      if (fake) then
        istat = nf_inq_varid(ncid,f%name,idvar)
        if (istat /= NF_NOERR) then
          write(char120,"('wrfdiag: Error getting id of field var ',a)")
     |      trim(f%name)
          call handle_ncerr(istat,char120)
        endif
        start_4d(1:3) = 1
        start_4d(4) = itime
        count_4d(1) = 1
        count_4d(2) = 1
        count_4d(3) = 1
        count_4d(4) = 1
        fakevar(1,1,1) = 0.
        istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d,fakevar)
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_vara_double',
     +      ' for fake field var ',a,' itime=',i2)") f%name,itime
          call handle_ncerr(istat,char120)
        endif
        return
      endif
!
      if (istep > 0) then
        if (any(f%data /= spval)) then ! was defined

        else
          write(6,"(/,'>>> WARNING wrfdiag: data for ',
     |      'diagnostic field ',a,' is all spval --',/,4x,
     |      'data was apparently not defined by addfsech.')") 
     |      trim(f%name)
        endif
      else  ! istep == 0

      endif
!
! Get field id:
      istat = nf_inq_varid(ncid,f%name,idvar)
      if (istat /= NF_NOERR) then
        write(char80,"('wrfdiag: Error getting id of field var ',a)")
     |    trim(f%name)
        call handle_ncerr(istat,char80)
      endif
!
! Get info about the field: 
      istat = nf_inq_var(ncid,idvar,rdname,itype,ndims,iddims,natts)
      if (ndims /= 4) then
        write(6,"(/,'>>> WARNING wrfdiag: bad ndims=',i3,
     |    ' (every diagnostic should have 4 dimensions)')") ndims
      endif
!
! Get info about dimensions:
      do i=1,ndims
        istat = nf_inq_dim(ncid,iddims(i),rdname,idimsizes(i))
      enddo
!
! Check dimension sizes:
      if (.not.fakeflds) then
        if (.not.f%magnetic) then ! geographic grid
          if (idimsizes(1) /= nlon .or. idimsizes(2) /= nlat .or. 
     |        idimsizes(3) /= nlev) then
            write(6,"(/,'>>> WARNING wrfdiag: bad dimension sizes',
     |        ' for diagnostic ',a)") trim(rdname)
            write(6,"('  dim sizes=',4i4,' but should be ',
     |        'nlon,nlat,nlev=',4i4)") idimsizes,nlon,nlat,nlev
          endif
        else                      ! magnetic grid
          if (idimsizes(1) /= nmlon .or. idimsizes(2) /= nmlat .or. 
     |        idimsizes(3) /= nmlev) then
            write(6,"(/,'>>> WARNING wrfdiag: bad dimension sizes',
     |        ' for magnetic diagnostic field ',a)") trim(rdname)
            write(6,"('  dim sizes=',4i4,' but should be ',
     |        'nmlon,nmlat,nmlev=',4i4)") idimsizes,nmlon,nmlat,nmlev
          endif

        endif
      endif
!
! Define f3d from f%data (netcdf does not like the pointer f%data):
!     local: real :: f3d(nlon,nlat,nlev)
!     local: real :: f3dmag(nmlon,nmlat,nmlev)
!     flds_mod: allocate(fsech(nfsech)%data(zimxp,zjmx,zkmxp),stat=ier)
!               allocate(fsechmag(nfsech)%data(imaxmp,jmaxm,zkmxp+3),
!
      if (.not.f%magnetic) then ! geographic grid
        do i=1,nlon
          f3d(i,:,:) = f%data(i+2,:,:)
        enddo
      else                      ! magnetic grid
        do i=1,nmlon
          f3dmag(i,:,:) = f%data(i,:,:)
        enddo
      endif
!
! Write data to the file:
      start_4d(1:3) = 1
      start_4d(4) = itime
      count_4d(4) = 1
      if (.not.f%magnetic) then ! geographic grid
        count_4d(1) = nlon
        count_4d(2) = nlat
        count_4d(3) = nlev
        istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d,f3d)
        if (istat /= NF_NOERR) then
          write(char80,"(/,'>>> wrfdiag: error return from ',
     |      'nf_put_vara_double for geo field var ',a,' itime=',i2)") 
     |      trim(f%name),itime
          call handle_ncerr(istat,char80)
        endif
      else                      ! magnetic grid
        count_4d(1) = nmlon
        count_4d(2) = nmlat
        count_4d(3) = nmlev
        istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d,f3dmag)
        if (istat /= NF_NOERR) then
          write(char80,"(/,'>>> wrfdiag: error return from ',
     |      'nf_put_vara_double for mag field var ',a,' itime=',i2)") 
     |      trim(f%name),itime
          call handle_ncerr(istat,char80)
        endif
      endif
!
! If model has executed one or more time steps, the long_name and
! units of the diagnostic sech field may have been defined by the
! user-called addfsech routine:
!
      if (istep > 0) then
        if (len_trim(f%long_name) > 0 .or. len_trim(f%units) > 0) then
          istat = nf_redef(ncid) ! put dataset in define mode
          if (istat /= NF_NOERR) call handle_ncerr(istat,
     +      'Error return from nf_redef')
          if (len_trim(f%long_name) > 0) then
            istat = nf_put_att_text(ncid,idvar,"long_name",
     |        len_trim(f%long_name),f%long_name)
            if (istat /= NF_NOERR) then
              write(char80,"('Error defining long_name of diagnostic',
     |          ' variable ',a,': long_name = ',a)") 
     |          trim(f%name),trim(f%long_name)
              call handle_ncerr(istat,trim(char80))
            endif
          endif
          if (len_trim(f%units) > 0) then
            istat = nf_put_att_text(ncid,idvar,"units",
     |        len_trim(f%units),f%units)
            if (istat /= NF_NOERR) then
              write(char80,"('Error defining units of diagnostic',
     |          ' variable ',a,' units = ',a)") 
     |          trim(f%name),trim(f%units)
              call handle_ncerr(istat,trim(char80))
            endif
          endif
          istat = nf_enddef(ncid) ! end define mode
          if (istat /= NF_NOERR) call handle_ncerr(istat,
     |      'wrfdiag: Error return from nf_enddef')
        endif ! long_name or units defined
      endif   ! istep > 0
!
! Write min,max of each field:
!     if (istep /= 0) then
!       call fminmaxspv(f3d,nlon*nlat*nlev,fmin,fmax,spval)
!       write(6,"('wrfdiag: Wrote field ',a,' istep=',i3,
!    |    ' 3d min,max=',2e12.4)") f%name(1:8),istep,fmin,fmax
!     endif
      end subroutine wrfdiag
!-------------------------------------------------------------------
      subroutine wrfprog(ncid,name,itime,fake)
!
! Write prognostic field f from fg-array to current netcdf output 
! history file attached to ncid:
!
      use fields_module,only: findnx
      use init_module,only: istep
      include "params.h"
      include "fgcom.h"
      include "grid.h"
      include "dynphi.h"  ! for DYNPOT
!
! Args:
      integer,intent(in) :: ncid,itime
      character(len=*),intent(in) :: name
      logical,intent(in) :: fake
!
! Local:
      integer :: nx,nxk,i,k,j,idvar,istat
      real :: f3d(nlon,nlat,nlev),fmin,fmax 
      character(len=120) :: char120
      real :: fakevar(1,1,1)
!
! Fake means 3d fields are dimensioned (1,1,1):
! (see fakeflds in flds_mod.f)
      if (fake) then
        istat = nf_inq_varid(ncid,name,idvar)
        if (istat /= NF_NOERR) then
          write(char120,"('wrfprog: Error getting id of field var ',a)")
     |      trim(name)
          call handle_ncerr(istat,char120)
        endif
        start_4d(1:3) = 1
        start_4d(4) = itime
        count_4d(1) = 1
        count_4d(2) = 1
        count_4d(3) = 1
        count_4d(4) = 1
        fakevar(1,1,1) = 0.
        istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d,fakevar)
        if (istat /= NF_NOERR) then
          write(char120,"('Error return from nf_put_vara_double',
     +      ' for fake field var ',a,' itime=',i2)") name,itime
          call handle_ncerr(istat,char120)
        endif
        return
      endif
!
! Get "pointer" index to fg: 
      nx = findnx(name)
      if (nx < 0) write(6,"(/,'>>> wrfprog: bad return from findnx',
     |  ' for prognostic field ',a)") trim(name)
!
! Define local f3d from fg. Note all prognostics in fg should have
! kmaxp1==nlev levels (NPHIH is not written to netcdf histories).
! Also note longitudes 3-74 from fg are used.
! (electric potential is in dynpot)
!
      if (trim(name) /= "POTEN") then
        do j=1,nlat
          nxk = nx
          do k=1,nlev
            nxk = nxk+1
            do i=1,nlon
              f3d(i,j,k) = fg(i+2,nxk,j,ixtimec)
            enddo
          enddo
        enddo
      else
        f3d(1:nlon,1:nlat,:) = dynpot(1:nlon,1:nlat,:)
      endif
!
! Get field id:
      istat = nf_inq_varid(ncid,name,idvar)
      if (istat /= NF_NOERR) then
        write(char120,"('wrfprog: Error getting id of field var ',a)")
     |    trim(name)
        call handle_ncerr(istat,char120)
      endif
!
! Put data onto netcdf file (this is where the majority of output
!   i/o happens). 
! (Should probably make a nc_put_vara_real for SUN vs non-SUN)
!
      start_4d(1:3) = 1
      start_4d(4) = itime
      count_4d(1) = nlon
      count_4d(2) = nlat
      count_4d(3) = nlev
      count_4d(4) = 1
      istat = nf_put_vara_double(ncid,idvar,start_4d,count_4d,f3d)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_put_vara_double',
     +    ' for field var ',a,' itime=',i2)") name,itime
        call handle_ncerr(istat,char120)
      endif
!
! Write min,max of each field if not echoing source history:
!     if (istep /= 0) then
!       call fminmaxspv(f3d,nlon*nlat*nlev,fmin,fmax,spval)
!       write(6,"('wrfprog: Wrote field ',a,' istep=',i3,
!    |    ' 3d min,max=',2e12.4)") name(1:8),istep,fmin,fmax
!     endif
      end subroutine wrfprog
!-------------------------------------------------------------------
      subroutine nc_fileinfo(ncid,fileinfo)
!
! A netcdf file on ncid has been written and is ready to be
! closed and disposed. This routine builds a fileinfo string
! giving beginning and ending history dates, model version string,
! and "primary" or "secondary". This is added as a global
! file attribute, and is included in the mss comment field
! if the file is disposed to the mss, and is therefore accessible 
! via msls -x. Example fileinfo string:
!
! 96080  80, 0, 0 to 96084  84, 0, 0 by  720 tgcm22 primary
!
! Args:
      integer,intent(in) :: ncid
      character(len=80),intent(out) :: fileinfo
!
! Local:
      integer :: istat,idvmtime,idtime,ntimes,idvyear,idvday,
     |  mtime_beg(3),mtime_end(3),istart2(2),icount2(2),
     |  iyear_beg,iyear_end,iday_beg,iday_end,idelhist,
     |  ibegiyd,iendiyd
      character(len=16) :: mversion,htype
      character(len=80) :: contents_desc
!
! Init:
      fileinfo = ' '
!
! Get number of histories (times) on the file:
      istat = nf_inq_unlimdim(ncid,idtime) ! id of unlimited record var
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'nc_fileinfo: Error getting unlimited dimension id')
      istat = nf_inq_dimlen(ncid,idtime,ntimes) ! length of time var
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'nc_fileinfo: Error getting length of time record variable')
!
! Get id for mtime variable:
      istat = nf_inq_varid(ncid,"mtime",idvmtime) 
      if (istat /= NF_NOERR)
     |  call handle_ncerr(istat,
     |    'nc_fileinfo: Error getting id of variable mtime')
!
! Get beginning and ending mtimes:
      istart2(1) = 1
      istart2(2) = 1
      icount2(1) = 3
      icount2(2) = 1
      istat = nf_get_vara_int(ncid,idvmtime,istart2,icount2,
     |  mtime_beg)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'nc_fileinfo: Error getting beginning mtime')
      istart2(2) = ntimes
      istat = nf_get_vara_int(ncid,idvmtime,istart2,icount2,
     |  mtime_end)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'nc_fileinfo: Error getting ending mtime')
!
! Get beginning and ending year:
      istat = nf_inq_varid(ncid,"year",idvyear) 
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'nc_fileinfo: Error getting id of variable year')
      istat = nf_get_var1_int(ncid,idvyear,1,iyear_beg)
      istat = nf_get_var1_int(ncid,idvyear,ntimes,iyear_end)
!
! Get beginning and ending day:
      istat = nf_inq_varid(ncid,"day",idvday) 
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |    'nc_fileinfo: Error getting id of variable day')
      istat = nf_get_var1_int(ncid,idvday,1,iday_beg)
      istat = nf_get_var1_int(ncid,idvday,ntimes,iday_end)
!
! Get delta time between histories (minutes), model version,
! and history type:
      istat = nf_get_att_int (ncid,NF_GLOBAL,"delhist_mins",idelhist)
      mversion = ' '
      istat = nf_get_att_text(ncid,NF_GLOBAL,"model_version",mversion)
      htype = ' '
      istat = nf_get_att_text(ncid,NF_GLOBAL,"history_type",htype)
!
! Construct file info string:
      ibegiyd = (iyear_beg-iyear_beg/100*100)*1000+iday_beg
      iendiyd = (iyear_end-iyear_end/100*100)*1000+iday_end
      write(fileinfo,"(i5.5,' ',i3,' ',i2,' ',i2,' to ',
     |                 i5.5,' ',i3,' ',i2,' ',i2,' by ',i4,
     |  ' ',a,' ',a)") ibegiyd, mtime_beg,iendiyd, mtime_end, 
     |  idelhist,trim(mversion),trim(htype)
!
! Add file info string and description of file info global attributes:
      istat = nf_redef(ncid)
      istat = nf_put_att_text(ncid,NF_GLOBAL,"contents",
     |  len_trim(fileinfo),trim(fileinfo))
      write(contents_desc,"('yyddd day hour min',
     |  ' to yyddd day hour min by delta_mins')")
      istat = nf_put_att_text(ncid,NF_GLOBAL,"contents_desc",
     |  len_trim(contents_desc),trim(contents_desc))
      istat = nf_enddef(ncid)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     +  'nc_fileinfo: Error return from nf_enddef')
      end subroutine nc_fileinfo
!-------------------------------------------------------------------
      subroutine handle_ncerr(istat,msg)
!
! Handle a netcdf lib error:
!
      integer,intent(in) :: istat
      character(len=*),intent(in) :: msg
!
      write(6,"(/72('-'))")
      write(6,"('>>> Error from netcdf library:')")
      write(6,"(a)") trim(msg)
      write(6,"('istat=',i5)") istat
      write(6,"(a)") nf_strerror(istat)
      write(6,"(72('-')/)")
      return
      end subroutine handle_ncerr
!-------------------------------------------------------------------
!
! Remaining functions are netcdf wrappers for 32-bit vs 64-bit platforms.
!
!-------------------------------------------------------------------
      integer function nc_put_att_real(ncid,varid,name,len,vals)
!
! Args:
      integer,intent(in) :: ncid,varid,len
      character(len=*),intent(in) :: name
      real,intent(in) :: vals(len)
!
! Local:
      integer :: istat
!
#if defined(SUN)
      istat = nf_put_att_real(ncid,varid,name,NF_FLOAT,
     |  len,vals)
#else
      istat = nf_put_att_double(ncid,varid,name,NF_DOUBLE,
     |  len,vals)
#endif
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'nc_put_att_real: error from nf_put_att_real/double')
      nc_put_att_real = istat
      end function nc_put_att_real
!-------------------------------------------------------------------
      integer function nc_put_var_real(ncid,idvar,vals)
!
! Args:
      integer,intent(in) :: ncid,idvar
      real,intent(in) :: vals(*)
!
! Local:
      integer :: istat
!
#if defined(SUN)
      istat = nf_put_var_real(ncid,idvar,vals)
#else
      istat = nf_put_var_double(ncid,idvar,vals)
#endif
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'nc_put_var_real: error from nf_put_var_real/double')
      nc_put_var_real = istat
      end function nc_put_var_real
!-------------------------------------------------------------------
      integer function nc_get_var_real(ncid,idvar,vals)
!
! Args:
      integer,intent(in) :: ncid,idvar
      real,intent(out) :: vals(*)
!
! Local:
      integer :: istat
!
#if defined(SUN)
      istat = nf_get_var_real(ncid,idvar,vals)
#else
      istat = nf_get_var_double(ncid,idvar,vals)
#endif
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'nc_get_var_real: error from nf_get_var_real/double')
      nc_get_var_real = istat
      end function nc_get_var_real
!-------------------------------------------------------------------
      integer function nc_get_var1_real(ncid,idvar,index,val)
!
! Args:
      integer,intent(in) :: ncid,idvar,index
      real,intent(out) :: val
!
! Local:
      integer :: istat
!
#if defined(SUN)
      istat = nf_get_var1_real(ncid,idvar,index,val)
#else
      istat = nf_get_var1_double(ncid,idvar,index,val)
#endif
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'nc_get_var1_real: error from nf_get_var1_real/double')
      nc_get_var1_real = istat
      end function nc_get_var1_real
!-------------------------------------------------------------------
      integer function nc_put_var1_real(ncid,idvar,index,val,name)
!
! Args:
      integer,intent(in) :: ncid,idvar,index
      real,intent(in) :: val
      character(len=*),intent(in) :: name
!
! Local:
      integer :: istat
      character(len=80) :: msg
!
#if defined(SUN)
      istat = nf_put_var1_real(ncid,idvar,index,val)
#else
      istat = nf_put_var1_double(ncid,idvar,index,val)
#endif
      if (istat /= NF_NOERR) then
        write(msg,"('nc_put_var1_real: error from nf_put_var1_real',
     |    '/double for variable ',a)") trim(name)
        call handle_ncerr(istat,msg)
      endif
      nc_put_var1_real = istat
      end function nc_put_var1_real
      end module netcdf_module
