!**********************************************************************
! subroutine sb1out
!**********************************************************************
subroutine sb1out (jj, nn, hr, ws, wdir, o_unit)
use p1werm_def
use constants_def
! include 'p1const.inc'
use anemometer_def
use datetime_def
use c1gen_def
use m1sim_def
use m1geo_def
use m2geo_def
use w1clig_def
use gridmod
use soilmod
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 's1surf.inc'
! + + + LOCAL COMMON BLOCKS + + +
!
include 's1dbh.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 m, icsr, x, y
integer 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
! 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
write (o_unit,*) 'Date of run: ', datetimestr
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", SoilRegionData(1)%SoilLayerData(1)%fracRock, " (m^3/m^3)"
write (o_unit,fmt="(a,a,f5.2,a)") "Initial soil ", &
& "mass fraction in surface layer < 0.10 mm ", soilFracDiamLt10ic, " (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 ", SoilRegionData(1)%asoilPM10FracAbraSusp
write (o_unit,fmt="(a,f5.2,a)") &
& "Soil fraction PM10 in emitted suspension ", SoilRegionData(1)%asoilPM10FracEmitSusp
write (o_unit,fmt="(a,f5.2,a)") &
& "Soil fraction PM10 in saltation breakage suspension ",SoilRegionData(1)%asoilPM10FracSaltBrkSusp
write (o_unit,fmt="(a,f5.2,a)") &
& "Coefficient of abrasion of aggregates ", SoilRegionData(1)%acoefAbraAgg
write (o_unit,fmt="(a,f5.2,a)") &
& "Coefficient of abrasion of crust ", SoilRegionData(1)%acoefAbraCrst
!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)") (grid(i,j)%fricVelc, 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)") (grid(i,j)%thrsFricVelc, 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)") (grid(i,j)%thrsFricVelcTrap, 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)") (grid(i,j)%randRoug, 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)") (grid(i,j)%ridgHght, 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)") (grid(i,j)%soilLayrRock, 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)") (grid(i,j)%soilFracDiamLt1, 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)") (grid(i,j)%soilFracDiamLt10, 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)") (grid(i,j)%soilFracDiamLt84, 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)") (grid(i,j)%soilFracDiamLt200, 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)") (grid(i,j)%soilFracDiamLt84mn, 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)") (grid(i,j)%soilMassAvalAggSurf, 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)") (grid(i,j)%soilMassAvalDelt, 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)") (grid(i,j)%soilCrstThck, 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)") (grid(i,j)%soilCrstFrac, 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)") (grid(i,j)%soilLoosCovFrac, 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)") (grid(i,j)%soilLoosMass, 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) (soilLossTot(k,n),k=1,(imax-1),m)
! write (o_unit,22) (soilLossSusp(k,n),k=1,(imax-1),m)
! write (o_unit,23??) ((soilLossSusp(k,n)/(soilLossTot(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) (soilFracDiamLt1(k,n),k=1,(imax-1),m)
! write (o_unit,23) (soilFracDiamLt10(k,n),k=1,(imax-1),m)
! write (o_unit,24) (soilFracDiamLt84(k,n),k=1,(imax-1),m)
! write (o_unit,35) (soilFracDiamLt200(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, soilFracDiamLt10ic, SoilRegionData(1)%SoilLayerData(1)%fracRock !edit ljh 1-22-05
! write (o_unit,42) SoilRegionData(1)%acoefAbraAgg, SoilRegionData(1)%acoefAbraCrst, awzypt
! write (o_unit,10) SoilRegionData(1)%asoilPM10FracAbraSusp, SoilRegionData(1)%asoilPM10FracEmitSusp, SoilRegionData(1)%asoilPM10FracSaltBrkSusp
! 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, 'soilLossTot=', 20f6.2)
! 22 format (1x, 'soilLossSusp=', 20f6.2)
! 13 format (1x, 'soilFracDiamLt1= ', 20f7.4)
! 23 format (1x, 'soilFracDiamLt10= ', 20f7.3)
! 24 format (1x, 'soilFracDiamLt84= ', 20f7.3)
! 12 format (1x, 'soilLayrRock=', 20f7.3) !edit ljh 1-22-05
! 35 format (1x, 'soilFracDiamLt200=', 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,' soilFracDiamLt10ic =',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)
subroutine sb2out (jj, nn, hr, o_unit)
use p1werm_def
use constants_def
use m1sim_def
use m2geo_def
use gridmod
IMPLICIT NONE
!
! + + + 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
real 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'
!
! + + + 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)") (grid(i,j)%soilLossTot, 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)")(grid(i,j)%soilLossTot-grid(i,j)%soilLossSusp,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)") (grid(i,j)%soilLossSusp, 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)") (grid(i,j)%soilLossPM10, 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)") (grid(i,j)%randRoug, 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)") (grid(i,j)%ridgHght, 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)") (grid(i,j)%soilLayrRock, 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)") (grid(i,j)%soilFracDiamLt1, 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)") (grid(i,j)%soilFracDiamLt10, 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)") (grid(i,j)%soilFracDiamLt84, 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)") (grid(i,j)%soilFracDiamLt200, 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)") (grid(i,j)%soilFracDiamLt84mn, 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)") (grid(i,j)%soilMassAvalAggSurf, 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)") (grid(i,j)%soilMassAvalDelt, 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)") (grid(i,j)%soilCrstThck, 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)") (grid(i,j)%soilCrstFrac, 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)") (grid(i,j)%soilLoosCovFrac, 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)") (grid(i,j)%soilLoosMass, 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) + soilLossTot(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) + grid(i,j)%soilLossTot/(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) (grid(k,n)%soilLossTot,k=1,(imax-1),m)
! write (o_unit,22) (grid(k,n)%soilLossSusp,k=1,(imax-1),m)
! write (o_unit,23) ((grid(k,n)%soilLossSusp)/(grid(k,n)%soilLossTot+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, 'soilLossTot= ', 20f6.2)
! 22 format (1x, 'soilLossSusp= ', 20f6.2)
! 23 format (1x, 'soilLossSusp/soilLossTot=', 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)
use p1werm_def
use datetime_def
use m1geo_def
use gridmod
IMPLICIT NONE
! +++ PURPOSE +++
! To print output desired from standalone EROSION submodel
! +++ ARGUMENT DECLARATIONS +++
integer o_unit, o_E_unit
integer i, j
real asoilLossTot, asoilLossSusp, asoilLossPM10
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 'm1flag.inc'
! + + + LOCAL COMMON BLOCKS
! include 'erosion/e2erod.inc'
include 'erosion/m2geo.inc'
!
integer x, y
! ++++ ARGUMENT DEFINITIONS +++
! +++ SUBROUTINES CALLED+++
! plotout.for
! ++++ LOCAL VARIABLES +++
! +++ END SPECIFICATIONS +++
! Calculate Averages Crossing Borders
! top border
asoilLossTot = 0.0
asoilLossSusp = 0.0
asoilLossPM10 = 0.0
j = jmax
do 1 i = 1, imax-1
asoilLossTot = asoilLossTot + grid(i,j)%soilLossTot
asoilLossSusp = asoilLossSusp + grid(i,j)%soilLossSusp
asoilLossPM10 = asoilLossPM10 + grid(i,j)%soilLossPM10
1 continue
! calc. average at top border
topt = asoilLossTot/(imax-1)
topss = asoilLossSusp/(imax-1)
top10 = asoilLossPM10/(imax-1)
! bottom border
asoilLossTot = 0.0
asoilLossSusp = 0.0
asoilLossPM10 = 0.0
j = 0
do 2 i = 1, imax-1
asoilLossTot = asoilLossTot + grid(i,j)%soilLossTot
asoilLossSusp = asoilLossSusp + grid(i,j)%soilLossSusp
asoilLossPM10 = asoilLossPM10 + grid(i,j)%soilLossPM10
2 continue
! calc. average at bottom border
bott = asoilLossTot/(imax-1)
botss = asoilLossSusp/(imax-1)
bot10 = asoilLossPM10/(imax-1)
! right border
asoilLossTot = 0.0
asoilLossSusp = 0.0
asoilLossPM10 = 0.0
i = imax
do 3 j = 1, jmax-1
asoilLossTot = asoilLossTot + grid(i,j)%soilLossTot
asoilLossSusp = asoilLossSusp + grid(i,j)%soilLossSusp
asoilLossPM10 = asoilLossPM10 + grid(i,j)%soilLossPM10
3 continue
! calc. average at right border
ritt = asoilLossTot/(jmax-1)
ritss = asoilLossSusp/(jmax-1)
rit10 = asoilLossPM10/(jmax-1)
!
! left border
asoilLossTot = 0.0
asoilLossSusp = 0.0
asoilLossPM10 = 0.0
i = 0
do 4 j = 1, jmax-1
asoilLossTot = asoilLossTot + grid(i,j)%soilLossTot
asoilLossSusp = asoilLossSusp + grid(i,j)%soilLossSusp
asoilLossPM10 = asoilLossPM10 + grid(i,j)%soilLossPM10
4 continue
! calc. average at left border
lftt = asoilLossTot/(jmax-1)
lftss = asoilLossSusp/(jmax-1)
lft10 = asoilLossPM10/(jmax-1)
! calculate averages of inner grid points
asoilLossTot = 0.0
asoilLossSusp = 0.0
asoilLossPM10 = 0.0
do 5 j=1,jmax-1
do 5 i= 1, imax-1
asoilLossTot= asoilLossTot + grid(i,j)%soilLossTot
asoilLossSusp = asoilLossSusp + grid(i,j)%soilLossSusp
asoilLossPM10 = asoilLossPM10 + grid(i,j)%soilLossPM10
5 continue
tt = (imax-1)*(jmax-1)
asoilLossTot = asoilLossTot/tt
asoilLossSusp = asoilLossSusp/tt
asoilLossPM10 = asoilLossPM10/tt
! calculate comparision of boundary and interior losses
lx = amxsim(1,2) - amxsim(1,1)
ly = amxsim(2,2) - amxsim(2,1)
tot = asoilLossTot*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
write (o_unit,*) 'Date of run: ', datetimestr
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) (grid(i,jmax)%soilLossTot+grid(i,jmax)%soilLossSusp, i = 1, imax-1)
write (o_unit,10) (grid(i,0)%soilLossTot+grid(i,0)%soilLossSusp, i = 1, imax-1)
write (o_unit,10) (grid(imax,j)%soilLossTot+grid(imax,j)%soilLossSusp, j = 1, jmax-1)
write (o_unit,10) (grid(0,j)%soilLossTot+grid(0,j)%soilLossSusp, 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) (grid(i,jmax)%soilLossTot, i = 1, imax-1)
write (o_unit,10) (grid(i,0)%soilLossTot, i = 1, imax-1)
write (o_unit,10) (grid(imax,j)%soilLossTot, j = 1, jmax-1)
write (o_unit,10) (grid(0,j)%soilLossTot, 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) (grid(i,jmax)%soilLossSusp, i = 1, imax-1)
write (o_unit,10) (grid(i,0)%soilLossSusp, i = 1, imax-1)
write (o_unit,10) (grid(jmax,j)%soilLossSusp, j = 1, jmax-1)
write (o_unit,10) (grid(0,j)%soilLossSusp, 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) (grid(i,jmax)%soilLossPM10, i = 1, imax-1)
write (o_unit,11) (grid(i,0)%soilLossPM10, i = 1, imax-1)
write (o_unit,11) (grid(imax,j)%soilLossPM10, j = 1, jmax-1)
write (o_unit,11) (grid(0,j)%soilLossPM10, 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) (grid(i,j)%soilLossTot, 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) (grid(i,j)%soilLossTot-grid(i,j)%soilLossSusp, 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) (grid(i,j)%soilLossSusp, 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) (grid(i,j)%soilLossPM10, 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,*) ' soilLossTot soilLossSusp soilLossPM10'
write (o_unit,*) ' -----------------kg/m^2--------------------'
write (o_unit,15) asoilLossTot, asoilLossTot-asoilLossSusp, asoilLossSusp, asoilLossPM10
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) asoilLossTot, asoilLossTot-asoilLossSusp, asoilLossSusp, asoilLossPM10
17 format (' repeat of total, salt/creep, susp, PM10:', 3f12.4,f12.6)
! output formats
! 6 format (1x,' Passing Border Grid Cells - Total soilLossTot+soilLossSusp(kg/m)')
7 format (1x,' Passing Border Grid Cells - Salt/Creep soilLossTot(kg/m)')
8 format (1x,' Passing Border Grid Cells - Suspension soilLossSusp(kg/m)')
9 format (1x,' Passing Border Grid Cells - PM10 soilLossPM10(kg/m)')
! 50 format (1x,' Leaving Field Grid Cells - Total soilLossTot(kg/m^2)')
! 60 format (1x,' Leaving Field Grid Cells - Salt/Creep soilLossTot-soilLossSusp(kg/m&
! &^2)')
! 70 format (1x,' Leaving Field Grid Cells - Suspension soilLossSusp(kg/m^2) &
! &')
! 80 format (1x,' Leaving Field Grid Cells - PM10 soilLossPM10(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") &
& asoilLossTot, (asoilLossTot-asoilLossSusp), asoilLossSusp, asoilLossPM10
endif
! end of plot section
return
end