diff --git a/examples/acoustics_2d_radial/Makefile_abl b/examples/acoustics_2d_radial/Makefile_abl
new file mode 100644
index 0000000..e821b7c
--- /dev/null
+++ b/examples/acoustics_2d_radial/Makefile_abl
@@ -0,0 +1,60 @@
+
+# Makefile for Clawpack code in this directory.
+# This version only sets the local files and frequently changed
+# options, and then includes the standard makefile pointed to by CLAWMAKE.
+CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common
+
+# See the above file for details and a list of make options, or type
+#   make .help
+# at the unix prompt.
+
+
+# Adjust these variables if desired:
+# ----------------------------------
+
+CLAW_PKG = classic                  # Clawpack package to use
+EXE = xclaw                         # Executable to create
+SETRUN_FILE = setrun_abl.py         # File containing function to make data
+OUTDIR = _output_abl                # Directory for output
+SETPLOT_FILE = setplot_abl.py       # File containing function to set plots
+PLOTDIR = _plots_abl                # Directory for plots
+
+OVERWRITE ?= True                   # False ==> make a copy of OUTDIR first
+
+# Environment variable FC should be set to fortran compiler, e.g. gfortran
+
+# Compiler flags can be specified here or set as an environment variable
+FFLAGS ?=
+
+# ---------------------------------
+# package sources for this program:
+# ---------------------------------
+
+include $(CLAW)/classic/src/2d/Makefile.classic_2d
+
+# ---------------------------------------
+# package sources specifically to exclude
+# (i.e. if a custom replacement source
+#  under a different name is provided)
+# ---------------------------------------
+
+EXCLUDE_MODULES = \
+
+EXCLUDE_SOURCES = \
+
+
+# ---------------------------------
+# List of custom sources for this program:
+# ---------------------------------
+
+MODULES = \
+
+SOURCES = \
+  qinit.f \
+  setprob.f \
+  $(CLAW)/riemann/src/rpn2_acoustics.f90 \
+  $(CLAW)/riemann/src/rpt2_acoustics.f90
+
+#-------------------------------------------------------------------
+# Include Makefile containing standard definitions and make options:
+include $(CLAWMAKE)
diff --git a/examples/acoustics_2d_radial/README.rst b/examples/acoustics_2d_radial/README.rst
index ff8878a..c5cc8b2 100644
--- a/examples/acoustics_2d_radial/README.rst
+++ b/examples/acoustics_2d_radial/README.rst
@@ -10,3 +10,21 @@ compute the "true solution" and then setplot.py contains code to produce a
 scatter plot of the computed 2d pressure vs. distance from origin and compare
 with the 1d solution.
 
