diff --git a/components/elm/bld/namelist_files/namelist_defaults.xml b/components/elm/bld/namelist_files/namelist_defaults.xml
index 4f02fc9df8ac..52842251684c 100644
--- a/components/elm/bld/namelist_files/namelist_defaults.xml
+++ b/components/elm/bld/namelist_files/namelist_defaults.xml
@@ -2155,4 +2155,7 @@ this mask will have smb calculated over the entire global land surface
.false.
.false.
+
+.false.;
+
diff --git a/components/elm/bld/namelist_files/namelist_definition.xml b/components/elm/bld/namelist_files/namelist_definition.xml
index 1353f8b84b64..22843eddc634 100644
--- a/components/elm/bld/namelist_files/namelist_definition.xml
+++ b/components/elm/bld/namelist_files/namelist_definition.xml
@@ -2135,4 +2135,14 @@ not be called.
Default: 0
+
+
+
+
+
+Use polygonal tundra parameterizations for ice wedge polygon microtopography and inindation fraction
+Default: .false.
+
+
diff --git a/components/elm/src/biogeophys/ActiveLayerMod.F90 b/components/elm/src/biogeophys/ActiveLayerMod.F90
index 1526b6ecb342..c6133d37fc7f 100644
--- a/components/elm/src/biogeophys/ActiveLayerMod.F90
+++ b/components/elm/src/biogeophys/ActiveLayerMod.F90
@@ -71,16 +71,21 @@ subroutine alt_calc(num_soilc, filter_soilc, &
real(r8), dimension(nlevgrnd) :: melt_profile ! profile of melted excess ice
!-----------------------------------------------------------------------
+ ! RF NOTE: use of 1989 ALT in these parameterizations is somewhat of a placeholder used to compare against
+ ! ATS runs for NGEE Arctic Phase 3. It should ultimately be replaced by a more physically based threshold
+ ! later on.
associate( &
t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd)
alt => canopystate_vars%alt_col , & ! Output: [real(r8) (:) ] current depth of thaw
altmax => canopystate_vars%altmax_col , & ! Output: [real(r8) (:) ] maximum annual depth of thaw
altmax_lastyear => canopystate_vars%altmax_lastyear_col , & ! Output: [real(r8) (:) ] prior year maximum annual depth of thaw
+ altmax_1989 => canopystate_vars%altmax_1989_col , & ! Output: [real(r8) (:) ] maximum ALT in 1989
altmax_ever => canopystate_vars%altmax_ever_col , & ! Output: [real(r8) (:) ] maximum thaw depth since initialization
alt_indx => canopystate_vars%alt_indx_col , & ! Output: [integer (:) ] current depth of thaw
altmax_indx => canopystate_vars%altmax_indx_col , & ! Output: [integer (:) ] maximum annual depth of thaw
altmax_lastyear_indx => canopystate_vars%altmax_lastyear_indx_col , & ! Output: [integer (:) ] prior year maximum annual depth of thaw
+ altmax_1989_indx => canopystate_vars%altmax_1989_indx_col, & ! Output: [integer (:) ] index of maximum ALT in 1989
altmax_ever_indx => canopystate_vars%altmax_ever_indx_col, & ! Output: [integer (:) ] maximum thaw depth since initialization
excess_ice => col_ws%excess_ice , & ! Input: [real(r8) (:,:) ] depth variable excess ice content in soil column (-)
rmax => col_pp%iwp_microrel , & ! Output: [real(r8) (:) ] ice wedge polygon microtopographic relief (m)
@@ -154,7 +159,6 @@ subroutine alt_calc(num_soilc, filter_soilc, &
endif
endif
-
! if appropriate, update maximum annual active layer thickness
if (alt(c) > altmax(c)) then
altmax(c) = alt(c)
@@ -170,6 +174,13 @@ subroutine alt_calc(num_soilc, filter_soilc, &
endif
endif
+ ! special loop for if year = 1989, see above note regarding
+ ! replacing with more physically based mechanism.
+ if (year .eq. 1989) then
+ altmax_1989(c) = altmax(c)
+ altmax_1989_indx(c) = altmax_indx(c)
+ endif
+
! update subsidence based on change in ALT
! melt_profile stores the amount of excess_ice
! melted in this timestep.
@@ -178,15 +189,22 @@ subroutine alt_calc(num_soilc, filter_soilc, &
melt_profile(j) = 0.0_r8
else if (j .eq. k_frz) then
melt_profile(j) = excess_ice(c,j) + ((z2-alt(c))/(z2-z1))*excess_ice(c,j+1) ! TODO: check indices here!!
+ ! remove melted excess ice:
+ excess_ice(c,j) = 0._r8
+ excess_ice(c,j+1) = excess_ice(c,j+1)*(1._r8 - min(1._r8,(z2-alt(c))/(z2-z1)))
else
melt_profile(j) = excess_ice(c,j)
+ ! remove melted excess ice
+ excess_ice(c,j) = 0._r8
end if
! calculate subsidence at this layer:
melt_profile(j) = melt_profile(j) * dzsoi(j)
end do
! subsidence is integral of melt profile:
- subsidence(c) = subsidence(c) + sum(melt_profile)
+ if ((year .ge. 1989) .and. (altmax_ever(c) .ge. altmax_1989(c))) then
+ subsidence(c) = subsidence(c) + sum(melt_profile)
+ end if
! limit subsidence to 0.4 m
subsidence(c) = min(0.4_r8, subsidence(c))
diff --git a/components/elm/src/biogeophys/CanopyStateType.F90 b/components/elm/src/biogeophys/CanopyStateType.F90
index 3ce62408e999..803366c05000 100644
--- a/components/elm/src/biogeophys/CanopyStateType.F90
+++ b/components/elm/src/biogeophys/CanopyStateType.F90
@@ -10,6 +10,7 @@ module CanopyStateType
use elm_varcon , only : spval,ispval
use elm_varpar , only : nlevcan, nvegwcs
use elm_varctl , only : iulog, use_cn, use_fates, use_hydrstress, use_fates_sp
+ use elm_varctl , only : use_polygonal_tundra
use LandunitType , only : lun_pp
use ColumnType , only : col_pp
use VegetationType , only : veg_pp
@@ -62,9 +63,11 @@ module CanopyStateType
integer , pointer :: alt_indx_col (:) ! col current depth of thaw
real(r8) , pointer :: altmax_col (:) ! col maximum annual depth of thaw
real(r8) , pointer :: altmax_lastyear_col (:) ! col prior year maximum annual depth of thaw
- real(r8) , pointer :: altmax_ever_col (:) ! col maximum thaw depth from beginning of simulation
integer , pointer :: altmax_indx_col (:) ! col maximum annual depth of thaw
integer , pointer :: altmax_lastyear_indx_col (:) ! col prior year maximum annual depth of thaw
+ real(r8) , pointer :: altmax_1989_col (:) ! col maximum annual thaw depth in 1989 RF note: temporary for IM1!
+ real(r8) , pointer :: altmax_ever_col (:) ! col maximum thaw depth from beginning of simulation
+ integer , pointer :: altmax_1989_indx_col (:) ! col 1989 maximum thaw depth RF note: temporary for IM1!
integer , pointer :: altmax_ever_indx_col (:) ! col maximum thaw depth from beginning of simulation
real(r8) , pointer :: dewmx_patch (:) ! patch maximum allowed dew [mm]
@@ -144,11 +147,13 @@ subroutine InitAllocate(this, bounds )
allocate(this%alt_col (begc:endc)) ; this%alt_col (:) = spval;
allocate(this%altmax_col (begc:endc)) ; this%altmax_col (:) = spval
allocate(this%altmax_lastyear_col (begc:endc)) ; this%altmax_lastyear_col (:) = spval
+ allocate(this%altmax_1989_col (begc:endc)) ; this%altmax_1989_col (:) = spval
allocate(this%altmax_ever_col (begc:endc)) ; this%altmax_ever_col (:) = spval
allocate(this%alt_indx_col (begc:endc)) ; this%alt_indx_col (:) = huge(1)
allocate(this%altmax_indx_col (begc:endc)) ; this%altmax_indx_col (:) = huge(1)
allocate(this%altmax_lastyear_indx_col (begc:endc)) ; this%altmax_lastyear_indx_col (:) = huge(1)
- allocate(this%altmax_ever_indx_col (begc:endc)) ; this%altmax_ever_indx_col (:) = huge(1)
+ allocate(this%altmax_1989_indx_col (begc:endc)) ; this%altmax_1989_indx_col (:) = huge(1)
+ allocate(this%altmax_ever_indx_col (begc:endc)) ; this%altmax_ever_indx_col (:) = huge(1)
allocate(this%dewmx_patch (begp:endp)) ; this%dewmx_patch (:) = spval
allocate(this%dleaf_patch (begp:endp)) ; this%dleaf_patch (:) = spval
@@ -282,10 +287,17 @@ subroutine InitHistory(this, bounds)
avgflag='A', long_name='maximum prior year active layer thickness', &
ptr_col=this%altmax_lastyear_col)
- this%altmax_ever_col(begc:endc) = spval
- call hist_addfld1d (fname='ALTMAX_EVER', units='m', &
- avgflag='A', long_name='maximum ever active layer thickness throughout simulation', &
- ptr_col=this%altmax_ever_col)
+ if (use_polygonal_tundra) then
+ this%altmax_1989_col(begc:endc) = spval
+ call hist_addfld1d (fname='ALTMAX_1989', units='m', &
+ avgflag='A', long_name='maximum 1989 active layer thickness', &
+ ptr_col=this%altmax_1989_col)
+
+ this%altmax_ever_col(begc:endc) = spval
+ call hist_addfld1d (fname='ALTMAX_EVER', units='m', &
+ avgflag='A', long_name='maximum ever active layer thickness throughout simulation', &
+ ptr_col=this%altmax_ever_col)
+ end if
end if
! Allow active layer fields to be optionally output even if not running CN
@@ -306,10 +318,17 @@ subroutine InitHistory(this, bounds)
avgflag='A', long_name='maximum prior year active layer thickness', &
ptr_col=this%altmax_lastyear_col, default='inactive')
- this%altmax_ever_col(begc:endc) = spval
- call hist_addfld1d (fname='ALTMAX_EVER', units='m', &
- avgflag='A', long_name='maximum ever active layer thickness throughout simulation', &
- ptr_col=this%altmax_ever_col, default='inactive')
+ if (use_polygonal_tundra) then
+ this%altmax_1989_col(begc:endc) = spval
+ call hist_addfld1d (fname='ALTMAX_1989', units='m', &
+ avgflag='A', long_name='maximum 1989 active layer thickness', &
+ ptr_col=this%altmax_1989_col, default='inactive')
+
+ this%altmax_ever_col(begc:endc) = spval
+ call hist_addfld1d (fname='ALTMAX_EVER', units='m', &
+ avgflag='A', long_name='maximum ever active layer thickness throughout simulation', &
+ ptr_col=this%altmax_ever_col, default='inactive')
+ end if
end if
! Accumulated fields
@@ -516,10 +535,12 @@ subroutine InitCold(this, bounds)
this%alt_col(c) = 0._r8
this%altmax_col(c) = 0._r8
this%altmax_lastyear_col(c) = 0._r8
+ this%altmax_1989_col(c) = 0._r8
this%altmax_ever_col(c) = 0._r8
this%alt_indx_col(c) = 0
this%altmax_indx_col(c) = 0
this%altmax_lastyear_indx_col(c) = 0
+ this%altmax_1989_indx_col(c) = 0
this%altmax_ever_indx_col(c) = 0
end if
end do
@@ -588,19 +609,27 @@ subroutine Restart(this, bounds, ncid, flag)
call restartvar(ncid=ncid, flag=flag, varname='altmax_lastyear', xtype=ncd_double, &
dim1name='column', long_name='', units='', &
interpinic_flag='interp', readvar=readvar, data=this%altmax_lastyear_col)
- call restartvar(ncid=ncid, flag=flag, varname='altmax_ever', xtype=ncd_double, &
- dim1name='column', long_name='', units='', &
- interpinic_flag='interp', readvar=readvar, data=this%altmax_ever_col)
call restartvar(ncid=ncid, flag=flag, varname='altmax_indx', xtype=ncd_int, &
dim1name='column', long_name='', units='', &
interpinic_flag='interp', readvar=readvar, data=this%altmax_indx_col)
call restartvar(ncid=ncid, flag=flag, varname='altmax_lastyear_indx', xtype=ncd_int, &
dim1name='column', long_name='', units='', &
interpinic_flag='interp', readvar=readvar, data=this%altmax_lastyear_indx_col)
- call restartvar(ncid=ncid, flag=flag, varname='altmax_ever_indx', xtype=ncd_int, &
- dim1name='column', long_name='', units='', &
- interpinic_flag='interp', readvar=readvar, data=this%altmax_ever_indx_col)
- end if
+ if (use_polygonal_tundra) then
+ call restartvar(ncid=ncid, flag=flag, varname='altmax_1989_indx', xtype=ncd_double, &
+ dim1name='column', long_name='', units='', &
+ interpinic_flag='interp', readvar=readvar, data=this%altmax_1989_indx_col)
+ call restartvar(ncid=ncid, flag=flag, varname='altmax_ever_indx', xtype=ncd_int, &
+ dim1name='column', long_name='', units='', &
+ interpinic_flag='interp', readvar=readvar, data=this%altmax_ever_indx_col)
+ call restartvar(ncid=ncid, flag=flag, varname='altmax_1989', xtype=ncd_double, &
+ dim1name='column', long_name='', units='', &
+ interpinic_flag='interp', readvar=readvar, data=this%altmax_1989_col)
+ call restartvar(ncid=ncid, flag=flag, varname='altmax_ever', xtype=ncd_double, &
+ dim1name='column', long_name='', units='', &
+ interpinic_flag='interp', readvar=readvar, data=this%altmax_ever_col)
+ end if ! polygonal tundra
+ end if ! cn or fates
if ( use_hydrstress ) then
call restartvar(ncid=ncid, flag=flag, varname='vegwp', xtype=ncd_double, &
dim1name='pft', dim2name='vegwcs', switchdim=.true., &
diff --git a/components/elm/src/biogeophys/SoilHydrologyMod.F90 b/components/elm/src/biogeophys/SoilHydrologyMod.F90
index 0a7e76cd4a26..d235223ff348 100644
--- a/components/elm/src/biogeophys/SoilHydrologyMod.F90
+++ b/components/elm/src/biogeophys/SoilHydrologyMod.F90
@@ -272,7 +272,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
use elm_varcon , only : denh2o, denice, roverg, wimp, mu, tfrz
use elm_varcon , only : pondmx, watmin
use column_varcon , only : icol_roof, icol_road_imperv, icol_sunwall, icol_shadewall, icol_road_perv
- use landunit_varcon , only : istsoil, istcrop, ilowcenpoly
+ use landunit_varcon , only : istsoil, istcrop, ilowcenpoly, iflatcenpoly, ihighcenpoly
use elm_time_manager , only : get_step_size, get_nstep
use atm2lndType , only : atm2lnd_type ! land river two way coupling
use lnd2atmType , only : lnd2atm_type
diff --git a/components/elm/src/main/controlMod.F90 b/components/elm/src/main/controlMod.F90
index 51e0b2898417..b73eac7b9171 100755
--- a/components/elm/src/main/controlMod.F90
+++ b/components/elm/src/main/controlMod.F90
@@ -60,12 +60,12 @@ module controlMod
! !PUBLIC TYPES:
implicit none
save
-
+
! !PUBLIC MEMBER FUNCTIONS:
public :: control_setNL ! Set namelist filename
public :: control_init ! initial run control information
public :: control_print ! print run control information
-
+
! !PRIVATE TYPES:
character(len= 7) :: runtyp(4) ! run type
character(len=SHR_KIND_CL) :: NLFilename = 'lnd.stdin' ! Namelist filename
@@ -79,14 +79,14 @@ module controlMod
!------------------------------------------------------------------------
subroutine control_setNL( NLfile )
-
+
! !DESCRIPTION:
! Set the namelist filename to use
-
+
! !ARGUMENTS:
implicit none
character(len=*), intent(IN) :: NLFile ! Namelist filename
-
+
! !LOCAL VARIABLES:
character(len=32) :: subname = 'control_setNL' ! subroutine name
logical :: lexist ! File exists
@@ -111,19 +111,19 @@ end subroutine control_setNL
!------------------------------------------------------------------------
subroutine control_init( )
-
+
! !DESCRIPTION:
! Initialize CLM run control information
-
+
! !USES:
use elm_time_manager , only : set_timemgr_init, get_timemgr_defaults
use fileutils , only : getavu, relavu
use shr_string_mod , only : shr_string_getParentDir
use elm_interface_pflotranMod , only : elm_pf_readnl
use ELMBeTRNLMod , only : betr_readNL
-
+
implicit none
-
+
! !LOCAL VARIABLES:
character(len=32) :: starttype ! infodata start type
integer :: i,j,n ! loop indices
@@ -306,8 +306,8 @@ subroutine control_init( )
namelist /elm_inparm/ &
do_budgets, budget_inst, budget_daily, budget_month, &
budget_ann, budget_ltann, budget_ltend
-
- namelist /elm_inparm/ &
+
+ namelist /elm_inparm/ &
use_atm_downscaling_to_topunit, precip_downscaling_method
namelist /elm_inparm/ &
@@ -318,16 +318,19 @@ subroutine control_init( )
namelist /elm_mosart/ &
lnd_rof_coupling_nstep
-
+
+ namelist /elm_inparm/ &
+ snow_shape, snicar_atm_type, use_dust_snow_internal_mixing
+
namelist /elm_inparm/ &
- snow_shape, snicar_atm_type, use_dust_snow_internal_mixing
-
- namelist /elm_inparm/ &
use_modified_infil
namelist /elm_inparm/ &
use_fan, fan_mode, fan_to_bgc_veg, nh4_ads_coef
+ ! NGEE Arctic options
+ namelist /elm_inparm/ &
+ use_polygonal_tundra
! ----------------------------------------------------------------------
! Default values
! ----------------------------------------------------------------------
@@ -506,7 +509,7 @@ subroutine control_init( )
end if
end if
- if (use_lch4 .and. use_vertsoilc) then
+ if (use_lch4 .and. use_vertsoilc) then
anoxia = .true.
else
anoxia = .false.
@@ -551,7 +554,7 @@ subroutine control_init( )
if (use_lnd_rof_two_way) then
if (lnd_rof_coupling_nstep < 1) then
call endrun(msg=' ERROR: lnd_rof_coupling_nstep cannot be smaller than 1.'//&
- errMsg(__FILE__, __LINE__))
+ errMsg(__FILE__, __LINE__))
endif
endif
@@ -686,7 +689,7 @@ end subroutine control_init
!------------------------------------------------------------------------
subroutine control_spmd()
-
+
! !DESCRIPTION:
! Distribute namelist data all processors. All program i/o is
! funnelled through the master processor. Processor 0 either
@@ -695,13 +698,13 @@ subroutine control_spmd()
! all processors and writes it to disk.
!
! !USES:
-
+
use spmdMod, only : mpicom, MPI_CHARACTER, MPI_INTEGER, MPI_LOGICAL, MPI_REAL8
use elm_varpar, only : numrad
-
+
! !ARGUMENTS:
implicit none
-
+
! !LOCAL VARIABLES:
integer ier !error code
!-----------------------------------------------------------------------
@@ -874,7 +877,7 @@ subroutine control_spmd()
call mpi_bcast (more_vertlayers,1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (const_climate_hist, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_top_solar_rad, 1, MPI_LOGICAL, 0, mpicom, ier) ! TOP solar radiation parameterization
-
+
! glacier_mec variables
call mpi_bcast (create_glacier_mec_landunit, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (maxpatch_glcmec, 1, MPI_INTEGER, 0, mpicom, ier)
@@ -948,11 +951,11 @@ subroutine control_spmd()
! PETSc-based thermal model
call mpi_bcast (use_petsc_thermal_model, 1, MPI_LOGICAL, 0, mpicom, ier)
-
+
! Downscaling of atmospheric forcing to topounits
call mpi_bcast (use_atm_downscaling_to_topunit, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (precip_downscaling_method, len(precip_downscaling_method), MPI_CHARACTER, 0, mpicom, ier)
-
+
! soil erosion
call mpi_bcast (use_erosion, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (ero_ccycle , 1, MPI_LOGICAL, 0, mpicom, ier)
@@ -974,28 +977,30 @@ subroutine control_spmd()
call mpi_bcast (snow_shape, len(snow_shape), MPI_CHARACTER, 0, mpicom, ier)
call mpi_bcast (snicar_atm_type, len(snicar_atm_type), MPI_CHARACTER, 0, mpicom, ier)
call mpi_bcast (use_dust_snow_internal_mixing, 1, MPI_LOGICAL, 0, mpicom, ier)
-
+
call mpi_bcast (mpi_sync_nstep_freq, 1, MPI_INTEGER, 0, mpicom, ier)
-
+
! use modified infiltration scheme in surface water storage
call mpi_bcast (use_modified_infil, 1, MPI_LOGICAL, 0, mpicom, ier)
+ !NGEE Arctic options
+ call mpi_bcast (use_polygonal_tundra, 1, MPI_LOGICAL, 0, mpicom, ier)
end subroutine control_spmd
!------------------------------------------------------------------------
subroutine control_print ()
-
+
! !DESCRIPTION:
! Write out the clm namelist run control variables
-
+
! !USES:
-
+
use AllocationMod, only : suplnitro, suplnNon
use AllocationMod, only : suplphos, suplpNon
-
+
! !ARGUMENTS:
implicit none
-
+
! !LOCAL VARIABLES:
integer i !loop index
character(len=32) :: subname = 'control_print' ! subroutine name
@@ -1035,7 +1040,7 @@ subroutine control_print ()
write(iulog,*) ' precip_downscaling_method = ', precip_downscaling_method
write(iulog,*) 'input data files:'
write(iulog,*) ' PFT physiology and parameters file = ',trim(paramfile)
- write(iulog,*) ' Soil order dependent parameters file = ',trim(fsoilordercon)
+ write(iulog,*) ' Soil order dependent parameters file = ',trim(fsoilordercon)
if (fsurdat == ' ') then
write(iulog,*) ' fsurdat, surface dataset not set'
else
@@ -1056,13 +1061,13 @@ subroutine control_print ()
else
write(iulog,*) ' atm topographic data = ',trim(fatmtopo)
end if
-
+
if (use_top_solar_rad) then
write(iulog,*) ' use TOP solar radiation parameterization instead of PP'
else
write(iulog,*) ' use_top_solar_rad is False, so do not run TOP solar radiation parameterization'
end if
-
+
if (use_cn) then
if (suplnitro /= suplnNon)then
write(iulog,*) ' Supplemental Nitrogen mode is set to run over Patches: ', &
@@ -1176,9 +1181,9 @@ subroutine control_print ()
write(iulog,*) ' atm_gustiness = ', atm_gustiness
write(iulog,*) ' force_land_gustiness = ', force_land_gustiness
write(iulog,*) ' more vertical layers = ', more_vertlayers
-
+
write(iulog,*) ' Sub-grid topographic effects on solar radiation = ', use_top_solar_rad ! TOP solar radiation parameterization
-
+
if (nsrest == nsrContinue) then
write(iulog,*) 'restart warning:'
write(iulog,*) ' Namelist not checked for agreement with initial run.'
@@ -1241,9 +1246,9 @@ subroutine control_print ()
write(iulog,*) ' use_lnd_rof_two_way = ', use_lnd_rof_two_way
write(iulog,*) ' lnd_rof_coupling_nstep = ', lnd_rof_coupling_nstep
write(iulog,*) ' mpi_sync_nstep_freq = ', mpi_sync_nstep_freq
-
+
write(iulog,*) ' use_modified_infil = ', use_modified_infil
-
+
! FAN
write(iulog,*) ' use_fan = ', use_fan
@@ -1253,6 +1258,8 @@ subroutine control_print ()
write(iulog,*) ' fan_to_bgc_veg = ', fan_to_bgc_veg
end if
+ ! NGEE Arctic options
+ if (use_polygonal_tundra) write(iulog, *) ' use_polygonal_tundra =', use_polygonal_tundra
end subroutine control_print
end module controlMod
diff --git a/components/elm/src/main/elm_initializeMod.F90 b/components/elm/src/main/elm_initializeMod.F90
index 2a8079692379..e131a1938076 100755
--- a/components/elm/src/main/elm_initializeMod.F90
+++ b/components/elm/src/main/elm_initializeMod.F90
@@ -14,7 +14,7 @@ module elm_initializeMod
use elm_varctl , only : use_lch4, use_cn, use_voc, use_c13, use_c14
use elm_varctl , only : use_fates, use_betr, use_fates_sp, use_fan
use elm_varsur , only : wt_lunit, urban_valid, wt_nat_patch, wt_cft, wt_glc_mec, topo_glc_mec,firrig,f_surf,f_grd
- use elm_varsur , only : fert_cft, fert_p_cft
+ use elm_varsur , only : fert_cft, fert_p_cft, wt_polygon
use elm_varsur , only : wt_tunit, elv_tunit, slp_tunit,asp_tunit,num_tunit_per_grd
use perf_mod , only : t_startf, t_stopf
!use readParamsMod , only : readParameters
@@ -62,7 +62,7 @@ subroutine initialize1( )
use elm_varpar , only: update_pft_array_bounds
use elm_varpar , only: surfpft_lb, surfpft_ub
use elm_varcon , only: elm_varcon_init
- use landunit_varcon , only: landunit_varcon_init, max_lunit, istice_mec
+ use landunit_varcon , only: landunit_varcon_init, max_lunit, istice_mec, max_polygon
use column_varcon , only: col_itype_to_icemec_class
use elm_varctl , only: fsurdat, fatmlndfrc, flndtopo, fglcmask, noland, version
use pftvarcon , only: pftconrd
@@ -86,7 +86,7 @@ subroutine initialize1( )
use filterMod , only: allocFilters
use reweightMod , only: reweight_wrapup
use topounit_varcon , only: max_topounits, has_topounit, topounit_varcon_init
- use elm_varctl , only: use_top_solar_rad
+ use elm_varctl , only: use_top_solar_rad, use_polygonal_tundra
!
! !LOCAL VARIABLES:
integer :: ier ! error status
@@ -279,6 +279,11 @@ subroutine initialize1( )
allocate (wt_glc_mec (1,1,1))
allocate (topo_glc_mec(1,1,1))
endif
+ if (use_polygonal_tundra) then
+ allocate (wt_polygon (begg:endg,1:max_topounits, max_polygon))
+ else
+ allocate (wt_polygon (1,1,1)) ! RF-not sure this is needed
+ endif
allocate (wt_tunit (begg:endg,1:max_topounits ))
allocate (elv_tunit (begg:endg,1:max_topounits ))
@@ -421,6 +426,7 @@ subroutine initialize1( )
!deallocate (wt_lunit, wt_cft, wt_glc_mec)
deallocate (wt_cft, wt_glc_mec) !wt_lunit not deallocated because it is being used in CanopyHydrologyMod.F90
deallocate (wt_tunit, elv_tunit, slp_tunit, asp_tunit,num_tunit_per_grd)
+ deallocate (wt_polygon) ! RF - might be used elsewhere, not sure if we want to deallocate here.
call t_stopf('elm_init1')
! initialize glc_topo
diff --git a/components/elm/src/main/elm_varctl.F90 b/components/elm/src/main/elm_varctl.F90
index e1d26775a5dc..318ca506f313 100644
--- a/components/elm/src/main/elm_varctl.F90
+++ b/components/elm/src/main/elm_varctl.F90
@@ -374,7 +374,11 @@ module elm_varctl
logical, public :: fan_nh3_to_atm = .false.
logical, public :: fan_to_bgc_crop = .false.
logical, public :: fan_to_bgc_veg = .false.
-
+
+ !----------------------------------------------------------
+ ! NGEE Arctic parameterizations
+ !----------------------------------------------------------
+ logical, public :: use_polygonal_tundra = .false.
!----------------------------------------------------------
! VSFM switches
@@ -524,8 +528,8 @@ module elm_varctl
!----------------------------------------------------------
logical, public :: use_lnd_rof_two_way = .false.
integer, public :: lnd_rof_coupling_nstep = 0
-
-
+
+
!----------------------------------------------------------
! SNICAR-AD
!----------------------------------------------------------
@@ -537,7 +541,7 @@ module elm_varctl
! MPI syncing
!----------------------------------------------------------
integer, public :: mpi_sync_nstep_freq = 0
-
+
!----------------------------------------------------------
! Modified infiltration scheme in surface water storage
!----------------------------------------------------------
diff --git a/components/elm/src/main/elm_varsur.F90 b/components/elm/src/main/elm_varsur.F90
index 3423febed108..e0c846e69005 100755
--- a/components/elm/src/main/elm_varsur.F90
+++ b/components/elm/src/main/elm_varsur.F90
@@ -2,7 +2,7 @@ module elm_varsur
!-----------------------------------------------------------------------
! Module containing 2-d surface boundary data information
- ! surface boundary data, these are all "gdc" local
+ ! surface boundary data, these are all "gdc" local
! Note that some of these need to be pointers (as opposed to just allocatable arrays) to
! match the ncd_io interface; for consistency, we make them all pointers
!
@@ -14,7 +14,7 @@ module elm_varsur
save
!
! weight of each landunit on the topounit in the grid cell
- real(r8), pointer :: wt_lunit(:,:,:)
+ real(r8), pointer :: wt_lunit(:,:,:)
! whether we have valid urban data in each grid cell or (in each topounit)
logical , pointer :: urban_valid(:,:)
@@ -22,12 +22,18 @@ module elm_varsur
! for natural veg landunit, weight of each pft on the landunit (adds to 1.0 on the
! landunit for all all grid cells, even! those without any natural pft)
! (second dimension goes natpft_lb:natpft_ub)
- real(r8), pointer :: wt_nat_patch(:,:,:)
+ real(r8), pointer :: wt_nat_patch(:,:,:)
+
+ ! for polygonal tundra special case of natural veg landunit, weight of polygon type
+ ! with respect to the natural vegetation landunit (NOT the grid cell). these do
+ ! not need to sum to 1 to preserve possibility of polygonal tundra and non-polygonal
+ ! ground on the same grid cell, on the vegetated landunit.
+ real(r8), pointer :: wt_polygon(:,:,:) ! dims: clump, topounit, polygon type
! for crop landunit, weight of each cft on the landunit (adds to 1.0 on the
! landunit for all all grid cells, even those without any crop)
! (second dimension goes cft_lb:cft_ub)
- real(r8), pointer :: wt_cft(:,:,:)
+ real(r8), pointer :: wt_cft(:,:,:)
! for each cft on the crop landunit prescribe annual fertilizer
! landunit for all all grid cells, even those without any crop)
@@ -37,19 +43,19 @@ module elm_varsur
! for glc_mec landunits, weight of glacier in each elevation class (adds to 1.0 on the
! landunit for all grid cells, even those without any glacier)
- real(r8), pointer :: wt_glc_mec(:,:,:)
+ real(r8), pointer :: wt_glc_mec(:,:,:)
! subgrid glacier_mec sfc elevation
- real(r8), pointer :: topo_glc_mec(:,:,:)
-
+ real(r8), pointer :: topo_glc_mec(:,:,:)
+
! Topounit related poiters
real(r8), pointer :: num_tunit_per_grd(:) ! Topounit area fraction
real(r8), pointer :: wt_tunit(:,:) ! Topounit area fraction
- real(r8), pointer :: elv_tunit(:,:) ! Topounit elevation
- real(r8), pointer :: slp_tunit(:,:) ! Topounit slope
+ real(r8), pointer :: elv_tunit(:,:) ! Topounit elevation
+ real(r8), pointer :: slp_tunit(:,:) ! Topounit slope
integer, pointer :: asp_tunit(:,:) ! Topounit aspect
-
- real(r8),pointer :: firrig(:,:)
+
+ real(r8),pointer :: firrig(:,:)
real(r8),pointer :: f_surf(:,:) ! fraction of water withdraws from surfacewater
real(r8),pointer :: f_grd(:,:) ! fraction of water withdraws from groundwater
!-----------------------------------------------------------------------
diff --git a/components/elm/src/main/initSubgridMod.F90 b/components/elm/src/main/initSubgridMod.F90
index ffc15f3f4133..23810a9ada36 100644
--- a/components/elm/src/main/initSubgridMod.F90
+++ b/components/elm/src/main/initSubgridMod.F90
@@ -13,11 +13,11 @@ module initSubgridMod
use elm_varctl , only : iulog
use elm_varcon , only : namep, namec, namel, namet
use decompMod , only : bounds_type
- use GridcellType , only : grc_pp
+ use GridcellType , only : grc_pp
Use TopounitType , only : top_pp
- use LandunitType , only : lun_pp
- use ColumnType , only : col_pp
- use VegetationType , only : veg_pp
+ use LandunitType , only : lun_pp
+ use ColumnType , only : col_pp
+ use VegetationType , only : veg_pp
!
! !PUBLIC TYPES:
implicit none
@@ -29,33 +29,34 @@ module initSubgridMod
public :: elm_ptrs_check ! checks and writes out a summary of subgrid data
public :: add_topounit ! add an entry in the topounit-level arrays
public :: add_landunit ! add an entry in the landunit-level arrays
+ public :: add_polygon_landunit ! adds an entry in the landunit-level arrays for the special type of polygonal ground.
public :: add_column ! add an entry in the column-level arrays
public :: add_patch ! add an entry in the patch-level arrays
!
!-----------------------------------------------------------------------
contains
-
+
!------------------------------------------------------------------------------
subroutine elm_ptrs_compdown(bounds)
!
! !DESCRIPTION:
- ! Assumes the part of the subgrid pointing up has been set. Fills
+ ! Assumes the part of the subgrid pointing up has been set. Fills
! in the data pointing down. Up is p_c, p_l, p_g, c_l, c_g, and l_g.
!
! This algorithm assumes all indices besides grid cell are monotonically
! increasing. (Note that grid cell index is NOT monotonically increasing,
- ! hence we cannot set initial & final indices at the grid cell level -
+ ! hence we cannot set initial & final indices at the grid cell level -
! grc_pp%luni, grc_pp%lunf, etc.)
!
! Algorithm works as follows. The p, c, and l loops march through
! the full arrays (nump, numc, and numl) checking the "up" indexes.
- ! As soon as the "up" index of the current (p,c,l) cell changes relative
- ! to the previous (p,c,l) cell, the *i array will be set to point down
+ ! As soon as the "up" index of the current (p,c,l) cell changes relative
+ ! to the previous (p,c,l) cell, the *i array will be set to point down
! to that cell. The *f array follows the same logic, so it's always the
! last "up" index from the previous cell when an "up" index changes.
!
- ! For example, a case where p_c(1:4) = 1 and p_c(5:12) = 2. This
+ ! For example, a case where p_c(1:4) = 1 and p_c(5:12) = 2. This
! subroutine will set c_pi(1) = 1, c_pf(1) = 4, c_pi(2) = 5, c_pf(2) = 12.
!
! !USES
@@ -78,7 +79,7 @@ subroutine elm_ptrs_compdown(bounds)
!--- Loop p through full local begp:endp length
!--- Separately check the p_c, p_l, and p_g indexes for a change in
!--- the "up" index.
- !--- If there is a change, verify that the current c,l,g is within the
+ !--- If there is a change, verify that the current c,l,g is within the
!--- valid range, and set c_pi, l_pi, or g_pi to that current c,l,g
!--- Constantly update the c_pf, l_pf, and g_pf array. When the
!--- g, l, c index changes, the *_pf array will be set correctly
@@ -123,8 +124,8 @@ subroutine elm_ptrs_compdown(bounds)
lun_pp%colf(curl) = c
lun_pp%ncolumns(curl) = lun_pp%colf(curl) - lun_pp%coli(curl) + 1
enddo
-
- ! Gridcell down pointers to topounits are monotonic, so those can be done like the
+
+ ! Gridcell down pointers to topounits are monotonic, so those can be done like the
! previous monotonic down pointers
curg = 0
do t = bounds%begt,bounds%endt
@@ -142,7 +143,7 @@ subroutine elm_ptrs_compdown(bounds)
! Determine landunit_indices: indices into landunit-level arrays for each grid cell.
! Note that landunits not present in a given grid cell are set to ispval.
- ! Preliminary implementation of topounits: leave this unchanged, but will only work
+ ! Preliminary implementation of topounits: leave this unchanged, but will only work
! for max_topounits = 1
grc_pp%landunit_indices(:,bounds%begg:bounds%endg) = ispval
do l = bounds%begl,bounds%endl
@@ -184,7 +185,7 @@ subroutine elm_ptrs_compdown(bounds)
call endrun(decomp_index=l, elmlevel=namel, msg=errMsg(__FILE__, __LINE__))
end if
end do
-
+
end subroutine elm_ptrs_compdown
!------------------------------------------------------------------------------
@@ -221,7 +222,7 @@ subroutine elm_ptrs_check(bounds)
begp => bounds%begp, &
endp => bounds%endp &
)
-
+
if (masterproc) write(iulog,*) ' '
if (masterproc) write(iulog,*) '---elm_ptrs_check:'
@@ -383,7 +384,7 @@ subroutine elm_ptrs_check(bounds)
if (masterproc) write(iulog,*) ' '
end associate
-
+
end subroutine elm_ptrs_check
!-----------------------------------------------------------------------
@@ -396,7 +397,7 @@ subroutine add_topounit(ti, gi, wtgcell,elv, slp, asp,topo_ind,is_tpu_active)
!
! !ARGUMENTS:
integer , intent(inout) :: ti ! input value is index of last topounit added; output value is index of this newly-added topounit
- integer , intent(in) :: gi ! gridcell index on which this topounit should be placed
+ integer , intent(in) :: gi ! gridcell index on which this topounit should be placed
real(r8) , intent(in) :: wtgcell ! weight of the topounit relative to the gridcell
real(r8) , intent(in) :: elv ! topounit elevation
real(r8) , intent(in) :: slp ! topounit slope
@@ -415,9 +416,9 @@ subroutine add_topounit(ti, gi, wtgcell,elv, slp, asp,topo_ind,is_tpu_active)
top_pp%elevation(ti) = elv
top_pp%slope(ti) = slp
top_pp%aspect(ti) = asp
- top_pp%topo_grc_ind(ti) = topo_ind
+ top_pp%topo_grc_ind(ti) = topo_ind
top_pp%active(ti) = is_tpu_active
-
+
end subroutine add_topounit
!-----------------------------------------------------------------------
@@ -438,18 +439,18 @@ subroutine add_landunit(li, ti, ltype, wttopounit)
real(r8) , intent(in) :: wttopounit ! weight of the landunit relative to the topounit
!
! !LOCAL VARIABLES:
-
+
character(len=*), parameter :: subname = 'add_landunit'
!-----------------------------------------------------------------------
-
+
li = li + 1
lun_pp%topounit(li) = ti
lun_pp%gridcell(li) = top_pp%gridcell(ti)
-
+
lun_pp%wttopounit(li) = wttopounit
lun_pp%itype(li) = ltype
-
+
if (ltype == istsoil .or. ltype == istcrop) then
lun_pp%ifspecial(li) = .false.
else
@@ -476,6 +477,51 @@ subroutine add_landunit(li, ti, ltype, wttopounit)
end subroutine add_landunit
+!-----------------------------------------------------------------------
+ subroutine add_polygon_landunit(li, ti, ltype, wttopounit, polytype)
+ !
+ ! !DESCRIPTION:
+ ! Add an entry in the landunit-level arrays. li gives the index of the last landunit
+ ! added; the new landunit is added at li+1, and the li argument is incremented
+ ! accordingly.
+ !
+ ! This verison of add_landunit is specific to polygonal tundra.
+ !
+ ! !USES:
+ use landunit_varcon , only : istsoil, istcrop, istice_mec, istdlak, isturb_MIN, isturb_MAX
+ !
+ ! !ARGUMENTS:
+ integer , intent(inout) :: li ! input value is index of last landunit added; output value is index of this newly-added landunit
+ integer , intent(in) :: ti ! topounit index on which this landunit should be placed
+ integer , intent(in) :: ltype ! landunit type
+ real(r8) , intent(in) :: wttopounit ! weight of the landunit relative to the topounit
+ integer , intent(in) :: polytype ! defines the type of ice wedge polygon this landunit corresponds to
+ !
+ ! !LOCAL VARIABLES:
+
+ character(len=*), parameter :: subname = 'add_polygon_landunit'
+ !-----------------------------------------------------------------------
+
+ li = li + 1
+
+ lun_pp%topounit(li) = ti
+ lun_pp%gridcell(li) = top_pp%gridcell(ti)
+
+ lun_pp%wttopounit(li) = wttopounit
+ lun_pp%itype(li) = ltype
+
+ if (ltype == istsoil) then
+ lun_pp%ifspecial(li) = .false.
+ lun_pp%ispolygon(li) = .true.
+ lun_pp%polygontype(li) = polytype
+ else
+ write (iulog, *) "ERROR: attempting to assign polygonal tundra landunit to special or crop landunit type"
+ call endrun(msg=errMsg(__FILE__, __LINE__))
+ end if
+
+ end subroutine add_polygon_landunit
+
+
!-----------------------------------------------------------------------
subroutine add_column(ci, li, ctype, wtlunit)
!
@@ -499,10 +545,10 @@ subroutine add_column(ci, li, ctype, wtlunit)
col_pp%landunit(ci) = li
col_pp%topounit(ci) = lun_pp%topounit(li)
col_pp%gridcell(ci) = lun_pp%gridcell(li)
-
+
col_pp%wtlunit(ci) = wtlunit
col_pp%itype(ci) = ctype
-
+
end subroutine add_column
!-----------------------------------------------------------------------
@@ -526,17 +572,17 @@ subroutine add_patch(pi, ci, ptype, wtcol)
! !LOCAL VARIABLES:
integer :: li ! landunit index, for convenience
integer :: lb_offset ! offset between natpft_lb and 1
-
+
character(len=*), parameter :: subname = 'add_patch'
!-----------------------------------------------------------------------
-
+
pi = pi + 1
veg_pp%column(pi) = ci
veg_pp%landunit(pi) = col_pp%landunit(ci)
veg_pp%topounit(pi) = col_pp%topounit(ci)
veg_pp%gridcell(pi) = col_pp%gridcell(ci)
-
+
veg_pp%wtcol(pi) = wtcol
veg_pp%itype(pi) = ptype
diff --git a/components/elm/src/main/landunit_varcon.F90 b/components/elm/src/main/landunit_varcon.F90
index 837679cda643..46a1512a6912 100644
--- a/components/elm/src/main/landunit_varcon.F90
+++ b/components/elm/src/main/landunit_varcon.F90
@@ -12,7 +12,7 @@ module landunit_varcon
implicit none
save
private
-
+
!------------------------------------------------------------------
! Initialize landunit type constants
!------------------------------------------------------------------
@@ -35,20 +35,19 @@ module landunit_varcon
integer, parameter, public :: landunit_name_length = 40 ! max length of landunit names
character(len=landunit_name_length), public :: landunit_names(max_lunit) ! name of each landunit type
-
+
! land unit polygonal ground types
integer, parameter, public :: ilowcenpoly = 1 ! low-centered polygons
integer, parameter, public :: iflatcenpoly = 2 ! flat-centered polygons
integer, parameter, public :: ihighcenpoly = 3 ! high-centered polygons
-
- integer, parameter, public :: max_poly = 3 !maximum value that lun_pp%polygontype can have
- !(i.e., largest value in the above list)
+
+ integer, parameter, public :: max_polygon = 3 !maximum value that lun_pp%polygontype can have
+ !(i.e., largest value in the above list)
integer, parameter, public :: polygon_name_length = 40 ! max length of landunit names
- character(len=polygon_name_length), public :: polygon_names(max_poly) ! name of each landunit type
-
-
-
+ character(len=polygon_name_length), public :: polygon_names(max_polygon)! name of each landunit type
+
+
! parameters that depend on the above constants
integer, parameter, public :: numurbl = isturb_MAX - isturb_MIN + 1 ! number of urban landunits
@@ -57,7 +56,7 @@ module landunit_varcon
! !PUBLIC MEMBER FUNCTIONS:
public :: landunit_varcon_init ! initialize constants in this module
public :: landunit_is_special ! returns true if this is a special landunit
-
+
!
! !PRIVATE MEMBER FUNCTIONS:
private :: set_landunit_names ! set the landunit_names vector
@@ -65,7 +64,7 @@ module landunit_varcon
!-----------------------------------------------------------------------
contains
-
+
!-----------------------------------------------------------------------
subroutine landunit_varcon_init()
!
@@ -77,15 +76,15 @@ subroutine landunit_varcon_init()
! !ARGUMENTS:
!
! !LOCAL VARIABLES:
-
+
character(len=*), parameter :: subname = 'landunit_varcon_init'
!-----------------------------------------------------------------------
-
+
call set_landunit_names()
call set_polygon_names()
end subroutine landunit_varcon_init
-
+
!-----------------------------------------------------------------------
function landunit_is_special(ltype) result(is_special)
!
@@ -124,7 +123,7 @@ subroutine set_polygon_names
character(len=*), parameter :: not_set = 'NOT_SET'
character(len=*), parameter :: subname = 'set_polygon_names'
!-----------------------------------------------------------------------
-
+
polygon_names(:) = not_set
polygon_names(ilowcenpoly) = 'low_centered_polygons'
@@ -134,7 +133,7 @@ subroutine set_polygon_names
if (any(polygon_names == not_set)) then
call shr_sys_abort(trim(subname)//': Not all polygon names set')
end if
-
+
end subroutine set_polygon_names
!-----------------------------------------------------------------------
@@ -149,7 +148,7 @@ subroutine set_landunit_names
character(len=*), parameter :: not_set = 'NOT_SET'
character(len=*), parameter :: subname = 'set_landunit_names'
!-----------------------------------------------------------------------
-
+
landunit_names(:) = not_set
landunit_names(istsoil) = 'vegetated_or_bare_soil'
diff --git a/components/elm/src/main/surfrdMod.F90 b/components/elm/src/main/surfrdMod.F90
index e35f8320b30c..3f6c24a4d014 100755
--- a/components/elm/src/main/surfrdMod.F90
+++ b/components/elm/src/main/surfrdMod.F90
@@ -25,8 +25,8 @@ module surfrdMod
use decompMod , only : get_elmlevel_gsmap
! use spmdMod , only : iam ! rank on the land communicator
#endif
- use spmdMod
- use topounit_varcon , only : max_topounits, has_topounit
+ use spmdMod
+ use topounit_varcon , only : max_topounits, has_topounit
!
! !PUBLIC TYPES:
@@ -62,7 +62,7 @@ subroutine surfrd_get_globmask(filename, mask, ni, nj)
!
! !DESCRIPTION:
! Read the surface dataset grid related information:
- ! This is the first routine called by clm_initialize
+ ! This is the first routine called by clm_initialize
! NO DOMAIN DECOMPOSITION HAS BEEN SET YET
!
! !USES:
@@ -70,14 +70,14 @@ subroutine surfrd_get_globmask(filename, mask, ni, nj)
!
! !ARGUMENTS:
character(len=*), intent(in) :: filename ! grid filename
- integer , pointer :: mask(:) ! grid mask
+ integer , pointer :: mask(:) ! grid mask
integer , intent(out) :: ni, nj ! global grid sizes
!
! !LOCAL VARIABLES:
logical :: isgrid2d
integer :: dimid,varid ! netCDF id's
integer :: ns ! size of grid on file
- integer :: n,i,j ! index
+ integer :: n,i,j ! index
integer :: ier ! error status
type(file_desc_t) :: ncid ! netcdf id
type(var_desc_t) :: vardesc ! variable descriptor
@@ -157,7 +157,7 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename)
use elm_varctl, only : use_pflotran
!
! !ARGUMENTS:
- integer ,intent(in) :: begg, endg
+ integer ,intent(in) :: begg, endg
type(domain_type),intent(inout) :: ldomain ! domain to init
character(len=*) ,intent(in) :: filename ! grid filename
character(len=*) ,optional, intent(in) :: glcfilename ! glc mask filename
@@ -172,7 +172,7 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename)
integer :: dimid,varid ! netCDF id's
integer :: start(1), count(1) ! 1d lat/lon array sections
integer :: ier,ret ! error status
- logical :: readvar ! true => variable is on input file
+ logical :: readvar ! true => variable is on input file
logical :: isgrid2d ! true => file is 2d lat/lon
logical :: istype_domain ! true => input file is of type domain
real(r8), allocatable :: rdata2d(:,:) ! temporary
@@ -206,7 +206,7 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename)
! Determine dimensions
call ncd_inqfdims(ncid, isgrid2d, ni, nj, ns)
-
+
! pflotran:beg-----------------------------------------------
call ncd_inqdlen(ncid, dimid, nv, 'nv')
if (nv>0) then
@@ -221,10 +221,10 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename)
call domain_init(ldomain, isgrid2d=isgrid2d, ni=ni, nj=nj, nbeg=begg, nend=endg)
! Determine type of file - old style grid file or new style domain file
- call check_var(ncid=ncid, varname='LONGXY', vardesc=vardesc, readvar=readvar)
+ call check_var(ncid=ncid, varname='LONGXY', vardesc=vardesc, readvar=readvar)
if (readvar) istype_domain = .false.
- call check_var(ncid=ncid, varname='xc', vardesc=vardesc, readvar=readvar)
+ call check_var(ncid=ncid, varname='xc', vardesc=vardesc, readvar=readvar)
if (readvar) istype_domain = .true.
! Read in area, lon, lat
@@ -235,11 +235,11 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename)
! convert from radians**2 to km**2
ldomain%area = ldomain%area * (re**2)
if (.not. readvar) call endrun( msg=' ERROR: area NOT on file'//errMsg(__FILE__, __LINE__))
-
+
call ncd_io(ncid=ncid, varname= 'xc', flag='read', data=ldomain%lonc, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: xc NOT on file'//errMsg(__FILE__, __LINE__))
-
+
call ncd_io(ncid=ncid, varname= 'yc', flag='read', data=ldomain%latc, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: yc NOT on file'//errMsg(__FILE__, __LINE__))
@@ -315,11 +315,11 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename)
call ncd_io(ncid=ncid, varname= 'AREA', flag='read', data=ldomain%area, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: AREA NOT on file'//errMsg(__FILE__, __LINE__))
-
+
call ncd_io(ncid=ncid, varname= 'LONGXY', flag='read', data=ldomain%lonc, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: LONGXY NOT on file'//errMsg(__FILE__, __LINE__))
-
+
call ncd_io(ncid=ncid, varname= 'LATIXY', flag='read', data=ldomain%latc, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: LATIXY NOT on file'//errMsg(__FILE__, __LINE__))
@@ -341,7 +341,7 @@ subroutine surfrd_get_grid(begg, endg, ldomain, filename, glcfilename)
! pflotran:end-----------------------------------------------
end if
-
+
! let lat1d/lon1d data available for all grid-types, if coupled with PFLOTRAN.
if (isgrid2d .or. use_pflotran) then
allocate(rdata2d(ni,nj), lon1d(ni), lat1d(nj))
@@ -525,7 +525,7 @@ subroutine surfrd_get_topo(domain,filename)
use domainMod , only : domain_type
use fileutils , only : getfil
use GridcellType, only : grc_pp
-
+
!
! !ARGUMENTS:
type(domain_type),intent(inout) :: domain ! domain to init
@@ -564,7 +564,7 @@ subroutine surfrd_get_topo(domain,filename)
domain%ns,ns
call endrun(msg=errMsg(__FILE__, __LINE__))
endif
-
+
beg = domain%nbeg
end = domain%nend
@@ -630,7 +630,7 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat)
use topounit_varcon, only : max_topounits, has_topounit
!
! !ARGUMENTS:
- integer, intent(in) :: begg, endg
+ integer, intent(in) :: begg, endg
type(domain_type),intent(in) :: ldomain ! land domain
character(len=*), intent(in) :: lfsurdat ! surface dataset filename
!
@@ -645,7 +645,7 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat)
real(r8) :: rmaxlon,rmaxlat ! local min/max vars
type(file_desc_t) :: ncid ! netcdf id
logical :: istype_domain ! true => input file is of type domain
- logical :: isgrid2d ! true => intut grid is 2d
+ logical :: isgrid2d ! true => intut grid is 2d
character(len=32) :: subname = 'surfrd_get_data' ! subroutine name
!-----------------------------------------------------------------------
@@ -670,21 +670,21 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat)
call ncd_io(ncid=ncid, varname= 'PFTDATA_MASK', flag='read', data=ldomain%pftm, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: pftm NOT on surface dataset'//errMsg(__FILE__, __LINE__))
-
- !! Read the actual number of topounits per grid
+
+ !! Read the actual number of topounits per grid
!call check_var(ncid=ncid, varname='topoPerGrid', vardesc=vardesc, readvar=readvar)
!if (readvar) then
! call ncd_io(ncid=ncid, varname= 'topoPerGrid', flag='read', data=ldomain%num_tunits_per_grd, &
! dim1name=grlnd, readvar=readvar)
- !endif
-
+ !endif
+
! Check if fsurdat grid is "close" to fatmlndfrc grid, exit if lats/lon > 0.001
- call check_var(ncid=ncid, varname='xc', vardesc=vardesc, readvar=readvar)
+ call check_var(ncid=ncid, varname='xc', vardesc=vardesc, readvar=readvar)
if (readvar) then
istype_domain = .true.
else
- call check_var(ncid=ncid, varname='LONGXY', vardesc=vardesc, readvar=readvar)
+ call check_var(ncid=ncid, varname='LONGXY', vardesc=vardesc, readvar=readvar)
if (readvar) then
istype_domain = .false.
else
@@ -736,7 +736,7 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat)
! Obtain special landunit info
call surfrd_special(begg, endg, ncid, ldomain%ns,ldomain%num_tunits_per_grd)
-
+
! Obtain vegetated landunit info
call surfrd_veg_all(begg, endg, ncid, ldomain%ns,ldomain%num_tunits_per_grd)
@@ -768,9 +768,9 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
use topounit_varcon , only : max_topounits, has_topounit
!
! !ARGUMENTS:
- integer , intent(in) :: begg, endg
+ integer , intent(in) :: begg, endg
type(file_desc_t), intent(inout) :: ncid ! netcdf id
- integer , intent(in) :: ns ! domain size
+ integer , intent(in) :: ns ! domain size
integer , intent(in) :: ntpu(:) ! Number of topounits per grid
!
! !LOCAL VARIABLES:
@@ -781,24 +781,24 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
integer :: nindx ! temporary for error check
integer :: ier ! error status
logical :: readvar
-
+
real(r8),pointer :: pctgla_old(:) ! percent of grid cell is glacier
real(r8),pointer :: pctlak_old(:) ! percent of grid cell is lake
real(r8),pointer :: pctwet_old(:) ! percent of grid cell is wetland
real(r8),pointer :: pcturb_old(:,:) ! percent of grid cell is urbanized
- integer ,pointer :: urban_region_id_old(:)
+ integer ,pointer :: urban_region_id_old(:)
real(r8),pointer :: pctglc_mec_tot_old(:) ! percent of grid cell is glacier (sum over classes)
real(r8),pointer :: pcturb_tot_old(:) ! percent of grid cell is urban (sum over density classes)
real(r8),pointer :: pctspec_old(:) ! percent of spec lunits wrt gcell
-
+
real(r8),pointer :: pctgla(:,:) ! percent of grid cell is glacier
real(r8),pointer :: pctlak(:,:) ! percent of grid cell is lake
real(r8),pointer :: pctwet(:,:) ! percent of grid cell is wetland
real(r8),pointer :: pcturb(:,:,:) ! percent of grid cell is urbanized
- integer ,pointer :: urban_region_id(:,:)
+ integer ,pointer :: urban_region_id(:,:)
real(r8),pointer :: pctglc_mec_tot(:,:) ! percent of grid cell is glacier (sum over classes)
real(r8),pointer :: pcturb_tot(:,:) ! percent of grid cell is urban (sum over density classes)
- real(r8),pointer :: pctspec(:,:) ! percent of spec lunits wrt gcell
+ real(r8),pointer :: pctspec(:,:) ! percent of spec lunits wrt gcell
integer :: dens_index ! urban density index
character(len=32) :: subname = 'surfrd_special' ! subroutine name
real(r8) closelat,closelon
@@ -813,7 +813,7 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
allocate(urban_region_id(begg:endg,1:max_topounits))
allocate(pctglc_mec_tot(begg:endg,1:max_topounits))
allocate(pctspec(begg:endg,1:max_topounits))
-
+
call check_dim(ncid, 'nlevsoi', nlevsoifl)
! Obtain non-grid surface properties of surface dataset other than percent pft
@@ -832,7 +832,7 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
! Read urban info
if (nlevurb == 0) then
- ! If PCT_URBAN is not multi-density then set pcturb to zero
+ ! If PCT_URBAN is not multi-density then set pcturb to zero
pcturb = 0._r8
urban_valid(begg:endg,:) = .false.
write(iulog,*)'PCT_URBAN is not multi-density, pcturb set to 0'
@@ -895,7 +895,7 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
pctglc_mec_tot(:,:) = 0._r8
pctspec = pctwet + pctlak + pcturb_tot + pctgla
-
+
endif
! Error check: glacier, lake, wetland, urban sum must be less than 100
@@ -904,11 +904,11 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
do nl = begg,endg
ti = (nl - begg) + 1
if (.not. has_topounit) then
- tm = max_topounits
+ tm = max_topounits
else
tm = ntpu(ti)
end if
- do t = 1, tm
+ do t = 1, tm
if (pctspec(nl,t) > 100._r8+1.e-04_r8) then
found = .true.
nindx = nl
@@ -938,9 +938,9 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
dens_index = n - isturb_MIN + 1
wt_lunit(nl,:,n) = pcturb(nl,:,dens_index) / 100._r8
end do
-
+
end do
-
+
! Obtain firrig and surface/grnd irrigation fraction
if (firrig_data) then
call ncd_io(ncid=ncid, varname='FIRRIG', flag='read', data=firrig, &
@@ -954,13 +954,13 @@ subroutine surfrd_special(begg, endg, ncid, ns,ntpu)
call ncd_io(ncid=ncid, varname='FGRD', flag='read', data=f_grd, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( trim(subname)//' ERROR: FGRD NOT on surfdata file' )
-
+
else
firrig(:,:) = 0.7_r8
f_surf(:,:) = 1.0_r8
f_grd(:,:) = 0.0_r8
end if
-
+
call CheckUrban(begg, endg, pcturb(begg:endg,:,:), subname,ntpu)
deallocate(pctgla,pctlak,pctwet,pcturb,pcturb_tot,urban_region_id,pctglc_mec_tot,pctspec)
@@ -1008,7 +1008,7 @@ subroutine surfrd_cftformat( ncid, begg, endg, wt_cft, fert_cft, fert_p_cft, cft
call ncd_io(ncid=ncid, varname='PCT_CFT', flag='read', data=wt_cft, &
dim1name=grlnd, readvar=readvar)
- if (.not. readvar) call endrun( msg=' ERROR: PCT_CFT NOT on surfdata file'//errMsg(__FILE__, __LINE__))
+ if (.not. readvar) call endrun( msg=' ERROR: PCT_CFT NOT on surfdata file'//errMsg(__FILE__, __LINE__))
if ( cft_size > 0 )then
call ncd_io(ncid=ncid, varname='NFERT', flag='read', data=fert_cft, &
@@ -1129,16 +1129,16 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
! Determine weight arrays for non-dynamic landuse mode
!
! !USES:
- use elm_varctl , only : create_crop_landunit, use_fates
+ use elm_varctl , only : create_crop_landunit, use_fates, use_polygonal_tundra
use elm_varctl , only : irrigate
use elm_varpar , only : surfpft_lb, surfpft_ub, surfpft_size, cft_lb, cft_ub, cft_size
use elm_varpar , only : crop_prog
- use elm_varsur , only : wt_lunit, wt_nat_patch, wt_cft, fert_cft, fert_p_cft
- use landunit_varcon , only : istsoil, istcrop
+ use elm_varsur , only : wt_lunit, wt_nat_patch, wt_cft, fert_cft, fert_p_cft, wt_polygon
+ use landunit_varcon , only : istsoil, istcrop, ilowcenpoly, iflatcenpoly, ihighcenpoly
use pftvarcon , only : nc3crop, nc3irrig, npcropmin
use pftvarcon , only : ncorn, ncornirrig, nsoybean, nsoybeanirrig
use pftvarcon , only : nscereal, nscerealirrig, nwcereal, nwcerealirrig
- use pftvarcon , only : ncassava, ncotton, nfoddergrass, noilpalm, nograins, nrapeseed
+ use pftvarcon , only : ncassava, ncotton, nfoddergrass, noilpalm, nograins, nrapeseed
use pftvarcon , only : nrice, nrtubers, nsugarcane, nmiscanthus, nswitchgrass, nwillow, npoplar
use pftvarcon , only : ncassavairrig, ncottonirrig, nfoddergrassirrig, noilpalmirrig, nograinsirrig, nrapeseedirrig
use pftvarcon , only : nriceirrig, nrtubersirrig, nsugarcaneirrig, nmiscanthusirrig, nswitchgrassirrig, nwillowirrig, npoplarirrig
@@ -1153,7 +1153,7 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
! !LOCAL VARIABLES:
integer :: nl, t ! index
integer :: dimid,varid ! netCDF id's
- integer :: ier ! error status
+ integer :: ier ! error status
integer :: cftsize ! size of CFT's
logical :: readvar ! is variable on dataset
logical :: cft_dim_exists ! does the dimension 'cft' exist on the dataset?
@@ -1173,15 +1173,35 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
call ncd_io(ncid=ncid, varname='PCT_NATVEG', flag='read', data=arrayl, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: PCT_NATVEG NOT on surfdata file'//errMsg(__FILE__, __LINE__))
- wt_lunit(begg:endg,1:max_topounits,istsoil) = arrayl(begg:endg,1:max_topounits)
+ wt_lunit(begg:endg,1:max_topounits,istsoil) = arrayl(begg:endg,1:max_topounits)
+
+ if (use_polygonal_tundra) then
+ call ncd_io(ncid=ncid, varname='PCT_HCP', flag='read', data=arrayl, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( msg=' ERROR: use_polygonal_tundra = .true., but PCT_HCP NOT on surfdata file'//errMsg(__FILE__, __LINE__))
+ wt_polygon(begg:endg,1:max_topounits,ihighcenpoly) = arrayl(begg:endg,1:max_topounits)
+
+ call ncd_io(ncid=ncid, varname='PCT_FCP', flag='read', data=arrayl, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( msg=' ERROR: use_polygonal_tundra = .true., but PCT_FCP NOT on surfdata file'//errMsg(__FILE__, __LINE__))
+ wt_polygon(begg:endg,1:max_topounits,iflatcenpoly) = arrayl(begg:endg,1:max_topounits)
+
+ call ncd_io(ncid=ncid, varname='PCT_LCP', flag='read', data=arrayl, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( msg=' ERROR: use_polygonal_tundra = .true., but PCT_LCP NOT on surfdata file'//errMsg(__FILE__, __LINE__))
+ wt_polygon(begg:endg,1:max_topounits,ilowcenpoly) = arrayl(begg:endg,1:max_topounits)
+ endif
+
+ ! add two other types
call ncd_io(ncid=ncid, varname='PCT_CROP', flag='read', data=arrayl, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: PCT_CROP NOT on surfdata file'//errMsg(__FILE__, __LINE__))
- wt_lunit(begg:endg,1:max_topounits,istcrop) = arrayl(begg:endg,1:max_topounits)
+ wt_lunit(begg:endg,1:max_topounits,istcrop) = arrayl(begg:endg,1:max_topounits)
+
deallocate(arrayl)
-
+
! Check the file format for CFT's and handle accordingly
call ncd_inqdid(ncid, 'cft', dimid, cft_dim_exists)
if ( cft_dim_exists .and. create_crop_landunit ) then
@@ -1214,7 +1234,7 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
else
call endrun( msg=' ERROR: New format surface datasets require create_crop_landunit TRUE'//errMsg(__FILE__, __LINE__))
end if
-
+
else
! PFTs contain the crops but create_crop_landunit = .true.
@@ -1245,12 +1265,15 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
end if
wt_lunit(begg:endg,:,istsoil) = wt_lunit(begg:endg,:,istsoil) / 100._r8
wt_lunit(begg:endg,:,istcrop) = wt_lunit(begg:endg,:,istcrop) / 100._r8
+ wt_polygon(begg:endg,:,ihighcenpoly) = wt_polygon(begg:endg,:,ihighcenpoly) / 100._r8
+ wt_polygon(begg:endg,:,iflatcenpoly) = wt_polygon(begg:endg,:,iflatcenpoly) / 100._r8
+ wt_polygon(begg:endg,:,ilowcenpoly) = wt_polygon(begg:endg,:,ilowcenpoly) / 100._r8
wt_nat_patch(begg:endg,:,:) = wt_nat_patch(begg:endg,:,:) / 100._r8
!call check_sums_equal_1_3d(wt_nat_patch, begg, 'wt_nat_patch', subname,ntpu)
call check_sums_equal_1_3d(wt_nat_patch, begg, 'wt_nat_patch', subname)
! If no irrigation, merge irrigated CFTs with rainfed
-
+
if (crop_prog .and. .not. irrigate) then
if (masterproc) then
write(iulog,*) trim(subname)//' crop=.T. and irrigate=.F., so merging irrigated pfts with rainfed'
@@ -1495,7 +1518,7 @@ subroutine surfrd_get_grid_conn(filename, cellsOnCell, edgesOnCell, &
call ncd_pio_closefile(ncid)
end subroutine surfrd_get_grid_conn
-
+
!-----------------------------------------------------------------------------------------------------
subroutine surfrd_topounit_data(begg, endg, lfsurdat)
!
@@ -1506,14 +1529,14 @@ subroutine surfrd_topounit_data(begg, endg, lfsurdat)
use ncdio_pio , only : file_desc_t, var_desc_t, ncd_pio_openfile, ncd_pio_closefile
use ncdio_pio , only : ncd_io, check_var, ncd_inqfdims, check_dim, ncd_inqdid, ncd_inqdlen
use elm_varctl , only: fsurdat
- use fileutils , only : getfil
+ use fileutils , only : getfil
use GridcellType , only : grc_pp
use elm_varsur , only : wt_tunit, elv_tunit, slp_tunit, asp_tunit,num_tunit_per_grd
use topounit_varcon , only : max_topounits, has_topounit
-
+
!
! !ARGUMENTS:
- integer , intent(in) :: begg, endg
+ integer , intent(in) :: begg, endg
character(len=*), intent(in) :: lfsurdat ! surface dataset filename
!
@@ -1521,23 +1544,23 @@ subroutine surfrd_topounit_data(begg, endg, lfsurdat)
type(var_desc_t) :: vardesc
integer :: n,t ! indices
character(len=256):: locfn ! local file name
-
+
logical :: readvar
integer :: dimid
type(file_desc_t) :: ncid ! netcdf id
-
+
real(r8),pointer :: maxTopoElv(:) ! Maximum topounit elevation
integer ,pointer :: numTopoPerGrid(:) ! Number of topounits per grid
real(r8),pointer :: TopounitFracArea(:,:) ! Topounit fractional area
real(r8) ,pointer :: TopounitElv(:,:) ! Topounit elevation
- real(r8),pointer :: TopounitSlope(:,:) ! Topounit slope
+ real(r8),pointer :: TopounitSlope(:,:) ! Topounit slope
integer ,pointer :: TopounitAspect(:,:) ! Topounit aspect
integer ,pointer :: num_topo_per_grid(:) ! Topounit aspect
real(r8),pointer :: GridElevation(:) ! Topounit aspect
! integer ,pointer :: TopounitIndices(:,:) ! Topounit indices in each grid
-
+
character(len=32) :: subname = 'surfrd_topounit_data' ! subroutine name
-!-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
allocate(maxTopoElv(begg:endg))
allocate(GridElevation(begg:endg))
@@ -1548,11 +1571,11 @@ subroutine surfrd_topounit_data(begg, endg, lfsurdat)
allocate(TopounitAspect(begg:endg,max_topounits))
allocate(num_topo_per_grid(begg:endg))
! allocate(TopounitIndices(begg:endg,max_topounits))
-
+
! Read surface data
call getfil( lfsurdat, locfn, 0 )
call ncd_pio_openfile (ncid, trim(locfn), 0)
-
+
!call check_dim(ncid, 'nlevsoi', nlevsoifl)
call check_var(ncid=ncid, varname='MaxTopounitElv', vardesc=vardesc, readvar=readvar)
if (readvar) then
@@ -1589,36 +1612,36 @@ subroutine surfrd_topounit_data(begg, endg, lfsurdat)
call ncd_io(ncid=ncid, varname='TopounitAspect', flag='read', data=TopounitAspect, &
dim1name=grlnd, readvar=readvar)
endif
-
+
call check_var(ncid=ncid, varname='topoPerGrid', vardesc=vardesc, readvar=readvar)
if (readvar) then
call ncd_io(ncid=ncid, varname='topoPerGrid', flag='read', data=num_topo_per_grid, &
dim1name=grlnd, readvar=readvar)
endif
-
+
call check_var(ncid=ncid, varname='TOPO2', vardesc=vardesc, readvar=readvar)
if (readvar) then
call ncd_io(ncid=ncid, varname='TOPO2', flag='read', data=GridElevation, &
dim1name=grlnd, readvar=readvar)
endif
if (readvar) then
- do n = begg,endg
- grc_pp%MaxElevation(n) = maxTopoElv(n)
- grc_pp%elevation(n) = GridElevation(n)
+ do n = begg,endg
+ grc_pp%MaxElevation(n) = maxTopoElv(n)
+ grc_pp%elevation(n) = GridElevation(n)
num_tunit_per_grd(n) = num_topo_per_grid(n)
grc_pp%ntopounits(n) = numTopoPerGrid(n)
do t = 1, max_topounits
wt_tunit(n,t) = TopounitFracArea(n,t)
elv_tunit(n,t) = TopounitElv(n,t)
! slp_tunit(n,t) = TopounitSlope(n,t)
- ! asp_tunit(n,t) = TopounitAspect(n,t)
+ ! asp_tunit(n,t) = TopounitAspect(n,t)
end do
- end do
- endif
+ end do
+ endif
deallocate(maxTopoElv,TopounitFracArea,TopounitElv,TopounitSlope,TopounitAspect,GridElevation)
-
+
call ncd_pio_closefile(ncid)
-
+
end subroutine surfrd_topounit_data
@@ -1677,7 +1700,7 @@ subroutine surfrd_get_topo_for_solar_rad(domain,filename)
domain%ns,ns
call endrun()
endif
-
+
beg = domain%nbeg
end = domain%nend
@@ -1734,18 +1757,18 @@ subroutine surfrd_fates_nocropmod( ncid, begg, endg )
! and will use a mapping table to connect its own pft and cft
! definitions to those it finds in the surface file.
!--------------------------------------------------------------------------
-
+
! !USES:
use elm_varsur , only : wt_nat_patch, wt_lunit
use elm_varpar , only : cft_size, surfpft_lb, surfpft_ub
use landunit_varcon , only : istsoil, istcrop
use topounit_varcon , only : max_topounits
-
+
! !ARGUMENTS:
implicit none
type(file_desc_t), intent(inout) :: ncid ! netcdf id
integer , intent(in) :: begg, endg
-
+
!
! !LOCAL VARIABLES:
logical :: readvar ! is variable on dataset
@@ -1753,9 +1776,9 @@ subroutine surfrd_fates_nocropmod( ncid, begg, endg )
real(r8),pointer :: array3d_cft(:,:,:) ! local array
integer :: g,p,t
integer :: cft_dimlen,surfpft_dimlen,dimid
-
+
character(len=32) :: subname = 'surfrd_fates_nocropmod'! subroutine name
-
+
call ncd_inqdlen(ncid, dimid, cft_dimlen, 'cft')
call ncd_inqdlen(ncid, dimid, surfpft_dimlen, 'natpft')
@@ -1764,14 +1787,14 @@ subroutine surfrd_fates_nocropmod( ncid, begg, endg )
call endrun( msg=' ERROR: PCT+CFT dimlen does not match array size for wt_nat_patch when fates is on'&
//errMsg(__FILE__, __LINE__))
end if
-
+
allocate( array3d_cft(begg:endg,1:max_topounits,1:cft_dimlen) )
allocate( array3d_pft(begg:endg,1:max_topounits,1:surfpft_dimlen) )
-
+
call ncd_io(ncid=ncid, varname='PCT_CFT', flag='read', data=array3d_cft, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: PCT_CFT NOT on surfdata file'//errMsg(__FILE__, __LINE__))
-
+
call ncd_io(ncid=ncid, varname='PCT_NAT_PFT', flag='read', data=array3d_pft, &
dim1name=grlnd, readvar=readvar)
if (.not. readvar) call endrun( msg=' ERROR: PCT_NAT_PFT NOT on surfdata file'//errMsg(__FILE__, __LINE__))
@@ -1781,7 +1804,7 @@ subroutine surfrd_fates_nocropmod( ncid, begg, endg )
wt_nat_patch(begg:,:,0:surfpft_dimlen-1) = array3d_pft(begg:,:,:)
wt_nat_patch(begg:,:,surfpft_dimlen:surfpft_dimlen+cft_dimlen-1) = array3d_cft(begg:,:,:)
-
+
do g = begg, endg
do t = 1, max_topounits
if ( wt_lunit(g,t,istcrop) > 0.0_r8 )then
@@ -1797,10 +1820,10 @@ subroutine surfrd_fates_nocropmod( ncid, begg, endg )
end if
end do
end do
-
+
deallocate(array3d_cft,array3d_pft)
-
-
+
+
end subroutine surfrd_fates_nocropmod
-
+
end module surfrdMod