module read_ncfile
!
! Read netcdf input file.
!
  use netcdf
  use shr_kind_mod ,only: r8 => shr_kind_r8            ! 8-byte reals
  use params       ,only: iulog,finit,fillvalue,pi,dtr,kbotdyn,zbotdyn
  use my_mpi       ,only: lev0,lev1,lon0,lon1,lat0,lat1
  use namelist     ,only: ctpoten_user => ctpoten
  use geogrid      ,only: &
    ntime,nlon,nlat,nlev,nilev,time,glat,glon,ylatg,ylong,&
    dlatg,dlong,zlev,zilev,lev_sequence,lon_sequence,     &
    jspole,jnpole,nlonp1,nlonp2,nlatp1
  use maggrid      ,only: zmlev,nmlev

  implicit none
  save

  integer :: ncid                ! netcdf file id
  integer,parameter :: nf=8      ! number of input fields
  character(len=8) :: fnames(nf) ! field names as on input file
!
! Inputs to edynamo on 3d distributed geographic subdomains 
! (lev0:lev1,lon0:lon1,lat0:lat1)
! These are read from the input history file for each time
!
  real(r8),allocatable,dimension(:,:,:),save :: & 
    tn,    & ! neutral temperature (deg K)
    un,    & ! neutral zonal wind (cm/s)
    vn,    & ! neutral meridional wind (cm/s)
    omega, & ! vertical motion (/sec)
    zpot,  & ! geopotential height (cm)
    barm,  & ! mean molecular weight
    ped,   & ! pedersen conductivity
    hall     ! hall conductivity
  real(r8) :: year,ut,ctpoten
  integer  :: iday
!
! TIMEGCM variables needed in output module:
!
  integer :: mtime(3) ! TIMEGCM model time (day,hour,min)
  real(r8),allocatable,dimension(:,:) :: & ! (nlon,nlat)
    tlbc,  & ! TN lower boundary
    ulbc,  & ! UN lower boundary
    vlbc     ! VN lower boundary
  real(r8) :: lbc ! bottom interface level (scalar)
!
! WACCM variables needed on output file:
!
  integer :: &
    date    ,& ! current date (YYYYMMDD)
    datesec ,& ! current seconds of current date
    nsteph  ,& ! current timestep
    ndcur   ,& ! current day (from base day)
    nscur      ! current seconds of current day

  real(r8),parameter :: eps = 1.e-6_r8

  private eps

  contains
!-----------------------------------------------------------------------
  subroutine get_geogrid(ncfile,model)
!
! Open netcdf input file, and get dimensions and coordinate variables
! of 3d global geographic grid, and time. This is called once per run, 
! and works for either TGCM or WACCM files.
!
! Args:
    character(len=*),intent(in) :: &
      ncfile,     & ! file path
      model         ! model name
!
! Local:
    integer :: istat,idunlim,id
    integer :: idv_lon,idv_lat,idv_lev,idv_ilev,idv_time
    character(len=1024) :: msg ! message for error handler
    character(len=NF90_MAX_NAME) :: varname
!
! Open netcdf history file.
!
    istat = nf90_open(ncfile,NF90_NOWRITE,ncid)
    if (istat /= NF90_NOERR) then
      write(msg,"('Error opening file ',a)") trim(ncfile)
      call handle_ncerr(istat,trim(msg),1)
    else
      write(iulog,"(/,'Opened ',a,' file ',a)") model,trim(ncfile)
!     write(iulog,"(/,'Opened file ',a)") trim(ncfile)
    endif
!
! Get number of times on the file (length of unlimited variable):
!
    istat = nf90_inq_dimid(ncid,'time',idunlim)
    istat = nf90_inquire_dimension(ncid,idunlim,varname,ntime)
!
! Get grid dimensions from the file:
!
    istat = nf90_inq_dimid(ncid,'lon',id)
    istat = nf90_inquire_dimension(ncid,id,varname,nlon)
    istat = nf90_inq_dimid(ncid,'lat',id)
    istat = nf90_inquire_dimension(ncid,id,varname,nlat)
    istat = nf90_inq_dimid(ncid,'lev',id)
    istat = nf90_inquire_dimension(ncid,id,varname,nlev)
    istat = nf90_inq_dimid(ncid,'ilev',id)
    istat = nf90_inquire_dimension(ncid,id,varname,nilev)
