c$Author: jt $ c$Date: 2002-03-11 20:40:25 $ c$Revision: 1.8 $ c$Source: /weru/cvs/weps/weps.src/soil/sinit.for,v $ subroutine sinit (daysim, & bhtsmx, & bhrwc, & bsfom, bszlyt, bslay, & bsfsan, bsfsil, bsfcla, & bszrgh, & bszrr, & bsfcce, bsfcec, * cump, dcump, * bsk4d, * bhtmx0, bhrwc0,szlyd, * bszrr0, bszrh0, bseags, * bseagm, bseagmn, bseagmx, * bslmin, bslmax, * rain, snow, sprink, * bwtdav, bhzirr, bszrho, * bm0irr, bhzsmt, bszrro, * bsdsblk) c + + + PURPOSE + + + c soil submodel for the Wind Erosion Research Model. c update the SOIL (SURFACE: roughness, ridges, crust, and erodible material, c and the LAYERS: aggregate size distribution, agg stability, and density). c for more details on equations and processes, see SOIL SUBMODEL TECHNICAL c DESCRIPTION. c + + + KEY WORDS + + + c wind erosion, soil processes, surface process, layer process c + + + GLOBAL COMMON BLOCKS + + + *$noereference include 'p1werm.inc' include 'wpath.inc' include 'm1subr.inc' include 'm1sim.inc' include 'm1flag.inc' include 'w1clig.inc' *$reference c + + + ARGUMENT DECLARATIONS + + + real & bhtsmx(mnsz), & bhrwc(mnsz), & bsfom(1:mnsz), bszlyt(mnsz), & bsfsan(1:mnsz), bsfsil(1:mnsz), bsfcla(1:mnsz), & bszrgh, & bszrr, & bsfcce(1:mnsz), bsfcec(1:mnsz), * cump, dcump, bsk4d(mnsz) real bhtmx0(mnsz) real bhrwc0(mnsz) real szlyd(0:mnsz) real bszrr0, bszrh0, bseags(0:mnsz) real bseagm(mnsz),bseagmn(mnsz),bseagmx(mnsz) real bslmin(mnsz),bslmax(mnsz) real rain, snow, sprink real bwtdav, bhzirr, bszrho, bhzsmt real bszrro real bsdsblk(mnsz) c integer daysim, bslay, bm0irr c + + + LOCAL VARIABLES + + + c the 0 at the end of bhtmx0, bhrwc0, bszrr0, bszrh0 refer to c prior day values of: c max temperature, soil water content, random roughnes & ridge height c bszrro , bszrho are right-after-tillage real sf84m(mnsz),sf84sd(mnsz),scecr real tsfom, tsfcce, tsfsacl integer ldx C Retain the values of these variables for the next day c + + + LOCAL DEFINITIONS + + + c c daysim - an index for the day of simulation. c bhrwc - soil water content for today, kg/kg. c bhrwc0 - soil water content for yesterday. mass basis kg/kg. c bhtmx0 - layer maximum temperature of yesterday. in C c bhtsmx - layer maximum temperature of today in C. c bsfcce - soil fraction calcium carbonate equivalent c bsfcec - soil cation exchange capacity (cmol/kg) c bsfcla - layer fraction of clay. c bsfom - layer fraction of organic matter. c bsfsan - layer fraction of sand. c bsfsil - layer fraction of silt. c bslay - number of soil layers c bszlyt - layer thickness, mm. c bszrgh - ridge height, mm. c bszrho - original ridge height, after tillage, mm. c bszrh0 - prior day ridge height, mm c bszrr0 - prior day random roughness, mm c bszrr - random roughness height, mm c bszrro - original random roughness height, after tillage, mm c cf4m - slope coefficient for process effects on agg geom mean c cump - cumulative (rain + sprinkler + snow-melt) to current c day from day 1 or time of last tillage c dcump - c bsk4d - drying process coef. to calc. aggregate stability c scecr - ratio of clay fraction cation exchange capacity c to percent clay c se - relative aggregate stability with partial update c se1 - relative aggreage stability after SOIL update c szlyd - depth to bottom of each soil layer, mm c + + + FUNCTIONS CALLED + + + integer lentrm c + + + SUBROUTINES CALLED + + + c + + + OUTPUT FORMATS + + + c c + + + END SPECIFICATIONS + + + c check for first day or tillage if (daysim .eq. 1) then bszrr0 = - 1.0 bszrh0 = - 1.0 do ldx = 1, bslay bhtmx0(ldx) = bhtsmx(ldx) bhrwc0(ldx) = bhrwc(ldx) end do endif c if ((bszrr0.ne.bszrr).or.(bszrh0.ne.bszrgh)) then c store initial or after tillage surface roughness bszrro = bszrr bszrho = bszrgh c szlyd(0) = 0.0 szlyd(1) = bszlyt(1) c set cumulative precip to zero cump = 0.0 c store/calculate initial layer values do 10 ldx = 1, bslay c store initial water content & yesterday's temperature bhtmx0(ldx) = bhtsmx(ldx) bhrwc0(ldx) = bhrwc(ldx) c calc. mean, min, and max agg. stability c (eq. S-26, S-27, S-28, *** sb S-25,26,27) if (bsfcla(ldx) .gt. 0.5) then bseagm(ldx) = 2.73 else bseagm(ldx) = 0.83+15.7*bsfcla(ldx)-23.8*bsfcla(ldx)**2. endif bseagmn(ldx) = bseagm(ldx) - 2*(0.16)*bseagm(ldx) bseagmx(ldx) = bseagm(ldx) + 2*(0.16)*bseagm(ldx) c calc. mean and standard deviation of fraction agg. < 0.84 mm c (eq. S-42, S-43, *** sb S-37, S-38) C *** sf84m(ldx) = 0.2902 + 0.31 * bsfsan(ldx) + 0.17 * bsfsil(ldx) + C *** & 0.0033*bsfsan(ldx)/bsfcla(ldx) - 4.66*bsfom(ldx) - 0.95*bsfcce(ldx) C *** debugging fix, 1st try C *** sf84m(ldx) = 0.2902 + 0.31 * bsfsan(ldx) + 0.17 * bsfsil(ldx) + C *** & 0.0033*bsfsan(ldx)/bsfcla(ldx) - 4.66*bsfom(ldx) C eodf C *** debugging fix, 2nd try C clamping upper limits on variables to keep them from forcing sf84m negative C 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 C *** eodf sf84sd(ldx) = (0.41 - 0.22*bsfsan(ldx))*sf84m(ldx) C *** write(*,*) ' sf84m(ldx), sf84sd(ldx) ', ldx, sf84m(ldx), sf84sd(ldx) c 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' C *** write(*,*) 'bslmin(ldx),bslmax(ldx)',ldx,bslmin(ldx),bslmax(ldx) c calc. depth to bottom of each layer if ( ldx .gt. 1) szlyd(ldx) = szlyd(ldx-1) + bszlyt(ldx) c 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 C *** remove calculation of cbd; replace with original cbd calc from inpsub C *** sdbko(ldx) = bsdsblk(ldx) C ***c calc. consolidated bulk density (eq. S-52) C *** sdbko(ldx) = 1.514 + 0.25*bsfsan(ldx) - C *** * 13.*bsfsan(ldx)*bsfom(ldx) -6.0*bsfcla(ldx)* C *** * bsfom(ldx) - 0.48*bsfcla(ldx)*scecr C *** if (sdbko(ldx) .gt. 2.2) sdbko(ldx) = 2.2 C *** if (sdbko(ldx) .lt. 0.5) sdbko(ldx) = 0.5 C ***c calc. increase in consolidated bulk density with depth C ***c (note the depths change slightly with time, but current C ***c update only occurs with tillage.) C ***c (eq. S-54) C *** debugging fix C *** sdbko(ldx) = sdbko(ldx)*(0.975+ 1.931/ C *** * (1+exp(-(szlyd(ldx)-506.8)/118.5))) C *** eodf c 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 C23456789*23456789*23456789*23456789*23456789*23456789*23456789*2345 c c set stability freezing and wetting process coefficients: c c initialize rain and snow variables rain = 0.0 snow = 0.0 c Determine if precip if (awzdpt .gt. 0.0) then c Determine if precip. is rain or snow c (note HYDROLOGY may do this in future using ELS results) if (awtdav .ge. 0.0) rain = awzdpt if (awtdav .lt. 0.0) snow = awzdpt endif c if (bm0irr .eq. 1) then sprink = bhzirr else sprink = 0.0 endif c Calc. daily and cumulative (rain + sprinkler irrigation + snowmelt) dcump = rain + sprink + bhzsmt c (note: cump not used in calc., but useful as ouptput for validation) cump = cump + dcump end