!
      module ncep_module
!
! Handle NCEP lower boundaries.
!
      use params_module,only: nlon,nlonp4,glon1,nlat,dlat,dlon,nlevp1
      use cons_module,only: grav,gask,atm_amu
      use input_module,only: ncep_ncfile
      use nchist_module,only: handle_ncerr,nc_open
      use addfld_module,only: addfld
      implicit none
#ifdef SUN
#include <netcdf3.inc>
#else
#include <netcdf.inc>
#endif
!
      integer,parameter :: nw=12, mxdays=30*365
      integer,parameter :: npr=3
      real :: zncep(nlonp4,nlat,npr),tncep(nlonp4,nlat)
!
      contains
!-----------------------------------------------------------------------
      subroutine ncep_bndry(istep,mtime,iyear,iday,isecs)
!
! Get ncep Z lower boundaries at 10 mb from netcdf file, and
!   calculate 10 mb T boundaries from Z.
!
! This routine is called at every timestep before the main latitude
!   scan. Data are interpolated to current model time. Z at 3 pressure
!   levels (30, 10, and 5 mb) is read from an external file. T is 
!   calculated using Z at all three pressures.
! Output arrays zncep(nlonp4,nlat,npr),tncep(nlonp4,nlat) are module data
!   
! Args:
      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,k,j,n,iyd,ixyd,ih,im,is,iyd4(4),ipos,idatsec,ier,
     |  ixyd0
      integer,save :: ncid,ndays,iyd0,iyd1,iydp=999999
      integer,save :: iyds(mxdays) ! yyddd's from data file
      real,save :: glon(nlonp4)
      real,save :: z4(nlonp4,nlat,npr,4) ! npr==3 (30, 10, and 5 mb) 
!                                       ! data for 4 days 
      real :: press(npr)=(/30.,10.,5./), sdif(npr)
! To save z and t on secondary histories (redundant in pressure): 
      real :: fmin,fmax,zncepik(nlonp4,nlevp1),tncepik(nlonp4,nlevp1),
     |  fmin1,fmax1
!
! External:
      integer,external :: inciyd
!
! iyd is current year-day yyyyddd (e.g., 2001080 for day 80 of 2001)
      iyd = iyear*1000+iday
      call isec2hms(isecs,ih,im,is)
!     write(6,"('ncep_bndry: istep=',i4,' iyear=',i4,' iday=',i4,
!    |  ' isecs=',i5,' iyd=',i8,' hr:min:sec=',i2,':',i2,':',i2 )") 
!    |  istep,iyear,iday,isecs,iyd,ih,im,is
!
! If first step, get ncep data file and read iyds(ndays).
! (ndays and iyds(ndays) are returned by getncepfile)
      if (istep==1) then
        call getncepfile(ncid,ndays,iyds)
!       write(6,"('ncep_bndry: ndays=',i5,' days (iyds)=',/,(8i8))") 
!    |    ndays,iyds(1:ndays) 
        iyd0 = iyds(1)     ! first day on the file
        iyd1 = iyds(ndays) ! last day on the file
        write(6,"(/72('-'))")
        write(6,"('NCEP_BNDRY:')")
        write(6,"('Lower boundary for heights (10mb) will be obtained ',
     +    ' at each time step from')")
        write(6,"('  NCEP data file ',a)") trim(ncep_ncfile)
        write(6,"('Interpolation starts at year-day yyyyddd ',i8,
     |    ' hr:min:sec ',i2,':',i2,':',i2)") iyd,ih,im,is
!
        do i=1,nlonp4
          glon(i) = glon1+(i-1)*dlon
        enddo
!
! Get data for 4-day interpolation (current day is day 2 of 4):
        ixyd0 = getixyd(iyds,ndays,iyd) ! index of current day in iyds(ndays)
        do i=1,4
          ixyd = ixyd0+i-2
          if (ixyd <= 0.or.ixyd > ndays) then
            write(6,"('>>> ncep_bndry step 1: bad index to iyds: ixyd0='
     |        ,i6,' i=',i2,' ixyd=',i6,' iyd=',i7,' iyd0=',i7,' iyd1=',
     |        i7)") ixyd0,i,ixyd,iyd,iyd0,iyd1
            call shutdown('ncep iyd')
          endif
          iyd4(i) = iyds(ixyd)
! subroutine getncepday(ncid,iyd,ixyd,glonx,z)
          call getncepday(ncid,iyd4(i),ixyd,glon,z4(1,1,1,i))
        enddo

        write(6,"('ncep_bndry step 1: iyd=',i8,' iyd4=',4i8)") iyd,iyd4

      else ! istep > 1
