!
      module ncep_reanalysis_module
!
! Read T,Z lower boundary (10 mb) data from ncep reanalysis data file.
! The file is specified with namelist input parameter "ncep_reanalysis"
! (see input.F)
! See hao:/home/tgcm/ncep_reanalysis for procedure to acquire the data 
! (NCEP/NCAR dataset ds090.0), and make the 10 mb ncep_reanalysis netcdf 
! file. 
!
      use params_module,only: nlon,nlonp4,glon1,nlat,dlat,dlon,nlevp1
      use cons_module,only: grav,gask,atm_amu
      use nchist_module,only: handle_ncerr,nc_open
      use addfld_module,only: addfld
      implicit none
      integer :: ncid,ncep_nlev,ncep_ntime 
#include <netcdf.inc>
!
! Year, day of year, and hour from ncep file:
      integer,allocatable,dimension(:) :: ncep_year,ncep_day,ncep_hour
!
! Output of this module is T,Z at 10 mb interpolated to current time.
      real,dimension(nlonp4,nlat) :: z_ncep, t_ncep
      contains
!-----------------------------------------------------------------------
      subroutine ncep_reanalysis_bndry(ncfile,istep,mtime,iyear,iday,
     |  isecs)
!
! Called from advance. Return global z_ncep and t_ncep at current time.
!
      implicit none
!
! Args:
      character(len=*),intent(in) :: ncfile
      integer,intent(in) :: 
     |  istep,   ! model time step (1->nstep)
     |  iyear,   ! current model year (4-digit)
     |  iday,    ! current model calendar day (1-367)
     |  isecs,   ! current model time in seconds
     |  mtime(4) ! model time (day,hr,min,sec), for info only
!
! Local:
      integer :: i
      real,dimension(nlon,nlat) :: tncep0,tncep,tncep1,
     |                             zncep0,zncep,zncep1
      real :: hour
!
      hour = float(isecs)/3600.
!
! If first step, open and read coord vars of ncep data file:
      if (istep==1) then
        write(6,"('ncep_reanalysis_bndry: ncfile=',a)") ncfile
        call getncepfile(ncfile)
      endif
!
! Get ncep data at current model time:
      call getnceptime(iyear,iday,hour)

      end subroutine ncep_reanalysis_bndry
!-----------------------------------------------------------------------
      subroutine getncepfile(ncfile)
      implicit none
!
! Args:
      character(len=*),intent(in) :: ncfile
!
! Local:
      integer :: istat,ndims,idunlim,ngatts,nvars
      integer :: nlat_rd,nlon_rd
      integer :: id_lat,id_lon,id_time,id_lev
      integer :: idv_lev,idv_year,idv_day,idv_hour
      integer :: istart1(1),icount1(1)
      character(len=80) :: dskfile
      real,allocatable :: ncep_levs(:)
