c$Author: fredfox $ c$Date: 2001-08-15 16:42:47 $ c$Revision: 1.43.2.4 $ c$Source: /weru/cvs/weps/weps.src/main/weps.for,v $ 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 C EDIT HISTORY C 05-Feb-99 wjr changed all open's to calls to fopenk C 19-Feb-99 wjr changed call to mfinit removing unit # C 19-Feb-99 wjr changed call to manage removing unit # C 03-Mar-99 wjr changed calls to hydro & soil to wrapper calls C 04-Mar-99 wjr changed call to crop to wrapper call C 04-Mar-99 wjr put plot and inithydr into separate files & C removed all hydro includes C 10-Mar-99 wjr did lots of stuff -- put init code into separate C files, put total loops into sep file, open into ... C created soilinit, cliginit, ... C c#include c external continue_hdl c external common_handler c + + + GLOBAL COMMON BLOCKS + + + include 'p1werm.inc' include 'wpath.inc' include 'p1unconv.inc' include 'm1subr.inc' include 'm1sim.inc' include 'm1geo.inc' include 'm1flag.inc' include 'm1dbug.inc' include 's1layr.inc' include 's1layd.inc' include 's1sgeo.inc' include 's1phys.inc' include 'c1info.inc' include 'c1gen.inc' include 'w1clig.inc' include 'w1wind.inc' include 'w1pavg.inc' include 'file.fi' include 'b1glob.inc' include 'd1glob.inc' include 'c1glob.inc' include 'c1db1.inc' include 'h1hydro.inc' include 'timer.fi' C *** remove after taking out debugging hack include 's1dbh.inc' include 'h1db1.inc' include 's1surf.inc' include 's1agg.inc' C *** 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 'decomp/decomp.inc' c + + + LOCAL VARIABLES + + + character *(64) bd_date !declarations for build date function character *(64) rel_ver !declarations for release/version function logical first C integer cd, cm, cy, ccd, ccm, ccy, cwd, cwm, cwy, integer cd, cm, cy, i ndiy, C *** i i, isr, ioc, l, j, k, i isr, index, jindex, kindex, i o_unit, i simyrs, i yrsim 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 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 daysim - This variable holds the total current days of simulation. 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 isr - This variable holds the subregion index. c iy - Starting year of simulation run (used in management). c index - This variable is a general index variable to access arrays 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 period - 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 + + + 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 erodinit - Erosion initialization routines c erosion - 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 initialization subroutines c sdbug - Prints global variables before and after call to SOIL c soil - Soil submodel c asdini - aggregate size distribution initialization c bpools - prints many biomass pool components (for debugging) 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 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 = ? - 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 = 'erod.out' - erosion output file o_unit = 20 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 32 = 'plot.out' - plotting output file c 40 = temporary file to hold accounting region erosion values c + + + DATA INITIALIZATIONS + + + C *** data xmin /0.0/, xmax /0.0/, ymin /0.0/, ymax /0.0/ data yrsim /0/ data first /.true./ c + + + INPUT FORMATS + + + c + + + OUTPUT FORMATS + + + c format for header of plot file 2050 format (1x,'day ','|','mon ','|',' yr ','|',' operation ', & '|',' tot_loss ','|',' suspen ','|',' pm10 ','|', & ' max_wind ','|',' dir_wind ','|',' precip ','|','Surf_H2O', & ' |',' ridge_ht ','|',' ridge_or ','|',' r_rough ','|', & ' gmd ','|',' cr_fract ','|',' ag_stab ','|',' lai ', & '|',' sai ','|','crop_yield','|',' st_mass ','|', & ' can_cov ', & '|',' fl_cov%','|',' st_cov% ','|', & 'loose_mass','|','loose_frac','|',' bulk_den') c format for header of erosion grid file 2055 format (1x,' d',1x,' mo',1x,' yr',7x,'dir',5x,'i',3x,'j', & 5x,'total',6x,'slt+crp',4x,'susp') c format for header of output file 2060 format (1x,'d',1x,'mo',1x,'yr|',' operation', 31x, & '| tot_loss|',4x,'st_dev|',2x,'crp+salt|', & 3x,'cs_loss|',4x,'cs_dep|', & 'cslos_area|','csdep_area|',4x, & 'suspen|', & 6x,'pm10|',7x,'cs1|',7x,'cs2|',7x,'cs3|',7x,'cs4|', & 7x,'ss1|',7x,'ss2|',7x,'ss3|',7x,'ss4|', & 7x,'pm1|',7x,'pm2|',7x,'pm3|',7x,'pm4|', & 4x,'precip|',2x,'w_energy|',3x,'dry_idx|',1x,'l_can_cov|', & 'l_sil_area|',1x,'l_st_mass|',2x,'d_fl_cov|', & 2x,'d_st_sil|',1x,'d_fl_mass|',1x,'d_st_mass|', & 'b_f_fl_cov|','b_f_st_sil|','b_m_fl_cov|', & 'b_m_st_sil|',4x,'rdg_or|',4x,'rdg_ht|',4x,'rdg_sp|', & 8x,'rr') c + + + END SPECIFICATIONS + + + C Before anything else is done, have WEPS output to stderr(?) C the build date and release/version information. write(6,*) write(6,*) 'WEPS 1.0 Version: ', trim(rel_ver()), & ' Built on date: ', trim (bd_date()) write(6,*) call timer(0,TIMSTART) call fopenk(luolog,'logfil.txt','unknown') c open input files and read run files call input C c temporarily initialize old random roughness aslrrc(1) = 10. as0rrk(1) = 0.9 call openfils c Initialize the management file and rotation counters call mfinit(1, tinfil, maxper) call asdini() call erodinit am0gdf = .true. C *** debugging write c write(*,*) ' weps 1: ', asdblk(1:7,1) C *** eodw operat =' ' if((am0hfl.gt.0).or.(am0sfl.gt.0).or.(am0tfl.gt.0) & .or.(am0cfl.gt.0).or.(am0dfl.gt.0).or.(am0efl.gt.0)) then if (am0ifl .eqv. .true.) write (luoplt, 2050) endif c if (am0ifl .eqv. .true.) write (7, 2055) if (am0ifl .eqv. .true.) write (41, 2060) 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). c these variables were used before they were given a value c and did not have another place where initialization was appropriate do 10 isr=1,nsubr call decoinit(isr) abdstm(isr) = 0.0 abffcv(isr) = 0.0 abfscv(isr) = 0.0 acxstm(isr) = 0.0 acxstmrep(isr) = 0.0 acrbc(isr) = 1 accovfact(isr) = 0.0 do index=1,mnbpls adrbc(index,isr) = 1 end do c Initialize the water holding capacity variable call hydrinit(isr) c initialize soil depth to bottom of layers (mm) from layer thickness (mm) c and initialize applied NO3 call soilinit(isr) 10 continue c just another needed initialization do index = 1,mnopdat do jindex = 1,mnmon do kindex = 1,mnryr outday(index, jindex, kindex) = 0 end do end do end do c read cligen header - remove when user interface is installed call cliginit 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 caldatw (cd, cm, cy) c determine number of days in the year ndiy = 365 if (isleap(cy) .eqv. .true.) ndiy = 366 call getcli(cd, cm, cy) C call getwin(cd, cm, cy) C c print current day of simulation to screen periodically daysim = daysim + 1 if ((cm .eq. 1) .and. (cd .eq. 1)) then yrsim = yrsim + 1 simyrs = (ly - iy + 1) - maxper if(am0sif) write(*,*) 'Initializing Surface' write(*,*) 'Year', yrsim, ' of', (ly - iy + 1) end if c write(*,*) 'weps:yr,mon,day',yrsim,cm,cd c do subregions do isr = 1, nsubr c if (am0jd.eq.ijday+1) call dbgdmp(daysim, isr) c if (am0jd.eq.ljday) call dbgdmp(daysim, isr) am0csr = isr c write(*,*) "Start callhydr" call callhydr(daysim, isr) !call HYDROLOGY submodel c write(*,*) "Start manage" call manage (isr, cd, cm, cy,iy,lopday,lopmon,lopyr) !MANAGEment (tillage) submodel c write(*,*) "Start updres" call updres(isr) !update decomp residue pools c write(*,*) "Start callsoil" call callsoil(daysim, isr) !SOIL submodel c write(*,*) "Start callcrop" if (am0cgf) call callcrop(isr) !CROP submodel if between planting and killing harvest c write(*,*) "Start decomp" call decomp(isr) !DECOMPosition submodel c write(*,*) "Start updres" call updres(isr) c write(*,*) "Start sumbio" call sumbio(isr) ! sum live and dead biomass end do ! end of sub-region loop c skip next section until surface is initialized if (.not.am0sif) then c if wind is great enough, call erosion if (awudmx .gt. 8.) then c write(*,*) "Start calcwu" call calcwu c write(*,*) "Start erosion" call erosion (o_unit ) end if c do bookkeeping for output of general report c write(*,*) "Start bookeep" call bookeep (cd, cm, cy) end if c set initialization flag to .false. after first day am0ifl = .false. if( (am0hfl.gt.0).or.(am0sfl.gt.0).or.(am0tfl.gt.0) & .or.(am0cfl.gt.0).or.(am0dfl.gt.0).or.(am0efl.gt.0) ) then c make operation name available for plot output if (am0sif .eqv. .true.) then if ((lopday .eq. cd) .and. (lopmon .eq. cm) .and. & (lopyr .eq. amnryr(1))) then operat = opname else operat = & ' ' end if total = 0.0 suspen = 0.0 pmten = 0.0 end if c print to plot data file call plotdata endif c write decomposition biomass pool amounts to files if (am0dfl.gt.0) call bpools(cd,cm,cy,1) c check for end of first rotation (surface initialization) 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 C close (luiwin) close (luicli) close (luiwsd) close (unit = 40) close (luoplt) c output the simulation general report c call genrep c call wsum c call caltot(maxper, simyrs, first) call caltot(simyrs, first) write (*,*) 'The WEPS simulation run is finished' call timer(0,TIMSTOP) call timer(0,TIMPRINT) stop 'The WEPS simulation run is finished ' end