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 '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? C This isn't sufficient if the lai and sai (abrcd) are zero - LEW C if (abzht(icsr) .gt. 0.0 ) then if ( (abzht(icsr) .gt. 0.0) .and. & (abrsai(icsr) .gt. 0.0) .and. & (abrlai(icsr) .gt. 0.0) .and. & (abrcd(icsr) .gt. 0.0) ) then print *, 'abzht,abrsai,abrlai,abrcd: ', & abzht(icsr),abrsai(icsr),abrlai(icsr),abrcd c Compute aerodynamic roughness of canopy (wzzov) call sbdzov i (abrcd(icsr), abzht(icsr),wzzo, o wzzov) c Compute friction velocity above (wusv) and below canopy (wus) call sbwusv i (abrcd(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 ? C This isn't sufficient if the lai and sai (abrcd) are zero - LEW C if (abzht(icsr) .gt. 0.0 ) then if ( (abzht(icsr) .gt. 0.0) .and. & (abrsai(icsr) .gt. 0.0) .and. & (abrlai(icsr) .gt. 0.0) .and. & (abrcd(icsr) .gt. 0.0) ) then print *, 'abzht,abrsai,abrlai,abrcd: ', & abzht(icsr),abrsai(icsr),abrlai(icsr),abrcd call sbdzov i (abrcd(icsr), abzht(icsr), wzzo, o wzzov) call sbwusv i (abrcd(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++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++