c********************************************************************** c EROSION submodel Beta 98.12 c********************************************************************** subroutine erosion ( i o_unit) c c +++ PURPOSE +++ c subroutine erosion serves as the control subroutine and calls other c subroutines in the EROSION submodel. c overall the EROSION submodel: c - determine if friction velocity exceeds threshold. c - calculates soil loss/deposition, suspension, PM-10 on grid c - updates soil & plant variables changed by erosion. c - prints output to file. c c ****See suggested user modifications at start of code c c +++ ARGUMENT DECLARATIONS +++ integer o_unit c +++ ARGUMENT DEFINITIONS +++ c o_unit - Unit number for output file (passed to sbout) c c c +++ PARAMETER +++ c c Minimum snow depth to prevent erosion real snodep, tim parameter (snodep = 20.0, tim = 3600*24) c c + + + GLOBAL COMMON BLOCKS + + + c compiler instr. - no warn of unreferenced symbols in include files *$noreference include 'p1werm.inc' include 'm1subr.inc' include 'p1const.inc' include 'b1glob.inc' c include 'b1geom.inc' include 'm1geo.inc' include 'w1wind.inc' include 's1dbh.inc' include 's1phys.inc' include 's1agg.inc' include 's1surf.inc' include 's1sgeo.inc' include 'h1db1.inc' include 'm1flag.inc' include 'm1sim.inc' include 'timer.fi' *$reference c c + + + LOCAL COMMON BLOCKS + + + *$noreference include 'erosion/e2grid.inc' include 'erosion/e3grid.inc' include 'erosion/m2geo.inc' include 'erosion/s2sgeo.inc' include 'erosion/s2agg.inc' include 'erosion/s2surf.inc' *$reference c c +++ LOCAL VARIABLES +++ integer i,j,wustfl, icsr, outfl integer nhill, n real*4 wuref, a real*4 wzzo, wzzov, wus, wust, wusp real*4 wr, sfd84(mnsub), time real*4 prev_dir c c + + + LOCAL VARIABLE DEFINITIONS + + + c i - c wustfl - c icsr - c outfl - c nhill - c wuref - c a - c wzzo - c wus - c wust - c wusp - c wr - c prev_dir - c c c +++ SAVES +++ c c +++ SUBROUTINES CALLED +++ c sbzo c sbwus c sbwust c sbinit c sbdirini c sbigrd c sbwind c sberod c c sbout c sbsfdi c c +++ FUNCTION DECLARATIONS +++ c c +++ END SPECIFICATIONS +++ c c start general erosion timer call timer(TIMEROS,TIMSTART) c initialize wind direction array for subhourly values c presently, there is only one direction for each day c erosion stand alone may want to input wind speed c and direction as subhourly pairs. If so, this can be disabled. do 3 i =1, ntstep awdir(i) = awadir 3 continue c ****Important notes to users of this submodel c ****flag to enable subdaily output for validation c turn off, i.e. set = 0 in WEPS outfl = 0 c c set a: defined as the ratio wus/wust a = 0 c c*****this check cannot be done if subdaily wind is used FAF do 20 icsr=1, nsubr c If snow depth > 20 mm in all subregions, then no erosion if (ahzsnd(icsr) .le. snodep) then c Have insufficient snow depth c c calc if daily max friction vel. exceeds threshold in any subregions c without hill and barrier effects c c ^^^ tmp out c write (*,*) 'call to sbzo ' c ^^^ tmp initialization for dike spacing asxdks(icsr) = 0.0 c Compute Zo (wzzo) of surface call sbzo i (awadir, asargo(icsr), aszrgh(icsr), asxrgs(icsr),aslrr(icsr), i wzoflg, abrlai(icsr), abrsai(icsr), abzht(icsr),asxdks(icsr), o wzzo, wzzov, awzzo) c c ^^^ tmp out c write(*,*) 'call to sbwus' c write(*,*) 'anemht awzzo awudmx ' c write(*,*) anemht, awzzo, awudmx c write(*,*) 'wzzov abrlai abrsai icsr' c write(*,*) wzzov, abrlai(icsr), abrsai(icsr),icsr c Compute soil surface friction velocity (wus) call sbwus i (anemht, awzzo, awudmx, wzzov, abrlai(icsr), i abrsai(icsr), o wus) c c ^^^ tmp out c write (*,*) 'call to sbsfdi ' c Calculate soil clod fraction less than 0.84 mm diameter c calc soil mass < 0.84 mm call sbsfdi i(aslagm(1,icsr), as0ags(1,icsr), aslagn(1,icsr), aslagx(1,icsr), i 0.84, o sfd84(icsr) ) c c c ^^^ tmp out c write(*,*) 'call to sbwust' c Compute friction velocity threshold for entrainment (wust) and c transport friction velocity threshold (wusp) call sbwust i (sfd84(icsr), asdagd(1,icsr), asfcr(icsr), asvroc(1,icsr), i asflos(icsr),abffcv(icsr),wzzo, ahrwc0(12,icsr), ahrwcw(1,icsr), o wust, wusp) c This is useful only for multiple subregions. c Checks to find maximum ratio between surface friction velocity c and friction velocity threshold among all subregion surfaces. c We have erosion if (wus/wust .gt. 1.0) - for flat fields only a = amax1 (a, wus/wust) c endif 20 continue c^^^ tmp output c write (*,*) 'erosion out after 1st loop using max ws' c write (*,*) 'wzzo, awzzo, wzoflg ' c write (*,*) wzzo, awzzo, wzoflg c write (*,*) 'wus, wust, wusp ' c write (*,*) wus, wust, wusp c write (*,*) c ^^^ end temp output c Some of Hagen's placeholder code for hills c (it adjusts the threshold ratio cutoff for erosion) c following st. is tmp. until sbhill is available nhill = 0 c if (nhill .eq. 0) then wr = 1.0 else wr = 0.7 endif c c If no erosion, then bail out of erosion submodel if (a .le. wr) then C write (*,*) 'Bailing out because there is no erosion' go to 100 endif c c update grid surface conditions call sbinit c calc. sweep direction based on wind direction for sberod prev_dir = awdir(1)+1.0 !make sure different to force calculation call sbdirini( awdir(1), prev_dir ) c set flag on to update threshold fric. vel. on grid wustfl = 1 c set ref wind speed to daily max wuref = awudmx c c step thru each periodic wind speed do 40 i =1, ntstep c check for erodible wind speed if (awu(i) .lt. 5.0) go to 40 c calc. sweep direction based on wind direction for sberod call sbdirini( awdir(i), prev_dir ) c if ((wustfl .eq. 1) .or. i (wustfl .eq. 0 .and. a*awu(i)/wuref .gt. wr)) then c force calculation to 15 minute time steps: c useful to allow enough updates of surface n = max(1,96/ntstep) c c calculate the time step in seconds time = tim/(n*ntstep) c c start the inner loop time step do 30 j = 1, n c prepare to update a and wus. a = 0. c updates the fric. vel and threshold fric. vel on grid c and calc. max. for a = wus/wust c this subroutine calls sbzo and sbwus c^^^ tmp out c write (*,*) 'call to sbwind' c stop general timer and start sbwind timer call timer(TIMEROS,TIMSTOP) call timer(TIMSBWIND,TIMSTART) call sbwind i (wustfl, awu(i), awdir(i), ntstep, i, o a) wuref = awu(i) wr = 1 c stop sbwind timer and start general timer call timer(TIMSBWIND,TIMSTOP) call timer(TIMEROS,TIMSTART) c^^^ tmpout c write (*,*) 'EROSION OUT: j, wustfl', j, wustfl c write (*,*) ' a, wuref, awu(i)', a, wuref, awu(i) c if (outfl .eq. 1) then c call sb1out c i (anemht, wzoflg, kbr, j, awu(i), o_unit) c c endif c^^^ end tmpout if (a .gt. wr) then wustfl = 1 if (outfl .eq. 1) then call sb1out i (anemht, wzoflg, kbr, j, awu(i), o_unit) endif c^^^ tmp c write(*,*) 'Call to sberod' c write(*,*) c stop gneral timer and start sberod timer call timer(TIMEROS,TIMSTOP) call timer(TIMSBEROD,TIMSTART) call sberod i(time ) call timer(TIMSBEROD,TIMSTOP) call timer(TIMEROS,TIMSTART) if (outfl .eq. 1) then call sb2out i (anemht, wzoflg, awu(i), kbr, o_unit) endif else c exit if all fric vel below threshold if (a*awudmx/wuref .le. wr) go to 50 wustfl = 0 endif 30 continue c endif 40 continue ! call saeinp ! output daily erosion stuff 50 continue c c Average surface conditions and update c c Purpose: update global variables changed by erosion at end of day c not used partly because windgen does not correlate days. c c Average the surface varibles across each subregion c currently only a single subregion c do 80 icsr = 1,1 c set global values to zero, for use in suming c aszcr(icsr) = 0.0 c asfcr(icsr) = 0.0 c asmlos(icsr) = 0.0 c aszrgh(icsr) = 0.0 c aslrr(icsr) = 0.0 c sfd84(icsr) = 0.0 c c sum across subregions c do 70 j = 1, jmax-1 c do 60 i = 1, imax-1 c aszcr(icsr) = aszcr(icsr) + szcr(i,j) c asfcr(icsr) = asfcr(icsr) + sfcr(i,j) c asmlos(icsr) = asmlos(icsr) + sflos(i,j) c aszrgh(icsr) = aszrgh(icsr) + szrgh(i,j) c aslrr(icsr) = aslrr(icsr) + slrr(i,j) c sfd84(icsr) = sfd84(icsr) + sf84(i,j) c 60 continue c 70 continue c c find averages from summed values c aszcr(icsr) = aszcr(icsr)/ ((jmax-1)*(imax-1)) c asfcr(icsr) = asfcr(icsr)/ ((jmax-1)*(imax-1)) c asmlos(icsr) = asmlos(icsr)/((jmax-1)*(imax-1)) c aszrgh(icsr) = aszrgh(icsr)/((jmax-1)*(imax-1)) c aslrr(icsr) = aslrr(icsr)/ ((jmax-1)*(imax-1)) c sfd84(icsr) = sfd84(icsr)/ ((jmax-1)*(imax-1)) c c update asd c aslagm(1,icsr) = exp(3.78 - 7.21*sfd84(icsr)) c as0ags(1,icsr) = 1/(0.0203 + 0.00193*aslagm(1,icsr) c & + 0.074/aslagm(1,icsr)**0.5) c c 80 continue 100 continue c call timer(TIMEROS,TIMSTOP) return end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++