!*==sbbr.spg  processed by SPAG 6.70Rc at 15:33 on 10 Dec 2012
!*------------------   SPAG Configuration Options   --------------------
!*--0323,76 000101,-1 000000102011332010100002000000210211210,136 10 --
!*--1100000011112111000000000000,10,10,10,10,10,10,900,100 200000000 --
!*--000000010000000000000,72,72 73,42,38,33 00011112110000100000000  --
!*----------------------------------------------------------------------
!$Author: joelevin $
!$Date: 2011-03-24 11:33:26 -0500 (Thu, 24 Mar 2011) $
!$Revision: 11724 $
!$HeadURL: https://eweru-dev1.eweru.ksu.edu/svn/code/weps1/branches/WEPS_F90_update/weps.src/src/lib_erosion/sbbr.for $
!**********************************************************************
!     subroutine sbbar       3/20/98
!**********************************************************************
      subroutine sbbr
      use i_p1werm
      use i_m1geo
      use i_p1const
      use i_e2grid
      use i_e3grid
      use i_m2geo
      use s_ang
      use s_fu
      use s_slen
      implicit none
!*--SBBR25
!
!*** Start of declarations rewritten by SPAG
!
! Local variables
!
      real :: a,aa,alpha,b,c,ca,ceta,dmin,lx,ly,tmpa,tmpw0,w,waa
      real,dimension(8) :: awar
      integer :: i,j,k,n,tmpv
!
!*** End of declarations rewritten by SPAG
!
! ^^^ note must add wlen to call for tst version
!     + + + PURPOSE + + +
!     to calculate the fraction of open field friction velocity
!     in shelter for 8 cardinal wind directions 0, 45, ...315
!     at all interior nodes
!
!     to assign the 8 fractons calculated at each node to as 3-d
!     array (W0br(i,j,k)) for all nodes inside sim. region.
 
!     + + + KEYWORDS + + +
!     wlen
 
!     + + + PARAMETERS + + +
 
!     + + + ARGUMENT DEFINITIONS  + + +
!     wlen = windward(-) and leeward(+) distances from barrier
!            calculated for each wind direction;
!  ^^^ used only in test version
 
!     + + + GLOBAL COMMON BLOCKS + + +
 
 
!     + + + LOCAL COMMON BLOCKS + + +
 
!     + + + ARGUMENT DECLARATIONS + + +
! ^^^ used only in test version
!      real wlen(0:mngdpt, 0:mngdpt, 8)
 
!     + + + LOCAL VARIABLES + + +
 
 
 
 
!     + + + LOCAL VARIABLE DEFINITIONS + + +
 
!     i, j, k, n = do-loop indices
!     a,c = lengths to barrier ends from cell x(i,j) (m)
!     b = barrier lenght (m)
!     aa, ca = angles opposite lenths a and c.
!     waa = angle between barrier and wind direction
!     alpha = clockwise angle from d (0 deg) to a
!     ceta  = clockwise angle from d (0 deg) to c
!     dmin = minimum distance from barrier to x(i,j)
!     w = length from barrier to x(i,j) along wind direction
!     tmpa = temporary (intermediate) angle calculation
 
!     + + + SUBROUTINES CALLED + + +
!     none
 
!     + + + FUNCTION DECLARATIONS + + +
 
!     + + + DATA INITIALIZATIONS + + +
 
!     initialize variable (fraction of open field fric. vel)
      do i = 0,imax
         do j = 0,jmax
            do k = 1,8
               w0br(i,j,k) = 1.0
            end do
         end do
      end do
 
!     + + + END SPECIFICATIONS + + +
 
      ! calc. 8 cardinal wind directions relative to the simulation
      ! region using a relative y-axis orientation of 0 deg for region.
 
      do i = 1,8
         awar(i) = (i-1)*45 - amasim
         if (awar(i)<0) then
            awar(i) = awar(i) + 360
         else if (awar(i)>360) then
            awar(i) = awar(i) - 360
         end if
 
        ! convert relative wind angles to radians
         awar(i) = awar(i)*pi/180
      end do
 
      ! update interior nodes
      do i = 1,imax - 1
         do j = 1,jmax - 1
          !*2   calculate distance to grid points(maybe offset from origin)
            lx = (i-0.5)*ix + amxsim(1,1)
            ly = (j-0.5)*jy + amxsim(2,1)
! ^^^tmp
!      if (i .eq. 1 .and. j .eq. jmax-1)then
!       write (*,*) 'sbbr output at 1, jmax-1 lx=', lx, 'ly=', ly
!       write (*,*) 'ix=',ix, 'jy=', jy, 'imax=', imax, 'jmax=',jmax
!      endif
! ^^^ end tmp
          ! barrier sweep
            do n = 1,nbr
            ! find distances a, c from x(i,j) to barrier ends
               a = slen(lx,ly,amxbr(1,1,n),amxbr(2,1,n))
               c = slen(lx,ly,amxbr(1,2,n),amxbr(2,2,n))
 
            ! find barrier length
               b = slen(amxbr(1,1,n),amxbr(2,1,n),amxbr(1,2,n),         &
                 & &amxbr(2,2,n))
 
            ! for a grid cell on a barrier calc. fric. vel. fraction
               if ((a+c)<=b) then
                  do k = 1,8
                     w = 0
                     tmpw0 = w0br(i,j,k)
                     w0br(i,j,k) = fu(w,amzbr(n),ampbr(n),amxbrw(n))
                     w0br(i,j,k) = min(tmpw0,w0br(i,j,k))
