Skip to content

Add options to core_init_atmosphere for P-X LSM #1010

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jherwehe wants to merge 6 commits into
base: develop
Choose a base branch
from
Open

Add options to core_init_atmosphere for P-X LSM #1010

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
799126f
Add options to core_init_atmosphere for P-X LSM
jherwehe Dec 9, 2022
7c3b0a4
Updates to resolve code conflicts
jherwehe Jan 10, 2023
64d3f66
* In the "MPAS-A Add options to core_init_atmosphere for P-X LSM":
ldfowler58 Nov 1, 2024
a46f6fb
* In ./src/core_init_atmosphere, removed sourcecode related to the im…
ldfowler58 Feb 9, 2025
13b7a96
* In ./src/core_init_atmosphere, removed sourcecode related to the in…
ldfowler58 Feb 9, 2025
2f6dde9
* In ./src/core_init_atmosphere/mpas_init_atm_surface.F, removed the …
ldfowler58 Mar 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
115 changes: 86 additions & 29 deletions src/core_atmosphere/physics/mpas_atmphys_initialize_real.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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).
Expand Down Expand Up @@ -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

Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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/), &
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading