c file: 'bookeep.for' subroutine bookeep (cd, cm, cy, xmax, xmin, ymax, ymin) c + + + PURPOSE + + + c This program keeps track of output statistics of the Wind c Erosion Prediction System. c Its purpose is to sum and average output by time and c management operation periods. c author: John Tatarko c version: 1 c + + + KEY WORDS + + + c wind, erosion, hydrology, tillage, soil, crop, decomposition c management, output, report c + + + GLOBAL COMMON BLOCKS + + + *$noereference include 'p1werm.inc' include 'p1unconv.inc' include 'm1subr.inc' include 'm1sim.inc' include 'm1geo.inc' include 'm1flag.inc' include 'm1dbug.inc' include 's1layr.inc' include 's1surf.inc' include 's1phys.inc' include 's1agg.inc' include 's1dbh.inc' include 's1dbc.inc' include 's1layd.inc' include 's1sgeo.inc' include 'c1info.inc' include 'c1gen.inc' include 'c1db1.inc' include 'c1db2.inc' include 'c1db3.inc' include 'c1db4.inc' include 'c1db5.inc' include 'c1glob.inc' include 'c1geom.inc' include 'd1mass.inc' include 'd1geom.inc' include 'd1glob.inc' include 'b1geom.inc' include 'b1glob.inc' include 'b1mass.inc' include 'w1clig.inc' include 'w1wind.inc' include 'w1pavg.inc' include 'h1et.inc' include 'h1hydro.inc' include 'h1scs.inc' include 'h1db1.inc' include 'h1temp.inc' c + + + LOCAL COMMON BLOCKS + + + include 'main/main.inc' include 'main/output.inc' include 'manage/man.inc' include 'manage/oper.inc' include 'manage/mproc.inc' include 'erosion/m2geo.inc' include 'erosion/e2erod.inc' *$reference c + + + LOCAL VARIABLES + + + integer cd, cm, cy, i i, j, k, i ndiy, ngdpt, rptdat, i xmax, xmin, ymax, ymin real avgetot, avgets, avget10, r cumpre, cumprery, r gx(0:30), gy(0:30), r maxegt, minegt, r sumegtot, sumegts, sumegt10, r tmp c + + + LOCAL DEFINITIONS + + + c am0*fl - These are switches for production of submodel c output, where the asterisk represents the first c letter of the submodel name. c am0gdf - c am0ifl - This variable is an initialization flag which is c set to .false. after the first simulation day. c avgetot - This varaible holds the average total erosion for the c user selected period. c avgets - This varaible holds the average erosion suspension for c the user selected period. c avget10 - This varaible holds the average PM-10 erosion for the c user selected period. c brcd - effective standing biomass c cd - The current day of simulation month. c cm - The current month of simulation year. c cy - The current year of simulation run. c daysim - This variable holds the total current days of simulation. c header - Dummy variable to read in character values which c are not used. c i - This variable is a counter for simulation loops. c id,im,iy - The initial day, month, and year of simulation. c ijday - This variable contains the initial julian day of c the simulation run. c ioc - This variable is set in read statement to zero if there c is no error and -1 if end-of-file c isr - This variable holds the subregion index. c l - This variable is an index on soil layers. c ld,lm,ly - The last day, month, and year of simulation. c ljday - This variable contains the last julian day of c the simulation run. c maxper - The maximum number of years in a rotation for all c subregions. c ndiy - The number of days in the year. c nsubr - This variable holds the total number of subregions. c ngdpt - This variable holds the total number of grig points in an c accounting region. c period - The number of years in a management rotation. This c variable is defined in (/inc/manage/man.inc) c subflg - This logical variable is used to read header information c in the subdaily wind file (if .true., read header). c usrid - This character variable is an identification string c to aid the user in identifying the simulation run. c usrloc - This character variable holds a location c description of the simulation site. c usrnam - This character variable holds the user name. c wflag - This logical variable controls output of warning c messages for mismatch of WINDGEN and simulation dates. c wcflag - This logical variable controls output of warning c messages for mismatch of CLIGEN and WINDGEN dates. c + + + SUBROUTINES CALLED + + + c + + + FUNCTIONS CALLED + + + integer lstday integer dayear c + + + DATA INITIALIZATIONS + + + if (am0oif .eqv. .true.) then if ((cd .eq. 1) .and. (cm .eq. 1)) then lopday = 31 lopmon = 12 c lopyr = maxper end if do 10 i = 1,32 do 20 j = 1,13 nyrs(j) = 1 do 30 k = 1,6 outyrs(k) = 0 opnam(i,j,k) = ' ' opimp(i,j,k) = ' ' opdir(i,j,k) = 0.0 ndays(i,j,k) = 1 tloss(i,j,k) = 0.0 std(i,j,k) = 0.0 sloss(i,j,k) = 0.0 pm10loss(i,j,k) = 0.0 cs1(i,j,k)= 0.0 cs2(i,j,k)= 0.0 cs3(i,j,k)= 0.0 cs4(i,j,k)= 0.0 ss1(i,j,k)= 0.0 ss2(i,j,k)= 0.0 ss3(i,j,k)= 0.0 ss4(i,j,k)= 0.0 pm1(i,j,k)= 0.0 pm2(i,j,k)= 0.0 pm3(i,j,k)= 0.0 pm4(i,j,k)= 0.0 precip(i,j,k) = 0.0 winde(i,j,k) = 0.0 drat(i,j,k) = 0.0 fldec(i,j,k) = 0.0 stdec(i,j,k) = 0.0 flmass(i,j,k) = 0.0 stmass(i,j,k) = 0.0 rdght(i,j,k) = 0.0 rdgsp(i,j,k) = 0.0 rr(i,j,k) = 0.0 30 continue 20 continue 10 continue end if am0oif = .false. c + + + INPUT FORMATS + + + c + + + OUTPUT FORMATS + + + 2000 format (2(i3),1x,i5,1x,f10.4) 2100 format (2(i3),1x,f10.4) 999 format (6(i5),2f10.2) c + + + END SPECIFICATIONS + + + c get erosion arrays for output c if subregion has been grided by EROSION, set up grid for output if(am0gdf .eqv. .true.) then c determine x,y coordinates for each grid point do 100 i = 0, imax do 100 j = 0, jmax if (i .eq. 0) then gx(i) = -(0.5 * ix) else gx(i) = gx(i-1) + ix end if if (j .eq. 0) then gy(j) = -(0.5 * jy) else gy(j) = gy(j-1) + jy end if c determine max and min index value for accounting c region (simulation region for WEPS-Lite) c based on gx and gy if (gx(i) .lt. amxar(1,1,1)) xmin = i + 1 if (gx(i) .le. amxar(1,2,1)) xmax = i if (gy(j) .lt. amxar(2,1,1)) ymin = j + 1 if (gy(j) .le. amxar(2,2,1)) ymax = j 100 continue ngdpt = ((xmax - xmin) + 1) * ((ymax - ymin) + 1) if (ngdpt .eq. 0) then write (*,*) ' the accounting region is too small, ', & ' please expand by ',(2*ix),' meters in the x-direction and by ' & ,(2*jy), ' meters in the y-direction and restart WEPS ' go to 600 end if am0gdf = .false. end if c sum values and get max and min from EROSION, then reset to 0.0 c if gridded, get max, min, and sums tmp = 0.0 if (am0eif .eqv. .false.) then do 200 i = xmin, xmax do 200 j = ymin, ymax sumegtot = sumegtot + egt(i,j) sumegts = sumegts + egtss(i,j) sumegt10 = sumegt10 + egt10(i,j) if (egt(i,j) .lt. maxegt) maxegt = egt(i,j) if (egt(i,j) .gt. minegt) minegt = egt(i,j) if (awudmx .gt. 8.) then write(7,2000)cd,cm,cy,awadir write(7,2100)i,j,egt(i,j) end if egt(i,j) = 0.0 egtss(i,j) = 0.0 egt10(i,j) = 0.0 200 continue end if if (minegt .eq. -1000000.) minegt = 0.0 if (maxegt .eq. 1000000.) maxegt = 0.0 300 if (ngdpt .eq. 0) ngdpt = 1 avgetot = sumegtot/ngdpt avgets = sumegts/ngdpt avget10 = sumegt10/ngdpt write(40,*) cd,cm,cy,avgetot,minegt,maxegt,avgets,avget10, & abffcv(1),abfscv(1),aszrgh(1),asxrgs(1),aslrr(1) c determine output for operation days and last day of the month c for each rotation year if ((lopday .eq. cd) .and. (lopmon .eq. cm) .and. & (lopyr .eq. amnryr(1))) then opnam(lopday,lopmon,lopyr) = opname opimp(lopday,lopmon,lopyr) = impl opdir(lopday,lopmon,lopyr) = tdir end if if (cd .eq. lstday(cm,cy)) then lopday = cd lopmon = cm outday(lopday,lopmon,lopyr) = 1 end if if ((lopday .eq. cd) .and. (lopmon .eq. cm) .and. & (lopyr .eq. amnryr(1))) then fldec(lopday,lopmon,lopyr) = fldec(lopday,lopmon,lopyr) & + abffcv(1) stdec(lopday,lopmon,lopyr) = stdec(lopday,lopmon,lopyr) & + abfscv(1) flmass(lopday,lopmon,lopyr) = flmass(lopday,lopmon,lopyr) & + abmres(1) stmass(lopday,lopmon,lopyr) = stmass(lopday,lopmon,lopyr) & + abmst(1) rdght(lopday,lopmon,lopyr) = rdght(lopday,lopmon,lopyr) & + aszrgh(1) rdgsp(lopday,lopmon,lopyr)=rdgsp(lopday,lopmon,lopyr)+asxrgs(1) rr(lopday,lopmon,lopyr) = rr(lopday,lopmon,lopyr) + aslrr(1) outday(lopday,lopmon,lopyr) = 1 end if ndays(lopday,lopmon,lopyr) = ndays(lopday,lopmon,lopyr) + 1 tloss(lopday,lopmon,lopyr) = tloss(lopday,lopmon,lopyr) + avgetot std(lopday,lopmon,lopyr) = 0.0 sloss(lopday,lopmon,lopyr) = sloss(lopday,lopmon,lopyr) + avgets pm10loss(lopday,lopmon,lopyr) = pm10loss(lopday,lopmon,lopyr) + & avget10 cs1(lopday,lopmon,lopyr) = 0.2 * tloss(lopday,lopmon,lopyr) cs2(lopday,lopmon,lopyr) = 0.4 * tloss(lopday,lopmon,lopyr) cs3(lopday,lopmon,lopyr) = 0.3 * tloss(lopday,lopmon,lopyr) cs4(lopday,lopmon,lopyr) = 0.1 * tloss(lopday,lopmon,lopyr) ss1(lopday,lopmon,lopyr) = 0.2 * sloss(lopday,lopmon,lopyr) ss2(lopday,lopmon,lopyr) = 0.4 * sloss(lopday,lopmon,lopyr) ss3(lopday,lopmon,lopyr) = 0.3 * sloss(lopday,lopmon,lopyr) ss4(lopday,lopmon,lopyr) = 0.1 * sloss(lopday,lopmon,lopyr) pm1(lopday,lopmon,lopyr) = 0.1 * pm10loss(lopday,lopmon,lopyr) pm2(lopday,lopmon,lopyr) = 0.2 * pm10loss(lopday,lopmon,lopyr) pm3(lopday,lopmon,lopyr) = 0.4 * pm10loss(lopday,lopmon,lopyr) pm4(lopday,lopmon,lopyr) = 0.3 * pm10loss(lopday,lopmon,lopyr) precip(lopday,lopmon,lopyr) = precip(lopday,lopmon,lopyr)+awzdpt drat(lopday,lopmon,lopyr) = drat(lopday,lopmon,lopyr) + ah0drat if (awudav .ge. 8.0) then winde(lopday,lopmon,lopyr) = (winde(lopday,lopmon,lopyr) & + (0.5*awdair*(awudav**2)*(awudav - 8.0))) & / 1000. end if c determine output of average values for the rotation year outday(32,13,lopyr) = 1 ndays(32,13,lopyr) = ndays(32,13,lopyr) + 1 tloss(32,13,lopyr) = tloss(32,13,lopyr) + avgetot sloss(32,13,lopyr) = sloss(32,13,lopyr) + avgets pm10loss(32,13,lopyr) = pm10loss(32,13,lopyr) + avget10 cs1(32,13,lopyr) = 0.2*tloss(32,13,lopyr) cs2(32,13,lopyr) = 0.4*tloss(32,13,lopyr) cs3(32,13,lopyr) = 0.3*tloss(32,13,lopyr) cs4(32,13,lopyr) = 0.1*tloss(32,13,lopyr) ss1(32,13,lopyr) = 0.2*sloss(32,13,lopyr) ss2(32,13,lopyr) = 0.4*sloss(32,13,lopyr) ss3(32,13,lopyr) = 0.3*sloss(32,13,lopyr) ss4(32,13,lopyr) = 0.1*sloss(32,13,lopyr) pm1(32,13,lopyr) = 0.2*pm10loss(32,13,lopyr) pm2(32,13,lopyr) = 0.4*pm10loss(32,13,lopyr) pm3(32,13,lopyr) = 0.3*pm10loss(32,13,lopyr) pm4(32,13,lopyr) = 0.1*pm10loss(32,13,lopyr) precip(32,13,lopyr) = precip(32,13,lopyr) + awzdpt drat(32,13,lopyr) = drat(32,13,lopyr)+ ah0drat if (awudav .ge. 8.0) then winde(32,13,lopyr) = (winde(32,13,lopyr) & + (0.5*awdair*(awudav**2)*(awudav - 8.0))) & / 1000. end if c determine output for months regardless of rotation year ndays(32,cm,6) = ndays(32,cm,6) + 1 outday(32,lopmon,6) = 1 tloss(32,cm,6) = tloss(32,cm,6) + avgetot sloss(32,cm,6) = sloss(32,cm,6) + avgets pm10loss(32,cm,6) = pm10loss(32,cm,6) + avget10 cs1(32,cm,6) = 0.2 * tloss(32,cm,6) cs2(32,lopmon,6) = 0.4 * tloss(32,lopmon,6) cs3(32,lopmon,6) = 0.3 * tloss(32,lopmon,6) cs4(32,lopmon,6) = 0.1 * tloss(32,lopmon,6) ss1(32,lopmon,6) = 0.2 * sloss(32,lopmon,6) ss2(32,lopmon,6) = 0.4 * sloss(32,lopmon,6) ss3(32,lopmon,6) = 0.3 * sloss(32,lopmon,6) ss4(32,lopmon,6) = 0.1 * sloss(32,lopmon,6) pm1(32,lopmon,6) = 0.2 * pm10loss(32,lopmon,6) pm2(32,lopmon,6) = 0.4 * pm10loss(32,lopmon,6) pm3(32,lopmon,6) = 0.3 * pm10loss(32,lopmon,6) pm4(32,lopmon,6) = 0.1 * pm10loss(32,lopmon,6) precip(32,cm,6) = precip(32,cm,6) + awzdpt drat(32,cm,6) = drat(32,cm,6) + ah0drat if (awudav .ge. 8.0) then winde(32,cm,6) = (winde(32,cm,6) & + (0.5*awdair*(awudav**2)*(awudav - 8.0))) & / 1000. end if if (cd .eq. lstday(cm,cy)) then fldec(32,cm,6) = fldec(32,cm,6) + abffcv(1) stdec(32,cm,6) = stdec(32,cm,6) + abfscv(1) flmass(32,cm,6) = flmass(32,cm,6) + abmres(1) stmass(32,cm,6) = stmass(32,cm,6) + abmst(1) rdght(32,cm,6) = rdght(32,cm,6) + aszrgh(1) rdgsp(32,cm,6) = rdgsp(32,cm,6) + asxrgs(1) rr(32,cm,6) = rr(32,cm,6) + aslrr(1) outday(32,cm,6) = 1 end if c determine output of average values for simulation ndays(32,13,6) = ndays(32,13,6) + 1 outday(32,13,6) = 1 tloss(32,13,6) = tloss(32,13,6) + avgetot sloss(32,13,6) = sloss(32,13,6) + avgets pm10loss(32,13,6) = pm10loss(32,13,6) + avget10 cs1(32,13,6) = 0.2 * tloss(32,13,6) cs2(32,13,6) = 0.4 * tloss(32,13,6) cs3(32,13,6) = 0.3 * tloss(32,13,6) cs4(32,13,6) = 0.1 * tloss(32,13,6) ss1(32,13,6) = 0.2 * sloss(32,13,6) ss2(32,13,6) = 0.4 * sloss(32,13,6) ss3(32,13,6) = 0.3 * sloss(32,13,6) ss4(32,13,6) = 0.1 * sloss(32,13,6) pm1(32,13,6) = 0.2 * pm10loss(32,13,6) pm2(32,13,6) = 0.4 * pm10loss(32,13,6) pm3(32,13,6) = 0.3 * pm10loss(32,13,6) pm4(32,13,6) = 0.1 * pm10loss(32,13,6) precip(32,13,6) = precip(32,13,6) + awzdpt drat(32,13,6) = drat(32,13,6) + ah0drat if (awudav .ge. 8.0) then winde(32,13,6) = (winde(32,13,6) & + (0.5*awdair*(awudav**2)*(awudav - 8.0))) & / 1000. end if if ((cd .eq. 1).and.(cm .eq. 1)) outyrs(lopyr) = outyrs(lopyr) + 1 outcnt = outcnt + 1 maxegt = 1000000.0 minegt = -1000000.0 sumegtot = 0. sumegts = 0. sumegt10 = 0. 500 return 600 end