!This file extracted from ASPEN rdmag.for, July 2003
!-------------------------------------------------------------------
      subroutine rdfield(imaxg,jmaxg,alatm,alonm)
!
! Read netcdf magnetic data file (mss /TGCM/data/magdat.nc)
! This file was written from the original cray-blocked file
!   /ECRIDLEY/ECR90/ECRMG6 by code in ~foster/tgcm/mkmag)
!
      use netcdf_module,only: nc_open,nc_close,handle_ncerr
      implicit none
      include "netcdf.inc"
!
! Args:
!
      integer,intent(in) :: imaxg,jmaxg
      real,intent(out) :: alatm(imaxg,jmaxg), alonm(imaxg,jmaxg)
! Local:
      character(len=*), parameter :: dskfile='TGCM.data.magdat.nc'
      integer :: zjmxp2
      integer :: istat,ncid
      integer :: ids1(1),ids2(2),ids3(3),ids4(4)
      integer :: id_zjmx,id_zimxp1,id_zjmxp2,id_dim2,id_dim3,
     |  id_imaxmp,id_jmaxm,id_dim4
      integer :: idv_alatm,idv_alonm,idv_xb,idv_yb,idv_zb,idv_bmod,
     |  idv_dmlat
      integer :: start_1d(1),count_1d(1),start_2d(2),count_2d(2),
     |  start_3d(3),count_3d(3),start_4d(4),count_4d(4)
      character(len=8) :: dimname
      real :: fmin,fmax
      real :: alattmp(imaxg,jmaxg+2), alontmp(imaxg,jmaxg+2)
!
      write(6,"(/,72('-'))")
      write(6,"('RDMAG: read magnetic field data file:')")
      zjmxp2=jmaxg+2
!
! Open the netcdf dataset:
      call nc_open(ncid,dskfile,'OLD','READ')
      if (ncid <= 0) then
        write(6,"(/,'>>> nc_rdmag: error opening netcdf mag data ',
     |    'file ',a)") trim(dskfile)
        stop 'nc_rdmag'
!     else
!       write(6,"('nc_rdmag: opened netcdf mag data file ',a,
!    |    ' ncid=',i8)") trim(dskfile),ncid
      endif
!
! Check dimensions:
      call checkdim(ncid,"zjmx"   ,jmaxg)
      call checkdim(ncid,"zimxp1" ,imaxg)
      call checkdim(ncid,"zjmxp2" ,zjmxp2)
!
! Read variables for fieldz.h:
! 2-d doubles (zimxp1,0:zjmxp1):
      start_2d(:) = 1
      count_2d(1) = imaxg
      count_2d(2) = zjmxp2
      call rd2dfld(ncid,'ALATM ',idv_alatm ,start_2d,count_2d,alattmp )
      call rd2dfld(ncid,'ALONM ',idv_alonm ,start_2d,count_2d,alontmp )
      !call rd2dfld(ncid,'XB    ',idv_xb    ,start_2d,count_2d,xb    )
      !call rd2dfld(ncid,'YB    ',idv_yb    ,start_2d,count_2d,yb    )
      !call rd2dfld(ncid,'ZB    ',idv_zb    ,start_2d,count_2d,zzb   )
      !call rd2dfld(ncid,'BMOD  ',idv_bmod  ,start_2d,count_2d,bmod  )

      alatm(:,1:jmaxg) = alattmp(:,2:jmaxg+1)
      alonm(:,1:jmaxg) = alontmp(:,2:jmaxg+1)
!
! Close the dataset:
      call nc_close(ncid)
      write(6,"('Completed read of magnetic field data file.')")
      write(6,"(72('-'),/)")
      end subroutine rdfield
!-------------------------------------------------------------------
      subroutine checkdim(ncid,dimname,iparam)
!
! Get length of dimension "dimname". If this length is not equal
! to iparam, stop with error message.
!
      use netcdf_module,only: handle_ncerr
      implicit none
      include "netcdf.inc"
!
! Args:
      integer,intent(in) :: ncid,iparam
      character(len=*),intent(in) :: dimname
!
! Local:
      integer :: istat,iddim,len
      character(len=80) :: char80
!
! Get dim id:
      istat = nf_inq_dimid(ncid,dimname,iddim)
      if (istat /= NF_NOERR) then
        write(char80,"('nc_rdmag: error getting dim id for ',
     |    a)") dimname
        call handle_ncerr(istat,char80)
      endif
!
! Get dim length:
      istat = nf_inq_dimlen(ncid,iddim,len)
      if (istat /= NF_NOERR) then
        write(char80,"('nc_rdmag: error getting length of ',
     |    'dimension ',a)") dimname
        call handle_ncerr(istat,char80)
      endif
!
! Compare with iparam:
      if (len /= iparam) then
        write(6,"(/,'>>> nc_rdmag: unexpected length for ',
     |    'dimension ',a)") dimname
        write(6,"('  length read = ',i3,' should be = ',i3)")
     |    len,iparam
        stop 'nc_rdmag'
      endif
      end subroutine checkdim
!-------------------------------------------------------------------
      subroutine rd2dfld(ncid,name,idvout,start_2d,count_2d,var)
      use netcdf_module,only: handle_ncerr
      implicit none
      include "netcdf.inc"
!
! Read 2-d double array from ncid to var:
!
! Args:
      integer,intent(in) :: ncid,start_2d(2),count_2d(2)
      character(len=*),intent(in) :: name
      integer,intent(out) :: idvout
      real,intent(out) :: var(count_2d(1),count_2d(2))
!
! Local:
      integer :: istat
      character(len=80) :: char80
      real :: fmin,fmax
!
      istat = nf_inq_varid(ncid,name,idvout)
      istat = nf_get_vara_double(ncid,idvout,start_2d,count_2d,var)
      write(char80,"('Error return from nf_get_vara_double for var',
     |  a)") name
      if (istat /= NF_NOERR) call handle_ncerr(istat,char80)
!     call fminmax(var,count_2d(1)*count_2d(2),fmin,fmax)
!     write(6,"('rd2dfld: ',a,' fmin,max=',2e12.4)") name,fmin,fmax
      end subroutine rd2dfld