      subroutine getvtau(vtaupak,vtaupal,nlon1,nlat,
     &                   irefyr,irefmn,iday,kount,ijpak,nuvt,gg)
C

C     * Dec 20, 2016 - J. Cole.    Modified code so it can read the data
C                                  as fully transient, repeated annual cycle
C                                  or repeating a particular month.
C
C     * April 17, 2015 - M.Lazare.     Cosmetic revision (same name kept)
C     *                                to have "keeptim" common block
C     *                                definition line in upper case,
C     *                                so that "myrssti" does not conflict
C     *                                with the CPP directive (compiler
C     *                                warnings generated).
C
c***********************************************************************
c Read monthly averaged Sato et al (1993) stratospheric volcanic aerosol
c optical depths from a file every mid month day.
c This is a multi year data set (currently for years 1850-2002) so this
c routine also depends on iyear, passed in via common /keeptim/.
c
c Larry Solheim Nov,2007
c***********************************************************************

      implicit real (a-h,o-z), integer (i-n) 

      !--- Reference year and month.
      integer irefyr, irefmn

      !--- pak is the interpolated value and pal is the target value
      real vtaupak(ijpak), vtaupal(ijpak)

      real gg(*)

      !--- local
      integer :: year_prev, year_next, mon_prev, mon_next

      !--- keeptime is required for iyear
      COMMON /KEEPTIM/ IYEAR,IMDH,MYRSSTI,ISAVDTS

      !--- mynode used to control i/o to stdout
      integer*4 mynode
      common /mpinfo/ mynode 

      !--- A list containing mid month days of the year
      integer, parameter, dimension(12) ::
     &  mm_doy = (/ 16,46,75,106,136,167,197,228,259,289,320,350 /)

      !--- A list containing mid month days for each month
      integer, parameter, dimension(12) ::
     &  mm_dom = (/ 16,14,16,15,16,15,16,16,15,16,15,16 /)
C----------------------------------------------------------------------

c.....Determine the year for which data is required

      if (irefyr .gt. 0) then
        !--- When irefyr > 0 read emissions from the data file for
        !--- a repeated annual cycle of year irefyr
        curr_year=irefyr
      else
        !--- Otherwise read time dependent emissions from the data
        !--- file using iyear to determine the current year
        curr_year=iyear
      endif

c.....check on iday
      if (iday.lt.1 .or. iday.gt.365) then
        if (mynode.eq.0) then
          write(6,*)'GETVTAU: IDAY is out of range. IDAY=',iday
        endif
        call xit("GETVTAU",-1)
      endif

c.....Determine previous and next year/month, relative to iday
      mon_prev=0
      if (irefmn .lt. 0) then
         if (iday.ge.1 .and. iday.lt.mm_doy(1)) then
            year_prev=curr_year-1
            year_next=curr_year
            mon_prev=12
            mon_next=1
         else if (iday.ge.mm_doy(12) .and. iday.le.365) then
            year_prev=curr_year
            year_next=curr_year+1
            mon_prev=12
            mon_next=1
         else
            do i=2,12
               if (iday.ge.mm_doy(i-1) .and. iday.lt.mm_doy(i)) then
                  year_prev=curr_year
                  year_next=curr_year
                  mon_prev=i-1
                  mon_next=i
                  exit
               endif
            enddo
         endif
      else
         !--- The month is not changing.
         mon_prev = irefmn
         mon_next = irefmn
      end if

      if (mon_prev.eq.0) then
        if (mynode.eq.0) then
          write(6,*)'GETVTAU: Problem setting year/month on IDAY=',iday
        endif
        call xit("GETVTAU",-2)
      endif

      if (kount.eq.0) then     !--- Starting up

        !--- Get field for previous mid-month day
        if (mynode.eq.0) then
          write(6,*)"GETVTAU: Reading volcanic optical depth values",
     &      " for year,month: ",year_prev,mon_prev," IDAY=",iday,
     &      " KOUNT=",kount
        endif
        ibuf2 = 1000000*year_prev+10000*mon_prev+100*mm_dom(mon_prev)
        rewind nuvt
        call getggbx(vtaupak,nc4to8("VTAU"),nuvt,
     &               nlon1,nlat,ibuf2,1,gg)

        !--- Get field for next mid-month day
        if (mynode.eq.0) then
          write(6,*)"GETVTAU: Reading volcanic optical depth values",
     &      " for year,month: ",year_next,mon_next," IDAY=",iday,
     &      " KOUNT=",kount
        endif
        ibuf2 = 1000000*year_next+10000*mon_next+100*mm_dom(mon_next)
        rewind nuvt
        call getggbx(vtaupal,nc4to8("VTAU"),nuvt,
     &               nlon1,nlat,ibuf2,1,gg)

        !--- Interpolate to current day
        day1=real(mm_doy(mon_prev)) 
        day2=real(iday)
        day3=real(mm_doy(mon_next))
        if(day2.lt.day1) day2=day2+365. 
        if(day3.lt.day2) day3=day3+365. 
        w1=(day2-day1)/(day3-day1)
        w2=(day3-day2)/(day3-day1)
        if (mynode.eq.0) then
          write(6,*)'GETVTAU: Interpolating between days ',
     &      mm_doy(mon_prev),' and ',mm_doy(mon_next),' to day ',iday,
     &      ' using weights ',w1,w2
        endif

        do i=1,ijpak
          vtaupak(i) = w1*vtaupal(i) + w2*vtaupak(i)
        enddo

      else     !--- Continue a run

        if (mynode.eq.0) then
          write(6,*)"GETVTAU: Reading volcanic optical depth values",
     &      " for year,month: ",year_next,mon_next," IDAY=",iday,
     &      " KOUNT=",kount
        endif
        !--- Get field for next mid-month day
        ibuf2 = 1000000*year_next+10000*mon_next+100*mm_dom(mon_next)
        rewind nuvt
        call getggbx(vtaupal,nc4to8("VTAU"),nuvt,
     &               nlon1,nlat,ibuf2,1,gg)
      endif

      return
      end
