      subroutine snohires (si,sm,glat1,glon1,list,lusaf,im,jm,nm)
!
!  inputs:  sm (eta land/sea mask)
!           glat (latitude array of eta grid)
!           glon (longitude array of eta grid)
!           list (unit number of printout)
!           lusaf (logical - true:use ims snow and usaf snow, 
!                            false: use ims snow only)
!           dtr   (degrees to radians conversion factor)
!
!  outputs: si (snowdepth and sea-ice on eta grid, see calling
!               subroutine cnsts for more details on si output)
!
      implicit real (a-h, o-z)
!
      real,intent(out) :: si(0:im+1,0:jm+1,nm)
      real,intent(in) :: sm(0:im+1,0:jm+1,nm)
      real,intent(in) :: glat1(0:im+1,0:jm+1,nm), glon1(0:im+1,0:jm+1,nm)
      logical,intent(inout) :: lusaf
      integer,intent(in) :: im,jm,nm
!
      integer::i,j,n
!
      real         afsno   (512,512)
      integer      mskaf   (512,512)
!
      real         scvh    (1024,1024)
      integer      mskscvh (1024,1024)
      real,parameter::dtr=3.1415926/180.
!
      real::glat(0:im+1,0:jm+1,nm), glon(0:im+1,0:jm+1,nm)
!
!  input units (45-nesdis daily snow/ice,
!               41-usaf daily snow/ice,
!               43-nesdis 1/16-bedient land/sea mask,
!               42-usaf   1/8-bedient land/sea mask)
! 
      data   insnoims/45/, insnoaf/41/, inmskaf/42/, inmskims/43/