!
! Advance the 4-day interpolates if necessary:
        ixyd0 = getixyd(iyds,ndays,iyd) ! index of current day in iyds(ndays)
        if (iyd > iydp) then
          do i=1,4
            ixyd = ixyd0+i-2
            if (ixyd <= 0.or.ixyd > ndays) then
              write(6,"('>>> ncep_bndry: bad index to iyds: ixyd0=',
     |          i6,' i=',i2,' ixyd=',i6,' iyd=',i7,' iyd0=',i7,' iyd1=',
     |          i7)") ixyd0,i,ixyd,iyd,iyd0,iyd1
              call shutdown('ncep iyd')
            endif
            iyd4(i) = iyds(ixyd)
            if (i <= 3) z4(:,:,:,i) = z4(:,:,:,i+1)
          enddo 

          write(6,"('ncep_bndry advance day: iyd=',i8,' iyd4=',4i8)") 
     |      iyd,iyd4

          ixyd = getixyd(iyds,ndays,iyd4(4))
          call getncepday(ncid,iyd4(4),ixyd,glon,z4(1,1,1,4))
        endif
      endif ! istep==1
!
! Do interpolation to current time, defining zncep:
      idatsec = 0
      ipos = 2
      if (iyd==iyd0) ipos = 1
      if (iyd==iyd1-1) ipos = 3
      do k=1,npr
        do j=1,nlat
          call timeterp(z4(1,j,k,1),z4(1,j,k,2),z4(1,j,k,3),z4(1,j,k,4),
     +      nlonp4,ipos,idatsec,isecs,zncep(1,j,k),0,ier)
        enddo
      enddo
!
! Calculate T at lower boundary from:
!  T = MBAR*g/R * dZ/ds
!  and
!  dZ/ds = ((ln(p2/p3))**2*(z2-z1) + (ln(p1/p2))**2*
!          (z3-z1))/(-ln(p1/p2)*ln(p21/p32)*ln(p31/p12))
!  Calculate differences between the three pressure levels
!  in terms of the pressure coordinate, s.
!
      do n = 1,3
        sdif(n) = alog(press(n)/press(1+mod(n,3)))
      enddo
      do j = 1,nlat
        do i = 1,nlonp4
          tncep(i,j) = atm_amu*grav/gask*(sdif(2)**2*
     |                (zncep(i,j,2)-zncep(i,j,1)) + sdif(1)**2*
     |                (zncep(i,j,3)-zncep(i,j,2))) /
     |                (-sdif(1)*sdif(2)*sdif(3))
        enddo
      enddo
!
! Report to stdout when new ncep dataset is obtained for new day:
      if (istep==1 .or. iyd > iydp) then
!     if (istep==1.or.mod(istep,10)==0 .or. iyd > iydp) then
        call fminmax(zncep(:,:,2),nlonp4*nlat,fmin,fmax)
        call fminmax(tncep,nlonp4*nlat,fmin1,fmax1)
!       write(6,"('NCEP at 10mb: istep=',i5,' mtime=',4i4,/,
!    |    '  Z min,max=',2e12.4,' T min,max=',2e12.4)")
!    |    istep,mtime,fmin,fmax,fmin1,fmax1
!       write(6,"('NCEP Z min,max=',2e12.4,' T min,max=',2e12.4)")
!    |    fmin,fmax,fmin1,fmax1
      endif
!
! Save output arrays zncep(nlonp4,nlat,npr),tncep(nlonp4,nlat) on
! secondary histories (redundant in pressure).
!
!     do j=1,nlat
!       do k=1,nlevp1 ! redundant in pressure for secondary histories
!         zncepik(:,k) = zncep(:,j,2) ! save at 10 mb level
!       enddo
!       call addfld('Z_NCEP',' ',' ',zncepik,'lon',1,nlonp4,
!    |    'lev',1,nlevp1,j)
!       do k=1,nlevp1 ! redundant in pressure for secondary histories
!         tncepik(:,k) = tncep(:,j) ! tncep is at 10 mb
!       enddo
!       call addfld('T_NCEP',' ',' ',tncepik,'lon',1,nlonp4,
!    |    'lev',1,nlevp1,j)
!     enddo
!
      iydp = iyd
      end subroutine ncep_bndry
!-----------------------------------------------------------------------
      subroutine getncepfile(ncid,ndays,iyds)
!
! Acquire ncep data file, verify dimensions, and read days array (yyddd)
!   from the file. Return ndays and iyds(ndays). File is left open in
!   read mode.
! This routine is called once per run, when istep==1.
!
! Args:
      integer,intent(out) :: ncid,ndays,iyds(mxdays)
