!$Author$ !$Date$ !$Revision$ !$HeadURL$ subroutine sinit (daysim, & & bhtsmx, bhrwc, bsfom, bszlyt, & & bslay, bsfsan, bsfsil, bsfcla, & & bszrgh, bszrr, bsfcce, bsfcec, & & cump, dcump, bsk4d, & & bhtmx0, bhrwc0, szlyd, & & bszrr0, bszrh0, & & bseagm, bseagmn, bseagmx, & & bslmin, bslmax, & & rain, snow, sprink, & & bhzirr, bszrho, & & bhlocirr, bhzsmt, bszrro, & & bsdsblk, bwzdpt, bwtdav, trigger) ! + + + PURPOSE + + + ! SOIL submodel for the Wind Erosion Prediction System model. ! daily initialization of soil properties ! (SURFACE: roughness, ridges, crust, and erodible material, ! and the LAYERS: aggregate size distribution, agg stability, and density). ! for more details on equations and processes, see SOIL SUBMODEL TECHNICAL ! DESCRIPTION. ! + + + KEY WORDS + + + ! wind erosion, soil processes, surface process, layer process ! + + + GLOBAL COMMON BLOCKS + + + include 'p1werm.inc' include 'wpath.inc' include 'm1subr.inc' include 'm1sim.inc' include 'm1flag.inc' ! + + + ARGUMENT DECLARATIONS + + + integer daysim real bhtsmx(mnsz), bhrwc(mnsz), bsfom(1:mnsz), bszlyt(mnsz) integer bslay real bsfsan(1:mnsz), bsfsil(1:mnsz), bsfcla(1:mnsz) real bszrgh, bszrr, bsfcce(1:mnsz), bsfcec(1:mnsz) real cump, dcump, bsk4d(mnsz) real bhtmx0(mnsz), bhrwc0(mnsz), szlyd(0:mnsz) real bszrr0, bszrh0 real bseagm(mnsz), bseagmn(mnsz), bseagmx(mnsz) real bslmin(mnsz),bslmax(mnsz) real rain, snow, sprink real bhzirr, bszrho real bhlocirr, bhzsmt, bszrro real bsdsblk(mnsz), bwzdpt, bwtdav integer trigger(bslay) ! + + + ARGUMENT DEFINITIONS + + + ! daysim - an index for the day of simulation. ! bhtsmx - layer maximum temperature of today in C. ! bhrwc - soil water content for today, kg/kg. ! bsfom - layer fraction of organic matter. ! bszlyt - layer thickness, mm. ! bslay - number of soil layers ! bsfsan - layer fraction of sand. ! bsfsil - layer fraction of silt. ! bsfcla - layer fraction of clay. ! bszrgh - ridge height, mm. ! bszrr - random roughness height, mm ! bsfcce - soil fraction calcium carbonate equivalent ! bsfcec - soil cation exchange capacity (cmol/kg) ! cump - cumulative (rain + sprinkler + snow-melt) to current ! day from day 1 or time of last tillage ! dcump - total rain + sprinkler + snow-melt for current day. ! bsk4d - drying process coef. to calc. aggregate stability ! bhtmx0 - layer maximum temperature of yesterday. in C ! bhrwc0 - soil water content for yesterday. mass basis kg/kg. ! szlyd - depth to bottom of each soil layer, mm ! bszrr0 - prior day random roughness, mm ! bszrh0 - prior day ridge height, mm ! bseagm - mean agg stability, ln(J/kg) ! bseagmn - minimum agg stability, ln(J/kg) ! bseagmx - maximum agg stability, ln(J/kg) ! bslmin - min value of aggregate gmd ! bslmax - max value of aggregate gmd ! rain - water added to soil as rain. ! snow - water equivalent added to soil surface as snow, mm. ! sprink - water added to soil as sprinkler irrigation, mm. ! bhzirr - irrigation water applied, mm/day. ! bszrho - ridge height right after tillage, mm. ! bhlocirr - location of irrigation application, mm. ! + means above the soil surface ! - means below the soil surface ! soil surface reference is the bottom of the furrow ! bhzsmt - snowmelt, mm/day. ! bszrro - random roughness height right after tillage, mm. ! bsdsblk - consolidated soil bulk density by layer, Mg/m^3 ! bwzdpt - rainfall depth (mm) ! bwtdav - Average daily air temperature (deg C) ! trigger - bitmapped integer showing the state of soil property change ! condition triggers for output into the layer detail file ! the 0 at the end of bhtmx0, bhrwc0, bszrr0, bszrh0 refer to ! prior day values of: ! max temperature, soil water content, random roughnes & ridge height ! bszrro , bszrho are right-after-tillage ! + + + LOCAL VARIABLES + + + real sf84m(mnsz), sf84sd(mnsz), scecr real tsfom, tsfcce, tsfsacl integer ldx ! + + + LOCAL DEFINITIONS + + + ! sf84m - mean of fraction agg. < 0.84 mm ! sf84sd - standard deviation of fraction agg. < 0.84 mm ! scecr - ratio of clay fraction cation exchange capacity ! to percent clay ! tsfom - temporary layer fraction of organic matter. ! tsfcce - temporary soil fraction calcium carbonate equivalent ! tsfsacl - temporary layer fraction of clay. ! ldx - layer index ! + + + FUNCTIONS CALLED + + + ! + + + SUBROUTINES CALLED + + + ! + + + END SPECIFICATIONS + + + ! check for first day if (daysim .eq. 1) then ! set up tillage check bszrr0 = - 1.0 bszrh0 = - 1.0 ! initialize previous day temperature and water content do ldx = 1, bslay bhtmx0(ldx) = bhtsmx(ldx) bhrwc0(ldx) = bhrwc(ldx) end do endif szlyd(0) = 0.0 do ldx = 1, bslay ! calc. depth to bottom of each layer szlyd(ldx) = szlyd(ldx-1) + bszlyt(ldx) ! zero out trigger condition array trigger(ldx) = 0 end do ! if tillage (or anything else outside of soil submodel) ! changed roughness or ridge height then update if ((bszrr0.ne.bszrr).or.(bszrh0.ne.bszrgh)) then ! store initial or after tillage surface roughness bszrro = bszrr bszrho = bszrgh ! set cumulative precip to zero cump = 0.0 ! store/calculate initial layer values do 10 ldx = 1, bslay ! store initial water content & yesterday's temperature bhtmx0(ldx) = bhtsmx(ldx) bhrwc0(ldx) = bhrwc(ldx) ! calc. mean, min, and max agg. stability ! (eq. S-26, S-27, S-28, *** sb S-25,26,27) if (bsfcla(ldx) .gt. 0.6) then bseagm(ldx) = 3.484 else bseagm(ldx) = -16.73 - 46.629*bsfcla(ldx)**2 & & + 23.514*bsfcla(ldx)**3 & & + 17.519*exp(bsfcla(ldx)) endif bseagmn(ldx) = bseagm(ldx) - 2*(0.16)*bseagm(ldx) bseagmx(ldx) = bseagm(ldx) + 2*(0.16)*bseagm(ldx) ! calc. mean and standard deviation of fraction agg. < 0.84 mm ! (eq. S-42, S-43, *** sb S-37, S-38) ! *** sf84m(ldx) = 0.2902 + 0.31 * bsfsan(ldx) + 0.17 * bsfsil(ldx) + ! *** & 0.0033*bsfsan(ldx)/bsfcla(ldx) - 4.66*bsfom(ldx) - 0.95*bsfcce(ldx) ! *** debugging fix, 1st try ! *** sf84m(ldx) = 0.2902 + 0.31 * bsfsan(ldx) + 0.17 * bsfsil(ldx) + ! *** & 0.0033*bsfsan(ldx)/bsfcla(ldx) - 4.66*bsfom(ldx) ! eodf ! *** debugging fix, 2nd try ! clamping upper limits on variables to keep them from forcing sf84m negative ! note that this needs correcting by a more robust regression equation tsfom = bsfom(ldx) if (tsfom.gt.0.03) tsfom = 0.03 tsfcce = bsfcce(ldx) if (tsfcce.gt.0.2) tsfcce = 0.2 if (bsfcla(ldx).eq.0) tsfsacl = 40. if (bsfcla(ldx).ne.0) tsfsacl = bsfsan(ldx) / bsfcla(ldx) if (tsfsacl.gt.40) tsfsacl = 40. sf84m(ldx) = 0.2902 + 0.31 * bsfsan(ldx) + 0.17 * & & bsfsil(ldx) + 0.0033*tsfsacl - 4.66*tsfom - 0.95*tsfcce ! *** eodf sf84sd(ldx) = (0.41 - 0.22*bsfsan(ldx))*sf84m(ldx) ! *** write(*,*) ' sf84m(ldx), sf84sd(ldx) ', ldx, sf84m(ldx), sf84sd(ldx) ! calc. min and max values of geom. mean agg. diameter (eq. S-45, S-46) bslmin(ldx) = exp(3.44 - 7.21*(sf84m(ldx) + 2.0 * & & sf84sd(ldx))) if (bslmin(ldx) .lt. 0.025) bslmin(ldx) = 0.025 bslmax(ldx) = exp(3.44 - 7.21*(sf84m(ldx) - 2.0 * & & sf84sd(ldx))) if (bslmax(ldx) .gt. 31.0) bslmax(ldx) = 31.0 if(bslmin(ldx).ge.bslmax(ldx)) write(*,*) 'sinit:min.gt.max' ! *** write(*,*) 'bslmin(ldx),bslmax(ldx)',ldx,bslmin(ldx),bslmax(ldx) ! calc. ratio of clay cation exchange capacity to percent clay (eq. S-53) scecr = (bsfcec(ldx) - bsfom(ldx) * (142. + 0.17 * & & szlyd(ldx)))/ (bsfcla(ldx) * 100.0 + 0.0001) if (scecr .lt. 0.15) scecr = 0.15 if (scecr .gt. 0.65) scecr = 0.65 ! *** remove calculation of cbd; replace with original cbd calc from inpsub ! *** sdbko(ldx) = bsdsblk(ldx) ! ***c calc. consolidated bulk density (eq. S-52) ! *** sdbko(ldx) = 1.514 + 0.25*bsfsan(ldx) - ! *** * 13.*bsfsan(ldx)*bsfom(ldx) -6.0*bsfcla(ldx)* ! *** * bsfom(ldx) - 0.48*bsfcla(ldx)*scecr ! *** if (sdbko(ldx) .gt. 2.2) sdbko(ldx) = 2.2 ! *** if (sdbko(ldx) .lt. 0.5) sdbko(ldx) = 0.5 ! ***c calc. increase in consolidated bulk density with depth ! ***c (note the depths change slightly with time, but current ! ***c update only occurs with tillage.) ! ***c (eq. S-54) ! *** debugging fix ! *** sdbko(ldx) = sdbko(ldx)*(0.975+ 1.931/ ! *** * (1+exp(-(szlyd(ldx)-506.8)/118.5))) ! *** eodf ! set stability drying process coefficient: bsk4d(ldx) = 0.46 - 0.23 * exp(-(szlyd(ldx-1) + & & (szlyd(ldx) - szlyd(ldx-1))/2.0)/88.57) 10 continue endif !23456789*23456789*23456789*23456789*23456789*23456789*23456789*2345 ! ! set stability freezing and wetting process coefficients: ! ! initialize rain and snow variables rain = 0.0 snow = 0.0 ! Determine if precip if (bwzdpt .gt. 0.0) then ! Determine if precip. is rain or snow ! (note HYDROLOGY may do this in future using ELS results) if (bwtdav .ge. 0.0) rain = bwzdpt if (bwtdav .lt. 0.0) snow = bwzdpt endif ! add irrigation to cumulative precipitation based on application ! height with respect to ridge height if (bhlocirr .ge. bszrgh ) then ! irrigation is being applied from above ridge height ! add full amount to degrade ridge height and random roughness sprink = bhzirr else if (bhlocirr .ge. 0.0 ) then ! irrigation application is below ridge height ! partially include reducing degradation (like furrow irrigation) sprink = bhzirr * bhlocirr / bszrgh else ! irrigation application underground ! no degradation of ridge height or randowm roughness sprink = 0.0 endif ! Calc. daily and cumulative (rain + sprinkler irrigation + snowmelt) dcump = rain + sprink + bhzsmt ! (note: cump not used in calc., but useful as ouptput for validation) cump = cump + dcump end