!
! Allocate coordinate variables:
!
    allocate(glon(nlon))
    allocate(glat(nlat))
    allocate(zlev(nlev))
    allocate(zilev(nilev))
    allocate(time(ntime))
!
! Read coordinate vars:
!
    istat = nf90_inq_varid(ncid,'lon',idv_lon)
    istat = nf90_get_var(ncid,idv_lon,glon,(/1/),(/nlon/))
    istat = nf90_inq_varid(ncid,'lat',idv_lat)
    istat = nf90_get_var(ncid,idv_lat,glat,(/1/),(/nlat/))
    istat = nf90_inq_varid(ncid,'lev',idv_lev)
    istat = nf90_get_var(ncid,idv_lev,zlev,(/1/),(/nlev/))
    istat = nf90_inq_varid(ncid,'ilev',idv_ilev)
    istat = nf90_get_var(ncid,idv_ilev,zilev,(/1/),(/nilev/))

    istat = nf90_inq_varid(ncid,'time',idv_time)
    istat = nf90_get_var(ncid,idv_time,time,(/1/),(/ntime/))
!
! Model-independent sizes/dimensions:
    nlonp1 = nlon+1
    nlonp2 = nlon+2
    nlatp1 = nlat+1

!   write(iulog,"('Read ',a,' input file: nlon=' ,i4,' glon=' ,/,(6f13.4))") model,nlon,glon
!   write(iulog,"('Read ',a,' input file: nlat=' ,i4,' glat=' ,/,(6f13.4))") model,nlat,glat
!   write(iulog,"('Read ',a,' input file: nlev=' ,i4,' zlev=' ,/,(6e13.4))") model,nlev,zlev
!   write(iulog,"('Read ',a,' input file: nilev=',i4,' zilev=',/,(6e13.4))") model,nilev,zilev
!   write(iulog,"('Read ',a,' input file: ntime=',i4,' time =',/,(6e13.4))") model,ntime,time

    write(iulog,"('get_geogrid for ',a,' ntime=',i4,' nlon=',i4,' nlat=',i4,' nlev=',i4,' nilev=',i4)") &
      trim(model),ntime,nlon,nlat,nlev,nilev

  end subroutine get_geogrid
!-----------------------------------------------------------------------
  subroutine read_waccm(ncfile,itime)
!
! Read WACCM history file inputs at time index itime (itime=1,ntime).
! File is open and in data mode under ncid when this routine is called.
!
! Note dimension and coordinate variables have been read by subroutine
! get_geogrid (at itime==1), and are stored in geogrid.F90 
! (either model, and in native model format).  Coord vars have been 
! allocated and defined by get_geogrid.
!
! Args:
    character(len=*),intent(in) :: ncfile ! netcdf file id
    integer,intent(in) :: itime
!
! Local:
    integer :: i,j,k,n,id,istat,start(4),count(4)
    integer :: idate,iyear,idatesec,month,day
    integer :: ivar1(1),ctpoten3(1,1,1)
    real(r8),allocatable,save :: f3d(:,:,:) ! global field (lev,lat,lon) 
    real(r8) :: real8,fmin,fmax
    logical,save :: first = .true.
!
! External:
    integer,external :: julian_day ! function in util.F90
!
! Input field names as they appear in the input  file:
!
    fnames = (/'T       ','U       ','V       ','OMEGA   ','Z3      ',&
               'MBARV   ','SIGMAPED','SIGMAHAL'/)

!   write(iulog,"('read_waccm: itime=',i4,' ncfile=',a)") itime,trim(ncfile)
!
! Read waccm date (YYYYMMDD), and convert to year and julian day of year,
! defining module data year and iday:
!
    istat = nf90_inq_varid(ncid,'date',id)
    istat = nf90_get_var(ncid,id,ivar1,(/itime/),(/1/))
    idate = ivar1(1)                        ! YYYYMMDAY
    date  = idate                           ! save to echo to output
    iyear = idate/10000                     ! 4-digit integer year
    year  = real(iyear)                     ! real year (module data)
    month = (idate-iyear*10000)/100         ! 2-digit month
    day   = ((idate-iyear*10000)-month*100) ! 2-digit day of month
    iday  = julian_day(month,day)           ! julian day (module data)
