!
      module readhist_module
!
! Read a TIEGCM/TIMEGCM history, defining a history structure.
!
      implicit none
#include "netcdf.inc"
#include "hist.h"
!
      real,parameter    :: real_unread = 99999.
      integer,parameter :: int_unread = 99999
      character(len=mxlen) :: input_file
      contains
!-----------------------------------------------------------------------
      subroutine readhist(infile,mtime,h,iprint,error)
      implicit none
!
! Read history model time mtime from netcdf history file, and return 
!   history structure h.
!
! Args:
! History file to read from:
      character(len=*),intent(in) :: infile
!
! h: history structure to return (read from infile into h):
      type(history),intent(out) :: h 
!
! mtime(3): model time (day,hr,min) of history to be read:
      integer,intent(inout) :: mtime(3)
      integer,intent(in) :: iprint    ! 0/1 flag for print to stdout
      integer,intent(out) :: error    ! error return
!
! Local:
      integer,parameter :: mxdims=20 ! assume <= mxdims dimensions
      integer :: i,id,ncid,istat,ier,itime
      integer :: ndims,nvars,ngatts,idunlim,itype,natts,len,
     |  iddims(mxdims),mtime_read(3)
      character(len=80) :: varname
      real :: var22(2,2),var2(2),var10(10)
      integer :: istart2(2),icount2(2)
      character(len=8) :: type
      logical :: exists
!
! Check existence of disk file:
      error = 0
      inquire(file=infile,exist=exists)
      if (.not.exists) then
        write(6,"('>>> Could not find file ',a)") trim(infile)
        error = 1
        return
      endif
!
! Open history file (assume netcdf):
      istat = nf_open(infile,NF_NOWRITE,ncid)
      if (istat /= NF_NOERR) then
        write(6,"(/,'>>> Error opening file ',a)") trim(infile)
        write(6,"('    May not be a valid netcdf file.',/)")
        error = 1
        return
      else
        if (iprint > 0)
     |    write(6,"('readhist: Opened netcdf file ',a)") trim(infile)
      endif
      h%history_file = trim(infile)
      input_file = trim(infile) ! module data, used by writehist
!
! Find desired mtime history on the file:
      call findhist(ncid,mtime,itime,ier)
      if (ier /= 0) then
        error = 1
        return
      endif
!
! Print this even if iprint==0:
      write(6,"('Reading history ',3i4,' (itime=',i3,') from file ',a)")
     |  mtime,itime,trim(infile)
!
! Initialize history structure:
      call hist_init(h)
!
! Get number of histories and model times on the file:
      istat = nf_inq_unlimdim(ncid,id)         ! id of unlimited record var
      istat = nf_inq_dimlen(ncid,id,h%nhist)   ! number of histories
      istat = nf_inq_varid(ncid,"mtime",id)    ! model time id
      allocate(h%mtimes(3,h%nhist),stat=ier)
      istart2(1) = 1
      icount2(1) = 3
      icount2(2) = 1
      do i=1,h%nhist
        istart2(2) = i
        istat = nf_get_vara_int(ncid,id,istart2,icount2,
     |    mtime_read)
        h%mtimes(:,i) = mtime_read
      enddo 
!     write(6,"('readhist: h%nhist=',i3,' h%mtimes=',/,(3i4))")
!    |  h%nhist,h%mtimes
!
! Define vars already found above:
      h%ihist = itime  ! nth history on input file

      h%modeltime(1:3) = mtime
      h%modeltime(4) = 0
!
! Read global attributes:
!
      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "Conventions",h%conventions)
      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "label",h%label)
      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "create_date",h%create_date)
      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,
     |  "output_file",h%output_file)
      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "history_type",h%history_type)
      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "run_type",h%run_type)
      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "source_file",h%source_file)
      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "initial_file",h%initial_file)

      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "lev_to_hPa_method1",h%lev_to_hPa_method1)
      if (istat /= 0) write(h%lev_to_hPa_method1,
     |  "('p0*exp(-lev(k))')")

      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "lev_to_hPa_method2",h%lev_to_hPa_method2)
      if (istat /= 0) write(h%lev_to_hPa_method2,
     |  "('p0_model*1.e-3*exp(-lev(k))')")

      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "potential_model",h%potential_model)
      istat = nf_get_att_int(ncid,NF_GLOBAL,
     |  "source_mtime",h%source_mtime)
      istat = nf_get_att_int(ncid,NF_GLOBAL,
     |  "initial_mtime",h%initial_mtime)
      istat = nf_get_att_int(ncid,NF_GLOBAL,
     |  "tuv_lbc_intop",h%tuv_lbc_intop)
      istat = nf_get_att_int(ncid,NF_GLOBAL,
     |  "delhist_mins",h%delhist_mins)

      istat = nf_get_att_double(ncid,NF_GLOBAL,
     |  "missing_value",h%missing_value)

      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "contents",h%contents)
      istat = nf_get_att_text(ncid,NF_GLOBAL,
     |  "contents_desc",h%contents_desc)
!
! Get number of dims, vars, atts, and id of unlimited dimension:
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
!
! Get dimensions and allocate coordinate vars:
      do i=1,ndims
        istat = nf_inq_dim(ncid,i,varname,len)
        select case(trim(varname))
!
! Geographic dims:
          case('lon')
            h%nlon = len
            allocate(h%glon(h%nlon),stat=ier)
          case('lat')
            h%nlat = len
            allocate(h%glat(h%nlat),stat=ier)
          case('lev')
            h%nlev = len
            allocate(h%lev(h%nlev),stat=ier)
          case('ilev')
            h%nilev = len
            allocate(h%ilev(h%nilev),stat=ier)
!
! Magnetic dims:
          case('mlon')
            h%mlon = len
            allocate(h%gmlon(h%mlon),stat=ier)
          case('mlat')
            h%mlat = len
            allocate(h%gmlat(h%mlat),stat=ier)
          case('mlev')
            h%mlev = len
            allocate(h%gmlev(h%mlev),stat=ier)
          case('imlev')
            h%imlev = len
            allocate(h%gimlev(h%imlev),stat=ier)
!
! Char lens:
          case('datelen')
            h%datelen = len
          case('filelen')
            h%filelen = len
        end select
      enddo
!
! Get variables:
      h%nf3d_geo = 0
      h%nf3d_mag = 0
      h%nf2d_geo = 0
      h%nf2d_mag = 0
      do i=1,nvars
        istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts)
