subroutine blendwind(cgrid,glat,glon,iadate,iadate2,u,v,w,nx,ny,nhwrf)
!***********************************************************************
! abstract: blend together wind fields coming from different models    *
!           originally used to blend rap winds with hwrf winds in the  * 
!           event of a tropical storm                                  *
!                                                                      *
! program history log:                                                 *
!   2012-10-25  pondeca                                                *
!***********************************************************************
 use kinds, only: r_single,i_kind
 use constants, only: init_constants,init_constants_derived

 implicit none

!**************************************************************
! input variables:
! cgrid             - ndfd grid name
! nx,ny             - dimensions of ndfd grid
! glat(nx,ny)       - ndfd earth latitudes in deg
! glon(nx,ny)       - ndfd earth longitudes in deg 
!                     Range: ]-180 , +180]
! iadate(5)         - input analysis date (iyear,imonth,iday,ihour,iminute)
! iadate2(5)        - input cycle date for hwrf
! u(nx,ny),v(nx,ny) - original ndfd wind 
! w(nx,ny)          - original ndfd wind speed 
! nhwrf             - hwrf storm index number. note that there
!                     can be several hwrf storms at a given time 

! output variables
! u(nx,ny),v(nx,ny) - ndfd wind after blending with hwrf wind
!**************************************************************
!**************************************************************
 
!Declare passed variables
 character(60)       , intent(in)    :: cgrid
 integer(i_kind)     , intent(in)    :: nx,ny
 integer(i_kind)     , intent(in)    :: iadate(5)
 integer(i_kind)     , intent(in)    :: iadate2(5)
 integer(i_kind)     , intent(in)    :: nhwrf
 real(r_single)      , intent(in)    :: glat(nx,ny)
 real(r_single)      , intent(in)    :: glon(nx,ny)
 real(r_single)      , intent(inout) :: u(nx,ny)
 real(r_single)      , intent(inout) :: v(nx,ny)
 real(r_single)      , intent(inout) :: w(nx,ny)

!Declare local parameters
 logical,parameter::regional=.true.