!
!  set lower threshold for snodepth in tenths of inches (mininum
!  snow depth is 1.5" where there is snow coverage).
!
      data   idepth /15/
!
!**************************  begin execution ***********************
!
      glon=-glon1
      glat=glat1

!     specify the unit number of the printer.
      ioutuprt = list
!
!     specify parameters of the n.h. 1024x1024 ims 1/16-mesh grid
!
!     location of the pole:
      xpnmc8 = 513.
      ypnmc8 = 512.
!     n -- the number of grid intervals from the pole to the equator:
      ennmc8 = 16.0e0 * 31.2043316e0
!     the longitudinal rotation angle of the grid
      alnmc8 = 10.e0
!     the orientation west longitude of the grid:
      orient8 = 80.e0
!
!  specify parameters of the n.h. 512 x 512 grid type 88 (usaf grid).
!
!  location of pole:
      xpnmcaf = 257.e0
      ypnmcaf = 256.e0
!  grid mesh length at 60n in km (fyi, but not needed)
!.... xmeshl  = 47.625e0
!  number of grid intervals from pole to equator
      ennmcaf = 8.0 * 31.2043316e0
!  the longitude rotation angle of the grid
      alnmcaf = 10.e0
!  the orientation west longitude of the grid
      orientaf  = 80.e0
!
!  call to read nesdis/ims 1/16-bedient daily n.h. snow/ice via sno16get
!  
      write(6,*) 'calling sno16get'
      call sno16get(scvh,iyear,imonth,iday,insnoims,list)
      write(6,*) 'back from sno16get'
!test
!test      write(60,7701) iyear, imonth, iday
!test  7701 format('ims data : ', 3i2.2)
!test
!
! at present time the ims file has a 2-digit year.  add century:
! now iyear from sno16get is 4-digits so we don't need to do this
!     if (iyear .gt. 90) then
!        iyear = iyear + 1900
!     else
!        iyear = iyear + 2000
!     endif
!
! calculate julian date:
      nday = iw3jdn(iyear,imonth,iday)
      write (ioutuprt, 1) insnoims, iyear, imonth, iday, nday 
    1 format(1h ,' daily nesdis snow read in via unit no.=',i3/  &
             1h ,' for valid date of y-m-d =',i4,'-',i2,'-',i2/  &
             1h ,' which is julian day =', i10) 
!
!   read the nesdis-ims land-water mask
!                                          (sea=0,land=1)
      do 10 j = 1, 1024
         read(inmskims,'(80i1)') (mskscvh(i,j),i=1,1024)
 10   continue
!
      call printims (scvh, mskscvh)
!
!  print specified usaf snow processing flag
      write(ioutuprt,2211)   lusaf
 2211 format(//1h ,5x,'subroutine snohires                    '/  &
           1h ,1x,'  will try to process usaf snowcover: lusaf = ',l2)
!
     
      if (lusaf) then
!
!  read air force 1/8-bedient daily n.h. snow/ice via sno8get
!
      write(6,*) 'calling sno8get'
      call sno8get(afsno,iyear,imonth,iday,insnoaf,list,lusaf)
!test
      write(60,7702) iyear, imonth, iday
  7702 format('usaf data: ', 3i2.2)
!test
!
!  if i/o error encountered in sno8get, then lusaf is returned false
!  
        if (.not. lusaf) then
          write(ioutuprt,56)
   56     format(1h ,'warning: file err in usaf snow anal')
        else 
!  the usaf snow header has a 2-digit year.  add century:
!  now iyear from sno8get is 4-digits so we don't need to do this
!         if (iyear .gt. 90) then
!           iyear = iyear + 1900
!         else
!           iyear = iyear + 2000
!         endif
!
! calculate julian date:
          nafday = iw3jdn(iyear,imonth,iday)
!
          write(ioutuprt,2) insnoaf, iyear, imonth, iday, nafday 
    2     format(1h ,' daily usaf snow read in via unit no.=',i3/  &
                 1h ,' for valid date of y-m-d ='i4,'-',i2,'-',i2/  &
                 1h ,' which is julian day =', i10) 
!
           if (nafday .lt. nday-4) then
             lusaf = .false.
             write(ioutuprt,2644) nday,nafday
 2644        format(/1h ,'******  warning   ****** warning ******'/  &
                  1h ,'date of usaf snow anal more than 4 days old'/  &
                  1h ,' -  will fall-back to nesdis ims snow cover'/  &
                  1h ,43x,'   nesdis snow anal julday=',i7/  &
                  1h ,43x,'     usaf snow anal julday=',i7)
! 
           endif
        endif
      endif
!
!    note:  upon return from call sno8get above, usaf snow/ice field
!            has following physical ranges:
!         - values over sea points are 0.0 or 11.0 (sea-ice flag)
!         - values over land/coast are 0.0 or pos depth in meters
!         - snowdepth over land is actual, not water equivalent yet
!
!   read the usaf afgwc land/coast/sea mask
!                           (sea=1,land=2,coastal-land=4,offworld=9)
      read(inmskaf)  mskaf
!
!---------- i/o of primary input fields is complete -------------
!                        init radians to degrees
      rtd = 1./dtr
!
! idepth is criterion of snowdepth threshold in tenths of
! inches below which usaf snowdepth will be assumed zero.
! here convert idepth to meters
!
       depth = float(idepth) * 2.54e-3   
!
       if ( lusaf) then
         write(ioutuprt,2321) lusaf,idepth
 2321    format(1h //' usaf snodep anal will be used, lusaf=',l2/  &
          1h , 35x,'snodepth threshold (tenths of inches) =',i3)
         call printaf (afsno,mskaf)
      else
         write(ioutuprt,2322) lusaf
 2322    format(1h //' usaf snodep analysis will be ignored'/  &
                1h , 35x,'logical flag lusaf=',l2)
      endif
!
!----------initialize snow/ice arrays to zero on eta grid------------
!
      si = 0.0
!--------------------------------------------------------------------
!
! ****** now begin major loop over all eta grids and points *******
!
         do n=1,nm
         do j=1,jm
         do i=1,im

 !
!--------------- determine lat/lon of eta grid point -------------
!                    (here long will be east long)
!
      yylat = glat(i,j,n)*rtd
      xlong = 360. - glon(i,j,n)*rtd
    
!!!!!!TEST!!!!!!
!      if (xlong.gt.360)then
!      print *,i,j,n,xlong,yylat
!      endif
!!!!!!!!!!!!!TEST!!!!!!!!!!!!



!  where eta domain south of 22 n lat (including any s.h.),
!  we keep default zero snow/ice 
!
!      print *,i,j,n,xlong,yylat,glon(i,j,n),glat(i,j,n),glon1(i,j,n) &
!           ,glat1(i,j,n)
      if (yylat.lt.22.0e0) go to 4300
!----------------------------------------------------------------
!
!  determine location of eta point on the 1024 x 1024 nesdis/ims grid
!
         rm= ennmc8*cos(yylat * dtr) / (1.e0 + sin(yylat * dtr))
         rad = (xlong - alnmc8) * dtr
         x = xpnmc8 + rm * cos(rad)
         y = ypnmc8 + rm * sin(rad)
!
      is  = int(x)
      ip1 = is + 1
      js  = int(y)
      jp1 = js + 1
!
!--------------------------------------------------------------------
!.......first interpolate nesdis snow/ice as primary default......
!
!     the value of snow cover or sea ice on the 1024 x 1024
!     nesdis/ims grid is 1 for snow/ice points, 0 for snow/ice-free.
!     we utilize an 1024 x 1024 land mask
!     from ims to distinguish sea ice from snow points
!     (actually, the data we got from ims have different values for
!     snow and for ice.  but since we want to use as much of the original
!     program [written for the 1/2-mesh sab snow/ice] as possible, we give
!     snow cover and sea ice the same value (=1) and let the mask do the
!     job.
!
!  now use eta and ims land-sea masks to ensure only land
!  points are interpolated to land points (to determine snow)
!  and only sea points are interpolated to sea points (for ice)
!  (nesdis/ims land mask: sea=0,land=1, while the eta mask is the
!  other way round).
!
      iland = 1
      if( sm(i,j,n) .gt. 0.9 ) iland=0
!
      ipoint = nint(x)
      jpoint = nint(y)
!      print *,i,j,n,ipoint,jpoint
      if ( mskscvh(ipoint,jpoint) .eq. iland ) then
        si(i,j,n) = scvh(ipoint,jpoint)
        go to 3351
      endif
!
!  nearest neighbor not same sfc type, so use all 4 surrounding points
!
      kount = 0
!
      xratio = x - real(is)
      yratio = y - real(js)
!
      area11 = (1.0e0 - xratio) * (1.0e0 - yratio)
      area21 = xratio * (1.0e0 - yratio)
      area12 = (1.0e0 - xratio) * yratio
      area22 = xratio * yratio
!
      if( mskscvh(is, js) .eq. iland) then
         kount  = kount + 1
         area   = area11
         ipoint = is
         jpoint = js
      end if
!
      if( mskscvh(is, jp1) .eq. iland ) then
         kount = kount +1
         if (kount .eq. 1) then
            ipoint = is
            jpoint = jp1
         elseif (area12 .gt. area) then
            area   = area12
            ipoint = is
            jpoint = jp1
         end if
      end if
!
      if( mskscvh(ip1, js) .eq. iland ) then
         kount = kount + 1
         if (kount .eq. 1) then
            area   = area21
            ipoint = ip1
            jpoint = js
         elseif (area21 .gt. area) then
            area   = area21
            ipoint = ip1
            jpoint = js
         end if
      end if
!
      if( mskscvh(ip1, jp1) .eq. iland ) then
         kount = kount + 1
         if (kount .eq. 1) then
            area   = area22
            ipoint = ip1
            jpoint = jp1
         elseif (area22 .gt. area) then
            area   = area22
            ipoint = ip1
            jpoint = jp1
         end if
      end if
!
!     determine sno/ice using nearest neighbor with same sfc type 
!
      if(kount .gt. 0) then
         si(i,j,n) = scvh(ipoint,jpoint)
      else
!
!         no immediately surrounding points in the 1024 x 1024 field of
!         snow/ice have the same land-sea type as the eta point.  the
!         eta point may be small island or lake or small bay or pennin.
!         (invariably a small lake in eta grid)
!         so expand search radius and take first sfc type match
!
          ipoint = nint(x)
          jpoint = nint(y)
!
!!!          do 3346  ll=1,4
          do 3346  ll=1,40              !DRAGAN
           jpe = min (1024, jpoint+ll)
           jpb = max (1 , jpoint-ll)
           ipe = min (1024, ipoint+ll)
           ipb = max (1 , ipoint-ll)
!
             do 2346 mk=jpb,jpe
             do 2346 nk=ipb,ipe
               if (mskscvh(nk,mk) .eq. iland) then
               si(i,j,n) = scvh(nk,mk)
               go to 3351
               endif
 2346        continue
 3346     continue
!
!  no land/sea mask matches found, so 
!     a) north of 55n, we assign snow/ice irrespective of sfc type
!     b) south of 55n, we keep a priori zero default
!   (the "b" option best for warmer lats of u.s., where this condition 
!   is virtually always a small eta lake with no counterpart water 
!   nearby in the nesdis/ims grid, e.g., salt lake, where we must  
!   avoid getting sea-ice owing to surrounding snow cover)
!
       write (ioutuprt, 3347) i,j,yylat,xlong,ipoint,jpoint,iland
 3347   format(1h ,'**no ims msk match,',  &
!               'eta-i,j,elat,elon,ims-i,ims-j,lnd:',2i6,2f7.2,2i3,i2)
                'eta-i,j,elat,elon,ims-i,ims-j,lnd:',2i6,2f7.2,2i5,i2)
             if (yylat .ge. 55.0 ) then
 
               si(i,j,n) = scvh(ipoint,jpoint) 
             else
               si(i,j,n) = 0.0
             endif 
    
      endif