+Extrapolation BCs
+------------------
+
+The code is set up to use extrapolation boundary conditions at all
+boundaries. This does a reasonably good job of providing non-reflecting 
+boundaries, but there are some artifacts visible at later times.
+
+Absorbing boundary layer
+------------------------
+
+New cababilities have been added to Clawpack 5.5.0 to allow extending the
+computational domain with an aborbing boundary layer that does a better job
+of eliminating artificial reflections.  [Add more discussion.]
+
+To try this version::
+
+    make all -f Makefile_abl
+
diff --git a/examples/acoustics_2d_radial/setplot_abl.py b/examples/acoustics_2d_radial/setplot_abl.py
new file mode 100644
index 0000000..18ceea0
--- /dev/null
+++ b/examples/acoustics_2d_radial/setplot_abl.py
@@ -0,0 +1,149 @@
+
+"""
+Set up the plot figures, axes, and items to be done for each frame.
+
+This module is imported by the plotting routines and then the
+function setplot is called to set the plot parameters.
+
+"""
+
+from __future__ import absolute_import
+from __future__ import print_function
+import os
+if os.path.exists('./1drad/_output'):
+    qref_dir = os.path.abspath('./1drad/_output')
+else:
+    qref_dir = None
+    print("Directory ./1drad/_output not found")
+
+
+#--------------------------
+def setplot(plotdata=None):
+#--------------------------
+
+    """
+    Specify what is to be plotted at each frame.
+    Input:  plotdata, an instance of clawpack.visclaw.data.ClawPlotData.
+    Output: a modified version of plotdata.
+
+    """
+
+    if plotdata is None:
+        from clawpack.visclaw.data import ClawPlotData
+        plotdata = ClawPlotData()
+
+
+    from clawpack.visclaw import colormaps
+    from clawpack.clawutil.data import ClawData
+
+    # determine original computational domain
+    abldata = ClawData()
+    abldata.read(plotdata.outdir + '/abl.data', force=True)
+    clawdata = ClawData()
+    clawdata.read(plotdata.outdir + '/claw.data', force=True)
+    x1 = clawdata.lower[0] + abldata.depth_lower[0]
+    x2 = clawdata.upper[0] - abldata.depth_upper[0]
+    y1 = clawdata.lower[1] + abldata.depth_lower[1]
+    y2 = clawdata.upper[1] - abldata.depth_upper[1]
+
+    plotdata.clearfigures()  # clear any old figures,axes,items data
+
+
+    # Figure for pressure
+    # -------------------
+
+    plotfigure = plotdata.new_plotfigure(name='Pressure', figno=0)
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = 'auto'
+    plotaxes.ylimits = 'auto'
+    plotaxes.title = 'Pressure'
+    plotaxes.scaled = True      # so aspect ratio is 1
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
+    plotitem.plot_var = 0
+    plotitem.pcolor_cmap = colormaps.blue_yellow_red
+    plotitem.pcolor_cmin = -2.0
+    plotitem.pcolor_cmax = 2.0
+    plotitem.add_colorbar = True
+
+    def plot_original_domain(current_data):
+        from matplotlib.pyplot import gca, text
+        ax = gca()
+        x = [x1,x2,x2,x1,x1]
+        y = [y1,y1,y2,y2,y1]
+        ax.plot(x, y, '--k')
+        text(-0.6,1.05,'Absorbing Boundary Layer')
+
+    plotaxes.afteraxes = plot_original_domain
+
+    # Figure for scatter plot
+    # -----------------------
+
+    plotfigure = plotdata.new_plotfigure(name='scatter', figno=3)
+    plotfigure.show = (qref_dir is not None)  # don't plot if 1d solution is missing
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = [0,1.5]
+    plotaxes.ylimits = [-2.,4.]
+    plotaxes.title = 'Scatter plot'
+
+    # Set up for item on these axes: scatter of 2d data
+    plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data')
+
+    def p_vs_r(current_data):
+        # Return radius of each grid cell and p value in the cell
+        from numpy.ma import MaskedArray, masked_where
+        from pylab import sqrt
+        x = current_data.x
+        y = current_data.y
+        r = sqrt(x**2 + y**2)
+        r = masked_where(x < x1, r)
+        r = masked_where(x > x2, r)
+        r = masked_where(y < y1, r)
+        r = masked_where(y > y2, r)
+        q = current_data.q
+        p = MaskedArray(q[0,:,:], mask=r.mask)
+        return r,p
+
+    plotitem.map_2d_to_1d = p_vs_r
+    plotitem.plot_var = 0
+    plotitem.plotstyle = 'o'
+    plotitem.color = 'b'
+    plotitem.show = (qref_dir is not None)       # show on plot?
+
+    # Set up for item on these axes: 1d reference solution
+    plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
+    plotitem.outdir = qref_dir
+    plotitem.plot_var = 0
+    plotitem.plotstyle = '-'
+    plotitem.color = 'r'
+    plotitem.kwargs = {'linewidth': 2}
+    plotitem.show = True       # show on plot?
+
+    def make_legend(current_data):
+        import matplotlib.pyplot as plt
+        plt.legend(('2d data (interior only)', '1d reference solution'))
+
+    plotaxes.afteraxes = make_legend
+
+
+    # Parameters used only when creating html and/or latex hardcopy
+    # e.g., via clawpack.visclaw.frametools.printframes:
+
+    plotdata.printfigs = True                # print figures
+    plotdata.print_format = 'png'            # file format
+    plotdata.print_framenos = 'all'          # list of frames to print
+    plotdata.print_fignos = 'all'            # list of figures to print
+    plotdata.html = True                     # create html files of plots?
+    plotdata.html_homelink = '../README.html'   # pointer for top of index
+    plotdata.html_movie = 'JSAnimation'      # new style, or "4.x" for old style
+    plotdata.latex = True                    # create latex file of plots?
+    plotdata.latex_figsperline = 2           # layout of plots
+    plotdata.latex_framesperline = 1         # layout of plots
+    plotdata.latex_makepdf = False           # also run pdflatex?
+
+    return plotdata
diff --git a/examples/acoustics_2d_radial/setrun_abl.py b/examples/acoustics_2d_radial/setrun_abl.py
new file mode 100644
index 0000000..02568e4
--- /dev/null
+++ b/examples/acoustics_2d_radial/setrun_abl.py
@@ -0,0 +1,260 @@
+"""
+Module to set up run time parameters for Clawpack -- classic code.
+
+The values set in the function setrun are then written out to data files
+that will be read in by the Fortran code.
+
+"""
+
+from __future__ import absolute_import
+import os
+import numpy as np
+
+#------------------------------
+def setrun(claw_pkg='classic'):
+#------------------------------
+
+    """
+    Define the parameters used for running Clawpack.
+
+    INPUT:
+        claw_pkg expected to be "classic" for this setrun.
+
+    OUTPUT:
+        rundata - object of class ClawRunData
+
+    """
+
+    from clawpack.clawutil import data
+
+
+    assert claw_pkg.lower() == 'classic',  "Expected claw_pkg = 'classic'"
+
+    num_dim = 2
+    rundata = data.ClawRunData(claw_pkg, num_dim)
+
+
+    #------------------------------------------------------------------
+    # Problem-specific parameters to be written to setprob.data:
+    #------------------------------------------------------------------
+
+    probdata = rundata.new_UserData(name='probdata',fname='setprob.data')
+    probdata.add_param('rho',     1.,  'density of medium')
+    probdata.add_param('bulk',    4.,  'bulk modulus')
+
+    #------------------------------------------------------------------
+    # Standard Clawpack parameters to be written to claw.data:
+    #   (or to amrclaw.data for AMR)
+    #------------------------------------------------------------------
+
+    clawdata = rundata.clawdata  # initialized when rundata instantiated
+
+
+    # Set single grid parameters first.
+    # See below for AMR parameters.
+
+
+    # ---------------
+    # Spatial domain:
+    # ---------------
+
+    # Number of space dimensions:
+    clawdata.num_dim = num_dim
+
+    # Lower and upper edge of computational domain:
+    clawdata.lower[0] = -1.000000e+00          # xlower
+    clawdata.upper[0] = 1.000000e+00          # xupper
+    clawdata.lower[1] = -1.000000e+00          # ylower
+    clawdata.upper[1] = 1.000000e+00          # yupper
+
+    # Number of grid cells:
+    clawdata.num_cells[0] = 50      # mx
+    clawdata.num_cells[1] = 50      # my
+
+    # -------------------------
+    # Absorbing boundary layer:
+    # -------------------------
+    abldata = clawdata.abldata
+    # Select grid mapping for layer
+    #   'trigonometric', 'appelo_colonius', 'petersson_sjogreen', user defined #
+    abldata.abltype = 'trigonometric'
+    # Set depth of layer at each computation boundary
+    abldata.depth_lower[0] = 0.2
+    abldata.depth_upper[0] = 0.2
+    abldata.depth_lower[1] = 0.2
+    abldata.depth_upper[1] = 0.2
+    # Set additional parameters if needed
+    #   1) trigonometric - no additional parameters needed
+    #   2) appelo_colonius - epsilon, p, q (e.g. {'epsilon': 1e-4, 'p': 4, 'q': 4})
+    #   3) petterson_sjogreen - epsilon (e.g. {'epsilon': 1e-4})
+    abldata.parameters = {}
+
+    # ---------------
+    # Size of system:
+    # ---------------
+
+    # Number of equations in the system:
+    clawdata.num_eqn = 3
+
+    # Number of auxiliary variables in the aux array (initialized in setaux)
+    clawdata.num_aux = 2
+
+    # Index of aux array corresponding to capacity function, if there is one:
+    clawdata.capa_index = 0
+
+
+    # -------------
+    # Initial time:
+    # -------------
+
+    clawdata.t0 = 0.000000
+
+
+    # Restart from checkpoint file of a previous run?
+    # If restarting, t0 above should be from original run, and the
+    # restart_file 'fort.chkNNNNN' specified below should be in
+    # the OUTDIR indicated in Makefile.
+
+    clawdata.restart = False               # True to restart from prior results
+    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
+
+
+    # -------------
+    # Output times:
+    #--------------
+
+    # Specify at what times the results should be written to fort.q files.
+    # Note that the time integration stops after the final output time.
+
+    clawdata.output_style = 1
+
+    if clawdata.output_style==1:
+        # Output ntimes frames at equally spaced times up to tfinal:
+        # Can specify num_output_times = 0 for no output
+        clawdata.num_output_times = 20
+        clawdata.tfinal = 1.0
+        clawdata.output_t0 = True  # output at initial (or restart) time?
+
+    elif clawdata.output_style == 2:
+        # Specify a list or numpy array of output times:
+        # Include t0 if you want output at the initial time.
+        clawdata.output_times =  [0., 0.1]
+
+    elif clawdata.output_style == 3:
+        # Output every step_interval timesteps over total_steps timesteps:
+        clawdata.output_step_interval = 2
+        clawdata.total_steps = 4
+        clawdata.output_t0 = True  # output at initial (or restart) time?
+
+
+    clawdata.output_format = 'ascii'       # 'ascii', 'binary', 'netcdf'
+
+
+    # ---------------------------------------------------
+    # Verbosity of messages to screen during integration:
+    # ---------------------------------------------------
+
+    # The current t, dt, and cfl will be printed every time step
+    # at AMR levels <= verbosity.  Set verbosity = 0 for no printing.
+    #   (E.g. verbosity == 2 means print only on levels 1 and 2.)
+    clawdata.verbosity = 0
+
+
+
+    # --------------
+    # Time stepping:
+    # --------------
+
+    # if dt_variable==True:  variable time steps used based on cfl_desired,
+    # if dt_variable==False: fixed time steps dt = dt_initial always used.
+    clawdata.dt_variable = True
+
+    # Initial time step for variable dt.
+    # (If dt_variable==0 then dt=dt_initial for all steps)
+    clawdata.dt_initial = 1.00000e-02
+
+    # Max time step to be allowed if variable dt used:
+    clawdata.dt_max = 1.000000e+99
+
+    # Desired Courant number if variable dt used
+    clawdata.cfl_desired = 0.900000
+    # max Courant number to allow without retaking step with a smaller dt:
+    clawdata.cfl_max = 1.000000
+
+    # Maximum number of time steps to allow between output times:
+    clawdata.steps_max = 50000
+
+
+    # ------------------
+    # Method to be used:
+    # ------------------
+
+    # Order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
+    clawdata.order = 2
+
+    # Use dimensional splitting? (not yet available for AMR)
+    clawdata.dimensional_split = 'unsplit'
+
+    # For unsplit method, transverse_waves can be
+    #  0 or 'none'      ==> donor cell (only normal solver used)
+    #  1 or 'increment' ==> corner transport of waves
+    #  2 or 'all'       ==> corner transport of 2nd order corrections too
+    clawdata.transverse_waves = 2
+
+
+    # Number of waves in the Riemann solution:
+    clawdata.num_waves = 2
+
+    # List of limiters to use for each wave family:
+    # Required:  len(limiter) == num_waves
+    # Some options:
+    #   0 or 'none'     ==> no limiter (Lax-Wendroff)
+    #   1 or 'minmod'   ==> minmod
+    #   2 or 'superbee' ==> superbee
+    #   3 or 'vanleer'  ==> van Leer
+    #   4 or 'mc'       ==> MC limiter
+    clawdata.limiter = ['vanleer','vanleer']
+
+    clawdata.use_fwaves = False    # True ==> use f-wave version of algorithms
+
+    # Source terms splitting:
+    #   src_split == 0 or 'none'    ==> no source term (src routine never
+    #   called)
+    #   src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
+    #   src_split == 2 or 'strang'  ==> Strang (2nd order) splitting used,
+    #   not recommended.
+    clawdata.source_split = 0
+
+
+    # --------------------
+    # Boundary conditions:
+    # --------------------
+
+    # Number of ghost cells (usually 2)
+    clawdata.num_ghost = 2
+
+    # Choice of BCs at xlower and xupper:
+    #   0 or 'user'     => user specified (must modify bcNamr.f to use this
+    #   option)
+    #   1 or 'extrap'   => extrapolation (non-reflecting outflow)
+    #   2 or 'periodic' => periodic (must specify this at both boundaries)
+    #   3 or 'wall'     => solid wall for systems where q(2) is normal
+    #   velocity
+
+    clawdata.bc_lower[0] = 'extrap'   # at xlower
+    clawdata.bc_upper[0] = 'extrap'   # at xupper
+
+    clawdata.bc_lower[1] = 'extrap'   # at ylower
+    clawdata.bc_upper[1] = 'extrap'   # at yupper
+
+    return rundata
+
+    # end of function setrun
+    # ----------------------
+
+
+if __name__ == '__main__':
+    # Set up run-time parameters and write all data files.
+    import sys
+    rundata = setrun(*sys.argv[1:])
+    rundata.write()
diff --git a/src/2d/Makefile.classic_2d b/src/2d/Makefile.classic_2d
index 37bc02b..c0be700 100644
--- a/src/2d/Makefile.classic_2d
+++ b/src/2d/Makefile.classic_2d
@@ -1,6 +1,11 @@
 #get the directory of this makefile
 LIB:=$(CLAW)/classic/src/2d
 