!       write(6,"('i=',i4,' nvars=',i4,' varname=',a,' ndims=',i3)") 
!    |    i,nvars,trim(varname),ndims
        select case(trim(varname))
          case('lon')
            istat = nf_get_var_double(ncid,i,h%glon)
          case('lat')
            istat = nf_get_var_double(ncid,i,h%glat)
          case('lev')
            istat = nf_get_var_double(ncid,i,h%lev)
          case('ilev')
            istat = nf_get_var_double(ncid,i,h%ilev)
          case('mlon')
            istat = nf_get_var_double(ncid,i,h%gmlon)
          case('mlat')
            istat = nf_get_var_double(ncid,i,h%gmlat)
          case('mlev')
            istat = nf_get_var_double(ncid,i,h%gmlev)
          case('imlev')
            istat = nf_get_var_double(ncid,i,h%gimlev)
          case('time')
            istat = nf_get_var1_double(ncid,i,itime,h%time)
          case('mtime') ! already got it from hist search above
          case('calendar_advance')
            istat = nf_get_var1_int(ncid,i,itime,h%calendar_advance)
          case('year')
            istat = nf_get_var1_int(ncid,i,itime,h%year)
          case('day')
            istat = nf_get_var1_int(ncid,i,itime,h%day)
          case('iter')
            istat = nf_get_var1_int(ncid,i,itime,h%iter)
          case('ut')
            istat = nf_get_var1_double(ncid,i,itime,h%ut)
          case('mag')
            istat = nf_get_var_double(ncid,i,var22)
            h%mag(:,:) = var22(:,:)
          case('dtide')
            istart2(1) = 1
            istart2(2) = h%ihist
            icount2(1) = 2
            icount2(2) = 1
            istat = nf_get_vara_double(ncid,i,istart2,icount2,var2)
            h%dtide(:) = var2(:)
          case('sdtide')
            istart2(1) = 1
            istart2(2) = h%ihist
            icount2(1) = 10
            icount2(2) = 1
            istat = nf_get_vara_double(ncid,i,istart2,icount2,var10)
            h%sdtide(:) = var10(:)

          case('f107d')
            istat = nf_get_var1_double(ncid,i,itime,h%f107d)
          case('f107a')
            istat = nf_get_var1_double(ncid,i,itime,h%f107a)
          case('hpower')
            istat = nf_get_var1_double(ncid,i,itime,h%hpower)
          case('ctpoten')
            istat = nf_get_var1_double(ncid,i,itime,h%ctpoten)

          case('bximf')
            istat = nf_get_var1_double(ncid,i,itime,h%bximf)
          case('byimf')
            istat = nf_get_var1_double(ncid,i,itime,h%byimf)
          case('bzimf')
            istat = nf_get_var1_double(ncid,i,itime,h%bzimf)
          case('swvel')
            istat = nf_get_var1_double(ncid,i,itime,h%swvel)
          case('swden')
            istat = nf_get_var1_double(ncid,i,itime,h%swden)

          case('colfac')
            istat = nf_get_var1_double(ncid,i,itime,h%colfac)

          case('alfa30')
            istat = nf_get_var1_double(ncid,i,itime,h%alfa30)
          case('e30')
            istat = nf_get_var1_double(ncid,i,itime,h%e30)
          case('ed2')
            istat = nf_get_var1_double(ncid,i,itime,h%ed2)
          case('alfad2')
            istat = nf_get_var1_double(ncid,i,itime,h%alfad2)

! al,e1,e2,h1,h2,alfac,ec,alfad,ed, ! auroral parameters ("new" hist)
          case('al')
            istat = nf_get_var1_double(ncid,i,itime,h%al)
          case('e1')
            istat = nf_get_var1_double(ncid,i,itime,h%e1)
          case('e2')
            istat = nf_get_var1_double(ncid,i,itime,h%e2)
          case('h1')
            istat = nf_get_var1_double(ncid,i,itime,h%h1)
          case('h2')
            istat = nf_get_var1_double(ncid,i,itime,h%h2)
          case('alfac')
            istat = nf_get_var1_double(ncid,i,itime,h%alfac)
          case('ec')
            istat = nf_get_var1_double(ncid,i,itime,h%ec)
          case('alfad')
            istat = nf_get_var1_double(ncid,i,itime,h%alfad)
          case('ed')
            istat = nf_get_var1_double(ncid,i,itime,h%ed)

          case('p0')
            istat = nf_get_var_double(ncid,i,h%p0)
          case('p0_model')
            istat = nf_get_var_double(ncid,i,h%p0_model)
          case('grav')
            istat = nf_get_var_double(ncid,i,h%grav)
          case('LBC')
            istat = nf_get_var_double(ncid,i,h%lbc)
          case('timestep')
            istat = nf_get_var_int(ncid,i,h%step)
          case('nflds')
!           no longer used (old histories)
          case('ncep')
            istat = nf_get_var1_int(ncid,i,itime,h%ncep)
          case('gpi')
            istat = nf_get_var1_int(ncid,i,itime,h%gpi)