!


 3351 continue
! 
      if (.not. lusaf) go to 4300
!
!
!-------------- begin usaf snow/ice interpolation------------------
!  if current usaf snow/ice anal was successfully read (lusaf=true), 
!  add the usaf information as follows:
!
!  the quality of nesdis/ims sea-ice cover and snow coverage (especially
!  in areas with small amounts of snow) are better than the usaf data.
!  we do not use usaf sea-ice cover data at all.  over land, the
!  presence/absence of snow is determined by the ims snow coverage, i.e., 
!  at a given location, if the ims data indicate no snow, then we assume 
!  there is no snow, no matter what the usaf data say.  if the ims data
!  show snow but the usaf data have no snow or less than 1.5" of snow at
!  the location, we assume there is a 1.5" of snow.  if both the ims data
!  and usaf data indicate snow, and the usaf snow depth is more than 1.5",
!  then we use the usaf snow depth.
!
!     on the 512 x 512 usaf grid, the data values are 0.0 for
!     no snow or ice, 11.0 for ice points, and positive but
!     less than 11.0 for snow points.
!
!--------------------------------------if eta sea point, skip usaf anl
      if ( sm(i,j,n) .gt. 0.9 ) go to 4300
!
! if the ims data indicate no snow on this point, also skip usaf anl:
      if (si(i,j,n).lt.0.001) go to 4300
