program drive_ncinfo
  !
  ! Requirements:
  ! base state file;
  ! 1 to N input files;
  ! one auxiliary grid information file.
  ! All in netCDF format.
  !
  ! Argument one is "history" or "singlevar" - REQUIRED
  !
  ! Argument two is the "base" ocn file, ie, the one from which
  !  sea level change is to be computed
  ! If arg0 is "history", then
  !  Arguments 3 to N are the ocean files for which SLC should be computed (up to 1000 files)
  ! If arg0 is "singlevar", then
  !  Argument  3 is the TEMP data
  !  Argument  4 is the SALT data
  !  Argument  5 is the SSH data
  !  Argument  6 is the RHO data
  !
  ! Example history input file:
  ! history
  ! B06.62.ocn.1000-1100Ave.nc
  ! B06.62.ocn.1000.nc
  ! B06.62.ocn.1001.nc
  ! [...]
  ! B06.62.ocn.1099.nc
  ! B06.62.ocn.1100.nc
  !
  ! Example singlevar input file:
  ! singlevar
  ! b40.1850.track1.1deg.006.pop.h.0800-ave-1299.nc
  ! b40.20th.track1.1deg.008.pop.h.TEMP.192001-192912.nc
  ! b40.20th.track1.1deg.008.pop.h.SALT.192001-192912.nc
  ! b40.20th.track1.1deg.008.pop.h.SSH.192001-192912.nc
  ! b40.20th.track1.1deg.008.pop.h.RHO.192001-192912.nc
  implicit none
  !
  call get_ncinfo
end program drive_ncinfo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_ncinfo
  !
  ! netCDF I/O subroutines
  !
  use counters_netcdf_jfl
  use interfaces_netcdf_jfl
  use definitions_netcdf_jfl
  use kinds_mod
  !
  implicit none
  !
  integer :: i,j,k,l,n,ier,length,io_status,ifile,nfile
  integer :: cnt_hist,cnt_temp,cnt_salt,cnt_ssh,cnt_rho
  integer :: cnt_base,cnt_grid,cnt_slc
  character(len=132) :: base_ocean_dataset
  character(len=132) :: temp_dataset,salt_dataset,rho_dataset,ssh_dataset
  character(len=132),dimension(1000) :: ocn_history
  !
  exists = .false.
  ishist = .false.
  issvar = .false.
  ispcm1 = .false.
  isccsm3 = .false.
  isccsm4 = .false.
  !
  read (*,*) dataset_type
  if (dataset_type == "history") then
     write(*,*) 'History file format data.'
     ishist = .true.
  elseif (dataset_type == "singlevar") then
     write(*,*) 'Single variable format data.'
     issvar = .true.
  else
     write(*,*) 'Unknown dataset format. DIE!'
     call abort()
  endif
  !
  read (*,*) base_ocean_dataset
  call check_exists(base_ocean_dataset)
  write(*,'(a,'' is base ocean state.'')') trim(base_ocean_dataset)
  !
  if (ishist) then
     io_status = 0
     ifile = 1
     do while (io_status == 0)
        read(*,*,iostat=io_status) ocn_history(ifile)
        if (io_status == 0) then
           call check_exists(ocn_history(ifile))
           ifile = ifile + 1
        endif
     enddo
     nfile = ifile - 1
  endif
  if (issvar) then
     read(*,*) temp_dataset
     read(*,*) salt_dataset
     read(*,*) ssh_dataset
     read(*,*) rho_dataset
     call check_exists(temp_dataset)
     call check_exists(salt_dataset)
     call check_exists(ssh_dataset)
     call check_exists(rho_dataset)
  endif
  !
  ! Get base ocean data
  !
  call open_cdf(ncid_base,trim(base_ocean_dataset),.true.)
  call get_dims(ncid_base)
  call get_vars(ncid_base)
  !
  ! Get base ocean dimensions
  !
  do n=1,dim_counter
     length = len_trim(dim_info(n)%name)
     if((dim_info(n)%name(:length).eq.'lat').or.(dim_info(n)%name(:length).eq.'nlat')) then
        nlat_base = dim_info(n)%length
     elseif((dim_info(n)%name(:length).eq.'lon').or.(dim_info(n)%name(:length).eq.'nlon')) then
        nlon_base = dim_info(n)%length
     elseif((dim_info(n)%name(:length).eq.'lev').or.(dim_info(n)%name(:length).eq.'z_t')) then
        nlev_base = dim_info(n)%length
     elseif(dim_info(n)%name(:length).eq.'time') then
        ntim_base = dim_info(n)%length
     endif
  enddo
  !
  if (nlev_base.eq.32) ispcm1  = .true.
  if (nlev_base.eq.40) isccsm3 = .true.
  if (nlev_base.eq.60) isccsm4 = .true.
  !
  ! If history format, step through ocean datasets
  !
  if (ishist) then
     do ifile = 1,nfile
        slc_dataset = 'DSL.'//trim(ocn_history(ifile))
        dim_counter = 0
        var_counter = 0
        call open_cdf(ncid_hist,trim(ocn_history(ifile)),.true.)
        call get_dims(ncid_hist)
        call get_vars(ncid_hist)
        !
        ! Get ocean dataset dimensions
        do n=1,dim_counter
           length = len_trim(dim_info(n)%name)
           if((dim_info(n)%name(:length).eq.'lat').or.(dim_info(n)%name(:length).eq.'nlat')) then
              nlat_hist = dim_info(n)%length
           elseif((dim_info(n)%name(:length).eq.'lon').or.(dim_info(n)%name(:length).eq.'nlon')) then
              nlon_hist = dim_info(n)%length
           elseif((dim_info(n)%name(:length).eq.'lev').or.(dim_info(n)%name(:length).eq.'z_t')) then
              nlev_hist = dim_info(n)%length
           elseif(dim_info(n)%name(:length).eq.'time') then
              ntim_hist = dim_info(n)%length
           endif
        enddo
        !
        ! Get coordinate variables
        !
        if (ispcm1) then
           call read_var(ncid_hist,'lon',ocn_lon)
           call read_var(ncid_hist,'lat',ocn_lat)
        else
           do i=1,nlon_hist
              ocn_lon(i)=i
           end do
           do i=1,nlat_hist
              ocn_lat(i)=i
           end do
        endif

     end do
