c file: 'bookeep.for' subroutine bookeep (cd, cm, cy) 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 + + + 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 'c1glob.inc' include 'd1glob.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' c + + + LOCAL VARIABLES + + + integer cd, cm, cy, i i, j, k, idx, i ndiy, ngdpt, numglosscs, numgdepcs, i rptdat, i xmax, xmin, ymax, ymin real avgetot, avgetcs, avgets, avget10, r avglosscs, avgdepcs, r apm1,apm2,apm3,apm4, r ass1,ass2,ass3,ass4, r acs1,acs2,acs3,acs4, r cslossarea, csdeparea, r cellarea, cumpre, cumprery, r gx(0:mngdpt), gy(0:mngdpt), r xl,yl, r maxegt, minegt, r sumegtot, sumegtcs, sumegts, sumegt10, r sumlosscs, sumdepcs, 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 amxar(xy,pt,ar) c - Coordinates of the accounting region where: c xy = 1 for x coordinate c xy = 2 for y coordinate c pt = 1 for first point, ie. lower left corner of region c pt = 2 for second point, ie. upper right corner of region c ar = number of the accounting region c eg. ar = 1 if only one acciunting region c Note that two numbers (ie x,y) for amxar are needed to c describe each point if the accounting region 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 avglosscs - The average creep + salatation loss only (deposition not included) c avgdepcs - The average creep + salatation deposition only (loss not included) 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 cellarea - The area of an individual grid cell. c cslossarea- The area of accounting region with creep+saltation loss. c csdeparea - The area of accounting region with deposition. 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 imax - Number of grid intervals in x direction on EROSION grid c ix - grid interval in x-direction(m) c jmax - Number of grid intervals in y direction on EROSION grid c jy - grid interval in y-direction (m) 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 xmax - max index value (numger of grid intervals) in x direction c for accounting region (= simulation region for WEPS1.0) c xmin - min index value (number of grid intervals) in x direction c for accounting region (= simulation region for WEPS1.0) c ymax - max index value (numger of grid intervals) in y direction c for accounting region (= simulation region for WEPS1.0) c ymin - min index value (number of grid intervals) in y direction c for accounting region (= simulation region for WEPS1.0) c xl - length of the accounting region in the x direction c yl - length of the accounting region in the y direction c + + + SUBROUTINES CALLED + + + c + + + FUNCTIONS CALLED + + + integer lstday integer dayear real biodrag c + + + DATA INITIALIZATIONS + + + data sumegtot, sumegtcs, sumegts, sumegt10 /4*0./ data sumlosscs, sumdepcs /2*0./ data maxegt, minegt /-1E9, 1e9/ data xmin /0.0/, xmax /0.0/, ymin /0.0/, ymax /0.0/ save xmin, xmax, ymin, ymax, ngdpt, cellarea 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,16 outyrs(k) = 0 opnam(i,j,k) = ' ' ndays(i,j,k) = 1 teros(i,j,k) = 0.0 sumx(i,j,k) = 0.0 sumxsq(i,j,k) = 0.0 std(i,j,k) = 0.0 cseros(i,j,k) = 0.0 seros(i,j,k) = 0.0 pm10eros(i,j,k) = 0.0 csloss(i,j,k) = 0.0 csdep(i,j,k) = 0.0 sumcslarea(i,j,k) = 0.0 sumcsdarea(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 ccovtot(i,j,k) = 0.0 cstsilh(i,j,k) = 0.0 cabmass(i,j,k) = 0.0 dflcov(i,j,k) = 0.0 dstsilh(i,j,k) = 0.0 dflmass(i,j,k) = 0.0 dstmass(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 rdgdir(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 numglosscs = 0 numgdepcs = 0 c + + + INPUT FORMATS + + + c + + + OUTPUT FORMATS + + + 2000 format (2(1x,i3),1x,i4,1x,f10.1,2(1x,i3),3(1x,f10.4)) 2100 format (2(i3),1x,2f10.4) 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 (i = x direction) and y (j = y direction) coordinates c for each grid point do 105 i = 0, imax do 100 j = 0, jmax if (i .eq. 0) then gx(i) = (-(0.5 * ix)) + amxar(1,1,1) else gx(i) = gx(i-1) + ix end if if (j .eq. 0) then gy(j) = (-(0.5 * jy)) + amxar(2,1,1) else gy(j) = gy(j-1) + jy end if c determine max and min index value for accounting c region (simulation region for WEPS1.0) c based on gx and gy c write(*,*)i,j,'gx',gx(i),' gy',gy(j) 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 105 continue c calculate the accounting region length in x direction xl = amxar(1,2,1) - amxar(1,1,1) c calculate the accounting region length in y direction yl = amxar(2,2,1) - amxar(2,1,1) c calculate area of a grid cell cellarea = (xl/imax)*(yl/jmax) c print accounting region information c write(*,*)'amxar(1,1,1)',amxar(1,1,1),' amxar(2,1,1)',amxar(2,1,1) c write(*,*)'amxar(1,2,1)',amxar(1,2,1),' amxar(2,2,1)',amxar(2,2,1) c write(*,*)'amxsim(1,1)',amxsim(1,1),' amxsim(1,2)',amxsim(1,2) c write(40,*) 'ix',ix,' jy',jy,' xl',xl,' yl',yl c write(40,*) 'xmax',xmax,' xmin',xmin,' ymax',ymax,' ymin',ymin c write(40,*) 'imax', imax, ' jmax', jmax c write(6,*) 'ix',ix,' jy',jy,' xl',xl,' yl',yl c write(6,*) 'xmax',xmax,' xmin',xmin,' ymax',ymax,' ymin',ymin c write(6,*) 'imax', imax, ' jmax', jmax ngdpt = ((xmax - xmin) + 1) * ((ymax - ymin) + 1) c write(*,*)'ngdpt',ngdpt 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. c else c ngdpt = 0 c cellarea = 0.0 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 if (am0eif .eqv. .false.) then do 205 i = xmin, xmax do 200 j = ymin, ymax egtcs(i,j) = egt(i,j) - egtss(i,j) sumegtot = sumegtot + egt(i,j) sumegtcs = sumegtcs + egtcs(i,j) sumegts = sumegts + egtss(i,j) sumegt10 = sumegt10 + egt10(i,j) c sum creep/salt loss and get number of grid points with loss if (egtcs(i,j) .lt. 0.0) then sumlosscs = sumlosscs + egtcs(i,j) numglosscs = numglosscs + 1 c write(*,*) 'sumlosscs',sumlosscs,' numg',numglosscs endif c sum creep/salt deposition and get number of grid points with deposition if (egtcs(i,j) .gt. 0.0) then sumdepcs = sumdepcs + egtcs(i,j) numgdepcs = numgdepcs + 1 c write(*,*) 'sumdepcs',sumdepcs,' numg',numgdepcs endif c determine the max and min soil loss from a grid point if (egt(i,j) .lt. maxegt) maxegt = egt(i,j) if (egt(i,j) .gt. minegt) minegt = egt(i,j) c if (egt(i,j) .ne. 0.0) then c write(7,2000)cd,cm,cy,awadir,i,j,egt(i,j), c & egtcs(i,j),egtss(i,j) c end if c if ((egt(i,j) .gt. 0.0).or.(egtcs(i,j) .gt. 0.0)) then c write(*,2000)'*****',cd,cm,cy,awadir,i,j,egt(i,j), c & egtcs(i,j),egtss(i,j) c end if 200 continue 205 continue c sum boundry loss - where extreme index values represent grid c points outside the accounting region c Note!!! c For the boundary case, egt is the creep+saltation not total soil loss do 201 i = 0, imax tcs1 = tcs1 + egt(i,0) tss1 = tss1 + egtss(i,0) tpm1 = tpm1 + egt10(i,0) c write(*,2100)i,0,egt(i,0), tcs1 tcs3 = tcs3 + egt(i,jmax) tss3 = tss3 + egtss(i,jmax) tpm3 = tpm3 + egt10(i,jmax) c write(*,*)'***',i,jmax,egt(i,jmax),tcs3,egtss(i,jmax) 201 continue do 202 j=0,jmax tcs2 = tcs2 + egt(0,j) tss2 = tss2 + egtss(0,j) tpm2 = tpm2 + egt10(0,j) c write(*,2100)0,j,egtcs(0,j), tcs2 tcs4 = tcs4 + egt(imax,j) tss4 = tss4 + egtss(imax,j) tpm4 = tpm4 + egt10(imax,j) c write(*,2100)imax,j,egtcs(imax,j),tcs4 202 continue c calculate average mass passing each boundary acs1 = tcs1 / (imax-1) acs3 = tcs3 / (imax-1) acs2 = tcs2 / (jmax-1) acs4 = tcs4 / (jmax-1) ass1 = tss1 / (imax-1) ass3 = tss3 / (imax-1) ass2 = tss2 / (jmax-1) ass4 = tss4 / (jmax-1) apm1 = tpm1 / (imax-1) apm3 = tpm3 / (imax-1) apm2 = tpm2 / (jmax-1) apm4 = tpm4 / (jmax-1) c reinitialize so that they will not be carried over to the next day tcs1 = 0. tss1 = 0. tpm1 = 0. tcs3 = 0. tss3 = 0. tpm3 = 0. tcs2 = 0. tss2 = 0. tpm2 = 0. tcs4 = 0. tss4 = 0. tpm4 = 0. do 215 i = 0, imax do 210 j = 0, jmax egt(i,j) = 0.0 egtcs(i,j) = 0.0 egtss(i,j) = 0.0 egt10(i,j) = 0.0 210 continue 215 continue end if if (minegt .eq. -1000000.) minegt = 0.0 if (maxegt .eq. 1000000.) maxegt = 0.0 if (ngdpt .eq. 0) ngdpt = 1 avgetot = sumegtot/ngdpt total=avgetot c write(*,*) cd, cm, cy, awadir, avgetot avgetcs = sumegtcs/ngdpt avgets = sumegts/ngdpt suspen = avgets avget10 = sumegt10/ngdpt pmten = avget10 c determine average loss and deposition and area of each if (numglosscs .eq. 0) then avglosscs = 0.0 cslossarea = 0.0 else avglosscs = sumlosscs / numglosscs cslossarea = cellarea * numglosscs / ngdpt end if if (numgdepcs .eq. 0) then avgdepcs = 0.0 csdeparea = 0.0 else avgdepcs = sumdepcs / numgdepcs csdeparea = cellarea * numgdepcs / ngdpt end if c write(*,*) 'sum', sumlosscs, ' num',numglosscs c & , ' avglosscs 1',avglosscs c if ((cd.eq.9).and.(cm.eq.5).and.(cy.eq.3.)) then c write(40,*) cd,cm,cy, 'ngdpt',ngdpt c write(40,*) 'sumegtot',sumegtot c write(40,*)'tot ',avgetot*xl*yl c write(40,*)'boundary ',(acs1+acs3+ass1+ass3) c & +(acs2+acs4+ass2+ass4) c write(40,*)'avgetot ',avgetot, c & (acs1+acs3+ass1+ass3)+(acs2+acs4+ass2+ass4) c write(40,*)'avgecs ',avgetcs,(acs1+acs2+acs3+acs4) c write(40,*)'avgess ',avgets,(ass1+ass2+ass3+ass4) c end if c write(40,*) cd,cm,cy,avgetot,minegt,maxegt,avgets,avget10, c & abffcv(1),abfscv(1),aszrgh(1),asxrgs(1),aslrr(1) c determine output for operation days and the 15 qnd last day of c the month for each rotation year c - first determine operation day if ((lopday .eq. cd) .and. (lopmon .eq. cm) .and. & (lopyr .eq. amnryr(1))) then opnam(lopday,lopmon,lopyr) = opname 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 leap year add the 29th to the 28th if((cm .eq. 2) .and. (cd .eq. 29)) then lopday = cd - 1 else lopday = cd end if lopmon = cm outday(lopday,lopmon,lopyr) = 1 end if if ((lopday .eq. cd) .and. (lopmon .eq. cm) .and. & (lopyr .eq. amnryr(1))) then ccovtot(lopday,lopmon,lopyr) = ccovtot(lopday,lopmon,lopyr) & + acftcv(1) cstsilh(lopday,lopmon,lopyr) = cstsilh(lopday,lopmon,lopyr) & + biodrag(acrlai(1), acrsai(1)) cabmass(lopday,lopmon,lopyr) = cabmass(lopday,lopmon,lopyr) & + acmst(1) c the first subscript is the dead biomass pool c this will need to be updated if number of pools change tmp = 0.0 do 340 idx = 1, mnbpls tmp = tmp + adftcv(idx,1) * (1.0-tmp) 340 continue dflcov(lopday,lopmon,lopyr) = dflcov(lopday,lopmon,lopyr)+tmp do 350 idx = 1, mnbpls dstsilh(lopday,lopmon,lopyr) = dstsilh(lopday,lopmon,lopyr) & + adrsai(idx,1) dflmass(lopday,lopmon,lopyr) = dflmass(lopday,lopmon,lopyr) & + admf(idx,1) dstmass(lopday,lopmon,lopyr) = dstmass(lopday,lopmon,lopyr) & + admst(idx,1) 350 continue 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) rdgdir(lopday,lopmon,lopyr) = rdgdir(lopday,lopmon,lopyr) & + asargo(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 c if leap year add the 29th to the 28th if ((cd .eq. 29) .and. (cm .eq. 2).and.(lopyr.eq.amnryr(1))) then ccovtot(28,2,lopyr) = ccovtot(28,2,mnryr) & + acftcv(1) cstsilh(28,2,lopyr) = cstsilh(28,2,mnryr) & + biodrag(acrlai(1), acrsai(1)) cabmass(28,2,lopyr) = cabmass(28,2,lopyr) & + acmst(1) c the first subscript is the dead biomass pool c this will need to be updated if number of pools change tmp = 0.0 do 354 idx = 1, mnbpls tmp = tmp + adftcv(idx,1) * (1.0-tmp) 354 continue dflcov(28,2,lopyr) = dflcov(28,2,lopyr)+tmp do 355 idx = 1, mnbpls dstsilh(28,2,lopyr) = dstsilh(28,2,lopyr) & + adrsai(idx,1) dflmass(28,2,lopyr) = dflmass(28,2,lopyr) & + admf(idx,1) dstmass(28,2,lopyr) = dstmass(28,2,lopyr) & + admst(idx,1) 355 continue fldec(28,2,lopyr) = fldec(28,2,lopyr) + abftcv(1) stdec(28,2,lopyr) = stdec(28,2,lopyr) & + biodrag(abrlai(1),abrsai(1)) flmass(28,2,lopyr) = flmass(28,2,lopyr) + abmf(1) stmass(28,2,lopyr) = stmass(28,2,lopyr) + abmst(1) rdgdir(28,2,lopyr) = rdgdir(28,2,lopyr) + asargo(1) rdght(28,2,lopyr) = rdght(28,2,lopyr) + aszrgh(1) rdgsp(28,2,lopyr)=rdgsp(28,2,lopyr)+asxrgs(1) rr(28,2,lopyr) = rr(28,2,lopyr) + aslrr(1) end if ndays(lopday,lopmon,lopyr) = ndays(lopday,lopmon,lopyr) + 1 teros(lopday,lopmon,lopyr) = teros(lopday,lopmon,lopyr) + avgetot sumxsq(lopday,lopmon,lopyr) = sumxsq(lopday,lopmon,lopyr) & + avgetot**2. sumx(lopday,lopmon,lopyr) = sumx(lopday,lopmon,lopyr) + avgetot cseros(lopday,lopmon,lopyr) = cseros(lopday,lopmon,lopyr)+avgetcs c sum average loss and deposition and area of each csloss(lopday,lopmon,lopyr) = csloss(lopday,lopmon,lopyr) & + avglosscs c write(*,*) 'csloss',csloss(lopday,lopmon,lopyr), c &' avglosscs',avglosscs csdep(lopday,lopmon,lopyr) = csdep(lopday,lopmon,lopyr)+avgdepcs sumcslarea(lopday,lopmon,lopyr) = sumcslarea(lopday,lopmon,lopyr) & + cslossarea sumcsdarea(lopday,lopmon,lopyr) = sumcsdarea(lopday,lopmon,lopyr) & + csdeparea seros(lopday,lopmon,lopyr) = seros(lopday,lopmon,lopyr) + avgets pm10eros(lopday,lopmon,lopyr) = pm10eros(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 (awudmx .gt. 8.0) then do 400 i = 1,ntstep if (awu(i) .gt. 8.0) then winde(lopday,lopmon,lopyr) = winde(lopday,lopmon,lopyr) & + (0.5 * awdair * (awu(i)**2.) * (awu(i) - 8.0)) & * (86400./ntstep) * 0.000001 end if 400 continue end if !write(6,*) awudmx, maxval(awu), (awu(i),i=1,ntstep) !write(6,*) '1 awudmx vs. maxval: ', awudmx, maxval(awu) c determine output of average values for the rotation year outday(32,13,lopyr) = 1 ndays(32,13,lopyr) = ndays(32,13,lopyr) + 1 teros(32,13,lopyr) = teros(32,13,lopyr) + avgetot sumxsq(32,13,lopyr) = sumxsq(32,13,lopyr) + avgetot**2. sumx(32,13,lopyr) = sumx(32,13,lopyr) + avgetot cseros(32,13,lopyr) = cseros(32,13,lopyr) + avgetcs c sum average loss and deposition and area of each csloss(32,13,lopyr) = csloss(32,13,lopyr) + avglosscs csdep(32,13,lopyr) = csdep(32,13,lopyr)+avgdepcs sumcslarea(32,13,lopyr) = sumcslarea(32,13,lopyr) + cslossarea sumcsdarea(32,13,lopyr) = sumcsdarea(32,13,lopyr) + csdeparea seros(32,13,lopyr) = seros(32,13,lopyr) + avgets pm10eros(32,13,lopyr) = pm10eros(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 (awudmx .gt. 8.0) then do 410 i = 1,ntstep if (awu(i) .gt. 8.0) then winde(32,13,lopyr) = winde(32,13,lopyr) & + (0.5 * awdair * (awu(i)**2.) * (awu(i) - 8.0)) & * (86400./ntstep) * 0.000001 end if 410 continue end if c write(6,*) 'awudmx vs. maxval: ', awudmx, maxval(awu) c determine output for months regardless of rotation year ndays(32,cm,16) = ndays(32,cm,16) + 1 outday(32,lopmon,16) = 1 teros(32,cm,16) = teros(32,cm,16) + avgetot sumxsq(32,cm,16) = sumxsq(32,cm,16) + avgetot**2. sumx(32,cm,16) = sumx(32,cm,16) + avgetot cseros(32,cm,16) = cseros(32,cm,16) + avgetcs c sum average loss and deposition and area of each csloss(32,cm,16) = csloss(32,cm,16) + avglosscs csdep(32,cm,16) = csdep(32,cm,16)+avgdepcs sumcslarea(32,cm,16) = sumcslarea(32,cm,16) + cslossarea sumcsdarea(32,cm,16) = sumcsdarea(32,cm,16) + csdeparea seros(32,cm,16) = seros(32,cm,16) + avgets pm10eros(32,cm,16) = pm10eros(32,cm,16) + avget10 cs1(32,cm,16) = cs1(32,cm,16) + acs1 cs2(32,cm,16) = cs2(32,cm,16) + acs2 cs3(32,cm,16) = cs3(32,cm,16) + acs3 cs4(32,cm,16) = cs4(32,cm,16) + acs4 ss1(32,cm,16) = ss1(32,cm,16) + ass1 ss2(32,cm,16) = ss2(32,cm,16) + ass2 ss3(32,cm,16) = ss3(32,cm,16) + ass3 ss4(32,cm,16) = ss4(32,cm,16) + ass4 pm1(32,cm,16) = pm1(32,cm,16) + apm1 pm2(32,cm,16) = pm2(32,cm,16) + apm2 pm3(32,cm,16) = pm3(32,cm,16) + apm3 pm4(32,cm,16) = pm4(32,cm,16) + apm4 precip(32,cm,16) = precip(32,cm,16) + awzdpt drat(32,cm,16) = drat(32,cm,16) + ah0drat if (awudmx .gt. 8.0) then do 420 i = 1,ntstep if (awu(i) .gt. 8.0) then winde(32,lopmon,16) = winde(32,lopmon,16) & + (0.5 * awdair * (awu(i)**2.) * (awu(i) - 8.0)) & * (86400./ntstep) * 0.000001 end if 420 continue end if if ((cd .eq. lstday(cm,cy)) .or. ((cd .eq.29).and.(cm.eq.2))) then ccovtot(32,cm,16) = ccovtot(28,2,mnryr) & + acftcv(1) cstsilh(32,cm,16) = cstsilh(28,2,mnryr) & + biodrag(acrlai(1), acrsai(1)) cabmass(32,cm,16) = cabmass(32,cm,16) & + acmst(1) c the first subscript is the dead biomass pool c this will need to be updated if number of pools change tmp = 0.0 do 424 idx = 1, mnbpls tmp = tmp + adftcv(idx,1) * (1.0-tmp) 424 continue dflcov(32,cm,16) = dflcov(32,cm,16)+tmp do 425 idx = 1, mnbpls dstsilh(32,cm,16) = dstsilh(32,cm,16) & + adrsai(idx,1) dflmass(32,cm,16) = dflmass(32,cm,16) & + admf(idx,1) dstmass(32,cm,16) = dstmass(32,cm,16) & + admst(idx,1) 425 continue fldec(32,cm,16) = fldec(32,cm,16) + abftcv(1) stdec(32,cm,16) = stdec(32,cm,16) + & biodrag(abrlai(1),abrsai(1)) flmass(32,cm,16) = flmass(32,cm,16) + abmf(1) stmass(32,cm,16) = stmass(32,cm,16) + abmst(1) rdgdir(32,cm,16) = rdgdir(32,cm,16) + asargo(1) rdght(32,cm,16) = rdght(32,cm,16) + aszrgh(1) rdgsp(32,cm,16) = rdgsp(32,cm,16) + asxrgs(1) rr(32,cm,16) = rr(32,cm,16) + aslrr(1) outday(32,cm,16) = 1 end if c determine output of average values for simulation ndays(32,13,16) = ndays(32,13,16) + 1 outday(32,13,16) = 1 teros(32,13,16) = teros(32,13,16) + avgetot sumxsq(32,13,16) = sumxsq(32,13,16) + avgetot**2. sumx(32,13,16) = sumx(32,13,16) + avgetot cseros(32,13,16) = cseros(32,13,16) + avgetcs c sum average loss and deposition and area of each csloss(32,13,16) = csloss(32,13,16) + avglosscs csdep(32,13,16) = csdep(32,13,16) + avgdepcs sumcslarea(32,13,16) = sumcslarea(32,13,16) + cslossarea sumcsdarea(32,13,16) = sumcsdarea(32,13,16) + csdeparea c write(*,*) 'loss', csloss(32,13,16),sumcslarea(32,13,16) c write(*,*) 'deposition', csdep(32,13,16) ,sumcsdarea(32,13,16) seros(32,13,16) = seros(32,13,16) + avgets pm10eros(32,13,16) = pm10eros(32,13,16) + avget10 cs1(32,13,16) = cs1(32,13,16) + acs1 cs2(32,13,16) = cs2(32,13,16) + acs2 cs3(32,13,16) = cs3(32,13,16) + acs3 cs4(32,13,16) = cs4(32,13,16) + acs4 ss1(32,13,16) = ss1(32,13,16) + ass1 ss2(32,13,16) = ss2(32,13,16) + ass2 ss3(32,13,16) = ss3(32,13,16) + ass3 ss4(32,13,16) = ss4(32,13,16) + ass4 pm1(32,13,16) = pm1(32,13,16) + apm1 pm2(32,13,16) = pm2(32,13,16) + apm2 pm3(32,13,16) = pm3(32,13,16) + apm3 pm4(32,13,16) = pm4(32,13,16) + apm4 precip(32,13,16) = precip(32,13,16) + awzdpt drat(32,13,16) = drat(32,13,16) + ah0drat if (awudmx .gt. 8.0) then do 430 i = 1,ntstep if (awu(i) .gt. 8.0) then winde(32,13,16) = winde(32,13,16) & + (0.5 * awdair * (awu(i)**2.) * (awu(i) - 8.0)) & * (86400./ntstep) * 0.000001 end if 430 continue 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. sumegtcs = 0. sumegts = 0. sumegt10 = 0. sumlosscs = 0. sumdepcs = 0. numglosscs = 0 numgdepcs = 0 C *** write(*,*) 'bookeep end' 500 return 600 end