! ^^^ used only in test version
!        wlen(i,j,k) = 0
 
                  end do
                  cycle
               end if
 
            ! find angles from x(i,j) to barrier ends
               tmpa = ((b*b+c*c-a*a)/(2*b*c))
               if (tmpa<-1) then
                  tmpa = -1
               else if (tmpa>1) then
                  tmpa = 1
               end if
               aa = acos(tmpa)
               tmpa = ((a*a+b*b-c*c)/(2*a*b))
               if (tmpa<-1) then
                  tmpa = -1
               else if (tmpa>1) then
                  tmpa = 1
               end if
               ca = acos(tmpa)
 
            ! find minimum distance to barrier
               if (aa>90*pi/180.OR.ca>90*pi/180) then
                  dmin = min(a,c)
               else
                  dmin = b*sin(aa)*sin(ca)/sin(aa+ca)
               end if
 
!     is x(i,j) within 35*barrier height of barrier?
! ^^^tmp
!      write (*,*) 'sbbr.for output'
!      write (*,*) 'dmin=', dmin, 'i=',i,'j=',j
! ^^^ end tmp
 
               if (dmin<=35*amzbr(n)) then
              ! Find angles alpa and ceta
                  alpha = ang(amxbr(1,1,n),amxbr(2,1,n),lx,ly)
                  ceta = ang(amxbr(1,2,n),amxbr(2,2,n),lx,ly)
              ! sweep relative wind directions
                  do k = 1,8
                ! if grid cell downwind from barrier make tmpv = 1
                     tmpv = 0
                     if (abs(alpha-ceta)<pi) then
                        if ((awar(k)>=min(alpha,ceta))&.AND.            &
                          & (awar(k)<max(alpha,ceta))) tmpv = 1
                     else if (abs(alpha-ceta)>=pi) then
                        if (awar(k)<=min(alpha,ceta)) tmpv = 1
                        if ((awar(k)>=max(alpha,ceta))&.AND.            &
                          & (awar(k)<=2*pi)) tmpv = 1
                     end if
 
! ^^^tmp used for test only
!      if (i .eq. 34 .and. j .eq. 1) then
!      write(*,*) ' barrier tmp output'
!      write (*,*) 'a=',a, ' b=', b, ' c=',c
!       write(*,*) 'aa=', aa, ' ca=',ca
!       write(*,*) 'dmin=', dmin
!       write(*,*) 'alpha=', alpha, '  ceta=', ceta
!       write(*,*) 'tmpv=', tmpv
!       write(*,*) 'awar(K)=', awar(k), ' k=', k
!      endif
! ^^^ end tmp
!
                ! calc. opposite side 'w'
                     if (tmpv>0) then
                  ! calc. barrier to x(i,j) distance 'w' along wind vector
                  ! calc. opposite side 'w'
                        tmpa = abs(awar(k)-alpha)
                        if (tmpa>pi) tmpa = 2*pi - tmpa
                        waa = pi - tmpa - ca
                        waa = max(waa,0.0001)
                                          !edit LH 3-27-07 prevent zero
                  ! by sin rule calc w
                        w = sin(ca)*a/sin(waa)
 
                  ! calc. min fraction unsheltered friction vel.
                        tmpw0 = w0br(i,j,k)
                        w0br(i,j,k) = fu(w,amzbr(n),ampbr(n),amxbrw(n))
                        w0br(i,j,k) = min(tmpw0,w0br(i,j,k))
 
! ^^^^  tmp used only in test version
!            wlen(i,j,k) = w
 
                     else if (dmin<5*amzbr(n)) then
                  ! is x(i,j) within 5H of barrier
                  ! calc. frac. friction vel. for other directions
                        tmpw0 = w0br(i,j,k)
                        w0br(i,j,k) = fu(-dmin,amzbr(n),ampbr(n),amxbrw(&
                                    & n))
                        w0br(i,j,k) = min(tmpw0,w0br(i,j,k))
 
! ^^^     tmp used only in test version
!              wlen(i,j,k) = -dmin
 
                     end if
                  end do
               end if
            end do
         end do
      end do
