SUBROUTINE onebox(thybp,ofil,orec,lwindo,x000,y000,tstep,ftmax,pstep,slack,pnmax,pqmin)
!
!  sink = grid of moisture source/sink in back/forward trajectory [output]
!
!  thybp = hybrid precip time series at the grid point x000,y000
!  lwindo = number of days over which parcels are launched
!  x000, y000 = indicies for grid point of launch site
!  tstep = number of time steps per reanalysis data interval
!          (negative for back trajectories)
!  ftmax = # reanalysis time steps in period
!  stday = start day (earliest day of lwindo for forward traj.
!                      latest day of lwindo for back traj)
!  pstep = # of precipitation steps in pentad 
!
   USE forcing_mod
   USE pointers_mod
!
   IMPLICIT NONE
!
   INTEGER       :: ofil       ! Output file number
   INTEGER       :: orec       ! Output record number
   INTEGER       :: lwindo     ! number of days over which parcels are launched
   INTEGER       :: dura       ! # reanal timesteps in launch window
   INTEGER       :: x000,y000  ! x & y grid box for parcel launches
   INTEGER       :: it1,it2    ! Start and end parcel launch window in reanal data
   INTEGER       :: ip1,ip2    ! Start and stop intervals for hourly precip data
   INTEGER       :: ilon,ilat  ! x & y grid box for source contribution
   INTEGER       :: ip         ! precip loop counter
   INTEGER       :: iz         ! sigma level loop counter
   INTEGER       :: islev      ! sigma level loop counter
   INTEGER       :: pact       ! number of active parcels
   INTEGER       :: pn         ! parcel loop counter
   INTEGER       :: pcnt       ! number of parcels launched
!                   
   INTEGER       :: i0,i1      ! Indices for instantaneous (sigma level) data
   INTEGER       :: l0,l1,l2   ! Indices for time averaged (flux) data
   INTEGER       :: p0,p1,p2   ! Indices for hybrid precipitation data
   INTEGER       :: ftmax      ! # reanalysis time steps in period
   INTEGER       :: pstep      ! # of precipitation steps in pentad
   INTEGER       :: pnmax      ! Total number of parcels to launch
   INTEGER       :: icc        ! Timestep count for model integration
   INTEGER       :: madtt      ! Number of model time steps in a precip data time step
   INTEGER       :: tid        ! Thread number
   REAL          :: jay        ! Disaggregation counter
   REAL          :: miss       ! Reanal missing data value
   REAL          :: iitt       ! max time steps to track parcels
   REAL          :: itt        ! Time counter
   REAL          :: ssum       ! Tally of sourced parcels
   REAL          :: xstart     ! Longitude of first reanal grid box
   REAL          :: dx,dy      ! rean grid spacing
   REAL          :: tirean     ! reanal data interval (hours)
   REAL          :: dint       ! reanal data interval (seconds)
   REAL          :: tpd        ! reanal data interval (times per day)
   REAL          :: ce         ! circumference of earth
   REAL          :: idt        ! fraction of rean data intervals per trajectory intervals
   REAL          :: slack      ! Maximum number of days outside (before) pentad to track parcels
   REAL          :: etot       ! total precip during pentad
   REAL          :: eint       ! precip represented by each parcel
   REAL          :: lon000,lat000 ! Lon and lat of launch grid box (center)
   REAL          :: tmax,tmin  ! time limits for interpolation in pentad
   REAL          :: enow       ! precip this time step
   REAL          :: ecnt       ! precip unaccounted for by parcels (remainder)
   REAL          :: rndlon,rndlat ! quasi-random offsets for parcel launches
   REAL          :: rndpw      ! quasi-random level for parcel launches (0-1)
   REAL          :: player(izmax) ! PW in each sigma layer
   REAL          :: pwp        ! total PW in column (lowest 16 layers)
   REAL          :: pwat       ! quasi-random level for parcel launches (units of PW)
   REAL          :: pwq        ! PW stepping up through layers
   REAL          :: L00        ! potential temperature for parcel release
   REAL          :: lon1,lat1  ! location after advection
   REAL          :: lev1       ! level after advection
   REAL          :: tstep      ! number of time steps per reanalysis data interval
                               ! (negative for back trajectories)
   REAL          :: stday      ! start day (latest day of lwindo for back traj)
   REAL          :: pmiss      ! precip missing data value
   REAL          :: offset     ! account for averaging period in reanal data
   REAL          :: totint     ! total time steps required to track parcels
   REAL          :: ipp         ! Real time step for hybrid precip data
   REAL          :: fac0,fac1   ! Weighting factors for time interp of sigma-level data
   REAL          :: fal0,fal1   ! Weighting factors for time interp of flux data
   REAL          :: fap0,fap1,fap2   ! Weighting factors for time interp of precip data
   REAL          :: lon0 (pnmax)
   REAL          :: lat0 (pnmax)
   REAL          :: vol  (pnmax) 
   REAL          :: lev0 (pnmax)
   REAL          :: sink(ixmax,iymax) ! Grid of evaporation source
   REAL          :: frac        ! Fraction of parcel contributed by evaporation at a time step 
   REAL          :: pqmin       ! Minimum threshold for parcel water content (%)
   REAL          :: psum,p      ! Total precip for the pentad
   REAL          :: denom, numer
