!>\class Erosion !!@author wjr !> Main entry point to erosion submodel. !! !! @param minErosWindSped Minimum erosive wind speed (m/s) to evaluate for erosion loss subroutine erosion (minErosWindSped) ! +++ 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 use p1werm_def use constants_def use file_def use m1flag_def use p1erode_def use anemometer_def use threshold_def use m1sim_def use m1subr_def use m1geo_def use gridmod use soilmod use sweep_init IMPLICIT NONE ! +++ ARGUMENT DECLARATIONS +++ real minErosWindSped !Minimum erosive wind speed (m/s) !to evaluate for erosion loss ! +++ ARGUMENT DEFINITIONS +++ ! +++ PARAMETER +++ !No erosion when snow depth >= 20mm real, parameter :: SNODEP = 20.0 !Minimum snow depth to prevent erosion ! + + + GLOBAL COMMON BLOCKS + + + include 'c1gen.inc' ! done include 'b1glob.inc' include 'c1glob.inc' include 'd1glob.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 'wpath.inc' include 'timer.inc' ! include 'command.inc' include 'main/main.inc' ! + + + LOCAL COMMON BLOCKS + + + ! +++ LOCAL VARIABLES +++ integer i,j,threFricVelUpdFlg, subRegX, outFlg ! integer nhill, n integer n integer day, mon, yr, hourX real refWindSped, fricThrsVelcRatio, tsFricThrsVelcRatio ! real fricThrsVelcRatio_preros(ntstep) real aeroRougRidg, aeroRougRand, aeroRougSurf, aeroRougCnpy real fricVelc, thrsFricVelc, thrsFricVelcTrap, bioDragCoef, fricVelcMod real minErosRatio, sfd84(maxNumSubRegs), time real sina, prev_dir real hr, sub_ntstep, hrs real thrsFricVelcMaxTs, thrsFricVelcMinTs, maxErosEngTs real fricVelcAnemLoc, fricVelcRandRoug, fricVelcRidgRoug, fricVelcBioDrag real thrsFricVelcBareSoil, thrsFricVelcSurfCov, thrsFricVelcSurfWet real thrsFricVelcSurfAgg, surfFracNoEmit ! + + + LOCAL VARIABLE DEFINITIONS + + + ! i,j - index ! threFricVelUpdFlg - flag to update threshold friction velocity ! subRegX - index of current subregion (now only one) ! outFlg - flag to call sb1out & sb2out for subdaily info ! not used, not sure why ! nhill - number of hills ! not used ! hourX - hour index (for referencing surface water content) ! refWindSped - reference wind speed (m/s) ! fricThrsVelcRatio - ratio of friction vel. to threshold friction vel. ! tsFricThrsVelcRatio - ratio for this timestep ! erosTsFricThrsVelcRatio - ratio of friction vel. to threshold friction vel. ! by timestep (max all subregions) before erosion starts ! aeroRougRidg - aerodynamic roughness of ridge (mm) ! aeroRougRand - aerodynamic roughness of random roughness (mm) ! aeroRougSurf - soil surface aerodynamic roughness (mm) ! aeroRougCnpy - above canopy aerodynamic roughness (mm) ! fricVelc - friction velocity (m/s) ! thrsFricVelc - threshold friction velocity (m/s) ! thrsFricVelcTrap - threshold friction velocity for trapping (m/s) ! minErosRatio - 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 ! tsFricThrsVelcRatio - updated ratio of rusust ! sfd84(maxNumSubRegs) - soil fraction with diameter < 0.84 mm ! time - time (s) ! thrsFricVelcMaxTs - est. of max. fricVelc on grid at ntstep ! thrsFricVelcMinTs - est. of min thrsFricVelc on grid at ntstep ! maxErosEngTs - est. of relative max erosive engergy at ntstep ! fricVelcAnemLoc - Friction velocity at the anemometer location ! fricVelcRandRoug - Friction velocity adjusted for site random roughness only ! fricVelcRidgRoug - Friction velocity adjusted for site ridge roughness only ! fricVelcBioDrag - Friction velocity adjusted for site biodrag only ! thrsFricVelcBareSoil - bare soil threshold friction velocity ! thrsFricVelcSurfCov - surface cover addition to bare soil threshold friction velocity ! thrsFricVelcSurfWet - surface wetness addition to bare soil threshold friction velocity ! thrsFricVelcSurfAgg - aggregate density addition to bare soil threshold friction velocity ! +++ SUBROUTINES CALLED +++ ! calcAeroRoug ! calcSurfThrsFricVelc ! initSubr ! calcWind ! zeroGrid ! updGridPts ! calcSoilLossGain ! calcSoilMassFrac ! +++ FUNCTION DECLARATIONS +++ real calcSoilMassFrac real calcFricVelc ! +++ END SPECIFICATIONS +++ ! call testGrid(grid, MaxGridPtsDir) ! call testGrid(grid) ! start general erosion timer call timer(TIMEROS,TIMSTART) hr = 0.0 sub_ntstep = 0.0 hrs = 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 calcEmit 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 calcEmit (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 avgWindDir(i) = windDir ! erosTsFricThrsVelcRatio(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 outFlg = 0 ! set ratio: defined as the ratio fricVelc/thrsFricVelc fricThrsVelcRatio = 0 !*****this check cannot be done if subhourly wind is used FAF do 20 subRegX=1, numSubr ! If snow depth > 20 mm in all subregions, then no erosion if (asoilSurfSnowDept(subRegX) .le. SNODEP) then ! Have insufficient snow depth ne_snowdepth(subRegX) = 0 ! calc if daily max friction vel. exceeds threshold in any ! subregions without hill and barrier effects ! calc. ridge spacing parallel the wind if (SoilRegionData(subRegX)%orgRidgHght > 5.0) then sina = abs(sin(degtorad*abs(avgWindDir(subRegX) - SoilRegionData(subRegX)%orgRidgOrient))) sina = abs(sin(degtorad*abs(avgWindDir(subRegX) - orgRidgOrient(subRegX)))) sina = max(0.10, sina) ridgSpacParaWind(subRegX) = SoilRegionData(subRegX)%orgRidgSpac/sina ! (E-61) if (SoilRegionData(subRegX)%orgDikeSpac > SoilRegionData(subRegX)%orgRidgSpac/3.) then ridgSpacParaWind(subRegX) = amin1(ridgSpacParaWind(subRegX), SoilRegionData(subRegX)%orgDikeSpac) endif else ridgSpacParaWind(subRegX) = 1000 endif ! Compute Zo (aeroRougSurf) of surface call calcAeroRoug & (ridgSpacParaWind(subRegX), SoilRegionData(subRegX)%orgRidgHght, SoilRegionData(subRegX)%orgRandRoug, & aeroAnemFlg, adrlaitot(subRegX), adrsaitot(subRegX), abioMassHght(subRegX), & acrlai(subRegX), acrsai(subRegX), aczht(subRegX), & cropRowSpac(subRegX), cropFurrFlg(subRegX), aeroRougRidg, aeroRougRand, & aeroRougSurf, aeroRougCnpy, aeroRougStatHght, bioDragCoef) ! Calculate soil clod fraction less than 0.84 mm diameter ! calc soil mass < 0.84 mm sfd84(subRegX) = calcSoilMassFrac( SoilRegionData(subRegX)%SoilLayerData(1)%aggGMD, & SoilRegionData(subRegX)%SoilLayerData(1)%aggGSD, SoilRegionData(subRegX)%SoilLayerData(1)%aggMinSiz, & SoilRegionData(subRegX)%SoilLayerData(1)%aggMaxSiz, 0.84) ! save the initial sfd84 value soilFracDiamLt84ic = sfd84(subRegX) soilFracDiamLt84ic = min (0.9999, max(soilFracDiamLt84ic,0.0001)) ! edit ljh 1-23-05 do i=1, ntstep ! find hour index (1-24) hourX = int(i*23.75/ntstep) + 1 ! (comparison) anemometer location surface friction velocity fricVelcAnemLoc = calcFricVelc( anemHght, aeroRougStatHght, avgWindSped(i), & aeroRougStatHght, 0.0 ) ! (comparison) site random roughness surface friction velocity fricVelcRandRoug = calcFricVelc( anemHght, aeroRougStatHght, avgWindSped(i), & aeroRougRand, 0.0) ! (comparison) site ridge (pattern) roughness surface friction velocity fricVelcRidgRoug = calcFricVelc( anemHght, aeroRougStatHght, avgWindSped(i), & aeroRougRidg, 0.0) ! (comparison) site biodrag surface friction velocity fricVelcBioDrag = calcFricVelc( anemHght, aeroRougStatHght, avgWindSped(i), & aeroRougStatHght, bioDragCoef) ! Compute soil surface friction velocity (wus) fricVelc = calcFricVelc( anemHght, aeroRougStatHght, avgWindSped(i), & aeroRougCnpy, bioDragCoef) ! Compute friction velocity threshold for entrainment (thrsFricVelc) and ! transport friction velocity threshold (thrsFricVelcTrap) grid(1,1)%soilMassAvalDelt = 0.0 grid(1,1)%soilMassAvalAggSurf = 0.0 grid(1,1)%soilMassAvalAggSurfmx= 0.0 call calcSurfThrsFricVelc( sfd84(subRegX), SoilRegionData(subRegX)%SoilLayerData(1)%aggDen, & & SoilRegionData(subRegX)%asoilCrstFrac, SoilRegionData(subRegX)%SoilLayerData(1)%fracRock, & & SoilRegionData(subRegX)%asoilLoosCovFrac,abioFracFlatCovr(subRegX), & & aeroRougSurf, asoilSurfWatrCont(hourX,subRegX), asoilLayrWiltPt(1,subRegX), fricVelc, soilFracDiamLt84ic, & & SoilRegionData(subRegX)%SoilLayerData(1)%fracRock, thrsFricVelc, thrsFricVelcTrap, & & fricVelcMod, grid(1,1)%soilFracDiamLt84mn, grid(1,1)%soilMassAvalAggSurf, grid(1,1)%soilMassAvalAggSurfmx, & & thrsFricVelcBareSoil, thrsFricVelcSurfCov, thrsFricVelcSurfWet, thrsFricVelcSurfAgg, surfFracNoEmit) ! Checks to find maximum ratio between surface friction velocity ! and friction velocity threshold among all time steps and subregion surfaces. ! We have erosion if (fricVelc/thrsFricVelc .gt. 1.0) - for flat fields only ! erosTsFricThrsVelcRatio(i) = max( erosTsFricThrsVelcRatio(i), fricVelc/thrsFricVelc ) ! fricThrsVelcRatio = max( fricThrsVelcRatio, erosTsFricThrsVelcRatio(i) ) if( fricVelc/thrsFricVelc .gt. fricThrsVelcRatio ) then ! set new maximum fricThrsVelcRatio = fricVelc/thrsFricVelc ! set reporting values for the new maximum (this is as close to erosion as we will get) ne_wus_anemom(subRegX) = fricVelcAnemLoc ne_wus_random(subRegX) = fricVelcRandRoug ne_wus_ridge(subRegX) = fricVelcRidgRoug ne_wus_biodrag(subRegX) = fricVelcBioDrag ne_wus(subRegX) = fricVelc ne_bare(subRegX) = thrsFricVelcBareSoil ne_flat_cov(subRegX) = thrsFricVelcSurfCov ne_surf_wet(subRegX) = thrsFricVelcSurfWet ne_ag_den(subRegX) = thrsFricVelcSurfAgg ne_wust(subRegX) = thrsFricVelc ne_sfd84(subRegX) = sfd84(subRegX) ne_asvroc(subRegX) = SoilRegionData(subRegX)%SoilLayerData(1)%fracRock ne_wzzo(subRegX) = aeroRougSurf ne_sfcv(subRegX) = surfFracNoEmit end if end do else ! snow prevented erosion in this subregion ne_snowdepth(subRegX) = 1 ne_wus_anemom(subRegX) = 0 ne_wus_random(subRegX) = 0 ne_wus_ridge(subRegX) = 0 ne_wus_biodrag(subRegX) = 0 ne_wus(subRegX) = 0 ne_bare(subRegX) = 1 ne_flat_cov(subRegX) = 1 ne_surf_wet(subRegX) = 1 ne_ag_den(subRegX) = 1 ne_wust(subRegX) = 1 ne_sfd84(subRegX) = 0 ne_asvroc(subRegX) = 0 ne_wzzo(subRegX) = 0 ne_sfcv(subRegX) = 0 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 minErosRatio = 1.0 ! else ! minErosRatio = 0.7 ! endif ! Check wind ration if (fricThrsVelcRatio .le. minErosRatio) then ! no erosion, set plot.out flag ne_erosion = 0 ! exit out of erosion submodel go to 100 endif ! entering erosion submodel, set plot.out flag ne_erosion = 1 ! initSubr 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 initSubr ! calc. sweep direction based on wind direction for calcSoilLossGain prev_dir = avgWindDir(1)+ 1.0 !make different to force calculation call calcWind( avgWindDir(1), prev_dir ) ! set flag on to update threshold fric. vel. on grid threFricVelUpdFlg = 1 ! set ref wind speed to daily max refWindSped = dailyMaxWindSped ! step thru each periodic wind speed do 41 i =1, ntstep ! check for erodible wind speed if (avgWindSped(i) .lt. minErosWindSped) then hr = hr + (24.0/ntstep) !No sub_ntstep's (no erosion calculated for this 'ntstep' go to 40 endif ! calc. sweep direction based on wind direction for calcSoilLossGain ! only needed if one reads hourly wind directions for input ! call calcWind( avgWindDir(i), prev_dir ) !tsFricThrsVelcRatio = erosTsFricThrsVelcRatio(i) This change sabotaged code logic tsFricThrsVelcRatio = fricThrsVelcRatio*avgWindSped(i)/refWindSped !if (threFricVelUpdFlg < 1) then ! no erosion aka surface updating has occurred yet !if( erosTsFricThrsVelcRatio(i) .gt. minErosRatio ) then ! erosion will occur, updated surface requires full grid calculation ! threFricVelUpdFlg = 1 if (tsFricThrsVelcRatio .le. minErosRatio) then hr = hr + (24.0/ntstep) !No sub_ntstep's (no erosion calculated for this 'ntstep' go to 40 ! skip since we do no erosion (still print calcEmit) 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 = secperday/(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 ! if update interval is 10 seconds: ! ntstep = 24, then n = 360 ! ntstep = 96, then n = 90 ! if update interval is 6 seconds: ! ntstep = 24, then n = 600 ! ntstep = 96, then n = 150 ! if update interval is 1 second: ! ntstep = 24, then n = 3600 ! ntstep = 96, then n = 900 ! 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 = secperday/(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 (avgWindSped(i) .lt. 15.0) then ! if (tsFricThrsVelcRatio < 1.1 ) then ! n = n - 2 ! elseif (tsFricThrsVelcRatio > 1.4) then ! n = n*2 ! endif ! elseif (tsFricThrsVelcRatio >1.4) then ! n = n*4 ! else ! n = n*2 ! endif ! ! est. n as fn. of approx max. erosive energy 9-6-06 LH thrsFricVelcMaxTs = 0.06*avgWindSped(i) thrsFricVelcMinTs = thrsFricVelcMaxTs/tsFricThrsVelcRatio maxErosEngTs = thrsFricVelcMaxTs*thrsFricVelcMaxTs*(thrsFricVelcMaxTs-thrsFricVelcMinTs) n = nint((0.5 + 4.6*maxErosEngTs)*n) ! insure at least 1 time step n = max(1,n) endif ! calculate the time step in seconds time = secperday/(n*ntstep) sub_ntstep = (24.0/ntstep)/n !fraction of hr == sub_ntstep ! start the inner loop time step j = 1 hrs = hr !Let hrs be sub_ntstep interval hr and keep current ntstep hr do while (j <= n) ! prepare to update fricThrsVelcRatio and fricVelc. ! note: when fricThrsVelcRatio= <0.1, updSurf does not calculate. fricThrsVelcRatio = 0.2 ! stop general timer and start updGridPts timer call timer(TIMEROS,TIMSTOP) call timer(TIMSBWIND,TIMSTART) ! updates the fric. vel and threshold fric. vel on grid ! and calc. max. for fricThrsVelcRatio = fricVelc/thrsFricVelc ! this subroutine calls calcAeroRoug and calcFricVelc ! call updGridPts (threFricVelUpdFlg, avgWindSped(i), avgWindDir(i), & call updGridPts (threFricVelUpdFlg, avgWindSped(i), & & ntstep, i, fricThrsVelcRatio) refWindSped = avgWindSped(i) minErosRatio = 1 ! stop updGridPts timer and start general timer call timer(TIMSBWIND,TIMSTOP) call timer(TIMEROS,TIMSTART) if (fricThrsVelcRatio .gt. minErosRatio) then ! erosion will occur this time step ! threFricVelUpdFlg = 1 if (btest(am0efl,3)) then call sb1out (j, n, hrs, avgWindSped(i), avgWindDir(i),& & luo_sgrd) endif ! stop gneral timer and start calcSoilLossGain timer call timer(TIMEROS,TIMSTOP) call timer(TIMSBEROD,TIMSTART) call calcSoilLossGain (time,SURF_UPD_FLG) call timer(TIMSBEROD,TIMSTOP) call timer(TIMEROS,TIMSTART) hrs = hrs + sub_ntstep !Compute end-of-period time in fraction of hours if (btest(am0efl,3)) then ! call sb2out (j, n, hrs, avgWindSped(i), avgWindDir(i),& call sb2out (j, n, hrs, luo_sgrd) endif j = j + 1 else ! print out initial state, even if we never call sberode() if (btest(am0efl,3).and.(j .eq. 1).and.(i .eq. 1)) then ! call updGridPts (threFricVelUpdFlg,avgWindSped(i),avgWindDir(i), & call updGridPts (threFricVelUpdFlg,avgWindSped(i), & & ntstep,i,fricThrsVelcRatio) refWindSped = avgWindSped(i) call sb1out (j, n, hrs, avgWindSped(i), avgWindDir(i),& & luo_sgrd) endif ! set to get out of inner loop and go to next wind speed - threFricVelUpdFlg = 0 j = n + 1 endif ! If we are ready to leave the sub_ntstep loop, update "hr" and go if (j .eq. (n+1)) then hr = hr + (24.0/ntstep) ! skip to next "nstep" hr - don't care if all sub_ntsteps don't equal one ntstep end if enddo 40 continue if (btest(am0efl,2)) then ! write(0,*) 'i is:', i, 'hr is:', hr ! Note that we use "hr" not "hrs" here so we report the end of the "ntstep" hr period call calcEmit (luo_emit, avgWindSped(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 !> Calculates aerodynamic roughness parameter, which is aeroRougSurf if there is no standing biomass or !! aeroRougCnpy if standing biomass is present. !! !! @param ridgSpacParaWind - row/dike spacing parallel the wind (mm) !! @param ridgHght - ridge height (mm) !! @param aeroAnemFlg - flag=0 - anemometer at station !! flag=1 - anemometer at field !! @param randRoug - random roughness (mm) !! @param resiLAI - residue leaf area index (total)(m2/m2) !! @param resiSAI - residue stem area index (total)(m2/m2) !! @param resiAvgHght - composite average residue height (m) !! @param cropLAI - crop leaf area index (m2/m2) !! @param cropSAI - crop stem area index (m2/m2) !! @param cropHght - crop height (m) !! @param cropRowSpac - crop row spacing (m) !! @param cropFurrFlg - flag=0 - crop planted in furrow bottom !! flag=1 - crop planted on ridge top !! @param aeroRougSurf - aerodynamic roughness of surface below canopy (mm) !! @param aeroRougRidg - aerodynamic roughness of ridge !! @param aeroRougRand - aerodynamic roughness of random roughness !! @param aeroRougCnpy - aerodynamic roughness length of canopy (mm) !! @param aeroRougStatHght - aerodynamic roughness at anemom. site (mm) !! @param bioDragCoef - biomass drag coefficient !! !! uses biodrag() !********************************************************************** ! subroutine calcAeroRoug (was sbzo) !********************************************************************** subroutine calcAeroRoug (ridgSpacParaWind, ridgHght, randRoug, & & aeroAnemFlg, resiLAI, resiSAI, resiAvgHght, & & cropLAI, cropSAI, cropHght, & & cropRowSpac, cropFurrFlg, aeroRougRidg, aeroRougRand, & & aeroRougSurf, aeroRougCnpy, aeroRougStatHght, bioDragCoef) ! ! +++ PURPOSE +++ ! ! Calc. aerodynamic roughness parm., aeroRougSurf, with no standing biomass ! aeroRougSurf is used by calcSurfThrsFricVelc ! ! Calc. aerodynamic roughness parm. as aeroRougCnpy, if standing biomass ! else let aeroRougCnpy = aeroRougSurf ! aeroRougCnpy is used by calcFricVelc ! Calc. biomass drag coef.(bioDragCoef) in function biodrag ! bioDragCoef is also used by calcFricVelc ! ! set anem aero. roughness and field roughness equal when anem. at ! the field site, ie. aeroAnemFlg = 1 ! to calculate aerodynamic roughness ! of vegetation canopy. ! References: Hagen, L.J., L. Lyles. 1988. Estimating Small Grain Equivalents of Shrub Dominated ! Rangelands for Wind Erosion Control. Trans ASAE 31(3):769-775 ! Armbrust, D., J.D. Bilbro, Jr. 1997. Relating Plant Canopy Characteristics to Soil Transport ! Capacity by Wind. Agronomy Journal 89:157-162 ! use p1erode_def ! +++ ARGUMENT DECLARATIONS +++ real ridgSpacParaWind, ridgHght, randRoug integer aeroAnemFlg real resiLAI, resiSAI, resiAvgHght real cropLAI, cropSAI, cropHght, cropRowSpac integer cropFurrFlg real aeroRougRidg, aeroRougRand, aeroRougSurf, aeroRougCnpy, aeroRougStatHght, bioDragCoef ! ! +++ ARGUMENT DEFINITIONS +++ ! ridgSpacParaWind - row/dike spacing parallel the wind (mm) ! ridgHght - ridge height (mm) ! aeroAnemFlg - flag=0 - anemometer at station ! flag=1 - anemometer at field ! randRoug - random roughness (mm) ! resiLAI - residue leaf area index (total)(m2/m2) ! resiSAI - residue stem area index (total)(m2/m2) ! resiAvgHght - composite average residue height (m) ! cropLAI - crop leaf area index (m2/m2) ! cropSAI - crop stem area index (m2/m2) ! cropHght - crop height (m) ! cropRowSpac - crop row spacing (m) ! cropFurrFlg - flag=0 - crop planted in furrow bottom ! flag=1 - crop planted on ridge top ! aeroRougSurf - aerodynamic roughness of surface below canopy (mm) ! aeroRougRidg - aerodynamic roughness of ridge ! aeroRougRand - aerodynamic roughness of random roughness ! aeroRougCnpy - aerodynamic roughness length of canopy (mm) ! aeroRougStatHght - aerodynamic roughness at anemom. site (mm) ! bioDragCoef - biomass drag coefficient ! ! +++ FUNCTIONS CALLED real biodrag ! +++ LOCAL VARIABLES +++ real ridgHghtSpacRatio, bioMassHght ! ! +++ LOCAL VARIABLE DEFINITIONS +++ ! ridgHghtSpacRatio - ratio of ridge height to parallel ridge spacing ! bioMassHght - biomass height (mm) ! ! +++ INCLUDE FILES+++ ! include 'erosion/p1erode.inc' !specify min/max aerodynamic roughness values include 'p1unconv.inc' ! +++ PARAMETERS +++ ! parameter(pid180 = 3.14159/180) ! pid180- radians per degree ! ! +++ END SPECIFICATIONS +++ !Note: in BLOCK.FOR !aeroAnemFlg should be set to 1 and anemomht changed ! if the anemomenter is at the field site to ! obtain correct values from calcAeroRoug !calc. for ridge aerodynamic roughness if (ridgHght .gt. 5.0) then ridgHghtSpacRatio = ridgHght / ridgSpacParaWind ! winds are never continually normal to ridges, so restrict ridgHghtSpacRatio. ridgHghtSpacRatio = min(0.20,ridgHghtSpacRatio) aeroRougRidg = ridgHght * 1/(-64.1+135.5*ridgHghtSpacRatio+(20.84/sqrt(ridgHghtSpacRatio))) ! (EQ-60) else aeroRougRidg = 0. endif !calculation for random aerodynamic roughness aeroRougRand = randRoug*0.3 ! (EQ-62) !set upper and lower limits on aerodynamic roughness aeroRougRand = min(WZZO_MAX, aeroRougRand) ! RR <= ~100.0mm aeroRougRand = max(aeroRougRand, WZZO_MIN) ! RR >= ~1.67mm !estimate combined ridge and random aerodynamic roughness !(later- no data sets at present) !chose the largest of the two. aeroRougSurf = max (aeroRougRidg, aeroRougRand) !calculate aerodynamic roughness of vegetation, if present ! calculate "effective" biomass drag coefficient ! new function for effective biomass drag coef. bioDragCoef = biodrag( resiLAI, resiSAI, cropLAI, cropSAI, cropFurrFlg, & & cropRowSpac, cropHght, ridgHght ) ! convert biomass height to mm bioMassHght = resiAvgHght * mtomm ! calculate roughness length of canopy ( in mm) if (bioDragCoef .gt. 0.1) then aeroRougCnpy = bioMassHght * 1/(17.27-(1.254*alog(bioDragCoef)/bioDragCoef)-(3.714/bioDragCoef)) !(E-68) else if( (bioMassHght .gt. 5.0) .and. (bioDragCoef .gt. 0.001) ) then ! aeroRougCnpy = bioMassHght*exp(alog(aeroRougSurf/bioMassHght) + (alog(0.11*bioMassHght/aeroRougSurf) & ! & * alog(bioDragCoef/0.01))/2.3) caused Simon's instability aeroRougCnpy = bioMassHght*(aeroRougSurf/bioMassHght+((0.11-aeroRougSurf/bioMassHght)/4.60517)*alog(bioDragCoef/0.001)) !(E-69) else aeroRougCnpy = 0.0 endif ! choose the maximum of canopy or surface roughness aeroRougCnpy = max(aeroRougCnpy, aeroRougSurf) ! if (aeroAnemFlg .eq. 1) then ! anemom. in field set aeroRougStatHght to aeroRougCnpy aeroRougStatHght = aeroRougCnpy endif !^tmp out ! write(*,*) 'calcAeroRoug out' ! write(*,*) 'aeroRougRidg aeroRougRand aeroRougCnpy bioDragCoef bioMassHght' ! write(*,*) aeroRougRidg, aeroRougRand, aeroRougCnpy, bioDragCoef, bioMassHght !^ end tmp return end !> Calculate threshold soil surface friction velocity !! as a function of ag size dist., aerodynamic roughness, !! crust, rock, & flat biomass cover,and soil surface wetness !! !! @param sfd84 - soil mass fraction in surface layer < 0.84 mm !! @param soilLayrAggDen - aggregate density (Mg/m^3) !! @param soilCrstFrac - fraction of crust cover. !! @param soilLayrRock - updated surface vol. rock > 2.0 mm (m^3/m^3). !! @param soilLoosCovFrac - soil fraction loose material cover on crust (m^3/m^3) !! @param bioFracFlatCovr - biomass fraction of flat cover (m^2/m^2) !! @param aeroRougSurf - aerodynamic roughness length of surface below canopy(mm) !! @param soilSurfWatrCont - soil water content on mass basis (at surface) (kg/kg). !! @param soilLayrWiltPt - soil water content on mass basis, at -1.5 MPa (kg/kg) !! @param fricVelc - Soil surface friction velocity (m/s) !! @param soilFracDiamLt84ic - surface soil fraction <0.84 mm initial condition !! @param asoilLayrRock - initial surface soil volume roc fraction !! @param soilMassAvalDelt - mobile soil mass change from erosion of aggregated !! @param surface (kg/m^2) (+)= gain, (-)= loss from surface. !! @param thrsFricVelc - friction velocity theshold for en (m/s) !! @param thrsFricVelcTrap - friction velocity threshold of tp and trans. cap.(m/s) !! @param soilFracDiamLt84mn - surface soil fraction <0.84 mm where thrsFricVelc= fricVelc of ag.sfc. !! @param soilMassAvalAggSurf - mobile soil reservoir of initial aggregated sfc.(kg/m^2) !! @param soilMassAvalAggSurfmx - max mobile soil reservoir of aggregateed sfc.(kg/m^2) !! @param thrsFricVelcBareSoil - bare soil threshold friction velocity !! @param thrsFricVelcSurfCov - surface cover addition to bare soil threshold friction velocity !! @param thrsFricVelcSurfWet - surface wetness addition to bare soil threshold friction velocity !! @param thrsFricVelcSurfAgg - aggregate density addition to bare soil threshold friction velocity !! @param surfFracNoEmit - fraction bare surface that does not emit !! @param fricVelcMod - calculate threshold frictional velocity of bare, smooth surface !********************************************************************** ! subroutine calcSurfThrsFricVelc (was sbwust) !********************************************************************** subroutine calcSurfThrsFricVelc (soilFracDiamLt84, soilLayrAggDen, soilCrstFrac, soilLayrRock, soilLoosCovFrac, & & bioFracFlatCovr, & ! & aeroRougSurf, soilSurfWatrCont, soilLayrWiltPt, fricVelc, soilFracDiamLt84ic, asoilLayrRock, soilMassAvalDelt, & & aeroRougSurf, soilSurfWatrCont, soilLayrWiltPt, fricVelc, soilFracDiamLt84ic, asoilLayrRock, & & thrsFricVelc, thrsFricVelcTrap, fricVelcMod, soilFracDiamLt84mn, soilMassAvalAggSurf, soilMassAvalAggSurfmx, & & thrsFricVelcBareSoil, thrsFricVelcSurfCov, thrsFricVelcSurfWet, thrsFricVelcSurfAgg, surfFracNoEmit) ! + + + PURPOSE + + + ! Calculate threshold soil surface friction velocity ! as a function of ag size dist., aerodynamic roughness, ! crust, rock, & flat biomass cover,and soil surface wetness ! + + + ARGUMENT DECLARATIONS + + + real soilFracDiamLt84, soilLayrAggDen, soilCrstFrac, soilLayrRock, soilLoosCovFrac, bioFracFlatCovr ! real aeroRougSurf, soilSurfWatrCont, soilLayrWiltPt, fricVelc, soilFracDiamLt84ic, asoilLayrRock, soilMassAvalDelt real aeroRougSurf, soilSurfWatrCont, soilLayrWiltPt, fricVelc, soilFracDiamLt84ic, asoilLayrRock real thrsFricVelc, thrsFricVelcTrap, soilFracDiamLt84mn, soilMassAvalAggSurf real soilMassAvalAggSurfmx, fricVelcMod real thrsFricVelcBareSoil, thrsFricVelcSurfCov, thrsFricVelcSurfWet, thrsFricVelcSurfAgg, surfFracNoEmit ! + + + ARGUMENT DEFINITIONS + + + ! soilFracDiamLt84 - soil mass fraction in surface layer < 0.84 mm ! soilLayrAggDen - aggregate density (Mg/m^3) ! soilCrstFrac - fraction of crust cover. ! soilLayrRock - updated surface vol. rock > 2.0 mm (m^3/m^3). ! soilLoosCovFrac - soil fraction loose material cover on crust (m^3/m^3) ! bioFracFlatCovr - biomass fraction of flat cover (m^2/m^2) ! aeroRougSurf - aerodynamic roughness length of surface below canopy(mm) ! soilSurfWatrCont - soil water content on mass basis (at surface) (kg/kg). ! soilLayrWiltPt - soil water content on mass basis, at -1.5 MPa (kg/kg) ! fricVelc - Soil surface friction velocity (m/s) ! soilFracDiamLt84ic - surface soil fraction <0.84 mm initial condition ! asoilLayrRock - initial surface soil volume roc fraction ! soilMassAvalDelt - mobile soil mass change from erosion of aggregated ! surface (kg/m^2) (+)= gain, (-)= loss from surface. ! thrsFricVelc - friction velocity theshold for en (m/s) ! thrsFricVelcTrap - friction velocity threshold of tp and trans. cap.(m/s) ! soilFracDiamLt84mn - surface soil fraction <0.84 mm where thrsFricVelc= fricVelc of ag.sfc. ! soilMassAvalAggSurf - mobile soil reservoir of initial aggregated sfc.(kg/m^2) ! soilMassAvalAggSurfmx - max mobile soil reservoir of aggregateed sfc.(kg/m^2) ! thrsFricVelcBareSoil - bare soil threshold friction velocity ! thrsFricVelcSurfCov - surface cover addition to bare soil threshold friction velocity ! thrsFricVelcSurfWet - surface wetness addition to bare soil threshold friction velocity ! thrsFricVelcSurfAgg - aggregate density addition to bare soil threshold friction velocity ! surfFracNoEmit - fraction bare surface that does not emit ! + + + LOCAL VARIABLES + + + real b1, b2, wet_rat ! + + + END SPECIFICATIONS + + + ! calculate threshold (fricVelcMod) of bare, smooth surface with soilFracDiamLt84ic ! for use in updSurf edit 6-27-06 LH b1 = -0.179 + 0.225*(log(1.5))**0.891 ! approx -0.078337 !(E-44) b2 = 0.3 + 0.06*0.5**1.2 ! approx 0.3261165 !(E-45) fricVelcMod =1.7-1.35*exp(-b1-b2*((1-soilFracDiamLt84ic)*(1-asoilLayrRock)+asoilLayrRock)**2) !(E-83) ! calc fraction bare surface that does not emit surfFracNoEmit = ((1 - soilCrstFrac)*(1 - soilFracDiamLt84) + soilCrstFrac & & - soilCrstFrac*soilLoosCovFrac)*(1 - soilLayrRock) + soilLayrRock ! to avoid a zero value surfFracNoEmit = surfFracNoEmit + 0.0001 ! check for total cover. ! if (surfFracNoEmit < 1.0) then ! calculate bare surface static threshold friction velocity b1 = -0.179 + 0.225*(alog(1 + aeroRougSurf))**0.891 !(E-84,E-44) b2 = 0.3 + 0.06*aeroRougSurf**1.2 !(E-85,E-45) thrsFricVelcBareSoil = 1.7 - 1.35*exp(-b1-b2*surfFracNoEmit**2) thrsFricVelcTrap = 1.7 - 1.35*exp(-b1-b2*0.4**2) ! else ! thrsFricVelcBareSoil = 1.85 ! thrsFricVelcTrap = 1.80 ! endif ! edit 07-17-01 ! calc change in threshold with flat cover if (bioFracFlatCovr .gt. 0) then thrsFricVelcSurfCov = (1 - exp(-1.2*bioFracFlatCovr))*(exp(-0.3*surfFracNoEmit)) else thrsFricVelcSurfCov = 0. endif ! calc change in threshold vel with wetness wet_rat = soilSurfWatrCont / soilLayrWiltPt ! if ( wet_rat .gt. 0.3) then ! if ( wet_rat .gt. 0.25) then ! triggers at 1/4 of the wilting pt if ( wet_rat .ge. 0.0) then ! thrsFricVelcSurfWet = 0.48 * wet_rat ! this modified curve closely matches the previous linear realationship ! in the 0.25 to 0.8 range, which is where the measured data are. ! (make sure the reference to the measured data is in the Tech Doc.) ! It is apparently not possible to measure threshold effects below ! 0.25 wetness ratio so whether the relationship should go smoothly ! through zero for a wetness ratio below 0.25 has not been determined. !!!!!! this function has a singularity and goes negative for ! values of wet_rat > 1.307674238 ! thrsFricVelcSurfWet = 1.0 / (11.906541 - 10.41204 * wet_rat**0.5) !(E-88,E-116) note commented out !!!! this function is better behaved thrsFricVelcSurfWet = 0.58 * (exp(wet_rat) - 1.0 - 0.7*wet_rat*wet_rat) else thrsFricVelcSurfWet = 0. endif ! After consultation with LH, it was decided that the adjustment ! to the friction velocity terms due to Agg. Density would be ! retained in the erosion code, but WEPS would default it to a ! value of 1.8 at this time. The standalone code would then ! be able to modify that value if desired (with 1.8 being the ! suggested default value) - Mar. 15, 2006 - LEW !correct for ag density, (use constant soilLayrAggDen=1.8, 5/28/03 LH) ! thrsFricVelcSurfAgg = -0.05275 !adjustment value if soilLayrAggDen == 1.8 thrsFricVelcSurfAgg = 0.3*(sqrt(soilLayrAggDen/2.65)-1.0) ! calc final static threshold friction velocity thrsFricVelc = thrsFricVelcBareSoil + thrsFricVelcSurfCov + thrsFricVelcSurfWet + thrsFricVelcSurfAgg thrsFricVelcTrap = thrsFricVelcTrap + thrsFricVelcSurfCov + thrsFricVelcSurfWet + thrsFricVelcSurfAgg ! calculate: soilMassAvalAggSurfmx; update: soilMassAvalAggSurf, soilFracDiamLt84mn call updSurf (fricVelc, thrsFricVelc, fricVelcMod, soilFracDiamLt84ic, asoilLayrRock, & & soilMassAvalAggSurfmx, soilMassAvalAggSurf, soilFracDiamLt84mn, soilFracDiamLt84) return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !> Calculates friction velocity, given station anemometer height, surface roughness, wind speed, and aerodynamic roughness. !! !! @param anemHght - parameter, anemometer height of input wind speed (m). !! @param aeroRougStatHght - parameter, surface aerodynamic roughness at input wind speed location (mm). !! @param avgWindSped - input wind speed driving EROSION submodel (m/s). !! @param aeroRougCnpy - subregion aerodynamic roughness (mm). !! @param bioDragCoef - biomass drag coefficient \f$ WU_{*}=WUF(\frac{WZO}{WZZO})^{0.067} \f$ !! @returns fricVelc - subregion soil surface friction velocity (m/s) i.e. below canopy, if one exists. !! !! Friction velocity at the weather stations !! (E-70) \f$ WUF=\frac{0.4WU}{ln(\frac{WZ}{WZZ0})} \f$ !! !! Translate friction velocity at station to friction velocity at the field !! (E-71) \f$ WU_{*}=WUF(\frac{WZO}{WZZO})^{0.067} \f$ !! !! If biomass is present translate friction velocity at station to friction velocity at the field !! (E-72) \f$ WU_{*v}=WUF(\frac{WZO_{v}}{WZZO})^{0.067} \f$ !! !! Adjust for above to below canopy when standing biomass is present !! (E-73) \f$WU_{*}=WU_{*v}(0.86exp(\frac{-BR_{cd}}{0.0298})+0.025exp(\frac{-BR_{cd}}{0.356}))\f$ !********************************************************************** ! subroutine calcFricVelc (sbwus) !********************************************************************** function calcFricVelc (anemHght, aeroRougStatHght, avgWindSped, & & aeroRougCnpy, bioDragCoef) result (fricVelc) ! ! +++ PURPOSE +++ ! To calculate subregion, friction velocity, given station ! anemometer height, surface roughness, wind speed; and subregion ! aerodynamic roughness. ! ! if standing biomass present, then calculate friction velocity ! at surface below the canopy (fricVelc). ! ! +++ ARGUMENT DECLARATIONS +++ real anemHght, aeroRougStatHght, avgWindSped, aeroRougCnpy, bioDragCoef, fricVelc ! ! +++ ARGUMENT DEFINITIONS +++ ! ! anemHght - anemometer height of input wind speed (m). ! aeroRougStatHght - surface aerodynamic roughness at input wind speed location (mm). ! avgWindSped - input wind speed driving EROSION submodel (m/s). ! aeroRougCnpy - subregion aerodynamic roughness (mm). ! bioDragCoef - biomass drag coefficient ! fricVelc - subregion soil surface friction velocity (m/s) i.e. below canopy, if one exists. ! ! +++ LOCAL VARIABLES +++ ! ! fricVelcStat - friction velocity at the weather station ! fricVelcAbovBio - friction velocity above biomass real fricVelcStat, fricVelcAbovBio ! ! +++ END SPECIFICATIONS +++ ! note: in BLOCK.FOR aeroAnemFlg should be set to 1 and anemomht ! set to correct height if anemometer is at field site ! to obtain correct values from calcFricVelc or read as ! input data in stand-alone EROSION. ! ! Calc station (input wind speed location) friction velocity fricVelcStat = avgWindSped*0.4/alog(anemHght*1000./aeroRougStatHght) !! Eq. E-70 ! ! calc subregion friction velocity fricVelc = fricVelcStat * (aeroRougCnpy/aeroRougStatHght)**0.067 !! Eq. E-71 & E-72 ! ! if standing biomass, calculate wus below canopy if (bioDragCoef .gt. 0.0001 ) then fricVelcAbovBio = fricVelc ! ! calculate friction velocity below canopy if( bioDragCoef.gt.2.56) then !check to avoid underflow fricVelc = fricVelcAbovBio * 0.25*exp(-bioDragCoef/0.356) else fricVelc = fricVelcAbovBio*(0.86*exp(-bioDragCoef/0.0298)+0.25*exp(-bioDragCoef/0.356)) !! Eq. E-73 endif fricVelc = min(fricVelc,fricVelcAbovBio) endif ! return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !> to update aeroRougSurf at each grid point; !! To update soil friction velocity on each grid point !! and modify it for barriers and hills; !! To initialize en thresh. and cp thresh. fr. velocites on grid; !! To calculate max ratios of friction velocity to threshold !! friction velocity !! !! @param intstep - current index of ntstep thru time !! @param ntstep - max. no. of time steps in day !! @param fricThrsVelcRatio - max ratio of friction velocity to thresh. friction vel. ! !! @param wind_dir - direction of wind (degrees from north) !! @param threFricVelUpdFlg - dummy !! @param avgWindSped - dummy !*********************************************************************** !* subroutine updGridPts (was sbwind) !*********************************************************************** ! subroutine updGridPts (threFricVelUpdFlg,avgWindSped, wind_dir, ntstep, intstep, fricThrsVelcRatio) subroutine updGridPts (threFricVelUpdFlg,avgWindSped, ntstep, intstep, fricThrsVelcRatio) ! ! +++ PURPOSE +++ ! to update aeroRougSurf at each grid point; ! To update soil friction velocity on each grid point ! and modify it for barriers and hills; ! To initialize en thresh. and cp thresh. fr. velocites on grid; ! To calculate max ratios of friction velocity to threshold ! friction velocity use p1werm_def use constants_def ! include 'p1const.inc' use anemometer_def use m1geo_def use m2geo_def use gridmod use soilmod IMPLICIT NONE ! ! +++ ARGUMENT DECLARATIONS +++ integer threFricVelUpdFlg,intstep, ntstep ! real avgWindSped, fricThrsVelcRatio, wind_dir real avgWindSped, fricThrsVelcRatio ! ! +++ ARGUMENT DEFINITIONS +++ ! intstep - current index of ntstep thru time ! ntstep - max. no. of time steps in day ! fricThrsVelcRatio - max ratio of friction velocity to thresh. friction vel. ! wind_dir - direction of wind (degrees from north) ! ! + + + GLOBAL COMMON BLOCKS + + + include 'c1gen.inc' include 'b1glob.inc' include 'c1glob.inc' include 'd1glob.inc' include 'h1db1.inc' include 's1agg.inc' include 's1dbh.inc' include 's1sgeo.inc' ! ! + + + GLOBAL VARIABLES USED + + + ! imax - no. grid intervals in x-direction. ! jmax - no. grid intervals in y-direction. ! subRegX - index of current subregion. ! thrsFricVelcTrap - subregion soil threshold friction vel. trans. cap. (m/s) ! thrsFricVelc - threshold fr. vel. for en. at grid points ! fricVelc - soil friction velocity at grid points corrected for ! hills and barriers (m/s). ! ! + + + LOCAL COMMON BLOCKS + + + ! include 'erosion/w2wind.inc' ! include 'erosion/e2grid.inc' include 'erosion/e3grid.inc' ! include 'erosion/s2agg.inc' ! include 'erosion/s2sgeo.inc' ! include 'erosion/s2surf.inc' ! ! +++ LOCAL VARIABLES +++ integer i,j, subRegX,k real aeroRougRidg, aeroRougRand, aeroRougSurf, aeroRougCnpy real at, rintstep, bioDragCoef real thrsFricVelcBareSoil, thrsFricVelcSurfCov, thrsFricVelcSurfWet, thrsFricVelcSurfAgg, surfFracNoEmit ! ! +++ SUBROUTINES CALLED ! calcAeroRoug ! calcFricVelc ! calcSurfThrsFricVelc real calcFricVelc ! ! + + + END SPECIFICATIONS + + + ! assign subregion index, currently only one subRegX = 1 fricThrsVelcRatio = 0.1 ! loop through grid interior to update do 40 i = 1, imax-1 do 30 j = 1, jmax-1 ! update aerodynamic roughness ! ^^^ tmp out ! write (*,*) 'in updGridPts, call to calcAeroRoug' call calcAeroRoug & & (ridgSpacParaWind(subRegX), grid(i,j)%ridgHght, grid(i,j)%randRoug, & & aeroAnemFlg, adrlaitot(subRegX), adrsaitot(subRegX), abioMassHght(subRegX), & & acrlai(subRegX), acrsai(subRegX), aczht(subRegX), & & cropRowSpac(subRegX), cropFurrFlg(subRegX), aeroRougRidg, aeroRougRand, & & aeroRougSurf, aeroRougCnpy, aeroRougStatHght, bioDragCoef) ! ^^^ tmp out ! write (*,*) 'in updGridPts, call to calcFricVelc' ! update surface (below canopy) friction velocity grid(i,j)%fricVelc = calcFricVelc (anemHght, aeroRougStatHght, avgWindSped, & & aeroRougCnpy, bioDragCoef) ! correct friction velocity for hills ! if (nhill .ne. 0 ) then ! grid(i,j)%fricVelc = grid(i,j)%fricVelc * hillRedu(i,j,kbr) ! endif ! ! correct friction velocity for barriers if (nbr .ne. 0 ) then grid(i,j)%fricVelc = grid(i,j)%fricVelc * grid(i,j)%barrRedu(kbr) endif if (threFricVelUpdFlg .eq. 1) then ! update threshold friction velocities ! calculate hour k for surface water content rintstep = intstep ! k = aint(rintstep*23.75/ntstep) + 1 k = int(rintstep*23.75/ntstep) + 1 ! I think this is what is intended aint returns real, int returns integer - LEW call calcSurfThrsFricVelc & & (grid(i,j)%soilFracDiamLt84, SoilRegionData(subRegX)%SoilLayerData(1)%aggDen, grid(i,j)%soilCrstFrac, & & grid(i,j)%soilLayrRock, & & grid(i,j)%soilLoosCovFrac, abioFracFlatCovr(subRegX),aeroRougSurf, asoilSurfWatrCont(k,subRegX), & & asoilLayrWiltPt(1,subRegX), & ! & grid(i,j)%fricVelc, soilFracDiamLt84ic, SoilRegionData(subRegX)%SoilLayerData(1)%fracRock, grid(i,j)%soilMassAvalDelt, & & grid(i,j)%fricVelc, soilFracDiamLt84ic, SoilRegionData(subRegX)%SoilLayerData(1)%fracRock, & & grid(i,j)%thrsFricVelc, grid(i,j)%thrsFricVelcTrap, fricVelcModGlob, grid(i,j)%soilFracDiamLt84mn, & & grid(i,j)%soilMassAvalAggSurf, & & grid(i,j)%soilMassAvalAggSurfmx, thrsFricVelcBareSoil, thrsFricVelcSurfCov, thrsFricVelcSurfWet, & & thrsFricVelcSurfAgg, surfFracNoEmit) ! ^^^ tmp out ! if( grid(i,j)%thrsFricVelc .le. 0.0 ) then ! write(*,*) B"updGridPts: i,j", i, j ! write(*,*) "updGridPts: fricVelc(i,j), grid(i,j)%thrsFricVelc, fricThrsVelcRatio", & ! & fricVelc(i,j), grid(i,j)%thrsFricVelc, fricThrsVelcRatio ! write(*,*) "grid(i,j)%soilFracDiamLt84 = ", grid(i,j)%soilFracDiamLt84 ! write(*,*) "SoilRegionData(subRegX)%SoilLayerData(1)%aggDen", SoilRegionData(subRegX)%SoilLayerData(1)%aggDen ! write(*,*) "grid(i,j)%soilCrstFrac", grid(i,j)%soilCrstFrac ! write(*,*) "grid(i,j)%soilLayrRock", grid(i,j)%soilLayrRock ! write(*,*) "grid(i,j)%soilLoosCovFrac ",grid(i,j)%soilLoosCovFrac ! write(*,*) "abioFracFlatCovr(subRegX)", abioFracFlatCovr(subRegX) ! write(*,*) "aeroRougSurf", aeroRougSurf ! write(*,*) "asoilSurfWatrCont(k,subRegX)", asoilSurfWatrCont(k,subRegX) ! write(*,*) "asoilLayrWiltPt(1,subRegX)", asoilLayrWiltPt(1,subRegX) ! write(*,*) "fricVelc(i,j)", fricVelc(i,j) ! write(*,*) "soilFracDiamLt84ic", soilFracDiamLt84ic ! write(*,*) "fricThrsVelcRatio", fricThrsVelcRatio ! write(*,*) "SoilRegionData(subRegX)%SoilLayerData(1)%fracRock", SoilRegionData(subRegX)%SoilLayerData(1)%fracRock ! write(*,*) "grid(i,j)%soilMassAvalDelt", grid(i,j)%soilMassAvalDelt ! write(*,*) "grid(i,j)%thrsFricVelc", grid(i,j)%thrsFricVelc ! write(*,*) "grid(i,j)%thrsFricVelcTrap", grid(i,j)%thrsFricVelcTrap ! write(*,*) "grid(i,j)%soilFracDiamLt84mn", grid(i,j)%soilFracDiamLt84mn ! write(*,*) "grid(i,j)%soilMassAvalAggSurf", grid(i,j)%soilMassAvalAggSurf ! stop ! end if endif at = grid(i,j)%fricVelc/grid(i,j)%thrsFricVelc fricThrsVelcRatio = amax1(fricThrsVelcRatio, at) 30 continue 40 continue ! write (*,*) 'at exit updGridPts fricThrsVelcRatio =', fricThrsVelcRatio ! write (*,*) ' fricVelc(3,3), thrsFricVelc(3,3)', fricVelc(3,3), thrsFricVelc(3,3) return end !++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !> To calc loss/dep of saltation/creep, susp. and PM-10 at cells !! To call sumOutput to calc. qo, qsso, q10o for each cell !! To calc. deposition in the boundary cells of sim. region !! To update the threshold friction velocity as the loose material depletes upwind and increases downwind !! !! @param time - time interval (seconds) !! @param flg - surface update flag (1=on, 0=off) !********************************************************************** ! subroutine calcSoilLossGain (was sberod) !********************************************************************** subroutine calcSoilLossGain (time,flg) ! To calc loss/dep of saltation/creep, susp. and PM-10 at cells ! To call sumOutput to calc. qo, qsso, q10o for each cell ! To calc. deposition in the boundary cells of sim. region ! To update the threshold friction velocity as the loose material ! depletes upwind and increases downwind ! use p1werm_def use m2geo_def use gridmod use soilmod IMPLICIT NONE ! +++ ARGUMENT DECLARATIONS +++ real time integer flg ! ! +++ ARGUMENT DEFINITIONS +++ ! time - time interval (seconds) ! flg - surface update flag (1=on, 0=off) ! + + + GLOBAL COMMON BLOCKS + + + include 's1agg.inc' include 's1sgeo.inc' include 's1surf.inc' include 's1dbh.inc' include 'b1glob.inc' include 'h1db1.inc' include 'timer.inc' include 'w1clig.inc' ! ! + + + LOCAL COMMON BLOCKS + + + ! include 'erosion/e2grid.inc' include 'erosion/e3grid.inc' ! include 'erosion/s2agg.inc' ! include 'erosion/s2surf.inc' ! include 'erosion/s2sgeo.inc' ! include 'erosion/e2erod.inc' ! include 'erosion/w2wind.inc' ! ! +++ PARAMETERS +++ ! ! +++ LOCAL VARIABLES +++ integer i, j, subRegX !* real lx,aa,bb,dd,la,lb,ld,ly ! real qi, qssi, q10i, qo, qsso, q10o, eg, egss, eg10 real qx(0:maxGridPtsDir, 0:maxGridPtsDir) real qy(0:maxGridPtsDir, 0:maxGridPtsDir) real qssx(0:maxGridPtsDir, 0:maxGridPtsDir), & & qssy(0:maxGridPtsDir, 0:maxGridPtsDir) real q10x(0:maxGridPtsDir, 0:maxGridPtsDir), & & q10y(0:maxGridPtsDir, 0:maxGridPtsDir) ! ! +++ END SPECIFICATIONS +++ ! ! set initial conditions to zero do 50 j = 0, jmax do 45 i = 0, imax qx(i,j) = 0. qy(i,j) = 0. qssx(i,j) = 0. qssy(i,j) = 0. q10x(i,j) = 0. q10y(i,j) = 0. 45 continue 50 continue ! ! set a correction term ! cc = (jy - ix)/ix !* set field length ! lx = ix/(abs(sin_awa)+0.001) ! if (lx .gt. max(ix,jy))then ! ! lx = max(ix,jy) ! endif ! grid length (lx): revised by LH 9-22-00 if (abs(tan_awa) .le. (ix/jy)) then la = jy lb = abs(tan_awa*jy) ld = abs(jy/cos_awa) else ld = abs(ix/sin_awa) lb = ix la = sqrt(ld*ld - lb*lb) endif lx = ld*(1.0 - 0.292893*la*lb/(ix*jy)) ly = ix*jy/lx ! ^^^ tmp out ! write (*,*) ' output from calcSoilLossGain' ! write (*,*) 'la=', la, 'lb=', lb,'ld=', ld ! write (*,*) 'lx=', lx,'ly=',ly,'ix=',ix,'jy=', jy ! write (*,*) '-----------------------------' ! ! update interior grid cells: do 110 i = i1, i2, i3 do 100 j = i4, i5, i6 ! if(ke .eq. 1) then ! calculate input discharge qi = (qx(i-i3,j)*jy + qy(i,j-i6)*ix)/ly qssi = (qssx(i-i3,j)*jy + qssy(i,j-i6)*ix)/ly q10i = (q10x(i-i3,j)*jy + q10y(i,j-i6)*ix)/ly ! calc. output discharge subRegX = grid(i,j)%subrReg !^^^ tmp out ! if (j .eq. 1 .and. i .eq. 1) then ! write (*,*) 'out from calcSoilLossGain line 131' ! write (*,*) 'imax jmax subRegX ridgSpacParaWind' ! write (*,*) imax, jmax, subRegX, ridgSpacParaWind(subRegX) ! write (*,*) 'fricVelc thrsFricVelc thrsFricVelcTrap soilFracDiamLt1 grid(i,j)%soilFracDiamLt10 soilFracDiamLt84' ! write (*,300) grid(i,j)%fricVelc,grid(i,j)%thrsFricVelc, grid(i,j)%thrsFricVelcTrap, grid(i,j)%soilFracDiamLt1, ! & grid(i,j)%soilFracDiamLt10(i,j), grid(i,j)%soilFracDiamLt84 ! write (*,*) ! 300 format(1x, 10f8.3) ! endif !^^^ end tmp out call timer(TIMSBEROD,TIMSTOP) call timer(TIMSBQOUT,TIMSTART) call sumOutput (flg, & & grid(i,j)%fricVelc, grid(i,j)%thrsFricVelc, grid(i,j)%thrsFricVelcTrap, grid(i,j)%soilFracDiamLt10, & & grid(i,j)%soilFracDiamLt84, & & grid(i,j)%soilFracDiamLt200, grid(i,j)%soilCrstThck, grid(i,j)%soilCrstFrac, grid(i,j)%soilLoosCovFrac, & & grid(i,j)%soilLoosMass, & & grid(i,j)%ridgHght, SoilRegionData(subRegX)%orgRidgSpac, ridgSpacParaWind(subRegX), grid(i,j)%randRoug, & & SoilRegionData(subRegX)%SoilLayerData(1)%fracClay, SoilRegionData(subRegX)%SoilLayerData(1)%fracSand, & & SoilRegionData(subRegX)%SoilLayerData(1)%fracVeryFineSand, grid(i,j)%soilLayrRock, abrsai(subRegX), abioMassHght(subRegX), & !edit ljh 1-22-05 & abioFracFlatCovr(subRegX), time, & & SoilRegionData(subRegX)%acoefAbraAgg, SoilRegionData(subRegX)%acoefAbraCrst,SoilRegionData(subRegX)%asoilPM10FracAbraSusp, & & SoilRegionData(subRegX)%asoilPM10FracEmitSusp, SoilRegionData(subRegX)%asoilPM10FracSaltBrkSusp, & ! & lx, qi, qssi, q10i, i, j, imax, jmax, & & lx, qi, qssi, q10i, & ! & grid(i,j)%soilMassAvalAggSurf, grid(i,j)%soilMassAvalDelt, grid(i,j)%soilFracDiamLt84mn, soilFracDiamLt84ic, & & grid(i,j)%soilMassAvalDelt, grid(i,j)%soilFracDiamLt84mn, soilFracDiamLt84ic, & & soilFracDiamLt10ic, & !edit ljh 1-22-05 & SoilRegionData(subRegX)%SoilLayerData(1)%fracRock, grid(i,j)%soilMassAvalAggSurfmx, & & qo, qsso, q10o ) call timer(TIMSBQOUT,TIMSTOP) call timer(TIMSBEROD,TIMSTART) ! ! update output accumulation arrays ! soil loss is negative: eg = -time*(qo - qi)/lx egss = -time*(qsso - qssi)/lx eg10 = -time*(q10o - q10i)/lx grid(i,j)%soilLossTot = grid(i,j)%soilLossTot + eg + egss grid(i,j)%soilLossSusp = grid(i,j)%soilLossSusp + egss grid(i,j)%soilLossPM10 = grid(i,j)%soilLossPM10 + eg10 ! !* update discharge scalars aa = abs(-ix*cos_awa) bb = abs(-jy*sin_awa) dd = abs(aa)+abs(bb) ! qx(i,j) = qo*ly*bb/(jy*dd) qy(i,j) = qo*ly*aa/(ix*dd) qssx(i,j) = qsso*ly*bb/(jy*dd) qssy(i,j) = qsso*ly*aa/(ix*dd) q10x(i,j) = q10o*ly*bb/(jy*dd) q10y(i,j) = q10o*ly*aa/(ix*dd) if (i .eq. i2) then if (qx(i,j) .gt. 1.0e-10) then grid(i2+i3,j)%soilLossTot = grid(i2+i3,j)%soilLossTot + time*qx(i2, j) endif if (qssx(i,j) .gt. 1.0e-10) then grid(i2+i3,j)%soilLossSusp = grid(i2+i3,j)%soilLossSusp + time*qssx(i2, j) grid(i2+i3,j)%soilLossPM10 = grid(i2+i3,j)%soilLossPM10 + time*q10x(i2, j) endif endif if (j .eq. i5) then if (qy(i,j) .gt. 1.0e-10) then grid(i,i5+i6)%soilLossTot = grid(i,i5+i6)%soilLossTot + time*qy(i,i5) endif if (qssy(i,j) .gt. 1.0e-10) then grid(i,i5+i6)%soilLossSusp = grid(i,i5+i6)%soilLossSusp + time*qssy(i,i5) grid(i,i5+i6)%soilLossPM10 = grid(i,i5+i6)%soilLossPM10 + time*q10y(i,i5) endif endif 100 continue 110 continue continue return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !> to calculate grid size and spacing for EROSION. !! grid size assumes outer points are outside field boundary !! a distance ix/2 !! to calculate number of grid points for EROSION. !! A max 'interior' square grid of 29X29 is assigned-no barriers !! A max 'interior' rectangular grid of 59X59 is assigned barriers !! to assign subregion index no. to each grid point. !********************************************************************** ! subroutine setupGrid (was sbgrid) !********************************************************************** subroutine setupGrid ! ! +++ PURPOSE +++ ! to calculate grid size and spacing for EROSION. ! grid size assumes outer points are outside field boundary ! a distance ix/2 ! to calculate number of grid points for EROSION. ! A max 'interior' square grid of 29X29 is assigned-no barriers ! A max 'interior' rectangular grid of 59X59 is assigned barriers ! to assign subregion index no. to each grid point. ! use p1werm_def use m1geo_def use m2geo_def use m1subr_def use gridmod IMPLICIT NONE ! +++ ARGUMENT DECLARATION +++ ! ! +++ LOCAL DEFINITIONS +++ ! imax - no. grid intervals in x-direction ! jmax - no. grid intervals in y-direction. ! ix - grid interval in x-direction (m) ! jy - grid interval in y-direction (m) ! dxmin - minimum grid interval (m) ! subrReg - current subr. index at grid point i,j ! subRegX - same as subrReg but not an array ! i,j - do loop indexes ! ! + + + GLOBAL COMMON BLOCKS + + + ! ! + + + LOCAL COMMON BLOCKS + + + ! include 'erosion/e2grid.inc' ! ! +++ LOCAL VARIABLES +++ integer ngdpt integer subRegX, i, j real dxmin, lx, ly ! ! +++ END SPECIFICATIONS +++ ! ! set min grid spacing dxmin = MIN_GRID_SP ! set max no. of grid points with no barrier ngdpt = N_G_DPT ! barriers? if (nbr .gt. 0) then ! find shortest barrier to determine dxmin do 5 i=1,nbr if (amzbr(i) > 0.0) then !Check for zero height barriers dxmin = min(dxmin, 5.0*amzbr(i)) endif 5 continue ngdpt = B_G_DPT !default to this value if a barrier exists endif ! calculate max grid intervals ! calc. lx and ly sides of field lx = amxsim (1,2)-amxsim(1,1) ly = amxsim (2,2)-amxsim(2,1) ! !^^^tmp out ! write(*,*) 'tmp out from setupGrid, line 69' ! write (*,*) 'lx=', lx, 'ly=',ly ! ^^^end tmp ! ! increase grid points on large field ! if((lx .gt. 200) .or. (ly .gt. 200)) then ! ngdpt = B_G_DPT ! endif ! ! case where lx > ly if ( lx .gt. ly)then imax = int ( lx / dxmin) imax = min(imax,ngdpt) imax = max(imax,2) ! calculate spacing for square or with barriers a rectangular grid ix = lx / (imax - 1) if (nbr .gt. 0) then jmax = int (ly / dxmin) jmax = min(jmax, ngdpt) else ! jmax = anint(ly/ix) + 1 jmax = nint(ly/ix) + 1 ! aint returns real, nint returns integer endif jmax = max(jmax,2) jy = ly/(jmax - 1) ! case where lx = ly or lx < ly else jmax = int (ly / dxmin) jmax = min(jmax,ngdpt) jmax = max(jmax,2) jy = ly / (jmax - 1) if (nbr .gt. 0) then imax = int (lx / dxmin) imax = min(imax,ngdpt) else ! imax = anint(lx/jy) + 1 imax = int(lx/jy) + 1 ! I think this is what is intended aint returns real, int returns integer - LEW endif imax = max(imax,2) ix = lx/(imax-1) endif ! ! determine subregion of each grid point ! for a single subregion now subRegX = 1 do 20 j = 0, jmax do 10 i = 0, imax ! ! for multiple subregions ! do 5 subRegX = 1, numSubr ! for multiple subregions ! if (i*ix .lt. amxsr(1,1,subRegX) .or. i*ix .gt. amxsr(1,2,subRegX)) ! & go to 10 ! if (j*jy .lt. amxsr(2,1,subRegX) .or. j*jy .gt. amxsr(2,2,subRegX)) ! & go to 10 ! 5 continue ! grid(i,j)%subrReg = subRegX 10 continue 20 continue return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! (was sbaglos) !> calc. minimum erodible fraction (soilFracDiamLt84mn) needed to stop emission !! on aggregated portion of surface at current threshold friction velocity !! that occurs when soilMassAvalAggSurf is removed. !! calc. potential mobile soil reservoir (soilMassAvalAggSurfmx) for a smooth !! surface 0.8 m/s friction velocity and soilFracDiamLt84ic that armors clod surface. !! calc. available mobile reservoir (soilMassAvalAggSurf) of soilFracDiamLt84 for current surface !! based on wus-thrsFricVelc ratio for !! !! @param fricVelc - friction velocity (m/s) !! @param thrsFricVelc - friction velocity theshold for en (m/s) !! @param fricVelcMod - threhold friction velocity = fricVelc minus flat biomass and wetness effects (m/s) !! @param soilFracDiamLt84ic - surface soil fraction <0.84 mm initial condition !! @param asoilLayrRock - surface soil volume rock (m^3/m^3) !! @param soilMassAvalAggSurfmx - max mobile soil reservoir of aggregateed sfc.(kg/m^2) !! @param soilMassAvalAggSurf - potential mobile soil reservoir of aggregated sfc.(kg/m^2) !! @param soilFracDiamLt84mn - surface soil fraction <0.84 mm where thrsFricVelc= fricVelc of ag.sfc. !! @param soilFracDiamLt84 - soil mass fraction in surface layer < 0.84 mm subroutine updSurf (fricVelc, thrsFricVelc, fricVelcMod, soilFracDiamLt84ic, asoilLayrRock, & & soilMassAvalAggSurfmx, soilMassAvalAggSurf, soilFracDiamLt84mn, soilFracDiamLt84) ! + + + PURPOSE + + + ! calc. minimum erodible fraction (soilFracDiamLt84mn) needed to stop emission ! on aggregated portion of surface at current threshold friction velocity ! that occurs when soilMassAvalAggSurf is removed. ! calc. potential mobile soil reservoir (soilMassAvalAggSurfmx) for a smooth ! surface 0.8 m/s friction velocity and soilFracDiamLt84ic that armors clod surface. ! ! calc. available mobile reservoir (soilMassAvalAggSurf) of soilFracDiamLt84 for current surface ! based on wus-thrsFricVelc ratio for ! + + + ARGUMENT DECLARATIONS + + + real fricVelc, thrsFricVelc, fricVelcMod, soilFracDiamLt84ic, asoilLayrRock real soilMassAvalAggSurfmx, soilMassAvalAggSurf, soilFracDiamLt84mn, soilFracDiamLt84 ! + + + ARGUMENT DEFINITIONS + + + ! fricVelc - friction velocity (m/s) ! thrsFricVelc - friction velocity theshold for en (m/s) ! fricVelcMod - threhold friction velocity = fricVelc minus flat biomass and wetness effects (m/s) ! soilFracDiamLt84ic - surface soil fraction <0.84 mm initial condition ! asoilLayrRock - surface soil volume rock (m^3/m^3) ! soilMassAvalAggSurfmx - max mobile soil reservoir of aggregateed sfc.(kg/m^2) ! soilMassAvalAggSurf - potential mobile soil reservoir of aggregated sfc.(kg/m^2) ! soilFracDiamLt84mn - surface soil fraction <0.84 mm where thrsFricVelc= fricVelc of ag.sfc. ! soilFracDiamLt84 - soil mass fraction in surface layer < 0.84 mm ! + + + END SPECIFICATIONS + + + ! edit LH 6-26-06 ! calc. max mobile soil at fricVelc = 0.75 m/s for bare, smooth surface ! soilFracDiamLt84 is assumed = 0 after this mass removal. soilMassAvalAggSurfmx = exp(2.708 - 7.603*((1-soilFracDiamLt84ic)*(1-asoilLayrRock) + asoilLayrRock)) ! reduce max mobile soil for roughness, cover, etc. soilMassAvalAggSurf = soilMassAvalAggSurfmx * (fricVelc - thrsFricVelc) / (0.75 - amin1(fricVelcMod,thrsFricVelc)) soilMassAvalAggSurf = max(0.0, soilMassAvalAggSurf) !soilMassAvalAggSurf = min(soilMassAvalAggSurfmx, soilMassAvalAggSurf) ! find soilFracDiamLt84 when soilMassAvalAggSurf has been removed from control volume mass(denominator). ! in sumOutput emission from cloddy surface goes to zero at soilFracDiamLt84mn if ((soilMassAvalAggSurf .eq. 0.0) .or. (soilFracDiamLt84 .eq. 0.0)) then soilFracDiamLt84mn = soilFracDiamLt84 !elseif (soilMassAvalAggSurfmx .le. soilMassAvalAggSurf) then ! soilFracDiamLt84mn = 0.05 else soilFracDiamLt84mn=(soilMassAvalAggSurfmx - soilMassAvalAggSurf)/(soilMassAvalAggSurfmx/(soilFracDiamLt84ic & & * (1.001-asoilLayrRock))) soilFracDiamLt84mn = amax1(0.0, soilFracDiamLt84mn) endif return end !> To calc the emissions for each time step of the input wind speed !! The emissions for EPA are the suspension component !! with units kg m-2 s-1. !! To write out a file in the format: !! 12 blank col, yr, mo, day, hr, soucename, emissionrate !! !********************************************************************** ! subroutine calcEmit (was sbemit) !********************************************************************** subroutine calcEmit (ounit, ws, hhr) ! To calc the emissions for each time step of the input wind speed ! The emissions for EPA are the suspension component ! with units kg m-2 s-1. ! To write out a file in the format: ! 12 blank col, yr, mo, day, hr, soucename, emissionrate ! ! Instructions & logic: ! To get ntstep period emissions output on erosion days: ! user sets am0efl = 3 in WEPS configuration screen ! subroutine openfils creates output file emit.out ! EROSION calls calcEmit to write heading in emit.out file, ! & sets am0efl to 98, then calls calcEmit ! to print (hourly) Weps emissions on erosion days. ! or ! user sets ae0efl (print flg)=4 in stand_alone input file ! EROSION opens emit.out file, calls calcEmit to write headings ! & sets ae0efl to 99, then calls calcEmit ! to print period emissions for an erosion day. ! use p1werm_def use constants_def use file_def ! include 'file.inc' use datetime_def use m2geo_def use m1sim_def use gridmod ! ! +++ ARGUMENT DECLARATIONS +++ integer ounit !Unit number for detail grid erosion real ws, hhr ! +++ ARGUMENT DEFINITIONS +++ ! ! ! + + + GLOBAL COMMON BLOCKS + + + ! include 'm1flag.inc' ! ! + + + LOCAL COMMON BLOCKS + + + ! include 'erosion/e2erod.inc' ! ! ! +++ PARAMETERS +++ ! ! +++ LOCAL VARIABLES +++ integer initflg save initflg integer prev_erosion_jday save prev_erosion_jday integer j,i integer yr, mo, da save yr, mo, da real tims, asoilLossTotp, asoilLossSuspp, asoilLossPM10p save tims, asoilLossTotp, asoilLossSuspp, asoilLossPM10p ! real hr ! save hr real asoilLossTot, asoilLossSusp, asoilLossPM10 real emittot, emitss, emit10, tt ! +++ LOCAL VARIABLE DEFINITIONS +++ ! ! ! +++ OUTPUT FORMATS +++ ! 100 format (1x,' yr mo day hr ws emission (kg m-2 s-1)') 110 format (22x,' total salt/creep susp PM10') 120 format (1x,3(i4),F7.3,F6.2, 1x,4(F11.8)) ! ! +++ END SPECIFICATIONS +++ ! ! set initial conditions if (initflg .eq. 0) then initflg = initflg + 1 prev_erosion_jday = am0jd - 1 !init to previous day ! asoilLossTotp = 0.0 ! asoilLossSuspp = 0.0 ! asoilLossPM10p = 0.0 tims = 3600*24/ntstep !seconds in each emission period call caldatw (da, mo, yr) !Set day, month and year write(0,*) 'First ntstep is: ', ntstep, tims, tims/3600 ! hr = 0 write (ounit,*) 'calcEmit output' ! write (ounit,*) 'Suspended emissions < 0.10 mm dia.' write (ounit,*) ! Print date of Run write (ounit,*) 'Date of run: ', datetimestr write (ounit,*) write (ounit,100) write (ounit,110) write (ounit,*) endif ! else !init prev erosion hr values to zero if this is new erosion day if (prev_erosion_jday .ne. am0jd) then prev_erosion_jday = am0jd asoilLossTotp = 0.0 asoilLossSuspp = 0.0 asoilLossPM10p = 0.0 endif write(0,*) 'Subsequent ntstep is: ', ntstep, tims, tims/3600 ! if (hr .ge. 24) then ! hr = tims/3600 call caldatw (da, mo, yr) ! else ! hr = hr + tims/3600 ! endif ! calculate averages of inner grid points asoilLossTot = 0.0 asoilLossSusp = 0.0 asoilLossPM10 = 0.0 ! do j=1,jmax-1 do i= 1, imax-1 asoilLossTot= asoilLossTot + grid(i,j)%soilLossTot asoilLossSusp = asoilLossSusp + grid(i,j)%soilLossSusp asoilLossPM10 = asoilLossPM10 + grid(i,j)%soilLossPM10 enddo enddo tt = (imax-1)*(jmax-1) asoilLossTot = - asoilLossTot/tt ! change signs to positive=emission asoilLossSusp = - asoilLossSusp/tt asoilLossPM10 = - asoilLossPM10/tt ! Commented out so that emission results don't get summed from hr to hr ! when erosion stops. - LEW 5/26/06 ! Hourly (or peroid) emission rate (kg m-2 s-1) ! if (asoilLossSuspp .gt. asoilLossSusp) then !main set soilLossTot arrays to zero ! asoilLossTotp = 0.0 ! asoilLossSuspp = 0.0 ! asoilLossPM10p = 0.0 ! endif ! emittot = (asoilLossTot - asoilLossTotp)/tims emitss = (asoilLossSusp - asoilLossSuspp)/tims emit10 = (asoilLossPM10 - asoilLossPM10p)/tims ! ! Save prior hour average emission asoilLossTotp = asoilLossTot asoilLossSuspp = asoilLossSusp asoilLossPM10p = asoilLossPM10 ! ! Write to emit.out file write (ounit,120) yr, mo, da, hhr, ws, & & emittot, emittot-emitss, emitss, emit10 ! ! endif return end !> Calc. wind angle on the sim. region !! Calc. sweep sequence for update of grid cells !! calc. ridge spacing parallel the wind !! !! @param wind_dir - direction of the wind in degrees from north !! @param prev_dir - previously computed direction of the wind !********************************************************************** ! subroutine calcWind - initializes grid/wind direction relationships (was sbdirini) !********************************************************************** subroutine calcWind(wind_dir, prev_dir) ! +++ purpose +++ ! Calc. wind angle on the sim. region ! Calc. sweep sequence for update of grid cells ! calc. ridge spacing parallel the wind ! use p1werm_def use m1geo_def use m2geo_def use m1subr_def use gridmod IMPLICIT NONE ! +++ ARGUMENT DECLARATIONS +++ real wind_dir real prev_dir ! ! +++ ARGUMENT DEFINITIONS +++ ! wind_dir - direction of the wind in degrees from north ! prev_dir - previously computed direction of the wind ! ! +++ PARAMETER +++ real pid180 parameter(pid180 = 3.14159/180.) ! ! + + + GLOBAL COMMON BLOCKS + + + include 's1sgeo.inc' include 'w1wind.inc' ! ! + + + LOCAL COMMON BLOCKS + + + include 'erosion/e3grid.inc' ! include 'erosion/s2sgeo.inc' ! ! + + + LOCAL VARIABLES + + + ! integer subRegX, i, j, ke ! integer i, j, ke integer ke ! real sfd1(maxNumSubRegs), sfd10(maxNumSubRegs), & ! & sfd84(maxNumSubRegs), sfd200(maxNumSubRegs) ! ! + + + LOCAL VARIABLE DEFINITIONS + + + ! subRegX = index of current subregion ! ! + + + SUBROUTINES CALLED + + + ! calcSoilMassFrac ! ! + + + END SPECIFICATION + + + ! ! check and do not calculate if done on last entry if( wind_dir.eq.prev_dir ) return prev_dir = wind_dir ! ! calc wind angle relative to the field Y-axis (+, - 45 deg. range) awa = wind_dir - amasim if (awa .lt. 0.0 ) awa = awa + 360.0 if (awa .gt. 360.0) awa = awa - 360.0 sin_awa = sin(awa*pid180) cos_awa = cos(awa*pid180) tan_awa = tan(awa*pid180) ! ! find wind quadrant relative to sim region & select sweep sequence ! If (awa .ge. 0.0 .and. awa .lt. 90.0) then i1 = imax - 1 i2 = 1 i3 = -1 i4 = jmax - 1 i5 = 1 i6 = -1 ke = 1 ! elseif (awa .ge. 90.0 .and. awa .lt. 180.0) then i1 = imax - 1 i2 = 1 i3 = -1 i4 = 1 i5 = jmax - 1 i6 = 1 ke = 1 ! elseif (awa .ge. 180.0 .and. awa .lt. 270.0) then i1 = 1 i2 = imax - 1 i3 = 1 i4 = 1 i5 = jmax - 1 i6 = 1 ke = 1 ! else i1 = 1 i2 = imax - 1 i3 = 1 i4 = jmax - 1 i5 = 1 i6 = -1 ke = 1 endif ! ! determine barrier influence direction index (kbr) ! if (wind_dir .ge. 337.5 .or. wind_dir .lt. 22.5) then kbr = 1 elseif (wind_dir .lt. 67.5) then kbr = 2 elseif (wind_dir.lt. 112.5) then kbr = 3 elseif (wind_dir .lt. 157.5) then kbr = 4 elseif (wind_dir .lt. 202.5) then kbr = 5 elseif (wind_dir.lt. 247.5) then kbr = 6 elseif (wind_dir .lt. 292.5) then kbr = 7 else kbr = 8 endif ! return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !> Calculates abrasion coefficients and PM10 fractions in sources of suspended soil !! !! @param soilLayrAggSta = aggreg. stability [Ln(J/Kg)] !! @param soilCrstStab = crust stabitlity [Ln(J/Kg)] !! @param soilLayrClay = soil surface fraction clay !! @param soilLayrSand = soil surface faction sand !! @param awzypt = annual average precipitation (mm) !! @param coefAbraAgg = coefficent of abrasion of aggregates (1/m) !! @param coefAbraCrst = coefficient of abrasion of crust (1/m) !! @param soilPM10FracAbraSusp = fraction pm10 in abraded supension size soil !! @param soilPM10FracEmitSusp = fraction pm10 in emitted suspension size soil !! @param soilPM10FracSaltBrkSusp = fraction pm10 in breakage from saltion size soil !********************************************************************** ! subroutine calcPM10 ver 1 3-4-2005 (was sbpm10) !********************************************************************** subroutine calcPM10 & & (soilLayrAggSta, soilCrstStab, soilLayrClay, soilLayrSand, awzypt, & & coefAbraAgg, coefAbraCrst, soilPM10FracAbraSusp, soilPM10FracEmitSusp, soilPM10FracSaltBrkSusp) ! ! + + + PURPOSE + + + ! Calculates abrasion coefficients and PM10 fractions in ! sources of suspended soil ! ! + + + ARGUMENT DECLARTAIONS + + + real soilLayrAggSta, soilCrstStab, soilLayrClay, soilLayrSand, awzypt real coefAbraAgg, coefAbraCrst,soilPM10FracAbraSusp, soilPM10FracEmitSusp, soilPM10FracSaltBrkSusp ! ! + + + ARGUMENT DEFINITIONS + + + ! soilLayrAggSta = aggreg. stability [Ln(J/Kg)] ! soilCrstStab = crust stabitlity [Ln(J/Kg)] ! soilLayrClay = soil surface fraction clay ! soilLayrSand = soil surface faction sand ! awzypt = annual average precipitation (mm) ! coefAbraAgg = coefficent of abrasion of aggregates (1/m) ! coefAbraCrst = coefficient of abrasion of crust (1/m) ! soilPM10FracAbraSusp = fraction pm10 in abraded supension size soil ! soilPM10FracEmitSusp = fraction pm10 in emitted suspension size soil ! soilPM10FracSaltBrkSusp = fraction pm10 in breakage from saltion size soil ! ! + + + LOCAL VARIABLES + + + real ratio, cla, lnc, soilLayrSilt, ppt ! ! + + + LOCAL VARIABLE DEFIINTIONS + + + ! ratio = silt/clay^2 ! cla = fraction clay with restricted range ! lnc = alog (cla) ! soilLayrSilt = soil fraction silt ! ppt = annual average precip restricted to 100-800 mm ! + + + END SPECIFICATIONS + + + ! ! calc. abrasion coefficients coefAbraAgg = exp(-2.07-0.077*soilLayrAggSta**2.5-0.119*alog(soilLayrAggSta)) !(E-102) if (soilCrstStab .eq. 0.0) then ! No crust stability value specified coefAbraCrst = 0.0 write(0,*) "Warning: Crust stability value is set to zero" else coefAbraCrst = exp(-2.07-0.077*soilCrstStab**2.5-0.119*alog(soilCrstStab)) !(E-103) end if ! ! calc. pm10 fractions in suspended soil soilPM10FracAbraSusp = 0.0116 + 0.00025/(coefAbraAgg+0.001) !(E-106) ! soilLayrSilt = 1 - soilLayrSand - soilLayrClay ratio = soilLayrSilt/(soilLayrClay + 0.0001)**2 ratio = min(300.0, ratio) soilPM10FracEmitSusp = 0.0067 + 0.0000487*ratio - 0.0000044*awzypt !(E-98) ! cla = min(0.42,max(0.017,soilLayrClay)) !restrict clay range lnc = alog(cla) ppt = min(800.0,max(100.0,awzypt)) !restrict precip range soilPM10FracSaltBrkSusp = -0.201-(0.52+(0.422+(0.1395+0.0156*lnc)*lnc)*lnc)*lnc & !(E-109) & + 0.131*exp(-ppt/175.6) ! return end !> Input subregion values of variables from other submodels to the grid points of the erosion submodel which erosion changes !! Initialize output grid array !! Calc. soil fraction of 4 dia. from asd, & rr shelter angles !********************************************************************** ! subroutine initSubr (was sbinit) !********************************************************************** subroutine initSubr ! +++ purpose +++ ! Input subregion values of variables from other submodels ! to the grid points of the erosion submodel which erosion changes ! Initialize output grid array ! Calc. soil fraction of 4 dia. from asd, & rr shelter angles ! use p1werm_def use p1erode_def use m2geo_def use m1subr_def use w1clig_def use gridmod use soilmod IMPLICIT NONE ! +++ ARGUMENT DECLARATIONS +++ ! ! +++ ARGUMENT DEFINITIONS +++ ! ! +++ PARAMETER +++ ! ! + + + GLOBAL COMMON BLOCKS + + + include 's1phys.inc' include 's1agg.inc' include 's1dbh.inc' include 's1surf.inc' include 's1sgeo.inc' include 'b1glob.inc' ! + + + LOCAL COMMON BLOCKS + + + ! include 'erosion/p1erode.inc' ! include 'erosion/e2grid.inc' ! include 'erosion/s2agg.inc' ! include 'erosion/s2surf.inc' ! include 'erosion/s2sgeo.inc' ! include 'erosion/e2erod.inc' ! ! ! + + + LOCAL VARIABLES + + + integer icsr, i, j real sfd1(maxNumSubRegs), sfd10(maxNumSubRegs), & & sfd84(maxNumSubRegs), sfd200(maxNumSubRegs) ! + + + LOCAL VARIABLE DEFINITIONS + + + ! icsr = index of current subregion ! ! + + + SUBROUTINES CALLED + + + real calcSoilMassFrac ! calcPM10 ! + + + END SPECIFICATION + + + ! ! calculate abrasion and pm10 parameters edit LH 3-4-05 do 5 icsr = 1, numSubr call calcPM10 & & (SoilRegionData(icsr)%SoilLayerData(1)%aggSta, SoilRegionData(icsr)%asoilCrstStab, & & SoilRegionData(icsr)%SoilLayerData(1)%fracClay, SoilRegionData(icsr)%SoilLayerData(1)%fracSand, & & awzypt, SoilRegionData(icsr)%acoefAbraAgg, SoilRegionData(icsr)%acoefAbraCrst, & & SoilRegionData(icsr)%asoilPM10FracAbraSusp, SoilRegionData(icsr)%asoilPM10FracEmitSusp, & & SoilRegionData(icsr)%asoilPM10FracSaltBrkSusp) ! ! calculate fraction less than diameter from asd ! determine current subregion ! do 5 icsr = 1, numSubr sfd1(icsr) = calcSoilMassFrac & & (SoilRegionData(icsr)%SoilLayerData(1)%aggGMD,SoilRegionData(icsr)%SoilLayerData(1)%aggGSD, & & SoilRegionData(icsr)%SoilLayerData(1)%aggMinSiz, SoilRegionData(icsr)%SoilLayerData(1)%aggMaxSiz, 0.01) sfd10(icsr) = calcSoilMassFrac & & (SoilRegionData(icsr)%SoilLayerData(1)%aggGMD, SoilRegionData(icsr)%SoilLayerData(1)%aggGSD, & & SoilRegionData(icsr)%SoilLayerData(1)%aggMinSiz, SoilRegionData(icsr)%SoilLayerData(1)%aggMaxSiz, 0.1) sfd84(icsr) = calcSoilMassFrac & & (SoilRegionData(icsr)%SoilLayerData(1)%aggGMD, SoilRegionData(icsr)%SoilLayerData(1)%aggGSD, & & SoilRegionData(icsr)%SoilLayerData(1)%aggMinSiz, SoilRegionData(icsr)%SoilLayerData(1)%aggMaxSiz, 0.84) ! store initial soilFracDiamLt84 soilFracDiamLt84ic = sfd84(icsr) soilFracDiamLt84ic = min(0.9999, max(soilFracDiamLt84ic,0.0001)) !set limits ! store initial soilFracDiamLt10 soilFracDiamLt10ic = sfd10(icsr) ! !^^^ tmp out ! write (*,*) ! write (*,*) 'initSubr out' ! write (*,*) 'SoilRegionData(icsr)%SoilLayerData(1)%aggGMD SoilRegionData(icsr)%SoilLayerData(1)%aggGSD& ! & SoilRegionData(icsr)%SoilLayerData(1)%aggMinSiz SoilRegionData(icsr)%SoilLayerData(1)%aggMaxSiz sfd84', & ! & SoilRegionData(icsr)%SoilLayerData(1)%aggGMD, SoilRegionData(icsr)%SoilLayerData(1)%aggGSD, & ! & SoilRegionData(icsr)%SoilLayerData(1)%aggMinSiz, SoilRegionData(icsr)%SoilLayerData(1)%aggMaxSiz(1),sfd84(icsr) ! write (*,*) !^^^ tmp end sfd200(icsr) = calcSoilMassFrac & & (SoilRegionData(icsr)%SoilLayerData(1)%aggGMD, SoilRegionData(icsr)%SoilLayerData(1)%aggGSD, & & SoilRegionData(icsr)%SoilLayerData(1)%aggMinSiz, SoilRegionData(icsr)%SoilLayerData(1)%aggMaxSiz, 2.0) 5 continue ! do 20 j = 0, jmax do 10 i = 0, imax ! determine subregion (at present only 1 subregion) ! input variables to grid cells icsr = grid(i,j)%subrReg ! soilFracDiamLt1 (i,j) = sfd1(icsr) grid(i,j)%soilFracDiamLt10 = sfd10(icsr) grid(i,j)%soilFracDiamLt84 = sfd84(icsr) grid(i,j)%soilFracDiamLt200 = sfd200(icsr) ! edit ljh - 1-22-04 grid(i,j)%soilLayrRock = SoilRegionData(icsr)%SoilLayerData(1)%fracRock ! if ifc has surface rock, 1st index maybe 0. ! grid(i,j)%soilCrstThck = SoilRegionData(icsr)%asoilCrstThck grid(i,j)%soilCrstFrac = SoilRegionData(icsr)%asoilCrstFrac grid(i,j)%soilLoosMass = SoilRegionData(icsr)%asoilLoosMass grid(i,j)%soilLoosCovFrac = SoilRegionData(icsr)%asoilLoosCovFrac ! grid(i,j)%ridgHght = SoilRegionData(icsr)%orgRidgHght !initialize RR values for each grid cell grid(i,j)%randRoug = SoilRegionData(icsr)%orgRandRoug if (grid(i,j)%randRoug < SLRR_MIN) then grid(i,j)%randRoug = SLRR_MIN else if (grid(i,j)%randRoug > SLRR_MAX) then grid(i,j)%randRoug = SLRR_MAX endif grid(i,j)%soilMassAvalDelt = 0.0 grid(i,j)%soilMassAvalAggSurf = 0.0 grid(i,j)%soilMassAvalAggSurfmx = 0.0 grid(i,j)%soilFracDiamLt84mn = 0.0 ! ! initialize output array- now in zeroGrid ! grid(i,j)%soilLossTot = 0 ! grid(i,j)%soilLossSusp = 0 ! grid(i,j)%soilLossPM10 = 0 ! 10 continue 20 continue ! return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !> Calculates the fraction of open field friction velocity !! in shelter for 8 cardinal wind directions 0, 45, ...315 !! at all interior nodes. !! Assigns the 8 fractons calculated at each node to as 3-d !! array (barrRedu(i,j,k)) for all nodes inside sim. region. !********************************************************************** ! subroutine calcOpenFldFricVelc 3/20/98 (was sbbr) !********************************************************************** subroutine calcOpenFldFricVelc ! ^^^ note must add wlen to call for tst version ! + + + PURPOSE + + + ! to calculate the fraction of open field friction velocity ! in shelter for 8 cardinal wind directions 0, 45, ...315 ! at all interior nodes ! ! to assign the 8 fractons calculated at each node to as 3-d ! array (barrRedu(i,j,k)) for all nodes inside sim. region. use p1werm_def use constants_def use gridmod use m1geo_def use m2geo_def IMPLICIT NONE ! + + + KEYWORDS + + + ! wlen ! + + + PARAMETERS + + + ! + + + ARGUMENT DEFINITIONS + + + ! wlen = windward(-) and leeward(+) distances from barrier ! calculated for each wind direction; ! ^^^ used only in test version ! + + + GLOBAL COMMON BLOCKS + + + ! + + + LOCAL COMMON BLOCKS + + + ! include 'erosion/e2grid.inc' include 'erosion/e3grid.inc' ! + + + ARGUMENT DECLARATIONS + + + ! ^^^ used only in test version ! real wlen(0:mngdpt, 0:mngdpt, 8) ! + + + LOCAL VARIABLES + + + integer i, j, k, n, tmpv real lx, ly real awar(8), a, b, c, aa, ca, waa real alpha, ceta, dmin, w, tmpa real tmpw0 ! + + + LOCAL VARIABLE DEFINITIONS + + + ! i, j, k, n = do-loop indices ! a,c = lengths to barrier ends from cell x(i,j) (m) ! b = barrier lenght (m) ! aa, ca = angles opposite lenths a and c. ! waa = angle between barrier and wind direction ! alpha = clockwise angle from d (0 deg) to a ! ceta = clockwise angle from d (0 deg) to c ! dmin = minimum distance from barrier to x(i,j) ! w = length from barrier to x(i,j) along wind direction ! tmpa = temporary (intermediate) angle calculation ! + + + SUBROUTINES CALLED + + + ! none ! + + + FUNCTION DECLARATIONS + + + real adjFricVelcForBarr, calcDist, calcAngToBarr ! + + + DATA INITIALIZATIONS + + + ! initialize variable (fraction of open field fric. vel) do 10 i = 0, imax do 10 j = 0, jmax do 10 k = 1, 8 grid(i,j)%barrRedu(k) = 1.0 10 continue ! + + + END SPECIFICATIONS + + + ! calc. 8 cardinal wind directions relative to the simulation ! region using a relative y-axis orientation of 0 deg for region. do 20 i = 1,8 awar(i) = (i-1)*45 - amasim if (awar(i) .lt. 0) then awar(i) = awar(i) + 360 elseif (awar(i) .gt. 360) then awar(i) = awar(i) - 360 endif ! convert relative wind angles to radians awar(i) = awar(i)*degtorad 20 continue ! update interior nodes do 90 i = 1, imax-1 do 90 j = 1, jmax-1 !*2 calculate distance to grid points(maybe offset from origin) lx = (i-0.5)*ix + amxsim(1,1) ly = (j-0.5)*jy + amxsim(2,1) ! ^^^tmp ! if (i .eq. 1 .and. j .eq. jmax-1)then ! write (*,*) 'calcOpenFldFricVelc output at 1, jmax-1 lx=', lx, 'ly=', ly ! write (*,*) 'ix=',ix, 'jy=', jy, 'imax=', imax, 'jmax=',jmax ! endif ! ^^^ end tmp ! barrier sweep do 70 n = 1, nbr ! find distances a, c from x(i,j) to barrier ends a = calcDist(lx, ly, amxbr(1,1,n), amxbr(2,1,n)) c = calcDist(lx, ly, amxbr(1,2,n), amxbr(2,2,n)) ! find barrier length b = calcDist(amxbr(1,1,n), amxbr(2,1,n), amxbr(1,2,n), & & amxbr(2,2,n)) ! for a grid cell on a barrier calc. fric. vel. fraction if ((a+c) .le. b) then do 50 k =1,8 w = 0 tmpw0 = grid(i,j)%barrRedu(k) grid(i,j)%barrRedu(k) = adjFricVelcForBarr(w, amzbr(n), ampbr(n), amxbrw(n)) grid(i,j)%barrRedu(k) = min (tmpw0, grid(i,j)%barrRedu(k)) ! ^^^ used only in test version ! wlen(i,j,k) = 0 50 continue go to 70 endif ! find angles from x(i,j) to barrier ends tmpa = ((b*b+c*c-a*a)/(2*b*c)) if (tmpa .lt. -1) then tmpa = -1 elseif (tmpa .gt. 1) then tmpa = 1 endif aa = acos(tmpa) tmpa = ((a*a+b*b-c*c)/(2*a*b)) if (tmpa .lt. -1) then tmpa = -1 elseif (tmpa .gt. 1) then tmpa = 1 endif ca = acos(tmpa) ! find minimum distance to barrier if (aa .gt. 90*degtorad .or. ca .gt. 90*degtorad) then dmin = min(a,c) else dmin = b*sin(aa)*sin(ca)/sin(aa+ca) endif ! is x(i,j) within 35*barrier height of barrier? ! ^^^tmp ! write (*,*) 'calcOpenFldFricVelc.for output' ! write (*,*) 'dmin=', dmin, 'i=',i,'j=',j ! ^^^ end tmp if (dmin .le. 35*amzbr(n)) then ! Find angles alpa and ceta alpha = calcAngToBarr(amxbr(1,1,n),amxbr(2,1,n),lx,ly) ceta = calcAngToBarr(amxbr(1,2,n),amxbr(2,2,n),lx,ly) ! sweep relative wind directions do 60 k = 1, 8 ! if grid cell downwind from barrier make tmpv = 1 tmpv = 0 if (abs(alpha - ceta) .lt. pi) then if((awar(k) .ge. min(alpha,ceta)) & & .and. (awar(k) .lt. max(alpha, ceta))) then tmpv = 1 endif elseif (abs(alpha-ceta) .ge. pi) then if(awar(k) .le. min(alpha,ceta)) then tmpv = 1 endif if((awar(k) .ge. max(alpha,ceta)) & & .and. (awar(k) .le. 2*pi)) then tmpv = 1 endif endif ! ^^^tmp used for test only ! if (i .eq. 34 .and. j .eq. 1) then ! write(*,*) ' barrier tmp output' ! write (*,*) 'a=',a, ' b=', b, ' c=',c ! write(*,*) 'aa=', aa, ' ca=',ca ! write(*,*) 'dmin=', dmin ! write(*,*) 'alpha=', alpha, ' ceta=', ceta ! write(*,*) 'tmpv=', tmpv ! write(*,*) 'awar(K)=', awar(k), ' k=', k ! endif ! ^^^ end tmp ! ! calc. opposite side 'w' if (tmpv .gt. 0) then ! calc. barrier to x(i,j) distance 'w' along wind vector ! calc. opposite side 'w' tmpa = abs(awar(k) - alpha) if (tmpa .gt. pi) then tmpa = 2*pi - tmpa endif waa = pi- tmpa - ca waa = max(waa, 0.0001) !edit LH 3-27-07 prevent zero ! by sin rule calc w w = sin(ca)*a/sin(waa) ! calc. min fraction unsheltered friction vel. tmpw0 = grid(i,j)%barrRedu(k) grid(i,j)%barrRedu(k) = adjFricVelcForBarr(w, amzbr(n),ampbr(n), amxbrw(n)) grid(i,j)%barrRedu(k) = min(tmpw0, grid(i,j)%barrRedu(k)) ! ^^^^ tmp used only in test version ! wlen(i,j,k) = w elseif (dmin .lt. 5*amzbr(n)) then ! is x(i,j) within 5H of barrier ! calc. frac. friction vel. for other directions tmpw0 = grid(i,j)%barrRedu(k) grid(i,j)%barrRedu(k) = adjFricVelcForBarr(-dmin, amzbr(n), ampbr(n), amxbrw(n)) grid(i,j)%barrRedu(k) = min( tmpw0, grid(i,j)%barrRedu(k)) ! ^^^ tmp used only in test version ! wlen(i,j,k) = -dmin endif 60 continue endif 70 continue 90 continue ! ^^^tmp output ! n = imax/15 ! n = max(1,n) ! write (*,*) ! write (*,*) 'BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB' ! write (*,*) 'output from calcOpenFldFricVelc.for, barrRedu(i,j,k)' ! write (*,*) ' i index values' ! write (*,310) (i, i=0,imax,n) ! do 220 k = 1,8 ! do 220 j = jmax,0,-1 ! write (*,300) (j, (barrRedu(i,j,k),i=0,imax,n)) ! 220 continue ! write (*,*) 'BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB' ! 300 format (1x, i3,1x, 30(f4.2,1x)) ! 310 format (1x, i8, 30i5) ! ^^^ end tmp end ! ********************************************************************* !> Calculates distance between two points given xy coordinates !! !! @param x1 - x coordinate of first point !! @param y1 - y coordinate of first point !! @param x2 - x coordinate of second point !! @param y2 - y coordinate of second point !! @return dist - distance between the points real function calcDist(x1,y1,x2,y2) result (dist) real x1, y1, x2, y2 dist = sqrt((x1-x2)**2+(y1-y2)**2) return end ! ********************************************************************* !> Calculates angle of vector from starting point to destination point !! !! @param x1 - x cooridinate of destination point (barrier) !! @param y1 - y cooridinate of destination point (barrier) !! @param x2 - x cooridinate of starting point (target cell) !! @param y2 - y cooridinate of starting point (target cell) real function calcAngToBarr(x1,y1,x2,y2) result (ang) real y1, x1, y2, x2, tempa,pi pi = 3.1415927 if (x1 .ne. x2) then tempa = atan (abs((y1-y2)/(x1-x2))) if ((y1 .eq. y2).and. (x1 > x2)) then !horizontal line east ang = pi/2.0 elseif ((y1 .eq. y2).and. (x1 < x2))then!horiz. line west ang = 1.5*pi elseif ((y1 > y2).and.(x1 > x2)) then !quad 1 ang = pi/2.0 - tempa elseif ((y1 < y2).and. (x1 > x2)) then !quad 2 ang = pi/2.0 + tempa elseif ((y1 < y2).and. (x1 < x2)) then !quad 3 ang = 1.5*pi - tempa elseif (( y1 > y2).and. (x1 < x2)) then !quad 4 ang = 1.5*pi + tempa endif elseif (y1 < y2) then ! vertical line to south ang = pi else ang = 2*pi ! vertical line to north endif return end ! ********************************************************************* !> Calculates reduction in friction velocity due to upwind barriers !! (ranges: porosity 0 to 0.9, distance: -5*barrHght to 50*barrHght) !! !! @param distToBarr - distance from barrier !! @param barrHght - barrier height !! @param effPor - effective porosity !! @param barrWid - barrier width !! @returns adjFrac - frictional velocity multipler (range 0.0 - 1.0) real function adjFricVelcForBarr(distToBarr, barrHght, effPor, barrWid) result (adjFrac) real distToBarr, barrHght, effPor, barrWid real a, b, c, d, adjDistToBarr, adjBarrWid, pb ! scale distance & width by barrier height adjDistToBarr = distToBarr/barrHght adjBarrWid = barrWid/barrHght ! increase effective porosity with barrier width pb = effPor + (1 - exp(-0.5*adjBarrWid))*0.3*(1-effPor) ! (E-74) ! calculate coef. as fn of porosity a = 0.008-0.17*pb+0.17*pb**1.05 ! (E-76) b = 1.35*exp(-0.5*pb**0.2) ! (E-77) c = 10*(1-0.5*pb) ! (E-78) d = 3 - pb ! (E-79) ! calc. frac. of fric. vel. adjFrac = 1 - exp(-a*adjDistToBarr**2) + b*exp(-0.003*(adjDistToBarr+c)**d) ! (E-75) ! Cap adjFrac at 1.0 if (adjFrac > 1.0) then adjFrac = 1.0 endif return end !> Sets the grid output arrays to zero !********************************************************************** ! subroutine zeroGrid (was sbigrd) !********************************************************************** subroutine zeroGrid use p1werm_def use gridmod use m2geo_def ! ! + + + PURPOSE + + + ! To set the grid output arrays to zero ! ! + + + ARGUMENT DECLARATIONS ! ! + + + GLOBAL COMON BLOCKS + + + ! compiler instr.- no warn of unreferenced symbols in include files ! ! + + + LOCAL COMMON BLOCKS + + + ! include 'erosion/e2erod.inc' ! ! + + + ARGUMENT DEFINITIONS + + + ! ! + + + SUBROUTINES CALLED + + + ! ! + + + LOCAL VARIABLES + + + integer i, j ! ! + + + END SPECIFICATIONS + + + ! do j = 0, jmax do i = 0, imax grid(i,j)%soilLossTot = 0 grid(i,j)%soilLossSusp = 0 grid(i,j)%soilLossPM10 = 0 enddo enddo return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !> Calculates the saltation/creep, suspension, and PM-10 discharge !! from each control volume !! !! @param fricVelc = friction velocity (m/s) !! @param thrsFricVelc = friction velocity threhold at emission (m/s) !! @param thrsFricVelcTrap = friction velocity threshold at transport cap.(m/s) !! @param soilFracDiamLt10 = soil fractions less than 0.1 !! @param soilFracDiamLt84 = soil fractions less than 0.84 mm !! @param soilFracDiamLt200 = soil fractions less than 2.0 mm !! @param soilFracDiamLt84ic = soil surface fraction less than 0.84 initially !! @param soilFracDiamLt10ic = soil surface fraction less than 0.10 initially !! @param asoilLayrRock = soil surface volume rock at start of event !! @param soilCrstFrac = soil fraction area crusted !! @param soilLoosCovFrac = soil fraction area of loose soil (only on crust) !! @param sarrc = soil angle r. roughness weibull c parm !! @param ridgHght = soil ridge height (mm) !! @param ridgSpacParaWind = soil ridge spacing parallel wind direction (mm) !! @param sxrg = soil ridge spacing (mm) !! @param soilLayrAggSta = soil aggregate stability (ln(J/kg)) !! @param soilCrstStab = soil crust stability (ln(j/kg)) !! @param soilLayrClay = soil fraction clay by mass !! @param soilLayrSand = soil fraction sand by mass !! @param soilLayrVeryFineSand = soil fraction very fine sand by mass (0.05-0.1 mm) !! @param soilLayrRock = soil fraction rock >2.0 mm by volume !! @param brsai = biomass stem area index !! @param bioMassHght = biomass height (m) !! @param bioFracFlatCovr = biomass fraction flat cover !! @param coefAbraAgg = coeficient of aggregate abrasion (1/m) !! @param coefAbraCrst = coefficient of crust abrasion (1/m) !! @param soilPM10FracAbraSusp = soil fraction of pm10 in abraded suspension size !! @param soilPM10FracEmitSusp = soil fraction of pm10 in emitted suspension size !! @param soilPM10FracSaltBrkSusp = soil fraction of pm10 in saltion/creep breakage !! @param lx = node spacing in x-direction (m) !! @param qi, qssi, q10i = input to C.V. of saltion, creep, pm-10 (kg/m*2) !! @param qo, qsso, q10o = output from C.V. of saltion, creep, pm-10 (kg/m*s) !! @param szc = tmp variable for prior crust thickness (mm) !! @param sz = tmp variable for roughness (mm) !! @param crlos = factor that decreases loose cover with roughness !! @param soilMassAvalDelt = change in loose mass on aggregated sfc. (kg/m^2) !! @param dmcld = change in clod mass on aggregated sfc. (kg/m^2) !! @param szv = change in height based on volume change (mm) !! @param sacd = specific area of clods per unit volume (mm^2/mm^3) !! @param flg = surface update flag, 1=update surface, 0=do not update surface !! @param soilCrstThck = consolidated crust thickness !! @param soilLoosMass = Amount of loose erodible material on crusted fraction of soil surface (kg/m^2) !! @param orgRidgSpac = Ridge spacing at beginning of day !! @param randRoug = soil random roughness (mm) !! @param time = time interval (seconds) !! @param i,j = grid cell indices !! @param imax = no. grid intervals in x-direction !! @param jmax = no. grid intervals in y-direction. !! @param soilMassAvalAggSurf = mobile soil reservoir of initial aggregated sfc.(kg/m^2) !! @param soilFracDiamLt84mn = soil fraction <0.84 mm where thrsFricVelc= fricVelc of ag.sfc. !! @param soilMassAvalAggSurfmx = max mobile soil reservoir of aggregateed sfc.(kg/m^2) !********************************************************************** ! subroutine sumOutput ver 9/2/98 file name sumOutput.for (was sbqout) !********************************************************************** subroutine sumOutput & & (flg, fricVelc, thrsFricVelc, thrsFricVelcTrap, soilFracDiamLt10, soilFracDiamLt84, & & soilFracDiamLt200, soilCrstThck, soilCrstFrac, soilLoosCovFrac, soilLoosMass, & & ridgHght, orgRidgSpac, ridgSpacParaWind, randRoug, & & soilLayrClay, soilLayrSand, & & soilLayrVeryFineSand, soilLayrRock, brsai, bioMassHght, & & bioFracFlatCovr, time, & & coefAbraAgg, coefAbraCrst, soilPM10FracAbraSusp, soilPM10FracEmitSusp, soilPM10FracSaltBrkSusp, & ! & lx, qi, qssi, q10i, i, j, imax, jmax, & & lx, qi, qssi, q10i, & ! & soilMassAvalAggSurf, soilMassAvalDelt, soilFracDiamLt84mn, soilFracDiamLt84ic, soilFracDiamLt10ic, & & soilMassAvalDelt, soilFracDiamLt84mn, soilFracDiamLt84ic, soilFracDiamLt10ic, & & asoilLayrRock,soilMassAvalAggSurfmx, & & qo, qsso, q10o ) ! ! +++ PURPOSE +++ ! calculate the saltation/creep, suspension, and PM-10 discharge ! from each control volume use p1erode_def ! ! +++ ARGUMENT DECLARATIONS +++ integer flg !Surface update flag (1=on, 0=off) real fricVelc, thrsFricVelc, thrsFricVelcTrap, soilFracDiamLt10, soilFracDiamLt84 real soilFracDiamLt200, soilCrstThck, soilCrstFrac, soilLoosCovFrac, soilLoosMass real ridgHght, orgRidgSpac, ridgSpacParaWind, randRoug real soilLayrClay, soilLayrSand,soilLayrVeryFineSand real soilLayrRock, brsai, bioMassHght, bioFracFlatCovr, time real coefAbraAgg, coefAbraCrst, soilPM10FracAbraSusp, soilPM10FracEmitSusp real soilPM10FracSaltBrkSusp real lx, qi, qssi, q10i, qo, qsso, q10o ! real soilMassAvalAggSurf, soilMassAvalDelt, soilFracDiamLt84mn, soilFracDiamLt84ic real soilMassAvalDelt, soilFracDiamLt84mn, soilFracDiamLt84ic real soilFracDiamLt10ic, asoilLayrRock, soilMassAvalAggSurfmx ! integer i, j, imax, jmax ! include 'erosion/p1erode.inc' !Parameters defining scope of erosion submodel ! +++ ARGUMENT DEFINITIONS +++ ! fricVelc = friction velocity (m/s) ! thrsFricVelc = friction velocity threhold at emission (m/s) ! thrsFricVelcTrap = friction velocity threshold at transport cap.(m/s) ! soilFracDiamLt10 = soil fractions less than 0.1 mm ! soilFracDiamLt84 = soil fractions less than 0.84 mm ! soilFracDiamLt200 = soil fractions less than 2.0 mm ! soilFracDiamLt84ic = soil surface fraction less than 0.84 initially ! soilFracDiamLt10ic = soil surface fraction less than 0.10 initially ! asoilLayrRock = soil surface volume rock at start of event ! soilCrstFrac = soil fraction area crusted ! soilLoosCovFrac = soil fraction area of loose soil (only on crust) ! sarrc = soil angle r. roughness weibull c parm ! ridgHght = soil ridge height (mm) ! ridgSpacParaWind = soil ridge spacing parallel wind direction (mm) ! sxrg = soil ridge spacing (mm) ! soilLayrAggSta = soil aggregate stability (ln(J/kg)) ! soilCrstStab = soil crust stability (ln(j/kg)) ! soilLayrClay = soil fraction clay by mass ! soilLayrSand = soil fraction sand by mass ! soilLayrVeryFineSand = soil fraction very fine sand by mass (0.05-0.1 mm) ! soilLayrRock = soil fraction rock >2.0 mm by volume ! brsai = biomass stem area index ! bioMassHght = biomass height (m) ! bioFracFlatCovr = biomass fraction flat cover ! coefAbraAgg = coeficient of aggregate abrasion (1/m) ! coefAbraCrst = coefficient of crust abrasion (1/m) ! soilPM10FracAbraSusp = soil fraction of pm10 in abraded suspension size ! soilPM10FracEmitSusp = soil fraction of pm10 in emitted suspension size ! soilPM10FracSaltBrkSusp = soil fraction of pm10 in saltion/creep breakage ! lx = node spacing in x-direction (m) ! qi, qssi, q10i = input to C.V. of saltion, creep, pm-10 (kg/m*2) ! qo, qsso, q10o = output from C.V. of saltion, creep, pm-10 (kg/m*s) ! szc = tmp variable for prior crust thickness (mm) ! sz = tmp variable for roughness (mm) ! crlos = factor that decreases loose cover with roughness ! soilMassAvalDelt = change in loose mass on aggregated sfc. (kg/m^2) ! dmcld = change in clod mass on aggregated sfc. (kg/m^2) ! szv = change in height based on volume change (mm) ! sacd = specific area of clods per unit volume (mm^2/mm^3) ! ! +++ PARAMETERS +++ real cs, cmp, ctf, cdp, c10dp parameter (cs=0.3, cmp=0.0001, ctf=1.2, cdp=0.02, & & c10dp=0.001) ! +++ PARAMETER DEFINITIONS +++ ! cs = saltation transport coef. (kg*s^2/m^4) ! cmp = mixing parameter coef. ! ct = ridge trapping coef. ! ctf = ridge trapping fraction ! cdp = deposition coef. for saltation ! c10dp = deposition coef. for pm10 ! ! +++ LOCAL VARIABLES +++ ! real fracd real cen, ceno, renb, renv, qen, qcp, ci real sfa12, surfFracNoEmit, sargc, sarrc, sfsn real fan, fanag, fancr, fancan, suspSoilFracFromEmit real sfssan, cbk,a,b,c,f,g,h,k real s, p, t1,t2, ct real cm, tmp, srrg, dmtlos, fdm real szc, sz, crlos real szv, qitmp, dmt, smasstot ! real lnslp ! integer m, n ! +++ LOCAL VARIABLE DEFINITIONS +++ ! fracd = dynamic threshold adjustment - added per LH's suggestions ! cen = coef. of emission (1/m) ! ceno = Coef. of emission for bare, loose, erodible sfc. ! renb = reduction in emission on bare soil by cover and rough. ! sargc = soil angle of ridge shelter weibull 'c' (deg.) ! sarrc = soil angle of random roughness weibull 'c'(deg.) ! sfa12 = soil fraction area with shelter angle > 12 deg. ! surfFracNoEmit = soil fraction clod & crust cover ! srrg = ratio of ridge spacing to ridge spacing parallel wind direction ! renv = reductin in emission on bare soil by flat biomass ! qen = transport capacity using emission threshold (kg/m*s) ! qcp = transport capacity using trans. cap. threshold (kg/m*s) ! ci = coef. of saltation interception by biomass stems (1/m) edit LH 8-29-05 ! a2 = intermediate calc. parameter ! fan = fraction saltation impacting clods & crust ! fanag = fraction saltation impacting clods ! fancr = fraction saltation impacting crust ! fancan = product of fan & can ie. effective abrasion coef. (1/m) ! suspSoilFracFromEmit = soil fraction to suspension from emission ! sfssan = soil fraction from abrasion to emission ! a,b,c = composite parameters for saltation discharge soln eq. ! f, g = composite parameters for suspension discharge soln eq. ! h, k = composite parameter for pm-10 discharge eq. ! s, p = intermediate calc. param. for saltation discharge soln. ! t1,t2 = intermediate calc. param. for salt. and susp. soln eq. ! tmp = intermdediate calc. param. for transport cap. soln eq. ! dmt = soil loss total in time-step (+ =depostion) (kg/m^2) ! dmtlos = soil loss/dep. for each time-step (+ = deposition) ! fdm = increase of dmtlos mass to agg. reservoir when loss ! exceeds soilLoosMass from crust ! lnslp = ln of soilMassAvalAggSurf_soilFracDiamLt84 slope, intermediate variable ! smasstot = total soil reservior mass (kg/m^2) ! +++ END SPECIFICATIONS +++ ! ! calc. fraction of area with shelter > 12 degrees sarrc = 2.3 * sqrt (randRoug) !(E-91) if ((ridgSpacParaWind .gt. 10.0) .and. (ridgHght .gt. 1.0)) then sargc = 65.4*(ridgHght/ridgSpacParaWind)**0.65 !(E-93) sfa12 =(1.- exp(-(12./sargc)**0.77))*(exp(-(12./sarrc)**0.77)) & !(E-92,E-94) & + exp(-(12./sargc)**0.77) else sfa12 = exp(-(12./sarrc)**0.77) endif ! ! ! The following changes were made based upon Larry Hagen's email ! me on Fri, Jan. 31, 2003 to represent dynamic threshold instead ! assuming a static threshold - LEW ! calc. change to dynamic transport ! edit ljh 1-24-05 fracd = 0.05 * (1.0 - sfa12) !(E-88) ! ! calc. transport cap. for emission edit 5/30/01 LH qen = cs * fricVelc * fricVelc * (fricVelc - (thrsFricVelc-fracd)) ! calc. transport cap. for trapping qcp = cs * fricVelc * fricVelc * (fricVelc - (thrsFricVelcTrap-fracd)) if (qcp .lt. 0.0) then qcp = 0.001 endif ! ! store tmp qi qitmp = qi ! test for trap region with both saltation & suspension deposition if (qen < 0.0001) then !edit ljh 1-24-05 qen = 0.00001 qo = 0.0 qi = qo !no abrasion when no transport edit 8-30-06 ct = 1.0 !traps all incoming saltation discharge fancan = 0.0 qsso = qssi*exp(-cdp*lx) !(E-110) q10o = q10i*exp(-c10dp*lx) !(E-111) go to 90 ! test for saltation deposition region with suspension emission ! We have set qi to qen to calc. suspension and abrasion ! and then used the real qi value to calc. the deposition. ! 12-5-2000 LH) ! elseif (qen .le. qi) then qi = qen endif ! ! Calc. emission params for saltation/creep ! emission coef. for bare, loose, erodible soil (1/m) ceno = 0.06 ! fraction emission coef. reduction by flat biomass renv = 0.075 + 0.934*exp(-bioFracFlatCovr/0.149) !(E-90) ! fraction emission coef. reduction by roughness and area not emitting if( soilFracDiamLt84 .le. soilFracDiamLt84mn) then ! edit LH 6-27-06 surfFracNoEmit = ((1.-soilCrstFrac)*(1.- 0.) + soilCrstFrac - soilLoosCovFrac*soilCrstFrac)*(1.-soilLayrRock) & !(E-80) & + soilLayrRock else surfFracNoEmit = ((1.-soilCrstFrac)*(1.-soilFracDiamLt84) + soilCrstFrac - soilLoosCovFrac*soilCrstFrac) & & *(1.-soilLayrRock) + soilLayrRock endif ! ! edit 7-17-01 LH renb = (-0.051 + 1.051*exp(-surfFracNoEmit/0.33050512))*(1 - sfa12) !(E-96) mod cen = ceno*renv*renb !(E-89) ! if (cen .lt. 0.0001) cen = 0.0 ! soil fraction suspension in emitted soil; lower for loose on cr suspSoilFracFromEmit = (soilFracDiamLt10/(soilFracDiamLt200+0.001))*(1. - soilCrstFrac*0.8) !(E-97,E-17) mod ! soil mixing parameter for suspension cm = cmp*suspSoilFracFromEmit ! ! calc. trap params. for saltation/creep ! calc. coef of stem interception ! for plants above 1 mm height edited 3-29-01 LH if (bioMassHght .lt. 0.001) then ci = 0. else if (orgRidgSpac .gt. 10) then srrg = sqrt(orgRidgSpac/ridgSpacParaWind) else srrg = 1.0 endif ci = 0.005*srrg*(1 - exp(-0.5*brsai/bioMassHght))*(1 - exp(-50*bioMassHght)) !LH 8-29-05 !(E-115) endif ! ! Calc. abrasion params. for saltation/creep ! soil fraction saltation/creep in saltation sfsn = 1 ! ! fraction of saltation abrading clods, crust and loose fan = (exp(- 4.0*bioFracFlatCovr))*(exp(-3*soilLayrRock))*sfsn !(E-99) mod fan = amax1(fan,0.) ! fraction abrading aggregates fanag = (1. - soilFracDiamLt84)*(1. - soilCrstFrac)*fan !(E-100) ! fraction abrading crust fancr = (soilCrstFrac - soilLoosCovFrac*soilCrstFrac)*fan !(E-101) ! calc. coef. of abrasion ! coefAbraAgg = exp(-2.07-0.077*soilLayrAggSta**2.5-0.119*alog(soilLayrAggSta)) ! coefAbraCrst = exp(-2.07-0.077*soilCrstStab**2.5-0.119*alog(soilCrstStab)) ! sum fancan = fancr*coefAbraCrst + fanag*coefAbraAgg !test3 if (fancan .lt. 0.0001) fancan = 0.0001 ! calc fraction of abraded of suspension size sfssan = (0.4 - 4.83*soilLayrClay + 27.18*soilLayrClay**2 & !(E-105) & - 53.7*soilLayrClay**3 + 42.25*soilLayrClay**4 - 10.7*soilLayrClay**5) ! ! set upper limit on sfssan for high sand content sfssan = min(sfssan, (1.0 - (soilLayrSand-soilLayrVeryFineSand))) ! ! calc suspension breakage coef. for saltation creep ! edit 6-9-01 LH cbk = 0.11*coefAbraAgg*(1. - (soilLayrSand-soilLayrVeryFineSand))*sfsn !(E-107) mod ! ! calc total trap coef. if ((qcp .lt. qen) & & .and. ((randRoug .gt. 10.) .or. (ridgHght .gt. 50.))) then ct = (1. - suspSoilFracFromEmit)*cen*ctf*(qen - qcp)/qen !(E-113,E-114) mod else ct = 0. endif ! ! assemble composite params. for saltation/creep a = cen *(1. - suspSoilFracFromEmit)*qen !(E-8) b = (1. - sfssan)*fancan - (1. - suspSoilFracFromEmit)*cen - cbk - ci - ct !(E-9) c = (1. - sfssan)*fancan*1./qen !(E-10) ! ! solve for saltation/creep out ! collect variables s = sqrt(4.*a*c + b**2.) !(E-12) p = (-2.*c*qi + b)/s !(E-13) ! change p-values that are out of range ! math range restriction: (-1 < p < 1) edit 7-18-01 LH if (p .le. -1.) then t1 = -20 elseif (p .ge. 1.) then t1 = 20 else t1 = (s/2.0)*(-lx) + 0.5*alog((1. + p)/(1. - p)) endif ! calculate atanh qo = (s/(2.*c))*(-tanh(t1) + b/s) !(E-11) mod ! ! assemble composite params. for suspension f = suspSoilFracFromEmit*cen*qen*(1.0 - 0.1*brsai) !edit 8-29-05 for suspension interception !(E-22) mod ! added last term to intercept some ss LH edit 6-11-01, changed 8-29-05 g = (sfssan*fancan - suspSoilFracFromEmit*cen + cbk + cm)*(1 - 0.1*brsai) !(E-23) mod ! ! ^^^ tmp out ! set counter ! m = (imax-1)/8 ! m = max(1,m) !^^^ tmp output ! if (fricVelc .gt. 1.1) then ! if (j .eq. 10) then ! do 2 i=1,10 ! write (*,*) ! write (*,*) 'output from sumOutput: i=',i,'j=',j ! write (*,*) 'a=', a, 'b=', b, 'c=',c ! write (*,*) 'g=', g, 'f=', f ! write (*,*) 'cen=', cen, 'qen=', qen ! write (*,*) 'suspSoilFracFromEmit=', suspSoilFracFromEmit, 'sfssan=',sfssan, ! & 'fancan=',fancan ! write (*,*) 'cbk=',cbk, 'ci=',ci,'ct=',ct,'cm=',cm ! 2 continue ! endif ! endif ! ^^^ end tmp out ! ! trap to prevent over range of exp if ((s*lx) .gt. 40.0) then s = 40.0/lx endif ! solve for suspension out ! collect variables if( p .gt. 1.0 ) then t2 = alog(2.0) else if( p .lt. -1 ) then t2 = alog(exp(s*lx)*2.0) else t2 = alog(exp(s*lx)*(1.-p)+ p + 1.) end if qsso=qssi+(1./(2.*c))*((-g*s+g*b+2.*f*c)*lx+2.*g*(-alog(2.0)+t2)) !(E-24) mod ! ! assemble composite params. for PM-10 h = soilPM10FracEmitSusp*f !(E-28) mod k = soilPM10FracAbraSusp*sfssan*fancan - soilPM10FracEmitSusp*suspSoilFracFromEmit*cen + soilPM10FracSaltBrkSusp*cbk & & + soilPM10FracEmitSusp*cm !(E-29) ! ! solve for PM-10 out (similar to suspension out) q10o=q10i+(1./(2.*c))*((-k*s+k*b+2.*h*c)*lx+2.*k*(-alog(2.0)+t2)) !(E-30) mod ! SURFACE UPDATE EQUATIONS !now added as a commandline argument to tsterode - LEW ! execute this section if flag is set 90 if (flg .eq. 1) then ! ! calc change in loose surface soil for time step (+ = surface gain) ! terms: - total loss, + created by abrasion dmt = ((-(qo-qitmp) - (qsso-qssi))/lx)*time !(E-31) dmtlos = dmt & & + (fancan*qi)*time !edit LH 3-1-05 !(E-31) cont ! set initial to zero fdm = 0.0 ! coef. for cover of loose mass (same as SOIL eq. s-24,2-25) sz = amax1(ridgHght, 4.0*randRoug) !(E-35) crlos = exp(-0.08*sz**0.5) !(E-34) ! ! execute this section if crust present if (soilCrstFrac .gt. 0.01) then ! increase loss per unit area crust if non-emitting aggregate sfc. if (dmtlos .lt. 0 .and. soilFracDiamLt84 .le. soilFracDiamLt84mn) then !edit LH 2-28-05 dmtlos = dmtlos*(1.0/soilCrstFrac) endif ! ! update loose on crust ! edit LH 2-24-05 soilLoosMass = soilLoosMass + dmtlos !(E-32) if (soilLoosMass < 0) then ! loose removal exceeded fdm = soilLoosMass* soilCrstFrac/(1.0001-soilCrstFrac) !remove extra from ag. area soilLoosMass = 0.0 endif ! ! update cover of loose mass if (soilLoosMass > 0) then tmp = 3.5*soilLoosMass**1.5 !(E-33) if(tmp.gt.80.) then !test and prevent underflow condition soilLoosCovFrac = crlos else soilLoosCovFrac = (1.0 - exp(-tmp))*crlos endif else soilLoosCovFrac = 0 endif ! if (fancan > 0.0) then ! abrasion occurs edit LH 1-9-07 ! update fancr if (dmtlos > 0.0 )then fancr = (soilCrstFrac - soilLoosCovFrac*soilCrstFrac)*fan !edit LH 2-16-05 endif ! ! update crust thickness szc = amax1(0.001, soilCrstThck) soilCrstThck = soilCrstThck - (fancr*coefAbraCrst*qi/(1.4*soilCrstFrac))*time !edit LH 3-1-05 !(E-36) soilCrstThck = amax1(soilCrstThck, 0.00) ! ! update crust (consolidated) zone cover soilCrstFrac = soilCrstFrac*soilCrstThck/szc !(E-37) soilCrstFrac = amax1(0.0, soilCrstFrac) endif else soilLoosMass = 0. soilLoosCovFrac = 0. endif ! ! execute this section if clods are present if (soilCrstFrac .lt. 0.99) then ! ! change in loose mass on aggregated and rock surface,(+)=increase soilMassAvalDelt = soilMassAvalDelt + dmtlos + fdm !(E-38) mod ! update soilFracDiamLt84 when net loose soil gain on agg. & rock ! added smasstot to simplify soilFracDiamLt84 & soilFracDiamLt10 equations smasstot= soilMassAvalAggSurfmx/(soilFracDiamLt84ic*(1.001-asoilLayrRock)) !edit LH 3-26-07 if (soilMassAvalDelt .lt. 0.0 ) then soilFracDiamLt84 = (soilMassAvalAggSurfmx + soilMassAvalDelt)/smasstot else soilFracDiamLt84=(soilMassAvalAggSurfmx+soilMassAvalDelt)/(smasstot + soilMassAvalDelt) ! bad LH 3-24-07 soilFracDiamLt84 = min(soilFracDiamLt84,(soilMassAvalDelt/0.4)) endif ! ! ! set limits on soilFracDiamLt84 soilFracDiamLt84 = max(0.0, soilFracDiamLt84) soilFracDiamLt84 = min(0.9999, soilFracDiamLt84) ! ! update rock cover based on dmt (soil loss(-) or gain(+)) if (asoilLayrRock .gt. 0.0 .and. asoilLayrRock .lt. 1.0) then soilLayrRock = soilLayrRock - 7.5*((1-soilLayrRock)/(1-asoilLayrRock)) & !(E-52) mod & *(dmt/(1200*(1.001-soilLayrRock))) soilLayrRock = min( 1.0, max(soilLayrRock, 0.0) ) endif ! ! update soilFracDiamLt200 soilFracDiamLt200 = (2 - soilFracDiamLt84)*soilFracDiamLt84 !(E-53) soilFracDiamLt200 = amin1(soilFracDiamLt200, 1.0) soilFracDiamLt200 = amax1(0.001,soilFracDiamLt200) ! update soilFracDiamLt10 edit LH 3-9-05 ! if (dmtlos > 0 ) then ! soilFracDiamLt10 = soilFracDiamLt10 - soilFracDiamLt10*dmtlos !(E-54) ! soilFracDiamLt10 = amax1(.01,soilFracDiamLt10) ! endif if (soilMassAvalDelt < 0.0) then soilFracDiamLt10 = soilFracDiamLt10ic*soilFracDiamLt84/soilFracDiamLt84ic !(E-55) mod else soilFracDiamLt10 = soilFracDiamLt10ic*smasstot/(smasstot+soilMassAvalDelt) !edit 3-26-07 LH endif ! soilFracDiamLt10 = soilFracDiamLt84*0.01 ! check that cumulative distribution points are always rational if( soilFracDiamLt84 .ge. soilFracDiamLt200 ) then soilFracDiamLt84 = min( soilFracDiamLt84, 0.9999999*soilFracDiamLt200) !write(*,*) 'sumOutput: soilFracDiamLt84 >= soilFracDiamLt200' end if if( soilFracDiamLt10 .ge. soilFracDiamLt84 ) then soilFracDiamLt10 = min( soilFracDiamLt10, 0.9999999*soilFracDiamLt84) !write(*,*) 'sumOutput: soilFracDiamLt10 >= s84' end if ! fanag = (1-soilFracDiamLt84)*(1-soilCrstFrac)*fan endif ! ! update surface roughness ! ! rate of volume change caused by emission ! (used bulk densities of 1.2 for loose and 1.4 for crust ! to calc. mm depth per m^2 of area) ! szv= cen*(qen - qi)/1.2 !(E-56) ! ! if trapping then emission lowers roughness, ! i.e. it comes from highest areas. ! If (ct .gt. 0.) then ! szv= -szv ! endif ! calc. change in mean surface depth ! terms: emission, trapped, abrasion, & abrasion not emitted ! szv = (szv - ct*qi/1.2 - fancan*qi/1.4 & ! & -(1.0 - sfssan)*fancan*qi*qi/(qen*1.2))*time !(E-57) ! if (szv .gt. 0.0) then ! slow roughness increase edit LH 8-30-06 ! szv = szv*0.5 ! endif ! ! new trial roughness update edit LH 8-31-06 if (ct .gt. 0.0) then !trapping if (dmtlos .gt. 0.0) then !deposition (-roughness) szv = -2.0*dmtlos/1.2 else szv = dmtlos/1.2 !emission (-roughness) endif else if (dmtlos .gt. 0.0) then szv = - dmtlos/1.2 !deposition (-roughness) else szv = 0.5*dmtlos/1.2 ! emission (+ roughness) endif endif szv = szv-(2.0*(fancan*qi/1.4))*time ! abrasion (- roughness) ! ! update ridge height if (ridgHght .gt. 10.) then ridgHght = ridgHght + szv ! set lower limit on ridge height ridgHght = amax1(ridgHght, 0.0) endif ! update random roughness randRoug = randRoug + szv/4.0 ! set lower limit on randRoug randRoug = max(randRoug, SLRR_MIN) ! endif ! 100 continue qi = qitmp return end !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! (was sbsfdi) !> Calculates soil mass fraction (soilMassFrac) < diameter (aggDistInpDiam) given modified lognormal distribution parameters !! !! @param aggDistGeoMean - aggregate distribution geometric mean diameter (mm). !! @param aggDistGeoStd - aggregate distribution geometric standard deviation. !! @param aggDistLoLmt - aggregate distribution lower limit (mm). !! @param aggDistUpLmt - aggregate distribution upper limit (mm). !! @param aggDistInpDiam - soil diameter in distribution (mm) !! @param soilMassFrac - soil mass fraction < aggDistInpDiam real function calcSoilMassFrac (aggDistGeoMean, aggDistGeoStd, aggDistLoLmt, aggDistUpLmt, aggDistInpDiam) & result (soilMassFrac) ! +++ PURPOSE +++ ! calc soil mass fraction (soilMassFrac) < diameter (aggDistInpDiam) ! given modified lognormal distribution parameters ! +++ ARGUMENT DECLARATIONS +++ ! real aggDistGeoMean, aggDistGeoStd, aggDistLoLmt, aggDistUpLmt, aggDistInpDiam, soilMassFrac real aggDistGeoMean, aggDistGeoStd, aggDistLoLmt, aggDistUpLmt, aggDistInpDiam ! +++ ARGUMENT DEFINITIONS +++ ! aggDistGeoMean - aggregate distribution geometric mean diameter (mm). ! aggDistGeoStd - aggregate distribution geometric standard deviation. ! aggDistLoLmt - aggregate distribution lower limit (mm). ! aggDistUpLmt - aggregate distribution upper limit (mm). ! aggDistInpDiam - soil diameter in distribution (mm) ! soilMassFrac - soil mass fraction < aggDistInpDiam ! +++ LOCAL VARIABLES +++ real slt ! +++ FUNCTIONS CALLED+++ real erf ! +++ END SPECIFICATIONS +++ ! calc soil mass < aggDistInpDiam if (aggDistInpDiam .lt. aggDistUpLmt .and. aggDistInpDiam .gt. aggDistLoLmt) then slt = ((aggDistInpDiam - aggDistLoLmt)*(aggDistUpLmt - aggDistLoLmt))/((aggDistUpLmt - & aggDistInpDiam)*aggDistGeoMean) !(E-81) soilMassFrac = 0.5*(1 + erf(alog(slt)/(sqrt(2.0)*alog(aggDistGeoStd)))) !(E-82) elseif (aggDistInpDiam .ge. aggDistUpLmt) then soilMassFrac = 1.0 else soilMassFrac = 0.0 endif return end !> Set field zero plane displacement equal to Anemometer zero plane displacement when anemometer at the field site !********************************************************************** ! subroutine calcZeroPlanDisp (was sbzdisp) !********************************************************************** subroutine calcZeroPlanDisp (ridgHght, cropRowSpac, cropFurrFlg, aeroAnemFlg, & & resiLAI, resiSAI, resiAvgHght, cropLAI, cropSAI, cropHght, & & dispHghtWeatStat, zeroPlanDisp) ! +++ PURPOSE +++ ! Calc. zero plane displacement (mm) ! set field zero plane displacement equal to Anemometer zero plane ! displacement when anem. at the field site, ie. aeroAnemFlg = 1 ! using equation from: Raupach, M.R. 1994. Simplified Expressions for ! Vegetation Roughness Length and Zero-Plane Displacement as functions ! of Canopy Height and Area index. Boundary Layer meteorology 71:211-216. ! +++ ARGUMENT DECLATION +++ real ridgHght, cropRowSpac integer cropFurrFlg, aeroAnemFlg real resiLAI, resiSAI, resiAvgHght real cropLAI, cropSAI, cropHght real dispHghtWeatStat, zeroPlanDisp ! +++ ARGUMENT DEFINITIONS +++ ! ridgHght - ridge height (mm) ! cropRowSpac - crop row spacing (m) ! cropFurrFlg - flag=0 - crop planted in furrow bottom ! flag=1 - crop planted on ridge top ! aeroAnemFlg - flag=0 - weather measurements are from a distant station ! flag=1 - weather measurements are in field ! resiLAI - residue leaf area index (total)(m2/m2) ! resiSAI - residue stem area index (total)(m2/m2) ! resiAvgHght - composite average residue height (m) ! cropLAI - crop leaf area index (m2/m2) ! cropSAI - crop stem area index (m2/m2) ! cropHght - crop height (m) ! cropFurrFlg - flag=0 - crop planted in furrow bottom ! flag=1 - crop planted on ridge top ! dispHghtWeatStat - zero plane displacement at weather location (mm) ! zeroPlanDisp - zero plane displacement at location (mm) ! +++ FUNCTIONS CALLED real biodrag ! +++ LOCAL VARIABLES +++ real bioMassHght, bsai ! +++ LOCAL VARIABLE DEFINITIONS +++ ! bioMassHght - biomass height (mm) ! bsai - biomass silhouette area index ! +++ INCLUDE FILES+++ include 'p1unconv.inc' ! +++ PARAMETERS +++ ! +++ END SPECIFICATIONS +++ ! calculate "effective" biomass silhouette area index bsai = biodrag( resiLAI, resiSAI, cropLAI, cropSAI, cropFurrFlg, & & cropRowSpac, cropHght, ridgHght ) ! find maximum biomass height and convert to mm bioMassHght = max(resiAvgHght, cropHght) * mtomm ! use silhouette area index and biomass height to find zero plane displacement if( bsai .gt. 1.0e-10 ) then zeroPlanDisp = bioMassHght * (1.0 - (1.0 - exp( -(15.0*bsai)**0.5) ) & & / (15.0*bsai)**0.5) else zeroPlanDisp = 0.0 end if if (aeroAnemFlg .eq. 1) then ! anemom. in field set weather displacement to field displacement dispHghtWeatStat = zeroPlanDisp endif return end