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 'b1glob.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 'main/plot.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 apm1,apm2,apm3,apm4, r ass1,ass2,ass3,ass4, r acs1,acs2,acs3,acs4, r cumpre, cumprery, r gx(0:30), gy(0:30), r maxegt, minegt, r sumegtot, sumegts, sumegt10, r tmp, r tpm1,tpm2,tpm3,tpm4, r tss1,tss2,tss3,tss4, r tcs1,tcs2,tcs3,tcs4 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 tpm1-4 - Local variable which holds the total pm10 passing each c boundary each day. c tss1-4 - Local variable which holds the total suspension passing c each boundary each day. c tcs1-4 - Local variable which holds the total creep and saltation c passing each boundary each day. 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 real biodrag c + + + DATA INITIALIZATIONS + + + data sumegtot, sumegts, sumegt10 /3*0/ data maxegt, minegt /-1E9, 1e9/ if (am0oif) 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. tcs1 = 0.0 tcs2 = 0.0 tcs3 = 0.0 tcs4 = 0.0 tss1 = 0.0 tss2 = 0.0 tss3 = 0.0 tss4 = 0.0 tpm1 = 0.0 tpm2 = 0.0 tpm3 = 0.0 tpm4 = 0.0 c + + + INPUT FORMATS + + + c + + + OUTPUT FORMATS + + + 2000 format (2(i3),1x,i5,1x,f10.4) 2100 format (2(i3),1x,2f10.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 c write(7,2000)cd,cm,cy,awadir c 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 do 201 i=0,imax tcs1 = tcs1 + egt(i,0) tss1 = tss1 + egtss(i,0) tpm1 = tpm1 + egt10(i,0) c write(7,2100)i,0,egt(i,0), tcs1 tcs3 = tcs3 + egt(i,jmax) tss3 = tss3 + egtss(i,jmax) tpm3 = tpm3 + egt10(i,jmax) c write(7,2100)i,jmax,egt(i,jmax),tcs3 201 continue do 202 j=0,jmax tcs2 = tcs2 + egt(0,j) tss2 = tss2 + egtss(0,j) tpm2 = tpm2 + egt10(0,j) c write(7,2100)0,j,egt(0,j), tcs2 tcs4 = tcs4 + egt(imax,j) tss4 = tss4 + egtss(imax,j) tpm4 = tpm4 + egt10(imax,j) c write(7,2100)imax,j,egt(imax,j),tcs4 202 continue c calculate average mass passing each boundary c C *** write(*,*) 'bookeep 1' acs1 = tcs1 / ((imax-1)*ix) acs3 = tcs3 / ((imax-1)*ix) acs2 = tcs2 / ((jmax-1)*jy) acs4 = tcs4 / ((jmax-1)*jy) ass1 = tss1 / ((imax-1)*ix) ass3 = tss3 / ((imax-1)*ix) ass2 = tss2 / ((jmax-1)*jy) ass4 = tss4 / ((jmax-1)*jy) apm1 = tpm1 / ((imax-1)*ix) apm3 = tpm3 / ((imax-1)*ix) apm2 = tpm2 / ((jmax-1)*jy) apm4 = tpm4 / ((jmax-1)*jy) 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 total = avgetot avgets = sumegts/ngdpt suspen = avgets avget10 = sumegt10/ngdpt pmten = avget10 write(40,*) cd,cm,cy,avgetot,minegt,maxegt,avgets,avget10, & abffcv(1),abfscv(1),aszrgh(1),asxrgs(1),aslrr(1) C *** write(*,*) 'bookeep 2' 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 c opimp(lopday,lopmon,lopyr) = impl opdir(lopday,lopmon,lopyr) = tdir end if c make operation name available for plot output if ((lopday .eq. cd) .and. (lopmon .eq. cm) .and. & (lopyr .eq. amnryr(1))) then operat = opname else operat = ' ' end if if ((cd .eq. lstday(cm,cy)) .or. (cd .eq.15)) then c if((cm .eq. 2) .and. (cd .eq. 29)) then c lopday = cd - 1 c else lopday = cd c end if lopmon = cm outday(lopday,lopmon,lopyr) = 1 end if C *** write(*,*) 'bookeep 3' if ((lopday .eq. cd) .and. (lopmon .eq. cm) .and. & (lopyr .eq. amnryr(1))) then fldec(lopday,lopmon,lopyr) = fldec(lopday,lopmon,lopyr) & + abftcv(1) stdec(lopday,lopmon,lopyr) = stdec(lopday,lopmon,lopyr) & + biodrag(abrlai(1),abrsai(1)) flmass(lopday,lopmon,lopyr) = flmass(lopday,lopmon,lopyr) & + abmf(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) = cs1(lopday,lopmon,lopyr) + acs1 cs2(lopday,lopmon,lopyr) = cs2(lopday,lopmon,lopyr) + acs2 cs3(lopday,lopmon,lopyr) = cs3(lopday,lopmon,lopyr) + acs3 cs4(lopday,lopmon,lopyr) = cs4(lopday,lopmon,lopyr) + acs4 ss1(lopday,lopmon,lopyr) = ss1(lopday,lopmon,lopyr) + ass1 ss2(lopday,lopmon,lopyr) = ss2(lopday,lopmon,lopyr) + ass2 ss3(lopday,lopmon,lopyr) = ss3(lopday,lopmon,lopyr) + ass3 ss4(lopday,lopmon,lopyr) = ss4(lopday,lopmon,lopyr) + ass4 pm1(lopday,lopmon,lopyr) = pm1(lopday,lopmon,lopyr) + apm1 pm2(lopday,lopmon,lopyr) = pm2(lopday,lopmon,lopyr) + apm2 pm3(lopday,lopmon,lopyr) = pm3(lopday,lopmon,lopyr) + apm3 pm4(lopday,lopmon,lopyr) = pm4(lopday,lopmon,lopyr) + apm4 precip(lopday,lopmon,lopyr) = precip(lopday,lopmon,lopyr)+awzdpt C *** write(*,*) 'bookeep 4' 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) = cs1(32,13,lopyr) + acs1 cs2(32,13,lopyr) = cs2(32,13,lopyr) + acs2 cs3(32,13,lopyr) = cs3(32,13,lopyr) + acs3 cs4(32,13,lopyr) = cs4(32,13,lopyr) + acs4 ss1(32,13,lopyr) = ss1(32,13,lopyr) + ass1 ss2(32,13,lopyr) = ss2(32,13,lopyr) + ass2 ss3(32,13,lopyr) = ss3(32,13,lopyr) + ass3 ss4(32,13,lopyr) = ss4(32,13,lopyr) + ass4 pm1(32,13,lopyr) = pm1(32,13,lopyr) + apm1 pm2(32,13,lopyr) = pm2(32,13,lopyr) + apm2 pm3(32,13,lopyr) = pm3(32,13,lopyr) + apm3 pm4(32,13,lopyr) = pm4(32,13,lopyr) + apm4 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) = cs1(32,cm,6) + acs1 cs2(32,cm,6) = cs2(32,cm,6) + acs2 cs3(32,cm,6) = cs3(32,cm,6) + acs3 cs4(32,cm,6) = cs4(32,cm,6) + acs4 ss1(32,cm,6) = ss1(32,cm,6) + ass1 ss2(32,cm,6) = ss2(32,cm,6) + ass2 ss3(32,cm,6) = ss3(32,cm,6) + ass3 ss4(32,cm,6) = ss4(32,cm,6) + ass4 pm1(32,cm,6) = pm1(32,cm,6) + apm1 pm2(32,cm,6) = pm2(32,cm,6) + apm2 pm3(32,cm,6) = pm3(32,cm,6) + apm3 pm4(32,cm,6) = pm4(32,cm,6) + apm4 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) + abftcv(1) stdec(32,cm,6) = stdec(32,cm,6) + & biodrag(abrlai(1),abrsai(1)) flmass(32,cm,6) = flmass(32,cm,6) + abmf(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) = cs1(32,13,6) + acs1 cs2(32,13,6) = cs2(32,13,6) + acs2 cs3(32,13,6) = cs3(32,13,6) + acs3 cs4(32,13,6) = cs4(32,13,6) + acs4 ss1(32,13,6) = ss1(32,13,6) + ass1 ss2(32,13,6) = ss2(32,13,6) + ass2 ss3(32,13,6) = ss3(32,13,6) + ass3 ss4(32,13,6) = ss4(32,13,6) + ass4 pm1(32,13,6) = pm1(32,13,6) + apm1 pm2(32,13,6) = pm2(32,13,6) + apm2 pm3(32,13,6) = pm3(32,13,6) + apm3 pm4(32,13,6) = pm4(32,13,6) + apm4 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. C *** write(*,*) 'bookeep end' 500 return 600 end