!
   REAL          :: thybp(pstep)       ! hybp(:,ix0,iy0)      
   REAL          :: pseries(0:pstep+1) ! Will be set to same as thyby(pstep) between 1 : pstep
   REAL          :: rmadtt             ! REAL(madtt), -tstep/REAL(ipdiur/irtpd)
!
!>>>>  Set some basic constants
!
      dura = lwindo * irtpd           ! # reanal timesteps in launch window
      ssum = 0.0                      ! tally parcels
      miss = -9.99e+33                ! reanal missing value
      pmiss = 9999.0                  ! precip missing value
      xstart = 0.0                    ! Longitude of 1st reanal grid box
      dx = 1.875                      ! delta X (degrees) of reanal data
      dy = (gauss(iymax)-gauss(1))/REAL(iymax-1) ! mean delta Y of reanal data
!
      tirean = 24.0/REAL(irtpd)       ! reanal data interval (hours)
      dint = tirean * 60.0 * 60.0     ! reanal data interval (seconds)
      tpd = 24/(dint/3600.)           ! times per day for reanal data
      ce = 2 * re * pi                ! circumference of earth
      idt = 1.0/tstep                 ! fraction of rean data intervals per trajectory intervals
      iitt = 1.0
      if (idt < 0.0) iitt = REAL(ftmax) - 1.0 ! if back-trajectory,
                                              ! work backward from end of period
      itt = iitt
      offset = 0.5000000              ! account for averaging period in reanal data
      tid = ofil - 30
!
!>>>>> Zero out the buckets, choose terminus
!
  pseries = 0.0
  DO ip = 1, pstep
     pseries(ip) = thybp(ip)
  ENDDO
  sink = 0.0
  lon000 = (REAL(x000)-1.0)*dx+xstart   ! Longitude of terminus box : Terminus = Place where to start tracing
  lat000 = gauss(y000)                  ! Latitude  of terminus box
  if (idt < 0) then
    it2 = int(itt)                      ! End   parcel launch window in reanal data
    it1 = int(it2) - dura + 1           ! Start parcel launch window in reanal data
  else
    return
  endif
!  print *,'itt=',itt,', lon=',lon000,', lat=',lat000
!  print *,'it1=',it1,', it2=',it2
!
!>>>>> Start and stop intervals for hourly precip data
!
  if (idt < 0) then
    ip2 = pstep
    ip1 = pstep - ipdiur * lwindo + 1
  else
    return
  endif
! print *,'ip1=',ip1,', ip2=',ip2
!
!>>>> Calculate total precip into box during duration
!
  etot = SUM(pseries)
! print *,pseries
  WRITE (6,1000) tid,x000,y000,etot,dura
!
!>>>>> Calculate precip launch interval
!
  eint = etot / REAL(pnmax) ! Precipitation per parcel
  ecnt = 0.5 * eint
  pcnt = 0    ! Parcel count
!
!>>>>>
!>>>>> TIME LOOP
!>>>>>
!
  psum = 0.0
  rmadtt = -tstep/REAL(ipdiur/irtpd)
  madtt = INT(rmadtt)
  tmax = itt                          ! Give slack beyond duration
  tmin = 1.0                          ! Farthest back to go
! print *,'TMIN = ',tmin

  DO WHILE (itt >= tmin)              ! itt is in reanal timesteps,
                                      ! but counts in increments of model time steps

    IF (eint == 0.0) EXIT             ! If no precip, no need to trace
                                      ! Exit!!!!!!!!

    icc = - NINT(itt/ idt)            ! timestep count for model integration
    ipp =  (itt-REAL(it1)+1)*REAL(ipdiur/irtpd)    ! record count for precip data
    i0 = INT(itt-0.0001)              ! i0, i1 are indices for instantaneous (sigma) data 
    i1 = i0 + 1
    l0 = INT(itt+offset-0.0001)       ! l0, l1 are indices for t-avg (flux) data
    l1 = l0 + 1
    p0 = INT(ipp-0.0001)              ! p0, p1, p2 are indices for hybrid precip data
    p1 = p0 + 1
    p2 = p1 + 1
    fac1 = itt - REAL(i0) 
    fac0 = 1.0 - fac1
    fal1 = itt - REAL(l0) + offset
    fal0 = 1.0 - fal1
!   fap1 = ipp - REAL(p0) 
!   fap0 = 1.0 - fap1
    jay =  (1 - ipp + REAL(p0))*rmadtt + 1
    jay =  (ipp - REAL(p0))*rmadtt
    fap1 = (2.0*rmadtt-ABS(2.0*jay-madtt-1.0))/(rmadtt*2.0)
    fap0 = MAX(1.0-(jay*2.0+madtt-1.0)/(rmadtt*2.0),0.0)
    fap2 = MAX(1.0-((madtt+1.0-jay)*2.0+madtt-1.0)/(rmadtt*2.0),0.0)