!
! otherwise, set snow depth to be the lower threshold, 1.5":
      si(i,j,n) = depth
!
!-------------------------------------------------------------------
!  this is an eta land point, so apply usaf snowdepth anal
!
!  determine location of eta point on the 512 x 512 usaf grid
!
         rm= ennmcaf*cos(yylat * dtr) / (1.e0 + sin(yylat * dtr))
         rad = (xlong - alnmcaf) * dtr
         x = xpnmcaf + rm * cos(rad)
         y = ypnmcaf + rm * sin(rad)
!
      is  = int(x)
      ip1 = is + 1
      js  = int(y)
      jp1 = js + 1
!
!-----if outside of usaf grid domain (i.e. s.h.) we keep zero default--
!
      if ((is .lt. 1) .or. (is .gt. 511) .or. (js .lt. 1)  &
               .or. (js .gt. 511))  then
        go to 4300
      endif
!--------------------------------------------------------------------
!
!  now use eta and usaf land-sea mask to ensure only land points are 
!  interpolated to land points (to determine snow)
!   (usaf land mask: sea=1,land=2,coastal-land=4,offworld=9)
!   note: in reaching this stage, we have already insured we are on 
!   eta land point and in n.h., i.e. not mskaf=9 (i.e. not offworld)
!  
      iland = 2
