c file name: soil.for subroutine soil (daysim, bm0irr, bhzirr, bhzsmt, & bhtsmx, bhtsmn, & bhrwc, bhrwca, bhrwcw, & bsfom, bszlyt, slay, & bsfsan, bsfsil, bsfcla, & bsxrgs, bszrgh, bslrrc, & bszcr, bsfcr, bsecr, bsdcr, & bsmlos, bsflos, & bsdsbk, & bsdblk, bsdagd, & bslagm, bslagn, & bs0ags, bslagx, bseags, & bcftcv, icover) c + + + PURPOSE + + + c soil submodel for the Wind Erosion Research Model. c from compilation by L. HAGEN, and I. ELMINYAWI 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 + + + AUTHOR + + + c Imam Elminyawi 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 bm0irr, bhzirr, bhzsmt, & bhtsmx(mnsz), bhtsmn(mnsz), & bhrwc(mnsz), bhrwca(mnsz), bhrwcw(mnsz), & bsfom(0:mnsz), bszlyt(mnsz), & bsfsan(0:mnsz), bsfsil(0:mnsz), bsfcla(0:mnsz), & bsxrgs, bszrgh, & bslrrc, bs0rrk, & 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), & bcftcv integer daysim, slay, icover c + + + LOCAL VARIABLES + + + character proces*10 real k3los, k4ft , k4wd, k4fd, k5wd, k6, k6wd, k6ft, k6fd, k7wd real k1los, k2los, k8bk real d(3) , f(3) , act(3), bct(3), zmax(3) real dmax1 c bhtmx0 and bhrwcy refer to yesterday's values for c max temperature, and soil water content respectively. c b0, bsagm0, bsdsb0, bsagx0 are right-after-tillage values. real bhtmx0(mnsz) real bhrwcy(mnsz) real hrwcp(mnsz) real bszrh0, b0 real bszlym(mnsz) real agsave(mnsz) real bsagm0(0:mnsz) real bsagx0(0:mnsz) real bsdsbk(mnsz), bsdsbk1(mnsz), bsdsb0(mnsz) c-------------------------------------------------- c local variables integer i, l integer day, mo, yr real a, abct, agsmn, agsmx, alosem real b, bmin, bsand1, bslagx3 real c, c1, c2, cf4, covcff, covcf2, cump, cumsr real dt, dt1, depth, dsqr, dcump, dcumsr, dx real rain, rgcf1 c real sdmn, snow, sprink, sum0, sum1, sum2, sum22 real sdmn, snow, sprink, sum0, sum1, sum2 real w c + + + LOCAL DEFINITIONS + + + c a - a factor depending on silt fraction in the surface c layer, used in random roughness c abct - see eq. S30 c act - eq. S23. c agsave - soil agg. ave. stability. (internal variable). c agsmn - minimum stability of aggregates. c agsmx - maximum stability of aggregates. c alosem - max loose erodible material. c b - c b0 - initial random roughness scale parameter. c bcftcv - plant cover fraction. c bct - eq. S24. c bhrwc - soil water content for today, kg/kg. c bhrwcy - 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 bmin - minimum random roughness scale parameter (will be assigned c 2.5 deg). c bs0ags - aggregate geometric standard deviation. c bs0rrk - random roughness shape parameter. not changed at this stage. c bsagm0 - geometric mean diameter after last tillage. c bsagx0 - agg. max. diameter after last tillage. c bsand1 - fraction of sand less than 1mm is set to 0.1. c bsdagd - aggregate density. c bsdblk - current layer density may be different from bsdsbk. c bsdcr - crust density. Mg/m^2 c bsdsb0 - settled bulk density at 300 mm based on sand, clay, and c om. per layer. c bsdsbk - settled bulk density, modified from bsdsb0 by the layer's c depth. c bseags - agg stability, ln(J/kg). c bsecr - dry crust stability. c bsfcla - layer fraction of clay. c bsfcr - fraction of soil crust cover. m^2/m^2. c bsflos - surface fraction of loose material. 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 mean diameter for aggregates in each c layer mm. c bslagx - maximum diameter (that aggregate may reach) for each layer. c bslagx3 - (2*gmd*gsd) used in calculation of surface aggr. size. dist. c bslrrc - random roughness scale parameter (current one). c bsmlos - amount of loose material mg/m^2. c bsxrgs - ridge spacing. we have a relation between this and bszrgh. c bszcr - crust thickness. c bszlym - depth to the mid point of a layer thickness. c bszlyt - layer thickness, mm. c bszrgh - ridge height. c bszrh0 - original ridge height. c c - a factor depending on sand fraction in the surface c layer, used in random roughness. c c1 - equivalent to previous crust thickness. c c2 - equivalent to current crust thickness. c cf4 - a correction factor to the crust thickness as a c result of surface clay content. c covcf2 - a plant cover correction factor for ridge height c decrease as a result of rain. c covcff - a plant cover correction factor for random roughness c decrease as a result of rain. c cump - equivalent (rain + sprinkler + snow-melt) for current day. c cumsr - total equivalent (rain + sprinkler) for this day. 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 the day. c dcumsr - total rain + sprinkler for the day. c depth - depth to bottom of current soil layer - used to determine c if depth > 300 mm for settled bulk density calculation. c dsqr - square of aggregate diameter - used to get crust thickness. c dt - the sum of the inverse of 3 particle diameter categories c on the soil surface. c dt1 - inverse of dt. c dx - log of ratio of crust thickness of two consecutive days. c hrwcp - difference between bhrwc and wilting point bhrwcw. c i - loop index. c icover - an indicator of where the plant cover is. c 1,2,3 is for ridge, forrow, or random cover respectively. c any other integer is for no cover. c l - loop index on soil layers. c mo - current month of simulation for output. c proces - character variable which signifies the type of process. c rain - water added to soil as rain. c rgcf1 - correction factor for ridge scale, eq. (S8). c sdmn - minimum soil bulk density, eq, (S76). c slay - no. of soil layers c slay + 1 - no. of (layers + surface). c snow - water added to soil as snow melt. c sprink - water added to soil as sprinkler irrigation. c sum0 - sum over max crust thickness contributions from 3 c particle diameters. c sum1 - sum needed during calculation of crust thickness. c sum2 - " " " " " " " c sum22 - " " " " " " " c w - ratio of soil water content to wilting point. c yr - current year of simulation for output. c + + + SUBROUTINES CALLED + + + c mpoint - claculate depth to midpoint of soil layers given c layers thickness. c aggs - calculate the aggregate stability for the different c processes from eq. S(54), S(53). c aggss - calculates aggregate stability increase due to pressure c and wet dry cycle (eq. S(58), S(59)). c asd - calculate agg size distribution eq. S(67), and S(62) c (deaggregation process). c asdd - aggregate size distribution eq. S(68), and S(69) c (aggregation process) due to wet/dry c + + + FUNCTION DECLARATION + + + real setbd c + + + FUNCTIONS DEFINITION + + + c setbd - estimates settled soil bulk density from intrinsic c properties. see rawls (1983) soil science 135, 123-125. integer lentrm c + + + INPUT FORMATS + + + c + + + OUTPUT FORMATS + + + 2100 format(1x,2(i2,'/'),i4,' Subregion # ',i4) 2150 format(1x,' ridge crust ', & ' loose erod. matl. ') 2200 format(1x,' --------------------- ---------------------------', & ' -----------------') 2250 format(1x,' height rough. scale thick. cov. fract. stabil.', & ' amount frac. ') 2260 format(1x,' - mm - ------------ - mm - ---------- J/m^2', & ' Mg/m^2 ------ ') 2300 format(2x, 'layer #', 3x, 'agg stability', 3x, 'max dia', 4x, & 'g.m.d. ', 5x, ' g.s.d.', 5x, 'density') 2310 format(2x, ' ', 3x, ' -- J/m^2 --', 3x, ' - mm -', 4x, & ' - mm -', 6x, ' ----- ', 4x, ' Mg/m^2 ') 2350 format(1x, 7(f7.2,4x)) 2400 format(2x,i4,8x, 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 + + + END SPECIFICATIONS + + + c set initial values proces = ' ' c need to test for tillage here also ??? if (daysim .eq. 880) then daysim = 880 endif if (daysim .eq. 1) then b0 = bslrrc bszrh0 = bszrgh do 10 l = 1, slay bhtmx0(l) = bhtsmx(l) bhrwcy(l) = bhrwc(l) 10 continue bsdblk(0) = bsdblk(1) bsdagd(0) = bsdagd(1) bseags(0) = bseags(1) bslagx(0) = bslagx(1) bslagm(0) = bslagm(1) bslagn(0) = bslagn(1) bs0ags(0) = bs0ags(1) c bsfcla(0) = bsfcla(1) c bsfsil(0) = bsfsil(1) c bsfsan(0) = bsfsan(1) do 20 l = 0, slay bsagx0(l) = bslagx(l) bsagm0(l) = bslagm(l) 20 continue end if c bsand1 fraction of sand less than 1mm is set to 0.1 and could be changed. c bmin, k3los, k8bk: minimum values for the b random roughnss parameter, c a coeff for loose erodible material, and bulk density respectively. c dcumsr is rain and sprinkler, c dcump rain, sprinkler irrigation, and snow-melt. bsand1 = 0.1 bmin = 2.5 k3los = 5.0 k8bk = 0.6 snow = 0.0 rain = 0.0 sprink = 0.0 if (( awzdpt .gt. 0.0) .and. (awtdav .gt. 0.)) then rain = awzdpt else if ( awzdpt .gt. 0.0 .and. awtdav .lt. 0.0 ) snow = awzdpt endif if ( bm0irr .eq. 1.0 ) sprink = bhzirr dcumsr = sprink + rain dcump = dcumsr + bhzsmt c skip surface processes (go to 55), if (precip + snowmelt) = 0.0 if (dcump .eq. 0.0) go to 55 c eq. S(17), S(18) parameters for random roughness change, then parameters for c loose erodible material eq. S(45), S(43), where alosem is max loose material. a = 91.0836 + 765.77 * bsfsil(1) c = 0.5284 + 4.66 * bsfsan(1) - 0.48 * (bsfsan(1))**1.05 & -0.1537 * (bsfsan(1))**0.05 k1los = 1.0 + 6.0 * bsfcla(1) k2los = bsand1 / bsfsan(1) alosem = 2.02 * k2los * exp( -(bsfsan(1) - 1.39)**2 / 0.1685 ) c we skip updating ridges, if ridge (height) does not exist at this day. c rgcf1 is a correction factor for ridge scale eq. S(8) if ( bszrgh .eq. 0.0 ) go to 8 rgcf1 = (348.0 / bsxrgs)**0.3 c ridge and furrow dike height as affected by cumulative depth of snow- c melt, rainfall and sprinkler irrigation water. icover refers to one of c three kinds of covers. (on the ridge, on the forrow, or random) eq S(9)-(11) if (bcftcv .gt. 0.0) then if (icover - 2) 2,3,4 2 covcf2 = 1.0 - 0.8 * bcftcv go to 5 3 covcf2 = 1.0 - 0.8 * (bcftcv - 0.5) go to 5 4 covcf2 = 1.0 - 0.4 * bcftcv 5 continue else covcf2 = 1.0 endif c bszrgh, bszrh0 are the current and the started-with ridge hight. eq. S(5),S(7) if (bszrgh .lt. bszrh0) then cump = ( (bszrgh /bszrh0 /covcf2 / rgcf1 - 1.0) / 0.055)**2 else cump = 0.0 bszrh0 = bszrgh endif cump = cump + dcump bszrgh = covcf2 * bszrh0 * ( 1.0 - 0.055 * sqrt(cump) ) c this part is for the random roughness as affected by total water. cump. c we only update b ________________ c eq. S(19) for plant cover effect, then eq. S(16), and S(15) in this order. 8 covcff = 1.0 - 0.8 * bcftcv if (bslrrc .lt. b0) then cump = a*( -alog( ( bslrrc - bmin )/ ( b0 - bmin )/ covcff ) ) & **( 1.0 / 0.8121 / c ) else cump = 0.0 b0 = bslrrc endif cump = cump + dcump bslrrc = covcff * (b0 - bmin) * exp(-(cump / a)**(0.8121 * c) ) & + bmin c next segment is for the crust thickness, current value is (bszcr). S(25)-S(29) c1 = bszcr if (bszcr .gt. 20.0) go to 32 d(1) = bslagm(0) / bs0ags(0) d(2) = bslagm(0) d(3) = bslagm(0) * bs0ags(0) dt = 1.0 / d(1) + 1.0 / d(2) + 1.0 / d(3) dt1 = 1.0 / dt cf4 = 1.6 * (1.0 - exp( -10.0 * bsfcla(1) ) ) c eq. S(20) for different components of crust due to the 3 diameters category. sum0 = 0.0 do 25 i = 1, 3 zmax(i) = ( 1.25 + 0.47 * d(i) ) * cf4 sum0 = sum0 + zmax(i) 25 continue if (bszcr .lt. sum0) then sum1 = 0.0 sum2 = 0.0 do 31 i = 1,3 f(i) = dt1 / d(i) dsqr = d(i) * d(i) c eq. S(23), S(24) act(i) = -1.41+ 0.822 / d(i) + 1.633 / dsqr if (act(i) .lt. -0.02) go to 13 act(i)= -0.02 13 bct(i) = exp( -0.383 - 0.00012 / dsqr - 3.06 * exp(-d(i)) )/ cf4 sum1 = sum1 + act(i) * zmax(i) sum2 = sum2 + bct(i) * zmax(i) c sum22 = sum2 * alog(cump) c if ( (sum1 + sum2) .gt. 0.0 ) go to 31 c sum1 = 0.001 c sum2 = 0.001 31 continue c fract = f(3) / 2.0 + f(1) / 2.0 + f(2) c cump = exp( ( 3.0 * bszcr / dt - sum1) / sum2 ) + dcump cump = exp( ( bszcr / dt - sum1) / sum2 ) + dcump bszcr = 0.0 do 34 i=1,3 abct = act(i) + bct(i) * alog(cump) if (abct .lt. 0.0) abct = 0.0 if (abct .gt. 1.0) abct = 1.0 bszcr = bszcr + abct * zmax(i) 34 continue else cump = 0.0 bszcr = sum0 endif 32 c2 = bszcr c this segment is for crust cover fraction, we took the original value =0.0 c eq. S(33) ____________________ dx = alog(c2) - alog(c1) c bsfcr = 1.0 - dx * exp (fract * bszcr) c bsfcr = (f(1)-f(3))/(alog(d(3))-alog(d(1)))*dx+bsfcr bsfcr = (f(1)-f(3))/(alog(d(3))-alog(d(1)))*dx c loose erodible material eq. S(43), S(46) if ( dcumsr .eq. 0.0 .or. bsmlos .ge. alosem ) go to 33 cumsr = - k1los * alog( 1.0 - bsmlos / alosem ) + dcumsr bsmlos = alosem * ( 1.0 - exp(- cumsr / k1los) ) 33 continue c loose erodible material cover. bhzsmt may consolidate loose erodible material. c ______________________________ eq. S(47), S(48) if ( bsmlos .ge. alosem ) bsmlos = alosem if (bhzsmt .eq. 0.0) goto 35 bhzsmt = -k3los * alog( bsmlos / alosem ) + bhzsmt bsmlos = alosem * exp(- bhzsmt / k3los ) c we now calulate eq. S(49) after calculating the gain in the loose, c erodible material due to sprinkler and rainfall, and loss in the c amount of loose erodible material due to some consolidation due to bhzsmt. bsflos = 1.5 * bsmlos * ( 1.0 - 0.6 * ( bsmlos / alosem ) ) c aggregate size distribution of the zero layer eq. S(36)....S(41) c _____________________________________________ 35 if (bszcr .gt. bslagm(0)) then dmax1 = bslagm(0) * exp(0.4725) else dmax1 = bslagm(0) * bslagm(0) / bszcr endif bslagx3 = 2.0 * bslagm(0) * bs0ags(0) bsagx0(0)= sqrt(dmax1* bslagx3 ) bslagm(0)= sqrt( c2 * bsagx0(0) ) bs0ags(0)= 0.5 * bslagm(0) / c2 c ********************************************************************* c ********************************************************************* c layer processes from 1 - slay. we skip this whole part if it is the 1st day. c _____ 55 if ( daysim .lt. 2) go to 6800 c we calculate the depth to each layer midpoint from given thicknesses. call mpoint (bszlyt, bszlym, slay) c 1st eq. is a function by Skidmore, for the settled bulk density at 300 mm. c then bulk density eq. S(73), for layers above 300 mm depth. this is c not the current bulk density, but what it should be depth = 0. do 4000 l = 1, slay if (bslagm(l) .gt. bsagm0(l)) bsagm0(l) = bslagm(l) bsdsb0(l) = setbd ( bsfsan(l), bsfcla(l), bsfom(l)) depth = depth + bszlyt(l) if (depth .le. 300.) then bsdsbk(l) = bsdsb0(l) * ( 0.72 + 0.00092 * bszlym(l)) c bsdsbk(l) = bsdsb0(l) * ( 0.72 + 0.00092 * depth) else end if c coefficient for aggregate size distribution. for top layer in case of rain/ c sprinkler irrigation is greater than 0 eq. S(63) if ( dcumsr .eq. 0.0 .or. l .gt. 1 ) go to 75 c Fix the divide by zero if (bsagm0(1) .lt. 0.0001) then bsagm0(1) = 0.1 endif k6 =( ( bsagm0(1) - 1.5) / 0.004 / bsdsbk(1) / bsagm0(1) & + 50.0 * bsfcla(1) ) / ( 1.01 - bcftcv ) call asd (k6, bsagm0(1), bslagm(1), bslagn(1)) 75 continue c dry aggregate stability. bseags. eq. S(50), S(51), S(52) respectively. c _______________________ agsave(l) = 0.83 + 15.7* bsfcla(l) - 23.8 * bsfcla(l) * bsfcla(l) agsmx = 1.32 * agsave(l) agsmn = 0.68 * agsave(l) hrwcp(l) = bhrwc(l) - bhrwcw(l) c start identifying the different cycles freeze/thaw, freeze/dry, wet/dry. c next is freeze/thaw cycle, coefficient for agg. stability eq. S(55). if( bhtsmx(l) .gt. 0.0 .and. bhtsmn(l) .lt. 0.0 ) then w = bhrwc(l) / bhrwcw(l) if (w .lt. 1.5) then k4ft = 100.0 else k4ft = 1.0 / ( 0.0346 * w - 0.0417 ) endif call aggs (k4ft, bseags(l), agsmn, agsmx) c eq. S(65) for agg size dist. coefficient. k6ft=( 20.0+50.0*bsfcla(l) ) * bhrwca(l)/hrwcp(l) call asd ( k6ft, bsagm0(l), bslagm(l), bslagn(l)) if (l .eq. 2) proces = 'ft' else if (bhtsmx(l) .lt. 0.0 .and. bhtmx0(l) .lt. 0.0) then if ( bhrwcy(l) .le. bhrwc(l)) go to 95 c freeze/dry coefficient S(57), S(66) for stability and asd respectively. k4fd = bhrwca(l) / hrwcp(l) call aggs ( k4fd, bseags(l), agsmn, agsmx ) k6fd =(2.0 + 5.0 * bsfcla(l)) * bhrwca(l) / hrwcp(l) call asd ( k6fd, bsagm0(l), bslagm(l), bslagn(l)) if (l .eq. 2) proces = 'fd' 95 continue else if ( bhrwcy(l) .lt. bhrwc(l) ) then c wet/dry coefficient eq. S(56) for stability. k4wd = 2.0 * bhrwca(l) / hrwcp(l) + 4.0 * bsfcla(l) & + ( bszlym(l) / 5.0 ) **3.6 call aggs ( k4wd, bseags(l), agsmn, agsmx ) if (l .eq. 2) proces = 'wd' c increase aggregate stability due to saturation/pressure, eq. S(60), and S(61) if ( hrwcp(l) .gt. bhrwca(l) ) then k5wd = bhrwca(l) / hrwcp(l) else k5wd = bhrwca(l) / hrwcp(l) + bszlym(slay) / & bszlym(l) endif c this is a different subroutine for agg stability than aggs. call aggss (k5wd, bseags(l), agsmn, agsmx) c eq. S(65) wet/dry cycle again that works better at lower layers. eq. S(64) k6wd = ( 20.0 + 50.0 * bsfcla(l) & + ( bszlym(l)/5.0)**3.6 )* bhrwca(l)/hrwcp(l) call asd (k6wd, bsagm0(l), bslagm(l), bslagn(l)) c eq. S(70) for gmd for aggregation effect on geometric mean diameter k7wd = bhrwca(l) / hrwcp(l) + bszlym(slay)/bszlym(l) c a new subroutine different from asd. call asdd ( k7wd, bslagm(l), bsagx0(l) ) else endif endif endif 500 continue c how bs0ags and bslagx are linked to the calculated gmd. eq. S(71), S(72). c these are geometric standard deviation, and maximum diameter respectively. c the 0 at the end of these names refer to the value right after tillage. if (bslagm(l) .lt. 0.0) then c write(*,*) 'l=',l,' bslagm= ',bslagm(l),' resetting blsagm >0 c & ' bslagm(l)=1.0 end if bs0ags(l) = ( bslagm(l) / 0.3 )**0.73 c Fix the divide by zero if (bsagm0(l) .lt. 0.001) then bsagm0(l) = 0.1 endif bslagx(l) = bslagm(l) * bsagx0(l) / bsagm0(l) c if current bulk density < bsdsbk, due to tillage for example, then wetting/ c increase in water content will bring it back to its previous value. c sdmn the minimum soil bulk density eq. S(76). if ( bsdblk(l) .ge. bsdsbk(l) ) go to 3999 if ( bhrwc(l) .le. bhrwcy(l) ) go to 3999 bsdsbk1(l) = bsdblk(l) sdmn = bsdsbk(l) - 0.6 bsdblk(l)= 0.6 * (1.0 - exp(-bhrwc(l) / bhrwca(l) /k8bk)) + sdmn c aggregate density, eq. S(77). then new layer thickness if layer denllsity changes. c bsdagd(l) = bsdsb0(l) bszlyt(l) = bszlyt(l) * bsdsbk1(l) / bsdblk(l) 3999 bsdagd(l) = bsdsb0(l) 4000 continue c crust density, and stability from eq. S(78), and S(42) respectively. bsdcr = 0.576 + 0.603 * bsdsb0(1) bsecr = -1.0888 + 0.7975 * bseags(1) c we assign today's values as yesterday's values to be ready for new day. 6800 do 6801 l = 1, slay bhtmx0(l) = bhtsmx(l) bhrwcy(l) = bhrwc(l) bsagm0(l) = bslagm(l) bszrh0 = bszrgh C b0 = b b0 = bslrrc 6801 continue c + + + OUTPUT SECTION + + + if ((am0sfl .eq. 1)) then open(16,file= outp(1:lentrm(outp)) // 'soil.out', & access='sequential', status='unknown') call caldat(am0jd,day,mo,yr) write(16,2600) write(16,2100) day, mo, yr, am0csr write(16,2550) write(16,2150) write(16,2200) write(16,2250) write(16,2260) write(16,2350) bszrgh, bslrrc, bszcr, bsfcr, bsecr, bsmlos, & bsflos write(16,2450) write(16,2300) write(16,2310) do 6900 l=1, slay write (16,2400) l, bseags(l), bslagx(l), bslagm(l), & bs0ags(l), bsdblk(l) 6900 continue write(16,2500) end if c we reassign today's values in temprature and water content as yesterdays 9003 do 9004 l = 1, slay bhrwcy(l) = bhrwc(l) bhtmx0(l) = bhtsmx(l) 9004 continue return end