!   write(6,2000) itt,ipp,icc,i0,i1,fac0,fac1,l0,l1,fal0,fal1,p0,p1,p2,fap0,fap1,fap2
! 2000 format('ITT=',f7.2,', IPP=',F6.2,', ICC=',i4,', i0 i1=',2i3,', fac0 fac1=',2f7.4/35x,'l0 l1=',2i3,', fal0 fal1=',2f7.4/29x,'p0 p1 p2=',3i3,', fap0 fap1 fap2=',3f7.4)
!
!>>>>> Disaggregate precipitation in a mass-conserving fashion
!
    p = 0.0
    IF (ipp > 0.0) THEN
!   IF (ipp > 0.0) p = -(pseries(p0)*fap0+pseries(p1)*fap1)*adtt
      denom = (0.5*(pseries(p0)+pseries(p2))+3.0*pseries(p1))*rmadtt
      numer = 4.0*pseries(p1)  
!     print *,denom,pseries(p0),pseries(p1),pseries(p2),numer
!
!     My correction to formatting
!
!     IF (denom > EPSILON(denom)) &
!     p = (pseries(p0)*fap0+pseries(p1)*fap1+pseries(p2)*fap2) * numer / denom

      IF (denom > EPSILON(denom)) THEN
         p = (pseries(p0)*fap0+pseries(p1)*fap1+pseries(p2)*fap2) * numer / denom
      ENDIF
!
      psum = psum + p
!     print *,ipp,'PSUM=',psum,' NUMER=',numer,' DENOM=',denom/rmadtt
!
    ENDIF

!
!>>>>> Associate temporal subsets of reanalysis data
!
    lh0 => evap(:,:,i0)
    lh1 => evap(:,:,i1)
!
!p  if (itt<31.0) print *,'-1',ps0(100,21) 
    ps0 => psfc(:,:,i0)
    ps1 => psfc(:,:,i1)
!   print *,i0,i1,lh0(x000,y000),lh1(x000,y000)
!p  if (itt<31.0) print *,'-2',ps0(100,21) 
!
    q0 => qhum(:,:,:,i0)
    q1 => qhum(:,:,:,i1)
!
    t0 => tmpr(:,:,:,i0)
    t1 => tmpr(:,:,:,i1)
!
    u0 => uwnd(:,:,:,i0)
    u1 => uwnd(:,:,:,i1)
!
    v0 => vwnd(:,:,:,i0)
    v1 => vwnd(:,:,:,i1)
!    if (itt<75.0) print *,'-2',u0(151,69,8)
!
!>>>>> If all parcels not yet released
!
    IF (pcnt < pnmax) THEN
      enow = p
      ecnt = ecnt + enow
!     write(6,2010)enow,ecnt
 2010 format(' Now=',f9.5,', Ecnt=',f9.5)
!
!>>>>> Tag new parcels
!
      DO WHILE (ecnt > eint)
        pcnt = pcnt + 1
!       print *,pcnt
        rndlon = ((10*ecnt + itt - INT(10*ecnt+itt))-0.5) * dx
        rndlat = ((itt+ecnt - INT(itt+ecnt))-0.5) * dy
        ecnt = ecnt - eint             ! decrement running total of rain mass
!       print *,'RND',ecnt,eint,rndlon,rndlat
        lon0(pcnt) = lon000 + rndlon   ! Location of parcel launch
        lat0(pcnt) = lat000 + rndlat
        vol(pcnt) = 1.0                ! The bag is full to start
!       print *,pcnt,lon0(pcnt),lat0(pcnt)
!
!>>>>> First, find precipitable water for this location
!
        player = 0.0    ! PW is each sigma layer
        FORALL(iz=1:izmax) &
          player(iz) = dlev(iz) * 1e-3 &
                * (q0(x000,y000,iz)*ps0(x000,y000)*fac0 &
                + q1(x000,y000,iz)*ps1(x000,y000)*fac1)
!       print *,'PLayer=',player
        pwp = SUM(ABS(player))/grav  ! Total precipitable water (in lowest 16 layers)
!
!>>>>> Randomly select a humidity-weighted level
!
        rndpw = 1000.0 * ABS(rndlon*rndlat)
        rndpw = rndpw - INT(rndpw)
        pwat = pwp * rndpw
!       print*,'PW=',pwp,' rand=',rndpw
        player = 0.0
        DO islev = 1,izmax-1
!         print *,'islev=',islev,player(islev),q0(x000,y000,islev),ps0(x000,y000),q1(x000,y000,islev),ps1(x000,y000)
          FORALL(iz=1:islev) &
            player(iz) = dlev(iz) * 1e-3 &
                 * (q0(x000,y000,iz)*ps0(x000,y000)*fac0 &
                 + q1(x000,y000,iz)*ps1(x000,y000)*fac1)
          pwq = SUM(ABS(player))/grav
!         print*,pwq,pwat
          IF (pwq > pwat) EXIT
        END DO
!
!>>>>> Determine initial theta for parcel release
!
        L00 = (t0(x000,y000,islev)*fac0+t1(x000,y000,islev)*fac1) / &
              ((ps0(x000,y000)*fac0+ps1(x000,y000)*fac1)*lev(islev)*1e-8)**(2./7.)
        lev0(pcnt) = L00
