!**********************************************************************
! subroutine sb1out
!**********************************************************************
subroutine sb1out (jj, nn, hr, ws, wdir, o_unit)
implicit none
!
! + + + 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
! + + + ARGUEMENT DECLARATIONS + + +
real ws, wdir, hr
integer jj, nn, o_unit
!
! + + + ARGUMENT DEFINITIONS + + +
! anemHght =
! aeroRougStatHght =
! wzz0 =
! awu =
! wus =
! thrsFricVelc =
! o_unit= Unit number for output file
!
! + + + GLOBAL COMMON BLOCKS + + +
include 'p1werm.inc'
include 'h1db1.inc'
include 'b1glob.inc'
include 'c1gen.inc'
include 's1surf.inc'
include 'w1clig.inc'
!
! + + + LOCAL COMMON BLOCKS + + +
!
include 'p1const.inc'
include 'm1sim.inc'
include 's1dbh.inc'
include 'm1geo.inc'
include 'erosion/s2agg.inc'
include 'erosion/s2surf.inc'
include 'erosion/s2sgeo.inc'
include 'erosion/w2wind.inc'
include 'erosion/m2geo.inc'
include 'erosion/e2erod.inc'
include 'erosion/e3grid.inc'
!
! + + + LOCAL VARIABLES + + +
!
integer m, n, k, icsr, x, y
integer initflag, ipd, npd
save initflag, ipd, npd
integer yr, mo, da
real hhrr, tims
save yr, mo, da, hhrr, tims
integer i,j
integer :: dt(8)
character(len=3) :: mstring
common / datetime / dt, mstring
!
! 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 caldatw (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
212 format(1x,'Date of run: ',a3,' ',i2.2,', ',i4,' ', &
& i2.2,':',i2.2,':',i2.2)
write(o_unit,212) mstring, dt(3), dt(1), dt(5), dt(6), dt(7)
write(o_unit,*)
write (unit=o_unit,fmt="(a,f5.2,a2,a,i1)") &
& ' anemHght = ', anemHght, 'm', ' aeroAnemFlg = ', aeroAnemFlg
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,*)
write (unit=o_unit,fmt="(a,i1)") ' wind quadrant = ', kbr
write (o_unit,*)
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(x,y),x=1,2),y=1,2)
write (o_unit,*)
write (o_unit,*) "Surface properties"
write (o_unit,fmt="(a,f8.2,a)") &
& "Ridge spacing parallel to wind direction", ridgSpacParaWind(icsr), " (mm)"
write (o_unit,fmt="(a,f5.2,a)") &
& "Crop row spacing", cropRowSpac(icsr), " (mm)"
write (o_unit,fmt="(a,i2,a)") &
& "Crop seeding location relative to ridge", cropFurrFlg(icsr), &
& " (0 - furrow, 1 - ridge)"
write (o_unit,fmt="(a,f5.2,a)") &
& "Composite weighted average biomass height", abioMassHght(icsr), " (m)"
write (o_unit,fmt="(a,f5.2,a)") &
& "Biomass leaf area index", abrlai(icsr), " (m^2/m^2)"
write (o_unit,fmt="(a,f5.2,a)") &
& "Biomass stem area index", abrsai(icsr), " (m^2/m^2)"
write (o_unit,fmt="(a,f5.2,a)") &
& "Biomass flat cover", abioFracFlatCovr(icsr), " (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(x,y),x=1,2),y=1,2)
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 caldatw (da, mo, yr)
else
hhrr = hhrr + tims/3600
endif
call caldatw (da, mo, yr)
! 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", asoilLayrRock(1,1), " (m^3/m^3)"
write (o_unit,fmt="(a,a,f5.2,a)") "Initial soil ", &
& "mass fraction in surface layer < 0.10 mm ", soilFracDiaLt10ic, " (kg/kg)"
write (o_unit,fmt="(a,a,f5.2,a)") "Initial soil ", &
& "mass fraction in surface layer < 0.84 mm ", soilFracDiamLt84ic, " (kg/kg)"
write (o_unit,*) "PM10 emission properties"
write (o_unit,fmt="(a,f5.2,a)") &
& "Soil fraction PM10 in abraded suspension ", asoilPM10FracAbraSusp(1)
write (o_unit,fmt="(a,f5.2,a)") &
& "Soil fraction PM10 in emitted suspension ", asoilPM10FracEmitSusp(1)
write (o_unit,fmt="(a,f5.2,a)") &
& "Soil fraction PM10 in saltation breakage suspension ",asoilPM10FracSaltBrkSusp(1)
write (o_unit,fmt="(a,f5.2,a)") &
& "Coefficient of abrasion of aggregates ", acoefAbraAgg(1)
write (o_unit,fmt="(a,f5.2,a)") &
& "Coefficient of abrasion of crust ", acoefAbraCrst(1)
!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)") (fricVelc(i,j), 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)") (thrsFricVelc(i,j), 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)") (thrsFricVelcTrap(i,j), 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)") (randRoug(i,j), 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)") (ridgHght(i,j), 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)") (soilLayrRock(i,j), 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)") (soilFracDiaLt1(i,j), 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)") (soilFracDiaLt10(i,j), 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)") (soilFracDiamLt84(i,j), 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)") (soilFracDiaLt200(i,j), 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)") (soilFracDiamLt84mn(i,j), 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)") (soilMassAvalAggSurf(i,j), 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)") (soilMassAvalDelt(i,j), 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)") (soilCrstThck(i,j), 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)") (soilCrstFrac(i,j), 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)") (soilLoosCovFrac(i,j), 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)") (soilLoosMass(i,j), i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write (o_unit,*)
! write (o_unit,20) anemHght,aeroAnemFlg,kbr, jj, ws
! set output increment
! m = (imax - 1)/8
! m = max0(m,1)
! n = (jmax-1)/2
! n = max(n,1)
!
! write (o_unit,*) 'sb1out output'
! write (o_unit,*) 'for prior wind speed'
! write (o_unit,21) (egt(k,n),k=1,(imax-1),m)
! write (o_unit,22) (egtss(k,n),k=1,(imax-1),m)
! write (o_unit,23??) ((egtss(k,n)/(egt(k,n)+0.0001)),k=1,(imax-1),m)
! write (o_unit,*)
! write (o_unit,18) (k , k=1,(imax-1),m), n
! write (o_unit,13) (soilFracDiaLt1(k,n),k=1,(imax-1),m)
! write (o_unit,23) (soilFracDiaLt10(k,n),k=1,(imax-1),m)
! write (o_unit,24) (soilFracDiamLt84(k,n),k=1,(imax-1),m)
! write (o_unit,35) (soilFracDiaLt200(k,n),k=1,(imax-1),m)
! write (o_unit,12) (soilLayrRock(k,n),k=1,(imax-1),m)! edit ljh 1-22-05
! write (o_unit,36) (soilMassAvalDelt(k,n),k=1,(imax-1),m)
! write (o_unit,37) (soilMassAvalAggSurf(k,n),k=1,(imax-1),m)
! write (o_unit,43) (soilMassAvalAggSurfmx(k,n),k=1,(imax-1),m)
! write (o_unit,39) (soilFracDiamLt84mn(k,n),k=1,(imax-1),m)
! write (o_unit,40) soilFracDiamLt84ic, soilFracDiaLt10ic, asoilLayrRock(1,1) !edit ljh 1-22-05
! write (o_unit,42) acoefAbraAgg(1), acoefAbraCrst(1), awzypt
! write (o_unit,10) asoilPM10FracAbraSusp(1), asoilPM10FracEmitSusp(1), asoilPM10FracSaltBrkSusp(1)
! write (o_unit,25) (soilCrstThck(k,n),k=1,(imax-1),m)
! write (o_unit,26) (soilCrstFrac(k,n),k=1,(imax-1),m)
! write (o_unit,27) (soilLoosMass(k,n),k=1,(imax-1),m)
! write (o_unit,28) (soilLoosCovFrac(k,n),k=1,(imax-1),m)
! write (o_unit,29) (ridgHght(k,n),k=1,(imax-1),m)
! write (o_unit,30) (randRoug(k,n),k=1,(imax-1),m)
! write (o_unit,38) ridgSpacParaWind(icsr), abioMassHght(icsr), abrlai(icsr), &
! & abrsai(icsr), abioFracFlatCovr(icsr)
! write (o_unit,41) cropRowSpac(icsr), cropFurrFlg(icsr)
! write (o_unit,31) asoilLayrWiltPt(1,1), asoilSurfWatrCont(12,1)
! write (o_unit,32) (fricVelc(k,n),k=1,(imax-1),m)
! write (o_unit,33) (thrsFricVelcTrap(k,n),k=1,(imax-1),m)
! write (o_unit,34) (thrsFricVelc(k,n),k=1,(imax-1),m)
! write (o_unit,44) fricVelcMod
! write (o_unit,*)
! output formats
10 format (1x, 'soilPM10FracAbraSusp =',f6.3,' soilPM10FracEmitSusp =',f6.3,' soilPM10FracSaltBrkSusp =',f6.3)
! 15 format (1x, ' (m) (m/s) ')
18 format (1x, 'i..n,j', 3i6, 17i7)
20 format (1x, 'anemHght aeroAnemFlg kbr jj ws', &
& f6.0, 3i6, f6.2)
21 format (1x, 'egt=', 20f6.2)
22 format (1x, 'egtss=', 20f6.2)
13 format (1x, 'soilFracDiaLt1= ', 20f7.4)
23 format (1x, 'soilFracDiaLt10= ', 20f7.3)
24 format (1x, 'soilFracDiamLt84= ', 20f7.3)
12 format (1x, 'soilLayrRock=', 20f7.3) !edit ljh 1-22-05
35 format (1x, 'soilFracDiaLt200=', 20f7.3)
36 format (1x, 'soilMassAvalDelt=', 20f7.3)
37 format (1x, 'soilMassAvalAggSurf=',20f7.3)
43 format (1x, 'soilMassAvalAggSurfmx=',20f7.3)
39 format (1x, 'soilFracDiamLt84mn=',20f7.3)
40 format (1x, 'soilFracDiamLt84ic =',f4.2,' soilFracDiaLt10ic =',f4.2,' asoilLayrRock=',f4.2)
42 format (1x, 'coefAbraAgg =', f6.3,' coefAbraCrst = 'f6.3,' awzypt=',f6.0)
41 format (1x, 'cropRowSpac=', f6.2, ' cropFurrFlg=', i3)
25 format (1x, 'soilCrstThck= ', 20f7.2)
26 format (1x, 'soilCrstFrac= ', 20f7.3)
27 format (1x, 'soilLoosMass=', 20f7.3)
28 format (1x, 'soilLoosCovFrac=', 20f7.3)
29 format (1x, 'ridgHght=', 20f7.2)
30 format (1x, 'randRoug= ', 20f7.2)
38 format (1x, 'ridgSpacParaWind=', f6.0, ' abioMassHght=', f6.2, ' abrlai=', f4.2, &
& ' abrsai=', f5.3, ' abioFracFlatCovr=',f4.3)
31 format (1x, 'asoilLayrWiltPt=',f4.2,' asoilSurfWatrCont(icsr,12)=', f6.2)
32 format (1x, 'fricVelc= ', 20f7.3)
33 format (1x, 'thrsFricVelcTrap=', 20f7.3)
34 format (1x, 'thrsFricVelc=', 20f7.3)
44 format (1x, 'fricVelcMod=', f5.3)
return
end
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!**********************************************************************
! subroutine sb2out
!**********************************************************************
subroutine sb2out (jj, nn, hr, ws, wdir, o_unit)
!
! + + + PURPOSE + + +
! To print to file tst.out some key variables used in erosion
! use wind direction of 270 to see output along downwind direction
!
! + + + ARGUEMENT DECLARATIONS + + +
real ws, wdir, hr
integer jj, nn, o_unit
!
! + + + ARGUMENT DEFINITIONS + + +
! anemHght = anemometer height (m)
! aeroRougStatHght = aerodynamic roughness at anemometer (mm)
! wzz0 = aerodynamic roughness length (mm)
! avgWindSped = wind speed (m/s)
! fricVelc = friction velocity (m/s)
! thrsFricVelc = threshold friction velocity (m/s)
! aeroAnemFlg = flag to showing anemometer at field (1) or wx sta (0)
! o_unit= Unit number for output file
!
! + + + GLOBAL COMMON BLOCKS + + +
include 'p1werm.inc'
include 'h1db1.inc'
include 'p1const.inc'
include 'm1sim.inc'
!
! + + + LOCAL COMMON BLOCKS + + +
!
include 'erosion/s2agg.inc'
include 'erosion/s2surf.inc'
include 'erosion/s2sgeo.inc'
include 'erosion/w2wind.inc'
include 'erosion/m2geo.inc'
include 'erosion/e2erod.inc'
include 'erosion/e3grid.inc'
!
! + + + LOCAL VARIABLES + + +
real egavg(maxGridPtsDir)
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 caldatw (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 caldatw (da, mo, yr)
else
hhrr = hhrr + tims/3600
endif
call caldatw (da, mo, yr)
! 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)") (egt(i,j), 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)")(egt(i,j)-egtss(i,j),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)") (egtss(i,j), 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)") (egt10(i,j), 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)") (randRoug(i,j), 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)") (ridgHght(i,j), 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)") (soilLayrRock(i,j), 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)") (soilFracDiaLt1(i,j), 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)") (soilFracDiaLt10(i,j), 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)") (soilFracDiamLt84(i,j), 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)") (soilFracDiaLt200(i,j), 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)") (soilFracDiamLt84mn(i,j), 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)") (soilMassAvalAggSurf(i,j), 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)") (soilMassAvalDelt(i,j), 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)") (soilCrstThck(i,j), 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)") (soilCrstFrac(i,j), 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)") (soilLoosCovFrac(i,j), 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)") (soilLoosMass(i,j), i = 1, imax-1)
end do
write(o_unit,fmt="(' ')")
write (o_unit,*)
!! endif
! write (o_unit,20) anemHght,aeroAnemFlg,ws,kbr
!
! 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 j = 1, (imax-1), m
! do 4 k = 1, j
! egavg(j) = egavg(j) + egt(k,n)
! 4 continue
! egavg(j) = egavg(j)/j
! 5 continue
! changed 1-12-07 LH
do 5 i = 1,(imax-1)
!average over y-direction
do 4 j = 1, (jmax-1)
egavg(i) = egavg(i) + egt(i,j)/(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,*) 'sb2out output'
! write (o_unit,18) (k , k=1,(imax-1),m), n
! write (o_unit,21) (egt(k,n),k=1,(imax-1),m)
! write (o_unit,22) (egtss(k,n),k=1,(imax-1),m)
! write (o_unit,23) ((egtss(k,n)/(egt(k,n)+0.00001)), &
! & k=1,(imax-1),m)
! changed 3-12-07 LH
write (o_unit,35) (egavg(k), k=1,(imax-1))
! write (o_unit,36)
! write (o_unit,37) (soilMassAvalAggSurf(k,n),k=1,(imax-1),m)
write (o_unit,*) '----------------------------------------------'
! output formats
! 10 format (1x, 'anemHght aeroAnemFlg ws kbr')
18 format(1x, 'i..n,j', 20i6)
! 15 format (1x, ' (m) (m/s) ')
20 format (1x, 'anemHght aeroAnemFlg ws kbr', &
& f6.0, i6, f8.2, i6)
21 format (1x, 'egt= ', 20f6.2)
22 format (1x, 'egtss= ', 20f6.2)
23 format (1x, 'egtss/egt=', 20f6.2)
35 format (1x, 'egavg = ', 20f6.2)
36 format (1x, 'corrected for Sf12')
37 format (1x, 'soilMassAvalAggSurf=', 20f6.2)
24 format (1x, 'soilFracDiamLt84=', 20f6.2)
25 format (1x, 'soilCrstThck=', 20f6.2)
26 format (1x, 'soilCrstFrac=', 20f6.2)
27 format (1x, 'soilLoosMass=', 20f6.2)
28 format (1x, 'soilLoosCovFrac=', 20f6.2)
29 format (1x, 'ridgHght=', 20f6.2)
30 format (1x, 'randRoug=', 20f6.2)
31 format (1x, 'asoilSurfWatrCont(icsr,12)', f6.2)
32 format (1x, 'fricVelc=', 20f6.2)
33 format (1x, 'thrsFricVelcTrap=', 20f6.2)
34 format (1x, 'thrsFricVelc=', 20f6.2)
return
end
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!$Author: joelevin $
!$Date: 2011-03-24 11:33:26 -0500 (Thu, 24 Mar 2011) $
!$Revision: 11724 $
!$HeadURL: https://svn.weru.ksu.edu/weru/weps1/trunk/weps.src/src/lib_erosion/daily_erodout.for $
subroutine daily_erodout (o_unit, o_E_unit)
! +++ PURPOSE +++
! To print output desired from standalone EROSION submodel
! +++ ARGUMENT DECLARATIONS +++
integer o_unit, o_E_unit
integer i, j
real aegt, aegtss, aegt10
real tt, lx, ly
real topt,topss, top10, bott, botss, bot10, ritt, ritss, rit10
real lftt, lftss, lft10, tot, totbnd
integer initflag
save initflag
integer yr, mo, da
! + + + GLOBAL COMMON BLOCKS + + +
include 'p1werm.inc'
include 'm1geo.inc'
include 'm1flag.inc'
! + + + LOCAL COMMON BLOCKS
include 'erosion/e2erod.inc'
include 'erosion/m2geo.inc'
!
integer x, y
integer :: dt(8)
character(len=3) :: mstring
common / datetime / dt, mstring
! ++++ ARGUMENT DEFINITIONS +++
! +++ SUBROUTINES CALLED+++
! plotout.for
! ++++ LOCAL VARIABLES +++
! +++ END SPECIFICATIONS +++
! Calculate Averages Crossing Borders
! top border
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
j = jmax
do 1 i = 1, imax-1
aegt = aegt + egt(i,j)
aegtss = aegtss + egtss(i,j)
aegt10 = aegt10 + egt10(i,j)
1 continue
! calc. average at top border
topt = aegt/(imax-1)
topss = aegtss/(imax-1)
top10 = aegt10/(imax-1)
! bottom border
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
j = 0
do 2 i = 1, imax-1
aegt = aegt + egt(i,j)
aegtss = aegtss + egtss(i,j)
aegt10 = aegt10 + egt10(i,j)
2 continue
! calc. average at bottom border
bott = aegt/(imax-1)
botss = aegtss/(imax-1)
bot10 = aegt10/(imax-1)
! right border
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
i = imax
do 3 j = 1, jmax-1
aegt = aegt + egt(i,j)
aegtss = aegtss + egtss(i,j)
aegt10 = aegt10 + egt10(i,j)
3 continue
! calc. average at right border
ritt = aegt/(jmax-1)
ritss = aegtss/(jmax-1)
rit10 = aegt10/(jmax-1)
!
! left border
aegt = 0.0
aegtss = 0.0
aegt10 = 0.0
i = 0
do 4 j = 1, jmax-1
aegt = aegt + egt(i,j)
aegtss = aegtss + egtss(i,j)
aegt10 = aegt10 + egt10(i,j)
4 continue
! calc. average at left border
lftt = aegt/(jmax-1)
lftss = aegtss/(jmax-1)
lft10 = aegt10/(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 + egt(i,j)
aegtss = aegtss + egtss(i,j)
aegt10 = aegt10 + egt10(i,j)
5 continue
tt = (imax-1)*(jmax-1)
aegt = aegt/tt
aegtss = aegtss/tt
aegt10 = aegt10/tt
! calculate comparision of boundary and interior losses
lx = amxsim(1,2) - amxsim(1,1)
ly = amxsim(2,2) - amxsim(2,1)
tot = aegt*lx*ly
totbnd = (topt + bott + topss + botss)*lx + &
& (ritt + lftt + ritss + lftss)*ly
if (btest(am0efl,1)) then
if (initflag .eq. 0) then !write header to files
write (o_unit,*) 'Grid cell output from daily_erodout.for'
write (o_unit,*)
! Print date of Run
12 format(1x,'Date of run: ',a3,' ',i2.2,', ',i4,' ', &
& i2.2,':',i2.2,':',i2.2)
write (o_unit,12) mstring, dt(3), dt(1), dt(5), dt(6), dt(7)
write (o_unit,*)
write (o_unit,fmt="(1x,a,5f10.2)") "", &
& amasim,((amxsim(x,y),x=1,2),y=1,2)
write (o_unit,*)
write (o_unit,*) &
& 'Total grid size: (', imax+1,',', jmax+1, ') ',&
& 'Inner grid size: (', imax-1,',', jmax-1, ')'
write (o_unit,*)
initflag = initflag + 1
endif
call caldatw (da, mo, yr) !get day, month and year
write (UNIT=o_unit,FMT="(' ',i5,i3,i3,' ')",ADVANCE="YES") &
& yr, mo, da
write (o_unit,*)
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) (egt(i,jmax)+egtss(i,jmax), i = 1, imax-1)
write (o_unit,10) (egt(i,0)+egtss(i,0), i = 1, imax-1)
write (o_unit,10) (egt(imax,j)+egtss(imax,j), j = 1, jmax-1)
write (o_unit,10) (egt(0,j)+egtss(0,j), 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) (egt(i,jmax), i = 1, imax-1)
write (o_unit,10) (egt(i,0), i = 1, imax-1)
write (o_unit,10) (egt(imax,j), j = 1, jmax-1)
write (o_unit,10) (egt(0,j), 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) (egtss(i,jmax), i = 1, imax-1)
write (o_unit,10) (egtss(i,0), i = 1, imax-1)
write (o_unit,10) (egtss(imax,j), j = 1, jmax-1)
write (o_unit,10) (egtss(0,j), 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) (egt10(i,jmax), i = 1, imax-1)
write (o_unit,11) (egt10(i,0), i = 1, imax-1)
write (o_unit,11) (egt10(imax,j), j = 1, jmax-1)
write (o_unit,11) (egt10(0,j), 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) (egt(i,j), 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) (egt(i,j)-egtss(i,j), 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) (egtss(i,j), 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) (egt10(i,j), i = 1, imax-1)
49 continue
write (o_unit,fmt="(' ')")
write (o_unit,*)
write (o_unit,*) '**Averages - Field'
write (o_unit,*) ' Total salt/creep susp PM10 '
write (o_unit,*) ' egt egtss egt10'
write (o_unit,*) ' -----------------kg/m^2--------------------'
write (o_unit,15) aegt, aegt-aegtss, aegtss, aegt10
write (o_unit,*)
write (o_unit,*) '**Averages - Crossing Boundaries '
write (o_unit,*) 'Location Total Salt/Creep Susp PM10'
write (o_unit,*) '--------------------kg/m----------------------'
write (o_unit,21) topt+topss, topt, topss, top10
write (o_unit,22) bott+botss, bott, botss, bot10
write (o_unit,23) ritt+ritss, ritt, ritss, rit10
write (o_unit,24) lftt+lftss, lftt, lftss, lft10
write (o_unit,*)
write (o_unit,*) ' Comparision 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
!^^^tmp out
! write (o_unit,*) 'lx=', lx, 'ly=', ly,'tot',tot
!^^^end tmpout
! additional output statements for easy shell script parsing
write (o_unit,*)
! write deposition as positive numbers
write (o_unit,17) aegt, aegt-aegtss, aegtss, aegt10
17 format (' repeat of total, salt/creep, susp, PM10:', 3f12.4,f12.6)
! 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)')
50 format (1x,' Leaving Field Grid Cells - Total egt(kg/m^2)')
60 format (1x,' Leaving Field Grid Cells - Salt/Creep egt-egtss(kg/m&
&^2)')
70 format (1x,' Leaving Field Grid Cells - Suspension egtss(kg/m^2) &
&')
80 format (1x,' Leaving Field Grid Cells - PM10 egt10(kg/m^2)')
10 format (1x, 500f12.4)
11 format (1x, 500f12.6)
15 format (1x, 3(f12.4,2x), f12.6)
16 format (1x, 2(f13.4,2x),2x, f13.4)
21 format (1x, 'top ', 1x, 4(f9.2,1x))
22 format (1x, 'bottom', 1x, 4(f9.2,1x))
23 format (1x, 'right ', 1x, 4(f9.2,1x))
24 format (1x, 'left ', 1x, 4(f9.2,1x))
endif !if (btest(am0efl,1)) then
!Erosion summary - total, salt/creep, susp, pm10
!(deposition values are positive - erosion values are negative)
if (btest(am0efl,0)) then
call caldatw (da, mo, yr) !get day, month and year
write (UNIT=o_E_unit,FMT="(' ',i5,i3,i3,' ')",ADVANCE="NO") &
& yr, mo, da
write (UNIT=o_E_unit,FMT="(4(f12.6),' ')",ADVANCE="YES") &
& aegt, (aegt-aegtss), aegtss, aegt10
endif
! end of plot section
return
end