!
! Get ut (hours) from waccm datesec (seconds in current day):
!
    istat = nf90_inq_varid(ncid,'datesec',id)
    istat = nf90_get_var(ncid,id,ivar1,(/itime/),(/1/))
    idatesec = ivar1(1) ! will be echoed to output
    datesec  = idatesec
    ut = real(idatesec)/3600. ! module data

    mtime(1) = iday
    mtime(2) = idatesec/3600
    mtime(3) = (idatesec-mtime(2)*3600)/60
!
! Get variables from waccm file that need to be echoed to output file:
! nsteph, ndcur, nscur
!
    istat = nf90_inq_varid(ncid,'nsteph',id)
    istat = nf90_get_var(ncid,id,ivar1,(/itime/),(/1/))
    nsteph = ivar1(1)

    istat = nf90_inq_varid(ncid,'ndcur',id)
    istat = nf90_get_var(ncid,id,ivar1,(/itime/),(/1/))
    ndcur = ivar1(1)

    istat = nf90_inq_varid(ncid,'nscur',id)
    istat = nf90_get_var(ncid,id,ivar1,(/itime/),(/1/))
    nscur = ivar1(1)

!   write(iulog,"('read_waccm: iday=',i4,' idatesec=',i8,' mtime=',3i4,' nsteph=',i4,' ndcur=',i4,' nscur=',i8)") &
!     iday,idatesec,mtime,nsteph,ndcur,nscur
!
! If ctpoten was not provided by user via namelist, 
! then read from input file.
! (CTPOTEN from Joe's waccm-x file is redundant in lat,lon)
!
    if (ctpoten_user == fillvalue) then ! read from file
      istat = nf90_inq_varid(ncid,'CTPOTEN',id)
      istat = nf90_get_var(ncid,id,ctpoten3,(/itime,1,1/),(/1,1,1/))
      ctpoten = ctpoten3(1,1,1)
!     write(iulog,"('read_waccm: ctpoten as read from waccm file = ',f10.2)") ctpoten
    else                                ! use from namelist
      ctpoten = ctpoten_user
!     write(iulog,"('read_waccm: ctpoten as provided by user = ',f10.2)") ctpoten
    endif

!   write(iulog,"('read_waccm: date (YYYYMMDD)=',i10,' year=',f6.0,' iday=',i4,' ut=',f8.2)") &
!     idate,year,iday,ut
!
! Allocate input arrays on subdomains, and a 3d global field for reading:
! (first call only)
!
    if (first) then
      call alloc_inputs
      allocate(f3d(nlon,nlat,nlev))
      f3d = finit
!
! Shift waccm longitudes 180 degrees from 0->360 to -180->+180
! (There is no periodic point, so last coord is 177.5 (assuming dlon=2.5))
! See also 180 deg shift of waccm input arrays below, and the shift
! back to native WACCM sequence before writing to output file (output.F90)
!
      call shift_lon(glon,nlon,'-180to180',.true.)
!
! Reverse order of vertical coordinate to match TIMEGCM, i.e., so k=1 is 
! the bottom boundary, and k=nlev is the top boundary:
!
      call reverse_vec(zlev,nlev)
      call reverse_vec(zilev,nilev)
      call reverse_vec(zmlev,nmlev)
!
! j-index to geographic poles (waccm):
      jspole = 1
      jnpole = nlat
!
! Set horizontal geographic grid in radians (for apex code):
!
      allocate(ylatg(nlat))   ! waccm grid includes poles
      allocate(ylong(nlon+1)) ! single periodic point
      real8 = real(nlat) ; dlatg = pi/real8
      real8 = real(nlon) ; dlong = 2._r8*pi/real8
      ylatg(1)    = -pi/2._r8+eps ! south pole
      ylatg(nlat) =  pi/2._r8-eps ! north pole
      do j=2,nlat-1
        real8 = real(j-1)
        ylatg(j) = -0.5_r8*(pi-dlatg)+real8*dlatg
      enddo
      do i=1,nlon+1
        real8 = real(i-1)
        ylong(i) = -pi+real8*dlong
      enddo

!     write(iulog,"('read_waccm: nlat=',i4,' ylatg(1:nlat)=',/,(6e12.4))") nlat,ylatg
!     write(iulog,"('read_waccm: nlon=',i4,' ylong(1:nlon+1)=',/,(6e12.4))") nlon,ylong

    endif ! first call
!
! Read global input fields at current time:
!
    do n=1,nf
      istat = nf90_inq_varid(ncid,fnames(n),id)
      if (istat /= NF90_NOERR) then
        write(iulog,"('>>> Error getting var id for field ',a)") &
          fnames(n)
        cycle
      endif
!
! Read variable into 3d global array:
!
      start(1:3) = 1
      start(4) = itime
      count(1) = nlon
      count(2) = nlat
      count(3) = nlev
      count(4) = 1

      istat = nf90_get_var(ncid,id,f3d,start=start,count=count)
      if (istat /= NF90_NOERR) then
        write(iulog,"('>>> Error reading var ')") fnames(n)
        call handle_ncerr(istat,'Error from nf90_get_var',1)
      endif
!     write(iulog,"('Read field ',a,' from file: min,max=',2e12.4)") fnames(n),&
!       minval(f3d),maxval(f3d)
!
! Reverse order of vertical dimension (as in zlev above):
      do j=1,nlat
        do i=1,nlon
          call reverse_vec(f3d(i,j,:),nlev)
        enddo
!       write(iulog,"('Field ',a,' after reverse_vec: j=',i4,' min,max=',2e12.4)") &
!         fnames(n),j,minval(f3d(:,j,:)),maxval(f3d(:,j,:))
!
! Shift fields in longitude dimension by 180 degrees to match
! coord var shift from 0->360 to -180->+180 (see glon above)
! (this shifts from waccm-format to edynamo-format)
!
        do k=1,nlev ! nlev or nilev?
          call shift_lon(f3d(:,j,k),nlon,'-180to180',.false.)
        enddo
!       write(iulog,"('Field ',a,' after shift_lon: j=',i4,' min,max=',2e12.4)") &
!         fnames(n),j,minval(f3d(:,j,:)),maxval(f3d(:,j,:))
      enddo
!
! Convert some units:
      if (trim(fnames(n))=='U'.or.trim(fnames(n))=='V') then
        f3d = f3d*100.    ! convert winds from m/s to cm/s
        if (first) &
          write(iulog,"('Read WACCM field ',a,' 3d global min,max=',2e12.4,' (cm/s)')") &
            fnames(n),minval(f3d),maxval(f3d)
      elseif (trim(fnames(n))=='Z3') then
        f3d = f3d*100.    ! convert geopotential from m to cm
        if (first) &
          write(iulog,"('Read WACCM field ',a,' 3d global min,max=',2e12.4,' (cm)')") &
            fnames(n),minval(f3d),maxval(f3d)

! Need to check on OMEGA here (OMEGA:units = "Pa/s")

      else
        if (first) &
          write(iulog,"('Read WACCM field ',a,' 3d global min,max=',2e12.4)") &
          fnames(n),minval(f3d),maxval(f3d)
      endif
!
! Define fields on processor subdomain:
      do j=lat0,lat1
        do i=lon0,lon1
          if (fnames(n)=='T       ')    tn(:,i,j) = f3d(i,j,:) 
          if (fnames(n)=='U       ')    un(:,i,j) = f3d(i,j,:)
          if (fnames(n)=='V       ')    vn(:,i,j) = f3d(i,j,:)
          if (fnames(n)=='OMEGA   ') omega(:,i,j) = f3d(i,j,:)
          if (fnames(n)=='Z3      ')  zpot(:,i,j) = f3d(i,j,:)
          if (fnames(n)=='MBARV   ')  barm(:,i,j) = f3d(i,j,:)
          if (fnames(n)=='SIGMAPED')   ped(:,i,j) = f3d(i,j,:)
          if (fnames(n)=='SIGMAHAL')  hall(:,i,j) = f3d(i,j,:)
        enddo
      enddo
!
! Set kbotdyn (params.F90):
!
      if (fnames(n)=='Z3      '.and.first) then
        kbotdyn = 0
        do k=1,nlev
          fmin = minval(f3d(:,:,k)*1.e-5) ! minimum km
          fmax = maxval(f3d(:,:,k)*1.e-5) ! maximum km
          if (fmin < zbotdyn.and.fmax > zbotdyn) then
            kbotdyn = k
            exit
          endif
          if (fmin > zbotdyn.and.fmax > zbotdyn.and.k > 1) then
            kbotdyn = k-1
            exit
          endif
        enddo ! k=1,nlev
        if (kbotdyn > 0) then
          write(iulog,"('Set kbotdyn=',i4,' for WACCM: zlev(kbotdyn)=',e12.4,' zbotdyn=',f10.2,' Z min,max at kbotdyn=',2f10.2)") &
            kbotdyn,zlev(kbotdyn),zbotdyn,fmin,fmax
        else
          call shutdown('Could not set kbotdyn for WACCM file')
        endif
      endif ! Z3 and first call
    enddo ! n=1,nf
!
! Report min,max to stdout if first timestep:
!
    if (first) then
!
! Report to stdout:
      write(iulog,"(/,'WACCM file: nlev=',i4,' lev0,1=',2i4,' nlon=',i4,' lon0,1=',2i4,' nlat=',i4,' lat0,1=',2i4)") &
        nlev,lev0,lev1,nlon,lon0,lon1,nlat,lat0,lat1

      write(iulog,"('read_waccm: TN    subdomain min,max=',2e12.4)") minval(tn),   maxval(tn)
      write(iulog,"('read_waccm: UN    subdomain min,max=',2e12.4)") minval(un),   maxval(un)
      write(iulog,"('read_waccm: VN    subdomain min,max=',2e12.4)") minval(vn),   maxval(vn)
      write(iulog,"('read_waccm: OMEGA subdomain min,max=',2e12.4)") minval(omega),maxval(omega)
      write(iulog,"('read_waccm: ZPOT  subdomain min,max=',2e12.4)") minval(zpot), maxval(zpot)
      write(iulog,"('read_waccm: BARM  subdomain min,max=',2e12.4)") minval(barm), maxval(barm)
      write(iulog,"('read_waccm: PED   subdomain min,max=',2e12.4)") minval(ped),  maxval(ped)
      write(iulog,"('read_waccm: HALL  subdomain min,max=',2e12.4)") minval(hall), maxval(hall)
    endif

    first = .false.
  end subroutine read_waccm
!-----------------------------------------------------------------------
  subroutine read_tgcm(ncfile,itime)
!
! Args:
    character(len=*),intent(in) :: ncfile
    integer,intent(in) :: itime
!
! Local:
    integer :: i,ii,j,k,istat,id,n,start(4),count(4),lonbeg,lonend
    integer :: iyear,iyear1(1),iday1(1)
    real(r8),allocatable,save :: f3d(:,:,:) ! global field (lev,lat,lon) 
    real(r8) :: ut1(1),ctpoten1(1)
    real(r8) :: real8,fmin,fmax
    logical,save :: first = .true.

    fnames = (/'TN_DYN  ','UN_DYN  ','VN_DYN  ','W_DYN   ','Z_DYN   ',&
               'BARM_DYN','PED_DYN ','HAL_DYN '/)

    istat = nf90_inq_varid(ncid,'mtime',id)
    istat = nf90_get_var(ncid,id,mtime,(/1,itime/),(/3,1/))

    istat = nf90_inq_varid(ncid,'day',id)
    istat = nf90_get_var(ncid,id,iday1,(/itime/),(/1/))
    iday = iday1(1)

    istat = nf90_inq_varid(ncid,'ut',id)
    istat = nf90_get_var(ncid,id,ut1,(/itime/),(/1/))
    ut = ut1(1)

    istat = nf90_inq_varid(ncid,'year',id)
    istat = nf90_get_var(ncid,id,iyear1,(/itime/),(/1/))
    iyear = iyear1(1)
    year = real(iyear)
!
! If ctpoten was not provided by user via namelist, 
! then read from input file:
!
    if (ctpoten_user == fillvalue) then ! read from file
      istat = nf90_inq_varid(ncid,'ctpoten',id)
      istat = nf90_get_var(ncid,id,ctpoten1,(/itime/),(/1/))
      ctpoten = ctpoten1(1)
    else                                ! use from namelist
      ctpoten = ctpoten_user
    endif

!   write(iulog,"('read_tgcm: itime=',i4,' mtime=',3i4,' ut=',f8.3,' iday=',i4,' year=',f7.0)") &
!     itime,mtime,ut,iday,year
!
! Allocate a 3d global field for reading:
    if (first) then
      call alloc_inputs
      allocate(tlbc(nlon,nlat))     ! TN lower boundary
      allocate(ulbc(nlon,nlat))     ! UN lower boundary
      allocate(vlbc(nlon,nlat))     ! VN lower boundary
      allocate(f3d(nlon,nlat,nlev)) ! 3d global read array
      f3d = finit
!
! Set horizontal geographic grid in radians (for apex code):
      allocate(ylatg(0:nlat+1)) ! include poles for TIMEGCM
      allocate(ylong(nlon+1))   ! single periodic point
      real8 = real(nlat) ; dlatg = pi/real8
      real8 = real(nlon) ; dlong = 2._r8*pi/real8
      do j=1,nlat
        real8 = real(j-1)
        ylatg(j) = -0.5_r8*(pi-dlatg)+real8*dlatg
      enddo
      ylatg(0)      = -pi/2._r8+eps
      ylatg(nlat+1) =  pi/2._r8-eps
      do i=1,nlon+1
        real8 = real(i-1)
        ylong(i) = -pi+real8*dlong
      enddo
!
! j-index to geographic poles (timegcm):
      jspole = 0
      jnpole = nlat+1
 
!     write(iulog,"('read_tgcm: nlat=',i4,' ylatg(0:nlat+1)=',/,(6e12.4))") nlat,ylatg
!     write(iulog,"('read_tgcm: nlon=',i4,' ylong(1:nlon+1)=',/,(6e12.4))") nlon,ylong

    endif ! first call
!
! Read global input fields at current time:
!
    do n=1,nf
      istat = nf90_inq_varid(ncid,fnames(n),id)
      if (istat /= NF90_NOERR) then
        write(iulog,"('>>> Error getting var id for field ',a)") &
          fnames(n)
        cycle
      endif
!
! Read variable into 3d global array:
!
      start(1:3) = 1
      start(4) = itime

      count(1) = nlon
      count(2) = nlat
      count(3) = nlev
      count(4) = 1

      istat = nf90_get_var(ncid,id,f3d,start=start,count=count)
      if (istat /= NF90_NOERR) then
        write(iulog,"('>>> Error reading var ')") fnames(n)
        call handle_ncerr(istat,'Error from nf90_get_var',1)
      endif

      if (first) write(iulog,"('Read TIMEGCM field ',a,' 3d global min,max=',2e12.4)") &
          fnames(n),minval(f3d),maxval(f3d)
!
! Transfer to geographic subdomains:
!
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlon+4) lonend = lon1-2
      
      do j=lat0,lat1
        do i=lonbeg,lonend
          ii = i-2                                             ! body of longitude domain
          if (i <= 2)      ii = nlon-2+i  ! i=1,2 -> ii=71,72  ! west end periodic points
          if (i >= nlon+3) ii = i-nlon-2  ! i=75,76 -> ii=1,2  ! east end periodic points