!       if ((x000==49).and.(y000==22)) &
!       write(6,2020)pcnt,lon0(pcnt),lat0(pcnt),lev0(pcnt),lev(islev)
 2020   format('Parcel ',i4,' released at ',f7.2,',',f7.2,', Theta=',f7.2, &
               ', Sig level=',f6.0)
!        print *,ecnt,eint
      END DO
    END IF
!
!>>>>> For each active parcel
!
!   print *,pact
    pact = 0
    DO pn = 1,pcnt
!
!>>>>>   Is parcel sucked dry?
!
!     print *,'VOL',vol(pn),pn
      IF (vol(pn) > pqmin*0.01) THEN
        pact = pact + 1
!
!>>>>> Check for rain/ET
!
!       if (itt<75.0) print *,'Pre-EVAP',pn
!       if (itt<75.0) print *,ps0(100,21)
        CALL isevap(tid,pn,lon0(pn),lat0(pn),itt,idt,fal0,fal1,fac0,fac1, &
                  tpd,gauss,ixmax,iymax,izmax,dx,xstart,dlev,grav,pqmin, &
                  ssum,vol(pn),ilon,ilat,frac)
!       write(6,4000)tid,itt,ssum
 4000   format(i2,'> itt=',f7.2,', ssum=',f7.2)
!
!>>>>> Dump into bucket
!
        sink(ilon,ilat) = sink(ilon,ilat) + frac
!       if (itt<75.0) print *,'Post-EVAP'
!
!>>>>> Advect parcel
!
        lon1 = 0.0
        lat1 = 0.0
        CALL blowS(lon0(pn),lat0(pn),lev0(pn),lon1,lat1,lev1, &
                 itt,idt,dint,fac0,fac1,ce,dx,dy,gauss,ixmax,iymax,izmax,xstart,lev)
!        if (itt<75.0) print *,'Post-BLOW'
!V2        if ((lon1<182.0).or.(lon1>358.0).or.(lat1<2.0).or.(lat1>88.0)) vol(pn) = 0.1
        IF (lat1<-90.0) THEN
           lat1 = -180.0 - lat1
           lon1 = lon1 + 180.0
        ENDIF
        IF (lat1>90.0) THEN
           lat1 = 180.0 - lat1
           lon1 = lon1 + 180.0
        ENDIF
        IF (lon1<0.0) lon1 = lon1 + 360.0
        IF (lon1>360.0) lon1 = lon1 - 360.0
        lon0(pn) = lon1
        lat0(pn) = lat1
        lev0(pn) = lev1
      ENDIF
    END DO
!   print *,pact,'parcels active'
!
!>>>> Increment time
    itt = itt + idt
    IF (pcnt == pnmax .AND. pact == 0) EXIT

  END DO      !======== END TIME LOOP
!
!>>>>> Write stats
!
  totint = iitt - itt
  sink = sink * etot / pnmax
  write (6,1010) tid,x000,y000,ssum,pnmax,totint
  1000 format(i2,'> X=',i3,', Y=',i3,' Total precip:',f9.4,' during ',i3, &
              ' time steps.')
  1010 format(i2,'> X=',i3,', Y=',i3,':',f7.2,' of ',i4,' parcels sank during ',f8.2, &
              ' time steps.')
!
  sink(1,1) = 1.0
  sink(2,1) = etot
  sink(3,1) = pnmax
  sink(4,1) = ssum
  sink(5,1) = dura
!
  WRITE(ofil,REC=orec) sink
RETURN
!
!
!
!
END SUBROUTINE onebox
!
!
!
!
SUBROUTINE blowS(lonS,latS,levS,lon1,lat1,lev1, &
                 itt,idt,dint,fac0,fac1,ce,dx,dy,gauss,ixmax,iymax,izmax,xstart,lev)
!                          
  USE pointers_mod
! Using u0,u1,v0,v1 from pointers_mod
!
  IMPLICIT NONE
  REAL, INTENT(IN)      :: lonS,latS,levS       ! Lon, lat, lev of parcel at beginning of time step
  INTEGER, INTENT(IN)   :: ixmax,iymax,izmax    ! Dimensions of 3-D grid 
  REAL, INTENT(IN)      :: gauss(iymax)         ! Gaussian latitudes for reanalysis grid
  REAL, INTENT(IN)      :: lev(izmax)           ! Sigmas of reanalysis model levels
  REAL, INTENT(IN)      :: itt                  ! Reanal time step
  REAL, INTENT(IN)      :: idt                  ! fraction of rean data intervals per trajectory intervals
  REAL, INTENT(IN)      :: dint                 ! reanal data interval (seconds)
  REAL, INTENT(IN)      :: ce                   ! circumference of earth        
  REAL, INTENT(IN)      :: dx                   ! Grid width in X direction
  REAL, INTENT(IN)      :: dy                   ! Mean grid width in Y direction
  REAL, INTENT(IN)      :: xstart               ! Longitude for 1st reanalysis grid box
  REAL, INTENT(IN)      :: fac0,fac1            ! Weights for time interpolation
