program dat_to_nc
  use read_data
  implicit none
#include <netcdf.inc>
!
! 4/15/08, B. Foster:
!   Read .dat files containing data for Weimer 2005 model, and 
!   write equivalent netcdf file for reading by tiegcm. These
!   ascii text files were made from original idl .xdr files
!   provided by Dan Weimer, April, 2008.
!
! To compile on hao linux:
! pgf90 -o dat_to_nc -r8 read_data.F90 dat_to_nc.F90 -I /opt/local/include -L/opt/local/lib -lnetcdf
!
  integer,parameter :: n_schfits=d1_pot, n_alschfits=d2_pot
  real :: epot_schfits(n_schfits,csize), epot_alschfits(n_alschfits,csize)
  real :: bpot_schfits(n_schfits,csize), bpot_alschfits(n_alschfits,csize)
  real :: ex_epot(2),ex_bpot(2)
  integer,parameter :: n1_scha=d1_scha, n2_scha=d2_scha, n3_scha=d3_scha,&
    nex=2
  character(len=32) :: infile,outfile
  character(len=1024) :: bigchar
  integer :: istat,ncid,id_na,id_nb,id_nex,id_n1_scha,id_n2_scha,id_n3_scha,&
    id_maxk_scha,id_maxm_scha,id_csize,id_n_schfits,id_n_alschfits,id_maxm_pot,&
    id_maxl_pot
  integer :: ids1(1),idv_bndya,idv_bndyb,idv_ex_bndy,idv_th0s,idv_ab,&
    idv_ls,idv_ms,idv_ex_epot,idv_ex_bpot
  integer :: ids3(3),idv_allnkm
  integer :: ids2(2),idv_epot_schfits,idv_epot_alschfits,&
                     idv_bpot_schfits,idv_bpot_alschfits
!
! Read ascii data file W05scEpot.dat:
!   integer,parameter :: csize=28, d1_pot=15, d2_pot=18
!   integer :: ab(csize), ls(csize), ms(csize)
!   real :: alschfits(d2_pot,csize), schfits(d1_pot,csize), ex_pot(2)
!   integer :: maxl_pot,maxm_pot
!
  infile = ' '
  write(infile,"('W05scEpot.dat')")
  call read_potential(trim(infile))
  epot_schfits = schfits
  epot_alschfits = alschfits
  ex_epot = ex_pot
  write(6,"('dat_to_nc: epot_schfits min,max=',2e12.4)") &
    minval(epot_schfits),maxval(epot_schfits)
  write(6,"('dat_to_nc: epot_alschfits min,max=',2e12.4)") &
    minval(epot_alschfits),maxval(epot_alschfits)
!
! Read ascii data file W05scBpot.dat:
!
  infile = ' '
  write(infile,"('W05scBpot.dat')")
  call read_potential(trim(infile))
  bpot_schfits = schfits
  bpot_alschfits = alschfits
  bpot_schfits = schfits
  bpot_alschfits = alschfits
  ex_bpot = ex_pot
  write(6,"('dat_to_nc: bpot_schfits min,max=',2e12.4)") &
    minval(bpot_schfits),maxval(bpot_schfits)
  write(6,"('dat_to_nc: bpot_alschfits min,max=',2e12.4)") &
    minval(bpot_alschfits),maxval(bpot_alschfits)
!
! Read ascii data file SCHAtable.dat:
!   integer,parameter :: d1_scha=19, d2_scha=7, d3_scha=68
!   real :: allnkm(d1_scha,d2_scha,d3_scha)
!   integer :: maxk_scha,maxm_scha
!   real :: th0s(d3_scha)
!
  infile = ' '
  write(infile,"('SCHAtable.dat')")
  call read_schatable(trim(infile))
!
! Read ascii data file W05scBndy.dat:
! integer,parameter :: na=6, nb=7
! real :: bndya(na),bndyb(nb),ex_bndy(2)
!
  infile = ' '
  write(infile,"('W05scBndy.dat')")
  call read_bndy(trim(infile))
!
! Open new nc output file:
!
  outfile = ' '
  write(outfile,"('w05sc.nc')")
  istat = nf_create(trim(outfile),NF_CLOBBER,ncid)
  write(6,"(/,'Opened output file ',a)") trim(outfile)