+#list of common modules files for classic 2d codes
+COMMON_MODULES = \
+  $(LIB)/abl_module.f90
+
+
 #list of common source files for classic 2d codes
 COMMON_SOURCES = \
   $(LIB)/qinit.f \
diff --git a/src/2d/abl_module.f90 b/src/2d/abl_module.f90
new file mode 100644
index 0000000..304de6f
--- /dev/null
+++ b/src/2d/abl_module.f90
@@ -0,0 +1,278 @@
+module abl_module
+
+  implicit none
+  integer, parameter :: NO_ABL = 0
+  integer, parameter :: TRIG = 1
+  integer, parameter :: AC = 2 ! Appelo & Colonius
+  integer, parameter :: PS = 3 ! Petersson & Sjogreen
+
+  integer, parameter :: MAX_PARAM = 3
+
+  real(kind=8), parameter :: PIH = 2.d0*datan(1.d0)
+
+
+  integer :: abl_type = 0
+  integer :: p_AC, q_AC
+  real(kind=8) :: eps_AC, eps_PS
+  real(kind=8) :: xdepth(2), ydepth(2), xpos(2), ypos(2)
+  logical :: initialized = .false.
+
+contains
+
+  subroutine initialize()
+
+    implicit none
+
+    real(kind=8) :: xlower, xupper, ylower, yupper
+
+    ! Read in ABL data
+    ! open the unit with new routine from Clawpack 4.4 to skip over
+    ! comment lines starting with #:
+    call opendatafile(7, 'abl.data')
+
+    ! Read in parameters:
+    read(7,*) abl_type
+    if (abl_type == TRIG) then
+      print *, "ABL is trigonometric"
+    else if (abl_type == AC) then
+      read(7,*) eps_AC
+      read(7,*) p_AC
+      read(7,*) q_AC
+      print *, "ABL is Appelo/Colonius with (eps,p,q)=", eps_AC, p_AC, q_AC
+    else if (abl_type == PS) then
+      read(7,*) eps_PS
+      print *, "ABL is Petersson/Sjogreen with eps=", eps_PS
+    end if
+    read(7,*) xdepth(1), ydepth(1)
+    read(7,*) xdepth(2), ydepth(2)
+    close(7)
+
+    ! Read in grid parameters
+    ! open the unit with new routine from Clawpack 4.4 to skip over
+    ! comment lines starting with #:
+    call opendatafile(7, 'claw.data')
+
+    read(7,*) xlower ! num_dim
+    read(7,*) xlower, ylower ! lower
+    read(7,*) xupper, yupper ! upper
+    close(7)
+
+    ! Compute ABL positions
+    xpos(1) = xlower + xdepth(1)
+    xpos(2) = xupper - xdepth(2)
+    ypos(1) = ylower + ydepth(1)
+    ypos(2) = yupper - ydepth(2)
+
+    initialized = .true.
+
+  end subroutine initialize
+
+  subroutine set_scaling_factor(xlower, ylower, dx, dy, mx, my, num_ghost, &
+                                ablX_M, ablX_J, ablX_R, ablY_M, ablY_J, ablY_R, &
+                                i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper)
+
+    implicit none
+
+    integer, intent(in) :: mx, my, num_ghost
+    real(kind=8), intent(in) :: xlower, ylower, dx, dy
+    integer, intent(out) :: i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper
+    real(kind=8), intent(out) :: ablX_M(1-num_ghost:mx+num_ghost)
+    real(kind=8), intent(out) :: ablX_J(1-num_ghost:mx+num_ghost)
+    real(kind=8), intent(out) :: ablX_R(1-num_ghost:mx+num_ghost)
+    real(kind=8), intent(out) :: ablY_M(1-num_ghost:my+num_ghost)
+    real(kind=8), intent(out) :: ablY_J(1-num_ghost:my+num_ghost)
+    real(kind=8), intent(out) :: ablY_R(1-num_ghost:my+num_ghost)
+
+    integer :: i, j
+    real(kind=8) :: xcenter, ycenter, xedge, yedge
+
+    ! initialize module if needed
+    if (.not. initialized) then
+      call initialize()
+    end if
+
+    ! x direction
+
+    ! determine which indices are in the layer and loop over those
+    i_abl_lower = int((1.d0 + 1.d-14)*(xpos(1)-xlower)/dx)
+    i_abl_upper = 1 + int((1.d0 + 1.d-14)*(xpos(2)-xlower)/dx)
+
+    do i=1-num_ghost,i_abl_lower
+
+      xcenter = xlower + (i-0.5d0)*dx
+
+      ! Calculate inverse of Jacobian at cell edges
+      xedge = max(xcenter-0.5d0*dx, xpos(1)-xdepth(1))
+      xedge = min(xedge, xpos(2)+xdepth(2))
+      ablX_J(i) = 1.d0/(1.d0 + gp(xedge,1))
+
+      ! Calculate grid map scaling and ratio at cell centers
+      xcenter = max(xcenter, xpos(1)-xdepth(1)+0.5d0*dx)
+      xcenter = min(xcenter, xpos(2)+xdepth(2)-0.5d0*dx)
+      if (abl_type == TRIG) then
+        xedge = min(xedge, xpos(2)+xdepth(2)-dx)
+        ablX_M(i) = 1.d0/(1.d0 + (g(xedge+dx,1)-g(xedge,1))/dx)
+        ablX_R(i) = ablX_M(i)*(1.d0 + gp(xcenter,1))
+      else
+        ablX_M(i) = 1.d0/(1.d0 + gp(xcenter,1))
+        ablX_R(i) = 1.d0
+      end if
+
+    end do
+
+    do i=i_abl_upper,mx+num_ghost
+
+      xcenter = xlower + (i-0.5d0)*dx
+
+      ! Calculate inverse of Jacobian at cell edges
+      xedge = max(xcenter-0.5d0*dx, xpos(1)-xdepth(1))
+      xedge = min(xedge, xpos(2)+xdepth(2))
+      ablX_J(i) = 1.d0/(1.d0 + gp(xedge,1))
+
+      ! Calculate grid map scaling and ratio at cell centers
+      xcenter = max(xcenter, xpos(1)-xdepth(1)+0.5d0*dx)
+      xcenter = min(xcenter, xpos(2)+xdepth(2)-0.5d0*dx)
+      if (abl_type == TRIG) then
+        xedge = min(xedge, xpos(2)+xdepth(2)-dx)
+        ablX_M(i) = 1.d0/(1.d0 + (g(xedge+dx,1)-g(xedge,1))/dx)
+        ablX_R(i) = ablX_M(i)*(1.d0 + gp(xcenter,1))
+      else
+        ablX_M(i) = 1.d0/(1.d0 + gp(xcenter,1))
+        ablX_R(i) = 1.d0
+      end if
+
+    end do
+
+
+    ! y direction
+
+    ! determine which indices are in the layer and loop over those
+    j_abl_lower = int((1.d0 + 1.d-14)*(ypos(1)-ylower)/dy)
+    j_abl_upper = 1 + int((1.d0 + 1.d-14)*(ypos(2)-ylower)/dy)
+
+    do j=1-num_ghost,j_abl_lower
+
+      ycenter = ylower + (j-0.5d0)*dy
+
+      ! Calculate inverse of Jacobian at cell edges
+      yedge = max(ycenter-0.5d0*dy, ypos(1)-ydepth(1))
+      yedge = min(yedge, ypos(2)+ydepth(2))
+      ablY_J(j) = 1.d0/(1.d0 + gp(yedge,2))
+
+      ! Calculate grid map scaling and ratio at cell centers
+      ycenter = max(ycenter, ypos(1)-ydepth(1)+0.5d0*dy)
+      ycenter = min(ycenter, ypos(2)+ydepth(2)-0.5d0*dy)
+      if (abl_type == TRIG) then
+        yedge = min(yedge, ypos(2)+ydepth(2)-dy)
+        ablY_M(j) = 1.d0/(1.d0 + (g(yedge+dy,2)-g(yedge,2))/dy)
+        ablY_R(j) = ablY_M(j)*(1.d0 + gp(ycenter,2))
+      else
+        ablY_M(j) = 1.d0/(1.d0 + gp(ycenter,2))
+        ablY_R(j) = 1.d0
+      end if
+
+    end do
+
+    do j=j_abl_upper,my+num_ghost
+
+      ycenter = ylower + (j-0.5d0)*dy
+
+      ! Calculate inverse of Jacobian at cell edges
+      yedge = max(ycenter-0.5d0*dy, ypos(1)-ydepth(1))
+      yedge = min(yedge, ypos(2)+ydepth(2))
+      ablY_J(j) = 1.d0/(1.d0 + gp(yedge,2))
+
+      ! Calculate grid map scaling and ratio at cell centers
+      ycenter = max(ycenter, ypos(1)-ydepth(1)+0.5d0*dy)
+      ycenter = min(ycenter, ypos(2)+ydepth(2)-0.5d0*dy)
+      if (abl_type == TRIG) then
+        yedge = min(yedge, ypos(2)+ydepth(2)-dy)
+        ablY_M(j) = 1.d0/(1.d0 + (g(yedge+dy,2)-g(yedge,2))/dy)
+        ablY_R(j) = ablY_M(j)*(1.d0 + gp(ycenter,2))
+      else
+        ablY_M(j) = 1.d0/(1.d0 + gp(ycenter,2))
+        ablY_R(j) = 1.d0
+      end if
+
+    end do
+
+
+  end subroutine set_scaling_factor
+
+  function gp(z,direction)
+
+    implicit none
+
+    real(kind=8) :: gp, z
+    integer :: direction
+
+    real(kind=8) :: pos(2), depth(2), sig, val
+
+    gp = 0.d0
+    sig = 0.d0
+    if (direction == 1) then
+      pos = xpos
+      depth = xdepth
+    else if (direction == 2) then
+      pos = ypos
+      depth = ydepth
+    end if
+
+    if (z < pos(1)) then
+      val = (pos(1)-z)/depth(1)
+    else if (z > pos(2)) then
+      val = (z-pos(2))/depth(2)
+    else
+      return
+    end if
+
+    if (abl_type == TRIG) then
+      gp = dtan(PIH*val)**2
+    else if (abl_type == AC) then
+      sig = (1.0 - eps_AC)*(1.0 - (1.0 - val)**p_AC)**q_AC
+      gp = sig/(1.0-sig)
+    else if (abl_type == PS) then
+      sig = (1.0 - eps_PS)* &
+      val**6*(462.0 - 1980.0*val + 3465.0*val**2 - 3080.0*val**3 + 1386.0*val**4 - 252.0*val**5)
+      gp = sig/(1.0-sig)
+    end if
+
+  end function gp
+
+  function g(z,direction)
+
+    implicit none
+
+    real(kind=8) :: g, z
+    integer :: direction
+
+    real(kind=8) :: pos(2), depth(2), val, depth_scalar
+
+    g = 0.d0
+    if (direction == 1) then
+      pos = xpos
+      depth = xdepth
+    else if (direction == 2) then
+      pos = ypos
+      depth = ydepth
+    end if
+
+    if (z < pos(1)) then
+      val = (z-pos(1))/depth(1)
+      depth_scalar = depth(1)
+    else if (z > pos(2)) then
+      val = (z-pos(2))/depth(2)
+      depth_scalar = depth(2)
+    else
+      return
+    end if
+
+    if (abl_type == TRIG) then
+      g = depth_scalar*(dtan(PIH*val)/PIH - val)
+    end if
+
+  end function g
+
+
+
+end module abl_module
diff --git a/src/2d/claw2.f b/src/2d/claw2.f
index 3de21fa..32debca 100644
--- a/src/2d/claw2.f
+++ b/src/2d/claw2.f
@@ -661,7 +661,7 @@ subroutine claw2(meqn,mwaves,maux,mbc,mx,my,
 c             
              call step2(maxm,meqn,mwaves,maux,mbc,mx,my,
      &                  work(i0qwrk1),q,aux,
-     &                  dx,dy,dt,method,mthlim,cfl,
+     &                  xlower,ylower,dx,dy,dt,method,mthlim,cfl,
      &                  work(i0qadd),work(i0fadd),work(i0gadd),
      &                  work(i0q1d),work(i0dtdx1),work(i0dtdy1),
      &                  work(i0aux1),work(i0aux2),work(i0aux3),
diff --git a/src/2d/claw2ez.f b/src/2d/claw2ez.f
index 78441bf..9afc68a 100644
--- a/src/2d/claw2ez.f
+++ b/src/2d/claw2ez.f
@@ -11,6 +11,8 @@ subroutine claw2ez    ! No arguments
 c
 c     Authors: Randall J. LeVeque, Grady I. Lemoine
 c
+      use abl_module, only: initialize_abl=>initialize
+ 
       implicit double precision (a-h,o-z)
       external bc2,rpn2,rpt2,src2,b4step2
 
@@ -233,6 +235,9 @@ subroutine claw2ez    ! No arguments
          go to 900
       end if
 
+c     # initialize abl module
+      call initialize_abl()
+
 c
 c     # set initial conditions:
 c
diff --git a/src/2d/flux2.f90 b/src/2d/flux2.f90
index a805e41..10c6a6a 100644
--- a/src/2d/flux2.f90
+++ b/src/2d/flux2.f90
@@ -5,7 +5,7 @@
     subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
     q1d,dtdx1d,aux1,aux2,aux3,method,mthlim, &
     qadd,fadd,gadd,cfl1d,wave,s, &
-    amdq,apdq,cqxx,bmasdq,bpasdq,rpn2,rpt2,use_fwave)
+    amdq,apdq,cqxx,bmasdq,bpasdq,rpn2,rpt2,use_fwave,abl,i_abl_lower,i_abl_upper)
 !     =====================================================
 
 !     # Compute the modification to fluxes f and g that are generated by