!
  REAL, INTENT(OUT)     :: lon1,lat1,lev1       ! Lon, lat, lev of parcel at end of time step
!
  INTEGER               :: ilon, ilat      ! Indices on 2-D grid for nearest point.
  REAL                  :: slon, slat      ! Lon and lat of center of nearest point on 2-D grid.
  REAL                  :: sig,la0,lb0,fln0,fln1,flt0,flt1
  REAL                  :: d00,d01,d10,d11      ! Distances to 4 surrounding grid centers
  REAL                  :: denom                ! Norm. factor for distance weighting
  REAL                  :: uuu,vvv              ! Interp'd wind components for back advect
  REAL                  :: uuN,vvN              ! Interp'd wind components for forward advect
  REAL                  :: latN,lonN,levN       ! Interim lon, lat, lev 
  REAL                  :: dlat,dlon            ! Distance parcel advected in lat and lon
  INTEGER               :: isig,a0,b0,a1,b1
  REAL                  :: npole,spole          ! Northern and southern limits to avoid pole singularity
!
  npole = (gauss(iymax)+gauss(iymax-1))*0.5
  spole = (gauss(1)+gauss(2))*0.5
!
!>>>>> Isentropic advection using sigma-level data
!
  CALL sigma(lonS,latS,levS,levN,isig,sig, &
                 gauss,ixmax,iymax,izmax,dx,xstart,fac0,fac1,lev)
                   ! Determine sigma level of parcel
                   ! Returns levN, isig and sig
!    if (itt<75.0) print *,'    Post-SIGMA S'
!
!>>>>> Find nearest grid point SW of parcel
!
  CALL nearpt(lonS,latS,slon,slat,ilon,ilat, &
                  gauss,ixmax,iymax,dx,xstart)
! Why do this? Would not just give la0=slon? .. Unless slon > 360.0; tracing around the globe more than once
  la0 = MOD(slon+360.0,360.0)
  a0  = MOD(ilon+ixmax-1,ixmax)+1
  IF (slon > lonS) THEN
    la0 = MOD(la0-dx+360.0,360.0)
    a0  = MOD(a0-1+ixmax-1,ixmax)+1
  ENDIF         
  lb0 = slat
  b0 = ilat
  IF (slat > latS) THEN
    lb0 = lb0 - dy
    b0 = b0 - 1
!   print *,'slat',slat,latS,lb0,b0
  ENDIF
  IF (lb0 >= npole) THEN
    lb0 = 2 * npole - lb0
    b0  = 2 * iymax - 1 - b0
    la0 = MOD(la0+180.0,360.0)
    a0 = MOD(a0+ixmax/2-1,ixmax)+1
  ENDIF
  IF (lb0 <= spole) THEN
    lb0 = 2 * spole - lb0
    b0  = 3 - b0
    la0 = MOD(la0+180.0,360.0)
    a0 = MOD(a0+ixmax/2-1,ixmax)+1
!   print *,'SPOLE',ilat,b0,slat,lb0,a0
  ENDIF
!
!>>>> Set up interpolation weights and find NE point
!
  a1 = MOD(a0+1+ixmax-1,ixmax)+1
  fln1 = lonS-slon
  IF (fln1 < 0.0) fln1 = fln1 + dx
  fln0 = dx - fln1
  b1 = b0 + 1
  IF (b1 >= iymax) b1 = 2 * iymax - 1 - b1
  flt1 = latS-slat
  IF (flt1 < 0.0) flt1 = flt1 + dy
  flt0 = dy - flt1
!
!>>>> calculate distances for horizontal interpolation
!
  d00 = SQRT(fln1*fln1+flt1*flt1)+1e-8
  d01 = SQRT(fln1*fln1+flt0*flt0)+1e-8
  d10 = SQRT(fln0*fln0+flt1*flt1)+1e-8
  d11 = SQRT(fln0*fln0+flt0*flt0)+1e-8
  denom = 1/d00+1/d01+1/d10+1/d11
!p    if (itt<31.0) print *,'    Forward',denom,d00,d01,d10,d11
!
!>>>> space and time interp of winds calculated
!
!p    if (itt<31.0) print *,a0,a1,b0,b1,isig,fac0,fac1
!p    if (itt<31.0) print *,u0(a0,b0,isig)
!p    if (itt<31.0) print *,u0(a1,b0,isig)
!p    if (itt<31.0) print *,u0(a0,b1,isig)
!p    if (itt<75.0) print *,u0(a1,b1,isig)
  uuu = ((u0(a0,b0,isig)*fac0 + u1(a0,b0,isig)*fac1)/d00 &
      + (u0(a1,b0,isig)*fac0 + u1(a1,b0,isig)*fac1)/d10 &
      + (u0(a0,b1,isig)*fac0 + u1(a0,b1,isig)*fac1)/d01 &
      + (u0(a1,b1,isig)*fac0 + u1(a1,b1,isig)*fac1)/d11) / denom 
