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, rptdat, i xmax, xmin, ymax, ymin real avgetot, avgetcs, 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, sumegtcs, 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 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 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 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 + + + SUBROUTINES CALLED + + + c + + + FUNCTIONS CALLED + + + integer lstday integer dayear real biodrag c + + + DATA INITIALIZATIONS + + + data sumegtot, sumegtcs, sumegts, sumegt10 /4*0./ data maxegt, minegt /-1E9, 1e9/ data xmin /0.0/, xmax /0.0/, ymin /0.0/, ymax /0.0/ save xmin, xmax, ymin, ymax 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) = ' ' opdir(i,j,k) = 0.0 ndays(i,j,k) = 1 tloss(i,j,k) = 0.0 sumx(i,j,k) = 0.0 sumxsq(i,j,k) = 0.0 std(i,j,k) = 0.0 csloss(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 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 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(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),y (j = y direction) coordinates c for each grid point do 105 i = 0, imax !30 do 100 j = 0, jmax !2 if (i .eq. 0) then gx(i) = (-(0.5 * ix)) + amxar(1,1,1) ! ~ -25 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 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(*,*) 'ix',ix,' jy',jy c write(*,*) 'xmax',xmax,' xmin',xmin,' ymax',ymax,' ymin',ymin c write(*,*) '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. 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 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) 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 non erosion 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 = xmin, xmax do 210 j = ymin, ymax 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 300 if (ngdpt .eq. 0) ngdpt = 1 avgetot = sumegtot/ngdpt c write(*,*) cd, cm, cy, awadir, avgetot total = avgetot avgetcs = sumegtcs/ngdpt 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 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 opdir(lopday,lopmon,lopyr) = odir 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 dflcov(lopday,lopmon,lopyr) = dflcov(lopday,lopmon,lopyr) & + (adffcv(1,1) & + (1.0-adffcv(1,1)) * adffcv(2,1) & + (1.0-(1.0-adffcv(1,1)) * adffcv(2,1)) * adffcv(3,1)) 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) 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 dflcov(28,2,lopyr) = dflcov(28,2,lopyr) & + (adffcv(1,1) & + (1.0-adffcv(1,1)) * adffcv(2,1) & + (1.0-(1.0-adffcv(1,1)) * adffcv(2,1)) * adffcv(3,1)) 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) 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 tloss(lopday,lopmon,lopyr) = tloss(lopday,lopmon,lopyr) + avgetot sumxsq(lopday,lopmon,lopyr) = sumxsq(lopday,lopmon,lopyr) & + avgetot**2. sumx(lopday,lopmon,lopyr) = sumx(lopday,lopmon,lopyr) + avgetot csloss(lopday,lopmon,lopyr) = csloss(lopday,lopmon,lopyr)+avgetcs 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 (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 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 sumxsq(32,13,lopyr) = sumxsq(32,13,lopyr) + avgetot**2. sumx(32,13,lopyr) = sumx(32,13,lopyr) + avgetot csloss(32,13,lopyr) = csloss(32,13,lopyr) + avgetcs 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 (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 determine output for months regardless of rotation year ndays(32,cm,16) = ndays(32,cm,16) + 1 outday(32,lopmon,16) = 1 tloss(32,cm,16) = tloss(32,cm,16) + avgetot sumxsq(32,cm,16) = sumxsq(32,cm,16) + avgetot**2. sumx(32,cm,16) = sumx(32,cm,16) + avgetot csloss(32,cm,16) = csloss(32,cm,16) + avgetcs sloss(32,cm,16) = sloss(32,cm,16) + avgets pm10loss(32,cm,16) = pm10loss(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 dflcov(32,cm,16) = dflcov(32,cm,16) & + (adffcv(1,1) & + (1.0-adffcv(1,1)) * adffcv(2,1) & + (1.0-(1.0-adffcv(1,1)) * adffcv(2,1)) * adffcv(3,1)) 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) 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 tloss(32,13,16) = tloss(32,13,16) + avgetot sumxsq(32,13,16) = sumxsq(32,13,16) + avgetot**2. sumx(32,13,16) = sumx(32,13,16) + avgetot csloss(32,13,16) = csloss(32,13,16) + avgetcs sloss(32,13,16) = sloss(32,13,16) + avgets pm10loss(32,13,16) = pm10loss(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. C *** write(*,*) 'bookeep end' 500 return 600 end