c********************************************************************** c subroutine sberod c********************************************************************** subroutine sberod i(time o ) 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 To update the threshold friction velocity as the loose material c depletes upwind and increases downwind c c +++ ARGUMENT DECLARATIONS +++ real time c c +++ ARGUMENT DEFINITIONS +++ c time = time interval (seconds) c + + + GLOBAL COMMON BLOCKS + + + *$noreference include 'p1werm.inc' include 's1agg.inc' include 's1sgeo.inc' include 's1surf.inc' include 's1dbh.inc' include 'b1glob.inc' include 'm1sim.inc' include 'h1db1.inc' include 'timer.fi' c 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 +++ PARAMETERS +++ c C +++ LOCAL VARIABLES +++ integer i, j, icsr c* real lx,aa,bb,dd,la,lb,ld,ly c real*4 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 c wzzo, c c slagm(0:mngdpt, 0:mngdpt), s0ags(0:mngdpt, 0:mngdpt) c 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 c +++ END SPECIFICATIONS +++ c 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 c set a correction term c cc = (jy - ix)/ix c* set field length c lx = ix/(abs(sin_awa)+0.001) c if (lx .gt. max(ix,jy))then c c lx = max(ix,jy) c endif c grid length (lx): revised by LH 9-22-00 if (abs(tan_awa) .le. (ix/jy)) then la = jy lb = abs(tan_awa*jy) ld = abs(jy/cos_awa) else ld = abs(ix/sin_awa) lb = ix la = sqrt(ld*ld - lb*lb) endif lx = ld*(1.0 - 0.292893*la*lb/(ix*jy)) ly = ix*jy/lx c ^^^ tmp out c write (*,*) ' output from sberod' c write (*,*) 'la=', la, 'lb=', lb,'ld=', ld c write (*,*) 'lx=', lx,'ly=',ly,'ix=',ix,'jy=', jy c write (*,*) '-----------------------------' c 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 = (qx(i-i3,j)*jy + qy(i,j-i6)*ix)/ly qssi = (qssx(i-i3,j)*jy + qssy(i,j-i6)*ix)/ly q10i = (q10x(i-i3,j)*jy + q10y(i,j-i6)*ix)/ly c calc. output discharge icsr = csr(i,j) c^^^ tmp out c if (j .eq. 1 .and. i .eq. 1) then c write (*,*) 'out from sberod line 131' c write (*,*) 'imax jmax icsr sxprg' c write (*,*) imax, jmax, icsr, sxprg(icsr) c write (*,*) 'wus wust wusp sf1 sf10 sf84' c write (*,300) wus(i,j),wust(i,j), wusp(i,j), sf1(i,j), c & sf10(i,j), sf84(i,j) c write (*,*) c 300 format(1x, 10f8.3) c endif c^^^ end tmp out call timer(TIMSBEROD,TIMSTOP) call timer(TIMSBQOUT,TIMSTART) 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), asxrgs(icsr), sxprg(icsr), slrr(i,j), i aseags(1,icsr), asecr(icsr), asfcla(1,icsr), asfsan(1,icsr), i asfvfs(1,icsr),asvroc(1,icsr), abrsai(icsr), abzht(icsr), i abffcv(icsr), time, i lx, qi, qssi, q10i, i, j, imax, jmax, o qo, qsso, q10o ) call timer(TIMSBQOUT,TIMSTOP) call timer(TIMSBEROD,TIMSTART) c c update output accumulation arrays c soil loss is negative: eg = -time*(qo - qi)/lx egss = -time*(qsso - qssi)/lx eg10 = -time*(q10o - q10i)/lx 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 scalars aa = abs(-ix*cos_awa) bb = abs(-jy*sin_awa) dd = abs(aa)+abs(bb) c qx(i,j) = qo*ly*bb/(jy*dd) qy(i,j) = qo*ly*aa/(ix*dd) qssx(i,j) = qsso*ly*bb/(jy*dd) qssy(i,j) = qsso*ly*aa/(ix*dd) q10x(i,j) = q10o*ly*bb/(jy*dd) q10y(i,j) = q10o*ly*aa/(ix*dd) c ^^^tmp c if (i .eq. 1 .and. qy(i,j) .gt. 0.00001) then c write (*,*) 'tmp out sberod line 172' c write (*,*) 'i=',i,'j=',j,'aa=',aa, 'bb=',bb,'qy=',qy(i,j) c write (*,*) 'i2=',i2,'i3=', i3,'i5=', i5, 'i6=', i6 c endif c* c qx(i,j) = -qo*sin_awa c qy(i,j) = -(qo + (qo - qi)*cc)*cos_awa c qssx(i,j) = -qsso*sin_awa c qssy(i,j) = -(qsso + (qsso -qssi)*cc)*cos_awa c q10x(i,j) = -q10o*sin_awa c q10y(i,j) = -(q10o + (q10o - q10i)*cc)*cos_awa 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 c if (i .eq. i2) then if (qx(i,j) .gt. 1.0e-10) then egt(i2+i3, j) = egt(i2+i3,j) & + time*qx(i2, j) endif if (qssx(i,j) .gt. 1.0e-10) then egtss(i2+i3,j) = egtss(i2+i3,j) & + time*qssx(i2, j) egt10(i2+i3, j) = egt10(i2+i3,j) & + time*q10x(i2, j) endif endif if (j .eq. i5) then if (qy(i,j) .gt. 1.0e-10) then egt(i, i5+i6) = egt(i,i5+i6) & + time*qy(i,i5) endif if (qssy(i,j) .gt. 1.0e-10) then egtss(i,i5+i6) = egtss(i,i5+i6) & + time*qssy(i,i5) egt10(i,i5+i6) = egt10(i,i5+i6) & + time*q10y(i,i5) endif endif c c ^^^ tmp output c if (i .eq. 1 .and. j .eq. 1) then c write (*,*) c write (*,*) 'out at sberod 1,1 call sbwust sf84=', sf84(i,j) c c write (*,*) c i slagm(i,j), s0ags(i,j), aslagn(1,icsr), aslagx(1,icsr), c i ' slagm, s0ags, aslagn,aslagx' c write(*,*) c i asdagd(1,icsr), sfcr(i,j), smlos(i,j), sflos(i,j), c i ' asdagd, sfcr, smlos, sflos' c Write (*,*) c i abffcv(icsr), wzzo, ahrwc0(12,icsr), ahrwcw(1,icsr), c i ' abffcv, wzzo, ahrwc0, ahrwcw' c write (*,*) c i szrgh(i,j), slrr(i,j), wust(i,j), wusp(i,j), c i ' szrgh, slrr, wust, wusp' c c endif c c c if (i .eq. (imax-1)/2 .and. j .eq. 1) then c write (*,*) 'out at sberod (imax-1)/2,1 call c i sbwust sf84=', sf84((imax-1)/2,j) c write (*,*) c i slagm(i,j), s0ags(i,j), aslagn(1,icsr), aslagx(1,icsr), c i ' slagm, s0ags, aslagn,aslagx' c write(*,*) c i asdagd(1,icsr), sfcr(i,j), smlos(i,j), sflos(i,j), c i ' asdagd, sfcr, smlos, sflos' c Write (*,*) c i abffcv(icsr), wzzo, ahrwc0(12,icsr), ahrwcw(1,icsr), c i ' abffcv, wzzo, ahrwc0, ahrwcw' c write (*,*) c i szrgh(i,j), slrr(i,j),wust(i,j), wusp(i,j), c i ' szrgh, slrr, wust, wusp' c c endif c if (i .eq. imax-1 .and. j .eq.1) then c write (*,*) 'out at sberod imax-1,1 call sbwust sf84=', c i sf84(imax-1,j) c write (*,*) c i slagm(i,j), s0ags(i,j), aslagn(1,icsr), aslagx(1,icsr), c i ' slagm, s0ags, aslagn,aslagx' c write(*,*) c i asdagd(1,icsr), sfcr(i,j), smlos(i,j), sflos(i,j), c i ' asdagd, sfcr, smlos, sflos' c Write (*,*) c i abffcv(icsr), wzzo, ahrwc0(12,icsr), ahrwcw(1,icsr), c i ' abffcv, wzzo, ahrwc0, ahrwcw' c write (*,*) c i szrgh(i,J), slrr(i,j), wust(i,j), wusp(i,j), c i ' szrgh, slrr, wust, wusp' c c endif c^^^ end tmp out 100 continue c c 110 continue c c temp code for output c 210 continue return end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++