! ! ! Program to simulate wind-blown PM10 emissions using WEPS' standalone ! EROSION submodel, WRF meteorology (via MCIP), and soil database ! "synthesize" by Guanglong Feng from the STATSGO database. ! program erosion_driver_v6a ! Description: ! ! ! Method: ! ! ! Requirements: ! ! ! Current code owner: Serena H. Chung ! ! History: ! ! Version Date Comment ! ------- ------- ------- ! 0.1 12/2008 ! 0.2b 09/2009 nmap as a dimension ! 0.2c 10/2009 module erosion_default; ! use ZRUF from METCRO2D or 10 mm for aerodynamic roughness ! ???SOIM1 for soil moisture in soil layer ! 0.2e 11/2009 correction on very fine sand fraction, asfvfs ! 0.2f 11/2009 correct failure to have output emissions include all mapunits ! 0.2g 11/2009 using zo_default ! 0.2h 10/2010 for AIRPACT3 domain in 12-km, force skipping the ! some western part. ! v4 reserved for fullweps ! 0.5a 11/2011 use MODIS LAI as an on/off switch to run EROSION ! 0.6a 01/2011 use MODIS NDVI (not LAI) as an on/off switch to run EROSION ! use Cropland Data Layer instead of BELD3 for cropland fraction ! wind speed adjustment factor, soil moisture ! ! Code Description: ! Language: Fortran 90, free form ! ! ! Declarations: ! Modules used: use soil_list use erosion_default implicit none ! Local scalar and array parameters !...netCDF library include 'netcdf.inc' !...Models-3 I/O API library include 'PARMS3.EXT' include 'IODECL3.EXT' include 'FDESC3.EXT' !...Name of this program character (len=16), parameter :: progname='erosion_drive_v6a' ! Common blocks !...Global common blocks for WEPS erosion code include 'p1werm.inc' include 'c1gen.inc' include 'p1const.inc' include 'm1sim.inc' include 'm1flag.inc' include 'm1subr.inc' include 'm1geo.inc' include 'b1glob.inc' include 'c1glob.inc' include 'd1glob.inc' include 's1layr.inc' include 's1phys.inc' include 's1dbh.inc' include 's1agg.inc' include 's1surf.inc' include 's1sgeo.inc' include 'h1db1.inc' include 'w1wind.inc' include 'w1pavg.inc' ! ! include 'file.fi' include 'erosion/m2geo.inc' include 'w1clig.inc' !Requires yrly average precip..... ! include 'util/misc/f2kcli.fi' !declarations for f2k commandline functions include 'command.inc' !...Local common blocks for WEPS erosion code include 'erosion/p1erode.inc' !obtain access to SURF_UPD_FLG variable ! Local scalars !..file names character (len=namlen3) :: infilename, outfilename !...For Models-3 I/O API file ID? integer :: logdev !...Date and time of data to read from and write to Models-3 I/O API !...files integer :: readDate, readTime, writeDate, writeTime integer :: runlen, stdate, sttime !...Time step in input file integer :: tstep !...Dimension variables integer :: nx, ny, ntime !...Grid dimension [=] m and area [=] m2 real :: dx, dy, area !...for using netcdf routines integer :: ncid, rcode, ndims, nvars, ngatts,unlimdimid integer :: dimid_lon, dimid_lat, dimid_layer,dimid_mapunit integer :: nlon,nlat,nlayer_total,nmapunit_max character (len = vname_char) :: varname !...input real :: wspdfac ! wind speed adjustment factor real :: soilm ! surface-soil moisture content real :: ustartb ! static threshold friction velocity for bare soil !...Counters, indices, misc integer :: i,j, k, itime, itime_start, itime_end, iday, s, nday integer :: iprop, ivar, imap !...soil data file name character(len=120) :: soil_file !...Misc character (len=namlen3*40) :: varlist character (len=80) :: mesg, vardesc integer :: status ! Local arrays real :: emis_hourly(1:24) real :: all_thickness(nlayer_max) real, allocatable :: tmp1d(:), tmp2d(:,:), tmp3d(:,:,:), tmp4d(:,:,:,:) real, allocatable :: WSPD10(:,:,:), WDIR10(:,:,:), DENS(:,:,:) real, allocatable :: ZRUF(:,:), SOIM1(:,:,:) real, allocatable :: soil_data(:,:,:,:,:), cropfrac_data(:,:) real, allocatable :: EMIS_PM10(:,:,:) real, allocatable :: area_by_grid(:,:,:) real, allocatable :: ndvi_data(:,:) integer, allocatable :: nmapunit_by_grid(:,:), nlayer_by_grid(:,:,:) ! External functions !...get environment variable as integer or real integer, external :: ENVINT real , external :: ENVREAL ! Additional WEPS stuff integer :: julday !utility date function real :: min_erosion_awu !Minimum erosiove wind speed (m/s) !to evaluate for erosion loss integer :: dt(8) character(len=3) :: mstring common / datetime / dt, mstring integer :: xplot character (len=112) :: xcharin(30) real :: xin(30) common /plot/ xplot, xcharin, xin integer :: init_sbemit_from_grid common /init_sbemit/ init_sbemit_from_grid ! End header --------------------------------------------------------- ! -------------------------------------------------------------------- ! Get simulation time information from environment variables ! (using Models-3 I/O API library) ! -------------------------------------------------------------------- ! Get start date of simulation vardesc = 'Scenario Starting Date (YYYYDDD)' stdate = envint( 'ERO_STDATE', vardesc, stdate, status ) if ( status .ne. 0 ) write( logdev, '(5x, a)' ) vardesc ! Get start time of simulation vardesc = 'Scenario Starting Time (HHMMSS)' sttime = envint( 'ERO_STTIME', vardesc, sttime, status ) if ( status .ne. 0 ) write( logdev, '(5x, a)' ) vardesc ! Get during of simulation (in hours) vardesc = 'Scenario Run Duration (HHMMSS)' runlen = envint( 'ERO_RUNLEN', vardesc, runlen, status ) if ( status .ne. 0 ) write( logdev, '(5x, a)' ) vardesc ! Get timestep of simulation/output ! (currently requires hourly ) vardesc = 'Scenario Output Time Step (HHMMSS)' tstep = envint( 'ERO_TSTEP', vardesc, tstep, status ) if ( status .ne. 0 ) write( logdev, '(5x, a)' ) vardesc if ( tstep /= 010000 ) then stop 'Program currently requires tstep=010000' end if ntime = runlen/tstep+1 if ( mod(ntime-1,24) /= 0 ) then stop 'ntime - 1 must be divisible by 24 ' end if nday = (ntime-1)/24 ! get input parameters from environment variables wspdfac = ENVREAL( 'WSPDFAC' , 'windspeed adjustment factor' , 1.0 , status ) soilm = ENVREAL( 'SOILM' , 'surface soil water content' , 0.01, status ) ustartb = ENVREAL( 'USTARTB' , 'static threshold friction vel', 0.2 , status ) ! obtain the name of the soil data file from environment variable call getenv ("SOIL_FILE",soil_file) ! -------------------------------------------------------------------- ! Open MCIP METCRO2D and METCRO3D files to get domain/grid ! information and to get wind speed, wind direction, and ! air density (from MM5 or WRF via MCIP) ! -------------------------------------------------------------------- ! Initialization of the Models-3 I/O API file logdev = init3() ! initialization returns unit # for log ! Open METCRO2D file for reading if ( .not. open3( 'MET_CRO_2D', fsread3, trim(progname) ) ) then mesg = "Could not open file MET_CRO_2D for input" call m3exit ( trim(progname), 0, 0, mesg, 2 ) end if ! Get file description if ( .not. desc3('MET_CRO_2D') ) then mesg = "Could not get file descriptionf for "//trim(infilename) call m3exit ( trim(progname), 0, 0, mesg, 2 ) end if ! Save dimension information, etc !... number of rows and columns nx = ncols3d ny = nrows3d !...time step tstep = tstep3d ! HHMMSS !...grid cell size in meters dx = XCELL3D dy = YCELL3D area = dx*dy ! meter^2 ! Open METCRO3D file for reading; ! This file contains meteorology if ( .not. open3( 'MET_CRO_3D', fsread3, trim(progname) ) ) then mesg = "Could not open file MET_CRO_3D for input" call m3exit ( trim(progname), 0, 0, mesg, 2 ) end if ! Read in meteorology data one time step a time readDate = stdate readTime = sttime allocate ( WSPD10(nx,ny,ntime), WDIR10(nx,ny,ntime), DENS(nx,ny,ntime) ) allocate ( ZRUF(nx,ny), SOIM1(nx,ny,ntime) ) allocate (tmp2d(nx,ny)) ! surface roughness length [=] m ! (assume time-independent, thus only read in the first time step) if ( read3( 'MET_CRO_2D', 'ZRUF', 1,readDate,readTime & , tmp2d) ) then ZRUF(1:nx,1:ny) = tmp2d(1:nx,1:ny) else mesg = "Could not read from MET_CRO_2D variable ZRUF" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if do itime = 1, ntime ! wind speed at 10 m [=] m/s WSPD10 if ( read3( 'MET_CRO_2D', 'WSPD10', 1,readDate,readTime & , tmp2d) ) then WSPD10(1:nx,1:ny,itime) = tmp2d(1:nx,1:ny) else mesg = "Could not read from MET_CRO_2D variable WSPD10" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! wind direction at 10 m [=] degrees WDIR10 if ( read3( 'MET_CRO_2D', 'WDIR10', 1,readDate,readTime & , tmp2d) ) then WDIR10(1:nx,1:ny,itime) = tmp2d(1:nx,1:ny) else mesg = "Could not read from MET_CRO_2D variable WDIR10" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if !!$ ! volumetric soil moisture in top cm [=] m3/m3 SOIM1 !!$ if ( read3( 'MET_CRO_2D', 'SOIM1', 1,readDate,readTime & !!$ , tmp2d) ) then !!$ SOIM1(1:nx,1:ny,itime) = tmp2d(1:nx,1:ny) !!$ else !!$ mesg = "Could not read from MET_CRO_2D variable SOIM1" !!$ call m3exit( trim(progname), 0, 0, mesg, 2 ) !!$ end if ! air density in the first layer [=] kg/m3 if ( read3( 'MET_CRO_3D', 'DENS', 1,readDate,readTime & , tmp2d) ) then DENS(1:nx,1:ny,itime) = tmp2d(1:nx,1:ny) else mesg = "Could not read from MET_CRO_3D variable DENS" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! increment to the next time step call nextime ( readDate, readTime, tstep ) end do ! itime = 1, ntime ! Close MCIP files if ( .not. close3('MET_CRO_2D') ) then mesg = "Could not close file MET_CRO_2D" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if if ( .not. close3('MET_CRO_3D') ) then mesg = "Could not close file MET_CRO_3D" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if deallocate (tmp2d) ! -------------------------------------------------------------------- ! Read in soil database ! -------------------------------------------------------------------- ! open soil database file for reading rcode = nf_open(soil_file,nf_nowrite,ncid) if ( rcode .ne. nf_noerr ) then write(6,*)'Error out of nf_open for file '//trim(soil_file) write(6,*)nf_strerror(rcode) stop end if ! inquire for some basic "header" information from soil data file rcode = nf_inq( ncid, ndims, nvars, ngatts, unlimdimid ) write(6,*) write(6,*) 'Output of nf_inq ' write(6,*) ' ndims = ', ndims write(6,*) ' nvars = ', nvars write(6,*) ' ngatts = ', ngatts write(6,*) ' unlimdimid = ', unlimdimid write(6,*) ! get dimension size of mapunit !...get dimension id for mapunit rcode = nf_inq_dimid( ncid,'mapunit', dimid_mapunit ) !...get dimension size for mapunit rcode=nf_inq_dimlen( ncid, dimid_mapunit, nmapunit_max ) ! do some consistency check on dimension size rcode = nf_inq_dimid( ncid,'lon', dimid_lon ) rcode=nf_inq_dimlen( ncid, dimid_lon, nlon) rcode = nf_inq_dimid( ncid,'lat', dimid_lat ) rcode=nf_inq_dimlen( ncid, dimid_lat, nlat) rcode = nf_inq_dimid( ncid,'layer', dimid_layer ) rcode=nf_inq_dimlen( ncid, dimid_layer, nlayer_total) if ( nlon /= nx ) then write(*,*)'nx, nlon = ', nx, nlon, ' not equal' stop else if ( nlat /= ny ) then write(*,*)'ny, nlat = ', ny, nlat, ' not equal' stop else if ( nlayer_total /= nlayer_max ) then write(*,*)'nlayer_max,nlayer_total = ',nlayer_max,nlayer_total,' not equal' stop end if ! read in soil data allocate( soil_data(nx,ny,nlayer_max,nmapunit_max,nsoilprop) ) !...first, read in soil properties at are in the STATGSO soil !..data base or are derived from the STATSGO soil database loop_soilprop: do iprop = 1, nsoilprop rcode = nf_inq_varname ( ncid, iprop, varname ) ! some consistency check if (trim(varname) /= trim(soilpropname(iprop)) ) then write(*,*)'when iprop = ',iprop write(*,*)'varname = ',trim(varname) write(*,*)'soilpropname = ', trim(soilpropname(iprop)) write(*,*)'something wrong is order of variables' stop end if if ( iprop == iTHICK ) then allocate (tmp1d(nlayer_max)) rcode = nf_get_var_real (ncid,iprop,tmp1d) all_thickness(1:nlayer_max)=tmp1d(1:nlayer_max) deallocate (tmp1d) else if ( (iprop==iWEG) .or. (iprop==iSTEEPNESS) ) then allocate ( tmp3d(nx,ny,nmapunit_max) ) rcode=nf_get_var_real ( ncid, iprop, tmp3d ) soil_data(1:nx,1:ny,1:nlayer_max,1:nmapunit_max,iprop ) = 0.0 soil_data(1:nx,1:ny,1 ,1:nmapunit_max,iprop ) = tmp3d deallocate ( tmp3d ) else allocate ( tmp4d (nx,ny,nlayer_max,nmapunit_max ) ) rcode=nf_get_var_real ( ncid, iprop, tmp4d ) if ( rcode /= nf_noerr ) then write(6,*)'Error out of nf_put_var_get ' write(6,*)nf_strerror(rcode) stop end if soil_data(1:nx,1:ny,1:nlayer_max,1:nmapunit_max,iprop ) = tmp4d deallocate (tmp4d) end if end do loop_soilprop ! read in area by of each map unit in each grid iprop = iArea rcode = nf_inq_varname ( ncid, iprop, varname ) !...some consistency check if (trim(varname) /= trim(soilpropname(iprop)) ) then write(*,*)'when iprop = ',iprop write(*,*)'varname = ',trim(varname) write(*,*)'soilpropname = ', trim(soilpropname(iprop)) write(*,*)'something wrong is order of variables' stop end if allocate ( area_by_grid(nx,ny,nmapunit_max) ) rcode=nf_get_var_real ( ncid, iprop, area_by_grid ) ! read in number of map units in each grid iprop = iNMapUnit rcode = nf_inq_varname ( ncid, iprop, varname ) !...some consistency check if (trim(varname) /= trim(soilpropname(iprop)) ) then write(*,*)'when iprop = ',iprop write(*,*)'varname = ',trim(varname) write(*,*)'soilpropname = ', trim(soilpropname(iprop)) write(*,*)'something wrong is order of variables' stop end if allocate ( nmapunit_by_grid(nx,ny) ) rcode=nf_get_var_int ( ncid, iprop, nmapunit_by_grid ) ! read in number of layers in each map unit of each each grid iprop = iNLayer rcode = nf_inq_varname ( ncid, iprop, varname ) !...some consistency check if (trim(varname) /= trim(soilpropname(iprop)) ) then write(*,*)'when iprop = ',iprop write(*,*)'varname = ',trim(varname) write(*,*)'soilpropname = ', trim(soilpropname(iprop)) write(*,*)'something wrong is order of variables' stop end if allocate ( nlayer_by_grid(nx,ny,nmapunit_max) ) rcode=nf_get_var_int ( ncid, iprop, nlayer_by_grid ) ! Put layer thickness into soil_data array soil_data(1:nx,1:ny,1:nlayer_max,1:nmapunit_max,iTHICK) = 0.0 do i = 1, nx do j = 1, ny do imap = 1, nmapunit_by_grid(i,j) do k = 1, nlayer_by_grid(i,j,imap) soil_data(i,j,k,imap,iTHICK) = all_thickness(k) end do end do end do end do ! Close soil data base file rcode = nf_close(ncid) ! -------------------------------------------------------------------- ! Read in Cropland Fraction from Cropland Data Layer data ! -------------------------------------------------------------------- ! Open CDL data file file for reading if ( .not. open3( 'CROP_FILE', fsread3, trim(progname) ) ) then mesg = "Could not open file CROP_FILE for input" call m3exit ( trim(progname), 0, 0, mesg, 2 ) end if ! read in cropland fraction data allocate( cropfrac_data(nx,ny)) if ( .not. read3( 'CROP_FILE', 'cropfrac', 1,0,0 & , cropfrac_data ) ) then mesg = "Could not read from file CROP_DATA" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! Close CDL database file if ( .not. close3('CROP_FILE') ) then mesg = "Could not close file CROP_FILE" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! -------------------------------------------------------------------- ! Read in MODIS NDVI data ! -------------------------------------------------------------------- ! Open MODIS NDVI data file file for reading if ( .not. open3( 'NDVI_DATA', fsread3, trim(progname) ) ) then mesg = "Could not open file NDVI_DATA for input" call m3exit ( trim(progname), 0, 0, mesg, 2 ) end if ! read in MODIS NDVI allocate( ndvi_data(nx,ny)) if ( .not. read3( 'NDVI_DATA', 'NDVI', 1,0,0 & , ndvi_data ) ) then mesg = "Could not read from file NDVI_DATA" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! Close MODIS NDVI database file if ( .not. close3('NDVI_DATA') ) then mesg = "Could not close file NDVI_DATA" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! -------------------------------------------------------------------- ! Read in crop management database ! -------------------------------------------------------------------- ! -------------------------------------------------------------------- ! Call WEPS's erosion submodule to calculate dust emissions ! -------------------------------------------------------------------- allocate(EMIS_PM10(nx,ny,ntime)) EMIS_PM10(:,:,:) = 0.0 i_loop: do i = 1, nx j_loop: do j = 1, ny ! run EROSION only if MODIS NDVI value is < 0.15 if ( ndvi_data(i,j) > 0.15 ) cycle j_loop if ( ndvi_data(i,j) < -0.1 ) cycle j_loop ! run EROSION only if cropfrac_data fraction > 1% ! ( cropfrac_data(:,:) is the fraction of the grid cell ! that is cropland if ( cropfrac_data(i,j) < 0.01 ) cycle j_loop ! For AIRPACT3 domain in 12-km, force skipping the ! some western part. if ( nx == 95 .and. ny == 95 ) then if (i <= 22) cycle j_loop end if ! run EROSION only if there exit soil mapunits ! in the grid with WEG (wind erodibility group) ! not in category 8 if ( nmapunit_by_grid(i,j) == 0 ) cycle j_loop ! loop over number of map units imap_loop: do imap = 1, nmapunit_by_grid(i,j) ! for some map unit, # of soil layers is 0 (why?); ! skip those map units if ( nlayer_by_grid(i,j,imap) == 0 ) cycle imap_loop ! Loop one day at a time, using 1 hour timestep day_loop: do iday = 1, nday ! Input into WEPS's standalone erosion module ! +++ INITIALIZATION STUFF +++ !......determine date of Run call date_and_time(values=dt) !......determine month of year select case (dt(2)) case (1); mstring = "Jan" case (2); mstring = "Feb" case (3); mstring = "Mar" case (4); mstring = "Apr" case (5); mstring = "May" case (6); mstring = "Jun" case (7); mstring = "Jul" case (8); mstring = "Aug" case (9); mstring = "Sep" case (10); mstring = "Oct" case (11); mstring = "Nov" case (12); mstring = "Dec" case default; mstring = "???" end select !...... MISC awzypt = 300.0 !Requires yrly average precip. (Hagen's best estimate for us to use) saeinp_daysim = 0 !initialize WEPS variables saeinp_jday = 0 !not used in the standalone submodel min_erosion_awu = 5.0 !default minimum erosive wind speed xgdpt = 0 !default grid spacing values if not specified ygdpt = 0 !on the commandline erod_interval = 0 !do not overide default surface updating interval SURF_UPD_FLG = 1 !enable erosion submodel surface updating by default !......Set the day, month and year for the "single day event" !......and get the "Julian Date" for it. !......(not applicable here, but do it anyway) am0jd = julday(1,1,0001) !......"initialization" flag !......(must be true for standalone EROSION runs am0eif = .true. !......"print" flag --- overwrite later????? am0efl = 3 !......to allow for initialization EROSION's subroutine SBEMIT init_sbemit_from_grid = 0 ! +++ SIMULATION REGION +++ !......Coordinates of simulation region amxsim(1,1) = 0. ! x1 amxsim(2,1) = 0 ! y1 amxsim(1,2) = dx ! x2 amxsim(2,2) = dy ! y2 !......Simulation orientation angle (degrees from North) amasim = 0.0 ! +++ ACCOUNTING REGION +++ !......Number of accounting regions (must always be 1 for now) !......and their coordinates (must match the simulation region) nacctr = 1 amxar(1,1,1) = amxsim(1,1) amxar(2,1,1) = amxsim(2,1) amxar(1,2,1) = amxsim(1,2) amxar(2,2,1) = amxsim(2,2) ! +++ BARRIERS +++ !......Number of barriers nbr = 0 ! +++ SUBREGIONS +++ !......Number of subregions (must always be 1 for now) !......and their coordinates (must match the simulation region) nsubr = 1 amxsr(1,1,1) = amxsim(1,1) amxsr(2,1,1) = amxsim(2,1) amxsr(1,2,1) = amxsim(1,2) amxsr(2,2,1) = amxsim(2,2) !......For the rest of the input , set for subregion 1 s = 1 ! +++ BIOMASS +++ !......"Overall" height of all standing biomass, !.......both crop and residue vegetation (meters) abzht(s) = abzht_default !......Average height of growing crop (meters) aczht(s) = aczht_default !......Growing crop stem area index (m2/m2) acrsai(s) = acrsai_default !......Growing crop leaf area index (m2/m2) acrlai(s) = acrlai_default !......Residue stem area index (m2/m2) adrsaitot(s) = adrsaitot_default !......Residue leaf area index (m2/m2) adrlaitot(s) = adrlaitot_default !......Growing crop row spacing (meters) !......(use value of 0.0 if not planted in rows, !......e.g. broadcast seeded) acxrow(s) = acxrow_default !......Specify seed location (0=furrow,1=ridge) !......(value doesn't matter if no rideges exist) ac0rg(s) = ac0rg_default !......Flat biomas cover (m2/m2) abffcv(s) = abffcv_default ! +++ SOIL ++++ !......Number of soil layers nslay(s) = 0 layer_loop1: do k = 1,nlayer_max if ( abs(soil_data(i,j,k,imap,iTHICK)) < 0.0001 ) then nslay(s) = k-1 exit layer_loop1 end if end do layer_loop1 if ( nslay(s) > 0 ) then layer_loop2 : do k = 1, nslay(s) ! layer thickness (mm) aszlyt(k,s) = soil_data(i,j,k,imap,iTHICK)*1000. ! bulk density (g/cm3) asdblk(k,s) = soil_data(i,j,k,imap,iBDM) ! fraction of sand content (g/g) asfsan(k,s) = soil_data(i,j,k,imap,iSAND)*0.01 ! fraction of very fine sand in soil layer ! estimate from Feng and Sharratt, ! Earth Surface Processes and Landforms ! 34, 1461-1468 (2009) ! (11/13/2009) asfvfs(k,s) = soil_data(i,j,k,imap,iSAND)*0.66*0.01 ! fraction of silt content (g/g) asfsil(k,s) = soil_data(i,j,k,imap,iSILT)*0.01 ! fraction of clay content (g/g) asfcla(k,s) = soil_data(i,j,k,imap,iCLAY)*0.01 ! rock volume in soil layer (m3/m3) asvroc(k,s) = asvroc_default ! average aggregate density of soil layer (Mg/m3) asdagd(k,s) = soil_data(i,j,k,imap,iSAGd) ! average dry aggregate stability (ln(J/kg) aseags(k,s) = soil_data(i,j,k,imap,iSAGs) ! geometric mean diameter of aggregates (mm) aslagm(k,s) = soil_data(i,j,k,imap,iAGMD) ! 0.4 ! minimum aggreagate size (mm) aslagn(k,s) = 0.01 ! maximum aggregate size (mm) aslagx(k,s) = soil_data(i,j,k,imap,iAMAX) ! geometric stand deviation of aggregate size as0ags(k,s) = soil_data(i,j,k,imap,iAGSD) end do layer_loop2 end if ! +++ SOIL SURFACE +++ !......surface crust fraction (m2/m2) asfcr(s) = asfcr_default !......surface crust thickness (mm) aszcr(s) = aszcr_default !......fraction of crusted surface with loose material on top !......of crust (m2/m2) asflos(s) = asflos_default !......mass of loose material on top of crust (kg/m2) asmlos(s) = asmlos_default !......density of soil crust (Mg/m2) asdcr(s) = asdcr_default !......dry crust stability (ln(J/kg)) asecr(s) = asecr_default !......allmaras random roughnes (mm) aslrr(s) = aslrr_default !......ridge height (mm) aszrgh(s) = aszrgh_default !......ridge spacing (mm) asxrgs(s) = asxrgs_default !......ridge width (mm) asxrgw(s) = asxrgw_default !......ridge orientation (degrees) asargo(s) = asargo_default !......dike spacing (if no dikes, specify 0.0) (mm) asxdks(s) = asxdks_default ! ++++ HYDROLOGY ++++ !......snow depth (mm) ahzsnd(s) = 0.0 do k = 1, nslay(s) ! wilting point water content of soild layer (g/g) ahrwcw(k,s) = soil_data(i,j,k,imap,iPWP) ! current water content of soil layer (g/g) ahrwca(k,s) = soil_data(i,j,k,imap,iPWP)*0.1 !SOIM1(i,j,1) ! all zeros end do !......surface layer water content (g/g) !......hourly basis ahrwc0(:,s) = soilm !ahrwc0_default ! +++ WEATHER ++++ !......number of intervals per day to run EROSION !......(must be 24 for now) ntstep = 24 !......time counter itime_start = (iday-1)*ntstep+1 itime_end = iday*ntstep !......air density (kg/m3) awdair = sum(DENS(i,j,itime_start:itime_end))/ntstep !......wind direction, measured clockwise from north (degrees) awdir(1:ntstep) = WDIR10(i,j,itime_start:itime_end) !......anemometer height (m) anemht = 10.0 !......aerodynamic roughness at anemometer site (mm) awzzo = zo_default ! ZRUF(i,j)*1000. ! not really at anemometer site !......zo location flag !......0 - at weather location - zo is constant !......1 - on field location - zo various based on field surface wzoflg = 0 !......wind speed awu(1:24) = WSPD10(i,j,itime_start:itime_end)*wspdfac !.........maximum windspeed awudmx = 0.0 do itime = 1,ntstep awudmx = max (awudmx,awu(itime)) end do !.......flag for output for plotting xplot = 0 xin(1:30) = 0.0 ! +++ CALL EROSION +++ !......initialize erosion code, create grid, etc: !......(must come after sim field size, & no. subr specified) call erodinit !...... emis_hourly(1:ntstep) = 0.0 ! ntstep has to be 24 for now call erosion (min_erosion_awu,emis_hourly) ! +++ SAVE EROSION OUTPUT !...... PM10 emissions in g/s EMIS_PM10(i,j,itime_start:itime_end)= & EMIS_PM10(i,j,itime_start:itime_end) & + emis_hourly(:)*area_by_grid(i,j,imap)*cropfrac_data(i,j)*1000. end do day_loop end do imap_loop end do j_loop end do i_loop ! -------------------------------------------------------------------- ! Output dust emissions in Models-3 I/O API file format ! -------------------------------------------------------------------- ! "Template" for the output Models-3 I/O API file nvars3d = 1 vname3d(1) = 'DUSTPM10' vdesc3d(1) = 'DUSTPM10' vtype3d(1) = m3real units3d(1) = 'g/s' !!$ vname3d(2) = 'CROPFRAC' !!$ vdesc3d(2) = 'Fraction of cropland assumed' !!$ vtype3d(2) = m3real !!$ units3d(2) = 'fraction' fdesc3d(:) = "" fdesc3d(1) = "created on aeolus.wsu.edu by /home/schung/WEPS/src/erosion_v6a.exe" tstep3d = tstep ! Open the new file if ( .not. open3( 'OUTFILE', FSCREA3, trim(progname) ) ) then mesg = 'Could not open file OUTFILE for output' call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! Put list of variables in global attribute of netCDF file do ivar = 1, nvars3d i = ivar*namlen3 varlist(i-15:i)= vname3d(ivar)(1:namlen3) end do if ( .not. wrattc ( 'OUTFILE', ALLVAR3, "VAR-LIST" , varlist(1:namlen3*nvars3d) )) then mesg = 'Could not write attribute VAR-LIST' call m3exit( trim(progname), 0, 0, mesg, 2 ) end if writeDate = stdate writeTime = sttime do itime = 1, ntime ! Write out to file ivar = 1 if ( .not. ( write3('OUTFILE', trim(vname3d(ivar)),writeDate,writeTime,EMIS_PM10(:,:,itime ) ) ) ) then mesg = "Could not write to OUTFILE variable "//trim(vname3d(ivar)) call m3exit( trim(progname), 0, 0, mesg, 2 ) end if !!$ ! Cropfrac_data fraction is actually time-invariant, but what the hack... !!$ ivar = 2 !!$ if ( .not. ( write3('OUTFILE', trim(vname3d(ivar)),writeDate,writeTime,cropfrac_data(:,:) ) ) ) then !!$ mesg = "Could not write to OUTFILE variable "//trim(vname3d(ivar)) !!$ call m3exit( trim(progname), 0, 0, mesg, 2 ) !!$ end if call nextime ( writeDate , writeTime , tstep ) ! an I/O API function end do mesg ='Exiting Models-3 I/O API' call m3exit( trim(progname), 0, 0, mesg, 0 ) write(*,*) write(*,*)'Program Successful !!! ' write(*,*) end program erosion_driver_v6a