c********************************************************************** c EROSION submodel Beta 98.12 c********************************************************************** subroutine erosion ( i o_unit) c c +++ PURPOSE +++ c subroutine sberos 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 parameter (snodep = 20.0) 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' 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' *$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 real*4 wuref, a real*4 wzzo, wzzov, wusv, wus, wust, wusp real*4 wr, aicsr, sfd84(mnsub) c c + + + LOCAL VARIABLE DEFINITIONS + + + c i = c j = c wustfl= c icsr = c outfl = c nhill = c wuref = c a = c wzzo = c wzzov = c wusv = c wus = c wust = c wusp = c wr = c aicsr = c c +++ SAVES +++ c c +++ SUBROUTINES CALLED +++ c sbdzo c sbdzov c sbwus c sbwusv c sbwust c sbinit c sbigrd c sbwind c sberod c c sbout c sbsfdi c +++ FUNCTION DECLARATIONS +++ c +++ END SPECIFICATIONS +++ 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 = 1 c ****set erosion soil loss output arrays to zero. c if user wants to accumulate beyond one day comment out call call sbigrd c set a: defined as the ratio wus/wust a = 0 c If snow depth > 20 mm in all subregions, then no erosion do 20 icsr=1, nsubr if (ahzsnd(icsr) .le. snodep) then c Have insufficient snow depth c calc if daily max friction vel. exceeds threshold in any subregions c without hill and barrier effects c Compute Zo (wzzo) of bare soil surface (below canopy) call sbdzo i (awadir, asargo(icsr), aszrgh(icsr), asxrgs(icsr), aslrr(icsr), o wzzo) c Check if a canopy is present? if (abzht(icsr) .gt. 0 ) then c Compute aerodynamic roughness of canopy (wzzov) call sbdzov i (abrlai(icsr), abrsai(icsr), abzht(icsr),wzzo, o wzzov) c Compute friction velocity above (wusv) and below canopy (wus) call sbwusv i (abrlai(icsr), abrsai(icsr), anemht, awzzo, awudmx, wzzov, o wusv, wus) else c Compute soil surface friction velocity (wus) call sbwus i (anemht, awzzo, awudmx, wzzo, o wus) endif c Compute friction velocity threshold for entrainment (wust) and c transport friction velocity threshold (wusp) call sbwust i (aslagm(1,icsr),as0ags(1,icsr),aslagn(1,icsr),aslagx(1,icsr), i asdagd(1,icsr), asfcr(icsr), asvroc(1,icsr), asflos(icsr), i 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 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 if (nhill .eq. 0) then wr = 1.0 else wr = 0.7 endif 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 c (may later want to modify sbinit to preserve c surface gradients on grid caused by erosion LH) call sbinit c 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 calc. 'a' (wus/wust) for current wind speed if ((a*awu(i)/wuref) .le. wr) go to 40 c c prepare to update a and wus. a = 0. c assign friction vel (wus) and threshold (wust) to each grid point c modified for barriers and hills c do 30 icsr = 1, nsubr if (wustfl .eq. 1) then call sbdzo i (awadir, asargo(icsr), aszrgh(icsr), asxrgs(icsr), aslrr(icsr), o wzzo ) c c canopy present ? if (abzht (icsr) .gt. 0) then call sbdzov i (abrlai(icsr), abrsai(icsr), abzht(icsr), wzzo, o wzzov) call sbwusv i (abrlai(icsr), abrsai(icsr), anemht, awzzo, awu(i), wzzov, o wusv, wus) c else call sbwus i ( anemht, awzzo, awu(i), wzzo, o wus) endif call sbwust i (aslagm(1,icsr), as0ags(1,icsr), aslagn(1,icsr), aslagx(1,icsr), i asdagd(1,icsr), asfcr(icsr), asvroc(1,icsr), asflos(icsr), i abffcv(icsr), wzzo, ahrwc0(12,icsr), ahrwcw(1,icsr), o wust, wusp) else c update only the fric. w/o updating sfc. using preceeding code wus = wus*awu(i)/wuref endif c update the grid fric. vel. and if wustfl=1, threshold fric. vel c also compute aicrs (wus/wust) for each subregion call sbwind i (wus, wust, wusp, wustfl, o aicsr) a = amax1(a, aicsr) 30 continue c turn off update wustfl = 0 c set ref. wind speed to latest awu wuref = awu(i) wr=1 c sim. region. < theshold at max wind speed? true= exit submodel if (a*awudmx/awu(i) .lt. 1. ) goto 100 c sim. region < threshold at current wind speed? 35 if (a .le. 1.) then go to 40 else c^^^ tmp C write(*,*) 'Call sberod' call sberod i (wustfl) endif c test for subdaily output of fric. vel. etc. if (outfl .eq. 1) then call sbout i (awzzo, wzzo, awu(i), wus, wust, wusp, kbr, o_unit) endif 40 continue c Average surface conditions and update c Purpose: update global variables changed by erosion at end of day c Average the surface varibles across each subregion c currently only a single subregion do 80 icsr = 1,1 c set global values to zero, for use in suming aszcr(icsr) = 0.0 asfcr(icsr) = 0.0 asmlos(icsr) = 0.0 aszrgh(icsr) = 0.0 aslrr(icsr) = 0.0 sfd84(icsr) = 0.0 c sum across subregions do 70 j = 1, jmax-1 do 60 i = 1, imax-1 aszcr(icsr) = aszcr(icsr) + szcr(i,j) asfcr(icsr) = asfcr(icsr) + sfcr(i,j) asmlos(icsr) = asmlos(icsr) + sflos(i,j) aszrgh(icsr) = aszrgh(icsr) + szrgh(i,j) aslrr(icsr) = aslrr(icsr) + slrr(i,j) sfd84(icsr) = sfd84(icsr) + sf84(i,j) 60 continue 70 continue c find averages from summed values aszcr(icsr) = aszcr(icsr)/ ((jmax-1)*(imax-1)) asfcr(icsr) = asfcr(icsr)/ ((jmax-1)*(imax-1)) asmlos(icsr) = asmlos(icsr)/((jmax-1)*(imax-1)) aszrgh(icsr) = aszrgh(icsr)/((jmax-1)*(imax-1)) aslrr(icsr) = aslrr(icsr)/ ((jmax-1)*(imax-1)) sfd84(icsr) = sfd84(icsr)/ ((jmax-1)*(imax-1)) c update asd aslagm(1,icsr) = exp(3.78 - 7.21*sfd84(icsr)) as0ags(1,icsr) = 1/(0.0203 + 0.00193*aslagm(1,icsr) & + 0.074/aslagm(1,icsr)**0.5) 80 continue 100 continue return end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++