c********************************************************************** c subroutine sberod c********************************************************************** subroutine sberod i (wustfl) c To calc loss/dep of saltation/creep, susp. and PM-10 at cells c To call sbqout to calc. qo, qsso, q10o for each cell c To calc. deposition in the boundary cells of sim. region c +++ ARGUMENT DEFINITIONS +++ c ix = c jy = c i1,i2,i3 = c i4,i5,i6 = c ke = c sf1 = c sf10 = c sf84 = c sf200 = c sfa12 = c sfcla = c svroc = c +++ LOCAL VARIABLE DEFINITIONS +++ c cc = c i, j = c qi,qssi, q10i = c qo,qsso, q10o = c eg = c egss = c eg10 = c egt = c egtss = c egt10 = c qx = c qy = c + + + GLOBAL COMMON BLOCKS + + + *$noreference include 'p1werm.inc' include 's1agg.inc' include 's1surf.inc' include 's1dbh.inc' include 'b1glob.inc' include 'm1sim.inc' c + + + LOCAL COMMON BLOCKS + + + include 'erosion/m2geo.inc' include 'erosion/e2grid.inc' include 'erosion/e3grid.inc' include 'erosion/s2agg.inc' include 'erosion/s2surf.inc' include 'erosion/s2sgeo.inc' include 'erosion/e2erod.inc' include 'erosion/w2wind.inc' c *$reference c c +++ ARGUMENT DECLARATIONS +++ integer wustfl, tim real*4 pid180, time c c +++ PARAMETERS +++ parameter(pid180 = 3.14159/180., tim = 24*3600) C +++ LOCAL VARIABLES +++ integer i, j, icsr real*4 cc, qi, qssi, q10i, qo, qsso, q10o, eg, egss, eg10, c qx(0:mngdpt, 0:mngdpt), c qy(0:mngdpt, 0:mngdpt), c qssx(0:mngdpt, 0:mngdpt), qssy(0:mngdpt, 0:mngdpt), c q10x(0:mngdpt, 0:mngdpt), q10y(0:mngdpt, 0:mngdpt) c +++ END SPECIFICATIONS +++ c set initial conditions to zero do 50 j = 0, jmax do 45 i = 0, imax qx(i,j) = 0. qy(i,j) = 0. qssx(i,j) = 0. qssy(i,j) = 0. q10x(i,j) = 0. q10y(i,j) = 0. 45 continue 50 continue c set a correction term cc = (jy - ix)/ix c calc. time step time = tim/ntstep c update interior grid cells: do 110 i = i1, i2, i3 do 100 j = i4, i5, i6 c if(ke .eq. 1) then c calculate input discharge qi = sqrt(qx(i-i3,j)**2 + qy(i,j-i6)**2) qssi = sqrt(qssx(i-i3,j)**2 + qssy(i,j-i6)**2) q10i = sqrt(q10x(i-i3,j)**2 + q10y(i,j-i6)**2) c calc. output discharge icsr = csr(i,j) call sbqout i (wus(i,j), wust(i,j), wusp(i,j), sf1(i,j), sf10(i,j), sf84(i,j), i sf200(i,j), szcr(i,j), sfcr(i,j), sflos(i,j), smlos(i,j), i szrgh(i,j), sxprg(icsr), sarrc(i,j), slrr(i,j), i aseags(1,icsr), asecr(icsr), asfcla(1,icsr), asfsan(1,icsr), i asvroc(1,icsr), abrsai(icsr), abzht(icsr), abffcv(icsr), time, i ix, qi, qssi, q10i, wustfl, o qo, qsso, q10o ) c update output accumulation arrays c soil loss is negative: eg = -time*(qo - qi)/ix egss = -time*(qsso - qssi)/ix eg10 = -time*(q10o - q10i)/ix egt(i,j) = egt (i,j) + eg + egss egtss(i,j) = egtss(i,j) + egss egt10(i,j) = egt10(i,j) + eg10 c c update discharge vectors qx(i,j) = -qo*sin(awa*pid180) qy(i,j) = -(qo + (qo - qi)*cc)*cos(awa*pid180) qssx(i,j) = -qsso*sin(awa*pid180) qssy(i,j) = -(qsso + (qsso -qssi)*cc)*cos(awa*pid180) q10x(i,j) = -q10o*sin(awa*pid180) q10y(i,j) = -(q10o + (q10o - q10i)*cc)*cos(awa*pid180) c c update salt/creep, suspension & pm-10 crossing boundary c note the units are kg/m and different than interior cells and c the meaning also differs. c egt = salt/creep discharge (not total) c egtss = suspension discharge c egt10 = pm-10 discharge c c calculate scalar discharge crossing borders, ie. (qx/qo)**2*qo c C Check for divide by zero errors here (qo, qsso, and q10o) if (i .eq. i2) then if (qo .gt. 1.0e-10) then egt(i2+i3, j) = egt(i2+i3,j) & + time*abs(qx(i2, j)**2.0/qo) endif if (qsso .gt. 1.0e-10) then egtss(i2+i3,j) = egtss(i2+i3,j) & + time*abs(qssx(i2, j)**2.0/qsso) egt10(i2+i3, j) = egt10(i2+i3,j) & + time*abs(q10x(i2, j)**2.0/q10o) endif endif if (j .eq. i5) then if (qo .gt. 1.0e-10) then egt(i, i5+i6) = egt(i,i5+i6) & + time*abs(qy(i, i5)**2.0/qo) endif if (qsso .gt. 1.0e-10) then egtss(i,i5+i6) = egtss(i,i5+i6) & + time*abs(qssy(i, i5)**2.0/qsso) egt10(i,i5+i6) = egt10(i,i5+i6) & + time*abs(q10y(i, i5)**2.0/q10o) endif endif c else c calc input discharge c qi = sqrt(qx(j-i6,i)**2 + qy(j,i-i3)**2) c qssi = sqrt(qssx(j-i6,i)**2 + qssy(j,i-i3)**2) c q10i = sqrt(q10x(j-i6,i)**2 + q10y(j,i-i3)**2) c c calc. output discharge c icsr = csr(j,i) c call sbqout c i (wus(j,i),wust(j,i), wusp(j,i), sf1(j,i), sf10(j,i), sf84(j,i), c i sf200(j,i), szcr(i,j), sfcr(j,i), sflos(j,i), smlos(i,j), c i szrgh(j,i), sxprg(icsr), sarrc(j,i), slrr(i,j), c i aseags(1,icsr), asecr(icsr), asfcla(1,icsr), asfsan(1,icsr), c i asvroc(1,icsr), abrsai(icsr), abzht(icsr), abffcv(icsr), c i ix, qi, qssi, q10i, wustfl, time c o qo, qsso, q10o ) c update output accumulation arrays c soil loss is negative: c eg = -time*(qo - qi)/ix c egss = -time*(qsso - qssi)/ix c eg10 = -time*(q10o - q10i)/ix c egt(j,i) = egt (j,i) + eg + egss c egtss(j,i) = egtss(j,i) + egss c egt10(j,i) = egt10(j,i) + eg10 c c update discharge vectors c qx(j,i) = -qo*sin(awa*pid180) c qy(j,i) = -(qo + (qo - qi)*cc)*cos(awa*pid180) c qssx(j,i) = -qsso*sin(awa*pid180) c qssy(i,j) = -(qsso + (qsso -qssi)*cc)*cos(awa*pid180) c q10x(j,i) = -q10o*sin(awa*pid180) c q10y(j,i) = -(q10o + (q10o - q10i)*cc)*cos(awa*pid180) c c calc salt/creep, suspension & pm-10 crossing boundary c c if (i .eq. i2) egt(j, i+i3) = egt(j, i+i3) + c & time*abs(qy(j,i2)/jy) c c endif 100 continue c c Calc. salt/creep deposition in boundary cells c if (ke .eq. 1) then c egt(i, i5+i6) = egt(i, i5+i6) + time* abs(qy(i, i5)/jy) c else c egt(i5+i6, i) = egt(i5+i6, i) + time*abs(qx(i5, i)/ix) c endif c 110 continue c c temp code for output c 210 continue return end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++