!
      dskfile = ' '
      call getfile(ncfile,dskfile)
      call nc_open(ncid,ncfile,'OLD','READ')
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
!     write(6,"('getncepfile: ndims=',i3,' nvars=',i3,' ngatts=',i3,
!    |  ' idunlim=',i3)") ndims,nvars,ngatts,idunlim
!
! Get and verify latitude dimension. 
! (read nlat_rd, compare with nlat, which is in params.h):
      istat = nf_inq_dimid(ncid,'lat',id_lat)
      istat = nf_inq_dimlen(ncid,id_lat,nlat_rd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting latitude dimension')
      if (nlat_rd /= nlat) then
        write(6,"(/,'>>> getncepfile: bad nlat_rd=',i3,
     |    ' -- should be nlat=',i3)") nlat_rd,nlat
        call shutdown('getncepfile')
      endif
!
! Get and verify longitude dimension. 
! (read nlon_rd, compare with nlon, which is in params.h):
      istat = nf_inq_dimid(ncid,'lon',id_lon)
      istat = nf_inq_dimlen(ncid,id_lon,nlon_rd)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting longitude dimension')
      if (nlon_rd /= nlon) then
        write(6,"(/,'>>> getncepfile: bad nlon_rd=',i3,
     |    ' -- should be nlon=',i3)") nlon_rd,nlon
        call shutdown('getncepfile')
      endif
!
! Get levels dimension (should be 1, i.e., only 10 mb level):
      istat = nf_inq_dimid(ncid,'lev',id_lev)
      istat = nf_inq_dimlen(ncid,id_lev,ncep_nlev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting levels dimension')
!
! Verify that first level on file = 10 mb.
      istat = nf_inq_varid(ncid,"lev",idv_lev) ! day var id
      allocate(ncep_levs(ncep_nlev),stat=istat)
      istart1(1) = 1
      icount1(1) = ncep_nlev
      istat = nf_get_vara_double(ncid,idv_lev,istart1,icount1,
     |  ncep_levs)
      write(6,"('getncepfile: ncep_nlev=',i3,' ncep_levs(1)=',f8.2)")
     |  ncep_nlev,ncep_levs(1)         
      if (ncep_levs(1) /= 10.) then
        write(6,"('>>> ncep_levs(1)=',f8.2,' -- expected 10. mb')")
     |    ncep_levs(1)
        call shutdown('getncepfile')
      endif
      deallocate(ncep_levs)
!
! Get number of times (value of unlimited dimension) on the file:
      istat = nf_inq_unlimdim(ncid,id_time) ! id of unlimited record var
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting unlimited dimension id')
      istat = nf_inq_dimlen(ncid,id_time,ncep_ntime) ! length of time var
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting length of day record variable')
        call shutdown('getncepfile')
      else
!       write(6,"('getncepfile: ncep_ntime=',i3)") ncep_ntime
      endif
!
! Read vars year,day_of_year,hour:
!     integer,allocatable,dimension(:) :: ncep_year,ncep_day,ncep_hour
!
      allocate(ncep_year(ncep_ntime),stat=istat)
      allocate(ncep_day (ncep_ntime),stat=istat)
      allocate(ncep_hour(ncep_ntime),stat=istat)
      istart1(1) = 1
      icount1(1) = ncep_ntime
!
! Year:
      istat = nf_inq_varid(ncid,"year",idv_year)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting id of var year')
      istat = nf_get_vara_int(ncid,idv_year,istart1,icount1,
     |  ncep_year)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error reading var year')
!     write(6,"('getncepfile: year=',/,(10i6))") ncep_year
!
! Day (of year):
      istat = nf_inq_varid(ncid,"day",idv_day)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting id of var day (day of year)')
      istat = nf_get_vara_int(ncid,idv_day,istart1,icount1,
     |  ncep_day)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error reading var day (day of year)')
!     write(6,"('getncepfile: day of year=',/,(10i6))") ncep_day
!
! Hour:
      istat = nf_inq_varid(ncid,"hour",idv_hour)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting id of var hour')
      istat = nf_get_vara_int(ncid,idv_hour,istart1,icount1,
     |  ncep_hour)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error reading var hour')
!     write(6,"('getncepfile: hour=',/,(10i6))") ncep_hour

      end subroutine getncepfile
!-----------------------------------------------------------------------
      subroutine getnceptime(iyear,iday,hour)
      implicit none
!
! Args:
      integer,intent(in) :: iyear,iday
      real,intent(in) :: hour
!     real,intent(in),dimension(nlon,nlat) :: tncep,zncep
!
! Local:
      integer :: i,itncep0,itncep1,ncephr1
!
!     write(6,"('getnceptime: searching for iyear=',i4,' iday=',i3,
!    |  ' hour=',f8.3)") iyear,iday,hour
      itncep0 = 0
      itncep1 = 0
!
! Check first time (if exact match, interp will not be necessary):
      if (ncep_year(1)==iyear.and.ncep_day(1)==iday.and.
     |  float(ncep_hour(1))==hour) then
        itncep0 = 1
      endif
!
! Search interior times:
      if (itncep0==0) then
        searchloop: do i=1,ncep_ntime-1
!
! Check for exact match at i:
          if (ncep_year(i)==iyear.and.ncep_day(i)==iday.and.
     |      float(ncep_hour(i))==hour) then
            itncep0 = i
            exit searchloop
          endif
!
! Check for exact match at i+1:
          if (ncep_year(i+1)==iyear.and.ncep_day(i+1)==iday.and.
     |      float(ncep_hour(i+1))==hour) then
            itncep0 = i+1 
            exit searchloop
          endif
!
! Check for bracket between i and i+1 (interp will be necessary):
! Note if ncep_hour(i) is 18, then ncep_hour(i+1) is 0, but need
! to test between 18 and 24 (not 18 and 0), see ncephr1:
          if (itncep0==0) then
            ncephr1 = ncep_hour(i+1)
            if (ncep_hour(i)==18) ncephr1 = 24
            if (ncep_year(i)==iyear.and.ncep_day(i)==iday.and.
     |        hour > float(ncep_hour(i)).and.
     |        hour < float(ncephr1)) then
              itncep0 = i
              itncep1 = i+1
              exit searchloop
            endif
          endif
        enddo searchloop
      endif
!
! Check for exact match at last time:
      if (itncep0==0) then
        if (ncep_year(ncep_ntime)==iyear.and.
     |      ncep_day(ncep_ntime)==iday.and.
     |      hour==ncep_hour(ncep_ntime)) then
          itncep0 = ncep_ntime 
        endif
      endif
!
! Time not found:
      if (itncep0==0) then
        write(6,"('>>> getnceptime: could not find iyear=',i4,
     |    ' iday=',i3,' hour=',f8.3)") iyear,iday,hour
        call shutdown('getnceptime')
      endif
!
!     write(6,"(/,'getnceptime: found iyear,iday,hour=',i5,i3,f8.3,
!    |  ' between itncep0,1=',2i4)") iyear,iday,hour,itncep0,itncep1
!     write(6,"('ncep_year,day,hour at itncep0=',3i4)")
!    |  ncep_year(itncep0),ncep_day(itncep0),ncep_hour(itncep0)
!     if (itncep1 > 0) 
!    |  write(6,"('ncep_year,day,hour at itncep1=',3i4)")
!    |    ncep_year(itncep1),ncep_day(itncep1),ncep_hour(itncep1)
!
! Read ncep data and do interpolation if necessary:
! (4th arg, ncep_hour, is ignored if itncep1==0)
      if (itncep1==0) then ! exact time match at itncep0
        call getncepdata(itncep0,itncep1,float(ncep_hour(itncep0)),
     |    float(itncep1),hour)
      else                 ! will interpolate between it0,it1
        ncephr1 = ncep_hour(itncep1)
        if (ncephr1==0) ncephr1 = 24
        call getncepdata(itncep0,itncep1,float(ncep_hour(itncep0)),
     |    float(ncephr1),hour)
      endif
      end subroutine getnceptime
!-----------------------------------------------------------------------
      subroutine getncepdata(it0,it1,hour0,hour1,hour)
!
! Read ncep data at time it0. If it1 > 0, then read ncep data
! at both it0 and it1, and interpolate to hour.
!
      implicit none
!
! Args:
      integer,intent(in) :: it0,it1
      real,intent(in) :: hour0,hour1,hour
!
! Local: 
      integer :: i,j,k
      real,dimension(nlon,nlat) :: tncep0,tncep1,zncep0,zncep1
      real,dimension(nlonp4,nlevp1) :: tncepik,zncepik
!
! External (util.F):
      real,external :: time_interp
!
! Read data at it0:
      call rdncepdata(it0,tncep0,zncep0)

!     write(6,"('getncepdata: it0,1=',2i4,' hour=',f7.3,
!    |  ' tncep0 min,max=',2e12.4,' zncep0 min,max=',2e12.4)") 
!    |  it0,it1,hour,minval(tncep0),maxval(tncep0),
!    |               minval(zncep0),maxval(zncep0)
!
! If it1 == 0, then no need to interpolate: 
! (exact time match was found at it0):
      if (it1==0) then
        t_ncep(3:nlon+2,:) = tncep0(:,:)
        z_ncep(3:nlon+2,:) = zncep0(:,:)
!
! If it1 > 0, then interpolate to needed hour:
      elseif (it1 > 0) then
        call rdncepdata(it1,tncep1,zncep1)
!
! Double check interpolates:
        if (hour0 > hour1) then
          write(6,"('>>> getncepdata: hour0 must be < hour1: hour0=',
     |      f8.3,' hour1=',f8.3)")
          call shutdown('getncepdata')
        elseif (hour < hour0.or.hour > hour1) then
          write(6,"('>>> getncepdata: hour must be between hour0,1:',
     |      ' hour0,1=',2f8.3,' hour=',f8.3)") hour0,hour1,hour
          call shutdown('getncepdata')
        endif
!
! Interpolate to hour:
! Util.F: real function time_interp(f0,f1,time0,time1,time)
!
        do j=1,nlat
          do i=3,nlon+2
            t_ncep(i,j) = time_interp(tncep0(i-2,j),tncep1(i-2,j),
     |        hour0,hour1,hour)
            z_ncep(i,j) = time_interp(zncep0(i-2,j),zncep1(i-2,j),
     |        hour0,hour1,hour)
          enddo
        enddo
      endif ! it1 > 0
!
! 1/26/07: This conversion is now included when the ncep_reanalysis file 
!          is written by the external code in hao:/home/tgcm/ncep_reanalysis.
! Convert Z from m to cm:
!     z_ncep = z_ncep*100.
!
! Periodic points:
! lons 1,2 <- nlon+1,nlon+2, and nlonp+3,nlon+4 <- 3,4
!
      t_ncep(1:2,:) = t_ncep(nlon+1:nlon+2,:)
      t_ncep(nlon+3:nlon+4,:) = t_ncep(3:4,:)
      z_ncep(1:2,:) = z_ncep(nlon+1:nlon+2,:)
      z_ncep(nlon+3:nlon+4,:) = z_ncep(3:4,:)
!
!     write(6,"('getncepdata: it0,1=',2i4,' hour0=',f7.3,' hour=',
!    |  f7.3,' hour1=',f7.3,' tncep0 mnmx=',2e12.4,' tncep1 mnmx=',
!    |  2e12.4)") it0,it1,hour0,hour,hour1,
!    |  minval(tncep0),maxval(tncep0),minval(tncep1),maxval(tncep1)
!     write(6,"('getncepdata: it0,1=',2i4,' hour0=',f7.3,' hour=',
!    |  f7.3,' hour1=',f7.3,' t_ncep mnmx=',2e12.4)")
!    |  it0,it1,hour0,hour,hour1,minval(t_ncep),maxval(t_ncep)
!
!     write(6,"('getncepdata: it0,1=',2i4,' hour0=',f7.3,' hour=',
!    |  f7.3,' hour1=',f7.3,' zncep0 mnmx=',2e12.4,' zncep1 mnmx=',
!    |  2e12.4)") it0,it1,hour0,hour,hour1,
!    |  minval(zncep0),maxval(zncep0),minval(zncep1),maxval(zncep1)
!     write(6,"('getncepdata: it0,1=',2i4,' hour0=',f7.3,' hour=',
!    |  f7.3,' hour1=',f7.3,' z_ncep mnmx=',2e12.4)")
!    |  it0,it1,hour0,hour,hour1,minval(z_ncep),maxval(z_ncep)
! 
!     call addfld('NCEP_T','NCEP reanalysis: TN at 10 mb','K',
!    |  t_ncep,'lon',1,nlonp4,'lat',1,nlat,0)
!     call addfld('NCEP_Z','NCEP reanalysis: Z at 10 mb','cm',
!    |  z_ncep,'lon',1,nlonp4,'lat',1,nlat,0)
!
! real,dimension(nlonp4,nlevp1) :: tncepik,zncepik
! Save on sec hist redundantly in vertical for tgcmproc_f90:
!
      do j=1,nlat
        do i=3,nlon+2
          tncepik(i,:) = t_ncep(i,j)
          zncepik(i,:) = z_ncep(i,j)
        enddo
        do k=1,nlevp1
          tncepik(1:2,k) = tncepik(nlon+1:nlon+2,j)
          tncepik(nlon+3:nlon+4,k) = tncepik(3:4,j)
          zncepik(1:2,k) = zncepik(nlon+1:nlon+2,j)
          zncepik(nlon+3:nlon+4,k) = zncepik(3:4,j)
        enddo
!
        call addfld('NCEP_T','NCEP reanalysis: TN at 10 mb','K',
     |    tncepik,'lon',1,nlonp4,'lev',1,nlevp1,j)
        call addfld('NCEP_Z','NCEP reanalysis: Z at 10 mb','cm',
     |    zncepik,'lon',1,nlonp4,'lev',1,nlevp1,j)
      enddo
!
      end subroutine getncepdata
!-----------------------------------------------------------------------
      subroutine rdncepdata(it,tncep,zncep)
!
! Read ncep T and Z at time index it.
!
      implicit none
!
! Args:
      integer,intent(in) :: it
      real,dimension(nlon,nlat),intent(out) :: tncep,zncep
!
! Local:
      integer :: istat,start_4d(4),count_4d(4),idv_t,idv_z
      character(len=120) :: char120
!
      start_4d(1:3) = 1
      start_4d(4) = it
      count_4d(1) = nlon
      count_4d(2) = nlat
      count_4d(3) = ncep_nlev
      count_4d(4) = 1
!
! Read T:
      istat = nf_inq_varid(ncid,"TMP_2_ISBL_10",idv_t)
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting var id of ncep TMP_2_ISBL_10')
        call shutdown('rdncepdata')
      endif
      istat = nf_get_vara_double(ncid,idv_t,start_4d,count_4d,tncep)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_get_vara_double',
     |    ' for field tncep: it=',i4)") it
        call handle_ncerr(istat,char120)
      else
!       write(6,"('rdncepdata: read tncep: it=',i4,' min,max=',2e12.4)") 
!    |    it,minval(tncep),maxval(tncep)
      endif
!
! Read Z:
      istat = nf_inq_varid(ncid,"HGT_2_ISBL_10",idv_z)
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting var id of ncep HGT_2_ISBL_10')
        call shutdown('rdncepdata')
      endif
      istat = nf_get_vara_double(ncid,idv_z,start_4d,count_4d,zncep)
      if (istat /= NF_NOERR) then
        write(char120,"('Error return from nf_get_vara_double',
     +    ' for field zncep')")
        call handle_ncerr(istat,char120)
      else
!       write(6,"('rdncepdata: read zncep: it=',i4,' min,max=',2e12.4)") 
!    |    it,minval(zncep),maxval(zncep)
      endif
      end subroutine rdncepdata
!-----------------------------------------------------------------------
      end module ncep_reanalysis_module