!
! Local:
      character(len=80) :: dskfile
      integer,parameter :: mxdims=4
      integer :: istat,ndims,iddims(mxdims),idunlim,ngatts,nvars
      integer :: id_lat,id_nw,id_lev,id_day,idv_day,nlat,nwave,nlev
      integer :: istart1(1),icount1(1)
!
! Get file from mss:
      dskfile = ' '
      call getfile(ncep_ncfile,dskfile)
      call nc_open(ncid,dskfile,'OLD','READ')
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
!
! Get and verify latitude dimension. 
! (read nlat, compare with nlat, which is in params.h):
      istat = nf_inq_dimid(ncid,'lat',id_lat)
      istat = nf_inq_dimlen(ncid,id_lat,nlat)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting latitude dimension')
      if (nlat /= nlat) then
        write(6,"(/,'>>> getncepfile: bad nlat=',i3,
     |    ' -- should be nlat=',i3)") nlat,nlat
        call shutdown('getncepfile')
      endif
!
! Get and verify nw dimension. 
! (read nwave, compare with nw, which is parameter in this module):
      istat = nf_inq_dimid(ncid,'nw',id_nw)
      istat = nf_inq_dimlen(ncid,id_nw,nwave)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting nw dimension')
      if (nwave /= nw) then
        write(6,"(/,'>>> getncepfile: bad nwave=',i3,
     |    ' -- should be nw=',i3)") nwave,nw
        call shutdown('getncepfile')
      endif
!
! Get and verify lev dimension.
! (read nlev, compare with npr, which is parameter in ncep.h):
      istat = nf_inq_dimid(ncid,'lev',id_lev)
      istat = nf_inq_dimlen(ncid,id_lev,nlev)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'getncepfile: error getting lev dimension')
      if (nlev /= npr) then
        write(6,"(/,'>>> getncepfile: bad nlev=',i3,
     |    ' -- should be npr=',i3)") nlev,npr
        call shutdown('getncepfile')
      endif
!
! Get number of days (value of unlimited dimension) on the file:
      istat = nf_inq_unlimdim(ncid,id_day) ! 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_day,ndays) ! length of time var
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'Error getting length of day record variable')
        call shutdown('getncepfile')
      else
        if (ndays > mxdays) then
          write(6,"('>>> getncepfile: ndays=',i6,' mxdays=',i6)")
          write(6,"('    Please increase mxdays to at least ndays')")
          call shutdown('getncepfile')
        endif
      endif
!
! Get days (yyddd) array:
      istat = nf_inq_varid(ncid,"day",idv_day) ! day var id
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'getncepfile: error getting id of variable day')
        call shutdown('getncepfile')
      endif
      istart1(1) = 1
      icount1(1) = ndays
      iyds(:) = -1 ! init
!
! This can apparently take several minutes on the Cray when ndays
! is large (like > 7500). Is much faster on the ibm (a few secs).
! This can be seen with interactive ncdump -c commands on the
! file on the different machines.
!
      istat = nf_get_vara_int(ncid,idv_day,istart1,icount1,
     |  iyds(1:ndays))
      if (istat /= NF_NOERR) then
        call handle_ncerr(istat,
     |    'getncepfile: error getting day array')
        call shutdown('getncepfile')
      endif
      end subroutine getncepfile
!-----------------------------------------------------------------------
      integer function getixyd(iyds,niyds,iyd)
!
! Args:
      integer,intent(in) :: niyds,iyd
      integer,intent(in) :: iyds(niyds)
!
! Local:
      integer :: i
!
      getixyd = 0 ! index to requested iyd in iyds(niyds)
      if (iyd < iyds(1)) then ! requested day is before 1st day
        write(6,"('>>> getixyd: requested year-day ',i6,' is less',
     |    ' than first available year-day',i6)") iyd,iyds(1)
        return
      elseif (iyd > iyds(niyds)) then ! requested day is after last day
        write(6,"('>>> getixyd: requested year-day ',i6,' is ',
     |    'greater than last available year-day',i6)") iyd,iyds(niyds)
        return
      else ! requested day is within range of first and last days
        iydloop: do i=1,niyds
          if (iyds(i)==iyd) then
            getixyd = i
            exit iydloop
          endif 
        enddo iydloop
      endif
      end function getixyd
!-----------------------------------------------------------------------
      subroutine getncepday(ncid,iyd,ixyd,glonx,z)
      use params_module,only: glat
!
! Args:
      integer,intent(in) :: ncid,iyd,ixyd
      real,intent(in) :: glonx(nlonp4)
      real,intent(out) :: z(nlonp4,nlat,npr) ! output ncep geopotential
