c********************************************************************** c subroutine sbbar 3/20/98 c********************************************************************** subroutine sbbr c ^^^ note must add wlen to call for tst version c + + + PURPOSE + + + c to calculate the fraction of open field friction velocity c in shelter for 8 cardinal wind directions 0, 45, ...315 c at all interior nodes c c to assign the 8 fractons calculated at each node to as 3-d c array (W0br(i,j,k)) for all nodes inside sim. region. c c + + + KEYWORDS + + + c wlen c c + + + PARAMETERS + + + c c c + + + ARGUMENT DEFINITIONS + + + c wlen = windward(-) and leeward(+) distances from barrier c calculated for each wind direction; c ^^^ used only in test version c c + + + GLOBAL COMMON BLOCKS + + + *$noreference include 'p1werm.inc' include 'm1geo.inc' include 'p1const.inc' c c + + + LOCAL COMMON BLOCKS + + + include 'erosion/e2grid.inc' include 'erosion/e3grid.inc' include 'erosion/m2geo.inc' *$reference c c + + + ARGUMENT DECLARATIONS + + + c ^^^ used only in test version c real wlen(0:mngdpt, 0:mngdpt, 8) c c + + + LOCAL VARIABLES + + + c integer i, j, k, n, tmpv c real awar(8), a, b, c, d, e, f, aa, ca, waa real alpha, ceta, dmin, w, tmpa real tmpw0 c c + + + LOCAL VARIABLE DEFINITIONS + + + c c i, j, k, n = do-loop indices c a,c = lengths to barrier ends from cell x(i,j) (m) c b = barrier lenght (m) c d = distance from x(i,j) to top of sim. region,x(i,jmax) (m) c e,f = distance from barrier ends to x(i,jmax) (m) c aa, ca = angles opposite lenths a and c. c waa = angle between barrier and wind direction c alpha = clockwise angle from d (0 deg) to a c ceta = clockwise angle from d (0 deg) to c c dmin = minimum distance from barrier to x(i,j) c w = length from barrier to x(i,j) along wind direction c tmpa = temporary (intermediate) angle calculation c c c c + + + SUBROUTINES CALLED + + + c none c c + + + FUNCTION DECLARATIONS + + + real fu, slen c c + + + DATA INITIALIZATIONS + + + c c initialize variable (fraction of open field fric. vel) do 12 i = 0, imax do 11 j = 0, jmax do 10 k = 1, 8 w0br(i,j,k) = 1.0 10 continue 11 continue 12 continue c c + + + END SPECIFICATIONS + + + c c calc. 8 cardinal wind directions relative to the simulation c region using a relative y-axis orientation of 0 deg for region. c do 20 i = 1,8 awar(i) = (i-1)*45 - amasim if (awar(i) .lt. 0) then awar(i) = awar(i) + 360 elseif (awar(i) .gt. 360) then awar(i) = awar(i) - 360 endif c c convert relative wind angles to radians awar(i) = awar(i)*pi/180 c 20 continue c c update interior nodes do 91 i = 1, imax-1 do 90 j = 1, jmax-1 c barrier sweep do 70 n = 1, nbr c find distances a, c from x(i,j) to barrier ends a = slen(ix*i, jy*j, amxbr(1,1,n), amxbr(2,1,n)) c = slen(ix*i, jy*j, amxbr(1,2,n), amxbr(2,2,n)) c find barrier length b = slen(amxbr(1,1,n), amxbr(2,1,n), amxbr(1,2,n), & amxbr(2,2,n)) c c for a grid cell on a barrier calc. fric. vel. fraction if ((a+c) .le. b) then do 50 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)) c ^^^ used only in test version c wlen(i,j,k) = 0 50 continue go to 70 endif c c find angles from x(i,j) to barrier ends tmpa = ((b*b+c*c-a*a)/(2*b*c)) if (tmpa .lt. -1) then tmpa = -1 elseif (tmpa .gt. 1) then tmpa = 1 endif aa = acos(tmpa) tmpa = ((a*a+b*b-c*c)/(2*a*b)) if (tmpa .lt. -1) then tmpa = -1 elseif (tmpa .gt. 1) then tmpa = 1 endif ca = acos(tmpa) c c find minimum distance to barrier if (aa .gt. 90*pi/180 .or. ca .gt. 90*pi/180) then dmin = min(a,c) else dmin = b*sin(aa)*sin(ca)/sin(aa+ca) endif c c is x(i,j) within 35*barrier height of barrier? if (dmin .le. 35*amzbr(n)) then c find angle alpha to side 'a' from relative 0 deg. c find angle ceta to side 'c' from relative 0 deg. c distance from x(i,j) to north boundary ref pt. d = slen(ix*i, jy*j, ix*i, jy*jmax) c distance from barrier pt2 to north boundary ref pt. f = slen(amxbr(1,2,n), amxbr(2,2,n), ix*i, jy*jmax) c distance from barrier pt1 to north boundar ref pt. e = slen(amxbr(1,1,n), amxbr(2,1,n), ix*i, jy*jmax) alpha = acos((a*a+d*d-e*e)/(2*a*d)) ceta = acos((c*c+d*d-f*f)/(2*c*d)) if (ix*i .gt. amxbr(1,1,n)) then alpha = 2*pi - alpha endif if (ix*i .gt. amxbr(1,2,n)) then ceta = 2*pi - ceta endif else go to 70 endif c c sweep relative wind directions do 60 k = 1, 8 c if grid cell downwind from barrier make tmpv = 1 tmpv = 0 if (abs(alpha - ceta) .lt. pi) then if((awar(k) .ge. min(alpha,ceta)) & .and. (awar(k) .lt. max(alpha, ceta))) then tmpv = 1 endif elseif (abs(alpha-ceta) .ge. pi) then if(awar(k) .le. min(alpha,ceta)) then tmpv = 1 endif if ((awar(k) .ge. max(alpha,ceta)) .and. & (awar(k) .le. 2*pi)) then tmpv = 1 endif endif c ^^^tmp used for test only c if (i .eq. 24 .and. j .eq. 5) then c write(*,*) 'tmp output' c write (*,*) 'a=',a, ' b=', b, ' c=',c c write(*,*) 'aa=', aa, ' ca=',ca c write(*,*) 'dmin=', dmin c write(*,*) 'alpha=', alpha, ' ceta=', ceta c write(*,*) 'tmpv=', tmpv c write(*,*) 'awar(K)=', awar(k), ' k=', k c endif c ^^^ end tmp c if (tmpv .gt. 0) then c calc. barrier to x(i,j) distance 'w' along wind vector c calc. opposite side 'w' tmpa = abs(awar(k) - alpha) if (tmpa .gt. pi) then tmpa = 2*pi - tmpa endif waa = pi- tmpa - ca c by sin rule calc w w = sin(ca)*a/sin(waa) c c 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)) c ^^^^ tmp used only in test version c wlen(i,j,k) = w else c is x(i,j) within 5H of barrier if (dmin .lt. 5*amzbr(n)) then c 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)) c ^^^ tmp used only in test version c wlen(i,j,k) = -dmin endif endif 60 continue 70 continue 90 continue 91 continue c end c c function to calc. length real function slen(x1,y1,x2,y2) real x1, y1, x2, y2 slen = abs(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))) return end c c function to calc. fu (fraction upwind fric. velocity c near the barrier) c (ranges: porosity 0 to 0.9, distance: -5*zbr to 50*zbr) real function fu (xh, zbr, pbr, xbrw) real xh, zbr, pbr, xbrw real a, b, c, d, x, xw, pb c scale distance & width by barrier height x = xh/zbr xw = xbrw/zbr c increase effective porosity with barrier width pb = pbr + 0.02*xw c calculate coef. as fn of porosity a = 0.003+0.005*pb b = 1.35*exp(-0.5*pb**0.2) c = 10*(1-0.5*pb) d = 3 - pb c calc. frac. of fric. vel. fu = 1 - exp(-a*x**2) + b*exp(-0.003*(x+c)**d) return end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++