@@ -57,6 +57,8 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
 !f2py external rpn2, rpt2
 !f2py intent(callback) rpn2,rpt2
 
+    use abl_module, only: abl_type
+
     implicit double precision (a-h,o-z)
     integer num_aux
     external rpn2,rpt2
@@ -78,13 +80,16 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
     dimension     s(num_waves,1-num_ghost:maxm+num_ghost)
     dimension  wave(num_eqn, num_waves, 1-num_ghost:maxm+num_ghost)
 
+    dimension abl(1-num_ghost:maxm+num_ghost)
+    dimension dtdxave(1-num_ghost:mx+num_ghost)
+
     dimension method(7),mthlim(num_waves)
     logical :: limit, use_fwave
     common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
 
     limit = .false.
     do mw=1,num_waves
-        if (mthlim(mw) > 0) limit = .TRUE. 
+        if (mthlim(mw) > 0) limit = .TRUE.
     end do
 
 !     # initialize flux increments:
@@ -98,6 +103,22 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
         end do
     end do
 
+    ! Store dtdxave for second order corrections so that it may be scaled for
+    ! ABL if needed
+    do i = 2-num_ghost, mx+num_ghost
+      dtdxave(i) = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
+    end do
+
+    ! apply ABL if needed
+    if (abl_type /= 0) then
+        do i = 2-num_ghost, i_abl_lower
+            dtdxave(i) = dtdxave(i) * abl(i)
+        end do
+        do i = i_abl_upper, mx+num_ghost
+            dtdxave(i) = dtdxave(i) * abl(i)
+        end do
+    end if
+
     ! Unnecessary to init fadd to zero here because it only gets
     ! modified by the second-order corrections
 