!p    if (itt<75.0) print *,'    Post-uuu'
  vvv = ((v0(a0,b0,isig)*fac0 + v1(a0,b0,isig)*fac1)/d00 &
      + (v0(a1,b0,isig)*fac0 + v1(a1,b0,isig)*fac1)/d10 &
      + (v0(a0,b1,isig)*fac0 + v1(a0,b1,isig)*fac1)/d01 &
      + (v0(a1,b1,isig)*fac0 + v1(a1,b1,isig)*fac1)/d11) / denom
!    if (itt<75.0) print *,'    Post-TERP S'
!
!  Advect from intial location
!
  dlon = uuu * idt * dint * (360/ce)
  lonN = MOD(lonS+dlon+360.0,360.0)
  dlat = vvv * idt * dint * (360/ce)
  latN = latS + dlat
  IF (latN >= npole) THEN
    latN = 2 * npole - latN
    lonN = MOD(lonN+180.0,360.0)
  ENDIF
!p    if (itt<75.0) print *,'    Post-Pole 1'
  IF (latN <= spole) THEN
    latN = 2 * spole - latN
    lonN = MOD(lonN+180.0,360.0)
  ENDIF

!
!>>>>> latN and lonN are forward-interp location
!>>>>> Next, find U and V at this new location
!
! Post-inital advect location
!
  CALL sigma(lonN,latN,levN,lev1,isig,sig, &
                 gauss,ixmax,iymax,izmax,dx,xstart,fac0,fac1,lev)
                   ! Determine sigma level of parcel
                   ! Returns lev1, isig and sig
!    if (itt<75.0) print *,'    Post-SIGMA N',lev1,isig,sig
!
!>>>>> Find nearest grid point SW of parcel
!
  CALL nearpt(lonN,latN,slon,slat,ilon,ilat, &
                  gauss,ixmax,iymax,dx,xstart)
  la0 = MOD(slon+360.0,360.0)
  a0 = MOD(ilon+ixmax-1,ixmax)+1
  IF (slon > lonN) THEN
    la0 = MOD(la0-dx+360.0,360.0)
    a0 = MOD(a0-1+ixmax-1,ixmax)+1
  ENDIF         
  lb0 = slat
  b0 = ilat
  IF (slat > latN) THEN
    lb0 = lb0 - dy
    b0 = b0 - 1
  ENDIF
  IF (lb0 >= npole) THEN
    lb0 = 2 * npole - lb0
    b0  = 2 * iymax - 1 - b0
    la0 = MOD(la0+180.0,360.0)
    a0 = MOD(a0+ixmax/2-1,ixmax)+1
  ENDIF
  IF (lb0 <= spole) THEN
    lb0 = 2 * spole - lb0
    b0  = 3 - b0
    la0 = MOD(la0+180.0,360.0)
    a0 = MOD(a0+ixmax/2-1,ixmax)+1
!    print *,'SPOL2',ilat,b0,slat,lb0,a0
  ENDIF
!
!>>>> Set up interpolation weights and find NE point
!
  a1 = MOD(a0+1+ixmax-1,ixmax)+1
  fln1 = lonN-slon
  IF (fln1 < 0.0) fln1 = fln1 + dx
  fln0 = dx - fln1
  b1 = b0 + 1
  IF (b1 >= iymax) b1 = 2 * iymax - 1 - b1
  flt1 = latN-slat
  IF (flt1 < 0.0) flt1 = flt1 + dy
  flt0 = dy - flt1
!
!>>>> calculate distances for horizontal interpolation
!
  d00 = SQRT(fln1*fln1+flt1*flt1)+1e-8
  d01 = SQRT(fln1*fln1+flt0*flt0)+1e-8
  d10 = SQRT(fln0*fln0+flt1*flt1)+1e-8
  d11 = SQRT(fln0*fln0+flt0*flt0)+1e-8
  denom = 1/d00+1/d01+1/d10+1/d11
!
!>>>> space and time interp of winds calculated
!
  uuN = ((u0(a0,b0,isig)*fac0 + u1(a0,b0,isig)*fac1)/d00 &
      + (u0(a1,b0,isig)*fac0 + u1(a1,b0,isig)*fac1)/d10 &
      + (u0(a0,b1,isig)*fac0 + u1(a0,b1,isig)*fac1)/d01 &
      + (u0(a1,b1,isig)*fac0 + u1(a1,b1,isig)*fac1)/d11) / denom
!p    if (itt<75.0) print *,'    Between'
  vvN = ((v0(a0,b0,isig)*fac0 + v1(a0,b0,isig)*fac1)/d00 &
      + (v0(a1,b0,isig)*fac0 + v1(a1,b0,isig)*fac1)/d10 &
      + (v0(a0,b1,isig)*fac0 + v1(a0,b1,isig)*fac1)/d01 &
      + (v0(a1,b1,isig)*fac0 + v1(a1,b1,isig)*fac1)/d11) / denom
!
!>>>> Two U's and V's calculated - solve implicit calculation
!
!    if (itt<75.0) print *,'    IMPLICIT'
  dlon = 0.5*(uuu+uuN) * idt * dint * (360/ce)
  lon1 = MOD(lonS+dlon+360.0,360.0)
  dlat = 0.5*(vvv+vvN) * idt * dint * (360/ce)
  lat1 = latS + dlat
  IF (lat1 >= npole) THEN
    lat1 = 2 * npole - lat1
    lon1 = MOD(lon1+180.0,360.0)
  ENDIF
  IF (lat1 <= spole) THEN
    lat1 = 2 * spole - lat1
    lon1 = MOD(lon1+180.0,360.0)