! Subscript #3 of the array F3D has value 50 which is greater than the upper bound of 49
          if (fnames(n)=='TN_DYN  ')    tn(:,i,j) = f3d(ii,j,:) 
          if (fnames(n)=='UN_DYN  ')    un(:,i,j) = f3d(ii,j,:)
          if (fnames(n)=='VN_DYN  ')    vn(:,i,j) = f3d(ii,j,:)
          if (fnames(n)=='W_DYN   ') omega(:,i,j) = f3d(ii,j,:)
          if (fnames(n)=='Z_DYN   ')  zpot(:,i,j) = f3d(ii,j,:)
          if (fnames(n)=='BARM_DYN')  barm(:,i,j) = f3d(ii,j,:)
          if (fnames(n)=='PED_DYN ')   ped(:,i,j) = f3d(ii,j,:)
          if (fnames(n)=='HAL_DYN ')  hall(:,i,j) = f3d(ii,j,:)
        enddo
      enddo ! j=lat0,lat1
!
! Set kbotdyn (params.F90):
!
      if (fnames(n)=='Z_DYN   '.and.first) then
        kbotdyn = 0
        do k=1,nlev
          fmin = minval(f3d(:,:,k)*1.e-5) ! minimum km
          fmax = maxval(f3d(:,:,k)*1.e-5) ! maximum km
          if (fmin < zbotdyn.and.fmax > zbotdyn) then
            kbotdyn = k
            exit
          endif
          if (fmin > zbotdyn.and.fmax > zbotdyn.and.k > 1) then
            kbotdyn = k-1
            exit
          endif
        enddo ! k=1,nlev
        if (kbotdyn > 0) then
          write(iulog,"('Set kbotdyn=',i4,' for TIMEGCM: zlev(kbotdyn)=',e12.4,' zbotdyn=',f10.2,' Z min,max at kbotdyn=',2f10.2)") &
            kbotdyn,zlev(kbotdyn),zbotdyn,fmin,fmax
        else
          call shutdown('Could not set kbotdyn for TIMEGCM file')
        endif
      endif ! Z3 and first call
    enddo ! n=1,nf