@@ -113,10 +134,10 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
 !     # Set qadd for the donor-cell upwind method (Godunov)
     do i = 1, mx+1
         do m = 1, num_eqn    ! qadd(:,i-1) is still in cache from last cycle of outer loop
-            qadd(m,i-1) = qadd(m,i-1) - dtdx1d(i-1)*amdq(m,i)
+            qadd(m,i-1) = qadd(m,i-1) - amdq(m,i)
         end do
         do m = 1, num_eqn
-            qadd(m,i) = qadd(m,i) - dtdx1d(i)*apdq(m,i)
+            qadd(m,i) = qadd(m,i) - apdq(m,i)
         end do
     end do
 
@@ -147,14 +168,13 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
     if (use_fwave.eqv. .FALSE. ) then
         do i = 2-num_ghost, mx+num_ghost
         !            # modified in Version 4.3 to use average only in cqxx, not transverse
-            dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
             do m = 1, num_eqn
                 cqxx(m,i) = 0.d0
             end do
             do mw=1,num_waves    ! Traverse the wave array in memory-contiguous fashion
                 do m=1,num_eqn
                     cqxx(m,i) = cqxx(m,i) + dabs(s(mw,i)) &
-                    * (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i)
+                    * (1.d0 - dabs(s(mw,i))*dtdxave(i)) * wave(m,mw,i)
                 enddo
             enddo
             do m = 1, num_eqn
