!$Author$
!$Date$
!$Revision$
!$HeadURL$
module sae_in_out_mod
implicit none
! This module is for the creation of single or multiple stand alone erosion input files,
! depending on the command line switches given
type make_sae_in_out
integer :: jday ! the present day in julian days for output
integer :: simday ! the present day in simulation days for creation of file
integer :: maxday ! maximum simday number
character*256 :: fullpath ! the root path plus subdirectory if indicated by multiple files
end type make_sae_in_out
type(make_sae_in_out) :: mksaeinp
type(make_sae_in_out) :: mksaeout
! placed here for sharing back with hagen_plot_flag by daily_erodout
real :: aegt, aegtss, aegt10, aegt2_5
logical :: in_weps
contains
subroutine saeinp( luo_saeinp, subrsurf )
! +++ PURPOSE +++
! print out input file for stand alone erosion
! + + + Modules Used + + +
use datetime_mod, only: caldat
use file_io_mod, only: fopenk, makenamnum
use grid_mod, only: amxsim, amasim, xgdpt, ygdpt
use subregions_mod
use barriers_mod, only: barrier
use erosion_data_struct_defs, only: subregionsurfacestate, awzypt, awdair, anemht, awzzo, wzoflg, awadir, subday, ntstep
! +++ ARGUMENT DECLARATIONS +++
integer, intent(inout) :: luo_saeinp ! output unit number
type(subregionsurfacestate), dimension(:), intent(in) :: subrsurf ! subregion surface conditions (erosion specific set)
! +++ LOCAL VARIABLES +++
integer k,l, sr, ip
integer b, nbr
integer day, mon, yr
! + + + LOCAL VARIABLE DEFINITIONS + + +
! sr - index used in subregion loop
! ip - index to polygon coordinates
! +++ END SPECIFICATIONS +++
if( luo_saeinp .lt. 0 ) then
call fopenk (luo_saeinp, trim(mksaeinp%fullpath) // makenamnum('saeros', mksaeinp%simday, mksaeinp%maxday, '.in'),'unknown')
call caldat (mksaeinp%jday,day,mon,yr)
write(*,'(4(a,i0))') 'Made SWEEP input file D/M/Y: ', day,'/', mon,'/', yr,' simulation day: ', mksaeinp%simday
write(luo_saeinp,*) 'Version: 100 multiple subregions'
write(luo_saeinp,2101) day, mon, yr
2101 format('# WEPS erosion day mon yr',2(1x,i2),2x,i4)
else
write(luo_saeinp,*) ' REPORT OF INPUTS (read by erodin.for) '
end if
! print header info
write(luo_saeinp,1005)
1005 format ( &
'#',65('*'),/, &
'#',/, &
'# +++ PURPOSE +++',/, &
'#',/, &
'# File for input to standalone erosion submodel program (sweep)',/, &
'#',/, &
'# All lines beginning with a "#" character are assumed to',/, &
'# be comment lines and are skipped.',/, &
'#',/, &
'# +++ DEFINITIONS +++',/, &
'#',/, &
'# All comments prior to each line of data input',/, &
'# in this template input file have the following format:',/, &
'#',/, &
'# Variable_Name, Var_type, Text Definition',/, &
'#',/, &
'# where Var_type is: I = integer L = logical R = real',/, &
'#',/, &
'#',/, &
'# +++ DEBUG FLAG +++',/, &
'#',/, &
'# debugflg - debug flag for providing different levels of debug info',/, &
'# currently useful to debug/check input file data format',/, &
'#',/, &
'# value of 0 will print no debug information',/, &
'# value of 1 will print out and number all input file lines',/, &
'# value of 2 will print out and number all data input lines',/, &
'# value of 3 will do both 1 and 2')
write(luo_saeinp,*) '0'
write(luo_saeinp,2000)
2000 format( &
'#',/, &
'#',/, &
'# +++ INITIALIZATION +++',/, &
'#',/, &
'# am0eif, L, EROSION "initialization" flag',/, &
'# Must be set to .TRUE. for standalone erosion runs')
write(luo_saeinp,*) '.TRUE.'
write(luo_saeinp,2100)
2100 format ( &
'#',/, &
'#',/, &
'# +++ SIMULATION REGION +++',/, &
'#',/, &
'# amxsim(x,y), R, Simulation Region diagonal coordinates (meters)',/, &
'# Input (x,y) coordinates in this form: x1,y1 x2,y2',/, &
'# If orientation in 0 degrees, then x axis is East-West',/, &
'# and y axis is North-South.',/, &
'# Typical Range: 10.0 to 1600.0',/, &
'#',/, &
'# NOTE: Accounting region and Subregion coordinates',/, &
'# must also be set to the same values',/, &
'#')
write(luo_saeinp,*) amxsim(1)%x, amxsim(1)%y, amxsim(2)%x, amxsim(2)%y
write(luo_saeinp,2105)
2105 format( &
'#',/, &
'#',/, &
'# amasim, R, Simulation Region orientation angle (degrees from North)')
write(luo_saeinp,*) amasim
write(luo_saeinp,2107)
2107 format( &
'#',/, &
'# +++ NUMBER OF GRID POINTS +++',/, &
'#',/, &
'# xgdpt, ygdpt, I, Number of grid points in simulation region x and y directions',/, &
'# A value of 0 for both values uses the default internal grid.',/, &
'# A value of 0 for only one generates and error.')
write(luo_saeinp,*) xgdpt, ygdpt
write(luo_saeinp,2110)
2110 format('#',/, &
'# +++ ACCOUNTING REGIONS +++',/, &
'#',/, &
'# nacctr, I, Number of accounting regions')
write(luo_saeinp,*) size(acct_poly)
! loop through all accounting regions
do sr = 1, size(acct_poly)
write(luo_saeinp,2115)
2115 format('#',/, &
'# accounting region polygon, count, xy pairs (subregions_mod)',/,&
'# polygons can be:',/, &
'# - open (the last point connects to the first point)',/,&
'# - closed (last point is same as first point)',/, &
'# - complex (multiple closed polygons entered as one)')
! number of coordinate pairs in polygon
write(luo_saeinp,*) acct_poly(sr)%np
! the coordinate pairs
do ip = 1, acct_poly(sr)%np
write(luo_saeinp,*) acct_poly(sr)%points(ip)%x, acct_poly(sr)%points(ip)%y
end do
end do
! barriers
write(luo_saeinp,2120)
2120 format ('#',/, &
'# +++ BARRIERS +++',/, &
'#',/, &
'# nbr, I, Number of barriers (0-5) ')
nbr = size(barrier)
write(luo_saeinp,*) nbr
write(luo_saeinp,2122)
2122 format ('#',/, &
'# NOTE: Remaining BARRIER inputs are repeated for each barrier specified',/, &
'# If no barriers specified (nbr=0), then no BARRIER inputs will be here',/, &
'#',/, &
'# barrier(b)%np, I, number of points in barrier polyline',/, &
'# Inputs are repeated for each point specified',/, &
'# barrier(b)%points(n)%x, barrier(b)%points(n)%y, R, x,y coordinate pair',/, &
'# Example: 0.0 500.0 ',/, &
'# barrier(b)%points(n)%amzbr, R, Barrier height (m)',/, &
'# barrier(b)%points(n)%amxbrw, R, Barrier width (m)',/, &
'# barrier(b)%points(n)%ampbr, R, Barrier porosity (m^2/m^2)',/, &
'# Example: 1.2 2.0 0.50 ',/, &
'#')
do b = 1, nbr
write(luo_saeinp,2125)
2125 format('# ',/, &
'# barrier(b)%np, I, number of points in barrier polyline')
write(luo_saeinp,*) barrier(b)%np
do ip = 1, barrier(b)%np
write(luo_saeinp,2127)
2127 format('# barrier(b)%points(n)%x, barrier(b)%points(n)%y, R, x,y coordinate pair')
write(luo_saeinp,*) barrier(b)%points(ip)%x, barrier(b)%points(ip)%y
write (luo_saeinp,2130)
2130 format('# amzbr(b), amxbrw(b), ampbr(b)')
write (luo_saeinp,*) barrier(b)%param(ip)%amzbr, barrier(b)%param(ip)%amxbrw, barrier(b)%param(ip)%ampbr
end do
end do
! subregions
write(luo_saeinp,2135)
2135 format( &
'#',/, &
'# +++ SUBREGIONS +++',/, &
'#',/, &
'# nsubr, I, Number of subregions (at least 1)',/, &
'# NOTE: together all must cover simulation region')
write(luo_saeinp,*) size(subr_poly)
write(luo_saeinp,2137)
2137 format( &
'#',/, &
'# NOTE: Remaining SUBREGION inputs (BIOMASS, SOIL, and HYDROLOGY,',/, &
'# ie. variables defined by subregion) are repeated for nsubr',/, &
'# subregions specified')
! loop through all subregions
do sr = 1, size(subr_poly)
write(luo_saeinp,2138)
2138 format( &
'#',/, &
'# subregion polygon, count, xy pairs (subregions_mod)',/, &
'# polygons can be:',/, &
'# - open (the last point connects to the first point)',/,&
'# - closed (last point is same as first point)',/, &
'# - complex (multiple closed polygons entered as one)')
! number of coordinate pairs in polygon
write(luo_saeinp,*) subr_poly(sr)%np
! the coordinate pairs
do ip = 1, subr_poly(sr)%np
write(luo_saeinp,*) subr_poly(sr)%points(ip)%x, subr_poly(sr)%points(ip)%y
end do
write(luo_saeinp,2140)
2140 format('#',/, &
'# +++ BIOMASS +++',/, &
'#',/, &
'# subrsurf(s)adzht_ave, R, Average residue height (m)')
write(luo_saeinp,*) subrsurf(sr)%adzht_ave
! Changed above to use "average residue height" instead of "overall height" - LEW 1/26/06
! '# biotot(s)%zht_ave, R, Overall biomass height (m)')
! write(luo_saeinp,*) biotot(s)%zht_ave
write(luo_saeinp,2146)
2146 format('#',/, &
'# aczht(s), R, Crop height (m)')
write(luo_saeinp,*) subrsurf(sr)%aczht
write(luo_saeinp,2147)
2147 format('#',/, &
'# acrsai(s), R, Crop stem area index (m^2/m^2)',/, &
'# acrlai(s), R, Crop leaf area index (m^2/m^2)')
write(luo_saeinp,*) subrsurf(sr)%acrsai, subrsurf(sr)%acrlai
write(luo_saeinp,2148)
2148 format('#',/, &
'# subrsurf(s)%adrsaitot, R, Residue stem area index (m^2/m^2)',/, &
'# subrsurf(s)%adrlaitot, R, Residue leaf area index (m^2/m^2)')
write(luo_saeinp,*) subrsurf(sr)%adrsaitot, subrsurf(sr)%adrlaitot
write(luo_saeinp,2149)
2149 format('#',/, &
'# acxrow(s) Crop row spacing (m)',/, &
'# ac0rg(s) Crop seed placement (0 - furrow, 1 - ridge)')
write(luo_saeinp,*) subrsurf(sr)%acxrow, subrsurf(sr)%ac0rg
write(luo_saeinp,2150)
2150 format('#',/, &
'# These are not implemented within EROSION',/, &
'# abrsaz(h,s), R, (b1geom.inc) Biomass stem area index by ht (1/m)',/, &
'# (should be 5 values here when used)',/, &
'# abrlaz(h,s), R, (b1geom.inc) Biomass leaf area index by ht (1/m)',/, &
'# (should be 5 values here when used)',/, &
'#',/, &
'# Only abffcv(s) is currently implemented within EROSION',/, &
'# abffcv(s), R, (b1geom.inc) Flat biomass cover (m^2/m^2)',/, &
'# abfscv(s), R, (b1geom.inc) Standing biomass cover (m^2/m^2)',/, &
'# abftcv(s), R, (b1geom.inc) Total biomass cover (m^2/m^2)',/, &
'# (should be 3 values here when abffcv(s) and abfscv(s) are used)')
write(luo_saeinp,*) subrsurf(sr)%abffcv
! soil
write(luo_saeinp,2160)
2160 format('#',/, &
'# +++ SOIL +++',/, &
'#',/, &
'# nslay(s), I, (s1layr.inc) Number of soil layers (3-10)')
write(luo_saeinp,*) subrsurf(sr)%nslay
write(luo_saeinp,2165)
2165 format('#',/, &
'# NOTE: Remaining SOIL inputs are repeated for each layer specified',/, &
'#',/, &
'# aszlyt(l,s), R, (s1layr.inc) Soil layer thickness (mm)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%aszlyt, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2170)
2170 format('#',/, &
'# asdblk(l,s), R, (s1phys.inc) Soil layer bulk density (Mg/m^3')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%asdblk, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2175)
2175 format('#',/, &
'# asfsan(l,s),R,(s1dbh.inc) Soil layer sand content (Mg/Mg)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%asfsan, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2177)
2177 format('#',/, &
'# asfvfs(l,s), R, (s1dbh.inc) Soil layer very fine sand (Mg/Mg)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%asfvfs, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2180)
2180 format('#',/, &
'# asfsil(l,s),R,(s1dbh.inc) Soil layer silt content (Mg/Mg)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%asfsil, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2185)
2185 format('#',/, &
'# asfcla(l,s),R,(s1dbh.inc) Soil layer clay content (Mg/Mg)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%asfcla, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2190)
2190 format('#',/, &
'# asvroc(l,s), R, (s1dbh.inc) Soil layer rock volume (m^3/m^3)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%asvroc, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2195)
2195 format('#',/, &
'# asdagd(l,s),R,(s1agg.inc) Soil layer agg density (Mg/m^3)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%asdagd, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2200)
2200 format('#',/, &
'# aseags(l,s), R, (s1agg.inc) Soil layer agg stability ln(J/kg)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%aseags, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2205)
2205 format('#',/, &
'# aslagm(l,s), R, (s1agg.inc) Soil layer GMD (mm)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%aslagm, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2210)
2210 format('#',/, &
'# aslagn(l,s), R, (s1agg.inc) Soil layer minimum agg size (mm)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%aslagn, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2215)
2215 format('#',/, &
'# aslagx(l,s), R, (s1agg.inc) Soil layer maximum agg size (mm)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%aslagx, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2220)
2220 format('#',/, &
'# as0ags(l,s), R, (s1agg.inc) Soil layer GSD (mm/mm)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%as0ags, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2225)
2225 format('#',/, &
'# asfcr(s), R, Surface crust fraction (m^2/m^2)',/, &
'# aszcr(s), R, Surface crust thickness (mm)',/, &
'# asflos(s), R, Fraction of loose material on surface (m^2/m^2)',/, &
'# asmlos(s), R, Mass of loose material on crust (kg/m^2)',/, &
'# asdcr(s), R, Soil crust density (Mg/m^3)',/, &
'# asecr(s), R, Soil crust stability ln(J/kg)')
write(luo_saeinp,*) subrsurf(sr)%asfcr, subrsurf(sr)%aszcr, subrsurf(sr)%asflos, subrsurf(sr)%asmlos, &
subrsurf(sr)%asdcr, subrsurf(sr)%asecr
write(luo_saeinp,2230)
2230 format('#',/, &
'# aslrr(s), R, Allmaras random roughness (mm)')
write(luo_saeinp,*) subrsurf(sr)%aslrr
write(luo_saeinp,2235)
2235 format('#',/, &
'# aszrgh(s), R, Ridge height (mm)',/, &
'# asxrgs(s), R, Ridge spacing (mm)',/, &
'# asxrgw(s), R, Ridge width (mm)',/, &
'# asargo(s), R, Ridge orientation (deg)')
write(luo_saeinp,*) subrsurf(sr)%aszrgh, subrsurf(sr)%asxrgs, subrsurf(sr)%asxrgw, subrsurf(sr)%asargo
write(luo_saeinp,2240)
2240 format('#',/, &
'# asxdks(s), R, Dike spacing (mm)')
write(luo_saeinp,*) subrsurf(sr)%asxdks
write(luo_saeinp,2245)
! hydrology
2245 format('#',/, &
'# +++ HYDROLOGY +++',/, &
'#',/, &
'# ahzsnd(s), R, Snow depth (mm)')
write(luo_saeinp,*) subrsurf(sr)%ahzsnd
write(luo_saeinp,2250)
2250 format('#',/, &
'# ahrwcw(l,s), R, (h1db1.inc) Soil layer wilting point water content (Mg/Mg)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%ahrwcw, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2255)
2255 format('#',/, &
'# ahrwca(l,s), R, (h1db1.inc) Soil layer water content (Mg/Mg)')
write(luo_saeinp,*) (subrsurf(sr)%bsl(l)%ahrwca, l=1,subrsurf(sr)%nslay)
write(luo_saeinp,2260)
2260 format('#',/, &
'# ahrwc0(h,s), R, (h1db1.inc) Surface layer water content (Mg/Mg)',/, &
'# NOTE: the near surface water content is specified on an',/, &
'# hourly basis. We read in the hrly water content',/, &
'# on two lines, with 12 values in each line.')
write(luo_saeinp,22) (subrsurf(sr)%ahrwc0(l), l=1,12)
22 format(12(1x,f10.7))
write(luo_saeinp,22) (subrsurf(sr)%ahrwc0(l), l=13,24)
! end of subregion loop
end do
! weather
write(luo_saeinp,2270)
2270 format('#',/, &
'# NOTE: This is the end of the SUBREGION variables',/, &
'#',/, &
'# +++ WEATHER +++',/, &
'#',/, &
'# awzypt, R, Average annual precipitation (mm)')
write(luo_saeinp,*) awzypt
write(luo_saeinp,2273)
2273 format('#',/, &
'# awdair, R, Air density (kg/m^3)')
write(luo_saeinp,*) awdair
write(luo_saeinp,2275)
2275 format('#',/, &
'# awadir, R, Wind direction (deg)')
write(luo_saeinp,*) awadir
write(luo_saeinp,2280)
2280 format('#',/, &
'# ntstep, I, (local variable) Number of intervals/day to run EROSION')
write(luo_saeinp,*) ntstep
write(luo_saeinp,2285)
2285 format('#',/, &
'# anemht, R anemometer height (m)',/, &
'# awzzo, R aerodynamic roughness at anemometer site (mm)', /, &
'# wzoflg, I (global variable) zo location flag',/, &
'# (flag =0 - zo fixed at wx sta. location)',/, &
'# (flag = 1 - zo variable at field location)')
write(luo_saeinp,*) anemht, awzzo, wzoflg
write(luo_saeinp,2290)
2290 format('#',/, &
'# wflg, I, (local variable) Wind/Weibull flag',/, &
'# (0 - read in Weibull parameters, 1 - read in wind speeds)')
write(luo_saeinp,*) '1'
write(luo_saeinp,2295)
2295 format('#',/, &
'# NOTE: This is only present when the above (wflg=0)',/, &
'# wfcalm, R, (local variable) Fraction of time winds are calm (hr/hr)',/, &
'# wuc, R, (local variable) Weibull "c" factor (m/s)',/, &
'# w0k, R, (local variable) Weibull "k" factor (fraction)',/, &
'# 0.263 5.856 1.720', &
'#',/, &
'# NOTE: The remaining data is only present when (wflg=1)',/, &
'# wflg=1 uses standard input from windgen in WEPS.',/, &
'#',/, &
'# awu(i), R, (w1wind) Wind speed for (ntstep) intervals (m/s)',/, &
'#',/, &
'# I think I can read multiple lines with variable number of values',/, &
'# We will try and see - LEW Must use 6 values per line LH.',/, &
'#')
write(luo_saeinp,*) (subday(k)%awu, k=1,6)
write(luo_saeinp,*) (subday(k)%awu, k=7,12)
write(luo_saeinp,*) (subday(k)%awu, k=13,18)
write(luo_saeinp,*) (subday(k)%awu, k=19,24)
write(luo_saeinp,2300)
2300 format( '#' )
close(luo_saeinp)
end subroutine saeinp
subroutine daily_erodout( o_unit, o_E_unit, sgrd_u, input_filename, cellstate )
! +++ PURPOSE +++
! To print output desired from standalone EROSION submodel
use file_io_mod, only: fopenk, makenamnum
use datetime_mod, only: get_systime_string, caldat
use erosion_data_struct_defs, only: cellsurfacestate, am0efl
use grid_mod, only: imax, jmax, amasim, amxsim
! +++ ARGUMENT DECLARATIONS +++
integer, intent(inout) :: o_unit, o_E_unit, sgrd_u
character(len=*), intent(in) :: input_filename
type(cellsurfacestate), dimension(0:,0:), intent(in) :: cellstate ! initialized grid cell state values
! ++++ LOCAL VARIABLES +++
character(len=21) :: rundatetime
integer i, j
real tt, lx, ly
real topt,topss, top10, top2_5, bott, botss, bot10, bot2_5
real ritt, ritss, rit10, rit2_5, lftt, lftss, lft10, lft2_5
real tot, totbnd
integer yr, mon, day
! +++ END SPECIFICATIONS +++
! Calculate Averages Crossing Borders
! top border
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
aegt2_5 = 0.0
j = jmax
do 1 i = 1, imax-1
aegt = aegt + cellstate(i,j)%egt
aegtss = aegtss + cellstate(i,j)%egtss
aegt10 = aegt10 + cellstate(i,j)%egt10
aegt2_5 = aegt2_5 + cellstate(i,j)%egt2_5
1 continue
! calc. average at top border
topt = aegt/(imax-1)
topss = aegtss/(imax-1)
top10 = aegt10/(imax-1)
top2_5 = aegt2_5/(imax-1)
! bottom border
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
aegt2_5 = 0.0
j = 0
do 2 i = 1, imax-1
aegt = aegt + cellstate(i,j)%egt
aegtss = aegtss + cellstate(i,j)%egtss
aegt10 = aegt10 + cellstate(i,j)%egt10
aegt2_5 = aegt2_5 + cellstate(i,j)%egt2_5
2 continue
! calc. average at bottom border
bott = aegt/(imax-1)
botss = aegtss/(imax-1)
bot10 = aegt10/(imax-1)
bot2_5 = aegt2_5/(imax-1)
! right border
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
aegt2_5 = 0.0
i = imax
do 3 j = 1, jmax-1
aegt = aegt + cellstate(i,j)%egt
aegtss = aegtss + cellstate(i,j)%egtss
aegt10 = aegt10 + cellstate(i,j)%egt10
aegt2_5 = aegt2_5 + cellstate(i,j)%egt2_5
3 continue
! calc. average at right border
ritt = aegt/(jmax-1)
ritss = aegtss/(jmax-1)
rit10 = aegt10/(jmax-1)
rit2_5 = aegt2_5/(jmax-1)
! left border
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
aegt2_5 = 0.0
i = 0
do 4 j = 1, jmax-1
aegt = aegt + cellstate(i,j)%egt
aegtss = aegtss + cellstate(i,j)%egtss
aegt10 = aegt10 + cellstate(i,j)%egt10
aegt2_5 = aegt2_5 + cellstate(i,j)%egt2_5
4 continue
! calc. average at left border
lftt = aegt/(jmax-1)
lftss = aegtss/(jmax-1)
lft10 = aegt10/(jmax-1)
lft2_5 = aegt2_5/(jmax-1)
! calculate averages of inner grid points
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
do 5 j=1,jmax-1
do 5 i= 1, imax-1
aegt= aegt + cellstate(i,j)%egt
aegtss = aegtss + cellstate(i,j)%egtss
aegt10 = aegt10 + cellstate(i,j)%egt10
aegt2_5 = aegt2_5 + cellstate(i,j)%egt2_5
5 continue
tt = (imax-1)*(jmax-1)
aegt = aegt/tt
aegtss = aegtss/tt
aegt10 = aegt10/tt
aegt2_5 = aegt2_5/tt
! calculate comparison of boundary and interior losses
lx = amxsim(2)%x - amxsim(1)%x
ly = amxsim(2)%y - amxsim(1)%y
tot = aegt*lx*ly
totbnd = (topt + bott + topss + botss)*lx + (ritt + lftt + ritss + lftss)*ly
if (btest(am0efl,1)) then
if( o_unit .lt. 0 ) then
call fopenk (o_unit, trim(mksaeout%fullpath) // makenamnum('saeros', mksaeout%simday, mksaeout%maxday, '.egrd'),'unknown')
call caldat (mksaeout%jday,day,mon,yr)
write(*,'(4(a,i0))') 'Made Daily Erosion grid file for: ', day,'/', mon,'/', yr,' simulation day: ', mksaeout%simday
write(o_unit, "('# WEPS erosion day mon yr',2(1x,i2),2x,i4)") day, mon, yr
write (o_unit,*)
write (o_unit,*) 'Grid cell output from WEPS run'
write (o_unit,*)
else
! write header to files
write (o_unit,*)
write (o_unit,*)
write (o_unit,*) 'Grid cell output from SWEEP run'
write (o_unit,*)
end if
! Print date of Run
rundatetime = get_systime_string() ! with Lahey f95, had to assign to variable first
write(o_unit,"(1x,'Date of run: ',a21)") rundatetime
write(o_unit,*)
write(o_unit,fmt="(1x,a)") ""
write(o_unit,fmt="(1x,5f10.2)") amasim, amxsim(1)%x, amxsim(1)%y, amxsim(2)%x, amxsim(2)%y
write(o_unit,fmt="(1x,a)") ""
write(o_unit,*)
write (o_unit,*) 'Total grid size: (', imax+1,',', jmax+1, ') ',&
'Inner grid size: (', imax-1,',', jmax-1, ')'
write (o_unit,*)
write (o_unit,6)
write (o_unit,*) ' top(i=1,imax-1,j=jmax) ', &
'bottom(i=1,imax-1,j=0) ', &
'right(i=imax,j=1,jmax-1) ', &
'left(i=0, j=1,jmax-1) '
write (o_unit,10) (cellstate(i,jmax)%egt+cellstate(i,jmax)%egtss, i = 1, imax-1)
write (o_unit,10) (cellstate(i,0)%egt+cellstate(i,0)%egtss, i = 1, imax-1)
write (o_unit,10) (cellstate(imax,j)%egt+cellstate(imax,j)%egtss, j = 1, jmax-1)
write (o_unit,10) (cellstate(0,j)%egt+cellstate(0,j)%egtss, j = 1, jmax-1)
write (o_unit,*)
write (o_unit,7)
write (o_unit,*) ' top(i=1,imax-1,j=jmax) ', &
'bottom(i=1,imax-1,j=0) ', &
'right(i=imax,j=1,jmax-1) ', &
'left(i=0, j=1,jmax-1) '
write (o_unit,10) (cellstate(i,jmax)%egt, i = 1, imax-1)
write (o_unit,10) (cellstate(i,0)%egt, i = 1, imax-1)
write (o_unit,10) (cellstate(imax,j)%egt, j = 1, jmax-1)
write (o_unit,10) (cellstate(0,j)%egt, j = 1, jmax-1)
write (o_unit,*)
write (o_unit,8)
write (o_unit,*) ' top(i=1,imax-1,j=jmax) ', &
'bottom(i=1,imax-1,j=0) ', &
'right(i=imax,j=1,jmax-1) ', &
'left(i=0, j=1,jmax-1) '
write (o_unit,10) (cellstate(i,jmax)%egtss, i = 1, imax-1)
write (o_unit,10) (cellstate(i,0)%egtss, i = 1, imax-1)
write (o_unit,10) (cellstate(imax,j)%egtss, j = 1, jmax-1)
write (o_unit,10) (cellstate(0,j)%egtss, j = 1, jmax-1)
write (o_unit,*)
write (o_unit,9)
write (o_unit,*) ' top(i=1,imax-1,j=jmax) ', &
'bottom(i=1,imax-1,j=0) ', &
'right(i=imax,j=1,jmax-1) ', &
'left(i=0,j=1,jmax-1) '
write (o_unit,11) (cellstate(i,jmax)%egt10, i = 1, imax-1)
write (o_unit,11) (cellstate(i,0)%egt10, i = 1, imax-1)
write (o_unit,11) (cellstate(imax,j)%egt10, j = 1, jmax-1)
write (o_unit,11) (cellstate(0,j)%egt10, j = 1, jmax-1)
write (o_unit,*)
write (o_unit,12)
write (o_unit,*) ' top(i=1,imax-1,j=jmax) ', &
'bottom(i=1,imax-1,j=0) ', &
'right(i=imax,j=1,jmax-1) ', &
'left(i=0,j=1,jmax-1) '
write (o_unit,11) (cellstate(i,jmax)%egt2_5, i = 1, imax-1)
write (o_unit,11) (cellstate(i,0)%egt2_5, i = 1, imax-1)
write (o_unit,11) (cellstate(imax,j)%egt2_5, j = 1, jmax-1)
write (o_unit,11) (cellstate(0,j)%egt2_5, j = 1, jmax-1)
write (o_unit,*)
write (o_unit,fmt="(' | ',3('|',a))") 'Total Soil Loss', 'soil loss', '(kg/m^2)'
do 19 j = jmax-1, 1, -1
write (o_unit,10) (cellstate(i,j)%egt, i = 1, imax-1)
19 continue
write (o_unit,fmt="(' ')")
write (o_unit,*)
write (o_unit,fmt="(' | ',3('|',a))") 'Saltation/Creep Soil Loss', 'salt/creep soil loss', '(kg/m^2)'
do 29 j = jmax-1, 1, -1
write (o_unit,10) (cellstate(i,j)%egt-cellstate(i,j)%egtss, i = 1, imax-1)
29 continue
write (o_unit,fmt="(' ')")
write (o_unit,*)
write (o_unit,fmt="(' | ',3('|',a))") 'Suspension Soil Loss', 'suspension soil loss', '(kg/m^2)'
do 39 j = jmax-1, 1, -1
write (o_unit,10) (cellstate(i,j)%egtss, i = 1, imax-1)
39 continue
write (o_unit,fmt="(' ')")
write (o_unit,*)
write (o_unit,fmt="(' | ',3('|',a))") 'PM10 Soil Loss', 'PM10 soil loss', '(kg/m^2)'
do 49 j = jmax-1, 1, -1
write (o_unit,11) (cellstate(i,j)%egt10, i = 1, imax-1)
49 continue
write (o_unit,fmt="(' ')")
write (o_unit,*)
write (o_unit,fmt="(' | ',3('|',a))") 'PM2_5 Soil Loss', 'PM2_5 soil loss', '(kg/m^2)'
do 59 j = jmax-1, 1, -1
write (o_unit,11) (cellstate(i,j)%egt2_5, i = 1, imax-1)
59 continue
write (o_unit,fmt="(' ')")
write (o_unit,*)
write (o_unit,*) '**Averages - Field'
write (o_unit,*) ' Total salt/creep susp PM10 PM2.5'
write (o_unit,*) ' egt egtss egt10 egt2_5'
write (o_unit,*) ' -----------------------kg/m^2---------------------------'
write (o_unit,15) aegt, aegt-aegtss, aegtss, aegt10, aegt2_5
write (o_unit,*)
write (o_unit,*) '**Averages - Crossing Boundaries '
write (o_unit,*) 'Location Total Salt/Creep Susp PM10 PM2_5'
write (o_unit,*) '-------------------------kg/m---------------------------'
write (o_unit,21) topt+topss, topt, topss, top10, top2_5
write (o_unit,22) bott+botss, bott, botss, bot10, bot2_5
write (o_unit,23) ritt+ritss, ritt, ritss, rit10, rit2_5
write (o_unit,24) lftt+lftss, lftt, lftss, lft10, lft2_5
write (o_unit,*)
write (o_unit,*) ' Comparison of interior & boundary loss'
write (o_unit,*) ' interior boundary int/bnd ratio'
if( totbnd.gt.1.0e-9 ) then
write (o_unit,16) tot, totbnd, tot/totbnd
else
!Boundary loss near or equal to zero
write (o_unit,16) tot, totbnd, 1.0e-9
end if
! additional output statements for easy shell script parsing
write (o_unit,*)
! write losses as positive numbers
write (o_unit,17) -aegt, aegtss-aegt, -aegtss, -aegt10, -aegt2_5
17 format (' repeat of total, salt/creep, susp, PM10:', 3f12.4,2f12.6)
close(o_unit)
! output formats
6 format (1x,' Passing Border Grid Cells - Total egt+egtss(kg/m)')
7 format (1x,' Passing Border Grid Cells - Salt/Creep egt(kg/m)')
8 format (1x,' Passing Border Grid Cells - Suspension egtss(kg/m)')
9 format (1x,' Passing Border Grid Cells - PM10 egt10(kg/m)')
12 format (1x,' Passing Border Grid Cells - PM2_5 egt2_5(kg/m)')
10 format (1x, 500f12.4)
11 format (1x, 500f12.6)
15 format (1x, 3(f12.4,2x), 2(f12.6,2x))
16 format (1x, 2(f13.4,2x),2x, f13.4)
21 format (1x, 'top ', 1x, 5(f9.2,1x))
22 format (1x, 'bottom', 1x, 5(f9.2,1x))
23 format (1x, 'right ', 1x, 5(f9.2,1x))
24 format (1x, 'left ', 1x, 5(f9.2,1x))
end if !if (btest(am0efl,1)) then
!Erosion summary - total, salt/creep, susp, pm10
!(loss values are positive - deposition values are negative)
if (btest(am0efl,0)) then
if( in_weps ) then
call caldat (mksaeout%jday,day,mon,yr)
write(*,'(4(a,i0))') 'Made Daily Erosion summary file for: ', day,'/', mon,'/', yr,' simulation day: ', mksaeout%simday
write (UNIT=o_E_unit,FMT="(5(f12.6),' ')",ADVANCE="NO") -aegt, -(aegt-aegtss), -aegtss, -aegt10, -aegt2_5
write (UNIT=o_E_unit,FMT="('# WEPS erosion day mon yr',2(1x,i2),2x,i4)",ADVANCE="NO") day, mon, yr
write (UNIT=o_E_unit,FMT="(A)",ADVANCE="YES") ' (loss values are positive - deposition values are negative)'
else
write (UNIT=o_E_unit,FMT="(5(f12.6),' ')",ADVANCE="NO") -aegt, -(aegt-aegtss), -aegtss, -aegt10, -aegt2_5
write (UNIT=o_E_unit,FMT="(A)",ADVANCE="NO") trim(input_filename)
write (UNIT=o_E_unit,FMT="(A)",ADVANCE="YES") ' (loss values are positive - deposition values are negative)'
end if
end if
!Duplicate Erosion summary info for the *.sgrd file so "sweep" interface
! can display this info on graphical report window
if (btest(am0efl,3) .and. (sgrd_u .ge. 0) ) then
write (sgrd_u,*)
write (sgrd_u,*) '**Averages - Field'
write (sgrd_u,*) ' Total salt/creep susp PM10 PM2.5'
write (sgrd_u,*) ' egt egtss egt10 egt2_5'
write (sgrd_u,*) ' -----------------------kg/m^2---------------------------'
write (sgrd_u,15) aegt, aegt-aegtss, aegtss, aegt10
write (sgrd_u,*)
write (sgrd_u,*) '**Averages - Crossing Boundaries '
write (sgrd_u,*) 'Location Total Salt/Creep Susp PM10 PM2_5'
write (sgrd_u,*) '--------------------------kg/m--------------------------'
write (sgrd_u,21) topt+topss, topt, topss, top10, top2_5
write (sgrd_u,22) bott+botss, bott, botss, bot10, bot2_5
write (sgrd_u,23) ritt+ritss, ritt, ritss, rit10, rit2_5
write (sgrd_u,24) lftt+lftss, lftt, lftss, lft10, lft2_5
write (sgrd_u,*)
write (sgrd_u,*) ' Comparison of interior & boundary loss'
write (sgrd_u,*) ' interior boundary int/bnd ratio'
if( totbnd.gt.1.0e-9 ) then
write (sgrd_u,16) tot, totbnd, tot/totbnd
else
!Boundary loss near or equal to zero
write (sgrd_u,16) tot, totbnd, 1.0e-9
end if
if( sgrd_u .ge. 0 ) then
close(sgrd_u)
end if
end if
end subroutine daily_erodout
subroutine sb1out( jj, nn, hr, ws, wdir, o_unit, subrsurf, cellstate )
! + + + PURPOSE + + +
! To print to file tst.out some key variables used in erosion
! use wind dir of 270 for most to see output along wind direction
use datetime_mod, only: get_systime_string, caldat
use erosion_data_struct_defs, only: subregionsurfacestate, cellsurfacestate, awzypt, anemht, wzoflg, ntstep
use grid_mod, only: awa, imax, jmax, amasim, amxsim
! + + + ARGUEMENT DECLARATIONS + + +
real ws, wdir, hr
integer jj, nn, o_unit
type(subregionsurfacestate), intent(in) :: subrsurf ! subregion surface conditions (erosion specific set)
type(cellsurfacestate), dimension(0:,0:), intent(inout) :: cellstate ! initialized grid cell state values
! + + + ARGUMENT DEFINITIONS + + +
! o_unit= Unit number for output file
! + + + LOCAL VARIABLES + + +
character(len=21) :: rundatetime
!integer m, n, k
integer initflag, ipd, npd
save initflag, ipd, npd
integer yr, mo, da
real hhrr, tims
save yr, mo, da, hhrr, tims
integer i,j, kbr
! + + + END SPECIFICATIONS + + +
! output headings?
if (initflag .eq. 0) then
ipd = 0
npd = nn * ntstep
tims = 3600*24/ntstep ! seconds in each emission period
call caldat( mksaeout%jday, da, mo, yr) ! Set day, month and year
hhrr = 0 - tims/3600 !Pre-set hhrr so we get end of period times
write (o_unit,*)
write (o_unit,*) 'OUT PUT from sb1out'
write (o_unit,*)
! Print date of Run
rundatetime = get_systime_string() ! with Lahey f95, had to assign to variable first
write(o_unit,"(1x,'Date of run: ',a21)") rundatetime
write(o_unit,*)
write (unit=o_unit,fmt="(a,f5.2,a2,a,i1)") ' anemht = ', anemht, 'm', ' wzoflg = ', wzoflg
write (unit=o_unit,fmt="(a,f6.2,a4)") ' wind direction = ', wdir, 'deg'
write (unit=o_unit,fmt="(a,f6.2,a4)") ' wind direction relative to field orientation = ', awa, 'deg'
write (o_unit,*)
! this section inserted for compatibility with previous output version
if (wdir .ge. 337.5 .or. wdir .lt. 22.5) then
kbr = 1
elseif (wdir .lt. 67.5) then
kbr = 2
elseif (wdir.lt. 112.5) then
kbr = 3
elseif (wdir .lt. 157.5) then
kbr = 4
elseif (wdir .lt. 202.5) then
kbr = 5
elseif (wdir.lt. 247.5) then
kbr = 6
elseif (wdir .lt. 292.5) then
kbr = 7
else
kbr = 8
endif
write (unit=o_unit,fmt="(a,i1)") ' wind quadrant = ', kbr
write (o_unit,*)
! end inserted section
write (o_unit,*) 'orientation and dimensions of sim region'
write (o_unit,*) 'amasim(deg) amxsim - (x1,y1) (x2,y2)'
write(o_unit,fmt="(1x,5f8.2)") amasim, amxsim(1)%x, amxsim(1)%y, amxsim(2)%x, amxsim(2)%y
write (o_unit,*)
write (o_unit,*) "Surface properties"
write (o_unit,fmt="(a,f8.2,a)") "Ridge spacing parallel to wind direction", subrsurf%sxprg, " (mm)"
write (o_unit,fmt="(a,f5.2,a)") "Crop row spacing", subrsurf%acxrow, " (mm)"
write (o_unit,fmt="(a,i2,a)") "Crop seeding location relative to ridge", subrsurf%ac0rg, " (0 - furrow, 1 - ridge)"
write (o_unit,fmt="(a,f5.2,a)") "Composite weighted average biomass height", subrsurf%abzht, " (m)"
write (o_unit,fmt="(a,f5.2,a)") "Biomass leaf area index", subrsurf%abrlai, " (m^2/m^2)"
write (o_unit,fmt="(a,f5.2,a)") "Biomass stem area index", subrsurf%abrsai, " (m^2/m^2)"
write (o_unit,fmt="(a,f5.2,a)") "Biomass flat cover", subrsurf%abffcv, " (m^2/m^2)"
write (o_unit,fmt="(a,f8.2,a)") "Average yearly total precipitation ", awzypt, " (mm)"
write (o_unit,*)
write(o_unit,fmt="(1x,a)") ""
write(o_unit,fmt="(1x,5f10.2)")amasim, amxsim(1)%x, amxsim(1)%y, amxsim(2)%x, amxsim(2)%y
write(o_unit,fmt="(1x,a)") ""
write (o_unit,*)
initflag = 1 ! turn off heading output
endif
ipd = ipd + 1
if (hhrr .ge. 24) then
hhrr = tims/3600
call caldat( mksaeout%jday, da, mo, yr) ! Set day, month and year
else
hhrr = hhrr + tims/3600
endif
call caldat( mksaeout%jday, da, mo, yr) ! Set day, month and year
! write (o_unit, fmt="(a, 3(i3), f6.2, 4(i4)f7.2)") ' day mon yr hhrr upd_pd jj nn npd ',da,mo,yr,hr,ipd,jj,nn,npd,hhrr
write (o_unit, fmt="(a, i5, 2(i3), f7.3, 4(i4))") ' yr mon day hr upd_pd jj nn(subpd) npd (sbqout 1)', &
yr,mo, da, hr,ipd, jj,nn, npd
write (o_unit,*)
write (o_unit, fmt="(a, f5.2, 2(f7.2))") ' pd wind speed, dir and dir rel to field ', ws, wdir, awa
write (o_unit,*)
write (o_unit,*) "Surface layer properties"
write (o_unit,fmt="(a,f5.2,a)") "Surface course fragments", subrsurf%bsl(1)%asvroc, " (m^3/m^3)"
write (o_unit,fmt="(a,a,f5.2,a)") "Initial soil ", "mass fraction in surface layer < 0.10 mm ", subrsurf%sf10ic, " (kg/kg)"
write (o_unit,fmt="(a,a,f5.2,a)") "Initial soil ", "mass fraction in surface layer < 0.84 mm ", subrsurf%sf84ic, " (kg/kg)"
write (o_unit,*) "PM10 emission properties"
write (o_unit,fmt="(a,f5.2,a)") "Soil fraction PM10 in abraded suspension ", subrsurf%asf10an
write (o_unit,fmt="(a,f5.2,a)") "Soil fraction PM10 in emitted suspension ", subrsurf%asf10en
write (o_unit,fmt="(a,f5.2,a)") "Soil fraction PM10 in saltation breakage suspension ", subrsurf%asf10bk
write (o_unit,fmt="(a,f5.2,a)") "Coefficient of abrasion of aggregates ", subrsurf%acanag
write (o_unit,fmt="(a,f5.2,a)") "Coefficient of abrasion of crust ", subrsurf%acancr
!Grid cell data
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Surface Friction Velocity', 'friction velocity', '(m/s)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%wus, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Threshold Surface Friction Velocity', 'threshold friction velocity', '(m/s)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%wust, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Transport Threshold Surface Friction Velocity', 'transport threshold friction velocity', '(m/s)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%wusp, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write (o_unit,*)
!Grid Cell Surface properties
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Surface Random Roughness', 'random roughness', '(mm)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%slrr, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Surface Oriented Roughness', 'ridge height', '(mm)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%szrgh, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Surface Rock', 'surface volume rock fraction', '(m^3/m^3)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%svroc, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write (o_unit,*)
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size<0.01', 'mass fraction < 0.01 mm size', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf1, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size<0.1', 'mass fraction < 0.1 mm size', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf10, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size<0.84', 'mass fraction < 0.84 mm size', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf84, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size<2.0', 'mass fraction < 2.0 mm size', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf200, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size for u* to be the thresh. friction velocity', '"effective" mass fraction < 0.84 mm size', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf84mn, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Mobile soil removable from aggregated surface', 'mass removable', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%smaglos, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Change in mobile soil on aggregated surface', 'net mass change', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%dmlos, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
! Crust properties
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Consolidated crust thickness', 'crust thickness', '(mm)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%szcr, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Fraction of Surface covered with Crust','crust cover','(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sfcr, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Fraction of Crusted Surface covered with Loose Erodible Soil ', 'loose erodible material', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sflos, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Mass of Loose Erodible Soil on Crusted Surface', 'loose erodible material', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%smlos, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write (o_unit,*)
end subroutine sb1out
subroutine sb2out (jj, nn, hr, o_unit, cellstate)
! + + + PURPOSE + + +
! To print to file tst.out some key variables used in erosion
! use wind direction of 270 to see output along downwind direction
use datetime_mod, only: caldat
use erosion_data_struct_defs, only: cellsurfacestate, ntstep
use grid_mod, only: imax, jmax
! + + + ARGUEMENT DECLARATIONS + + +
real hr
integer jj, nn, o_unit
type(cellsurfacestate), dimension(0:,0:), intent(inout) :: cellstate ! initialized grid cell state values
! + + + ARGUMENT DEFINITIONS + + +
! o_unit= Unit number for output file
! + + + LOCAL VARIABLES + + +
real egavg(imax)
integer m, n, k, icsr
integer initflag, ipd, npd
save initflag, ipd, npd
integer yr, mo, da
real hhrr, tims
save yr, mo, da, hhrr, tims
integer i,j
! outflag = 0 - print heading output, 1 - no more heading
! + + + END SPECIFICATIONS + + +
! define index of current subregions
icsr = 1
! output headings?
if (initflag .eq. 0) then
ipd = 0
npd = nn * ntstep
tims = 3600*24/ntstep !seconds in each emission period
call caldat( mksaeout%jday, da, mo, yr) ! Set day, month and year
hhrr = 0 !Pre-set hhrr so we get start of period times
! write (o_unit,*)
! write (o_unit,*) 'OUT PUT from sb2out'
! write (o_unit,*)
initflag = 1 ! turn off heading output
endif
ipd = ipd + 1
if (hhrr .ge. 24) then
hhrr = tims/3600
call caldat( mksaeout%jday, da, mo, yr) ! Set day, month and year
else
hhrr = hhrr + tims/3600
endif
call caldat( mksaeout%jday, da, mo, yr) ! Set day, month and year
! write (o_unit, fmt="(a, 3(i3), f5.2, i3)") &
! ' day mon yr hhrr update_period ', da, mo, yr, hhrr, jj
write (o_unit, fmt="(a, i5, 2(i3), f7.3, 4(i4))") &
' yr mon day hr upd_pd jj nn(subpd) npd (sbqout 2)', &
yr,mo, da, hr,ipd, jj,nn, npd
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Cumulative Total Soil Loss', 'soil loss', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%egt, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Cumulative Saltation/Creep Soil Loss', 'salt/creep soil loss', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)")(cellstate(i,j)%egt-cellstate(i,j)%egtss,i=1,imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Cumulative Suspension Soil Loss', 'suspension soil loss', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%egtss, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Cumulative PM10 Soil Loss', 'PM10 soil loss', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.6)") (cellstate(i,j)%egt10, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Cumulative PM2_5 Soil Loss', 'PM2_5 soil loss', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.6)") (cellstate(i,j)%egt2_5, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
!! if (ipd .eq. npd) then
!Grid Cell Surface properties
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Surface Random Roughness (after)', 'random roughness', '(mm)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%slrr, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Surface Oriented Roughness (after)', 'ridge height', '(mm)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%szrgh, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Surface Rock (after)', 'surface volume rock fraction', '(m^3/m^3)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%svroc, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write (o_unit,*)
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size < 0.01','mass fraction < 0.01 mm size','(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf1, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size < 0.1', 'mass fraction < 0.1 mm size', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf10, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size < 0.84','mass fraction < 0.84 mm size','(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf84, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size < 2.0', 'mass fraction < 2.0 mm size', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf200, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Soil Agg. Size for u* to be the thresh. friction velocity (af)','"effective" mass fraction < 0.84 mm size', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sf84mn, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Mobile soil removable from aggregated surface (after)', 'mass removable', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%smaglos, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Change in mobile soil on aggregated surface (after)', 'net mass change', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%dmlos, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
! Crust properties
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Consolidated crust thickness (after)', 'crust thickness', '(mm)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%szcr, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Fraction of Surface covered with Crust (after)', 'crust cover','(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sfcr, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Fraction of Crusted Surface covered with Loose Erodible Soil(a)', 'loose erodible material', '(fract.)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%sflos, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write(o_unit,fmt="(' |',i5,2(i3),f7.3,3('|',a))") yr, mo, da, hr, &
'Mass of Loose Erodible Soil on Crusted Surface (after)', 'loose erodible material', '(kg/m^2)'
do j = jmax-1, 1, -1
write (o_unit, fmt="(500f12.4)") (cellstate(i,j)%smlos, i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write (o_unit,*)
!! endif
! set output increment
m = (imax - 1)/8
m = max0(m,1)
n = (jmax-1)/2
n = max(n,1)
! initialize avg erosion variable
do 3 j = 1, imax
egavg(j) = 0.0
3 continue
! calc. avg erosion over a given field length
do 5 i = 1,(imax-1)
!average over y-direction
do 4 j = 1, (jmax-1)
egavg(i) = egavg(i) + cellstate(i,j)%egt/(jmax-1)
4 continue
5 continue
!average over x-direction
do 6 i = 2, (imax-1)
egavg(i) = ((i-1)*(egavg(i-1))+egavg(i))/i
6 continue
write (o_unit,35) (egavg(k), k=1,(imax-1))
write (o_unit,*) '----------------------------------------------'
! output formats
35 format (1x, 'egavg = ', 20f6.2)
end subroutine sb2out
subroutine sbemit (ounit, ws, hhr, cellstate, first_emit)
! To calc the emissions for each time step of the input wind speed
! The emissions for EPA are the suspension component
! with units kg m-2 s-1.
! To write out a file in the format:
! 12 blank col, yr, mo, day, hr, soucename, emissionrate
!
! Instructions & logic:
! To get ntstep period emissions output on erosion days:
! user sets am0efl = 3 in WEPS configuration screen
! subroutine openfils creates output file emit.out
! EROSION calls sbemit to write heading in emit.out file,
! & sets am0efl to 98, then calls sbemit
! to print (hourly) Weps emissions on erosion days.
! or
! user sets ae0efl (print flg)=4 in stand_alone input file
! EROSION opens emit.out file, calls sbemit to write headings
! & sets ae0efl to 99, then calls sbemit
! to print period emissions for an erosion day.
use datetime_mod, only: get_systime_string, caldat
use erosion_data_struct_defs, only: cellsurfacestate, ntstep
use grid_mod, only: imax, jmax
! +++ ARGUMENT DECLARATIONS +++
integer ounit !Unit number for detail grid erosion
real ws, hhr
type(cellsurfacestate), dimension(0:,0:), intent(inout) :: cellstate ! initialized grid cell state values
logical, intent(inout) :: first_emit ! indicates entry of emit from erosion for the first time this day
! +++ LOCAL VARIABLES +++
character(len=21) :: rundatetime
integer initflg
save initflg
integer j,i
integer yr, mo, da
save yr, mo, da
real tims, aegtp, aegtssp, aegt10p, aegt2_5p
save tims, aegtp, aegtssp, aegt10p, aegt2_5p
! real hr
! save hr
real aegt, aegtss, aegt10, aegt2_5 ! these have local scope
real emittot, emitss, emit10, emit2_5, tt
! +++ OUTPUT FORMATS +++
100 format (1x,' yr mo day hr ws emission (kg m-2 s-1)')
110 format (22x,' total salt/creep susp PM10 PM2_5')
120 format (1x,3(i4),F7.3,F6.2, 1x,5(F11.8))
! +++ END SPECIFICATIONS +++
! set initial conditions
if (initflg .eq. 0) then
initflg = initflg + 1
tims = 3600*24/ntstep !seconds in each emission period
call caldat( mksaeout%jday, da, mo, yr) ! Set day, month and year
write(0,*) 'First ntstep is: ', ntstep, tims, tims/3600
write (ounit,*) 'SBEMIT output'
! write (ounit,*) 'Suspended emissions < 0.10 mm dia.'
write (ounit,*)
! Print date of Run
rundatetime = get_systime_string() ! with Lahey f95, had to assign to variable first
write(ounit,"(1x,'Date of run: ',a21)") rundatetime
write (ounit,*)
write (ounit,100)
write (ounit,110)
write (ounit,*)
endif
! init prev erosion hr values to zero if this is new erosion day
if( first_emit ) then
first_emit = .false.
aegtp = 0.0
aegtssp = 0.0
aegt10p = 0.0
aegt2_5p = 0.0
endif
!write(0,*) 'Subsequent ntstep is: ', ntstep, tims, tims/3600
call caldat( mksaeout%jday, da, mo, yr) ! Set day, month and year
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
aegt2_5 = 0.0
do j=1,jmax-1
do i= 1, imax-1
aegt= aegt + cellstate(i,j)%egt
aegtss = aegtss + cellstate(i,j)%egtss
aegt10 = aegt10 + cellstate(i,j)%egt10
aegt2_5 = aegt2_5 + cellstate(i,j)%egt2_5
enddo
enddo
tt = (imax-1)*(jmax-1)
aegt = - aegt/tt ! change signs to positive=emission
aegtss = - aegtss/tt
aegt10 = - aegt10/tt
aegt2_5 = - aegt2_5/tt
emittot = (aegt - aegtp)/tims
emitss = (aegtss - aegtssp)/tims
emit10 = (aegt10 - aegt10p)/tims
emit2_5 = (aegt2_5 - aegt2_5p)/tims
! Save prior hour average emission
aegtp = aegt
aegtssp = aegtss
aegt10p = aegt10
aegt2_5p = aegt2_5
! Write to emit.out file
write (ounit,120) yr, mo, da, hhr, ws, emittot, emittot-emitss, emitss, emit10, emit2_5
end subroutine sbemit
end module sae_in_out_mod