!
! Define dimensions:
!
  istat = nf_def_dim(ncid,"na",na,id_na)
  istat = nf_def_dim(ncid,"nb",nb,id_nb)
  istat = nf_def_dim(ncid,"nex",nex,id_nex)
  istat = nf_def_dim(ncid,"n1_scha",n1_scha,id_n1_scha)
  istat = nf_def_dim(ncid,"n2_scha",n2_scha,id_n2_scha)
  istat = nf_def_dim(ncid,"n3_scha",n3_scha,id_n3_scha)
  istat = nf_def_dim(ncid,"maxk_scha",maxk_scha,id_maxk_scha)
  istat = nf_def_dim(ncid,"maxm_scha",maxm_scha,id_maxm_scha)
  istat = nf_def_dim(ncid,"maxl_pot",maxl_pot,id_maxl_pot)
  istat = nf_def_dim(ncid,"maxm_pot",maxm_pot,id_maxm_pot)
  istat = nf_def_dim(ncid,"csize",csize,id_csize)
  istat = nf_def_dim(ncid,"n_schfits",n_schfits,id_n_schfits)
  istat = nf_def_dim(ncid,"n_alschfits",n_alschfits,id_n_alschfits)
!
! Define bndy vars:
!
  ids1(1) = id_na
  istat = nf_def_var(ncid,"bndya",NF_DOUBLE,1,ids1,idv_bndya)
  ids1(1) = id_nb
  istat = nf_def_var(ncid,"bndyb",NF_DOUBLE,1,ids1,idv_bndyb)
  ids1(1) = id_nex
  istat = nf_def_var(ncid,"ex_bndy",NF_DOUBLE,1,ids1,idv_ex_bndy)
!
! Define scha table vars:
!
  ids1(1) = id_n3_scha
  istat = nf_def_var(ncid,"th0s",NF_DOUBLE,1,ids1,idv_th0s)
  ids3(1) = id_n1_scha
  ids3(2) = id_n2_scha
  ids3(3) = id_n3_scha
  istat = nf_def_var(ncid,"allnkm",NF_DOUBLE,3,ids3,idv_allnkm)
!
! Define epot and bpot vars:
!
! int ab(csize)
!
  ids1(1) = id_csize
  istat = nf_def_var(ncid,"ab",NF_INT,1,ids1,idv_ab)
  write(bigchar,"('Denotes if coef. is for a sine or cosine term')")
  istat = nf_put_att_text(ncid,idv_ab,'description',len_trim(bigchar),&
    trim(bigchar))
!
! int ls(csize)
!
  istat = nf_def_var(ncid,"ls",NF_INT,1,ids1,idv_ls)
  write(bigchar,"('Array of L values')")
  istat = nf_put_att_text(ncid,idv_ls,'description',len_trim(bigchar),&
    trim(bigchar))
!
! int ms(csize)
!
  istat = nf_def_var(ncid,"ms",NF_INT,1,ids1,idv_ms)
  write(bigchar,"('Array of M values')")
  istat = nf_put_att_text(ncid,idv_ms,'description',len_trim(bigchar),&
    trim(bigchar))
!
! real ex_epot(nex),ex_bpot(nex)
!
  ids1(1) = id_nex
  istat = nf_def_var(ncid,"ex_epot",NF_DOUBLE,1,ids1,idv_ex_epot)
  write(bigchar,&
   "('Constants for an exponential equation for saturation effect in Electric Potential model')")
  istat = nf_put_att_text(ncid,idv_ex_epot,'description',len_trim(bigchar),&
    trim(bigchar))
!
  ids1(1) = id_nex
  istat = nf_def_var(ncid,"ex_bpot",NF_DOUBLE,1,ids1,idv_ex_bpot)
  write(bigchar,&
   "('Constants for an exponential equation for saturation effect in Magnetic Potential model')")
  istat = nf_put_att_text(ncid,idv_ex_bpot,'description',len_trim(bigchar),&
    trim(bigchar))
!
  ids2(1) = id_n_schfits
  ids2(2) = id_csize
  istat = nf_def_var(ncid,"epot_schfits",NF_DOUBLE,2,ids2,idv_epot_schfits)
  write(bigchar,"('Spherical Cap Harmonic coefficients for Electric Potential model')")
  istat = nf_put_att_text(ncid,idv_epot_schfits,"Description",&
    len_trim(bigchar),trim(bigchar))