@@ -162,8 +182,11 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
             end do
         enddo
     else
+        if (abl_type /= 0) then
+            print *, "ABL not implemented yet for F-wave splitting"
+            stop
+        end if
         do i = 2-num_ghost, mx+num_ghost
-            dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
             do m = 1, num_eqn
                 cqxx(m,i) = 0.d0
             end do
@@ -171,7 +194,7 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
                 do m=1,num_eqn
                 !                 # second order corrections:
                     cqxx(m,i) = cqxx(m,i) + dsign(1.d0,s(mw,i)) &
-                    * (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i)
+                    * (1.d0 - dabs(s(mw,i))*dtdxave(i)) * wave(m,mw,i)
                 enddo
             enddo
             do m = 1, num_eqn
@@ -214,7 +237,7 @@ subroutine flux2(ixy,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
             gadd(m,2,i-1) = -0.5d0*dtdx1d(i-1) * bpasdq(m,i)
         end do
     end do
-     
+
 
 !     # split the right-going flux difference into down-going and up-going:
     call rpt2(ixy,2,maxm,num_eqn,num_waves,num_aux,num_ghost,mx,q1d,q1d, &
diff --git a/src/2d/step2.f90 b/src/2d/step2.f90
index c4cd3e1..ce44ca9 100644
--- a/src/2d/step2.f90
+++ b/src/2d/step2.f90
@@ -1,6 +1,6 @@
 !     ==========================================================
     subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, &
-    qold,qnew,aux,dx,dy,dt,method,mthlim,cfl, &
+    qold,qnew,aux,xlower,ylower,dx,dy,dt,method,mthlim,cfl, &
     qadd,fadd,gadd,q1d,dtdx1d,dtdy1d, &
     aux1,aux2,aux3,work,mwork,use_fwave,rpn2,rpt2)
 !     ==========================================================
@@ -15,6 +15,7 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, &
 !     # fadd and gadd are used to return flux increments from flux2.
 !     # See the flux2 documentation for more information.
 
+    use abl_module, only: abl_type, set_scaling_factor
 
     implicit double precision (a-h,o-z)
     external :: rpn2,rpt2
@@ -38,6 +39,16 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, &
     logical ::          use_fwave
     double precision :: work(mwork)
 
+    ! added for abl implementation
+    integer :: i_abl_lower, i_abl_upper, j_abl_lower, j_abl_upper
+    double precision :: xlower, ylower, dtdxm, dtdxc, dtdxp, dtdyp, dtdyc, dtdym
+    double precision :: ablX_M(1-num_ghost:mx+num_ghost)
+    double precision :: ablX_J(1-num_ghost:mx+num_ghost)
+    double precision :: ablX_R(1-num_ghost:mx+num_ghost)
+    double precision :: ablY_M(1-num_ghost:my+num_ghost)
+    double precision :: ablY_J(1-num_ghost:my+num_ghost)
+    double precision :: ablY_R(1-num_ghost:my+num_ghost)
+
     common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
 
 !f2py intent(out) cfl
@@ -76,6 +87,14 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, &
     cfl = 0.d0
     dtdx = dt/dx
     dtdy = dt/dy
+    dtdxm = dtdx
+    dtdxc = dtdx
+    dtdxp = dtdx
+    dtdym = dtdy
+    dtdyc = dtdy
+    dtdyp = dtdy
+
+
 
     if (index_capa == 0) then
     !        # no capa array:
@@ -85,25 +104,36 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, &
         end do
     endif
 
+!   compute abl scaling factors if needed
+    if (abl_type /= 0) then
+        call set_scaling_factor(xlower,ylower,dx,dy,mx,my,num_ghost, &
+                                ablX_M,ablX_J,ablX_R,ablY_M,ablY_J,ablY_R, &
+                                i_abl_lower,i_abl_upper,j_abl_lower,j_abl_upper)
+    else
+        i_abl_lower = 0
+        j_abl_lower = 0
+        i_abl_upper = mx+1
+        j_abl_upper = my+1
+    end if
 
 !     # perform x-sweeps
 !     ==================
 
     do j = 0,my+1
-    
+
     !        # copy data along a slice into 1d arrays:
         do i = 1-num_ghost, mx+num_ghost
             do m=1,num_eqn
                 q1d(m,i) = qold(m,i,j)
             end do
         end do
-    
+
         if (index_capa > 0)  then
             do i = 1-num_ghost, mx+num_ghost
                 dtdx1d(i) = dtdx/aux(index_capa,i,j)
             end do
         endif
-    
+
         if (num_aux > 0)  then
             do i = 1-num_ghost, mx+num_ghost
                 do ma=1,num_aux
@@ -113,45 +143,88 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, &
                 end do
             end do
         endif
-    
+
     !     # Store the value of j along this slice in the common block
     !        # comxyt in case it is needed in the Riemann solver (for
     !        # variable coefficient problems)
         jcom = j
-    
+
+        ! adjust dtdy* if ABL is present
+        if (abl_type /= 0) then
+            if (j-1 <= j_abl_lower .or. j-1 >= j_abl_upper) then
+                dtdym = dtdy * ablY_R(j-1) * ablY_M(j-1)
+            else
+                dtdym = dtdy
+            end if
+            if (j <= j_abl_lower .or. j >= j_abl_upper) then
+                dtdyc = dtdy * ablY_R(j) * ablY_M(j)
+                dtdxc = dtdx * ablY_R(j)
+            else
+                dtdyc = dtdy
+                dtdxc = dtdx
+            end if
+            if (j+1 <= j_abl_lower .or. j+1 >= j_abl_upper) then
+                dtdyp = dtdy * ablY_R(j+1) * ablY_M(j+1)
+            else
+                dtdyp = dtdy
+            end if
+        end if
+
     !        # compute modifications fadd and gadd to fluxes along this slice:
         call flux2(1,maxm,num_eqn,num_waves,num_aux,num_ghost,mx, &
         q1d,dtdx1d,aux1,aux2,aux3,method,mthlim, &
         qadd,fadd,gadd,cfl1d, &
         work(i0wave),work(i0s),work(i0amdq),work(i0apdq), &
         work(i0cqxx),work(i0bmadq),work(i0bpadq),rpn2,rpt2, &
-        use_fwave)
+        use_fwave,ablX_J,i_abl_lower,i_abl_upper)
         cfl = dmax1(cfl,cfl1d)
