diff --git a/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F b/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F
index d6dc1bc0c1..facd88dd0b 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F
@@ -52,6 +52,13 @@ module mpas_atmphys_initialize_real
 ! * In subroutine physics_init_seaice, removed the initialization of isice_lu since it is now defined in
 !   Registry.xml and initialized in subroutine init_atm_static.
 !   Laura D. Fowler (laura@ucar.edu) / 2017-01-12.
+! * generalized soilcat (isltyp) to work with either soiltype_top or soiltype_bot data.
+!   Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2018-09-10.
+! * added fractional soil type (soiltypf) and land use (landusef) dependent on config_frac_landuse option;
+!   otherwise dominant land use and soil type will be used by default.
+!   Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2020-02-10.
+! * added initialization of lai_modis from MODIS climatological monthly mean LAI for use by PX LSM.
+!   Jerold A. Herwehe (herwehe.jerry@epa.gov) / 2020-02-20.
 ! * In subroutine physics_init_seaice, added the initialization of the annual maximum snow albedo over seaice
 !   points to 0.75.
 ! * Laura D. Fowler (laura@ucar.edu) / 2022-03-15).
@@ -88,9 +95,9 @@ subroutine physics_initialize_real(mesh, fg, dminfo, dims, configs)
  real(kind=RKIND),dimension(:,:),pointer:: albedo12m
 
  real(kind=RKIND),dimension(:),pointer:: seaice,xice,xland
- real(kind=RKIND),dimension(:),pointer:: vegfra,shdmin,shdmax
+ real(kind=RKIND),dimension(:),pointer:: vegfra,shdmin,shdmax,lai_modis
  real(kind=RKIND),dimension(:),pointer:: snow,snowc,snowh
- real(kind=RKIND),dimension(:,:),pointer:: greenfrac
+ real(kind=RKIND),dimension(:,:),pointer:: greenfrac,lai12m
 
  real(kind=RKIND),dimension(:),pointer:: skintemp,sst
 
@@ -116,9 +123,11 @@ subroutine physics_initialize_real(mesh, fg, dminfo, dims, configs)
  call mpas_pool_get_array(mesh, 'greenfrac', greenfrac)
  call mpas_pool_get_array(mesh, 'shdmin', shdmin)
  call mpas_pool_get_array(mesh, 'shdmax', shdmax)
+ call mpas_pool_get_array(mesh, 'lai12m', lai12m)
 
  call mpas_pool_get_array(fg, 'sfc_albbck', sfc_albbck)
  call mpas_pool_get_array(fg, 'vegfra', vegfra)
+ call mpas_pool_get_array(fg, 'lai_modis', lai_modis)
  call mpas_pool_get_array(fg, 'snow', snow)
  call mpas_pool_get_array(fg, 'snowc', snowc)
  call mpas_pool_get_array(fg, 'snowh', snowh)
@@ -157,10 +166,12 @@ subroutine physics_initialize_real(mesh, fg, dminfo, dims, configs)
     if(landmask(iCell) .eq. 0) sfc_albbck(iCell) = 0.08_RKIND
  enddo
 
-!initialization of the green-ness (vegetation) fraction: interpolation of the monthly values to
-!the initial date. get the min/max for each cell for the monthly green-ness fraction:
+!initialization of the green-ness (vegetation) fraction and leaf area index: 
+!interpolation of the monthly values to the initial date, and 
+!get the min/max for each cell for the monthly green-ness fraction:
  initial_date = trim(config_start_time)
  call monthly_interp_to_date(nCellsSolve,initial_date,greenfrac,vegfra)
+ call monthly_interp_to_date(nCellsSolve,initial_date,lai12m,lai_modis)
 
 !calculates the maximum and minimum green-ness (vegetation) fraction:
  call monthly_min_max(nCellsSolve,greenfrac,shdmin,shdmax)
@@ -304,9 +315,6 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
 
  call mpas_pool_get_config(configs, 'config_nsoillevels', config_nsoillevels)
 
- if(config_nsoillevels .ne. 4) &
-    call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.')
-
  do iCell = 1, nCellsSolve
     iSoil = 1
     zs_fg(iSoil,iCell) = 0.5_RKIND * dzs_fg(iSoil,iCell)
@@ -317,21 +325,45 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
     enddo
  enddo
 
- do iCell = 1, nCellsSolve
-    dzs(1,iCell) = 0.10_RKIND
-    dzs(2,iCell) = 0.30_RKIND
-    dzs(3,iCell) = 0.60_RKIND
-    dzs(4,iCell) = 1.00_RKIND
+ if(config_nsoillevels .eq. 4) then  ! for Noah lsm
+
+    do iCell = 1, nCellsSolve
+       dzs(1,iCell) = 0.10_RKIND
+       dzs(2,iCell) = 0.30_RKIND
+       dzs(3,iCell) = 0.60_RKIND
+       dzs(4,iCell) = 1.00_RKIND
+
+       iSoil = 1
+       zs(iSoil,iCell)  = 0.5_RKIND * dzs(iSoil,iCell)
+       do iSoil = 2, nSoilLevels
+          zs(iSoil,iCell) = zs(iSoil-1,iCell)              &
+                          + 0.5_RKIND * dzs(iSoil-1,iCell) &
+                          + 0.5_RKIND * dzs(iSoil,iCell)
+       enddo
 
-    iSoil = 1
-    zs(iSoil,iCell)  = 0.5_RKIND * dzs(iSoil,iCell)
-    do iSoil = 2, nSoilLevels
-       zs(iSoil,iCell) = zs(iSoil-1,iCell)              &
-                       + 0.5_RKIND * dzs(iSoil-1,iCell) &
-                       + 0.5_RKIND * dzs(iSoil,iCell)
     enddo
 
- enddo
+ elseif(config_nsoillevels .eq. 2) then  ! for PX lsm
+
+    do iCell = 1, nCellsSolve
+       dzs(1,iCell) = 0.01_RKIND
+       dzs(2,iCell) = 0.99_RKIND
+
+       iSoil = 1
+       zs(iSoil,iCell)  = 0.5_RKIND * dzs(iSoil,iCell)
+       do iSoil = 2, nSoilLevels
+          zs(iSoil,iCell) = zs(iSoil-1,iCell)              &
+                          + 0.5_RKIND * dzs(iSoil-1,iCell) &
+                          + 0.5_RKIND * dzs(iSoil,iCell)
+       enddo
+
+    enddo
+
+ else
+
+    call physics_error_fatal('NOAH lsm uses 4 and PX lsm uses 2 soil layers. Correct config_nsoillevels.')
+
+ endif
 
  end subroutine init_soil_layers_depth
 
@@ -409,8 +441,8 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
     call mpas_log_write('Error in interpolation of sm_fg to MPAS grid: num_sm = $i', messageType=MPAS_LOG_CRIT, intArgs=(/num_sm/))
  endif 
 
- if(config_nsoillevels .ne. 4) &
-    call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.')
+ if(config_nsoillevels .ne. 4 .and. config_nsoillevels .ne. 2) &
+    call physics_error_fatal('NOAH lsm uses 4 and PX lsm uses 2 soil layers. Correct config_nsoillevels.')
 
  if(.not.allocated(zhave)   ) allocate(zhave(nFGSoilLevels+2,nCellsSolve)   )
  if(.not.allocated(st_input)) allocate(st_input(nFGSoilLevels+2,nCellsSolve))
@@ -441,14 +473,14 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
 
  enddo
 
-!... interpolate the soil temperature, soil moisture, and soil liquid temperature to the four
-!    layers used in the NOAH land surface scheme:
+!... interpolate the soil temperature, soil moisture, and soil liquid temperature to either
+!    the four layers used in the NOAH land surface scheme or the two layers used in the PX LSM:
 
  do iCell = 1, nCellsSolve
 
     if(landmask(iCell) .eq. 1) then
 