!
  istat = nf_def_var(ncid,"bpot_schfits",NF_DOUBLE,2,ids2,idv_bpot_schfits)
  write(bigchar,"('Spherical Cap Harmonic coefficients for Magnetic Potential model')")
  istat = nf_put_att_text(ncid,idv_bpot_schfits,"Description",&
    len_trim(bigchar),trim(bigchar))
!
  ids2(1) = id_n_alschfits
  istat = nf_def_var(ncid,"epot_alschfits",NF_DOUBLE,2,ids2,idv_epot_alschfits)
  write(bigchar,"('Spherical Cap Harmonic coefficients for Electric Potential model',&
    &' with optional AL indices')")
  istat = nf_put_att_text(ncid,idv_epot_alschfits,"Description",&
    len_trim(bigchar),trim(bigchar))
!
  istat = nf_def_var(ncid,"bpot_alschfits",NF_DOUBLE,2,ids2,idv_bpot_alschfits)
  write(bigchar,"('Spherical Cap Harmonic coefficients for Magnetic Potential model',&
    &' with optional AL indices')")
  istat = nf_put_att_text(ncid,idv_bpot_alschfits,"Description",&
    len_trim(bigchar),trim(bigchar))
!
!
! Global file attributes:
!
  bigchar = ' '
  write(bigchar,"('Ben Foster (NCAR/HAO) foster@ucar.edu')")
  istat = nf_put_att_text(ncid,NF_GLOBAL,"Author",len_trim(bigchar),trim(bigchar))
!
  bigchar = ' '
  write(bigchar,"('April 15, 2008')")
  istat = nf_put_att_text(ncid,NF_GLOBAL,"Creation_date",len_trim(bigchar),&
    trim(bigchar))
!
  bigchar = ' '
  write(bigchar,&
    "('W05scBndy.xdr, SCHAtable.xdr, W05scEpot.xdr, W05scBpot.xdr')")
  istat = nf_put_att_text(ncid,NF_GLOBAL,"Original_IDL_save_files",&
    len_trim(bigchar),trim(bigchar))
!
  bigchar = ' '
  write(bigchar,&
    "('Data required by 2005 version of Electric and Magnetic Potential Models',&
     &' by Dan Weimer')")
  istat = nf_put_att_text(ncid,NF_GLOBAL,"Description",len_trim(bigchar),&
    trim(bigchar))
!
  bigchar = ' '
  write(bigchar,&
    "('Weimer, D. R., Predicting Surface Geomagnetic Variations Using Ionospheric',&
      &' Electrodynamic Models, J. Geophys. Res., 110, A12307, doi:10.1029/',&
      &'2005JA011270, 2005.')")
  istat = nf_put_att_text(ncid,NF_GLOBAL,"Reference",len_trim(bigchar),&
    trim(bigchar))
!
! Exit define mode:
!
  istat = nf_enddef(ncid)
!
! Write bndy data:
!
  istat = nf_put_var_double(ncid,idv_bndya,bndya)
  istat = nf_put_var_double(ncid,idv_bndyb,bndyb)
  istat = nf_put_var_double(ncid,idv_ex_bndy,ex_bndy)
!
! Write scha table data:
!
  istat = nf_put_var_double(ncid,idv_th0s,th0s)
  istat = nf_put_var_double(ncid,idv_allnkm,allnkm)
!
! Write epot and bpot data:
!
  istat = nf_put_var_int(ncid,idv_ab,ab)
  istat = nf_put_var_int(ncid,idv_ls,ls)
  istat = nf_put_var_int(ncid,idv_ms,ms)
  istat = nf_put_var_double(ncid,idv_ex_epot,ex_epot)
  istat = nf_put_var_double(ncid,idv_ex_bpot,ex_bpot)
  istat = nf_put_var_double(ncid,idv_epot_schfits,epot_schfits)
  istat = nf_put_var_double(ncid,idv_bpot_schfits,bpot_schfits)
  istat = nf_put_var_double(ncid,idv_epot_alschfits,epot_alschfits)
  istat = nf_put_var_double(ncid,idv_bpot_alschfits,bpot_alschfits)
!
! Close output file:
!
  istat = nf_close(ncid)
  write(6,"('Closed output file ',a)") trim(outfile)
end program dat_to_nc
!-----------------------------------------------------------------------
      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
!-----------------------------------------------------------------------