!   print *,'SPOLZ',spole,lat1,lon1
  ENDIF
!
END SUBROUTINE blowS
!
!
!
SUBROUTINE isevap(tid,pn,lon0,lat0,itt,idt,fal0,fal1,fac0,fac1, &
                  tpd,gauss,ixmax,iymax,izmax,dx,xstart,dlev,grav,pqmin, &
                  ssum,vol0,ilon,ilat,frac)
!
  USE pointers_mod
! Using lh0,lh1,ps0,ps1,q0,q1 from pointers_mod
  IMPLICIT NONE
!
  INTEGER, INTENT(IN)   :: tid         ! Thread number
  INTEGER, INTENT(IN)   :: pn          ! Parcel number
  REAL, INTENT(IN)      :: lon0,lat0   ! Exact parcel latitude and longitude
  REAL, INTENT(IN)      :: itt         ! Reanal time step
  REAL, INTENT(IN)      :: idt         ! fraction of rean data intervals per trajectory intervals
  REAL, INTENT(IN)      :: fal0,fal1   ! Weights for precip time interpolation
  REAL, INTENT(IN)      :: fac0,fac1   ! Weights for met time interpolation
  REAL, INTENT(IN)      :: tpd         ! times per day for reanal data
  REAL, INTENT(IN)      :: gauss(iymax)! Gaussian latitudes for reanalysis grid
  REAL, INTENT(IN)      :: dlev(izmax) ! Sigma layer thicknesses
  REAL, INTENT(IN)      :: grav        ! Gravity 
  REAL, INTENT(IN)      :: pqmin       ! Minimum threshold for parcel water content (%)
!
  INTEGER, INTENT(IN)   :: ixmax,iymax,izmax ! 3-D dimensions of reanalysis grid
  REAL, INTENT(IN)      :: dx                   ! Grid distance in X direction
  REAL, INTENT(IN)      :: xstart               ! Longitude for 1st reanalysis grid box
!
  REAL, INTENT(INOUT)   :: ssum        ! Total accumulated precip
  REAL, INTENT(INOUT)   :: vol0        ! volume remaining in parcel
!
  INTEGER, INTENT(OUT)  :: ilon, ilat  ! Indices on 2-D grid for nearest point.
  REAL, INTENT(OUT)     :: frac        ! Fraction contributed by evaporation from this point
!
  REAL                  :: slon, slat      ! Lon and lat of center of nearest point on 2-D grid.
  REAL                  :: player(izmax)
  REAL                  :: rip, eee, pwp, ratio, vol1, pfrac
  integer               :: ip, iz 
!
!>>>>> Depletes back-trajectory parcel by evaporation
! 
  pfrac = pqmin * 0.01
  vol1 = vol0
  frac = 0.0
!
  rip = -NINT((itt-1.0)/idt) + 1
  ip = INT(rip+0.5)            ! Hourly count for precip data
  CALL nearpt(lon0,lat0,slon,slat,ilon,ilat, &
                  gauss,ixmax,iymax,dx,xstart)
!
!>>>>> Is there evaporation?
!
!p  if (itt < 71.) print *,ilon,ilat
!p  if (itt < 71.) print *,fal0,fal1
!p  if (itt < 71.) print *,lh0(ilon,ilat),lh1(ilon,ilat),fal0,fal1
  eee = -(lh0(ilon,ilat)*fal0+lh1(ilon,ilat)*fal1)*0.03455 &
        * idt / tpd
!p  if (itt < 71.) print *,eee
  IF (eee <= 0.0) RETURN
  player=0.0
!p  if (itt < 31.) print *,q0(ilon,ilat,:)
!p  if (itt < 71.) print *,ps0(ilon,ilat)
!p  if (itt < 31.) print *,q1(ilon,ilat,:)
!p  if (itt < 75.) print *,ps1(ilon,ilat)
  FORALL(iz=1:izmax) &
     player(iz) = dlev(iz) * 1e-3 &
                * (q0(ilon,ilat,iz)*ps0(ilon,ilat)*fac0 &
                + q1(ilon,ilat,iz)*ps1(ilon,ilat)*fac1)
!p  if (itt < 75.) print *,player
  pwp = SUM(ABS(player))/grav
!
  ratio = eee/pwp            ! Fraction of PW contrib by evap 
  IF (vol1 - ratio < pfrac) ratio = vol1 -pfrac 
  frac = ratio * 1.0 / (1.0-pfrac)  
  vol1 = vol1 - ratio
  ssum = ssum + frac
!  WRITE(6,3000)pn,vol1,ratio,eee,pwp,ilon,ilat
 3000 format(i4,': Vol=',f7.4,', Ratio=',f7.4,', E=',f9.4,', PW=',f9.4, &
      ' at X=',i3,', Y=',i3)
  vol0 = vol1
!
END SUBROUTINE isevap
!
!
!