! integer :: gswmdi,gswmsdi,gswmnmdi,gswmnmsdi
          case('gswmdi')
            istat = nf_get_var1_int(ncid,i,itime,h%gswmdi)
          case('gswmsdi')
            istat = nf_get_var1_int(ncid,i,itime,h%gswmsdi)
          case('gswmnmdi')
            istat = nf_get_var1_int(ncid,i,itime,h%gswmnmdi)
          case('gswmnmsdi')
            istat = nf_get_var1_int(ncid,i,itime,h%gswmnmsdi)
          case('write_date')
            if (h%datelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%datelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%write_date)
            endif
          case('gpi_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%gpi_ncfile)
            endif
          case('ncep_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%ncep_ncfile)
            endif
          case('imf_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%imf_ncfile)
            endif
          case('ecmwf_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%ecmwf_ncfile)
            endif
          case('see_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%see_ncfile)
            endif
          case('gswm_mi_di_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%gswm_mi_di_ncfile)
            endif
          case('gswm_mi_sdi_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%gswm_mi_sdi_ncfile)
            endif
          case('gswm_nm_di_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%gswm_nm_di_ncfile)
            endif
          case('gswm_nm_sdi_ncfile')
            if (h%filelen /= int_unread) then
              istart2(1) = 1         ; istart2(2) = itime
              icount2(1) = h%filelen ; icount2(2) = 1
              istat = nf_get_vara_text(ncid,i,istart2,icount2,
     |          h%gswm_nm_sdi_ncfile)
            endif
          case default
!
! Determine if field is 3d or 2d, geo or mag:
!   do i=1,nvars
!     istat = nf_inq_var(ncid,i,varname,itype,ndims,iddims,natts)
!
            call ftype(h,ncid,varname,i,ndims,iddims,type) ! type is returned
            select case(trim(type))
              case('3dgeo')
                h%nf3d_geo = h%nf3d_geo+1
                h%f3d_geo(h%nf3d_geo)%short_name = varname
                call read_f3d(ncid,varname,h,'geo',itime)
              case('3dmag')
                h%nf3d_mag = h%nf3d_mag+1
                h%f3d_mag(h%nf3d_mag)%short_name = varname
                call read_f3d(ncid,varname,h,'mag',itime)
              case('2dgeo')
                h%nf2d_geo = h%nf2d_geo+1
                h%f2d_geo(h%nf2d_geo)%short_name = varname
                call read_f2d(ncid,varname,h,'geo',itime)
              case('2dmag')
                h%nf2d_mag = h%nf2d_mag+1
                h%f2d_mag(h%nf2d_mag)%short_name = varname
                call read_f2d(ncid,varname,h,'mag',itime)
              case('other')
                if (iprint > 0) write(6,"('Unknown variable not read:',
     |            a)") trim(varname)
              case default
                write(6,"('>>> Warning: unknown type=',a)") trim(type)
            end select ! case type
        end select ! case(trim(varname))
      enddo ! i=1,nvars
      if (h%nilev == int_unread) then                    ! is an "old" history
        h%old = .true.
        h%new = .false.
      else
        h%old = .false.
        h%new = .true.
      endif
!
! If old histories do not have certain variables, use missing value:
! (if there are others, they will be left real_unread or int_unread,
!  or empty strings)
!
      if (h%bximf == real_unread) h%bximf = h%missing_value
      if (h%byimf == real_unread) h%byimf = h%missing_value
      if (h%bzimf == real_unread) h%bzimf = h%missing_value
      if (h%swvel == real_unread) h%swvel = h%missing_value
      if (h%swden == real_unread) h%swden = h%missing_value
      if (h%al    == real_unread) h%al    = h%missing_value
      if (h%e1    == real_unread) h%e1    = h%missing_value
      if (h%e2    == real_unread) h%e2    = h%missing_value
      if (h%h1    == real_unread) h%h1    = h%missing_value
      if (h%h2    == real_unread) h%h2    = h%missing_value
      if (h%alfac == real_unread) h%alfac = h%missing_value
      if (h%ec    == real_unread) h%ec    = h%missing_value
      if (h%alfad == real_unread) h%alfad = h%missing_value
      if (h%ed    == real_unread) h%ed    = h%missing_value
      if (h%grav  == real_unread) h%grav  = h%missing_value
      if (h%p0_model == real_unread) h%p0_model = .0005
      if (h%gpi > 0 .and. len_trim(h%gpi_ncfile) == 0)
     |  write(h%gpi_ncfile,"('[yes]')")
      if (h%ncep > 0 .and. len_trim(h%ncep_ncfile) == 0)
     |  write(h%ncep_ncfile,"('[yes]')")
! integer :: gswmdi,gswmsdi,gswmnmdi,gswmnmsdi
      if (h%gswmdi > 0 .and. len_trim(h%gswm_mi_di_ncfile)==0)
     |  write(h%gswm_mi_di_ncfile,"('[yes]')")
      if (h%gswmsdi > 0 .and. len_trim(h%gswm_mi_sdi_ncfile)==0)
     |  write(h%gswm_mi_sdi_ncfile,"('[yes]')")
      if (h%gswmnmdi > 0 .and. len_trim(h%gswm_nm_di_ncfile)==0)
     |  write(h%gswm_nm_di_ncfile,"('[yes]')")
      if (h%gswmnmsdi > 0 .and. len_trim(h%gswm_nm_sdi_ncfile)==0)
     |  write(h%gswm_nm_sdi_ncfile,"('[yes]')")
!
      if (iprint > 0) call hist_print(h,infile)
      istat = nf_close(ncid)
      write(6,"('End subroutine readhist')")
      end subroutine readhist
!-----------------------------------------------------------------------
      subroutine findhist(ncid,mtime,itime,ier)
      implicit none
!
! Position file connected to ncid at history mtime.
! If history not found, return ier > 0.
!
! Args:
      integer,intent(inout) :: mtime(3)
      integer,intent(in) :: ncid
      integer,intent(out) :: itime,ier
!
! Local:
      integer :: i,istat,ntimes,mtime_read(3),idv_mtime,id_time
      integer :: istart2(2),icount2(2)
!
! Init:
      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')
!     write(6,"('findhist: There are ',i3,' histories on this file')")
!    |  ntimes
!
! 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,"('findhist: seeking ',3i4,' found ',3i4,'  i=',i3)")
!    |    mtime,mtime_read,i
        if (mtime(1) > 366) then
          itime = i
          mtime(:) = mtime_read(:)
          return
        endif
        if (all(mtime_read==mtime)) then
          itime = i
          return
        endif
      enddo
      if (itime==0) then
        write(6,"(/,'>>> findhist: mtime ',3i4,' NOT found.')") mtime
        ier = 1
        return
      else
!       write(6,"('findhist: Found history ',3i4,' at n=',i4)") 
!    |    mtime,itime
      endif
      end subroutine findhist
!-----------------------------------------------------------------------
      subroutine hist_init(h)
      implicit none
!
! Initialize a history structure:
!
      type(history),intent(out) :: h
!
! Global text attributes:
      h%Conventions = ' ' 
      h%label = ' ' 
      h%create_date = ' ' 
      h%logname = ' ' 
      h%host = ' ' 
      h%system = ' ' 
      h%model_name = ' ' 
      h%model_version = ' ' 
      h%output_file = ' ' 
      h%history_type = ' ' 
      h%run_type = ' ' 
      h%source_file = ' ' 
      h%initial_file = ' ' 
      h%lev_to_hPa_method1 = ' '
      h%lev_to_hPa_method2 = ' '
      h%potential_model = ' '
      h%contents = ' '
      h%contents_desc = ' '
!
! Integer global attributes:
      h%tuv_lbc_intop = int_unread
      h%delhist_mins = int_unread
      h%initial_mtime = int_unread
      h%source_mtime = int_unread
!
! Double global attribute:
      h%missing_value = real_unread
! 
      h%write_date = ' '
      h%ncep_ncfile = ' '
      h%gpi_ncfile = ' '
      h%see_ncfile = ' '
      h%imf_ncfile = ' '
      h%ecmwf_ncfile = ' '
      h%gswm_mi_di_ncfile = ' '
      h%gswm_mi_sdi_ncfile = ' '
      h%gswm_nm_di_ncfile = ' '
      h%gswm_nm_sdi_ncfile = ' '

      h%nlat = int_unread
      h%nlon = int_unread
      h%nlev = int_unread
      h%nilev = int_unread

      h%mlat = int_unread
      h%mlon = int_unread
      h%mlev = int_unread
      h%imlev = int_unread

      h%time = int_unread
      h%iter = int_unread
      h%year = int_unread
      h%day = int_unread
      h%calendar_advance = int_unread
      h%ut = real_unread
      h%step = int_unread
      h%ncep = int_unread
      h%gpi = int_unread
!
! Old histories: integer :: gswmdi = gswmsdi,gswmnmdi,gswmnmsdi
!
      h%gswmdi  = int_unread
      h%gswmsdi = int_unread
      h%gswmnmdi  = int_unread
      h%gswmnmsdi = int_unread

! These are init to zero in readhist before var loop
!     h%nf3d_geo = int_unread
!     h%nf3d_mag = int_unread
!     h%nf2d_geo = int_unread
!     h%nf2d_mag = int_unread

      h%dynamo = .false.
      h%coupled_ccm = .false.
      h%hpower = real_unread
      h%ctpoten = real_unread
      h%f107d = real_unread
      h%f107a = real_unread
      h%bximf = real_unread
      h%byimf = real_unread
      h%bzimf = real_unread
      h%mag(:,1) = real_unread
      h%mag(:,2) = real_unread
      h%dtide = real_unread
      h%sdtide = real_unread
      h%colfac = real_unread
      h%p0     = real_unread
      h%p0_model = real_unread
      h%alfa30 = real_unread
      h%e30    = real_unread
      h%alfad2 = real_unread
      h%ed2    = real_unread
      h%missing_value = real_unread
      h%swvel  = real_unread
      h%swden  = real_unread
      h%grav   = real_unread
      h%lbc    = real_unread

! al,e1,e2,h1,h2,alfac,ec,alfad,ed, ! auroral parameters ("new" hist)
      h%al = real_unread
      h%e1 = real_unread
      h%e2 = real_unread
      h%h1 = real_unread
      h%h2 = real_unread
      h%alfac = real_unread
      h%ec = real_unread
      h%alfad = real_unread
      h%ed = real_unread

      h%datelen = int_unread
      h%filelen = int_unread
      h%nhist = int_unread

!     if (associated(h%data3d_geo)) deallocate(h%data3d_geo)
!     if (associated(h%data3d_mag)) deallocate(h%data3d_mag)
!     if (associated(h%data2d_geo)) deallocate(h%data2d_geo)
!     if (associated(h%data2d_mag)) deallocate(h%data2d_mag)
!     if (associated(h%mtimes)) deallocate(h%mtimes)

      end subroutine hist_init
!-----------------------------------------------------------------------
      subroutine hist_print(h,file)
      implicit none
!
! Print info about a history structure:
!
! Args:
      type(history),intent(in) :: h
      character(len=*),intent(in) :: file
!
! Locals:
      integer :: i
      real :: fmin,fmax
!
      write(6,"(/,72('-'))")
      write(6,"('Reading from file ',a,/,'modeltime=',4i4,/)") 
     |  trim(file),h%modeltime
      write(6,"(2x,'model_name    = ',a)") trim(h%model_name)
      write(6,"(2x,'model_version = ',a)") trim(h%model_version)
      if (len_trim(h%write_date) > 0)
     |  write(6,"(2x,'write_date = ',a)") trim(h%write_date)
      if (len_trim(h%create_date) > 0)
     |  write(6,"(2x,'create_date = ',a)") trim(h%create_date)
      write(6,"(2x,'logname       = ',a)") trim(h%logname)
      write(6,"(2x,'host          = ',a)") trim(h%host)
      write(6,"(2x,'system        = ',a)") trim(h%system)
      write(6,"(2x,'run_type      = ',a)") trim(h%run_type)
      write(6,"(2x,'ihist         = ',i3,
     |  ' (nth history on history file)')") h%ihist
      if (h%delhmins /= int_unread)
     |  write(6,"(2x,'delhmins= ',i4,
     |    ' (delta minutes between histories)')") h%delhmins
      write(6,"(2x,'calendar year,day = ',i4,',',i3)") h%year,h%day
      if (h%calendar_advance /= int_unread) then
        if (h%calendar_advance <= 0) then
          write(6,"(2x,'calendar_advance=',i4,' (model is NOT being ',
     |      'advanced in calendar time)')") h%calendar_advance
        else
          write(6,"(2x,'calendar_advance=',i4,' (model IS being ',
     |      'advanced in calendar time)')") h%calendar_advance
        endif
      endif
      write(6,"(2x,'modeltime = ',i3,',',i2,',',i2,',',i2,
     |  ' (model time day,hour,minute,seconds)')") h%modeltime
      if (h%time <= 1440.) then
        write(6,"(2x,'time    = ',f10.2,' (minutes in current day)')") 
     |    h%time
      else
        write(6,"(2x,'time    = ',e12.4)") h%time
      endif
      write(6,"(2x,'ut      = ',f5.2,' (ut hours)')") h%ut
      write(6,"(2x,'step    = ',i4,' (time step in seconds)')")
     |  h%step
      write(6,"(2x,'iter    = ',i8,' (number of steps from 0,0,0)')")
     |  h%iter
      write(6,"(2x,'nlat    = ',i4,' (number of latitudes)')") h%nlat
      write(6,"(2x,'nlon    = ',i4,' (number of longitudes)')") h%nlon
      write(6,"(2x,'nlev    = ',i4,' (number of levels)')") h%nlev
      if (h%nilev /= int_unread)
     |  write(6,"(2x,'nilev   = ',i4,' (number of interface levels)')") 
     |    h%nilev
      write(6,"(2x,'mlat    = ',i4,' (number of mag latitudes)')") 
     |  h%mlat
      write(6,"(2x,'mlon    = ',i4,' (number of mag longitudes)')") 
     |  h%mlon
      write(6,"(2x,'mlev    = ',i4,' (number of mag levels)')") h%mlev
      if (h%imlev /= int_unread)
     |  write(6,"(2x,'imlev   = ',i4,' (number of mag interface levels)'
     |    )") h%imlev
      if (any(h%mag < 1.e36)) 
     |  write(6,"(2x,'mag     = ',4e10.2,' (magnetic pole coords)')") 
     |    h%mag
      write(6,"(2x,'dtide   = ',e8.1,' ',f5.1,
     |  ' (amp/phase of diurnal tide)')") h%dtide
      write(6,"(2x,'sdtide  = ',5e8.1,' ',5f5.1,/,
     |  4x,'(amp/phase of semi-diurnal tide)')") h%sdtide
      write(6,"(2x,'f107d   = ',f6.2,' (daily solar flux)')") h%f107d
      write(6,"(2x,'f107a   = ',f6.2,' (average solar flux)')") h%f107a
      write(6,"(2x,'hpower  = ',f5.2,' (Gw)')") h%hpower
      write(6,"(2x,'ctpoten = ',f5.2,' (Volts)')") h%ctpoten
      if (h%bximf /= real_unread) 
     |  write(6,"(2x,'bximf   = ',e12.2)") h%bximf
      if (h%byimf /= real_unread) 
     |  write(6,"(2x,'byimf   = ',e12.2)") h%byimf
      if (h%bzimf /= real_unread) 
     |  write(6,"(2x,'bzimf   = ',e12.2)") h%bzimf
      if (h%swvel /= real_unread) 
     |  write(6,"(2x,'swvel   = ',e12.2)") h%swvel
      if (h%swden /= real_unread) 
     |  write(6,"(2x,'swden   = ',e12.2)") h%swden

! al,e1,e2,h1,h2,alfac,ec,alfad,ed, ! auroral parameters ("new" hist)
      if (h%al /= real_unread) write(6,"(2x,'al  = ',e8.2)") h%al
      if (h%e1 /= real_unread) write(6,"(2x,'e1  = ',e8.2)") h%e1
      if (h%e2 /= real_unread) write(6,"(2x,'e2  = ',e8.2)") h%e2
      if (h%h1 /= real_unread) write(6,"(2x,'h1  = ',e8.2)") h%h1
      if (h%h2 /= real_unread) write(6,"(2x,'h2  = ',e8.2)") h%h2
      if (h%alfac /= real_unread) write(6,"(2x,'alfac  = ',e8.2)") 
     |  h%alfac
      if (h%alfad /= real_unread) write(6,"(2x,'alfad  = ',e8.2)") 
     |  h%alfad
      if (h%ec /= real_unread) write(6,"(2x,'ec  = ',e8.2)") h%ec
      if (h%ed /= real_unread) write(6,"(2x,'ed  = ',e8.2)") h%ed

! Auroral parmeters ("old" hist)
      if (h%alfa30 /= real_unread) 
     |  write(6,"(2x,'alfa30  = ',e8.2,' (KeV)')") h%alfa30
      if (h%e30 /= real_unread) 
     |  write(6,"(2x,'e30     = ',e8.2,' (ergs/cm2/s)')") h%e30
      if (h%alfad2 /= real_unread) 
     |  write(6,"(2x,'alfad2  = ',e8.2,' (KeV)')") h%alfad2
      if (h%ed2 /= real_unread) 
     |  write(6,"(2x,'ed2     = ',e8.2,' (ergs/cm2/s)')") h%ed2

      write(6,"(2x,'colfac  = ',f5.2)") h%colfac
      write(6,"(2x,'p0      = ',e8.2)") h%p0
      if (h%p0_model /= real_unread)
     |  write(6,"(2x,'p0_model = ',e12.4)") h%p0_model
      if (h%grav /= real_unread)
     |  write(6,"(2x,'grav = ',e12.4)") h%grav
      if (h%lbc /= real_unread)
     |  write(6,"(2x,'LBC = ',e12.4)") h%lbc

      if (h%ncep /= int_unread) 
     |  write(6,"(2x,'ncep    = ',i2,' (NMC/NCEP flag)')") h%ncep
      if (h%gpi /= int_unread) 
     |  write(6,"(2x,'gpi     = ',i2,' (GPI flag)')") h%gpi

! integer :: gswmdi,gswmsdi,gswmnmdi,gswmnmsdi
      if (h%gswmdi /= int_unread) 
     |  write(6,"(2x,'gswmdi        = ',i2,' (gswm flag)')") h%gswmdi
      if (h%gswmsdi /= int_unread) 
     |  write(6,"(2x,'gswmsdi       = ',i2,' (gswm flag)')") h%gswmsdi
      if (h%gswmnmdi /= int_unread) 
     |  write(6,"(2x,'gswmnmdi      = ',i2,' (gswm flag)')") h%gswmnmdi
      if (h%gswmnmsdi /= int_unread) 
     |  write(6,"(2x,'gswmnmsdi     = ',i2,' (gswm flag)')") h%gswmnmsdi

      if (h%missing_value /= real_unread) 
     |  write(6,"(2x,'missing_value = ',e12.4,' (global attribute)')") 
     |  h%missing_value

      if (len_trim(h%ncep_ncfile) > 0)
     |  write(6,"(2x,'ncep_ncfile = ',a)") trim(h%ncep_ncfile)
      if (len_trim(h%gpi_ncfile) > 0)
     |  write(6,"(2x,'gpi_ncfile = ',a)") trim(h%gpi_ncfile)
      if (len_trim(h%imf_ncfile) > 0)
     |  write(6,"(2x,'imf_ncfile = ',a)") trim(h%imf_ncfile)
      if (len_trim(h%ecmwf_ncfile) > 0)
     |  write(6,"(2x,'ecmwf_ncfile = ',a)") trim(h%ecmwf_ncfile)
      if (len_trim(h%see_ncfile) > 0)
     |  write(6,"(2x,'see_ncfile = ',a)") trim(h%see_ncfile)

      if (len_trim(h%gswm_mi_di_ncfile) > 0)
     |  write(6,"(2x,'gswm_mi_di_ncfile  = ',a)") 
     |    trim(h%gswm_mi_di_ncfile)
      if (len_trim(h%gswm_mi_sdi_ncfile) > 0)
     |  write(6,"(2x,'gswm_mi_sdi_ncfile = ',a)") 
     |    trim(h%gswm_mi_sdi_ncfile)

      if (len_trim(h%gswm_nm_di_ncfile) > 0)
     |  write(6,"(2x,'gswm_nm_di_ncfile  = ',a)") 
     |    trim(h%gswm_nm_di_ncfile)
      if (len_trim(h%gswm_nm_sdi_ncfile) > 0)
     |  write(6,"(2x,'gswm_nm_sdi_ncfile = ',a)") 
     |    trim(h%gswm_nm_sdi_ncfile)

      write(6,"(2x,'nf3d_geo = ',i3,' (number of 3d geo fields)')")
     |  h%nf3d_geo
      write(6,"(2x,'nf2d_geo = ',i3,' (number of 2d geo fields)')")
     |  h%nf2d_geo
      write(6,"(2x,'nf3d_mag = ',i3,' (number of 3d mag fields)')")
     |  h%nf3d_mag
      write(6,"(2x,'nf2d_mag = ',i3,' (number of 2d mag fields)')")
     |  h%nf2d_mag
!
      if (associated(h%glon))
     |  write(6,"(/,'nlon=',i3,' glon = ',/,(6e12.4))") h%nlon,h%glon 
      if (associated(h%glat))
     |  write(6,"(/,'nlat=',i3,' glat = ',/,(6e12.4))") h%nlat,h%glat 
      if (associated(h%lev))
     |  write(6,"(/,'nlev=',i3,' lev = ',/,(6e12.4))") h%nlev,h%lev 
      if (associated(h%ilev))
     |  write(6,"(/,'nilev=',i3,' ilev = ',/,(6e12.4))") h%nilev,h%ilev 

      if (associated(h%gmlon))
     |  write(6,"(/,'mlon=',i3,' gmlon = ',/,(6e12.4))") h%mlon,h%gmlon 
      if (associated(h%gmlat))
     |  write(6,"(/,'mlat=',i3,' gmlat = ',/,(6e12.4))") h%mlat,h%gmlat 
      if (associated(h%gmlev))
     |  write(6,"(/,'mlev=',i3,' gmlev = ',/,(6e12.4))") h%mlev,h%gmlev 
      if (associated(h%gimlev))
     |  write(6,"(/,'imlev=',i3,' gimlev= ',/,(6e12.4))")
     |    h%imlev,h%gimlev 
!
      if (h%new) then
        write(6,"(/,'  This is a ""new"" history.')")
      else
        write(6,"(/,'  This is an ""old"" history.')")
      endif
!
! Print min,max of 3d geo fields:
      if (h%nf3d_geo > 0) then
        write(6,"(/)")
        do i=1,h%nf3d_geo
          call fminmaxspv(h%f3d_geo(i)%data(:,:,:),h%nlon*h%nlat*h%nlev,
     |      fmin,fmax,h%f3d_geo(i)%missing_value)
          if (h%f3d_geo(i)%levdimname=='ilev') then
            write(6,"(2x,'3d geo Field ',a,' min,max=',2e12.4,
     |        ' (interface levels)')") h%f3d_geo(i)%short_name,fmin,fmax
          else
            write(6,"(2x,'3d geo Field ',a,' min,max=',2e12.4,
     |        ' (midpoint levels)')") h%f3d_geo(i)%short_name,fmin,fmax
          endif
        enddo
      endif
!
! Print min,max of 3d mag fields:
      if (h%nf3d_mag > 0) then
        write(6,"(/)")
        do i=1,h%nf3d_mag
          call fminmaxspv(h%f3d_mag(i)%data(:,:,:),h%mlon*h%mlat*h%mlev,
     |      fmin,fmax,h%f3d_mag(i)%missing_value)
          write(6,"(2x,'3d mag Field ',a,' min,max=',2e12.4)")
     |      h%f3d_mag(i)%short_name,fmin,fmax
        enddo
      endif
!
! Print min,max of 2d geo fields:
      if (h%nf2d_geo > 0) then
        write(6,"(/)")
        do i=1,h%nf2d_geo
          call fminmaxspv(h%f2d_geo(i)%data(:,:),h%nlon*h%nlat,
     |      fmin,fmax,h%f2d_geo(i)%missing_value)
          write(6,"(2x,'2d geo Field ',a,' min,max=',2e12.4)")
     |      h%f2d_geo(i)%short_name,fmin,fmax
        enddo
      endif
!
! Print min,max of 2d mag fields:
      if (h%nf2d_mag > 0) then
        write(6,"(/)")
        do i=1,h%nf2d_mag
          call fminmaxspv(h%f2d_mag(i)%data(:,:),h%mlon*h%mlat,
     |      fmin,fmax,h%f2d_mag(i)%missing_value)
          write(6,"(2x,'2d mag Field ',a,' min,max=',2e12.4)")
     |      h%f2d_mag(i)%short_name,fmin,fmax
        enddo
      endif
      write(6,"(72('-'),/)")
      end subroutine hist_print
!-----------------------------------------------------------------------
      subroutine handle_ncerr(istat,msg)
      implicit none
!
! 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
!-----------------------------------------------------------------------
      subroutine fminmax(f,n,fmin,fmax)
      implicit none
!
! Return min and max of fields f(n) (between -1.e36,+1.e36)
!
! Args: 
      integer,intent(in) :: n
      real,intent(in) :: f(n)
      real,intent(out) :: fmin,fmax
!
! Local:
      integer :: i
!
      fmin = 1.e38
      fmax = -1.e38
      do i=1,n
!       fmin = amin1(f(i),fmin)
!       fmax = amax1(f(i),fmax)
        fmin = min(f(i),fmin)
        fmax = max(f(i),fmax)
      enddo
      end subroutine fminmax
!-------------------------------------------------------------------
      subroutine fminmaxspv(f,n,fmin,fmax,spv)
      implicit none
!
! Return min and max of fields f(n) (between -1.e36,+1.e36)
! Ignore any f(i)==spv.
!
! Args:
      integer,intent(in) :: n
      real,intent(in) :: f(n),spv
      real,intent(out) :: fmin,fmax
!
! Local:
      integer :: i
!
      fmin = 1.e36
      fmax = -1.e36
      do i=1,n
        if (f(i) /= spv) then
          fmin = min(f(i),fmin)
          fmax = max(f(i),fmax)
        endif
      enddo
      end subroutine fminmaxspv
!-----------------------------------------------------------------------
      subroutine ftype(h,ncid,varname,id,ndims,iddims,type)
!
! Given a variable to be read, its number of dimensions ndims, and
! the dimension id's iddims(ndims), return type as one of the following:
!
!   '3dgeo' A 3d field (lon,lat,lev) on geographic grid
!   '3dmag' A 3d field (lon,lat,lev) on magnetic grid
!   '2dgeo' A 2d field (lon,lat) on geographic grid
!   '2dmag' A 2d field (lon,lat) on magnetic grid
!   'other' None of the above
!
      implicit none
!
! Args: 
      type(history),intent(in) :: h
      character(len=*),intent(in) :: varname
      integer,intent(in) :: id,ndims,ncid
      integer,intent(in) :: iddims(ndims)
      character(len=8),intent(out) :: type
!
! Local:
      integer :: istat,i,len,idims(ndims)
!
      type = 'other   ' ! init
      do i=1,ndims 
        istat = nf_inq_dimlen(ncid,iddims(i),idims(i)) ! dim lengths
      enddo
!
      if (ndims == 4) then ! (lon,lat,lev,time)
        if (idims(1)==h%nlon.and.idims(2)==h%nlat.and.idims(3)==h%nlev
     |    .and.idims(4)==h%nhist) then
          type = '3dgeo   '
        endif
        if (idims(1)==h%mlon.and.idims(2)==h%mlat.and.idims(3)==h%mlev
     |    .and.idims(4)==h%nhist) then
          type = '3dmag   '
        endif

!       write(6,"('ftype: varname=',a,' id=',i3,' ndims=',i3,' iddims=',
!    |    4i4,' idims=',4i4,' type=',a)") varname(1:8),id,ndims,iddims,
!    |    idims,type
 
      elseif (ndims == 3) then ! (lon,lat,time)
        if (idims(1)==h%nlon.and.idims(2)==h%nlat.and.idims(3)==h%nhist) 
     |    then
          type = '2dgeo   '
        endif
        if (idims(1)==h%mlon.and.idims(2)==h%mlat.and.idims(3)==h%nhist) 
     |    then
          type = '2dmag   '
        endif

!       write(6,"('ftype: varname=',a,' id=',i3,' ndims=',i3,' iddims=',
!    |    3i4,' idims=',3i4,' type=',a)") varname(1:8),id,ndims,iddims,
!    |    idims,type

      endif

      end subroutine ftype
!-----------------------------------------------------------------------
      subroutine read_f3d(ncid,fname,h,geomag,itime)
!
! Args:
      integer,intent(in) :: ncid,itime
      character(len=*),intent(in) :: fname,geomag
      type(history),intent(inout) :: h
!
! Local:
      integer,parameter :: mxdims=20 ! assume <= mxdims dimensions
      integer :: istat,id,itype,ndims,iddims(mxdims),natts,len,if3d
      character(len=80) :: varname,attname,levname
      character(len=shortname_len) :: short_name
      character(len=longname_len) :: long_name
      character(len=units_len) :: units
      type(field_3d) :: f3d
      integer :: start_4d(4),count_4d(4),ier
      real :: missing_value
      
      long_name = ' '
      units = ' '
      istat = nf_inq_varid(ncid,fname,id)

      istat = nf_inq_var(ncid,id,varname,itype,ndims,iddims,natts)
      istat = nf_inq_dimname(ncid,iddims(3),levname)
!     write(6,"('read_f3d: geomag=',a,' varname=',a,' ndims=',i3,
!    |  ' iddims=',4i4,' 3rd levname=',a)") geomag,trim(varname),ndims,
!    |  iddims(1:ndims),trim(levname)

      istat = nf_get_att_text(ncid,id,'long_name',long_name)
      istat = nf_get_att_text(ncid,id,'units',units)
      istat = nf_get_att_double(ncid,id,'missing_value',missing_value)
      if (istat /= 0) missing_value = real_unread

      start_4d(1:3) = 1
      start_4d(4) = itime
      count_4d(4) = 1

      if (geomag=='geo') then
        if3d = h%nf3d_geo
        call init_f3d(h%f3d_geo(if3d))
        h%f3d_geo(if3d)%long_name  = long_name
        h%f3d_geo(if3d)%short_name = fname
        h%f3d_geo(if3d)%units      = units
        h%f3d_geo(if3d)%missing_value = missing_value
        h%f3d_geo(if3d)%geo = .true.
        h%f3d_geo(if3d)%mag = .false.
        allocate(h%f3d_geo(if3d)%data(h%nlon,h%nlat,h%nlev),stat=ier)
        h%f3d_geo(if3d)%id1 = h%nlon
        h%f3d_geo(if3d)%id2 = h%nlat
        h%f3d_geo(if3d)%id3 = h%nlev
        count_4d(1) = h%nlon
        count_4d(2) = h%nlat
        count_4d(3) = h%nlev
        h%f3d_geo(if3d)%levdimname = levname
        if (trim(levname)=='ilev') then ! interfaces
          h%f3d_geo(if3d)%id3 = h%nilev
          count_4d(3) = h%nilev
        endif
        istat = nf_get_vara_double(ncid,id,start_4d,count_4d,
     |    h%f3d_geo(if3d)%data(:,:,:))
        h%f3d_geo(if3d)%type = ' '
        if (trim(fname) == 'O1'  .or. trim(fname) == 'O2'  .or. 
     |      trim(fname) == 'OX'  .or. trim(fname) == 'O3'  .or. 
     |      trim(fname) == 'N4S' .or. trim(fname) == 'NOZ' .or. 
     |      trim(fname) == 'CO'  .or. trim(fname) == 'CO2' .or. 
     |      trim(fname) == 'H2O' .or. trim(fname) == 'H2'  .or. 
     |      trim(fname) == 'HOX' .or. trim(fname) == 'OP'  .or. 
     |      trim(fname) == 'CH4' .or. trim(fname) == 'AR'  .or. 
     |      trim(fname) == 'HE'  .or. trim(fname) == 'NAT' .or. 
     |      trim(fname) == 'O21D'.or. trim(fname) == 'NO2' .or. 
     |      trim(fname) == 'NO'  .or. trim(fname) == 'OH'  .or. 
     |      trim(fname) == 'HO2' .or. trim(fname) == 'HO2' .or. 
     |      trim(fname) == 'H'   .or. trim(fname) == 'N2D' .or. 
     |      trim(fname) == 'O2P') h%f3d_geo(if3d)%type = 'DENSITY'

!       write(6,"('Read 3d ',a,' if3d=',i3,' Field ',a,':')") 
!    |    geomag,if3d,trim(h%f3d_geo(if3d)%short_name)
!       write(6,"('  long_name:  ',a)") trim(h%f3d_geo(if3d)%long_name)
!       write(6,"('  units:      ',a)") trim(h%f3d_geo(if3d)%units)
!       write(6,"('  levdimname: ',a)") trim(h%f3d_geo(if3d)%levdimname)
!       write(6,"('  missing value: ',e12.4)") 
!    |    h%f3d_geo(if3d)%missing_value
!       write(6,"('  data min,max=',2e12.4)") 
!    |    minval(h%f3d_geo(if3d)%data),
!    |    maxval(h%f3d_geo(if3d)%data)

      else                                    ! 3d mag field
        if3d = h%nf3d_mag
        call init_f3d(h%f3d_mag(if3d))
        h%f3d_mag(if3d)%long_name  = long_name
        h%f3d_mag(if3d)%short_name = fname
        h%f3d_mag(if3d)%units      = units
        h%f3d_mag(if3d)%geo = .false.
        h%f3d_mag(if3d)%mag = .true.
        allocate(h%f3d_mag(if3d)%data(h%mlon,h%mlat,h%mlev),stat=ier)
        h%f3d_mag(if3d)%id1 = h%mlon
        h%f3d_mag(if3d)%id2 = h%mlat
        h%f3d_mag(if3d)%id3 = h%mlev
        count_4d(1) = h%mlon
        count_4d(2) = h%mlat
        count_4d(3) = h%mlev
        h%f3d_mag(if3d)%levdimname = levname
        if (trim(levname)=='imlev') then ! interfaces
          h%f3d_mag(if3d)%id3 = h%imlev
          count_4d(3) = h%imlev
        endif
        istat = nf_get_vara_double(ncid,id,start_4d,count_4d,
     |    h%f3d_mag(if3d)%data(:,:,:))

!       write(6,"('Read 3d ',a,' if3d=',i3,' Field ',a,':')") 
!    |    geomag,if3d,trim(h%f3d_mag(if3d)%short_name)
!       write(6,"('  long_name:  ',a)") trim(h%f3d_mag(if3d)%long_name)
!       write(6,"('  units:      ',a)") trim(h%f3d_mag(if3d)%units)
!       write(6,"('  levdimname: ',a)") trim(h%f3d_mag(if3d)%levdimname)
!       write(6,"('  data min,max=',2e12.4)") 
!    |    minval(h%f3d_mag(if3d)%data),
!    |    maxval(h%f3d_mag(if3d)%data)

      endif
      end subroutine read_f3d
!-----------------------------------------------------------------------
      subroutine read_f2d(ncid,fname,h,geomag,itime)
!
! Args:
      integer,intent(in) :: ncid,itime
      character(len=*),intent(in) :: fname,geomag
      type(history),intent(inout) :: h
!
! Local:
      integer,parameter :: mxdims=20 ! assume <= mxdims dimensions
      integer :: istat,id,itype,ndims,iddims(mxdims),natts,len,if2d
      character(len=80) :: varname,attname
      character(len=shortname_len) :: short_name
      character(len=longname_len) :: long_name
      character(len=units_len) :: units
      integer :: start_3d(3),count_3d(3),ier
      real :: missing_value
      
      long_name = ' '
      units = ' '
      istat = nf_inq_varid(ncid,fname,id)
      istat = nf_get_att_text(ncid,id,'long_name',long_name)
      istat = nf_get_att_text(ncid,id,'units',units)
      istat = nf_get_att_double(ncid,id,'missing_value',missing_value)
      if (istat /= 0) missing_value = real_unread

      start_3d(1:2) = 1
      start_3d(3) = itime
      count_3d(3) = 1

      if (geomag=='geo') then
        if2d = h%nf2d_geo
        call init_f2d(h%f2d_geo(if2d))
        h%f2d_geo(if2d)%long_name  = long_name
        h%f2d_geo(if2d)%short_name = fname
        h%f2d_geo(if2d)%units      = units
        h%f2d_geo(if2d)%missing_value = missing_value
        h%f2d_geo(if2d)%geo = .true.
        h%f2d_geo(if2d)%mag = .false.
        allocate(h%f2d_geo(if2d)%data(h%nlon,h%nlat),stat=ier)
        h%f2d_geo(if2d)%id1 = h%nlon
        h%f2d_geo(if2d)%id2 = h%nlat
        count_3d(1) = h%nlon
        count_3d(2) = h%nlat
        istat = nf_get_vara_double(ncid,id,start_3d,count_3d,
     |    h%f2d_geo(if2d)%data(:,:))

!       write(6,"('Read 2d ',a,' if2d=',i3,' Field ',a,':')") 
!    |    geomag,if2d,trim(h%f2d_geo(if2d)%short_name)
!       write(6,"('  long_name: ',a)") trim(h%f2d_geo(if2d)%long_name)
!       write(6,"('  units:     ',a)") trim(h%f2d_geo(if2d)%units)
!       write(6,"('  missing_value: ',e12.4)") 
!    |    h%f2d_geo(if2d)%missing_value
!       write(6,"('  data min,max=',2e12.4)") 
!    |    minval(h%f2d_geo(if2d)%data),
!    |    maxval(h%f2d_geo(if2d)%data)

      else
        if2d = h%nf2d_mag
        call init_f2d(h%f2d_mag(if2d))
        h%f2d_mag(if2d)%long_name  = long_name
        h%f2d_mag(if2d)%short_name = fname
        h%f2d_mag(if2d)%units      = units
        h%f2d_mag(if2d)%geo = .false.
        h%f2d_mag(if2d)%mag = .true.
        allocate(h%f2d_mag(if2d)%data(h%mlon,h%mlat),stat=ier)
        h%f2d_geo(if2d)%id1 = h%mlon
        h%f2d_geo(if2d)%id2 = h%mlat
        count_3d(1) = h%mlon
        count_3d(2) = h%mlat
        istat = nf_get_vara_double(ncid,id,start_3d,count_3d,
     |    h%f2d_mag(if2d)%data(:,:))

!       write(6,"('Read 2d ',a,' if2d=',i3,' Field ',a,':')") 
!    |    geomag,if2d,trim(h%f2d_mag(if2d)%short_name)
!       write(6,"('  long_name: ',a)") trim(h%f2d_mag(if2d)%long_name)
!       write(6,"('  units:     ',a)") trim(h%f2d_mag(if2d)%units)
!       write(6,"('  missing_value: ',e12.4)") 
!    |    h%f2d_mag(if2d)%missing_value
!       write(6,"('  data min,max=',2e12.4)") 
!    |    minval(h%f2d_mag(if2d)%data),
!    |    maxval(h%f2d_mag(if2d)%data)

      endif
      end subroutine read_f2d
!-----------------------------------------------------------------------
      subroutine init_f2d(f)
      type(field_2d),intent(out) :: f
      f%short_name = ' '
      f%long_name = ' '
      f%units = ' '
      end subroutine init_f2d
!-----------------------------------------------------------------------
      subroutine init_f3d(f)
      type(field_3d),intent(out) :: f
      f%short_name = ' '
      f%long_name = ' '
      f%units = ' '
      f%levdimname = ' '
      end subroutine init_f3d
!-----------------------------------------------------------------------
      end module readhist_module