-    
+
     !        # update qnew by flux differencing.
     !        # (rather than maintaining arrays f and g for the total fluxes,
     !        # the modifications are used immediately to update qnew
     !        # in order to save storage.)
-    
+
         if (index_capa == 0) then
-        
+
         !            # no capa array.  Standard flux differencing:
-            do i=1,mx
+            do i=i_abl_lower+1,i_abl_upper-1
+                do m=1,num_eqn
+                    qnew(m,i,j) = qnew(m,i,j) + dtdxc * qadd(m,i) &
+                    - dtdxc * (fadd(m,i+1) - fadd(m,i)) &
+                    - dtdyc * (gadd(m,2,i) - gadd(m,1,i))
+                    qnew(m,i,j-1) = qnew(m,i,j-1) - dtdym * gadd(m,1,i)
+                    qnew(m,i,j+1) = qnew(m,i,j+1) + dtdyp * gadd(m,2,i)
+                end do
+            end do
+            ! these loops are skipped if no ABL present
+            do i=1,i_abl_lower
                 do m=1,num_eqn
-                    qnew(m,i,j) = qnew(m,i,j) + qadd(m,i) &
-                    - dtdx * (fadd(m,i+1) - fadd(m,i)) &
-                    - dtdy * (gadd(m,2,i) - gadd(m,1,i))
-                    qnew(m,i,j-1) = qnew(m,i,j-1) - dtdy * gadd(m,1,i)
-                    qnew(m,i,j+1) = qnew(m,i,j+1) + dtdy * gadd(m,2,i)
+                    qnew(m,i,j) = qnew(m,i,j) + ablX_M(i)*( dtdxc * qadd(m,i) &
+                    - dtdxc * (fadd(m,i+1) - fadd(m,i)) &
+                    - dtdyc * (gadd(m,2,i) - gadd(m,1,i)) )
+                    qnew(m,i,j-1) = qnew(m,i,j-1) - dtdym*ablX_M(i) * gadd(m,1,i)
+                    qnew(m,i,j+1) = qnew(m,i,j+1) + dtdyp*ablX_M(i) * gadd(m,2,i)
                 end do
             end do
-        
+            do i=i_abl_upper,mx
+                do m=1,num_eqn
+                    qnew(m,i,j) = qnew(m,i,j) + ablX_M(i)*( dtdxc * qadd(m,i) &
+                    - dtdxc * (fadd(m,i+1) - fadd(m,i)) &
+                    - dtdyc * (gadd(m,2,i) - gadd(m,1,i)) )
+                    qnew(m,i,j-1) = qnew(m,i,j-1) - dtdym*ablX_M(i) * gadd(m,1,i)
+                    qnew(m,i,j+1) = qnew(m,i,j+1) + dtdyp*ablX_M(i) * gadd(m,2,i)
+                end do
+            end do
+
         else
-        
         !            # with capa array.
+            if (abl_type /= 0) then
+                print *, "ABL not implemented yet for capa array"
+                stop
+            end if
             do i=1,mx
                 do m=1,num_eqn