SUBROUTINE sigma(lon0,lat0,lev0,lev1,isig,sig, &
                 gauss,ixmax,iymax,izmax,dx,xstart,fac0,fac1,lev)
!
!
  USE pointers_mod
! Using t0,t1,ps0,ps1 from pointers_mod
  IMPLICIT NONE
  INTEGER, INTENT(IN)   :: ixmax,iymax,izmax    ! Dimensions of 3-D grid 
  REAL, INTENT(IN)      :: lon0,lat0,lev0       ! Exact parcel longitude, latitude, and sigma level
  REAL, INTENT(IN)      :: gauss(iymax)         ! Gaussian latitudes for reanalysis grid
  REAL, INTENT(IN)      :: lev(izmax)           ! Sigmas of reanalysis model levels
  REAL, INTENT(IN)      :: dx                   ! Grid distance in X direction
  REAL, INTENT(IN)      :: xstart               ! Longitude for 1st reanalysis grid box
  REAL, INTENT(IN)      :: fac0,fac1            ! Weights for time interpolation
!
  REAL, INTENT(OUT)     :: lev1, sig 
  INTEGER, INTENT(OUT)  :: isig
!
  INTEGER               :: ilon, ilat      ! Indices on 2-D grid for nearest point.
  REAL                  :: slon, slat      ! Lon and lat of center of nearest point on 2-D grid.
  REAL                  :: prof(izmax)     ! Theta in P=P(:)
  REAL                  :: odel,del,sgn,olev,slev,tmin,l00
  INTEGER               :: iz

!
!>>>>> Determines theta and sigma level of parcel
!
  l00 = lev0
  CALL nearpt(lon0,lat0,slon,slat,ilon,ilat, &
                  gauss,ixmax,iymax,dx,xstart)
!
!>>>>> Determine surface pressure (?, Theta?) at boundary times
!
  FORALL(iz=1:izmax) &
    prof(iz) = t0(ilon,ilat,iz)*fac0/(ps0(ilon,ilat)*lev(iz)*1e-8)**(2./7.) &
             + t1(ilon,ilat,iz)*fac1/(ps1(ilon,ilat)*lev(iz)*1e-8)**(2./7.)
  tmin = minval(prof)            ! Lowest potential temp in profile
  l00 = max(l00,tmin+1)          ! Adjust parcel theta to 1K above lower than minimum in prof(:) if necessary
                                 ! prof(1) = surface
  lev1 = l00

  olev = 995.0
  odel = 0.0
!
!>>>>> Determine nearest sigma level for this theta
!
  DO iz=1,izmax-1
    isig = iz
    slev = lev(iz)
    del = l00 - prof(iz)
    sgn = del + odel
!p  IF (itt<31.0) print *,iz,slev,del,sgn,prof(iz)
    IF (del <= 0.0) THEN
      IF (sgn <= 0.0) THEN
        slev = olev
        isig = isig - 1
        EXIT
      ENDIF
    ENDIF
    olev = slev
    odel = del
  ENDDO 
  sig = slev
!
END SUBROUTINE sigma
!
!
!
!
SUBROUTINE nearpt(lon0,lat0,slon,slat,ilon,ilat, &
                  gauss,ixmax,iymax,dx,xstart)
! expects: gauss, ixmax, iymax, dx, xstart
!
  IMPLICIT NONE
  INTEGER, INTENT(IN)   :: ixmax, iymax    ! Dimensions of 2-D grid 
  REAL, INTENT(IN)      :: lon0, lat0      ! Exact parcel latitude and longitude
  REAL, INTENT(IN)      :: gauss(iymax)    ! Gaussian latitudes for reanalysis grid
  REAL, INTENT(IN)      :: dx              ! Grid distance in X direction
  REAL, INTENT(IN)      :: xstart          ! Longitude for 1st reanalysis grid box
  REAL, INTENT(OUT)     :: slon, slat      ! Lon and lat of center of nearest point on 2-D grid.
  INTEGER, INTENT(OUT)  :: ilon, ilat      ! Indices on 2-D grid for nearest point.
  REAL odel,del,sgn,olon,olat
  INTEGER iy,ix
!
  odel = 0.0
  olat = gauss(1)
!
!>>>>> Determine nearest latitude index for lat0
!
  DO iy=1,iymax
    ilat = iy
    slat = gauss(iy)
    del = lat0 - slat
    sgn = del + odel
    IF (del <= 0.0) THEN
      IF (sgn <= 0.0) THEN
        slat = olat
        ilat = iy - 1
      ENDIF
      EXIT
    ENDIF
    olat = slat
    odel = del
  ENDDO 
!
!>>>>> Determine nearest longitude index for lon0
!
  odel = 0.0
  olon = xstart
!
  DO ix=1,ixmax
    ilon = ix
    slon = xstart + dx * REAL(ix-1)
    del = lon0 - slon
    sgn = del + odel
    IF (del <= 0.0) THEN
      IF (sgn <= 0.0) THEN
        slon = olon
        ilon = ix - 1
      ENDIF
      EXIT
    ENDIF
    olon = slon
    odel = del
  ENDDO 
!
END SUBROUTINE nearpt