c file: 'weps.for' c program WEPS c + + + PURPOSE + + + c This is the MAIN program of the Wind Erosion Prediction System. c Its purpose is to control the sequence of events during a c simulation run. c author: John Tatarko c version: 92.10 c + + + KEY WORDS + + + c wind, erosion, hydrology, tillage, soil, c crop, decomposition, management c#include c external continue_hdl c external common_handler c + + + GLOBAL COMMON BLOCKS + + + *$noereference include 'p1werm.inc' include 'p1unconv.inc' c$INCLUDE:'p1strlen.inc' c$INCLUDE:'p1const.inc' include 'm1subr.inc' include 'm1sim.inc' include 'm1geo.inc' include 'm1flag.inc' include 'm1dbug.inc' c$INCLUDE:'m1eros.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' c$INCLUDE:'s1psd.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 'b1glob.inc' include 'b1mass.inc' include 'b1geom.inc' include 'w1clig.inc' include 'w1wind.inc' c$INCLUDE:'w1info.inc' include 'w1pavg.inc' c$INCLUDE:'h1et.inc' include 'h1hydro.inc' include 'h1scs.inc' include 'h1db1.inc' include 'h1temp.inc' c$INCLUDE:'asd.inc' c + + + LOCAL COMMON BLOCKS + + + include 'main/main.inc' include 'main/output.inc' include 'manage/man.inc' include 'erosion/m2geo.inc' include 'erosion/e2erod.inc' *$INCLUDE fsublib.fi *$reference c + + + LOCAL VARIABLES + + + character header*80 logical cflag, wflag, wcflag, first c*$INCLUDE integer.wat integer cd, cm, cy, ccd, ccm, ccy, cwd, cwm, cwy, c i diff, div, i ndiy, i i, isr, ioc, l, j, k, i simyrs, i xmax, xmin, ymax, ymin, yrsim real avgetot, avgets, avget10, brcd, cropcvr, dur, r grad, r maxegt, minegt, r sumegtot, sumegts, sumegt10, r tp, r wdir, wvel, r xmav c used for debugging under unix c integer ieeer c integer ieee_handler c character out*80 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 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 ccd - The current day of the CLIGEN file. c ccm - The current month of the CLIGEN file. c ccy - The current year of the CLIGEN file. c cflag - This logical variable controls output of warning c messages for mismatch of CLIGEN and simulation dates. 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 clifil - This variable holds the CLIGEN input file name. c cwd - The current day of the WINDGEN file. c cwm - The current month of the WINDGEN file. c cwy - The current year of the WINDGEN file. c daysim - This variable holds the total current days of simulation. c diff - This variable holds the number of simulation days. c div c grad - Global radiation (ly/day) as read in from CLIGEN. 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 iy - Starting year of simulation run (used in management). c l - This variable is an index on soil layers. c ld,lm,ly - The last day, month, and year of simulation, including c surface initialization. c ljday - This variable contains the last julian day of c the simulation run. c maxper - The maximum number of years in a rotation of 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 mperod - The number of years in a management rotation. This c variable is defined in (/inc/manage/man.inc) c simyrs - The number of years in a simulation run excluding the c years for surface initialization. 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 winfil - This variable holds the WINDGEN input file 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 calcwu - Subdaily wind speed generation c caldat - Converts julian day to day, month, and year (cd,cm,cy) c cdbug - Prints global variables before and after call to CROP c crop - Crop submodel c ddbug - Prints global variables before and after call to DECOMP c decomp - Decomposition submodel c erode - Erosion submodel c genrep - Output the general report c growhandles - Expands the number of handles allowed. c hdbug - Prints global variables before and after call to HYDRO c hydro - Hydrology submodel main program c input - Open files and perform input c manage - Management (tillage) submodel main program c mfinit - Management initiation subroutines c sdbug - Prints global variables before and after call to SOIL c soil - Soil submodel c + + + FUNCTIONS CALLED + + + integer dayear, julday logical isleap c + + + UNIT NUMBERS FOR INPUT/OUTPUT DEVICES + + + c * = screen and keyboard c 1 = simulation run file c 2 = general report output file c 3 = CLIGEN run file c 4 = WINDGEN run file c 5 = Reserved c 6 = Reserved - screen c 7 = Reserved c 8 = sub-daily wind speed data file c 9 = SOIL/HYDROLOGY run file c 10 = management (tillage) run file c 11 = decomp run file (not now used) c 12 = 'water.out' - hourly hydrology output file c 13 = 'temp.out' - daily soil temperature output file c 14 = 'hydro.out' - daily hydrology output file c 15 = 'manage.out' - management output file c 16 = 'soil.out' - soil output file c 17 = 'crop.out' - crop output file c 18 = 'dabove.out' - decomp above ground output file c 19 = 'dbelow.out' - decomp below ground output file c 20 = 'eros.out' - erosion output file c 21 = 'eegt' - period loss or deposition from grid point file c 22 = 'eegtss - period suspension loss from grid point file c 23 = 'eegt10' - period PM-10 loss from grid point file c 24 = hourly wind distribution output file (subday.out) c 25 = debug hydro c 26 = debug soil c 27 = debug crop c 28 = debug decomp c 29 = debug management c 31 = '../db/crop.db' - crop growth and decomposition parameter database c 32 = 'plot.out' - plotting output file c 40 = temporary file to hold accounting region erosion values c + + + DATA INITIALIZATIONS + + + avgetot = 0. avgets = 0. avget10 = 0. daysim = 0 outcnt = 0 sumegtot = 0.0 sumegts = 0.0 sumegt10 = 0.0 maxegt = 1000000.0 minegt = -1000000.0 xmin = 0.0 xmax = 0.0 ymin = 0.0 ymax = 0.0 yrsim = 0 c set initialization flags am0dif = .true. am0eif = .true. am0sif = .true. c no crop growing at start of simulation am0cgf = .false. am0ifl = .true. cflag = .false. subflg = .true. wflag = .false. wcflag = .false. c set grid flag until first gridding is done am0gdf = .false. c set output flag to initialize output arrays am0oif = .true. first = .true. c + + + INPUT FORMATS + + + 1020 format (a80) 1025 format (32x, f7.2) 1026 format (1x,12f6.2) 1030 format (2(2x,i2),1x,i4,1x,2f6.2,f5.2,1x,f6.2,3f7.2,f6.2,2f7.2) 1040 format (2(1x, i2), 1x, i4, 4f6.1) c + + + OUTPUT FORMATS + + + 2030 format (' warning !',28x,' day month year') 2035 format (' warning !',24x,' day month year') 2040 format (' current simulation date - ',i2,9x,i2,8x,i4, & /,' does not match current WINDGEN date - ',i2,9x,i2,8x,i4, & /) 2050 format (' current simulation date - ',i2,9x,i2,8x,i4, & /,' does not match current CLIGEN date - ',i2,9x,i2,8x,i4, & /) 2055 format (' current CLIGEN date - ',i2,9x,i2,8x,i4, & /,' is the end of file - rewinding to top of CLIGEN file', &/) 2060 format (' current CLIGEN date - ',i2,9x,i2,8x,i4, & /,' does not match current WINDGEN date - ',i2,9x,i2,8x,i4, & /) 2065 format (' current WINDGEN date - ',i2,9x,i2,8x,i4, & /,' is the end of file - rewinding to top of WINDGEN file', & /) 2070 format (1x,'day ',i6,' of ',i6) 2080 format (1x, 5f10.3) 2100 format (2(i3),i5,2(i3),1x,f10.4) c 2500 format (2(i2,'/'),i2'|',a3,1x,a20,a20,'|',4(f9.4,'|'), c & 12(f8.4,'|'),2(f6.2,'|'),f6.1,'|',7(f6.2,'|')) 2500 format (2(i2,'/'),i2'|',f3.0,1x,a20,a20,'|',4(f9.4,'|'), & 12(f8.4,'|'),2(f6.2,'|'),f6.1,'|',7(f6.2,'|')) 2520 format ('#',5x,(i1,'|'),41x,'N/A|',4(f9.4,'|'),12(f8.4,'|'), & 2(f6.2,'|'),f6.1,'|',7('N/A|')) 2540 format (2x,i2,'|',48x,'|',4(f9.4,'|'), & 12(f8.4,'|'),10(f6.2,'|')) 2560 format('*|',47x,'N/A|',4(f9.4,'|'),12(f8.4,'|'), &2(f6.2,'|'),f6.1,'|',7('N/A|')) c + + + END SPECIFICATIONS + + + c*$include growhnd.wat c call ieee_flags('clear', 'exception', 'all', out) c ieeer = ieee_handler('set','division',common_handler) c ieeer = ieee_handler('set','division',SIGFPE_ABORT) c if (ieeer .ne. 0) print*, 'could not set common_handler' c ieeer = ieee_handler('set','inexact',SIGFPE_IGNORE) c ieeer = ieee_handler('set','underflow',SIGFPE_IGNORE) c ieeer = ieee_handler('set','all',continue_hdl) *$ call growhandles (60) c open input files and read run files call input c temporarily initialize old random roughness aslrrc(1) = 10. as0rrk(1) = 0.9 c open crop and decomp parameter file open (unit = 31, file = '../db/crop.db',status='old') c open erosion output files open (unit = 20, file = '../out/erosion.out') open (unit = 21, file = '../out/eegt.out') c open (unit = 22, file = '../out/eegtss.out') c open (unit = 23, file = '../out/eegt10.out') open (unit = 7, file = 'maxmin.out') c open plot data file open(unit = 32, file = '../plot/plot.out', status = 'unknown') c open temporary file to hold accounting region erosion values open(unit = 40, file = 'eros.tmp', status = 'unknown') open(unit = 41, file = 'output.tmp', status = 'unknown') c Initialize the management file and rotation counters c Management will eventually open files with unit #'s "51" c through "50+nsubr" (not yet implemented - may be temporary c solution anyway) c munit = 50 do 4 isr = 1, nsubr c call mfinit(isr, munit+isr, tinfil) c Have management open same file each time with old unit # (10) c for now (can only handle one subregion yet anyway) call mfinit(isr, 10, tinfil) 4 continue c Likely that we will put all management data into memory c and only read and initialize everything here, looping through c each management file (one for each subregion). maxper = 1 do 5 isr = 1, nsubr call decini(isr) amnryr(isr) = 1 if (mperod(isr) .gt. maxper) maxper = mperod(isr) 5 continue c convert soil layer thickness (mm) to depth to bottom of layer (mm) c and initialize applied NO3 do 10 isr=1,nsubr aszlyd(1, isr) = aszlyt(1, isr) do 15 l = 2, nslay(isr) aszlyd(l,isr) = aszlyt(l,isr) + aszlyd(l-1, isr) 15 continue c temporarily initialize residue to bypass an error in crop-jt 7/22/94 c remove when problem is solved in crop do 16 l = 1, nslay(isr) admb(1,l,isr) = 0.001 16 continue asmno3(isr) = 0. 10 continue c read cligen header - remove when user interface is installed do 20 i = 1, 3 read(3,1020) header 20 continue read(3,1025) awtyav read(3,1026) (awtmav(i), i = 1,12) do 25 i = 1, 3 read(3,1020) header 25 continue c read windgen header do 30 i = 1, 7 read (4,1020) header 30 continue c calculate number of simulation days ijday = julday(id, im, iy) c add first rotation (for surface initialization) to last year ly = ly + maxper ljday = julday(ld, lm, ly) c begin simulation days do 50 am0jd = ijday,ljday tillday = .false. call caldat (am0jd, cd, cm, cy) c determine number of days in the year if (isleap(cy) .eqv. .true.) then ndiy = 366 else ndiy = 365 end if c read weather cligen and windgen daily data c at end of file rewind to top read(unit=3, fmt=1030, iostat=ioc) ccd,ccm,ccy,awzdpt,dur,tp, & xmav,awtdmx,awtdmn,grad,wvel,wdir,awtdpt if (ioc .eq. -1) then rewind 3 write(6,2035) write(6,2055) ccd, ccm, ccy do 31 i = 1, 8 read(3,1020) header 31 continue read(unit=3, fmt=1030) ccd,ccm,ccy,awzdpt,dur,tp, & xmav,awtdmx,awtdmn,grad,wvel,wdir,awtdpt ioc = 0 end if read(unit=4, fmt=1040, iostat=ioc) cwd,cwm,cwy,awadir,awudmx, & awudmn,awhrmx if (ioc .eq. -1) then rewind 4 write(6,2035) write(6,2065) cwd, cwm, cwy do 41 i = 1, 7 read (4,1020) header 41 continue read(4,1040) cwd, cwm, cwy, awadir, awudmx, awudmn, awhrmx ioc = 0 end if aweirr = grad * 0.04186 awudav = (awudmx + awudmn) / 2. c calculate air density from temperature and pressure awtdav = (awtdmx + awtdmn) / 2. awdair = 348.56 * (1.013-0.1183*(amzele/1000.) & + 0.0048 * (amzele/1000.)**2.) / (awtdav + 273.1) c test for date mismatch with weather files c print out warning on first mismatch only, skip thereafter if ((wflag .eqv. .true.) .and. (cflag .eqv. .true.) .and. & (wcflag .eqv. .true.)) go to 57 if (wflag .eqv. .true.) go to 55 if ((cd .ne. cwd) .or. (cm .ne. cwm) .or. (cy .ne. cwy)) then write (6,2030) write (6,2040) cd,cm,cy,cwd,cwm,cwy wflag = .true. end if 55 if (cflag .eqv. .true.) go to 56 if ((cd .ne. ccd) .or. (cm .ne. ccm) .or. (cy .ne. ccy)) then write (6,2030) write (6,2050) cd, cm, cy, ccd, ccm, ccy cflag = .true. end if 56 if (wcflag .eqv. .true.) go to 57 if ((ccd .ne. cwd) .or. (ccm .ne. cwm) .or. (ccy .ne. cwy)) & then write (6,2030) write (6,2060) ccd, ccm, ccy, cwd, cwm, cwy wcflag = .true. end if c print current day of simulation to screen periodically 57 daysim = daysim + 1 c diff = ljday - ijday + 1 c div = 100 c if (mod(daysim, div) .eq. 0) then c write(*,2070) daysim, diff c end if if ((cm .eq. 1) .and. (cd .eq. 1)) then yrsim = yrsim + 1 simyrs = (ly - iy + 1) - maxper if(am0sif .eqv. .true.) write(*,*) 'Initializing Surface' write(*,*) 'Year', yrsim, ' of', (ly - iy + 1) end if c do subregions 36 do 60 isr = 1, nsubr am0csr = isr c call HYDROLOGY submodel c ieeer = ieee_handler('set','underflow',SIGFPE_DEFAULT) if (am0hdb .eq. 1) then call hdbug(cd, cm, cy, isr, nslay(isr)) end if call hydro( nslay(isr), amrslp(isr), & acftcv(isr), acrlai(isr), & admres(isr), aczrtd(isr), ahfwsf(isr), & aszlyd(1, isr), asdblk(1, isr), & ahrwc(1, isr), ahrwcs(1, isr), & ahrwcf(1, isr), ahrwcw(1, isr), ahrwca(1,isr), & ah0cb(1,isr), aheaep(1,isr), & asfsan(1,isr), asfsil(1,isr), asfcla(1,isr), & ah0cng(isr), ah0cnp(isr), & ahzper(isr), ahzirr(isr), ahzrun(isr), & awudav, ahrsk(1, isr), & ahtsmx(1, isr), ahtsmn(1, isr), & ahrwc0(1, isr), daysim, & asfald(isr), asfalw(isr) ) if (am0hdb .eq. 1) then call hdbug(cd, cm, cy, isr, nslay(isr)) end if c if (am0tdb .eq. 1) then c write (29,*) '** Date ',cd, cm, cy c end if c call MANAGEment (tillage) submodel c set the flag am0hrvfl to zero indicating no harvest has c been done. When manage does a harvest this flag will be set to c 1 or 2 depending on the type of harvest, and will remain 1 or 2 c until the next simulation day (next call to manage) This flag c replaces am0nkfl. ANH 12/7/95 am0hrvfl = 0 C call manage (munit(sr)+isr, isr, cd, cm, cy,iy,lopday,lopmon,lopyr) call manage (10, isr, cd, cm, cy,iy,lopday,lopmon,lopyr) c call SOIL submodel if (am0sdb .eq. 1) then call sdbug(cd, cm, cy, isr, nslay(isr)) end if c calculate crop cover cropcvr = 1. - exp(-0.65*acrlai(isr)) call soil (daysim, am0irr(isr), ahzirr(isr), ahzsmt(isr), & ahtsmx(1,isr), ahtsmn(1,isr), & ahrwc(1,isr), ahrwca(1,isr), ahrwcw(1,isr), & asfom(1,isr), aszlyt(1,isr), nslay(isr), & asfsan(1,isr), asfsil(1,isr), asfcla(1,isr), & asxrgs(isr), aszrgh(isr), aslrr(isr), & aszcr(isr), asfcr(isr), asecr(isr), & asdcr(isr), asmlos(isr), asflos(isr), & asdblk(0,isr), asdagd(0,isr), & aslagm(0,isr), aslagn(0,isr), & as0ags(0,isr), aslagx(0,isr), aseags(0,isr), & cropcvr, ac0rg(isr)) if (am0sdb .eq. 1) then call sdbug(cd, cm, cy, isr, nslay(isr)) end if aszlyd(1, isr) = aszlyt(1, isr) do 70 l = 2, nslay(isr) aszlyd(l,isr) = aszlyt(l,isr) + aszlyd(l-1, isr) 70 continue c call CROP submodel if between planting and killing harvest if (am0cgf .eqv. .true.) then if (am0cdb .eq. 1) then call cdbug(cd, cm, cy, isr, nslay(isr)) end if call crop (nslay(isr), aszlyd(1,isr), asdblk(1,isr), & asfcce(1,isr), asfom(1,isr), asfcec(1,isr), asfsmb(1,isr), & asfcla(1,isr), as0ph(1,isr), asftan(1,isr), asftap(1,isr), & acrlai(isr), aczrtd(isr), acmst(isr), asmno3(isr), & ac0bn1(isr), ac0bn2(isr), ac0bn3(isr), & ac0bp1(isr), ac0bp2(isr), ac0bp3(isr), ac0ck(isr), & acrhi(isr), acehu0(isr), aczmxc(isr), ac0idc(isr), & acephu(isr), aczmrt(isr), actmin(isr), actopt(isr), & ac0fd1(1,isr), ac0fd2(1,isr), ac0fd1(2,isr), ac0fd2(2,isr), & ac0be1(1,isr), ac0be2(1,isr), ac0be1(2,isr),ac0be2(2,isr), & admb(1,1,isr), ac0alf(isr), ac0blf(isr), ac0clf(isr), & ac0dlf(isr), ac0arp(isr), ac0brp(isr), ac0crp(isr), & ac0drp(isr), ac0aht(isr), ac0bht(isr), ac0ssa(isr), & ac0ssb(isr), ac0sla(isr), ac0hue(isr), ac0lfe(isr)) if (am0cdb .eq. 1) then call cdbug(cd, cm, cy, isr, nslay(isr)) end if end if c call the DECOMPosition submodel if (am0ddb .eq. 1) then call ddbug(cd, cm, cy, isr, nslay(isr)) end if call decomp(isr) if (am0ddb .eq. 1) then call ddbug(cd, cm, cy, isr, nslay(isr)) end if if ((am0dfl .eq. 1).or.(am0dfl .eq. 2).or.(am0dfl .eq.3)) & then call decout (cd, cm, cy) end if c sum live and dead biomass abffcv(isr) = adffcv(isr) brcd = acrlai(isr)*(0.2-(0.15)* exp(-8.*acrlai(isr))) & +acrsai(isr) abfscv(isr) = adfscv(isr) + brcd c the next two statements were added since these should not be c allowed to go over 1. fix this jt 6/25/98 if (abfscv(isr) .gt. 1.) abfscv(isr) = 1. if (abffcv(isr) .gt. 1.) abffcv(isr) = 1. flmres(isr)=0. stmres(isr)=0. do 59 i=1,3 flmres(isr) = admf(i,isr) + flmres(isr) stmres(isr) = admst(i,isr) + stmres(isr) 59 continue abmst(isr) = acmst(isr) + stmres(isr) abmres(isr) = flmres(isr) 60 continue c skip next section until surface is initialized if (am0sif .eqv. .true.) goto 90 c if wind is great enough, call erosion if (awudmx .gt. 8.) then call calcwu call erosion end if c set initialization flag to .false. after first day am0ifl = .false. c print to plot data file write (32, 2080) awudmx, awzdpt, ahrwc0(1, 1), & aszrgh(1), acrlai(1) c do bookkeeping for output of general report call bookeep (cd, cm, cy, xmax, xmin, ymax, ymin) c check for end of first rotation (surface initialization) 90 do 95 isr = 1, nsubr c if last day of year, check for end of rotation if (dayear(cd,cm,cy) .eq. ndiy) then if (amnryr(isr) .eq. maxper) then c first time surface initilization is done set c write(*,*) amnryr(isr),maxper am0sif = .false. end if c if current rotation year is the last, reset to one if (mperod(isr) .eq. amnryr(isr)) then amnryr(isr) = 1 lopyr = amnryr(isr) else c if surface init done - increment rotation year counter amnryr(isr) = amnryr(isr) + 1 lopyr = amnryr(isr) end if end if 95 continue 50 continue close (unit = 3) close (unit = 4) close (unit = 40) c output the simulation general report c call genrep c call wsum c adjust output to print previous period values on current date c this next section is not the best way to do this - jt c****************************************************************** do 451 k = 1,5 do 461 j = 1,12 do 471 i = 1,31 if (outday(i,j,k) .eq. 1) then if (first .eqv. .true.) then ttmp1 = tloss(i,j,k) stmp1 = sloss(i,j,k) pm10tmp1 = pm10loss(i,j,k) cs1tmp1 = cs1(i,j,k) cs2tmp1 = cs2(i,j,k) cs3tmp1 = cs3(i,j,k) cs4tmp1 = cs4(i,j,k) ss1tmp1 = ss1(i,j,k) ss2tmp1 = ss2(i,j,k) ss3tmp1 = ss3(i,j,k) ss4tmp1 = ss4(i,j,k) pm1tmp1 = pm1(i,j,k) pm2tmp1 = pm2(i,j,k) pm3tmp1 = pm3(i,j,k) pm4tmp1 = pm4(i,j,k) ptmp1 = precip(i,j,k) dtmp1 = drat(i,j,k) wtmp1 = winde(i,j,k) tloss(i,j,k) = tloss(31,12,1) sloss(i,j,k) = sloss(31,12,1) pm10loss(i,j,k) = pm10loss(31,12,1) cs1(i,j,k) = cs1(31,12,1) cs2(i,j,k) = cs2(31,12,1) cs3(i,j,k) = cs3(31,12,1) cs4(i,j,k) = cs4(31,12,1) ss1(i,j,k) = ss1(31,12,1) ss2(i,j,k) = ss2(31,12,1) ss3(i,j,k) = ss3(31,12,1) ss4(i,j,k) = ss4(31,12,1) pm1(i,j,k) = pm1(31,12,1) pm2(i,j,k) = pm2(31,12,1) pm3(i,j,k) = pm3(31,12,1) pm4(i,j,k) = pm4(31,12,1) precip(i,j,k) = precip(31,12,1) drat(i,j,k) = drat(31,12,1) winde(i,j,k) = winde(31,12,1) else ttmp2 = tloss(i,j,k) stmp2 = sloss(i,j,k) pm10tmp2 = pm10loss(i,j,k) cs1tmp2 = cs1(i,j,k) cs2tmp2 = cs2(i,j,k) cs3tmp2 = cs3(i,j,k) cs4tmp2 = cs4(i,j,k) ss1tmp2 = ss1(i,j,k) ss2tmp2 = ss2(i,j,k) ss3tmp2 = ss3(i,j,k) ss4tmp2 = ss4(i,j,k) pm1tmp2 = pm1(i,j,k) pm2tmp2 = pm2(i,j,k) pm3tmp2 = pm3(i,j,k) pm4tmp2 = pm4(i,j,k) ptmp2 = precip(i,j,k) dtmp2 = drat(i,j,k) wtmp2 = winde(i,j,k) tloss(i,j,k) = ttmp1 sloss(i,j,k) = stmp1 pm10loss(i,j,k) = pm10tmp1 cs1(i,j,k) = cs1tmp1 cs2(i,j,k) = cs2tmp1 cs3(i,j,k) = cs3tmp1 cs4(i,j,k) = cs4tmp1 ss1(i,j,k) = ss1tmp1 ss2(i,j,k) = ss2tmp1 ss3(i,j,k) = ss3tmp1 ss4(i,j,k) = ss4tmp1 pm1(i,j,k) = pm1tmp1 pm2(i,j,k) = pm2tmp1 pm3(i,j,k) = pm3tmp1 pm4(i,j,k) = pm4tmp1 precip(i,j,k) = ptmp1 drat(i,j,k) = dtmp1 winde(i,j,k) = wtmp1 ttmp1 = ttmp2 stmp1 = stmp2 pm10tmp1 = pm10tmp2 cs1tmp1 = cs1tmp2 cs2tmp1 = cs2tmp2 cs3tmp1 = cs3tmp2 cs4tmp1 = cs4tmp2 ss1tmp1 = ss1tmp2 ss2tmp1 = ss2tmp2 ss3tmp1 = ss3tmp2 ss4tmp1 = ss4tmp2 pm1tmp1 = pm1tmp2 pm2tmp1 = pm2tmp2 pm3tmp1 = pm3tmp2 pm4tmp1 = pm4tmp2 ptmp1 = ptmp2 dtmp1 = dtmp2 wtmp1 = wtmp2 end if first = .false. end if 471 continue 461 continue 451 continue tloss(31,1,maxper) = tloss(31,1,maxper) - & (tloss(31,1,maxper)/(outyrs(maxper)+1)) sloss(31,1,maxper) = sloss(31,1,maxper) - & (sloss(31,1,maxper)/(outyrs(maxper)+1)) pm10loss(31,1,maxper) = pm10loss(31,1,maxper) - & (pm10loss(31,1,maxper)/(outyrs(maxper)+1)) cs1(31,1,maxper) = cs1(31,1,maxper) - (cs1(31,1,maxper)/ & (outyrs(maxper)+1)) cs2(31,1,maxper) = cs2(31,1,maxper) - (cs2(31,1,maxper)/ & (outyrs(maxper)+1)) cs3(31,1,maxper) = cs3(31,1,maxper) - (cs3(31,1,maxper)/ & (outyrs(maxper)+1)) cs4(31,1,maxper) = cs4(31,1,maxper) - (cs4(31,1,maxper)/ & (outyrs(maxper)+1)) ss1(31,1,maxper) = ss1(31,1,maxper) - (ss1(31,1,maxper)/ & (outyrs(maxper)+1)) ss2(31,1,maxper) = ss2(31,1,maxper) - (ss2(31,1,maxper)/ & (outyrs(maxper)+1)) ss3(31,1,maxper) = ss3(31,1,maxper) - (ss3(31,1,maxper)/ & (outyrs(maxper)+1)) ss4(31,1,maxper) = ss4(31,1,maxper) - (ss4(31,1,maxper)/ & (outyrs(maxper)+1)) pm1(31,1,maxper) = pm1(31,1,maxper) - (pm1(31,1,maxper)/ & (outyrs(maxper)+1)) pm2(31,1,maxper) = pm2(31,1,maxper) - (pm2(31,1,maxper)/ & (outyrs(maxper)+1)) pm3(31,1,maxper) = pm3(31,1,maxper) - (pm3(31,1,maxper)/ & (outyrs(maxper)+1)) pm4(31,1,maxper) = pm4(31,1,maxper) - (pm4(31,1,maxper)/ & (outyrs(maxper)+1)) precip(31,1,maxper) = precip(31,1,maxper) - & (precip(31,1,maxper)/(outyrs(maxper)+1)) drat(31,1,maxper) = drat(31,1,maxper) - & (drat(31,1,maxper)/(outyrs(maxper)+1)) winde(31,1,maxper) = winde(31,1,maxper) - & (winde(31,1,maxper)/(outyrs(maxper)+1)) c********************************************************************* do 450 k = 1,5 do 460 j = 1,13 do 470 i = 1,32 if (outday(i,j,k) .eq. 1) then if (i .eq. 32) then write(41,2520) k,tloss(i,j,k)/outyrs(k), & std(i,j,k), & sloss(i,j,k)/outyrs(k), & pm10loss(i,j,k)/outyrs(k), & cs1(i,j,k),cs2(i,j,k),cs3(i,j,k),cs4(i,j,k), & ss1(i,j,k),ss2(i,j,k),ss3(i,j,k),ss4(i,j,k), & pm1(i,j,k),pm2(i,j,k),pm3(i,j,k),pm4(i,j,k), & precip(i,j,k)/outyrs(k), & winde(i,j,k)/outyrs(k), & drat(i,j,k)/outyrs(k) else write(41,2500) i,j,k, c & opdir(i,j,k),opnam(i,j,k),opimp(i,j,k), & opdir(i,j,k),opnam(i,j,k), & tloss(i,j,k)/outyrs(k), & std(i,j,k), & sloss(i,j,k)/outyrs(k), & pm10loss(i,j,k)/outyrs(k), & cs1(i,j,k),cs2(i,j,k),cs3(i,j,k),cs4(i,j,k), & ss1(i,j,k),ss2(i,j,k),ss3(i,j,k),ss4(i,j,k), & pm1(i,j,k),pm2(i,j,k),pm3(i,j,k),pm4(i,j,k), & precip(i,j,k)/(outyrs(k)), & winde(i,j,k)/outyrs(k), & drat(i,j,k)/outyrs(k), & fldec(i,j,k)/outyrs(k), & stdec(i,j,k)/outyrs(k), & flmass(i,j,k)/outyrs(k), & stmass(i,j,k)/outyrs(k), & rdght(i,j,k)/outyrs(k), & rdgsp(i,j,k)/outyrs(k), & rr(i,j,k)/outyrs(k) end if end if 470 continue 460 continue 450 continue do 480 j = 1,12 write(41,2540) j, tloss(32,j,6)/simyrs,std(32,j,6), & sloss(32,j,6), & pm10loss(32,j,6), & cs1(32,j,6),cs2(32,j,6),cs3(32,j,6),cs4(32,j,6), & ss1(32,j,6),ss2(32,j,6),ss3(32,j,6),ss4(32,j,6), & pm1(32,j,6),pm2(32,j,6),pm3(32,j,6),pm4(32,j,6), & precip(32,j,6)/simyrs,winde(32,j,6)/simyrs, & drat(32,j,6)/simyrs, & fldec(32,j,6)/simyrs,stdec(32,j,6)/simyrs, & flmass(32,j,6)/simyrs,stmass(32,j,6)/simyrs, & rdght(32,j,6)/simyrs,rdgsp(32,j,6)/simyrs, & rr(32,j,6)/simyrs 480 continue write(41,2560) tloss(32,13,6)/simyrs,std(32,13,6), & sloss(32,13,6)/simyrs, pm10loss(32,13,6)/simyrs, & cs1(32,13,6),cs2(32,13,6),cs3(32,13,6),cs4(32,13,6), & ss1(32,13,6),ss2(32,13,6),ss3(32,13,6),ss4(32,13,6), & pm1(32,13,6),pm2(32,13,6),pm3(32,13,6),pm4(32,13,6), & precip(32,13,6)/simyrs,winde(32,13,6)/simyrs, & drat(32,13,6)/simyrs 500 stop 'The WEPS simulation run is finished' end