!
      ipoint = nint(x)
      jpoint = nint(y)
      if ( mskaf(ipoint,jpoint) .ge. iland) then
        si(i,j,n) = afsno(ipoint,jpoint)
        si(i,j,n) = amax1(depth,si(i,j,n))
        go to 4300
      endif
!
      kount = 0
!
      xratio = x - real(is)
      yratio = y - real(js)
!
      area11 = (1.0e0 - xratio) * (1.0e0 - yratio)
      area21 = xratio * (1.0e0 - yratio)
      area12 = (1.0e0 - xratio) * yratio
      area22 = xratio * yratio
!
      if( mskaf(is, js) .ge. iland) then
         kount  = kount + 1
         area   = area11
         ipoint = is
         jpoint = js
      end if
!
      if( mskaf(is, jp1) .ge. iland ) then
         kount = kount +1
         if (kount .eq. 1) then
            ipoint = is
            jpoint = jp1
         elseif (area12 .gt. area) then
            area   = area12
            ipoint = is
            jpoint = jp1
         end if
      end if
!
      if( mskaf(ip1, js) .ge. iland ) then
         kount = kount + 1
         if (kount .eq. 1) then
            area   = area21
            ipoint = ip1
            jpoint = js
         elseif (area21 .gt. area) then
            area   = area21
            ipoint = ip1
            jpoint = js
         end if
      end if
!
      if( mskaf(ip1, jp1) .ge. iland ) then
         kount = kount + 1
         if (kount .eq. 1) then
            area   = area22
            ipoint = ip1
            jpoint = jp1
         elseif (area22 .gt. area) then
            area   = area22
            ipoint = ip1
            jpoint = jp1
         end if
      end if
!
!     determine sno/ice considering the number of points surrounding
!     eta grid point with the same land-sea mask flag
!
      if (kount .gt. 0) then
          si(i,j,n) = afsno(ipoint,jpoint)
          si(i,j,n) = amax1(depth,si(i,j,n))
!
      else
!
!         no immediately surrounding points in the 512 x 512 field of
!         sno/ice have the same land-sea mask as the eta point.
!         the eta point may be an island or lake or bay or penninsula.
!         so expand search radius and take first mask flag match
!
          ipoint = nint(x)
          jpoint = nint(y)
!
!!!!!          do 7346  ll=1,7
         do 7346  ll=1,70            !DRAGAN  
           jpe = min (512, jpoint+ll)
           jpb = max (1 ,  jpoint-ll)
           ipe = min (512, ipoint+ll)
           ipb = max (1 ,  ipoint-ll)
!
             do 6346 mk=jpb,jpe
             do 6346 nk=ipb,ipe
               if (mskaf(nk,mk) .ge. iland) then
               si(i,j,n) = afsno(nk,mk)
               si(i,j,n) = amax1(depth,si(i,j,n))
               go to 4300
               endif
 6346        continue
 7346     continue
!
!   no land mask match found, so we print warning and stay
!   with earlier value determined from nesdis/ims anal, which 
!   we convert from cover flag (0,1) to default depth of 1.5".
!
          si(i,j,n) = depth
          write (ioutuprt, 7347) i,j,yylat,xlong,is,js,si(i,j,n)
 7347      format(1h ,'*** warning ***..no usaf land msk match ',  &
       ' at eta-i,j,elat,elon,usaf-i,j:',2i6,2f7.2,2i3/  &
       1h ,' default to nesdis value of ',f5.1,' times .10 m')
!
      endif
!     
!
!******************** end major eta grid point loop ****************
!
 4300 continue
         enddo
         enddo
         enddo
!
!  print sample of snow/ice on eta grid
!
      call printeta (si,sm,im,jm,nm)
!test call printyl(si,sm)

      return

      end