! ^^^tmp output
!      n = imax/15
!      n = max(1,n)
!      write (*,*)
!      write (*,*) 'BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB'
!      write (*,*) 'output from sbbr.for, w0br(i,j,k)'
!      write (*,*) '    i index values'
!      write (*,310) (i, i=0,imax,n)
!      do 220 k = 1,8
!      do 220 j = jmax,0,-1
!      write (*,300) (j, (w0br(i,j,k),i=0,imax,n))
!  220 continue
!      write (*,*) 'BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB'
!  300 format (1x, i3,1x, 30(f4.2,1x))
!  310 format (1x, i8, 30i5)
! ^^^ end tmp
 
      end subroutine sbbr
!*==slen.spg  processed by SPAG 6.70Rc at 15:33 on 10 Dec 2012
!*------------------   SPAG Configuration Options   --------------------
!*--0323,76 000101,-1 000000102011332010100002000000210211210,136 10 --
!*--1100000011112111000000000000,10,10,10,10,10,10,900,100 200000000 --
!*--000000010000000000000,72,72 73,42,38,33 00011112110000100000000  --
!*----------------------------------------------------------------------
 
!     function to calc. length
      function slen(x1,y1,x2,y2)
      implicit none
!*--SLEN275
!
!*** Start of declarations rewritten by SPAG
!
! Dummy arguments
!
      real :: x1,x2,y1,y2
      real :: slen
      intent (in) x1,x2,y1,y2
!
!*** End of declarations rewritten by SPAG
!
      slen = abs(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)))
      end function slen
!*==ang.spg  processed by SPAG 6.70Rc at 15:33 on 10 Dec 2012
!*------------------   SPAG Configuration Options   --------------------
!*--0323,76 000101,-1 000000102011332010100002000000210211210,136 10 --
!*--1100000011112111000000000000,10,10,10,10,10,10,900,100 200000000 --
!*--000000010000000000000,72,72 73,42,38,33 00011112110000100000000  --
!*----------------------------------------------------------------------
 
!     function to calc quadrant and then angles alpha and ceta
!     line direction point to br; ang clockwise from north
      function ang(xbr,ybr,xp,yp)
      implicit none
!*--ANG300
!
!*** Start of declarations rewritten by SPAG
!
! Dummy arguments
!
      real :: xbr,xp,ybr,yp
      real :: ang
      intent (in) xbr,xp,ybr,yp
!
! Local variables
!
      real :: pi,tempa
!
!*** End of declarations rewritten by SPAG
!
      pi = 3.1415927
      if (xbr/=xp) then
         tempa = atan(abs((ybr-yp)/(xbr-xp)))
         if ((ybr==yp).AND.(xbr>xp)) then         !horizontal line east
            ang = pi/2.0
         else if ((ybr==yp).AND.(xbr<xp)) then      !horiz. line west
            ang = 1.5*pi
         else if ((ybr>yp).AND.(xbr>xp)) then     !quad 1
            ang = pi/2.0 - tempa
         else if ((ybr<yp).AND.(xbr>xp)) then      !quad 2
            ang = pi/2.0 + tempa
         else if ((ybr<yp).AND.(xbr<xp)) then      !quad 3
            ang = 1.5*pi - tempa
         else if ((ybr>yp).AND.(xbr<xp)) then       !quad 4
            ang = 1.5*pi + tempa
         end if
      else if (ybr<yp) then     ! vertical line to south
         ang = pi
      else
         ang = 2*pi        ! vertical line to north
      end if
 
      end function ang
!*==fu.spg  processed by SPAG 6.70Rc at 15:33 on 10 Dec 2012
!*------------------   SPAG Configuration Options   --------------------
!*--0323,76 000101,-1 000000102011332010100002000000210211210,136 10 --
!*--1100000011112111000000000000,10,10,10,10,10,10,900,100 200000000 --
!*--000000010000000000000,72,72 73,42,38,33 00011112110000100000000  --
!*----------------------------------------------------------------------
 
!     function to calc. fu (fraction upwind fric. velocity
!     near the  barrier)
!     (ranges: porosity 0 to 0.9, distance: -5*zbr to 50*zbr)
      function fu(xh,zbr,pbr,xbrw)
      implicit none
!*--FU351
!
!*** Start of declarations rewritten by SPAG
!
! Dummy arguments
!
      real :: pbr,xbrw,xh,zbr
      real :: fu
      intent (in) pbr,xbrw,xh,zbr
!
! Local variables
!
      real :: a,b,c,d,pb,x,xw
!
!*** End of declarations rewritten by SPAG
!
!     scale distance & width by barrier height
      x = xh/zbr
      xw = xbrw/zbr
!     increase effective porosity with barrier width
      pb = pbr + (1-exp(-0.5*xw))*0.3*(1-pbr)
 
!     calculate coef. as fn of porosity
      a = 0.008 - 0.17*pb + 0.17*pb**1.05
      b = 1.35*exp(-0.5*pb**0.2)
      c = 10*(1-0.5*pb)
      d = 3 - pb
!     calc. frac. of fric. vel.
      fu = 1 - exp(-a*x**2) + b*exp(-0.003*(x+c)**d)
!     Cap fu at 1.0
      if (fu>1.0) fu = 1.0
      end function fu