!
! Read lower boundaries of t,u,v:
!
    istat = nf90_inq_varid(ncid,'LBC',id)
    if (istat /= NF90_NOERR) write(iulog,"('>>> Could not find LBC on input file.')")
    istat = nf90_get_var(ncid,id,lbc)

    istat = nf90_inq_varid(ncid,'TLBC',id)
    if (istat /= NF90_NOERR) write(iulog,"('>>> Could not find TLBC on input file.')")
    istat = nf90_get_var(ncid,id,tlbc,start=(/1,1,itime/),count=(/nlon,nlat,1/))

    istat = nf90_inq_varid(ncid,'ULBC',id)
    if (istat /= NF90_NOERR) write(iulog,"('>>> Could not find ULBC on input file.')")
    istat = nf90_get_var(ncid,id,ulbc,start=(/1,1,itime/),count=(/nlon,nlat,1/))

    istat = nf90_inq_varid(ncid,'VLBC',id)
    if (istat /= NF90_NOERR) write(iulog,"('>>> Could not find VLBC on input file.')")
    istat = nf90_get_var(ncid,id,vlbc,start=(/1,1,itime/),count=(/nlon,nlat,1/))
!
! Report min,max to stdout if first timestep:
!
    if (first) then
      write(iulog,"(/,'read_tgcm (first call): nlev=',i4,' lev0,1=',2i4,' nlon=',i4,' lon0,1=',2i4,' nlat=',i4,' lat0,1=',2i4)") &
        nlev,lev0,lev1,nlon,lon0,lon1,nlat,lat0,lat1

      write(iulog,"('read_tgcm: TN    subdomain min,max=',2e12.4)") minval(tn),   maxval(tn)
      write(iulog,"('read_tgcm: UN    subdomain min,max=',2e12.4)") minval(un),   maxval(un)
      write(iulog,"('read_tgcm: VN    subdomain min,max=',2e12.4)") minval(vn),   maxval(vn)
      write(iulog,"('read_tgcm: OMEGA subdomain min,max=',2e12.4)") minval(omega),maxval(omega)
      write(iulog,"('read_tgcm: ZPOT  subdomain min,max=',2e12.4)") minval(zpot), maxval(zpot)
      write(iulog,"('read_tgcm: BARM  subdomain min,max=',2e12.4)") minval(barm), maxval(barm)
      write(iulog,"('read_tgcm: PED   subdomain min,max=',2e12.4)") minval(ped),  maxval(ped)
      write(iulog,"('read_tgcm: HALL  subdomain min,max=',2e12.4)") minval(hall), maxval(hall)