-       noah: do iSoil = 1 , nSoilLevels
+       do iSoil = 1 , nSoilLevels
           input: do ifgSoil = 1 , nFGSoilLevels+2-1
              if(iCell .eq. 1) call mpas_log_write('$i $i $r $r $r', &
                                                   intArgs=(/iSoil,ifgSoil/), &
@@ -477,7 +509,7 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
              endif
           enddo input
           if(iCell.eq. 1) call mpas_log_write('')
-       enddo noah
+       enddo
 
     elseif(landmask(iCell) .eq. 0) then
 
@@ -604,14 +636,16 @@ subroutine physics_init_seaice(mesh, input, dims, configs)
  real(kind=RKIND):: xice_threshold
  real(kind=RKIND):: mid_point_depth
  real(kind=RKIND),pointer:: tsk_seaice_threshold
- real(kind=RKIND),dimension(:),pointer  :: vegfra
+ real(kind=RKIND),dimension(:),pointer  :: vegfra,lai_modis
  real(kind=RKIND),dimension(:),pointer  :: seaice,snoalb,xice
  real(kind=RKIND),dimension(:),pointer  :: skintemp,tmn,xland
  real(kind=RKIND),dimension(:,:),pointer:: tslb,smois,sh2o,smcrel
+ real(kind=RKIND),dimension(:,:),pointer:: landusef,soiltypf
 
- logical, pointer :: config_frac_seaice
+ logical, pointer :: config_frac_seaice,config_frac_landuse
  character(len=StrKIND),pointer:: config_landuse_data
- integer,pointer:: isice_lu
+ integer,pointer:: isice_lu,ismax_lu
+ integer:: i
 
 !note that this threshold is also defined in module_physics_vars.F.It is defined here to avoid
 !adding "use module_physics_vars" since this subroutine is only used for the initialization of
@@ -625,19 +659,24 @@ subroutine physics_init_seaice(mesh, input, dims, configs)
  call mpas_pool_get_config(configs, 'config_frac_seaice'         , config_frac_seaice)
  call mpas_pool_get_config(configs, 'config_tsk_seaice_threshold', tsk_seaice_threshold)
  call mpas_pool_get_config(configs, 'config_landuse_data'        , config_landuse_data)
+ call mpas_pool_get_config(configs, 'config_frac_landuse'        , config_frac_landuse)
 
  call mpas_pool_get_dimension(dims, 'nCellsSolve', nCellsSolve)
  call mpas_pool_get_dimension(dims, 'nSoilLevels', nSoilLevels)
 
  call mpas_pool_get_array(mesh, 'isice_lu', isice_lu)
+ call mpas_pool_get_array(mesh, 'ismax_lu', ismax_lu)
  call mpas_pool_get_array(mesh, 'landmask', landmask)
  call mpas_pool_get_array(mesh, 'lu_index', ivgtyp)
- call mpas_pool_get_array(mesh, 'soilcat_top', isltyp)
  call mpas_pool_get_array(mesh, 'snoalb', snoalb)
+ call mpas_pool_get_array(mesh, 'landusef', landusef)
+ call mpas_pool_get_array(mesh, 'soilcat', isltyp)
+ call mpas_pool_get_array(mesh, 'soiltypf', soiltypf)
 
  call mpas_pool_get_array(input, 'seaice', seaice)
  call mpas_pool_get_array(input, 'xice', xice)
  call mpas_pool_get_array(input, 'vegfra', vegfra)
+ call mpas_pool_get_array(input, 'lai_modis', lai_modis)
 
  call mpas_pool_get_array(input, 'skintemp', skintemp)
  call mpas_pool_get_array(input, 'tmn', tmn)
@@ -649,6 +688,7 @@ subroutine physics_init_seaice(mesh, input, dims, configs)
  call mpas_pool_get_array(input, 'smcrel', smcrel)
 
  call mpas_log_write('--- isice_lu   = $i', intArgs=(/isice_lu/))
+ call mpas_log_write('--- ismax_lu   = $i', intArgs=(/ismax_lu/))
 
 !assign the threshold value for xice as a function of config_frac_seaice:
  if(.not. config_frac_seaice) then
@@ -678,10 +718,27 @@ subroutine physics_init_seaice(mesh, input, dims, configs)
        !... sea-ice points are converted to land points:
        if(landmask(iCell) .eq. 0) tmn(iCell) = 271.4_RKIND
        ivgtyp(iCell) = isice_lu
+!       if(.not. config_frac_landuse) then
+          landusef(isice_lu,iCell) = 1.0  ! defines dominant land use only
+          do i = 1,isice_lu-1
+             landusef(i,iCell) = 0.0
+          end do
+          do i = isice_lu+1,ismax_lu
+             landusef(i,iCell) = 0.0
+          end do
+!       else
+          ! Placeholder for land conversion with fractional land use
+!       endif
        isltyp(iCell) = 16
        snoalb(iCell) = 0.75
        vegfra(iCell) = 0._RKIND
        xland(iCell)  = 1._RKIND
+       lai_modis(iCell) = 0._RKIND
+
+       soiltypf(16,iCell) = 1.0  ! Currently defines dominant soil type fraction only
+       do i = 1,15
+          soiltypf(i,iCell) = 0.0
+       end do
 
        !... recalculate the soil temperature and soil moisture:
        do iSoil = 1, nSoilLevels
diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index b107d16953..08d82b0ae1 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -36,6 +36,10 @@
                      description="The number of first-guess atmospheric layers"/>
                 <dim name="nFGSoilLevels" definition="namelist:config_nfgsoillevels"
                      description="The number of first-guess soil layers"/>
+                <dim name="nLandCat" definition="namelist:config_nlandcat"
+                     description="The number of land use categories used by the land-surface scheme"/>
+                <dim name="nSoilCat" definition="16"
+                     description="The number of soil type categories used by the land-surface scheme"/>
                 <dim name="nVertLevelsP1" definition="nVertLevels+1"
                      description="The number of atmospheric levels, always one more than the number of layers"/>
                 <dim name="nSoilComps" definition="8"
@@ -63,9 +67,9 @@
                                         5 = super-cell, \newline
                                         6 = mountain wave, \newline
                                         7 = real-data initial conditions from, e.g., GFS, \newline
-                                        8 = surface field (SST, sea-ice) update file for use with real-data simulations \newline
-                                        9 = lateral boundary conditions update file for use with real-data simulations \newline
-                                        13 = CAM-MPAS 3-d grid with specified topography and zeta levels"
+                                        8 = surface field (SST, sea-ice) update file for use with real-data simulations, \newline
+                                        9 = lateral boundary conditions update file for use with real-data simulations, \newline
+                                       13 = CAM-MPAS 3-d grid with specified topography and zeta levels"
                      possible_values="1 -- 9"/>
 
                 <nml_option name="config_calendar_type" type="character" default_value="gregorian" in_defaults="false"
@@ -127,6 +131,11 @@
                      description="The number of vertical soil levels in the first-guess dataset (case 7 only)"
                      possible_values="Positive integer values"/>
 
+                <nml_option name="config_nlandcat"              type="integer"       default_value="40"
+                     units="-"
+                     description="The number of land categories"
+                     possible_values="Positive integer values"/>
+
                 <nml_option name="config_gocartlevels"          type="integer"       default_value="30"
                      units="-"
                      description="The number of vertical gocart levels in the climatological gocart data set"
@@ -164,7 +173,12 @@
                 <nml_option name="config_landuse_data"          type="character"     default_value="MODIFIED_IGBP_MODIS_NOAH"
                      units="-"
                      description="The land use classification to use (case 7 only)"
-                     possible_values="`USGS' or `MODIFIED_IGBP_MODIS_NOAH'"/>
+                     possible_values="`USGS' or `MODIFIED_IGBP_MODIS_NOAH' or `NLCD40'"/>
+
+                <nml_option name="config_frac_landuse"          type="logical"       default_value="false"
+                     units="-"
+                     description="Whether to use dominant (false) or fractional (true) land use and soil categories (case 7 only)"
+                     possible_values="true or false"/>
 
                 <nml_option name="config_topo_data"             type="character"     default_value="GMTED2010"
                      units="-"
@@ -240,6 +254,96 @@
                      description="Whether to use the vertical layer profile that was developed for use in real-time TC experiments"
                      possible_values="true or false"/>
 
+                <nml_option name="config_EPA_vertical_grid"     type="logical"       default_value="false"
+                     units="-"
+                     description="Whether to use the vertical layer profile originally used at EPA"
+                     possible_values="true or false"/>
+
+                <nml_option name="config_spline_vertical_grid"  type="logical"       default_value="false"
+                     units="-"
+                     description="Whether to use a vertical layer profile based on a cubic spline fit"
+                     possible_values="true or false"/>
+
+                <nml_option name="config_spline_x1"             type="real"          default_value="0.0"
+                     units="meters"
+                     description="spline control point #1 for height dimension"
+                     possible_values="Real values greater than or equal to zero"/>
+
+                <nml_option name="config_spline_y1"             type="real"          default_value="20.0"
+                     units="meters"
+                     description="spline control point #1 for layer thickness"
+                     possible_values="Real values greater than zero"/>
+
+                <nml_option name="config_spline_x2"             type="real"          default_value="350.0"
+                     units="meters"
+                     description="spline control point #2 for height dimension"
+                     possible_values="Real values greater than or equal to zero"/>
+
+                <nml_option name="config_spline_y2"             type="real"          default_value="100.0"
+                     units="meters"
+                     description="spline control point #2 for layer thickness"
+                     possible_values="Real values greater than zero"/>
+
+                <nml_option name="config_spline_x3"             type="real"          default_value="800.0"
+                     units="meters"
+                     description="spline control point #3 for height dimension"
+                     possible_values="Real values greater than or equal to zero"/>
+
+                <nml_option name="config_spline_y3"             type="real"          default_value="300.0"
+                     units="meters"
+                     description="spline control point #3 for layer thickness"
+                     possible_values="Real values greater than zero"/>
+
+                <nml_option name="config_spline_x4"             type="real"          default_value="2500.0"
+                     units="meters"
+                     description="spline control point #4 for height dimension"
+                     possible_values="Real values greater than or equal to zero"/>
+
+                <nml_option name="config_spline_y4"             type="real"          default_value="600.0"
+                     units="meters"
+                     description="spline control point #4 for layer thickness"
+                     possible_values="Real values greater than zero"/>
+
+                <nml_option name="config_spline_x5"             type="real"          default_value="8000.0"
+                     units="meters"
+                     description="spline control point #5 for height dimension"
+                     possible_values="Real values greater than or equal to zero"/>
+
+                <nml_option name="config_spline_y5"             type="real"          default_value="600.0"
+                     units="meters"
+                     description="spline control point #5 for layer thickness"
+                     possible_values="Real values greater than zero"/>
+
+                <nml_option name="config_spline_x6"             type="real"          default_value="16000.0"
+                     units="meters"
+                     description="spline control point #6 for height dimension"
+                     possible_values="Real values greater than or equal to zero"/>
+
+                <nml_option name="config_spline_y6"             type="real"          default_value="700.0"
+                     units="meters"
+                     description="spline control point #6 for layer thickness"
+                     possible_values="Real values greater than zero"/>
+
+                <nml_option name="config_spline_x7"             type="real"          default_value="30000.0"
+                     units="meters"
+                     description="spline control point #7 for height dimension"
+                     possible_values="Real values greater than or equal to zero"/>
+
+                <nml_option name="config_spline_y7"             type="real"          default_value="2000.0"
+                     units="meters"
+                     description="spline control point #7 for layer thickness"
+                     possible_values="Real values greater than zero"/>
+
+                <nml_option name="config_spline_x8"             type="real"          default_value="50000.0"
+                     units="meters"
+                     description="spline control point #8 for height dimension"
+                     possible_values="Real values greater than or equal to zero"/>
+
+                <nml_option name="config_spline_y8"             type="real"          default_value="4000.0"
+                     units="meters"
+                     description="spline control point #8 for layer thickness"
+                     possible_values="Real values greater than zero"/>
+
                 <nml_option name="config_specified_zeta_levels" type="character"     default_value=""    in_defaults="false"
                      units="-"
                      description="The name of a text file containing a list of vertical coordinate (zeta) values in increasing order at layer interfaces."
@@ -424,13 +528,17 @@
                         <var name="mminlu" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="isice_lu" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="iswater_lu" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
+                        <var name="ismax_lu" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="ivgtyp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
+                        <var name="landusef" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="isltyp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
+                        <var name="soiltypf" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="snoalb" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="soiltemp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="greenfrac" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="shdmin" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="shdmax" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
+                        <var name="lai12m" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="albedo12m" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="soilcomp" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
                         <var name="soilcl1" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
@@ -529,13 +637,17 @@
                         <var name="mminlu"/>
                         <var name="isice_lu"/>
                         <var name="iswater_lu"/>
+                        <var name="ismax_lu"/>
                         <var name="ivgtyp"/>
+                        <var name="landusef"/>
                         <var name="isltyp"/>
+                        <var name="soiltypf"/>
                         <var name="snoalb"/>
                         <var name="soiltemp"/>
                         <var name="greenfrac"/>
                         <var name="shdmin"/>
                         <var name="shdmax"/>
+                        <var name="lai12m"/>
                         <var name="albedo12m"/>
                         <var name="soilcomp"/>
                         <var name="soilcl1"/>
@@ -595,6 +707,7 @@
                         <var name="seaice" packages="met_stage_out"/>
                         <var name="vegfra" packages="met_stage_out"/>
                         <var name="sfc_albbck" packages="met_stage_out"/>
+                        <var name="lai_modis" packages="met_stage_out"/>
                         <var name="xland" packages="met_stage_out"/>
                         <var name="rho" packages="met_stage_out"/>
                         <var name="theta" packages="met_stage_out"/>
@@ -796,11 +909,17 @@
                 <var name="ivgtyp" name_in_code="lu_index" type="integer"  dimensions="nCells" units="unitless"
                      description="dominant vegetation category"/>
 
+                <var name="landusef" type="real" dimensions="nLandCat nCells" units="unitless"
+                     description="grid cell fractional land use"/>
+
                 <var name="mminlu" type="text" dimensions="" units="unitless"
                      description="land use classification"/>
 
-                <var name="isltyp" name_in_code="soilcat_top"  type="integer"  dimensions="nCells" units="unitless"
-                     description="dominate soil category"/>
+                <var name="isltyp" name_in_code="soilcat"  type="integer"  dimensions="nCells" units="unitless"
+                     description="dominant soil category"/>
+
+                <var name="soiltypf" type="real" dimensions="nSoilCat nCells" units="unitless"
+                     description="grid cell fractional soil type"/>
 
                 <var name="snoalb" type="real" dimensions="nCells" units="unitless"
                      description="annual maximum snow albedo"/>
@@ -808,16 +927,19 @@
                 <var name="soiltemp" type="real" dimensions="nCells" units="K"
                      description="annual mean deep soil temperature"/>
 
-                <var name="greenfrac" type="real" dimensions="nMonths nCells" units="unitless"
-                     description="monthly-mean climatological greeness fraction"/>
+                <var name="greenfrac" type="real" dimensions="nMonths nCells" units="percent"
+                     description="monthly-mean climatological greenness fraction"/>
 
-                <var name="shdmin" type="real" dimensions="nCells" units="unitless"
+                <var name="shdmin" type="real" dimensions="nCells" units="percent"
                      description="minimum fractional coverage of annual green vegetation fraction"/>
 
-                <var name="shdmax" type="real" dimensions="nCells" units="unitless"
+                <var name="shdmax" type="real" dimensions="nCells" units="percent"
                      description="maximum fractional coverage of annual green vegetation fraction"/>
 
-                <var name="albedo12m" type="real" dimensions="nMonths nCells" units="unitless"
+                <var name="lai12m" type="real" dimensions="nMonths nCells" units="m^{2} m^{-2}"
+                     description="monthly-mean climatological leaf area index"/>
+
+                <var name="albedo12m" type="real" dimensions="nMonths nCells" units="percent"
                      description="monthly-mean climatological surface albedo"/>
 
                 <var name="isice_lu" type="integer" dimensions="" units="unitless" default_value="24"
@@ -826,6 +948,9 @@
                 <var name="iswater_lu" type="integer" dimensions="" units="unitless" default_value="16"
                      description="Index category for water"/>
 
+                <var name="ismax_lu" type="integer" dimensions="" units="unitless" default_value="24"
+                     description="Maximum number of land use categories"/>
+
                 <!-- SOIL COMPOSITION fields needed for the NOAH-MP land surface scheme -->
                 <var name="soilcomp" type="real" dimensions="nSoilComps nCells" units="percent"
                      description="soil composition needed as input in the NOAH-MP land surface model"
@@ -1218,7 +1343,7 @@
                      description="geopotential height vertically interpolated from first guess"
                      packages="met_stage_out"/>
 
-                <var name="vegfra" type="real" dimensions="nCells Time" units="unitless"
+                <var name="vegfra" type="real" dimensions="nCells Time" units="percent"
                      description="vegetation fraction"
                      packages="met_stage_out"/>
 
@@ -1226,6 +1351,10 @@
                      description="background surface albedo"
                      packages="met_stage_out"/>
 
+                <var name="lai_modis" type="real" dimensions="nCells Time" units="m^{2} m^{-2}"
+                     description="leaf area index"
+                     packages="met_stage_out"/>
+
                 <var name="xland" type="real" dimensions="nCells Time" units="unitless"
                      description="land-ocean mask (1=land including sea-ice ; 2=ocean)"
                      packages="met_stage_out"/>
diff --git a/src/core_init_atmosphere/mpas_init_atm_cases.F b/src/core_init_atmosphere/mpas_init_atm_cases.F
index 7d43b5ee83..f60ff5c729 100644
--- a/src/core_init_atmosphere/mpas_init_atm_cases.F
+++ b/src/core_init_atmosphere/mpas_init_atm_cases.F
@@ -2668,7 +2668,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
       real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
 
       real (kind=RKIND) :: target_z
-      integer :: iCell, iCell1, iCell2 , iEdge, i, k, nz, cell1, cell2
+      integer :: iCell, iCell1, iCell2 , iEdge, i, k, nz, itr, itrp, cell1, cell2
       integer, pointer :: nCellsSolve, nz1
       integer :: nInterpPoints, ndims
 
@@ -2715,7 +2715,7 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
 
       real (kind=RKIND) :: p0
 
-      real (kind=RKIND) :: etavs, ztemp, zd, zt, dz, str
+      real (kind=RKIND) :: etavs, ztemp, zd, zt, dz, str, grd, kfrac
 
       real (kind=RKIND) :: es, rs, xnutr, znut, rcv
 
@@ -2744,6 +2744,24 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
       real (kind=RKIND), pointer :: config_ztop
       logical, pointer :: config_tc_vertical_grid
       character (len=StrKIND), pointer :: config_specified_zeta_levels
+      logical, pointer :: config_spline_vertical_grid
+      real (kind=RKIND), pointer :: config_spline_x1
+      real (kind=RKIND), pointer :: config_spline_y1
+      real (kind=RKIND), pointer :: config_spline_x2
+      real (kind=RKIND), pointer :: config_spline_y2
+      real (kind=RKIND), pointer :: config_spline_x3
+      real (kind=RKIND), pointer :: config_spline_y3
+      real (kind=RKIND), pointer :: config_spline_x4
+      real (kind=RKIND), pointer :: config_spline_y4
+      real (kind=RKIND), pointer :: config_spline_x5
+      real (kind=RKIND), pointer :: config_spline_y5
+      real (kind=RKIND), pointer :: config_spline_x6
+      real (kind=RKIND), pointer :: config_spline_y6
+      real (kind=RKIND), pointer :: config_spline_x7
+      real (kind=RKIND), pointer :: config_spline_y7
+      real (kind=RKIND), pointer :: config_spline_x8
+      real (kind=RKIND), pointer :: config_spline_y8
+      logical, pointer :: config_EPA_vertical_grid
       logical, pointer :: config_use_spechumd
       integer, pointer :: config_nfglevels
       integer, pointer :: config_nfgsoillevels
@@ -2807,6 +2825,31 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
       real (kind=RKIND), dimension(:), pointer :: rh2
       real (kind=RKIND), dimension(:), pointer :: t2m
 
+      ! For cubic spline definition of air quality layer thickness profile
+      integer, parameter :: nprof = 8
+      real (kind=RKIND), dimension(nprof) :: &
+        xx,      &!< x-dimension control point for cubic spline interpolation
+        yy        !< y-dimension control point for cubic spline interpolation
+      real (kind=RKIND), parameter :: &
+        cadj = 0.3,     &!< layer thickness adjustment coefficient
+        ctol = 0.001     !< tolerance for layer fit to model top
+      real (kind=RKIND), dimension(nprof) :: &
+        ya,      &!< Thickness scaled to match model top
+        y2ndDer, &!< 2nd derivative at nodes of layer profile definition
+        acoef     !< Tridiagonal decomposition factor
+      real (kind=RKIND) :: &
+        thick,   &!< Interpolated layer thickness
+        zint,    &!< Height for interpolation
+        ptef,    &!< Profile top error factor
+        adj,     &!< Thickness adjustment factor
+        xRatio,  &!< X-dimension scale for interpolation
+        tvar,    &!< temporary variable
+        aa,      &!< temporary variable
+        bb,      &!< temporary variable
+        cc        !< temporary variable
+      real (kind=RKIND), dimension(nVertLevels+1) :: zwtemp
+      integer :: nlayerint, klo, khi, ncount
+
       character (len=StrKIND) :: errstring
       character (len=StrKIND), pointer :: config_interface_projection
 
@@ -2820,6 +2863,24 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
       call mpas_pool_get_config(configs, 'config_ztop', config_ztop)
       call mpas_pool_get_config(configs, 'config_tc_vertical_grid', config_tc_vertical_grid)
       call mpas_pool_get_config(configs, 'config_specified_zeta_levels', config_specified_zeta_levels)
+      call mpas_pool_get_config(configs, 'config_EPA_vertical_grid', config_EPA_vertical_grid)
+      call mpas_pool_get_config(configs, 'config_spline_vertical_grid', config_spline_vertical_grid)
+      call mpas_pool_get_config(configs, 'config_spline_x1', config_spline_x1)
+      call mpas_pool_get_config(configs, 'config_spline_y1', config_spline_y1)
+      call mpas_pool_get_config(configs, 'config_spline_x2', config_spline_x2)
+      call mpas_pool_get_config(configs, 'config_spline_y2', config_spline_y2)
+      call mpas_pool_get_config(configs, 'config_spline_x3', config_spline_x3)
+      call mpas_pool_get_config(configs, 'config_spline_y3', config_spline_y3)
+      call mpas_pool_get_config(configs, 'config_spline_x4', config_spline_x4)
+      call mpas_pool_get_config(configs, 'config_spline_y4', config_spline_y4)
+      call mpas_pool_get_config(configs, 'config_spline_x5', config_spline_x5)
+      call mpas_pool_get_config(configs, 'config_spline_y5', config_spline_y5)
+      call mpas_pool_get_config(configs, 'config_spline_x6', config_spline_x6)
+      call mpas_pool_get_config(configs, 'config_spline_y6', config_spline_y6)
+      call mpas_pool_get_config(configs, 'config_spline_x7', config_spline_x7)
+      call mpas_pool_get_config(configs, 'config_spline_y7', config_spline_y7)
+      call mpas_pool_get_config(configs, 'config_spline_x8', config_spline_x8)
+      call mpas_pool_get_config(configs, 'config_spline_y8', config_spline_y8)
       call mpas_pool_get_config(configs, 'config_use_spechumd', config_use_spechumd)
       call mpas_pool_get_config(configs, 'config_nfglevels', config_nfglevels)
       call mpas_pool_get_config(configs, 'config_nfgsoillevels', config_nfgsoillevels)
@@ -3113,6 +3174,127 @@ subroutine init_atm_case_gfs(block, mesh, nCells, nEdges, nVertLevels, fg, state
             if (k > 1) dzw(k-1) = zw(k)-zw(k-1)
          end do
 
+      else if (config_spline_vertical_grid) then
+
+         call mpas_log_write('Setting up vertical levels using a cubic spline fit')
+
+         zt = config_ztop
+         xx(1) = config_spline_x1
+         yy(1) = config_spline_y1
+         xx(2) = config_spline_x2
+         yy(2) = config_spline_y2
+         xx(3) = config_spline_x3
+         yy(3) = config_spline_y3
+         xx(4) = config_spline_x4
+         yy(4) = config_spline_y4
+         xx(5) = config_spline_x5
+         yy(5) = config_spline_y5
+         xx(6) = config_spline_x6
+         yy(6) = config_spline_y6
+         xx(7) = config_spline_x7
+         yy(7) = config_spline_y7
+         xx(8) = config_spline_x8
+         yy(8) = config_spline_y8
+
+! Set up cubic spline definition of layer thickness vertical profile
+         y2ndDer(1)=0.0
+         y2ndDer(nprof)=0.0
+         acoef(1)=0.0
+         ya(1:nprof) = yy(1:nprof)
+
+         do i=2,nprof-1
+            xRatio=(xx(i)-xx(i-1))/(xx(i+1)-xx(i-1))
+            tvar=1.0/(2.0+xRatio*y2ndDer(i-1))
+            y2ndDer(i)=tvar*(xRatio-1.0)
+            acoef(i) = tvar*(6.0*((ya(i+1)-ya(i))/(xx(i+1)-xx(i)) &
+                -(ya(i)-ya(i-1))/(xx(i)-xx(i-1)))/(xx(i+1)-xx(i-1)) &
+                -xRatio*acoef(i-1))
+         enddo
+
+         do i=nprof-1,1,-1
+            y2ndDer(i)=y2ndDer(i)*y2ndDer(i+1)+acoef(i)
+         enddo
+
+! Loop to set layer interface heights from air quality profile
+         ptef = 0.
+         zwtemp(1)= 0.
+         ncount = 0
+         do while(abs(1.0-ptef) > ctol)
+            zwtemp(2) = ya(1)
+            zint = ya(1)
+            nlayerint = 2
+            do while (nlayerint.lt.nz)
+               klo=1
+               khi=nprof
+               do while (khi-klo.gt.1)
+                  k=(khi+klo)/2
+                  if(xx(k).gt.zint)then
+                     khi=k
+                  else
+                     klo=k
+                  endif
+               enddo
+               cc=xx(khi)-xx(klo)
+               aa=(xx(khi)-zint)/cc
+               bb=(zint-xx(klo))/cc
+               thick=aa*ya(klo)+bb*ya(khi)+((aa**3-aa)*y2ndDer(klo)+ &
+                  (bb**3-bb)*y2ndDer(khi))*(cc**2)/6.
+               zint = zint+thick
+               nlayerint = nlayerint+1
+               zwtemp(nlayerint) = zint
+            enddo
+            if(ptef.eq.0.) then
+               ptef = zt/zint
+               if(ptef > 1.5) then
+                  call mpas_log_write('*** Significant layer inflation required ***')
+                  call mpas_log_write('***   You may want to use more layers    ***')
+                  call mpas_log_write('***      and/or a lower model top        ***')
+               else if(ptef < 0.6) then
+                  call mpas_log_write('*** Significant layer deflation required ***')
+                  call mpas_log_write('***   You may want to use fewer layers   ***')
+                  call mpas_log_write('***      and/or a higher model top       ***')
+               endif
+            else
+               ptef = zt/zint
+            endif
+            call mpas_log_write('Using first layer, ptef = $r $r',realArgs=(/ya(1),ptef/))
+            ncount = ncount+1
+            if(ncount > 20) then
+               call mpas_log_write('*** ERROR: layer adjustment failure ***')
+               call mpas_log_write('  *** Iteration did not converge ***')
+               call mpas_dmpar_global_abort('****************************************************************************', deferredAbort=.true.)
+               call mpas_dmpar_global_abort('ERROR: Aborted while setting up vertical levels for Air Quality applications', deferredAbort=.true.)
+               call mpas_dmpar_global_abort('****************************************************************************')
+            endif
+            adj = 1.0+(cadj*(ptef-1.0))
+            ya(1:nprof) = ya(1:nprof)*adj
+            y2ndDer(1:nprof) = y2ndDer(1:nprof)*adj
+         enddo
+
+! Now make minor adjustment to layers to fit within specified top
+         do k=1,nz
+            zw(k) = zwtemp(k)*ptef
+         enddo
+         call mpas_log_write('Layers for Air Quality applications generated')
+         call mpas_log_write('Number of layers and Model top = $i $r',intArgs=(/nz1/),realArgs=(/zt/))
+         call mpas_log_write('Adjusted bottom layer thickness (m) = $r',realArgs=(/zw(2)/))
+         call mpas_log_write('')
+
+      else if (config_EPA_vertical_grid) then
+
+         call mpas_log_write('Setting up vertical levels as in MPAS 4.0 and earlier, but')
+         call mpas_log_write('*** Using custom layer definition ***')
+
+         str = 3.
+         grd = 0.03
+         zt = config_ztop
+
+         do k=1,nz
+            kfrac = real(k-1)/real(nz1)
+            zw(k) = zt*((1-grd)*kfrac**str+grd*kfrac)
+            if (k > 1) dzw(k-1) = zw(k)-zw(k-1)
+         end do
+
       !
       ! Otherwise, use the vertical level configuration from MPAS v2.0
       !
@@ -6620,32 +6802,39 @@ subroutine physics_idealized_init(mesh, fg)
       type (mpas_pool_type), intent(inout) :: fg
       
       !local variables:
-      integer :: iCell, iMonth, iSoil
-      integer, pointer :: nCells, nSoilLevels, nMonths
-      integer, dimension(:), pointer :: landmask, lu_index, soilcat_top
-      real (kind=RKIND), dimension(:), pointer :: ter, xice, shdmin, shdmax, vegfra, sfc_albbck, xland, seaice
+      integer :: iCell, iMonth, iSoil, iLandCat, iSoilCat
+      integer, pointer :: nCells, nSoilLevels, nMonths, nLandCat, nSoilCat
+      integer, dimension(:), pointer :: landmask, lu_index, soilcat
+      real (kind=RKIND), dimension(:), pointer :: ter, xice, shdmin, shdmax, vegfra, sfc_albbck, xland, seaice, lai_modis
       real (kind=RKIND), dimension(:), pointer :: snow, snowc, snoalb, snowh, skintemp, sst, tmn
-      real (kind=RKIND), dimension(:,:), pointer :: tslb, smcrel, sh2o, smois, dzs, albedo12m, greenfrac
+      real (kind=RKIND), dimension(:,:), pointer :: landusef, soiltypf
+      real (kind=RKIND), dimension(:,:), pointer :: tslb, smcrel, sh2o, smois, dzs, albedo12m, greenfrac, lai12m
       
       !---------------------------------------------------------------------------------------------
 
       call mpas_pool_get_dimension(mesh, 'nCells', nCells)
       call mpas_pool_get_dimension(mesh, 'nSoilLevels', nSoilLevels)
       call mpas_pool_get_dimension(mesh, 'nMonths', nMonths)
+      call mpas_pool_get_dimension(mesh, 'nLandCat', nLandCat)
+      call mpas_pool_get_dimension(mesh, 'nSoilCat', nSoilCat)
 
       call mpas_pool_get_array(mesh, 'ter', ter)
       call mpas_pool_get_array(mesh, 'landmask', landmask)
       call mpas_pool_get_array(mesh, 'lu_index', lu_index)
-      call mpas_pool_get_array(mesh, 'soilcat_top', soilcat_top)
+      call mpas_pool_get_array(mesh, 'landusef', landusef)
+      call mpas_pool_get_array(mesh, 'soilcat', soilcat)
+      call mpas_pool_get_array(mesh, 'soiltypf', soiltypf)
       call mpas_pool_get_array(mesh, 'shdmin', shdmin)
       call mpas_pool_get_array(mesh, 'shdmax', shdmax)
       call mpas_pool_get_array(mesh, 'snoalb', snoalb)
       call mpas_pool_get_array(mesh, 'albedo12m', albedo12m)
       call mpas_pool_get_array(mesh, 'greenfrac', greenfrac)
+      call mpas_pool_get_array(mesh, 'lai12m', lai12m)
 
       call mpas_pool_get_array(fg, 'xice', xice)
       call mpas_pool_get_array(fg, 'vegfra', vegfra)
       call mpas_pool_get_array(fg, 'sfc_albbck', sfc_albbck)
+      call mpas_pool_get_array(fg, 'lai_modis', lai_modis)
       call mpas_pool_get_array(fg, 'xland', xland)
       call mpas_pool_get_array(fg, 'seaice', seaice)
       call mpas_pool_get_array(fg, 'snow', snow)
@@ -6671,11 +6860,12 @@ subroutine physics_idealized_init(mesh, fg)
          xice(iCell) = 0.0
          landmask(iCell) = 0
          lu_index(iCell) = 0
-         soilcat_top(iCell) = 0
+         soilcat(iCell) = 0
          shdmin(iCell) = 0.0
          shdmax(iCell) = 0.0
          vegfra(iCell) = 0.0
          sfc_albbck(iCell) = 0.0
+         lai_modis(iCell) = 0.0
          xland(iCell) = 0.0
          seaice(iCell) = 0.0
       
@@ -6698,11 +6888,22 @@ subroutine physics_idealized_init(mesh, fg)
             smois(iSoil,iCell) =   0.0
             dzs(iSoil,iCell) =   0.0
          end do
+
+         !fractional soil type:
+         do iSoilCat = 1, nSoilCat
+            soiltypf(iSoilCat,iCell) = 0.0
+         end do
+      
+         !fractional land use:
+         do iLandCat = 1, nLandCat
+            landusef(iLandCat,iCell) = 0.0
+         end do
       
-         !monthly climatological surface albedo and greeness fraction:
+         !monthly climatological surface albedo, greeness fraction, and leaf area index:
          do iMonth = 1, nMonths
             albedo12m(iMonth,iCell) = 0.08
             greenfrac(iMonth,iCell) = 0.0
+            lai12m(iMonth,iCell)    = 0.0
          end do
       
       end do
diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F
index b9d289d8b6..f7a8e13d01 100644
--- a/src/core_init_atmosphere/mpas_init_atm_static.F
+++ b/src/core_init_atmosphere/mpas_init_atm_static.F
@@ -88,11 +88,13 @@ end subroutine interp_accumulation_function
  ! 12.1.2.2 of the Fortran standard) and currently, only PGI does not support this, so
  ! use module level variables for now...
  !
+ real (kind=RKIND), dimension(:,:), pointer :: landusef
+ real (kind=RKIND), dimension(:,:), pointer :: soiltypf
  integer (kind=I8KIND), dimension(:), pointer :: ter_integer
  integer (kind=I8KIND), dimension(:,:), pointer :: soilcomp_int
  real (kind=RKIND) :: soilcomp_msgval = 255.0_RKIND   ! Modified later based on index file for soilcomp
  integer, dimension(:), pointer :: lu_index
- integer, dimension(:), pointer :: soilcat_top
+ integer, dimension(:), pointer :: soilcat
  integer, dimension(:), pointer :: nhs
  integer, dimension(:,:), allocatable:: ncat
  ! Landmask is used by the accumulation function for maxsnoalb and soilcomp,
@@ -153,6 +155,7 @@ subroutine init_atm_static(mesh, dims, configs)
 
  integer, pointer :: isice_lu, iswater_lu
  integer :: iswater_soil
+ integer, pointer :: ismax_lu
  integer, pointer :: nCells, nCellsSolve, nEdges, nVertices, maxEdges
  logical, pointer :: on_a_sphere
  real (kind=RKIND), pointer :: sphere_radius
@@ -177,12 +180,16 @@ subroutine init_atm_static(mesh, dims, configs)
  integer (kind=I8KIND), dimension(:), pointer :: snoalb_integer
  real (kind=RKIND), dimension(:), pointer :: shdmin, shdmax
  real (kind=RKIND), dimension(:,:), pointer :: greenfrac
+ real (kind=RKIND), dimension(:,:), pointer :: lai12m
+ integer (kind=I8KIND), dimension(:,:), pointer :: lai12m_int
  real (kind=RKIND), dimension(:,:), pointer :: albedo12m
  integer (kind=I8KIND), dimension(:,:), pointer :: albedo12m_int
+ real (kind=RKIND), dimension(:,:), pointer :: landusef
+ real (kind=RKIND), dimension(:,:), pointer :: soiltypf
  real (kind=RKIND) :: fillval
  real (kind=RKIND), pointer :: missing_value
  integer, dimension(:), pointer :: lu_index
- integer, dimension(:), pointer :: soilcat_top
+ integer, dimension(:), pointer :: soilcat
  integer, dimension(:), pointer :: landmask
  integer, dimension(:), pointer :: bdyMaskCell
  character(len=StrKIND), pointer :: mminlu
@@ -262,13 +269,17 @@ subroutine init_atm_static(mesh, dims, configs)
  call mpas_pool_get_array(mesh, 'verticesOnCell', verticesOnCell)
 
  call mpas_pool_get_array(mesh, 'lu_index', lu_index)
+ call mpas_pool_get_array(mesh, 'landusef', landusef)
  call mpas_pool_get_array(mesh, 'mminlu', mminlu)
  call mpas_pool_get_array(mesh, 'isice_lu', isice_lu)
  call mpas_pool_get_array(mesh, 'iswater_lu', iswater_lu)
- call mpas_pool_get_array(mesh, 'soilcat_top', soilcat_top)
+ call mpas_pool_get_array(mesh, 'ismax_lu', ismax_lu)
+ call mpas_pool_get_array(mesh, 'soilcat', soilcat)
+ call mpas_pool_get_array(mesh, 'soiltypf', soiltypf)
  call mpas_pool_get_array(mesh, 'landmask', landmask)
  call mpas_pool_get_array(mesh, 'snoalb', snoalb)
  call mpas_pool_get_array(mesh, 'greenfrac', greenfrac)
+ call mpas_pool_get_array(mesh, 'lai12m', lai12m)
  call mpas_pool_get_array(mesh, 'albedo12m', albedo12m)
  call mpas_pool_get_array(mesh, 'shdmin', shdmin)
  call mpas_pool_get_array(mesh, 'shdmax', shdmax)
@@ -352,11 +363,14 @@ subroutine init_atm_static(mesh, dims, configs)
        write(mminlu,'(a)') 'USGS'
     case('MODIFIED_IGBP_MODIS_NOAH')
        write(mminlu,'(a)') 'MODIFIED_IGBP_MODIS_NOAH'
+    case('NLCD40')
+       write(mminlu,'(a)') 'NLCD40'
     case default
          call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
          call mpas_log_write('Invalid land use dataset '''//trim(config_landuse_data) &
                                        //''' selected for config_landuse_data',                   messageType=MPAS_LOG_ERR)
-         call mpas_log_write('   Possible options are: ''USGS'', ''MODIFIED_IGBP_MODIS_NOAH''',   messageType=MPAS_LOG_ERR)
+         call mpas_log_write('   Possible options are: ''USGS'', ''MODIFIED_IGBP_MODIS_NOAH'', ', messageType=MPAS_LOG_ERR)
+         call mpas_log_write('                         or ''NLCD40'''                           , messageType=MPAS_LOG_ERR)
          call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
          call mpas_log_write('Please correct the namelist.', messageType=MPAS_LOG_CRIT)
  end select surface_input_select0
@@ -396,48 +410,80 @@ subroutine init_atm_static(mesh, dims, configs)
     case('MODIFIED_IGBP_MODIS_NOAH')
        call mpas_log_write('Using 20-class MODIS 30-arc-second land cover dataset')
        geog_sub_path = 'modis_landuse_20class_30s/'
+    case('NLCD40')
+       call mpas_log_write('Using 20-class MODIS 30-arc-second land cover dataset for most of globe')
+       geog_sub_path = 'modis_landuse_20class_30s/'
     case default
        call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
        call mpas_log_write('Invalid land use dataset '''//trim(config_landuse_data) &
                                      //''' selected for config_landuse_data',                   messageType=MPAS_LOG_ERR)
-       call mpas_log_write('   Possible options are: ''USGS'', ''MODIFIED_IGBP_MODIS_NOAH''',   messageType=MPAS_LOG_ERR)
+       call mpas_log_write('   Possible options are: ''USGS'', ''MODIFIED_IGBP_MODIS_NOAH'', ', messageType=MPAS_LOG_ERR)
+       call mpas_log_write('                         or ''NLCD40'''                           , messageType=MPAS_LOG_ERR)
        call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR)
        call mpas_log_write('Please correct the namelist.', messageType=MPAS_LOG_CRIT)
  end select surface_input_select1
 
  call mpas_log_write('--- start interpolate LU_INDEX')
- call interp_landuse(mesh, tree, trim(geog_data_path)//trim(geog_sub_path), isice_lu, iswater_lu, &
+ call interp_landuse(configs, mesh, tree, trim(geog_data_path)//trim(geog_sub_path), isice_lu, iswater_lu, ismax_lu, &
                      supersample_fac=supersample_fac_30s)
  call mpas_log_write('--- end interpolate LU_INDEX')
 
+ if (config_landuse_data == 'NLCD40') then
+    ! Incorporate NLCD fractional percentage land use for CONUS, with MODIS for rest of world
+    call mpas_log_write('Overlaying 40-class NLCD 9-arc-second land cover dataset for CONUS')
+    geog_sub_path = 'nlcd2011_ll_9s/'
+    call interp_landuse2(mesh, dims, configs, trim(config_geog_data_path)//trim(geog_sub_path), isice_lu, iswater_lu, ismax_lu)
+ end if
+
 !
-! Interpolate SOILCAT_TOP
+! Interpolate SOILCAT
 !
- geog_sub_path = 'soiltype_top_30s/'
+ surface_input_select2: select case(trim(config_landuse_data))
+    case('USGS')
+       geog_sub_path = 'soiltype_top_30s/'
+    case('MODIFIED_IGBP_MODIS_NOAH')
+       geog_sub_path = 'soiltype_top_30s/'
+    case('NLCD40')
+       geog_sub_path = 'soiltype_bot_30s/'
+    case default
+ end select surface_input_select2
 
- call mpas_log_write('--- start interpolate SOILCAT_TOP')
- call interp_soilcat(mesh, tree, trim(geog_data_path)//trim(geog_sub_path), iswater_soil, &
+ call mpas_log_write('--- start interpolate SOILCAT')
+ call interp_soilcat(configs, mesh, tree, trim(config_geog_data_path)//trim(geog_sub_path), iswater_soil, &
                      supersample_fac=supersample_fac_30s)
- call mpas_log_write('--- end interpolate SOILCAT_TOP')
+ call mpas_log_write('--- end interpolate SOILCAT')
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! KLUDGE TO FIX SOIL TYPE OVER ANTARCTICA
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- where (lu_index == isice_lu) soilcat_top = 16
+! where (lu_index == isice_lu) soilcat = 16
+ do iCell = 1,nCells
+    if (lu_index(iCell) == isice_lu) then
+       do i = 1,16
+          soiltypf(i,iCell) = 0.0
+       end do
+       soilcat(iCell) = 16
+       soiltypf(soilcat(iCell),iCell) = 1.0
+    end if
+ end do
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! CORRECT INCONSISTENT SOIL AND LAND USE DATA
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  do iCell = 1,nCells
     if (lu_index(iCell) == iswater_lu .or. &
-        soilcat_top(iCell) == iswater_soil) then
+        soilcat(iCell) == iswater_soil) then
         if (lu_index(iCell) /= iswater_lu) then
             call mpas_log_write('Turning lu_index into water at $i', intArgs=(/iCell/))
+            landusef(lu_index(iCell),iCell) = 0.0
             lu_index(iCell) = iswater_lu
+            landusef(lu_index(iCell),iCell) = 1.0  ! Currently defines dominant land use fraction only
+            soiltypf(soilcat(iCell),iCell) = 1.0   ! Currently defines dominant soil type fraction only
         end if
-        if (soilcat_top(iCell) /= iswater_soil) then
-            call mpas_log_write('Turning soilcat_top into water at $i', intArgs=(/iCell/))
-            soilcat_top(iCell) = iswater_soil
+        if (soilcat(iCell) /= iswater_soil) then
+            call mpas_log_write('Turning soilcat into water at $i', intArgs=(/iCell/))
+            soilcat(iCell) = iswater_soil
+            soiltypf(soilcat(iCell),iCell) = 1.0   ! Currently defines dominant soil type fraction only
         end if
     end if
  end do
@@ -943,6 +989,150 @@ subroutine init_atm_static(mesh, dims, configs)
  call mpas_log_write('--- end interpolate GREENFRAC')
 
 
+!
+! Interpolate LAI
+!
+
+ call mpas_log_write('Using MODIS LAI 30-arc-second data for climatological monthly leaf area index')
+
+ geog_sub_path = 'lai_modis_30s/'
+
+ ierr = mgr % init(trim(config_geog_data_path)//trim(geog_sub_path))
+ if (ierr /= 0) then
+    call mpas_log_write('Error occurred when initalizing the interpolation of monthly leaf area index (lai12m)', &
+                         messageType=MPAS_LOG_CRIT)
+ endif
+
+ call mpas_pool_get_config(mgr % pool, 'tile_bdr', tile_bdr)
+ call mpas_pool_get_config(mgr % pool, 'tile_x', tile_nx)
+ call mpas_pool_get_config(mgr % pool, 'tile_y', tile_ny)
+ call mpas_pool_get_config(mgr % pool, 'tile_z', tile_nz)
+ call mpas_pool_get_config(mgr % pool, 'missing_value', missing_value)
+ call mpas_pool_get_config(mgr % pool, 'scale_factor', scalefactor_ptr)
+ scalefactor = scalefactor_ptr
+
+ allocate(nhs(nCells))
+ allocate(lai12m_int(tile_nz, nCells))
+ nhs(:) = 0
+ lai12m(:,:) = 0.0
+ lai12m_int(:,:) = 0_I8KIND
+ fillval = 0.0
+
+ do iCell = 1, nCells
+     if (nhs(iCell) == 0) then
+         tile => null()
+         ierr = mgr % get_tile(latCell(iCell), lonCell(iCell), tile)
+         if (ierr /= 0 .or. .not. associated(tile)) then
+             call mpas_log_write('Could not get tile that contained cell $i', intArgs=(/iCell/), messageType=MPAS_LOG_CRIT)
+         end if
+
+         ierr = mgr % push_tile(tile)
+         if (ierr /= 0) then
+             call mpas_log_write("Error pushing this tile onto the stack: "//trim(tile % fname), messageType=MPAS_LOG_CRIT)
+         end if
+     end if
+
+     do while (.not. mgr % is_stack_empty())
+         tile => mgr % pop_tile()
+
+         if (tile % is_processed) then
+             cycle
+         end if
+
+         call mpas_log_write('Processing tile: '//trim(tile % fname))
+
+         all_pixels_mapped_to_halo_cells = .true.
+
+         do j = tile_bdr + 1, tile_ny + tile_bdr, 1
+             do i = tile_bdr + 1, tile_nx + tile_bdr, 1
+
+                 call mgr % tile_to_latlon(tile, j, i, lat_pt, lon_pt)
+                 call mpas_latlon_to_xyz(xPixel, yPixel, zPixel, sphere_radius, lat_pt, lon_pt)
+                 call mpas_kd_search(tree, (/xPixel, yPixel, zPixel/), res, max_distance=max_kdtree_distance2)
+
+                 if (.not. associated(res)) cycle
+
+                 !
+                 ! This field only matters for land cells, and for all but the outermost boundary cells,
+                 ! we can safely assume that the nearest model grid cell contains the pixel (else, a different
+                 ! cell would be nearest)
+                 !
+                 if (landMask(res % id) == 1 .and. bdyMaskCell(res % id) < nBdyLayers) then
+                     do k = 1, tile_nz
+                         if (tile % tile(i, j, k) == missing_value) then
+                             i8val = int(fillval, kind=I8KIND)
+                         else
+                             i8val = int(tile % tile(i,j,k), kind=I8KIND)
+                         end if
+                         lai12m_int(k, res % id) = lai12m_int(k, res % id) + i8val
+                     end do
+                     nhs(res % id) = nhs(res % id) + 1
+
+                     if (res % id <= nCellsSolve) then
+                         all_pixels_mapped_to_halo_cells = .false.
+                     end if
+
+                 ! For outermost land cells, additional work is needed to verify that the pixel
+                 ! actually lies within the nearest cell
+                 else if (landMask(res % id) == 1) then
+                     if (mpas_in_cell(xPixel, yPixel, zPixel, xCell(res % id), yCell(res % id), zCell(res % id), &
+                                       nEdgesOnCell(res % id), verticesOnCell(:,res % id), xVertex, yVertex, zVertex)) then
+                         do k = 1, tile_nz
+                             if (tile % tile(i, j, k) == missing_value) then
+                                 i8val = int(fillval, kind=I8KIND)
+                             else
+                                 i8val = int(tile % tile(i,j,k), kind=I8KIND)
+                             end if
+                             lai12m_int(k, res % id) = lai12m_int(k, res % id) + i8val
+                         end do
+                         nhs(res % id) = nhs(res % id) + 1
+
+                         if (res % id <= nCellsSolve) then
+                             all_pixels_mapped_to_halo_cells = .false.
+                         end if
+                     end if
+                 end if
+             end do
+         end do
+
+         tile % is_processed = .true.
+
+         if (.not. all_pixels_mapped_to_halo_cells) then
+             ierr = mgr % push_neighbors(tile)
+             if (ierr /= 0) then
+                 call mpas_log_write("Error pushing the tile neighbors of: "//trim(tile%fname), messageType=MPAS_LOG_CRIT)
+             end if
+         end if
+
+     end do
+ end do
+
+ do iCell = 1, nCells
+     ! For land points that have no overlap with valid data, and for water points,
+     ! just use the fill value...
+     if (nhs(iCell) == 0) then
+         lai12m(:,iCell) = fillval
+     else
+         lai12m(:,iCell) = real(real(lai12m_int(:,iCell), kind=R8KIND) / real(nhs(iCell), kind=R8KIND), kind=RKIND)
+         lai12m(:,iCell) = lai12m(:,iCell) * scalefactor
+     end if
+     if (lu_index(iCell) == isice_lu) then
+         lai12m(:,iCell) = fillval
+     end if
+ end do
+
+ deallocate(nhs)
+ deallocate(lai12m_int)
+
+ ierr = mgr % finalize()
+ if (ierr /= 0) then
+    call mpas_log_write('Error occurred when finalizing the interpolation of monthly leaf area index (lai12m)', &
+                         messageType=MPAS_LOG_CRIT)
+ endif
+
+ call mpas_log_write('--- end interpolate LAI')
+
+
 !
 ! Interpolate ALBEDO12M
 !
@@ -1769,16 +1959,18 @@ end subroutine categorical_interp_accumulation
  !> that isice and iswater are in the dataset's index file.
  !
  !-----------------------------------------------------------------------
- subroutine interp_landuse(mesh, kdtree, geog_data_path, isice_lu, iswater_lu, supersample_fac)
+ subroutine interp_landuse(configs, mesh, kdtree, geog_data_path, isice_lu, iswater_lu, ismax_lu, supersample_fac)
 
      implicit none
 
      ! Input variables
      type (mpas_pool_type), intent(inout) :: mesh
      type (mpas_kd_type), pointer, intent(in) :: kdtree
+     type (mpas_pool_type), intent(in) :: configs
      character (len=*), intent(in) :: geog_data_path
      integer, intent(out) :: isice_lu
      integer, intent(out) :: iswater_lu
+     integer, intent(out) :: ismax_lu
      integer, intent(in), optional :: supersample_fac
 
      ! Local variables
@@ -1787,9 +1979,12 @@ subroutine interp_landuse(mesh, kdtree, geog_data_path, isice_lu, iswater_lu, su
      integer, pointer :: isice_lu_ptr
      integer, pointer :: iswater_lu_ptr
 
+     logical, pointer :: config_frac_landuse
+
      real (kind=RKIND), pointer :: scalefactor
 
-     integer :: iCell
+     integer :: ncattot
+     integer :: i, iCell
      integer :: ierr
 
      ierr = mgr % init(trim(geog_data_path))
@@ -1798,8 +1993,11 @@ subroutine interp_landuse(mesh, kdtree, geog_data_path, isice_lu, iswater_lu, su
          return ! Program execution should not reach this statement since the preceding message is a critical error
      end if
 
+     call mpas_pool_get_config(configs, 'config_frac_landuse', config_frac_landuse)
+
      call mpas_pool_get_dimension(mesh, 'nCells', nCells)
      call mpas_pool_get_array(mesh, 'lu_index', lu_index)
+     call mpas_pool_get_array(mesh, 'landusef', landusef)
      call mpas_pool_get_config(mgr % pool, 'scale_factor', scalefactor)
      call mpas_pool_get_config(mgr % pool, 'category_min', category_min)
      call mpas_pool_get_config(mgr % pool, 'category_max', category_max)
@@ -1808,19 +2006,39 @@ subroutine interp_landuse(mesh, kdtree, geog_data_path, isice_lu, iswater_lu, su
 
      isice_lu = isice_lu_ptr
      iswater_lu = iswater_lu_ptr
+     ismax_lu = category_max
 
      allocate(ncat(category_min:category_max, nCells))
      ncat(:,:) = 0
+     landusef(:,:) = 0.0
 
      call init_atm_map_static_data(mesh, mgr, kdtree, categorical_interp_criteria, categorical_interp_accumulation, &
                                    supersample_fac=supersample_fac)
 
-     do iCell = 1, nCells
-         ! Because maxloc returns the location of the maximum value of an array as if the
-         ! starting index of the array is 1, and dataset categories do not necessarily start
-         ! at 1, we need to use category_min to ensure the correct category location is chosen.
-         lu_index(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
-     end do
+     if (.not. config_frac_landuse) then  ! compute dominant land use
+         do iCell = 1, nCells
+             ! Because maxloc returns the location of the maximum value of an array as if the
+             ! starting index of the array is 1, and dataset categories do not necessarily start
+             ! at 1, we need to use category_min to ensure the correct category location is chosen.
+             lu_index(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
+             landusef(lu_index(iCell),iCell) = 1.0  ! defines dominant land use only
+         end do
+     else                                 ! compute fractional land use
+         do iCell = 1, nCells
+             ! Note that lu_index currently defines dominant land use only for ivgtyp, where 
+             ! ivgtyp is used by Noah LSM (and PX LSM when landusef=0 everywhere).
+             lu_index(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
+             ncattot = 0
+             do i = 1,ismax_lu
+                 ncattot = ncattot + ncat(i,iCell)
+             end do
+             do i = 1,ismax_lu
+                 if (ncat(i,iCell) > 0) then
+                     landusef(i,iCell) = real(ncat(i,iCell))/ncattot  ! defines fractional land use for PX LSM
+                 end if
+             end do
+         end do
+     end if
      deallocate(ncat)
 
      ierr = mgr % finalize()
@@ -1834,6 +2052,172 @@ subroutine interp_landuse(mesh, kdtree, geog_data_path, isice_lu, iswater_lu, su
 
  end subroutine interp_landuse
 
+ !***********************************************************************
+ !
+ !  routine interp_landuse2
+ !
+ !> \brief   Interpolate NLCD land use for CONUS
+ !> \details
+ !> Interpolate NLCD land use by using the nlcd2011_ll_9s/ dataset. Then
+ !> overlay the interpolated NLCD data onto CONUS, leaving the previously-
+ !> interpolated MODIS land use in place over the rest of the globe.
+ !
+ !-----------------------------------------------------------------------
+ subroutine interp_landuse2(mesh, dims, configs, geog_data_path, isice_lu, iswater_lu, ismax_lu)
+
+     implicit none
+
+     ! Input variables
+     type (mpas_pool_type), intent(inout) :: mesh
+     type (mpas_pool_type), intent(in) :: dims
+     type (mpas_pool_type), intent(in) :: configs
+     character (len=*), intent(in) :: geog_data_path
+     integer, intent(inout) :: isice_lu
+     integer, intent(inout) :: iswater_lu
+     integer, intent(inout) :: ismax_lu
+
+     ! Local variables
+     type (mpas_geotile_mgr_type) :: mgr
+     character (len=StrKIND) :: fname
+     character (kind=c_char), dimension(StrKIND+1) :: c_fname
+     real (kind=c_float), dimension(:,:,:), pointer, contiguous :: rarray
+     real (kind=RKIND) :: scalefactor
+     real (kind=RKIND) :: dx, dy
+     real (kind=RKIND) :: lat_pt, lon_pt
+     real (kind=RKIND) :: start_lat, start_lon
+     real (kind=RKIND), dimension(:), pointer :: latCell, lonCell
+     integer, pointer :: nCells, maxEdges
+     integer, dimension(:), pointer :: nEdgesOnCell
+     integer, dimension(:,:), pointer :: cellsOnCell
+     integer (c_int) :: endian, isigned, istatus, wordsize
+     integer (c_int) :: nx, ny, nz
+     integer:: i,j,k
+     integer:: iCell,iPoint,iTileStart,iTileEnd,jTileStart,jTileEnd,iEnd,jEnd
+     type (c_ptr) :: rarray_ptr
+
+     logical, pointer :: config_frac_landuse
+
+     integer :: ncattot
+     integer :: ierr
+
+     ierr = mgr % init(trim(geog_data_path))
+     if (ierr /= 0) then
+         call mpas_log_write("Error occured initalizing interpolation for "//trim(geog_data_path), messageType=MPAS_LOG_CRIT)
+         return ! Program execution should not reach this statement since the preceding message is a critical error
+     end if
+
+     call mpas_pool_get_config(configs, 'config_frac_landuse', config_frac_landuse)
+
+     call mpas_pool_get_dimension(dims, 'nCells', nCells)
+     call mpas_pool_get_dimension(dims, 'maxEdges', maxEdges)
+     call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
+     call mpas_pool_get_array(mesh, 'cellsOnCell', cellsOnCell)
+     call mpas_pool_get_array(mesh, 'latCell', latCell)
+     call mpas_pool_get_array(mesh, 'lonCell', lonCell)
+     call mpas_pool_get_array(mesh, 'lu_index', lu_index)
+     call mpas_pool_get_array(mesh, 'landusef', landusef)
+     call mpas_pool_get_config(mgr % pool, 'category_min', category_min)
+
+     ismax_lu = 40
+     iswater_lu = 17 
+     isice_lu = 15
+     nx = 6000
+     ny = 3800
+     nz = 40
+     isigned  = 0
+     endian   = 0
+     wordsize = 2
+     scalefactor = 0.1
+     dx = 0.0025
+     dy = 0.0025
+     start_lat = 9.00125000
+     start_lon = -139.99875000
+     jEnd = 19001
+     iEnd = 30001
+     allocate(rarray(nx,ny,nz))
+     allocate(ncat(ismax_lu,nCells))
+     ncat(:,:) = 0
+
+     rarray_ptr = c_loc(rarray)
+
+     do jTileStart = 1,jEnd,ny
+         jTileEnd = jTileStart + ny - 1
+
+         do iTileStart = 1,iEnd,nx
+             iTileEnd = iTileStart + nx - 1
+             write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(geog_data_path), &
+                   iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+             call mpas_log_write(trim(fname))
+             call mpas_f_to_c_string(fname, c_fname)
+
+             call read_geogrid(c_fname,rarray_ptr,nx,ny,nz,isigned,endian, &
+                               wordsize,istatus)
+             call init_atm_check_read_error(istatus, fname)
+             rarray(:,:,:) = rarray(:,:,:) * real(scalefactor, kind=c_float)  ! rarray contains fractional % (0-100) land use by category
+
+             iPoint = 1
+             do k=1,nz
+                 do j=1,ny
+                     do i=1,nx
+                         ! Check if the NLCD dataset has zeros
+                         if (nint(rarray(i,j,k)) == 0) cycle
+
+                         lat_pt = start_lat + (jTileStart + j - 2) * dy
+                         lon_pt = start_lon + (iTileStart + i - 2) * dx
+                         lat_pt = lat_pt * PI / 180.0
+                         lon_pt = lon_pt * PI / 180.0
+
+                         iPoint = nearest_cell(lat_pt,lon_pt,iPoint,nCells,maxEdges, &
+                                               nEdgesOnCell,cellsOnCell, &
+                                               latCell,lonCell)
+                         ncat(k,iPoint) = ncat(k,iPoint) + 1
+                     end do
+                 end do
+             end do
+
+         end do
+     end do
+
+     if (.not. config_frac_landuse) then  ! compute dominant NLCD land use for CONUS
+         do iCell = 1, nCells
+             if (maxloc(ncat(:,iCell), dim=1) < 21) cycle  ! do not change MODIS land use outside CONUS
+             ! Because maxloc returns the location of the maximum value of an array as if the
+             ! starting index of the array is 1, and dataset categories do not necessarily start
+             ! at 1, we need to use category_min to ensure the correct category location is chosen.
+             lu_index(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
+             landusef(lu_index(iCell),iCell) = 1.0  ! defines dominant NLCD land use only
+         end do
+     else                                 ! compute fractional NLCD land use for CONUS
+         do iCell = 1, nCells
+             if (maxloc(ncat(:,iCell), dim=1) < 21) cycle  ! do not change MODIS land use outide CONUS
+             ! Note that lu_index currently defines dominant land use only for ivgtyp, where 
+             ! ivgtyp is used by Noah LSM (and PX LSM when landusef=0 everywhere).
+             lu_index(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
+             ncattot = 0
+             do i = 1,ismax_lu
+                 ncattot = ncattot + ncat(i,iCell)
+             end do
+             do i = 1,ismax_lu
+                 if (ncat(i,iCell) > 0) then
+                     do j = 1,ismax_lu
+                         landusef(j,iCell) = 0.0 ! clear MODIS landusef if NLCD (CONUS) cell
+                     end do
+                 end if
+             end do
+             do i = 1,ismax_lu
+                 if (ncat(i,iCell) > 0) then
+                     landusef(i,iCell) = real(ncat(i,iCell))/ncattot  ! defines fractional land use for PX LSM
+                 end if
+             end do
+         end do
+     end if
+     deallocate(rarray)
+     deallocate(ncat)
+
+     nullify(category_min)
+
+ end subroutine interp_landuse2
+
  !***********************************************************************
  !
  !  routine interp_soilcat
@@ -1852,13 +2236,14 @@ end subroutine interp_landuse
  !> iswater is present in the dataset's index file.
  !>
  !-----------------------------------------------------------------------
- subroutine interp_soilcat(mesh, kdtree, geog_data_path, iswater_soil, supersample_fac)
+ subroutine interp_soilcat(configs, mesh, kdtree, geog_data_path, iswater_soil, supersample_fac)
 
      implicit none
 
      ! Input variables
      type (mpas_pool_type), intent(inout) :: mesh
      type (mpas_kd_type), pointer, intent(in) :: kdtree
+     type (mpas_pool_type), intent(in) :: configs
      character (len=*), intent(in) :: geog_data_path
      integer, intent(out) :: iswater_soil
      integer, intent(in), optional :: supersample_fac
@@ -1870,7 +2255,10 @@ subroutine interp_soilcat(mesh, kdtree, geog_data_path, iswater_soil, supersampl
 
      real (kind=RKIND), pointer :: scalefactor
 
-     integer :: iCell
+     logical, pointer :: config_frac_landuse
+
+     integer :: ncattot
+     integer :: i, iCell
      integer :: ierr
 
      ierr = mgr % init(trim(geog_data_path))
@@ -1879,8 +2267,11 @@ subroutine interp_soilcat(mesh, kdtree, geog_data_path, iswater_soil, supersampl
          return
      end if
 
+     call mpas_pool_get_config(configs, 'config_frac_landuse', config_frac_landuse)
+
      call mpas_pool_get_dimension(mesh, 'nCells', nCells)
-     call mpas_pool_get_array(mesh, 'soilcat_top', soilcat_top)
+     call mpas_pool_get_array(mesh, 'soilcat', soilcat)
+     call mpas_pool_get_array(mesh, 'soiltypf', soiltypf)
      call mpas_pool_get_config(mgr % pool, 'scale_factor', scalefactor)
      call mpas_pool_get_config(mgr % pool, 'category_min', category_min)
      call mpas_pool_get_config(mgr % pool, 'category_max', category_max)
@@ -1890,16 +2281,40 @@ subroutine interp_soilcat(mesh, kdtree, geog_data_path, iswater_soil, supersampl
 
      allocate(ncat(category_min:category_max, nCells))
      ncat(:,:) = 0
+     soiltypf(:,:) = 0.0
 
      call init_atm_map_static_data(mesh, mgr, kdtree, categorical_interp_criteria, categorical_interp_accumulation, &
                                    supersample_fac=supersample_fac)
 
-     do iCell = 1, nCells
-        ! Because maxloc returns the location of the maximum value of an array as if the
-        ! starting index of the array is 1, and dataset categories do not necessarily start
-        ! at 1, we need to use category_min to ensure the correct category location is chosen.
-        soilcat_top(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
-     end do
+     if(.not. config_frac_landuse) then  ! compute dominant soil type for dominant land use
+         do iCell = 1, nCells
+             ! Because maxloc returns the location of the maximum value of an array as if the
+             ! starting index of the array is 1, and dataset categories do not necessarily start
+             ! at 1, we need to use category_min to ensure the correct category location is chosen.
+             soilcat(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min
+             soiltypf(soilcat(iCell),iCell) = 1.0  ! defines dominant soil type only
+         end do
+     else                                ! compute fractional soil type for fractional land use
+         do iCell = 1,nCells
+             soilcat(iCell) = 1
+             do i = 2,16
+                 if(ncat(i,iCell) > ncat(soilcat(iCell),iCell)) then
+                     soilcat(iCell) = i      ! currently defines dominant soil type only for isltyp
+                 end if
+             end do
+         end do
+         do iCell = 1,nCells
+             ncattot = 0
+             do i = 1,16
+                 ncattot = ncattot + ncat(i,iCell)
+             end do
+             do i = 1,16
+                 if(ncat(i,iCell) > 0) then
+                     soiltypf(i,iCell) = real(ncat(i,iCell))/ncattot  ! defines fractional soil type for PX LSM
+                 end if
+             end do
+         end do
+     end if
      deallocate(ncat)
 
      ierr = mgr % finalize()