c file name: soil3.for subroutine soil (daysim, bm0irr, bhzirr, bhzsmt, & bhtsmx, bhtsmn, & bhrwc, bhrwca, bhrwcw, & bsfom, bszlyt, bslay, & bsfsan, bsfsil, bsfcla, bsvroc, & bsxrgs, bszrgh, & bszrr, & bszcr, bsfcr, bsecr, bsdcr, & bsmlos, bsflos, & bsdblk, bsdagd, & bslagm, bslagn, & bs0ags, bslagx, bseags, & bbffcv, bbfscv, & bsfcce, bsfcec) 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 Edit History C 24-Feb-99 wjr changed szrgo to bszrho per LH instructions C c + + + CONTRIBUTORS to CODE + + + c Imam Elminyawi, Erik Monson, L. Hagen, Andy Hawkins, T. Zobeck 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 bhzirr, bhzsmt, & bhtsmx(mnsz), bhtsmn(mnsz), & bhrwc(mnsz), bhrwca(mnsz), bhrwcw(mnsz), & bsfom(1:mnsz), bszlyt(mnsz), & bsfsan(1:mnsz), bsfsil(1:mnsz), bsfcla(1:mnsz), & bsvroc(1:mnsz), & bsxrgs, bszrgh, & bszrr, & bszcr, bsfcr, bsecr, bsdcr, & bsmlos, bsflos, & bsdblk(0:mnsz), bsdagd(0:mnsz), & bslagm(0:mnsz), bslagn(0:mnsz), & bs0ags(0:mnsz), bslagx(0:mnsz), bseags(0:mnsz), & bbffcv, bbfscv, & bsfcce(1:mnsz), bsfcec(1:mnsz) c integer bm0irr, daysim, bslay 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 rain, snow, sprink real cump, cumpa real k4f, k4d(mnsz), k4w, c4m, c4p real cf1rg, cf2cov real arr, crr real bhtmx0(mnsz) real bhrwc0(mnsz) real hrwc0(mnsz), hrwc1(mnsz) real szlyd(0:mnsz) real bsmls0, cflos, sz real bszrr0, bszrh0 real se0, se, se1 real sdbko(mnsz), bsdbk0 real seagm(mnsz),seagmn(mnsz),seagmx(mnsz) real slmin(mnsz),slmax(mnsz) real sf84m(mnsz),sf84sd(mnsz),scecr real dcump real bszrro, bszrho integer i, j, nj integer day, mo, yr logical crpdbg C Declared here right now because it is specified in a write statement - LEW integer bm0til c + + + LOCAL DEFINITIONS + + + c c arr - regression coef. to calc. random roughness c bbffcv - biomass fraction flat cover c bbfscv - biomass fraction standing cover c bhrwc - soil water content for today, kg/kg. c bhrwc0 - soil water content for yesterday. mass basis kg/kg. c bhrwca - soil avaiable water content on mass basis kg water/kg soil. c bhrwcw - wilting point = 15 bar-grav. soil water content, kg/kg c bhtmx0 - layer maximum temperature of yesterday. in C c bhtsmn - layer minimum temperature of today in C. c bhtsmx - layer maximum temperature of today in C. c bhzirr - irrigation water applied, mm/day. c bhzsmt - snowmelt, mm/day. c bm0irr - a flag. 0 means no irrigation, 1 sprinkler, 2 furrow. c 2.5 deg). c bm0til - flag, 1 if tillage, 0 if no tillage c bs0ags - aggregate geometric standard deviation. c bsdagd - aggregate density. c bsdblk - current layer density may be different from bsdsbk. c bsdcr - crust density. Mg/m^2 c bsdbk0 - bulk density prior to update by SOIL, Mg/m^3 c bseags - agg stability, ln(J/kg). c bsfcce - soil fraction calcium carbonate equivalent c bsfcec - soil cation exchange capacity (cmol/kg) c bsecr - dry crust stability, ln(J/kg). c bsfcla - layer fraction of clay. c bsfcr - fraction of soil crust cover. m^2/m^2. c bsflos - surface cover fraction of loose material c crust area, m^2/m^2. c bsfom - layer fraction of organic matter. c bsfsan - layer fraction of sand. c bsfsil - layer fraction of silt. c bslagm - aggregate geometric mean diameter, mm. c bslagn - minimum geometric diameter for aggregates in each c layer, mm. c bslagx - maximum geometric diameter (that aggregate may reach) c for each layer, mm c bslay - number of soil layers c bsmlos - amount of loose material on crusted area, kg/m^2. c bsmls0 - prior value of bsmlos before update by SOIL, kg/m^2 c bsvroc - soil volume fraction of rock in each layer c bsxrgs - ridge spacing. we have a relation between this and bszrgh. c bszcr - crust thickness. 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 cf1rg - correction factor for ridge scale c cf4m - slope coefficient for process effects on agg geom mean c cf2cov - a plant cover correction factor for ridge height c and random roughness decrease as a result of rain. c cflos - correction factor for decease of fraction loose cover c area on crust caused by roughness c crr - regression coefficient for random roughness decrease c cump - cumulative (rain + sprinkler + snow-melt) to current c day from day 1 or time of last tillage c cumpa - apparent (rain + sprinkler + snow-metl) to current c day from time of last tillage c day - current day of simulation for output. c daysim - an index for the day of simulation. c dcump - total rain + sprinkler + snow-melt for current day. c hrwc0 - relative water content on prior day of each layer c hrwc1 - relative water content on current day of each layer c i,j - loop indexes. c k4d - drying process coef. to calc. aggregate stability c k4f - freezing process coef. to calc. aggregate stability c k4w - wetting process coef. to calc. aggregate stability c mo - current month of simulation for output. c nj - number of 5 mm water increments added in current day c rain - water added to soil as rain. c scecr - ratio of clay fraction cation exchange capacity c to percent clay c sdbko - consolidated soil bulk density by layer, Mg/m^3 c se - relative aggregate stability with partial update c se0 - relative aggregate stability prior to SOIL update c se1 - relative aggreage stability after SOIL update c szlyd - depth to bottom of each soil layer, mm c snow - water equivalent added to soil surface as snow, mm. c sprink - water added to soil as sprinkler irrigation, mm. c yr - current year of simulation for output. c + + + FUNCTIONS CALLED + + + integer lentrm c + + + SUBROUTINES CALLED + + + c caldat - input: julian day, output: day, mo, yr c + + + OUTPUT FORMATS + + + 2100 format(1x,2(i2,'/'),i4,' Subregion # ',i4) 2150 format (' tillage -', i3,' cumu. precip. -', f7.2, & ' daily precip. -',f7.2) 2200 format(1x,' -------- ridge -------- -------- crust -------', & ' --- l.e.m. ---') 2250 format(1x,' height rr cs angle thick. fract. stab.', & ' amount frac. ') 2260 format(1x,' - mm - ------ -------- - mm - ------', & ' J/m^2 Mg/m^2 ------ ') 2300 format(1x, '#',2x, 'thickness', 3x, 'agg stab', 4x, 'max dia',6x, & 'g.m.d. ', 4x, ' g.s.d.', 5x, 'density') 2310 format(4x, ' - mm -', 3x, ' - J/m^2 -', 3x, ' - mm -', 6x, & ' - mm -', 5x, ' ----- ', 4x, ' Mg/m^2 ') 2350 format(1x, 6(f7.2,2x)) 2400 format(i2, 3x, f6.1, 3x, 5(f7.2,5x)) 2450 format (2x,27('-'),' layer processes',26('-')) 2500 format (75('-')) 2550 format (1x,27('-'),' surface processes ',28('-')) 2600 format (' SOIL output for: ') c c + + + END SPECIFICATIONS + + + data crpdbg /.false./ C Initialized here right now because it is specified in a write statement - LEW bm0til = 0 c + + + INITIALIZATION SECTION + + + if (crpdbg) write (*,*) 'soil: start' c check for first day or tillage if (daysim .eq. 1) then bszrr0 = - 1.0 bszrh0 = - 1.0 endif c if ((bszrr0 .ne. bszrr) .or. (bszrh0 .ne. bszrgh)) then if (crpdbg) write (*,*) 'soil: beg 1 time init' 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 i = 1, bslay c store initial water content bhrwc0(i) = bhrwc(i) bhtmx0(i) = bhtsmx(i) c calc. mean, min, and max agg. stability (eq. S-26, S-27, S-28, *** sb S-25,26,27) if (crpdbg) write(*,*) 'bsfcla(i) ', bsfcla(i) seagm(i) = 0.83 + 15.7*bsfcla(i) - 23.8*bsfcla(i)**2. seagmn(i) = seagm(i) - 2*(0.16)*seagm(i) seagmx(i) = seagm(i) + 2*(0.16)*seagm(i) c calc. mean and standard deviation of fraction agg. < 0.84 mm c (eq. S-42, S-43) if (crpdbg) write(*,*) 'bsfsan(i) ', bsfsan(i), * ' bsfsil(i) ', bsfsil(i) if (crpdbg) write(*,*) 'bsfom(i) ', bsfom(i), * ' bsfcce(i) ', bsfcce(i) sf84m(i) = 0.2902 + 0.31 * bsfsan(i) + 0.17 * bsfsil(i) + & 0.0033*bsfsan(i)/bsfcla(i) - 4.66*bsfom(i) - 0.95*bsfcce(i) sf84sd(i) = (0.41 - 0.22*bsfsan(i))*sf84m(i) c calc. min and max values of geom. mean agg. diameter (eq. S-45, S-46) if (crpdbg) write(*,*) 'i ', i, ' bslay ', bslay if (crpdbg) write(*,*) 'sf84m(i) ', sf84m(i), * ' sf84sd(i) ', sf84sd(i) if (crpdbg) write(*,*) 'exp(3.44 - 7.21*(sf84m(i) + * 2*sf84sd(i))) ', * (3.44 - 7.21*(sf84m(i) + 2*sf84sd(i))) slmin(i) = exp(3.44 - 7.21*(sf84m(i) + 2.0 * sf84sd(i))) if (slmin(i) .lt. 0.025) slmin(i) = 0.025 slmax(i) = exp(3.44 - 7.21*(sf84m(i) - 2.0 * sf84sd(i))) if (slmax(i) .gt. 31.0) slmax(i) = 31.0 c calc. depth to bottom of each layer if ( i .gt. 1) szlyd(i) = szlyd(i-1) + bszlyt(i) c calc. ratio of clay cation exchange capacity to percent clay (eq. S-53) if (crpdbg) write (*,*) 'bsfcla(i) ', bsfcla(i) scecr = (bsfcec(i) - bsfom(i) * (142 + 0.17 * szlyd(i)))/ & (bsfcla(i) * 100.0 + 0.0001) if (scecr .lt. 0.15) scecr = 0.15 if (scecr .gt. 0.65) scecr = 0.65 c calc. consolidated bulk density (eq. S-52) sdbko(i) = 1.514 + 0.25*bsfsan(i) - 13*bsfsan(i)*bsfom(i) & -6.0*bsfcla(i)*bsfom(i) - 0.48*bsfcla(i)*scecr if (sdbko(i) .gt. 2.2) sdbko(i) = 2.2 if (sdbko(i) .lt. 0.5) sdbko(i) = 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) if (crpdbg) write (*,*) 'szlyd(i) ', szlyd(i) sdbko(i) = sdbko(i)*(0.975+ 1.931/(1+exp(-(szlyd(i)-506.8)/ & 118.5))) c c set stability drying process coefficient: k4d(i) = 0.23 +0.015*(szlyd(i-1)+(szlyd(i)-szlyd(i-1))/2.0)**0.5 10 continue endif if (crpdbg) write (*,*) 'soil: beg init 2' c c set stability freezing and wetting process coefficients: k4f = 1.4 k4w = 0.36 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 c c UPDATE SURFACE if (crpdbg) write (*,*) 'soil: beg upd surf' c skip surface processes (go to 55), if (rain+sprinkler+snowmelt=0) if (dcump .eq. 0.0) go to 55 c RIDGE SECTION: if (crpdbg) write (*,*) 'soil: beg ridge' c calculate biomass cover sheltering factor (eq. S-9 & S-10 combined) cf2cov = 1.0 - 0.6 * (bbfscv + (1.0 - bbfscv)*bbffcv) C c if ridge height is zero, skip ridge update if (bszrgh .ne. 0.0) then c calc. ridge scale factor (eq. S-8) cf1rg = (348.0 / bsxrgs)**0.3 c calculate apparent cum. precip. (eq. S-5) cumpa = ((1. - bszrgh/bszrho)/(0.034*cf1rg))**2. c update ridge height (eq. S-6) C *** debugging fix write(*,*) '********* cumpa=',cumpa,'dcump',dcump,bszrgh/bszrho c if ((cumpa + dcump * cf2cov*(1.0-bsvroc(i)))*cf1rg .ge. 0.) then write(*,*) 'soil3: ' ,(cumpa + dcump * cf2cov*(1.0-bsvroc(i))) * *cf1rg bszrgh = bszrho * (1.0 - 0.034 * sqrt(cumpa + dcump * cf2cov * & (1.0 - bsvroc(i)))*cf1rg) c else c bszrgh = bszrho c write(*,*) 'soil: debugging fix executed 2' c endif C *** end of debugging fix C c check to see that minimum bszrgh/bszrho > 0.05 if not then set c the ratio to 0.05. if ((bszrgh/bszrho) .lt. 0.05) bszrgh = 0.05 * bszrho endif C if (crpdbg) write (*,*) 'soil: beg ran rough' c RANDOM ROUGHNESS SECTION: c calc. reg. coefficients (eq. S-12, S-13) arr = 91.08 + 765.8 * bsfsil(1) crr = 0.53 + 4.66 * bsfsan(1) - 3.8 * bsfsan(1)**1.5 & - 1.22*(bsfsan(1))**0.5 c calc. apparent precip. (eq. S-11 is S-14 solved for a bare surface) cumpa = -arr * (alog(bszrr / bszrro)) * (1.0 / crr) c update random roughness (eq. S-14) C *** debugging fix c if ((cumpa + (1.0 - bsvroc(1)) * cf2cov * dcump)/arr. lt. 0.) then c bszrr = bszrro c write(*,*) 'soil: debugging fix executed 1' c else C *** end of debugging fix bszrr = bszrro*exp(-((cumpa + (1.0 - bsvroc(1))*cf2cov*dcump) & /arr)**crr) c endif if ( bszrr .lt. 2.0) bszrr = 2.0 c if (crpdbg) write (*,*) 'soil: beg crust' c CRUST SECTION: c calc. apparent precip. (eq. S-15) cumpa = (alog(1.0 - bszcr/7.6))/(0.0705 - 0.0687*bsfcla(1) & **0.146) c check for threshold precip. c ie. check to see if a H2O addition exceeding 10mm has been made C *** threshold is not noted for S-15, this test should go later if((cumpa + dcump) .lt. 10. ) go to 12 C *** if (crpdbg) write(*,*) 'soil: threshold precip' c calc. crust thickness (eq. S-16, *** sb S-15) bszcr = 7.6*(1.0 - exp(-(0.0705 + 0.0687*bsfcla(1)**0.146) & *(cumpa + dcump))) c c calc. apparent precip (eq. S-17 *** sb S-16) cumpa = -(alog(1.0 - bsfcr))/0.045 c calc. crust cover fraction (eq. S-18, *** sb S-17) bsfcr = 1.0 - exp(- 0.045*(cumpa + dcump)) c a c loose erodible material on crust c check for initial mass and set max loose mass (eq S-20, *** sb S-19) if ((bsmlos .eq. 0.0 ) .and. (bhzsmt .eq. 0.0)) then C *** bsmlos = 0.1*exp(-0.57 + 0.22*bsfsan(1)/bsfcla(1) & +7.0*bsfcce(1) - bsfom(1)) c set upper limit on loose mass (eq. S-21, *** sb S-20) if (bsmlos .gt. 3.0) bsmlos = 3.0 else c check if water is from snowmelt (eq. S-22, *** sb S-21,22) if (bhzsmt .gt. 0.0) then bsmls0 = bsmlos bsmlos = bsmlos * (1.0 - 0.1 * bhzsmt) write(*,*) 'soil: 2' if ((bsmlos/bsmls0) .lt. 0.1) bsmlos = bsmls0 * 0.1 else bsmlos = bsmlos * (1.0 - 0.0053 * dcump) endif endif c c fraction cover of loose erodible material (eq. S-24, S-25, sb S-23,24) sz = amax1(bszrr, bszrgh) write(*,*) 'soil3: bsmlos ', bsmlos cflos = sqrt(bsmlos)/(0.24*sz) if (cflos .gt. 1.0) cflos = 1.0 bsflos = (1.0 - exp(-3.5*bsmlos**1.5))*cflos c 12 continue c skip layer update on first simulation day 55 if (daysim .lt. 2) go to 800 c c + + + UPDATE LAYERS: + + + if (crpdbg) write (*,*) 'soil: beg upd lay' i = 1 do while ((szlyd(i) .le. 300.) .or. (i .eq. 1)) c c check for water content change or dry soil - then no changes if (bhrwc(i) .eq. bhrwc0(i)) go to 90 if ((bhrwc(i) .le. bhrwcw(i)) .and. & (bhrwc0(i) .le. bhrwcw(i))) go to 90 c c AGGREGATE STABILITY SECTION: c if (crpdbg) write (*,*) 'soil: beg agg sta' c calc. relative agg stability and water content for prior day se0 = (bseags(i) - seagmn(i))/(seagmx(i) - seagmn(i)) hrwc0(i) = (bhrwc0(i) - bhrwcw(i))/bhrwca(i) c c calc. relative water content for current day hrwc1(i) = (bhrwc(i) - bhrwcw(i)) / bhrwca(i) if (hrwc1(i) .lt. 0.0) hrwc1(i) = 0.0 c c check soil temp to determine process c c check if unfrozen then wet or dry if ((bhtmx0(i) .gt. 0.0) .and. (bhtsmn(i) .gt. 0.0)) go to 70 c c check for freeze after wet or dry if ((bhtsmx(i) .lt. 0.0) .and. (bhtsmn(i) .lt. 0.0) & .and. (bhtmx0(i) .gt. 0.0)) go to 70 c c check for thaw process alone if ((bhtmx0(i) .lt. 0.0) .and. (bhtsmx(i) .ge. 0.0) & .and. (bhtsmn(i) .ge. 0.0)) go to 60 c c check for freeze/thaw if ((bhtsmx(i) .gt. 0.0) .and. (bhtsmn(i) .lt. 0.0) & .and. (bhtmx0(i) .gt. 0.0)) then goto 50 c c check for drying while frozen elseif ((bhtsmn(i) .lt. 0.0) .and. (bhtmx0(i) .lt. 0) & .and. (bhrwc(i) .lt. bhrwc0(i))) then c update stability with drying while frozen se1 = se0 - (2.7 - se0) * k4f * (hrwc0(i) - hrwc1(i))/ & (2.5 - k4f * hrwc1(i)) go to 80 else go to 80 endif c freeze process with prior day water content 50 se = se0 * (1.0 - k4w * k4f * hrwc0(i))/(1.0 - k4w * hrwc0(i)) se0 = se + ((2.7 - se)/2.5)*k4f*hrwc0(i) c c thaw process with prior day water content 60 continue if (hrwc0(i) .gt. 1.0) then c soil puddling process when saturated se0 = 0.74/hrwc0(i) else se = se0 - (2.7 - se0)*k4f*hrwc0(i)/(2.5 - k4f * hrwc0(i)) se0 = se*(1.0 - k4d(i)*hrwc0(i))/(1.0 - k4d(i)*k4f* & hrwc0(i)) endif go to 70 c c check if wetting process 70 if (bhrwc(i) .gt. bhrwc0(i)) then c wetting process se1 = se0*(1.0 - k4w*hrwc1(i))/(1.0 - k4w*hrwc0(i)) else c drying process se1 = se0*(1.0 - k4d(i)*hrwc1(i))/(1.0 - k4d(i)*hrwc0(i)) endif c c check if freeze after w/d if ((bhtmx0(i) .gt. 0.0) .and. (bhtsmx(i) .lt. 0.0) & .and. (bhtsmn(i) .lt. 0.0)) then c freeze process with today water content se = se0*(1.0 - k4w*k4f*hrwc1(i))/(1.0 - k4w* & hrwc1(i)) se1 = se + ((2.7 - se)/2.5)*k4f*hrwc1(i) endif c c calc. new agg. stability (S-28) 80 bseags(i) = se1*(seagmx(i) - seagmn(i)) + seagmn(i) c c calc. new crust stability if (i .eq. 1) bsecr = bseags(i) c if (crpdbg) write (*,*) 'soil: beg asd' c ASD SECTION: c calc. slope parameter using prior geometric mean dia. write(*,*) ' soil3: slope param ', * (bslagm(i) - slmin(i))/(slmax(i) - slmin(i)), * bslagm(i), slmin(i), slmax(i) c4m = sqrt((bslagm(i) - slmin(i))/(slmax(i) - slmin(i))) & * (1.0/se0) if (c4m .lt. 0.5) c4m = 0.5 if (c4m .gt. 2.0) c4m = 2.0 c c calc. new geom. mean dia. bslagm(i) = slmin(i) + (slmax(i) - slmin(i))*(c4m*se1)**2. c restrict upper size if not frozen if ((bhtsmx(i) .ge. 0.0) .and. (bslagm(i) .gt. slmax(i))) then bslagm(i) = slmax(i) endif c c calc geom. standard deviation (eq. S-49) write (*,*) 'soil3: i, bslagm(i) ', i, bslagm(i) bs0ags(i) = 1.0 / (0.0203 + 0.00193 *bslagm(i) + & 0.074 / sqrt(bslagm(i))) c c calc. max. diameter for asd c4p = 1.52 * bslagm(i)**(-0.449) bslagx(i) = bslagm(i) * bs0ags(i)**c4p c if (crpdbg) write (*,*) 'soil: beg den' c DENSITY SECTION: c update bulk density c c store initial value of layer density bsdbk0 = bsdblk(i) c if ((dcump .gt. 0) .and. (bsdblk(i) .lt. 0.97 * sdbko(i))) then c update for water additions in 5 mm increments nj =nint(dcump/5.) do 85 j = 1,nj bsdblk(i) = bsdblk(i) + 0.75 * (1-(bsdblk(i)/ & (0.97 * sdbko(i))))**1.5 85 continue else c daily update density (except precip days) for other forces bsdblk(i) = bsdblk(i) + 0.01*(sdbko(i) - bsdblk(i)) endif c c update layer thickness bszlyt(i) = bszlyt(i)*bsdbk0/bsdblk(i) c c update aggregate density bsdagd(i) = sdbko(i) c 90 continue i = i + 1 c temporary fix to keep i from incrementing above the number of c soil layers - should be fixed permanently - jt 4/19/99 if (i .gt. bslay) go to 91 end do c c update crust density (S-58) 91 bsdcr = 0.576 + 0.603 * sdbko(1) c c Assign today's values to 'yesterday storage' 800 do 801 i = 1, bslay bhtmx0(i) = bhtsmx(i) bhrwc0(i) = bhrwc(i) 801 continue bszrr0 = bszrr bszrh0 = bszrgh c if (crpdbg) write (*,*) 'soil: beg out' c + + + OUTPUT SECTION + + + C if ((am0sfl .eq. 1)) then open(16,file= rootp(1:lentrm(rootp)) // 'soil.out', & access='sequential', & status='unknown') call caldat(am0jd,day,mo,yr) write(16,2600) write(16,2100) day, mo, yr, am0csr write(16,2150) bm0til, cump, dcump write(16,2550) write(16,2200) write(16,2250) write(16,2260) write(16,2350) bszrgh, bszcr, bsfcr, & bsecr, bsmlos, bsflos write(16,2450) write(16,2300) write(16,2310) c output new values by layer to the soil output file. i = 1 do while ((szlyd(i) .le. 300) .or. (i .eq. 1)) write (16,2400) i, bszlyt(i), bseags(i), bslagx(i), & bslagm(i), bs0ags(i), bsdblk(i) i = i + 1 end do write(16,2500) endif return end