!     write(iulog,"('read_tgcm: TLBC  global min,max=',2e12.4)") minval(tlbc), maxval(tlbc)
!     write(iulog,"('read_tgcm: ULBC  global min,max=',2e12.4)") minval(ulbc), maxval(ulbc)
!     write(iulog,"('read_tgcm: VLBC  global min,max=',2e12.4)") minval(vlbc), maxval(vlbc)
    endif

    first = .false.
  end subroutine read_tgcm
!-----------------------------------------------------------------------
  subroutine alloc_inputs
!
! Allocate input fields on geographic subdomain, and init:
!
!   write(iulog,"('allocate input fields on subdomains: lev0,1=',2i4,' lon0,1=',2i4,' lat0,1=',2i4)") &
!     lev0,lev1,lon0,lon1,lat0,lat1

    allocate(tn   (lev0:lev1,lon0:lon1,lat0:lat1))
    allocate(un   (lev0:lev1,lon0:lon1,lat0:lat1))
    allocate(vn   (lev0:lev1,lon0:lon1,lat0:lat1))
    allocate(omega(lev0:lev1,lon0:lon1,lat0:lat1))
    allocate(zpot (lev0:lev1,lon0:lon1,lat0:lat1))
    allocate(barm (lev0:lev1,lon0:lon1,lat0:lat1))
    allocate(ped  (lev0:lev1,lon0:lon1,lat0:lat1))
    allocate(hall (lev0:lev1,lon0:lon1,lat0:lat1))
    tn   =finit
    un   =finit
    vn   =finit 
    omega=finit 
    zpot =finit
    barm =finit
    ped  =finit 
    hall =finit
  end subroutine alloc_inputs
!-----------------------------------------------------------------------
end module read_ncfile