!Declare local variables
 real(r_single),allocatable,dimension(:,:,:):: ud2,vd2,rd2,wd2

 integer(i_kind) n,i,j
 integer(i_kind) nminanl,nmins2,nrecs
 real(r_single) rdiff
 logical hexist(2)
 character(1)  ch1
 character(80) fname
 character(2) clun

 print*,'in blendwind: cgrid,nx,ny=',trim(cgrid),nx,ny

 allocate(ud2(nx,ny,2))
 allocate(vd2(nx,ny,2))
 allocate(rd2(nx,ny,2))
 allocate(wd2(nx,ny,2))

 hexist(1:2)=.false.

 write (clun,'(i2.2)') nhwrf

 inquire(file='hwrfin1_'//clun,exist=hexist(1)) !hwrf outer grid
 inquire(file='hwrfin2_'//clun,exist=hexist(2)) !hwrf inner grid

 if (.not.hexist(1) .and. .not.hexist(2)) then 
    print*,'in blendwind: could not find hwrf storm number =', nhwrf
    return
 endif

 print*,'in blendwind: processing hwrf storm number ',nhwrf

 call init_constants_derived
 call init_constants(regional)

!convert analysis time and hwrf cycle time to minutes relative to fixed date
 call w3fs21(iadate,nminanl)
 call w3fs21(iadate2,nmins2)

 rdiff=real((nminanl-nmins2),r_single)
 nrecs=nint((rdiff+0.1_r_single)/60_r_single)+1

 print*,'in blendwind: analysis date,minutes ',iadate,nminanl
 print*,'in blendwind: hwrf date,minutes     ',iadate2,nmins2
 print*,'in blendwind: nrecs=',nrecs

 do n=1,2
    if (hexist(n))  then 
       write(ch1,'(i1.1)') n
       fname='hwrfin'//ch1//'_'//clun
       call hwrfwind_on_ndfd(cgrid,fname,ud2(1,1,n),vd2(1,1,n), & 
                             rd2(1,1,n),glon,glat,nrecs,nx,ny)  
    endif
 enddo

 print*,'in blendwind: before blending: umin,max=',minval(u),maxval(u)
 print*,'in blendwind: before blending: vmin,max=',minval(v),maxval(v)
  
 call blend0(cgrid,u,ud2,rd2,hexist,nx,ny,2)
 call blend0(cgrid,v,vd2,rd2,hexist,nx,ny,2)

 print*,'in blendwind: after blending: umin,max=',minval(u),maxval(u)
 print*,'in blendwind: after blending: vmin,max=',minval(v),maxval(v)

 do n=1,2
   do j=1,ny
     do i=1,nx
       wd2(i,j,n)=ud2(i,j,n)*ud2(i,j,n)+vd2(i,j,n)*vd2(i,j,n) 
       if ( wd2(i,j,n) > 0.) wd2(i,j,n)=sqrt(wd2(i,j,n))
     enddo
   enddo
 enddo

 call blend0(cgrid,w,wd2,rd2,hexist,nx,ny,2)

 deallocate(ud2)
 deallocate(vd2)
 deallocate(rd2)
 deallocate(wd2)
 end subroutine blendwind
!====================================================================
!====================================================================
 subroutine hwrfwind_on_ndfd(cgrid,fname,uout,vout,rout, & 
            glon,glat,nrecs,nx,ny)

 use kinds, only: r_kind,r_single,i_kind
 use constants, only: deg2rad,zero
 use generalized_transf, only: init_dims,init_general_transform_2, &
                               tll2xy,rotate_wind_xy2ll

 implicit none

!Declare passed variables

 character(60)   , intent(in)  :: cgrid
 character(80)   , intent(in)  :: fname
 integer(i_kind) , intent(in)  :: nx,ny
 integer(i_kind) , intent(in)  :: nrecs
 real(r_single)  , intent(in)  :: glon(nx,ny)
 real(r_single)  , intent(in)  :: glat(nx,ny)
 real(r_single)  , intent(out) :: uout(nx,ny)
 real(r_single)  , intent(out) :: vout(nx,ny)
 real(r_single)  , intent(out) :: rout(nx,ny)

!Declare local parameters:
 real(r_single),parameter::spval=-9999._r_single
 real(r_kind),parameter:: r360 = 360.0_r_kind
 integer(i_kind),parameter::lun=1
 integer(i_kind),parameter::npass=10

!Declare local variables
 integer(i_kind) nxh,nyh
 real(r_single),allocatable,dimension(:,:):: uh,vh,acph
 real(r_single),allocatable,dimension(:,:):: glonh,glath
 real(r_kind)  ,allocatable,dimension(:,:):: glonh8,glath8

 integer(i_kind) i,j,m,n
 real(r_single) xx,yy,xx1,yy1
 real(r_kind)   xx8,yy8
 real(r_kind)   u8,v8
 real(r_kind)   u08,v08
 real(r_kind)   rlon08,rlat08
 logical outside

 open(lun,file=trim(fname),form='unformatted')
 read(lun) nxh,nyh

 print*,'in hwrfwind_on_ndfd: input file is ',trim(fname)
 print*,'in hwrfwind_on_ndfd: nxh,nyh=',nxh,nyh

 allocate(uh(nxh,nyh))
 allocate(vh(nxh,nyh))
 allocate(glonh(nxh,nyh))
 allocate(glath(nxh,nyh))
 allocate(acph(nxh,nyh))

 do n=1,nrecs
    read(lun) uh
    read(lun) vh
    read(lun) glath
    read(lun) glonh
    read(lun) acph
 enddo
 close(lun)

 open(lun,file=trim(fname)//'.grads',form='unformatted')
 write(lun) uh
 write(lun) vh
 write(lun) glath
 write(lun) glonh
 close(lun)

 allocate(glonh8(nxh,nyh))
 allocate(glath8(nxh,nyh))

 glath8=real(glath*deg2rad,r_kind)
 glonh8=real(glonh*deg2rad,r_kind)

 print*,'in hwrfwind_on_ndfd: glonh8min,max=', & 
            minval(glonh8),maxval(glonh8)

 print*,'in hwrfwind_on_ndfd: glath8min,max=', & 
            minval(glath8),maxval(glath8)

 call init_dims(nxh,nyh)
 call init_general_transform_2(glath8,glonh8)

 do j=1,nyh
 do i=1,nxh
    u8=real(uh(i,j),r_kind)
    v8=real(vh(i,j),r_kind)
    rlon08=real(glonh(i,j),r_kind)
    if (rlon08 >= r360) rlon08=rlon08-r360
    if (rlon08  < zero) rlon08=rlon08+r360
    rlon08=rlon08*deg2rad
    rlat08=real(glath(i,j),r_kind)*deg2rad
    
    call tll2xy(rlon08,rlat08,xx8,yy8,outside)
    call rotate_wind_xy2ll(u8,v8,u08,v08,rlon08,rlat08,xx8,yy8)

    uh(i,j)=real(u08,r_single)
    vh(i,j)=real(v08,r_single)
 enddo
 enddo

 call ndfd_rel_hwrfwinds(cgrid,glonh,uh,vh,nxh,nyh)

 uout=0.
 vout=0.
 rout=spval

 do j=1,ny
 do i=1,nx
    rlon08=real(glon(i,j),r_kind)
    if (rlon08 >= r360) rlon08=rlon08-r360
    if (rlon08  < zero) rlon08=rlon08+r360
    rlon08=rlon08*deg2rad
    rlat08=real(glat(i,j),r_kind)*deg2rad

    call tll2xy(rlon08,rlat08,xx8,yy8,outside)

    if (.not. outside) then
       xx=real(xx8,r_single)
       yy=real(yy8,r_single)
       call bilinear_2d0(uh,nxh,nyh,uout(i,j),yy,xx)!Note the reverse order "yy,xx"
       call bilinear_2d0(vh,nxh,nyh,vout(i,j),yy,xx)!Note the reverse order "yy,xx"
       xx1=real(nxh,r_single)-xx
       yy1=real(nyh,r_single)-yy
       rout(i,j)=min(xx,xx1,yy,yy1)
    endif
 enddo
 enddo

 call smther_one_2(uout,1,nx,1,ny,npass)
 call smther_one_2(vout,1,nx,1,ny,npass)

 open(lun,file=trim(fname)//'.grads2',form='unformatted')
 write(lun) uout
 write(lun) vout
 close(lun)

 deallocate(uh)
 deallocate(vh)
 deallocate(glath)
 deallocate(glonh)
 deallocate(acph)
 deallocate(glath8)
 deallocate(glonh8)


 end subroutine hwrfwind_on_ndfd
!====================================================================
!====================================================================
 subroutine ndfd_rel_hwrfwinds(cgrid,glonh,uh,vh,nxh,nyh)

 use kinds, only: r_single,i_kind
 use constants, only: deg2rad
 implicit none

 include 'param.incl'

!Declare passed variables
 character(60)   , intent(in)    :: cgrid
 integer(i_kind) , intent(in)    :: nxh,nyh
 real(r_single)  , intent(in)    :: glonh(nxh,nyh)
 real(r_single)  , intent(inout) :: uh(nxh,nyh)
 real(r_single)  , intent(inout) :: vh(nxh,nyh)

!Declare local parameters:
 real(r_single),parameter:: s180 = 180.0_r_single
 real(r_single),parameter:: s360 = 360.0_r_single
 
!Declare local variables
 real(r_single) xn,elonv,ylon,ue,ve,u,v,sdeg2rad
 real(r_single) angle2,sinx2,cosx2
 integer(i_kind) i,j

 sdeg2rad=real(deg2rad,r_single)

 if (trim(cgrid)=='conus' .or. trim(cgrid)=='alaska' .or. &
  trim(cgrid)=='cohres' .or. trim(cgrid)=='akhres' .or. trim(cgrid)=='hrrr' .or. & 
  trim(cgrid)=='cohresext' .or. trim(cgrid)=='juneau') then

    if (trim(cgrid)=='conus') then 
        elonv=real(elonv_conus,r_single)
        xn=sin(elonv*sdeg2rad)
    endif

    if (trim(cgrid)=='alaska') then 
        elonv=real(elonv_alaska,r_single)
        xn=1._r_single
    endif

    if (trim(cgrid)=='cohres') then 
        elonv=real(elonv_cohres,r_single)
        xn=sin(elonv*sdeg2rad)
    endif

    if (trim(cgrid)=='akhres') then 
        elonv=real(elonv_akhres,r_single)
        xn=1._r_single
    endif

    if (trim(cgrid)=='hrrr') then 
        elonv=real(elonv_hrrr,r_single)
        xn=sin(elonv*sdeg2rad)
    endif

    if (trim(cgrid)=='cohresext') then 
        elonv=real(elonv_cohresext,r_single)
        xn=sin(elonv*sdeg2rad)
    endif

    if (trim(cgrid)=='juneau') then 
        elonv=real(elonv_juneau,r_single)
        xn=sin(elonv*sdeg2rad)
    endif

    do j=1,nyh
    do i=1,nxh
       ue=uh(i,j)
       ve=vh(i,j)
       ylon=glonh(i,j)
       if (ylon  < 0._r_single) ylon=ylon+s360
       angle2=xn*(ylon-elonv)*sdeg2rad
       sinx2 = sin(angle2)
       cosx2 = cos(angle2)
       uh(i,j)=cosx2*ue-sinx2*ve
       vh(i,j)=sinx2*ue+cosx2*ve
    enddo
    enddo
 
 endif
end subroutine ndfd_rel_hwrfwinds
!====================================================================
!====================================================================
 subroutine blend0(cgrid,field1,field2,rd2,hexist,nx,ny,nk)

 use kinds, only: r_single,i_kind
 implicit none

!Declare passed variables
 character(60)  ,intent(in)  :: cgrid
 integer(i_kind),intent(in)  ::nx,ny,nk
 real(r_single),intent(inout)::field1(nx,ny)
 real(r_single),intent(in)   ::field2(nx,ny,nk)
 real(r_single),intent(in)   ::rd2(nx,ny,nk)
 logical,intent(in)          ::hexist(nk)

!Declare local parameters
 real(r_single),parameter::spval=-9999._r_single
 real(r_single),parameter::epsilon=1.e-3_r_single

!Declare local variables
 real(r_single),allocatable,dimension(:,:):: fldaux
 real(r_single) alpha1,alpha2,diff
 real(r_single) s1,s2,rl1,rl2,amp1,amp2
 integer(i_kind) i,j,n
 integer(i_kind) i1,i2,j1,j2
 integer(i_kind) istart,iend,jstart,jend

 data s1/2._r_single/
 data s2/4._r_single/
 data rl1/7._r_single/
 data rl2/7._r_single/
 data amp1/1._r_single/
 data amp2/1._r_single/

 if (trim(cgrid)=='conus' .or. trim(cgrid)=='alaska') then
    rl2=4._r_single
 endif

 if (trim(cgrid)=='cohres' .or. trim(cgrid)=='cohresext' .or. & 
     trim(cgrid)=='akhres' .or. trim(cgrid)=='hrrr') then
    rl2=8._r_single 
 endif

 if (trim(cgrid)=='juneau') then
    rl2=13._r_single 
 endif 

 allocate(fldaux(nx,ny))

 if (hexist(1)) then !start with outer domain 
    do j=1,ny
    do i=1,nx
       diff=abs(rd2(i,j,1)-spval)
       if (diff < epsilon) then  
           fldaux(i,j)=field1(i,j)
         else
           alpha1=amp1*exp(-rd2(i,j,1)*rd2(i,j,1)/(s2*rl1*rl1))
           fldaux(i,j)=alpha1*field1(i,j)+(1.-alpha1)*field2(i,j,1)
        endif
    enddo
    enddo

    do j=1,ny
    do i=1,nx
       field1(i,j)=fldaux(i,j)
    enddo
    enddo
 endif

 if (hexist(2)) then !now blend with inner domain
    i1=+huge(i1)
    i2=-huge(i1)
    j1=+huge(j1)
    j2=-huge(j1)
    do j=1,ny
    do i=1,nx
       diff=abs(rd2(i,j,2)-spval)
       if (diff < epsilon) then  
           fldaux(i,j)=field1(i,j)
         else
           alpha2=amp2*exp(-rd2(i,j,2)*rd2(i,j,2)/(s2*rl2*rl2))
           fldaux(i,j)=alpha2*field1(i,j)+(1.-alpha2)*field2(i,j,2)
           i1=min(i1,i) ; i2=max(i2,i)
           j1=min(j1,j) ; j2=max(j2,j)
        endif
    enddo
    enddo
    print*,'in blendwind: i1,i2,j1,j2=',i1,i2,j1,j2

    do j=1,ny
    do i=1,nx
      field1(i,j)=fldaux(i,j)
    enddo
    enddo

!   if ( i1<nx .and. i2<nx .and. j1<ny .and. j2<ny .and.i2>i1 .and. j2>j1 ) then
!
!      print*,'in blendwind: 5-point smoothing'
!      istart=max(1,i1-20) ; iend=min(nx,i2+20)
!      jstart=max(1,j1-20) ; jend=min(ny,j2+20)
!
!      print*,'in blendwind: istart,iend,jstart,jend=',istart,iend,jstart,jend
!
!      do n=1,5
!        do j=jstart,jend
!        do i=istart,iend
!           fldaux(i,j)=0.5*field1(i,j)+0.125*(field1(i-1,j-1)+field1(i-1,j+1)+field1(i+1,j-1)+field1(i+1,j+1))
!!          fldaux(i,j)=0.25*(field1(i-1,j-1)+field1(i-1,j+1)+field1(i+1,j-1)+field1(i+1,j+1))
!        enddo
!        enddo
!
!        do j=1,ny
!        do i=1,nx
!          field1(i,j)=fldaux(i,j)
!        enddo
!        enddo
!      enddo
!   endif
 endif

 deallocate(fldaux)
end subroutine blend0
!==============================================================
!************************************************************
!------------------------------------------------------
      subroutine bilinear_2d0(rffcst,ix,jx,rfobs,xx,yy)
!
!      input: rffcst (model grid value.)
!      output: rfobs (interpolated value)
!      const: xx,yy (define coordinates in grid units
!                    of point for which interpolation is
!                    performed)
!
!
!     i+1,j |          | i+1,j+1
!         --+----------+---
!           |          | dym
!           |    *     + -
!           |   x,y    | dy
!           |          |
!         --+----+-----+---
!        i,j|<dx>|<dxm>| i,j+1
!*******************************************************************
      use kinds, only: r_single,i_kind

      implicit none

!declare passed variables
      integer(i_kind),intent(in):: ix,jx
      real(r_single),intent(in)::rffcst(ix,jx)
      real(r_single),intent(out):: rfobs

!declare local variables
      integer(i_kind) i,j
      real(r_single) xx,yy,dx,dy,dxm,dym

      i  = ifix(yy)
      j  = ifix(xx)
      if((i.ge.1) .and. (i.le.(ix-1)) .and. &
         (j.ge.1) .and. (j.le.(jx-1)) ) then
        dx = xx - float(j)
        dy = yy - float(i)
        dxm= 1.0-dx
        dym= 1.0-dy
          rfobs=dxm*(dym*rffcst(i,j)+dy*rffcst(i+1,j)) &
                   + dx *(dym*rffcst(i,j+1)+dy*rffcst(i+1,j+1))
      else
        rfobs=huge(rfobs)        
      endif

      return
      end
!------------------------------------------------------
subroutine smther_one_2(g1,is,ie,js,je,ns)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:   smther_one
! prgmmr: pondeca          org: np22                date: 2006-08-01
!
! abstract: apply 1-2-1 smoother in each direction of data slab
!
! program history log:
!   2006-08-01  pondeca
!
!   input argument list:
!    g1        - 2d array of field to be smoothed
!    is,ie     - first and last i values of data slab
!    js,je     - first and last j values of data slab
!    ns        - number of passes
!
!   output argument list:
!    g1        - smoothed 2d field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!

  use kinds,only: r_single,i_kind
  implicit none

  integer(i_kind)  is, ie, js, je
  integer(i_kind)  i,j,l,ip,im,jp,jm
  integer(i_kind), intent(in) :: ns

  real(r_single), dimension(is:ie, js:je), intent(inout) :: g1
                                   !  on input: original data slab
                                   !  on ouput: filtered data slab
 

  real(r_single), allocatable:: g2(:,:)

   allocate(g2(is:ie,js:je))
   do l=1,ns

     do j=js,je
      do i=is,ie
       ip=min(i+1,ie) ; im=max(is,i-1)
         g2(i,j)=.25*(g1(ip,j)+g1(im,j))+.5*g1(i,j)
      end do
     end do

     do i=is,ie
      do j=js,je
       jp=min(j+1,je) ; jm=max(js,j-1)
       g1(i,j)=.25*(g2(i,jp)+g2(i,jm))+.5*g2(i,j)
      end do
     end do

   end do
   deallocate(g2)

   return
end subroutine smther_one_2
!------------------------------------------------------