!$Author: hagen $ !$Date: 2007-03-13 19:21:05 $ !$Revision: 1.44.2.4 $ !$Source: /weru/cvs/weps/weps.src/erosion/erosion.for,v $ subroutine erosion ! +++ 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 variables changed by erosion. ! ****See user modifications at start of code to enable test ! output subroutines sb1out and sb2out ! +++ ARGUMENT DECLARATIONS +++ ! +++ ARGUMENT DEFINITIONS +++ ! +++ PARAMETER +++ real SNODEP !Minimum snow depth to prevent erosion parameter (SNODEP = 20.0) !No erosion when snow depth >= 20mm real PID180 parameter(PID180 = 3.14159/180.) ! + + + GLOBAL COMMON BLOCKS + + + include 'p1werm.inc' include 'c1gen.inc' include 'm1subr.inc' include 'p1const.inc' include 'b1glob.inc' include 'c1glob.inc' include 'd1glob.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 'file.fi' 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, hidx real wuref, rusust, rut, rusust_preros(ntstep) real wzzo, wzzov, wus, wust, wusp, brcd, wusto real wr, sfd84(mnsub), time real sina, prev_dir real hr, sub_ntstep real wuse, wuste, enge ! + + + 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 ! hidx - hour index (for referencing surface water content) ! wuref - reference wind speed (m/s) ! rusust - ratio of friction vel. to threshold friction vel. ! rut - ratio for this timestep ! rusust_preros - ratio of friction vel. to threshold friction vel. ! by timestep (max all subregions) before erosion starts ! 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) ! wuse - est. of max. wus on grid at ntstep ! wuste - est. of min wust on grid at ntstep ! enge - est. of relative max erosive engergy at ntstep ! +++ SUBROUTINES CALLED +++ ! sbzo ! sbwus ! sbwust ! sbinit ! sbdirini ! sbigrd ! sbwind ! sberod ! sbsfdi ! +++ FUNCTION DECLARATIONS +++ ! +++ END SPECIFICATIONS +++ ! start general erosion timer call timer(TIMEROS,TIMSTART) hr = 0.0 sub_ntstep = 0.0 ! 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 ! We need to ensure that sbemit only gets called once here to write the header ! multiple days in WEPS will mess this up. ! if (btest(am0efl,2)) then ! write(0,*) 'i is:', i, 'hr is:',hr,'should be header only here' ! call sbemit (luo_emit, awu(i), hr) !Should only write the header one time ! endif ! 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 i =1, ntstep awdir(i) = awadir rusust_preros(i) = 0.0 end do ! ****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, adrlaitot(icsr), adrsaitot(icsr), abzht(icsr), & & acrlai(icsr), acrsai(icsr), aczht(icsr), & & acxrow(icsr), ac0rg(icsr), & & wzzo, wzzov, awzzo, brcd) ! 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) sf84ic = min (0.9999, max(sf84ic,0.0001)) ! edit ljh 1-23-05 do i=1, ntstep ! find hour index (1-24) hidx = int(i*23.75/ntstep) + 1 ! Compute soil surface friction velocity (wus) call sbwus( anemht, awzzo, awu(i), wzzov, brcd, wus ) ! Compute friction velocity threshold for entrainment (wust) and ! transport friction velocity threshold (wusp) dmlos(1,1) = 0.0 smaglos(1,1) = 0.0 smaglosmx(1,1)= 0.0 call sbwust( sfd84(icsr), asdagd(1,icsr), asfcr(icsr), & & asvroc(1,icsr), asflos(icsr),abffcv(icsr), wzzo, & & ahrwc0(hidx,icsr), ahrwcw(1,icsr), wus, sf84ic, & & rusust, asvroc(1,icsr), dmlos(1,1), wust, wusp, & & wusto, sf84mn(1,1), smaglos(1,1),smaglosmx(1,1)) ! Checks to find maximum ratio between surface friction velocity ! and friction velocity threshold among all time steps and subregion surfaces. ! We have erosion if (wus/wust .gt. 1.0) - for flat fields only rusust_preros(i) = max( rusust_preros(i), wus/wust ) rusust = max( rusust, rusust_preros(i) ) end do 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 41 i =1, ntstep hr = i*24/ntstep - 1.0 !change LH 3-12-07 ! check for erodible wind speed if (awu(i) .lt. 5.0) then !change LH 3-12-07 ! hr = hr + (24.0/ntstep) !No sub_ntstep's (no erosion calculated for this 'ntstep' go to 41 endif ! 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_preros(i) This change sabotaged code logic rut = rusust*awu(i)/wuref !if (wustfl < 1) then ! no erosion aka surface updating has occurred yet !if( rusust_preros(i) .gt. wr ) then ! erosion will occur, updated surface requires full grid calculation ! wustfl = 1 If (rut .le. wr) then !change LH 3-12-07 ! hr = hr + (24.0/ntstep) !No sub_ntstep's (no erosion calculated for this 'ntstep' go to 41 ! no sbemit since we do no erosion endif !endif ! The following code determines how often the surface updating code ! is executed. The internal loop interval appears to be determined based ! upon the value of "ntstep". "ntstep" determines the base updating interval ! value used within WEPS and is currently defaulting to 24 (ie, updating ! once an hour). The internal loop then can be set to update the surface ! on a smaller interval, depending upon the value of "ntstep" and some ! other variables (probably to reduce runtime). ! "ntstep" is used in the standalone erosion submodel to determine the ! default reporting and updating interval as well as the number of wind ! speed values for the day when a daily Weibull wind speed distribution ! is provided as input. Currently this value can have a maximum of 96 ! in the code). Thus, we cannot just change the value of "ntstep" as ! an input to the standalone erosion submodel to modify the frequency ! of the surface updating without expanding array dimensions throughout ! entire WEPS code and changing the wind speed input information. ! Therefore, a commandline option has been added to the standalone erosion ! submodel code to allow us to completely overide the current code ! that determines the frequency of surface updating, if desired. If the ! value of the "surface update interval" variable is set to a value ! greater than 59 and the product of "ntstep" and "erod_interval" is ! evenly divisble into the number of seconds in a day, the surface updating ! will occur at the computed interval of: ! update_interval = SEC_PER_DAY/(ntstep * n) ! (seconds) ! if update interval is 15 minutes (900 seconds): ! ntstep = 24, then n = 4 ! ntstep = 96, then n = 1 ! if update interval is 1 minute: ! ntstep = 24, then n = 60 ! ntstep = 96, then n = 15 ! NOTE: "n" is the number of steps within a single "ntstep" to obtain ! the desired update interval. if (erod_interval > 0) then !overide the default update interval ! check for even divisibility done in tsterode main program n = SEC_PER_DAY/(erod_interval*ntstep) ! write(6,*) 'erod_interval and n',erod_interval, n else !default surface updating behavior ! 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 (awu(i) .lt. 15.0) then ! if (rut < 1.1 ) then ! n = n - 2 ! elseif (rut > 1.4) then ! n = n*2 ! endif ! elseif (rut >1.4) then ! n = n*4 ! else ! n = n*2 ! endif ! ! est. n as fn. of approx max. erosive energy 9-6-06 LH wuse = 0.06*awu(i) wuste = wuse/rut enge = wuse*wuse*(wuse-wuste) n = nint((0.5 + 4.6*enge)*n) ! insure at least 1 time step n = max(1,n) endif ! calculate the time step in seconds time = SEC_PER_DAY/(n*ntstep) sub_ntstep = (24.0/ntstep)/n !fraction of hr == sub_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 ! stop general timer and start sbwind timer call timer(TIMEROS,TIMSTOP) call timer(TIMSBWIND,TIMSTART) ! updates the fric. vel and threshold fric. vel on grid ! and calc. max. for rusust = wus/wust ! this subroutine calls sbzo and sbwus 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 ! erosion will occur this time step ! wustfl = 1 if (btest(am0efl,3)) then call sb1out (j, n, hr, awu(i), awdir(i), luo_sgrd) 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) hr = hr + sub_ntstep !Compute end-of-period time in fraction of hours if (btest(am0efl,3)) then call sb2out (j, n, hr, awu(i), awdir(i), luo_sgrd) endif j = j + 1 else ! end inner loop and go to next wind speed ! wustfl = 0 j = n + 1 ! hr = hr + sub_ntstep !Compute end-of-period time in fraction of hours ! We do this because we jump out of the inside loop (see j=n+1 statement above) !change LH 3-12-07 !hr = hr + (24.0/ntstep) !No sub_ntstep's (no erosion calculated for this 'ntstep' endif enddo if (btest(am0efl,2)) then ! write(0,*) 'i is:', i, 'hr is:', hr call sbemit (luo_emit, awu(i), hr) !Should only write data here endif 41 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