#include "dims.h"
!
      module ncep_module
!
! Handle NCEP lower boundaries.
!
      use cons_module,only: grav,gask,atm_amu
      use input_module,only: ncep,ncep_ncfile,tempdir
      use netcdf_module,only: handle_ncerr,nc_open,nc_close,
     |  nc_get_var_real
      implicit none
#include "params.h"
#include "ncep.h"
#include "netcdf.inc"
!
      integer,parameter :: nw=12, mxdays=30*365
!
      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(zimxp,zjmx,npr),tncep(zimxp,zjmx) are 
!   in shared common in ncep.h.
! 
!   
! 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
      integer,save :: ncid,ndays,iyd0,iyd1,iydp=999999
      integer,save :: iyds(mxdays) ! yyddd's from data file
      real,save :: glon(zimxp)
      real :: z4(zimxp,zjmx,npr,4) ! npr==3 (30, 10, and 5 mb)
      real :: press(npr)=(/30.,10.,5./), sdif(npr)
! To save z and t on secondary histories (redundant in pressure): 
      real :: fmin,fmax,zncepik(zimxp,zkmxp),tncepik(zimxp,zkmxp),
     |  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,zimxp
          glon(i) = -190.+(i-1)*5.
        enddo
!
! Get data for 4-day interpolation:
        do i=1,4
          iyd4(i) = inciyd(iyd,i-2)
          ixyd = getixyd(iyds,ndays,iyd4(i))
          if (ixyd==0) then
            write(6,"('>>> ncep_bndry: need day ',i5,' but could ',
     |        'not find it on ncep data file.')") iyd4(i)
            write(6,"('ndays=',i5,' iyds on file = ',/,(8i8))") 
     |        ndays,iyds(1:ndays)
            stop 'ncep iyd'
          endif
          call getncepday(ncid,iyd4(i),ixyd,glon,z4(1,1,1,i))
        enddo
      else ! istep > 1
!
! Advance the 4-day interpolates if necessary:
        if (iyd > iydp) then
          do i=1,4
            iyd4(i) = inciyd(iyd4(i),iyd-iydp) ! presumably iyd-iydp=1
            if (i <= 3) z4(:,:,:,i) = z4(:,:,:,i+1)
          enddo 
          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,zjmx
          call timeterp(z4(1,j,k,1),z4(1,j,k,2),z4(1,j,k,3),z4(1,j,k,4),
     +      zimxp,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,zjmx
        do i = 1,zimxp
          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),zimxp*zjmx,fmin,fmax)
        call fminmax(tncep,zimxp*zjmx,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
      endif
!
! Save output arrays zncep(zimxp,zjmx,npr),tncep(zimxp,zjmx) on
! secondary histories (redundant in pressure).
!
!     do j=1,zjmx
!       do k=1,zkmxp ! redundant in pressure for secondary histories
!         zncepik(:,k) = zncep(:,j,2) ! save at 10 mb level
!       enddo
!       call addfsech('Z_NCEP',' ',' ',zncepik,zimxp,zkmxp,zkmxp,j)
!
!       do k=1,zkmxp ! redundant in pressure for secondary histories
!         tncepik(:,k) = tncep(:,j) ! tncep is at 10 mb
!       enddo
!       call addfsech('T_NCEP',' ',' ',tncepik,zimxp,zkmxp,zkmxp,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:
      call mkdiskflnm(ncep_ncfile,dskfile)
      call getms(ncep_ncfile,dskfile,tempdir,' ')
      call nc_open(ncid,dskfile,'OLD','READ')
      istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim)
!
! Get and verify latitude dimension. 
! (read nlat, compare with zjmx, 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 /= zjmx) then
        write(6,"(/,'>>> getncepfile: bad nlat=',i3,
     |    ' -- should be zjmx=',i3)") nlat,zjmx
        stop '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
        stop '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
        stop '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')
        stop 'getncepfile'
      else
        if (ndays > mxdays) then
          write(6,"('>>> getncepfile: ndays=',i6,' mxdays=',i6)")
          write(6,"('    Please increase mxdays to at least ndays')")
          stop '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')
        stop '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')
        stop '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)
#include "grid.h"
!
! Args:
      integer,intent(in) :: ncid,iyd,ixyd
      real,intent(in) :: glonx(zimxp)
      real,intent(out) :: z(zimxp,zjmx,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(zjmx,npr), zs(zjmx,nw,npr), zc(zjmx,nw,npr)
      real :: fmin,fmax,pi,dtr
      real :: alfa,angle,wt1,wt,start(zimxp), stop(zimxp)
      character(len=80) :: char80
!
! Get zm(zjmx,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) = zjmx
      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(zjmx,nw,npr):
      istat = nf_inq_varid(ncid,"ZS",idv_zs)
      start4d(1:3) = 1
      start4d(4) = ixyd
      count4d(1) = zjmx
      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(zjmx,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,zjmx*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,:),zjmx*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,:),zjmx*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,zjmx
          do i = 1,zimxp
            z(i,j,k) = zm(j,k)
          enddo
        enddo
      enddo
      do k=1,npr
        do j = 1,zjmx
          angle = abs(glat(1) + (j-1)*dlat)*2.*alfa
          angle = amin1(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,zimxp
              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,zimxp
          start(i) = 25.*(z(i,1,k) + z(1+mod(i+35,zimx),1,k))
        enddo
        do j = 1,zjmx-1
          do i = 1,zimxp
            z(i,j,k) = 25.*(z(i,j,k) + z(i,j+1,k))
          enddo
        enddo
        do i = 1,zimxp
          stop(i) = 25.*(z(i,zjmx,k) + z(1+mod(i+35,zimx),zjmx,k))
        enddo
        do i = 1,zimxp
          z(i,zjmx,k) = stop(i)
        enddo
        do j = zjmx,2,-1
          do i = 1,zimxp
            z(i,j,k) = z(i,j,k) + z(i,j-1,k)
          enddo
        enddo
        do i = 1,zimxp
          z(i,1,k) = z(i,1,k) + start(i)
        enddo
      enddo
      end subroutine getncepday
      end module ncep_module