! ! ! Program to simulate wind-blown PM10 emissions using WEPS' standalone ! EROSION submodel, WRF meteorology (via MCIP), and soil database ! "synthesized" by Guanglong Feng from the STATSGO database. In this ! version, WEPS is run at a finer resolution than WRF. The emissions ! of WEPS at the finer grids are summed up to arrive at total emissions ! at WRF grids. ! program erosion_nest_v7a ! ! ! Description: ! ! ! Method: ! ! ! Requirements: ! ! ! Current code owner: Serena H. Chung ! ! History: ! ! Version Date Comment ! ------- ------- ------- ! 0.2c 10/2009 modified from erosion_nest.f90 ! use nmap as a dimension ! use ZRUF from METCRO2D or 10 mm for aerodynamic roughness ! 0.2d 11/2009 bug correction: needs to use met/transport grid ! dimension for field size, i.e., ! amxsim(1,2) = dx_met ! amxsim(2,2) = dy_met ! 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. ! ! 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 added as input from evironment variable ! 0.7a 02/2011 use AMSR-E Soil Moisture Product ! ! 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_nest_v7a' !...Minimum erodible fraction allowed to fun the erosion code real, parameter :: erodible_min= 0.01 ! 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 !...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_met, ny_met, nx_ero, ny_ero, ntime !...Grid dimension [=] m and area [=] m2 real :: dx_met, dy_met, area_met real :: dx_ero, dy_ero !...ratio of coarse adnn finer grisd integer :: dx_ratio, dy_ratio !...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 !...Counters, indices, misc integer :: i_start, j_start integer :: i_outer,j_outer, i_inner, j_inner, i, j integer :: 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(:,:,:,:,:) real, allocatable :: EMIS_PM10_met(:,:,:), EMIS_PM10_ero(:,:,:) real, allocatable :: area_by_grid(:,:,:) real, allocatable :: soilm_data(:,:) integer, allocatable :: nmapunit_by_grid(:,:), nlayer_by_grid(:,:,:) real, allocatable :: ndvi_data(:,:), cropfrac_data(:,:) ! External functions !...get environment variable as integer 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) ! -------------------------------------------------------------------- vardesc = 'Scenario Starting Date (YYYYDDD)' stdate = envint( 'ERO_STDATE', vardesc, stdate, status ) if ( status .ne. 0 ) write( logdev, '(5x, a)' ) vardesc vardesc = 'Scenario Starting Time (HHMMSS)' sttime = envint( 'ERO_STTIME', vardesc, sttime, status ) if ( status .ne. 0 ) write( logdev, '(5x, a)' ) vardesc vardesc = 'Scenario Run Duration (HHMMSS)' runlen = envint( 'ERO_RUNLEN', vardesc, runlen, status ) if ( status .ne. 0 ) write( logdev, '(5x, a)' ) vardesc 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 ) ! 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 MET_CRO_2D" call m3exit ( trim(progname), 0, 0, mesg, 2 ) end if ! Save dimension information, etc nx_met = ncols3d ny_met = nrows3d tstep = tstep3d ! HHMMSS dx_met = XCELL3D dy_met = YCELL3D area_met = dx_met*dy_met ! Open METCRO3D file for reading 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 readDate = stdate readTime = sttime allocate ( WSPD10(nx_met,ny_met,ntime), WDIR10(nx_met,ny_met,ntime), DENS(nx_met,ny_met,ntime) ) allocate ( ZRUF(nx_met,ny_met), SOIM1(nx_met,ny_met,ntime) ) allocate (tmp2d(nx_met,ny_met)) ! 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_met,1:ny_met) = tmp2d(1:nx_met,1:ny_met) 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_met,1:ny_met,itime) = tmp2d(1:nx_met,1:ny_met) 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_met,1:ny_met,itime) = tmp2d(1:nx_met,1:ny_met) 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_met,1:ny_met,itime) = tmp2d(1:nx_met,1:ny_met) !!$ 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_met,1:ny_met,itime) = tmp2d(1:nx_met,1:ny_met) 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_2D" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! -------------------------------------------------------------------- ! Open/Create the output dust emissions file of emissions ! in the met grid ! -------------------------------------------------------------------- ! Modify file "Template" of the output file from that of ! METCRO2D file for the output Models-3 I/O API file nvars3d = 1 vname3d(1) = 'DUSTPM10' vdesc3d(1) = 'DUSTPM10' vtype3d(1) = m3real units3d(1) = 'g/s' fdesc3d(:) = "" fdesc3d(1) = "created on aeolus.wsu.edu by /home/schung/WEPS/src/erosion_nest_v7a.exe" tstep3d = tstep ! Open the new file if ( .not. open3( 'OUTFILE1', FSCREA3, trim(progname) ) ) then mesg = 'Could not open file OUTFILE1 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_inner = ivar*namlen3 varlist(i_inner-15:i_inner)= vname3d(ivar)(1:namlen3) end do if ( .not. wrattc ( 'OUTFILE1', 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 ! -------------------------------------------------------------------- ! 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 ! Get file description and save dimension information if ( .not. desc3('CROP_FILE') ) then mesg = "Could not get file descriptionf for CROP_FILE" call m3exit ( trim(progname), 0, 0, mesg, 2 ) end if nx_ero = ncols3d ny_ero = nrows3d dx_ero = XCELL3D dy_ero = YCELL3D ! read in cropland fraction data allocate( cropfrac_data(nx_ero,ny_ero)) 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 ! -------------------------------------------------------------------- ! Open/Create the output dust emissions file of emissions ! the landuse data file grid ! -------------------------------------------------------------------- ! Modify file "Template" of the output file from that of ! landuse file for the output Models-3 I/O API file sdate3d = stdate stime3d = sttime tstep3d = tstep nvars3d = 1 vname3d(1) = 'DUSTPM10' vdesc3d(1) = 'DUSTPM10' vtype3d(1) = m3real units3d(1) = 'g/s' fdesc3d(:) = "" fdesc3d(1) = "created on aeolus.wsu.edu by /home/schung/WEPS/src/erosion_nest_v7a.exe" tstep3d = tstep ! Open the new file if ( .not. open3( 'OUTFILE2', FSCREA3, trim(progname) ) ) then mesg = 'Could not open file OUTFILE2 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_inner = ivar*namlen3 varlist(i_inner-15:i_inner)= vname3d(ivar)(1:namlen3) end do if ( .not. wrattc ( 'OUTFILE2', 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 ! -------------------------------------------------------------------- ! 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_ero ) then write(*,*)'nx_ero, nlon = ', nx_ero, nlon, ' not equal' stop else if ( nlat /= ny_ero ) then write(*,*)'ny_ero, nlat = ', ny_ero, 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_ero,ny_ero,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_ero,ny_ero,nmapunit_max) ) rcode=nf_get_var_real ( ncid, iprop, tmp3d ) soil_data(1:nx_ero,1:ny_ero,1:nlayer_max,1:nmapunit_max,iprop ) = 0.0 soil_data(1:nx_ero,1:ny_ero,1 ,1:nmapunit_max,iprop ) = tmp3d deallocate ( tmp3d ) else allocate ( tmp4d (nx_ero,ny_ero,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_ero,1:ny_ero,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_ero,ny_ero,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_ero,ny_ero) ) 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_ero,ny_ero,nmapunit_max) ) rcode=nf_get_var_int ( ncid, iprop, nlayer_by_grid ) ! Put layer thickness into soil_data array soil_data(1:nx_ero,1:ny_ero,1:nlayer_max,1:nmapunit_max,iTHICK) = 0.0 do i = 1, nx_ero do j = 1, ny_ero 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 i = -1; j = -1 ! to insure no mis-use of index later ! Close soil data base file rcode = nf_close(ncid) ! -------------------------------------------------------------------- ! 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_ero,ny_ero)) 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 AMSR-E Soil Moisture data ! -------------------------------------------------------------------- ! Open AMSR-E Soil Moisture data file file for reading if ( .not. open3( 'SOILM_DATA', fsread3, trim(progname) ) ) then mesg = "Could not open file SOILM_DATA for input" call m3exit ( trim(progname), 0, 0, mesg, 2 ) end if ! read in AMSR-E Soil Moisture allocate( soilm_data(nx_ero,ny_ero)) if ( .not. read3( 'SOILM_DATA', 'SoilM', 1,0,0 & , soilm_data ) ) then mesg = "Could not read from file SOILM_DATA" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! Close AMSR-E Soil Moisture database file if ( .not. close3('SOILM_DATA') ) then mesg = "Could not close file SOILM_DATA" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! -------------------------------------------------------------------- ! Call WEPS's erosion submodule to calculate dust emissions ! -------------------------------------------------------------------- dx_ratio = nx_ero/nx_met dy_ratio = ny_ero/ny_met allocate(EMIS_PM10_met(nx_met,ny_met,ntime)) allocate(EMIS_PM10_ero(nx_ero,ny_ero,ntime)) EMIS_PM10_met(:,:,:) = 0.0 EMIS_PM10_ero(:,:,:) = 0.0 i_loop_met: do i_outer = 1, nx_met j_loop_met: do j_outer = 1, ny_met ! For AIRPACT3 domain in 12-km, force skipping ! some western part. if ( nx_met == 95 .and. ny_met == 95 ) then if (i_outer <= 22) cycle j_loop_met end if i_start = dx_ratio*(i_outer-1)+1 j_start = dy_ratio*(j_outer-1)+1 i_loop_ero: do i_inner = i_start, i_start+dx_ratio-1 j_loop_ero: do j_inner = j_start, j_start+dy_ratio-1 ! run EROSION only if erodible fraction is greager than a minimum ! erodible fraction is the fraction of cropland as determined by ! Cropland Data Layer if (cropfrac_data(i_inner,j_inner) < erodible_min ) cycle j_loop_ero ! 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_inner,j_inner) == 0 ) cycle j_loop_ero ! for now, also skip grid cells in which there is no soil ! moisture data from AMSR-E, which seems to occur over ! non-erodible grid cells anyway if ( soilm_data(i_inner,j_inner) > 9998. ) cycle j_loop_ero ! run EROSION only if MODIS NDVI value is < 0.15 if ( ndvi_data(i_inner,j_inner) > 0.15 ) cycle j_loop_ero if ( ndvi_data(i_inner,j_inner) < -0.1 ) cycle j_loop_ero ! loop over number of map units imap_loop: do imap = 1, nmapunit_by_grid(i_inner,j_inner) ! for some map unit, # of soil layers is 0 (why?); ! skip those map units if ( nlayer_by_grid(i_inner,j_inner,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 !......(use the meteorological/transport grid dimension; 11/02/2009) amxsim(1,1) = 0. ! x1 amxsim(2,1) = 0 ! y1 amxsim(1,2) = dx_met amxsim(2,2) = dy_met !......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) = nlayer_by_grid(i_inner,j_inner,imap) layer_loop2 : do k = 1, nslay(s) ! layer thickness (mm) aszlyt(k,s) = soil_data(i_inner,j_inner,k,imap,iTHICK)*1000. ! bulk density (g/cm3) asdblk(k,s) = soil_data(i_inner,j_inner,k,imap,iBDM) ! fraction of sand content (g/g) asfsan(k,s) = soil_data(i_inner,j_inner,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_inner,j_inner,k,imap,iSAND)*0.66*0.01 ! fraction of silt content (g/g) asfsil(k,s) = soil_data(i_inner,j_inner,k,imap,iSILT)*0.01 ! fraction of clay content (g/g) asfcla(k,s) = soil_data(i_inner,j_inner,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_inner,j_inner,k,imap,iSAGd) ! average dry aggregate stability (ln(J/kg) aseags(k,s) = soil_data(i_inner,j_inner,k,imap,iSAGs) ! geometric mean diameter of aggregates (mm) aslagm(k,s) = soil_data(i_inner,j_inner,k,imap,iAGMD) ! 0.4 ! minimum aggreagate size (mm) aslagn(k,s) = 0.01 ! maximum aggregate size (mm) aslagx(k,s) = soil_data(i_inner,j_inner,k,imap,iAMAX) ! geometric stand deviation of aggregate size as0ags(k,s) = soil_data(i_inner,j_inner,k,imap,iAGSD) end do layer_loop2 ! +++ 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_inner,j_inner,k,imap,iPWP) ! current water content of soil layer (g/g) ahrwca(k,s) = soil_data(i_inner,j_inner,k,imap,iPWP)*0.1 !??? SOIM1(i_inner,j_inner,1) ! all zeros end do !......surface layer water content (g/g) !......hourly basis if ( asdblk(1,s) > 1.0 ) then ahrwc0(:,s) = soilm_data(i_inner,j_inner)/asdblk(1,s) else ahrwc0(:,s) = soilm_data(i_inner,j_inner) endif ! +++ 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_outer,j_outer,itime_start:itime_end))/ntstep !......wind direction, measured clockwise from north (degrees) awdir(1:ntstep) = WDIR10(i_outer,j_outer,itime_start:itime_end) !......anemometer height (m) anemht = 10.0 !......aerodynamic roughness at anemometer site (mm) awzzo = zo_default ! ZRUF(i_outer,j_outer)*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_outer,j_outer,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_ero(i_inner,j_inner,itime_start:itime_end) = & EMIS_PM10_ero(i_inner,j_inner,itime_start:itime_end) & + emis_hourly(:) & * cropfrac_data(i_inner,j_inner)*1000. & * area_by_grid(i_inner,j_inner,imap) EMIS_PM10_met(i_outer,j_outer,itime_start:itime_end) = & EMIS_PM10_met(i_outer,j_outer,itime_start:itime_end) & + emis_hourly(:) & * cropfrac_data(i_inner,j_inner)*1000. & * area_by_grid(i_inner,j_inner,imap) end do day_loop end do imap_loop end do j_loop_ero end do i_loop_ero end do j_loop_met end do i_loop_met ! -------------------------------------------------------------------- ! Output dust emissions in Models-3 I/O API file format ! -------------------------------------------------------------------- ! ( File was open/created earlier in the program ) writeDate = stdate writeTime = sttime do itime = 1, ntime ! Write out to files if ( .not. ( write3('OUTFILE1', 'DUSTPM10',writeDate,writeTime,EMIS_PM10_met(:,:,itime ) ) ) ) then mesg = "Could not write to OUTFILE1 variable DUSTPM10" call m3exit( trim(progname), 0, 0, mesg, 2 ) end if ! Write out to file if ( .not. ( write3('OUTFILE2', 'DUSTPM10',writeDate,writeTime,EMIS_PM10_ero(:,:,itime ) ) ) ) then mesg = "Could not write to OUTFILE2 variable DUSTPM10" 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_nest_v7a