!$Author: fredfox $ !$Date: 2002-12-04 14:59:34 $ !$Revision: 1.10.4.10 $ !$Source: /weru/cvs/weps/weps.src/erosion/erosion.for,v $ !********************************************************************** ! EROSION submodel Beta 98.12 !********************************************************************** subroutine erosion (o_unit) ! ! +++ PURPOSE +++ ! subroutine erosion is the control subroutine and calls other ! subroutines in the EROSION submodel. ! overall the EROSION submodel: ! - calculates ridge/dike spacing parallel the wind ! - determines if friction velocity exceeds threshold. ! - calculates soil loss/deposition, suspension, PM-10 on grid ! - updates soil & plant variables changed by erosion. ! - prints output to file. ! ! ****See user modifications at start of code to enable test ! output subroutines sb1out and sb2out ! ! +++ ARGUMENT DECLARATIONS +++ integer o_unit ! +++ ARGUMENT DEFINITIONS +++ ! o_unit - Unit number for output file (passed to sb1out) ! ! ! +++ PARAMETER +++ ! ! Minimum snow depth to prevent erosion real snodep, tim parameter (snodep = 20.0, tim = 3600*24) real pid180 parameter(pid180 = 3.14159/180.) ! ! + + + GLOBAL COMMON BLOCKS + + + 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' include 'wpath.inc' include 'timer.fi' include 'command.inc' include 'main/main.inc' ! ! + + + LOCAL COMMON BLOCKS + + + include 'erosion/p1erode.inc' !Needs the SURF_UPD_FLG variable include 'erosion/e2grid.inc' include 'erosion/e3grid.inc' include 'erosion/m2geo.inc' include 'erosion/s2sgeo.inc' include 'erosion/s2agg.inc' include 'erosion/s2surf.inc' ! ! +++ LOCAL VARIABLES +++ integer i,j,wustfl, icsr, outfl integer nhill, n integer day, mon, yr real wuref, rusust, rut real wzzo, wzzov, wus, wust, wusp real wr, sfd84(mnsub), time real sina, prev_dir ! ! + + + LOCAL VARIABLE DEFINITIONS + + + ! i,j - index ! wustfl - flag to update threshold friction velocity ! icsr - index of current subregion (now only one) ! outfl - flag to call sb1out & sb2out for subdaily info ! nhill - number of hills ! wuref - reference wind speed (m/s) ! rusust - ratio of friction vel. to threshold friction vel. ! wzzo - aerodynamic roughness (mm) ! wus - friction velocity (m/s) ! wust - threshold friction velocity (m/s) ! wusp - threshold friction velocity for trapping (m/s) ! wr - ratio to be exceeded fric. vel. to thresh. fric. vel. ! for erosion ! sina - sin of the acute angle between ridge and wind angles ! prev_dir - prior direction ! rut - updated ratio of rusust ! wzzo - aerodynamic roughness ! sfd84(mnsub)- soil fraction with diameter < 0.84 mm ! time - time (s) ! ! +++ SAVES +++ ! ! +++ SUBROUTINES CALLED +++ ! sbzo ! sbwus ! sbwust ! sbinit ! sbdirini ! sbigrd ! sbwind ! sberod ! sbsfdi ! ! +++ FUNCTION DECLARATIONS +++ ! +++ END SPECIFICATIONS +++ ! ! start general erosion timer call timer(TIMEROS,TIMSTART) ! code to output standalone erosion input file on specified date ! check for day of simulation for which you want a file created if( saeinp_daysim.gt.0 ) then if ((am0jd-ijday+1).eq.saeinp_daysim) then call caldat (am0jd,day,mon,yr) write(*,*) 'Stand alone erosion input file created D/M/Y',& & day,'/',mon,'/',yr,'simulation day',saeinp_daysim call saeinp ! output daily erosion stuff end if else if( saeinp_jday.gt.0 ) then if ((am0jd).eq.saeinp_jday) then call caldat (am0jd,day,mon,yr) write(*,*) 'Stand alone erosion input file created D/M/Y',& & day,'/',mon,'/',yr,'simulation day',am0jd-ijday+1 call saeinp ! output daily erosion stuff end if end if ! initialize wind direction array for subhourly values ! presently, there is only one direction for each day ! erosion stand alone may want to input wind speed ! and direction as subhourly pairs. If so, this can be disabled. do 3 i =1, ntstep awdir(i) = awadir 3 continue ! ****Important notes to users of this submodel ! ****flag to enable subdaily output for validation ! turn off, i.e. set = 0 in WEPS outfl = 0 ! ! set ratio: defined as the ratio wus/wust rusust = 0 ! !*****this check cannot be done if subdaily wind is used FAF do 20 icsr=1, nsubr ! If snow depth > 20 mm in all subregions, then no erosion if (ahzsnd(icsr) .le. snodep) then ! Have insufficient snow depth ! ! calc if daily max friction vel. exceeds threshold in any subregions ! without hill and barrier effects ! ! calc. ridge spacing parallel the wind if (aszrgh(icsr) > 5.0) then sina = abs(sin(pid180*abs(awdir(icsr) - asargo(icsr)))) sina = max(0.10, sina) sxprg(icsr) = asxrgs(icsr)/sina if (asxdks(icsr) > asxrgs(icsr)/3.) then sxprg(icsr) = amin1(sxprg(icsr), asxdks(icsr)) endif else sxprg(icsr) = 1000 endif ! ! Compute Zo (wzzo) of surface call sbzo & & (sxprg(icsr), aszrgh(icsr), aslrr(icsr), & & wzoflg, abrlai(icsr), abrsai(icsr), abzht(icsr), & & wzzo, wzzov, awzzo) ! ! Compute soil surface friction velocity (wus) call sbwus & & (anemht, awzzo, awudmx, wzzov, abrlai(icsr), & & abrsai(icsr), wus) ! ! Calculate soil clod fraction less than 0.84 mm diameter ! calc soil mass < 0.84 mm call sbsfdi & &(aslagm(1,icsr), as0ags(1,icsr), aslagn(1,icsr), aslagx(1,icsr), & & 0.84, sfd84(icsr) ) ! ! save the initial sf84 value sf84ic = sfd84(icsr) ! save the initial rock cover svrocic = asvroc(1,icsr) ! ! Compute friction velocity threshold for entrainment (wust) and ! transport friction velocity threshold (wusp) call sbwust & & (sfd84(icsr), asdagd(1,icsr), asfcr(icsr), asvroc(1,icsr), & & asflos(icsr),abffcv(icsr),wzzo, ahrwc0(12,icsr), ahrwcw(1,icsr), & & wus,sf84ic, rusust, svrocic, dmlos(1,1), & & wust, wusp, sf84mn(1,1), smaglos(1,1)) ! This is useful only for multiple subregions. ! Checks to find maximum ratio between surface friction velocity ! and friction velocity threshold among all subregion surfaces. ! We have erosion if (wus/wust .gt. 1.0) - for flat fields only rusust = amax1 (rusust, wus/wust) ! endif 20 continue ! Some placeholder code for hills ! (it adjusts the threshold ratio cutoff for erosion) ! following st. is tmp. until sbhill is available ! nhill = 0 ! ! if (nhill .eq. 0) then wr = 1.0 ! else ! wr = 0.7 ! endif ! ! If no erosion, then exit out of erosion submodel if (rusust .le. wr) then go to 100 endif ! ! sbinit calls sbsdfi to get sf< 0.01,0.1,0.84,2.0 mm ! and writes to grid, writes other var. to grid and ! zeros eros output arrays. call sbinit ! calc. sweep direction based on wind direction for sberod prev_dir = awdir(1)+ 1.0 !make different to force calculation call sbdirini( awdir(1), prev_dir ) ! set flag on to update threshold fric. vel. on grid wustfl = 1 ! set ref wind speed to daily max wuref = awudmx ! ! step thru each periodic wind speed do 40 i =1, ntstep ! check for erodible wind speed if (awu(i) .lt. 5.0) go to 40 ! calc. sweep direction based on wind direction for sberod ! only needed if one reads hourly wind directions for input ! call sbdirini( awdir(i), prev_dir ) ! rut = rusust*awu(i)/wuref if (wustfl < 1) then if(rut .gt. wr) then wustfl = 1 else go to 40 endif endif ! force calculation to 15 minute time steps: ! useful to allow enough updates of surface n = max(1,96/ntstep) ! modify time step to more or less than 15 minutes if (rut < 1.2 ) then n = n - 2 elseif (rut > 1.6) then n = n*2 endif ! insure at least 1 time step n = max(1,n) ! ! calculate the time step in seconds time = tim/(n*ntstep) ! ! start the inner loop time step j = 1 do while (j <= n) ! prepare to update rusust and wus. ! note: when rusust= <0.1, sbaglos does not calculate. rusust = 0.2 ! updates the fric. vel and threshold fric. vel on grid ! and calc. max. for rusust = wus/wust ! this subroutine calls sbzo and sbwus ! stop general timer and start sbwind timer call timer(TIMEROS,TIMSTOP) call timer(TIMSBWIND,TIMSTART) call sbwind (wustfl, awu(i), awdir(i), ntstep, i, rusust) wuref = awu(i) wr = 1 ! stop sbwind timer and start general timer call timer(TIMSBWIND,TIMSTOP) call timer(TIMEROS,TIMSTART) ! if (rusust .gt. wr) then wustfl = 1 if (outfl .eq. 1) then call sb1out (anemht, wzoflg, kbr, j, awu(i), o_unit) endif ! stop gneral timer and start sberod timer call timer(TIMEROS,TIMSTOP) call timer(TIMSBEROD,TIMSTART) call sberod (time,SURF_UPD_FLG) call timer(TIMSBEROD,TIMSTOP) call timer(TIMEROS,TIMSTART) if (outfl .eq. 1) then call sb2out (anemht, wzoflg, awu(i), kbr, o_unit) endif j = j + 1 else ! exit if all fric vel below threshold if (rusust*awudmx/wuref .le. wr) go to 50 ! end inner loop and go to next wind speed wustfl = 0 j = n + 1 endif enddo 40 continue ! call saeinp ! output daily erosion inputs for testing 50 continue ! ! Average surface conditions and update ! ! Purpose: update global variables changed by erosion at end of day ! not implemented partly because windgen does not correlate days. 100 continue ! call timer(TIMEROS,TIMSTOP) return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++