!
! Local:
      integer :: i,k,j,iw,istat,idv_zm,idv_zs,idv_zc
      integer :: start3d(3),count3d(3),start4d(4),count4d(4)
      real :: zm(nlat,npr), zs(nlat,nw,npr), zc(nlat,nw,npr)
      real :: fmin,fmax,pi,dtr
      real :: alfa,angle,wt1,wt,start(nlonp4), stop(nlonp4)
      character(len=80) :: char80
!
! Get zm(nlat,npr):
      istat = nf_inq_varid(ncid,"ZM",idv_zm)
      start3d(1:2) = 1
      start3d(3) = ixyd  ! time index on file for requested iyd
      count3d(1) = nlat
      count3d(2) = npr
      count3d(3) = 1
      istat = nf_get_vara_double(ncid,idv_zm,start3d,count3d,zm)
      if (istat /= NF_NOERR) then
        write(char80,"('Error return from nf_get_vara_double',
     +    ' for var zm iyd=',i8,' ixyd=',i8)") iyd,ixyd
        call handle_ncerr(istat,char80)
      endif
!
! Get zs(nlat,nw,npr):
      istat = nf_inq_varid(ncid,"ZS",idv_zs)
      start4d(1:3) = 1
      start4d(4) = ixyd
      count4d(1) = nlat
      count4d(2) = nw
      count4d(3) = npr
      count4d(4) = 1
      istat = nf_get_vara_double(ncid,idv_zs,start4d,count4d,zs)
      if (istat /= NF_NOERR) then
        write(char80,"('Error return from nf_get_vara_double',
     +    ' for var zs iyd=',i8,' ixyd=',i8)") iyd,ixyd
        call handle_ncerr(istat,char80)
      endif
!
! Get zc(nlat,nw,npr):
      istat = nf_inq_varid(ncid,"ZC",idv_zc)
      istat = nf_get_vara_double(ncid,idv_zc,start4d,count4d,zc)
      if (istat /= NF_NOERR) then
        write(char80,"('Error return from nf_get_vara_double',
     +    ' for var zc iyd=',i8,' ixyd=',i8)") iyd,ixyd
        call handle_ncerr(istat,char80)
      endif
!
! Report to stdout:
!     call fminmax(zm,nlat*npr,fmin,fmax)
!     write(6,"(/,'getncepday: iyd=',i8,' zm(nlat,npr) min,max=',
!    |  2e12.4)") iyd,fmin,fmax
!     do i=1,nw
!       call fminmax(zs(:,i,:),nlat*npr,fmin,fmax)
!       write(6,"('i=',i3,'  nw=',i3,' zs(nlat,i,npr) min,max=',
!    |    2e12.4)") i,nw,fmin,fmax
!       call fminmax(zc(:,i,:),nlat*npr,fmin,fmax)
!       write(6,"(14x,'zc(nlat,i,npr) min,max=',2e12.4)") fmin,fmax
!     enddo
!
! Do expansion to the desired longitudes:
      pi = 4.*atan(1.)
      dtr = pi/180.
      alfa = 3.
      do k=1,npr
        do j = 1,nlat
          do i = 1,nlonp4
            z(i,j,k) = zm(j,k)
          enddo
        enddo
      enddo
      do k=1,npr
        do j = 1,nlat
          angle = abs(glat(1) + float(j-1)*dlat)*2.*alfa
          angle = min(angle,180.)
          wt1 = (1.-cos(angle*dtr))/2.
          do iw = 1,4
            if (iw <= 2) then
              wt = 1.
            else
              wt = wt1
            endif
            do i = 1,nlonp4
              z(i,j,k) = z(i,j,k) + wt*(zs(j,iw,k) * sin(float(iw)*
     +          glonx(i)*dtr)+ zc(j,iw,k) * cos(float(iw)*glonx(i)*dtr))
            enddo
          enddo
        enddo
      enddo
! 
! (1:2:1) latitudinal smoothing and x 100 for cm
!
      do k=1,npr
        do i = 1,nlonp4
          start(i) = 25.*(z(i,1,k) + z(1+mod(i+35,nlon),1,k))
        enddo
        do j = 1,nlat-1
          do i = 1,nlonp4
            z(i,j,k) = 25.*(z(i,j,k) + z(i,j+1,k))
          enddo
        enddo
        do i = 1,nlonp4
          stop(i) = 25.*(z(i,nlat,k) + z(1+mod(i+35,nlon),nlat,k))
        enddo
        do i = 1,nlonp4
          z(i,nlat,k) = stop(i)
        enddo
        do j = nlat,2,-1
          do i = 1,nlonp4
            z(i,j,k) = z(i,j,k) + z(i,j-1,k)
          enddo
        enddo
        do i = 1,nlonp4
          z(i,1,k) = z(i,1,k) + start(i)
        enddo
      enddo
      end subroutine getncepday
      end module ncep_module