-                    qnew(m,i,j) = qnew(m,i,j) + qadd(m,i) &
+                    qnew(m,i,j) = qnew(m,i,j) + dtdx * qadd(m,i) &
                     - (dtdx * (fadd(m,i+1) - fadd(m,i)) &
                     +  dtdy * (gadd(m,2,i) - gadd(m,1,i))) &
                     / aux(index_capa,i,j)
@@ -169,22 +242,21 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, &
 !     # perform y sweeps
 !     ==================
 
-
     do i = 0, mx+1
-    
+
     !        # copy data along a slice into 1d arrays:
         do j = 1-num_ghost, my+num_ghost
             do m=1,num_eqn
                 q1d(m,j) = qold(m,i,j)
             end do
         end do
-    
+
         if (index_capa > 0)  then
             do j = 1-num_ghost, my+num_ghost
                 dtdy1d(j) = dtdy / aux(index_capa,i,j)
             end do
         endif
-    
+
         if (num_aux > 0)  then
             do j = 1-num_ghost, my+num_ghost
                 do ma=1,num_aux
@@ -194,47 +266,91 @@ subroutine step2(maxm,num_eqn,num_waves,num_aux,num_ghost,mx,my, &
                 end do
             end do
         endif
-    
-    
+
     !     # Store the value of i along this slice in the common block
     !        # comxyt in case it is needed in the Riemann solver (for
     !        # variable coefficient problems)
         icom = i
-    
+
+        ! adjust dtdx* if ABL is present
+        if (abl_type /= 0) then
+            if (i-1 <= i_abl_lower .or. i-1 >= i_abl_upper) then
+                dtdxm = dtdx * ablX_R(i-1) * ablX_M(i-1)
+            else
+                dtdxm = dtdx
+            end if
+            if (i <= i_abl_lower .or. i >= i_abl_upper) then
+                dtdxc = dtdx * ablX_R(i) * ablX_M(i)
+                dtdyc = dtdy * ablX_R(i)
+            else
+                dtdxc = dtdx
+                dtdyc = dtdy
+            end if
+            if (i+1 <= i_abl_lower .or. i+1 >= i_abl_upper) then
+                dtdxp = dtdx * ablX_R(i+1) * ablX_M(i+1)
+            else
+                dtdxp = dtdx
+            end if
+        end if
+
+
     !        # compute modifications fadd and gadd to fluxes along this slice:
         call flux2(2,maxm,num_eqn,num_waves,num_aux,num_ghost,my, &
         q1d,dtdy1d,aux1,aux2,aux3,method,mthlim, &
         qadd,fadd,gadd,cfl1d, &
         work(i0wave),work(i0s),work(i0amdq),work(i0apdq), &
         work(i0cqxx),work(i0bmadq),work(i0bpadq),rpn2,rpt2, &
-        use_fwave)
-    
+        use_fwave,ablY_J,j_abl_lower,j_abl_upper)
+
         cfl = dmax1(cfl,cfl1d)
-    
+
     !        # update qnew by flux differencing.
     !        # Note that the roles of fadd and gadd are reversed for
     !        # the y-sweeps -- fadd is the modification to g-fluxes and
     !        # gadd is the modification to f-fluxes to the left and right.
-    
+
         if (index_capa == 0) then
-        
+
         !            # no capa array.  Standard flux differencing:
-            do j=1,my
+            do j=j_abl_lower+1,j_abl_upper-1
+                do m=1,num_eqn
+                    qnew(m,i,j) = qnew(m,i,j) + dtdyc * qadd(m,j) &
+                    - dtdyc * (fadd(m,j+1) - fadd(m,j)) &
+                    - dtdxc * (gadd(m,2,j) - gadd(m,1,j))
+                    qnew(m,i-1,j) = qnew(m,i-1,j) - dtdxm * gadd(m,1,j)
+                    qnew(m,i+1,j) = qnew(m,i+1,j) + dtdxp * gadd(m,2,j)
+                end do
+            end do
+            ! these loops are skipped if no ABL present
+            do j=1,j_abl_lower
                 do m=1,num_eqn
-                    qnew(m,i,j) = qnew(m,i,j) + (qadd(m,j) &
-                    - dtdy * (fadd(m,j+1) - fadd(m,j)) &
-                    - dtdx * (gadd(m,2,j) - gadd(m,1,j)))
-                    qnew(m,i-1,j) = qnew(m,i-1,j) - dtdx * gadd(m,1,j)
-                    qnew(m,i+1,j) = qnew(m,i+1,j) + dtdx * gadd(m,2,j)
+                    qnew(m,i,j) = qnew(m,i,j) + ablY_M(j)*( dtdyc * qadd(m,j) &
+                    - dtdyc * (fadd(m,j+1) - fadd(m,j)) &
+                    - dtdxc * (gadd(m,2,j) - gadd(m,1,j)) )
+                    qnew(m,i-1,j) = qnew(m,i-1,j) - dtdxm*ablY_M(j) * gadd(m,1,j)
+                    qnew(m,i+1,j) = qnew(m,i+1,j) + dtdxp*ablY_M(j) * gadd(m,2,j)
                 end do
             end do
-        
+            do j=j_abl_upper,my
+                do m=1,num_eqn
+                    qnew(m,i,j) = qnew(m,i,j) + ablY_M(j)*( dtdyc * qadd(m,j) &
+                    - dtdyc * (fadd(m,j+1) - fadd(m,j)) &
+                    - dtdxc * (gadd(m,2,j) - gadd(m,1,j)) )
+                    qnew(m,i-1,j) = qnew(m,i-1,j) - dtdxm*ablY_M(j) * gadd(m,1,j)
+                    qnew(m,i+1,j) = qnew(m,i+1,j) + dtdxp*ablY_M(j) * gadd(m,2,j)
+                end do
+            end do
+
         else
-        
+
         !            # with capa array.
+            if (abl_type /= 0) then
+                print *, "ABL not implemented yet for capa array"
+                stop
+            end if
             do j=1,my
                 do m=1,num_eqn
-                    qnew(m,i,j) = qnew(m,i,j) + qadd(m,j) &
+                    qnew(m,i,j) = qnew(m,i,j) + dtdy * qadd(m,j) &
                     - (dtdy * (fadd(m,j+1) - fadd(m,j)) &
                     + dtdx * (gadd(m,2,j) - gadd(m,1,j))) &
                     / aux(index_capa,i,j)