Skip to content
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

mksurfdata_esmf: Refactor mksoiltex #2737

Open
wants to merge 3 commits into
base: b4b-dev
Choose a base branch
from
Open

mksurfdata_esmf: Refactor mksoiltex #2737

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
761e6c7
Add subroutine for filling output arrays.
samsrabin Aug 30, 2024
7d305c9
Put back commented-out organic_o option 2.
samsrabin Aug 30, 2024
9a2759f
Fixes to organic_o option 2.
samsrabin Sep 4, 2024
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
293 changes: 115 additions & 178 deletions tools/mksurfdata_esmf/src/mksoiltexMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,90 @@ module mksoiltexMod
contains
!=================================================================================

subroutine mksoiltex_i_to_o(no, lookup_index, n_scid, nlay, n_mapunits, name, val_neg_4, val_neg_other, data_i, data_o, data_modifier, organic_o)
!
! Arguments
integer, intent(in) :: no
integer, intent(in) :: lookup_index
integer, intent(in) :: n_scid, nlay, n_mapunits
character(*), intent(in) :: name
real(r4), intent(in) :: val_neg_4 ! Fallback value if read-in value is -4
real(r4), intent(in) :: val_neg_other ! Fallback value if read-in value is negative but not -4
real(r4), intent(in) :: data_i(:,:,:)
real(r4), intent(out) :: data_o(:,:)
real(r4), optional, intent(in) :: data_modifier(:,:,:)
real(r4), optional, intent(out) :: organic_o(:,:)
!
! Local variables
integer :: l
integer :: ier
logical :: is_neg_4
logical :: organic2
real(r4), allocatable :: organic_i(:,:,:)

! Apply data_modifer, if needed (organic_o OPTION 2)
organic2 = .false.
if (present(data_modifier)) then
if (.not. present(organic_o)) then
call shr_sys_abort('mksoiltex_i_to_o: data_modifier not needed if not providing organic_o')
else if (trim(name) /= 'orgc') then
call shr_sys_abort('mksoiltex_i_to_o: organic_o should only be provided along with orgc')
end if
organic2 = .true.
allocate(organic_i(nlay,n_scid,n_mapunits), stat=ier)
if (ier/=0) call shr_sys_abort()
organic_i = data_i * data_modifier
else if (present(organic_o)) then
call shr_sys_abort('mksoiltex_i_to_o: If providing organic_o, also provide data_modifier')
end if

! Fill first layer of output array with first positive value on SCID dim of input array
data_o(no,1) = data_i(1,1,lookup_index)
if (organic2) organic_o(no,1) = organic_i(1,1,lookup_index)
if (data_o(no,1) < 0._r4) then
do l = 2,n_scid
if (data_i(1,l,lookup_index) >= 0._r4) then
data_o(no,1) = data_i(1,l,lookup_index)
if (organic2) organic_o(no,1) = organic_i(1,l,lookup_index) ! TODO: This should probably be under a conditional looking at organic_i instead of data_i, but that's not how it was in the original
exit
end if
end do
end if

! Handle cases where no positive value was found
if (data_o(no,1) < 0._r4) then
is_neg_4 = int(data_o(no,1)) == -4 ! TODO: Rename. Maybe indicates sand dune?
if (is_neg_4) then
data_o(no,:) = val_neg_4
if (organic2) organic_o(no,:) = val_neg_4 ! TODO: This should probably be under conditionals looking at organic_o instead of data_o, but that's not how it was in the original
else
data_o(no,:) = val_neg_other
if (organic2) organic_o(no,:) = val_neg_other ! TODO: This should probably be under conditionals looking at organic_o instead of data_o, but that's not how it was in the original
end if
end if

! Ensure negative values are handled
! TODO: Remove as unnecessary?
! TODO: Also handle organic_o? Not handled in original
if (data_o(no,1) < 0._r4 .and. trim(name) /= "sand") then
write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index
call shr_sys_abort('could not find a value >= 0 for '//name)
end if

! Top soil layer is filled above. Here, we fill the other layers.
do l = 2,nlay
data_o(no,l) = data_i(l,1,lookup_index)
if (organic2) organic_o(no,l) = organic_i(l,1,lookup_index)

! If a layer is negative, fill it with the previous layer's value
if (data_o(no,l) < 0._r4) then
data_o(no,l) = data_o(no,l-1)
if (organic2) organic_o(no,l) = organic_o(no,l-1) ! TODO: This should probably be under a conditional looking at organic_o instead of data_o, but that's not how it was in the original
end if
end do

end subroutine mksoiltex_i_to_o

subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o, rc)
!
! make %sand, %clay, organic carbon content, coarse fragments, bulk density,
Expand Down Expand Up @@ -91,6 +175,7 @@ subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o
type(var_desc_t) :: pio_varid_bulk
type(var_desc_t) :: pio_varid_phaq
type(var_desc_t) :: pio_varid_organic
integer, parameter :: organic_o_option = 1
integer :: starts(3) ! starting indices for reading lookup table
integer :: counts(3) ! dimension counts for reading lookup table
integer :: srcTermProcessing_Value = 0
Expand All @@ -108,6 +193,10 @@ subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o
write(ndiag,'(a)') ' Input mesh/grid file is '//trim(file_mesh_i)
end if

if (organic_o_option < 1 .or. organic_o_option > 2) then
call shr_sys_abort('organic_o_option must be 1 or 2')
end if

! Determine ns_o and allocate output data
call ESMF_MeshGet(mesh_o, numOwnedElements=ns_o, rc=rc)
if (ChkErr(rc,__LINE__,u_FILE_u)) return
Expand Down Expand Up @@ -318,180 +407,26 @@ subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o
! Determine lookup_index
lookup_index = mapunit_lookup(mapunit_o(no))

! Determine top soil layer sand_o
! If less than 0 search within the SCID array for the first index
! that gives a value greater than or equal to 0
! Then determine the other soil layers
sand_o(no,1) = float(sand_i(1,1,lookup_index))
if (sand_o(no,1) < 0._r4) then
do l = 2,n_scid
if (float(sand_i(1,l,lookup_index)) >= 0._r4) then
sand_o(no,1) = float(sand_i(1,l,lookup_index))
exit
end if
end do
end if
if (sand_o(no,1) < 0._r4) then
if (int(sand_o(no,1)) == -4) then
sand_o(no,:) = 99._r4
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@samsrabin I haven't looked at this in a long while, but judging from sand_o + clay_o = 100 when these two start equal to -4, I will agree with your assessment that the flag -4 indicates sand dunes regardless of which variable we look at.

else
sand_o(no,:) = 43._r4
end if
end if
do l = 2,nlay
sand_o(no,l) = float(sand_i(l,1,lookup_index))
if (sand_o(no,l) < 0._r4) then
sand_o(no,l) = sand_o(no,l-1)
end if
end do

! Same algorithm for clay_o as for sand_o
clay_o(no,1) = float(clay_i(1,1,lookup_index))
if (clay_o(no,1) < 0._r4) then
do l = 2,n_scid
if (float(clay_i(1,l,lookup_index)) >= 0._r4) then
clay_o(no,1) = float(clay_i(1,l,lookup_index))
exit
end if
end do
end if
if (clay_o(no,1) < 0._r4) then
if (int(clay_o(no,1)) == -4) then
clay_o(no,:) = 1._r4
else
clay_o(no,:) = 18._r4
end if
end if
if (clay_o(no,1) < 0._r4) then
write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index
call shr_sys_abort('could not find a value >= 0 for clay_i')
end if
do l = 2,nlay
clay_o(no,l) = float(clay_i(l,1,lookup_index))
if (clay_o(no,l) < 0._r4) then
clay_o(no,l) = clay_o(no,l-1)
end if
end do

! Same algorithm for orgc_o as for sand_o
! organic_o OPTION 2 (commented out)
! Calculate from multiple input variables to get the output variable
orgc_o(no,1) = orgc_i(1,1,lookup_index)
! organic_o(no,1) = orgc_i(1,1,lookup_index) * bulk_i(1,1,lookup_index) * float(100 - cfrag_i(1,1,lookup_index)) * 0.01_r4 / 0.58_r4
if (orgc_o(no,1) < 0._r4) then
do l = 2,n_scid
if (orgc_i(1,l,lookup_index) >= 0._r4) then
orgc_o(no,1) = orgc_i(1,l,lookup_index)
! organic_o(no,1) = orgc_i(1,l,lookup_index) * bulk_i(1,l,lookup_index) * float(100 - cfrag_i(1,l,lookup_index)) * 0.01_r4 / 0.58_r4
exit
end if
end do
end if
if (orgc_o(no,1) < 0._r4) then
if (int(orgc_o(no,1)) == -4) then ! sand dunes
orgc_o(no,:) = 1._r4
! organic_o(no,:) = 1._r4
else
orgc_o(no,:) = 0._r4
! organic_o(no,:) = 0._r4
end if
end if
if (orgc_o(no,1) < 0._r4) then
write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index
call shr_sys_abort('could not find a value >= 0 for orgc_i')
end if
do l = 2,nlay
orgc_o(no,l) = orgc_i(l,1,lookup_index)
! organic_o(no,l) = orgc_i(l,1,lookup_index) * bulk_i(l,1,lookup_index) * float(100 - cfrag_i(l,1,lookup_index)) * 0.01_r4 / 0.58_r4
if (orgc_o(no,l) < 0._r4) then
orgc_o(no,l) = orgc_o(no,l-1)
! organic_o(no,l) = organic_o(no,l-1)
end if
end do

! Same algorithm for cfrag_o as for sand_o
cfrag_o(no,1) = float(cfrag_i(1,1,lookup_index))
if (cfrag_o(no,1) < 0._r4) then
do l = 2,n_scid
if (float(cfrag_i(1,l,lookup_index)) >= 0._r4) then
cfrag_o(no,1) = float(cfrag_i(1,l,lookup_index))
exit
end if
end do
end if
if (cfrag_o(no,1) < 0._r4) then
if (int(cfrag_o(no,1)) == -4) then ! sand dunes
cfrag_o(no,:) = 1._r4
else
cfrag_o(no,:) = 0._r4
end if
end if
if (cfrag_o(no,1) < 0._r4) then
write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index
call shr_sys_abort('could not find a value >= 0 for cfrag_i')
! Fill output arrays
call mksoiltex_i_to_o(no, lookup_index, n_scid, nlay, n_mapunits, &
"sand", 99._r4, 43._r4, float(sand_i(:,:,:)), sand_o)
call mksoiltex_i_to_o(no, lookup_index, n_scid, nlay, n_mapunits, &
"clay", 1._r4, 18._r4, float(clay_i(:,:,:)), clay_o)
call mksoiltex_i_to_o(no, lookup_index, n_scid, nlay, n_mapunits, &
"cfrag", 1._r4, 0._r4, float(cfrag_i(:,:,:)), cfrag_o)
call mksoiltex_i_to_o(no, lookup_index, n_scid, nlay, n_mapunits, &
"bulk", 1.5_r4, 1.5_r4, bulk_i, bulk_o) ! TODO: 1.5 ok for sand dunes and -7?
call mksoiltex_i_to_o(no, lookup_index, n_scid, nlay, n_mapunits, &
"phaq", 7._r4, 7._r4, phaq_i, phaq_o)
if (organic_o_option == 1) then
call mksoiltex_i_to_o(no, lookup_index, n_scid, nlay, n_mapunits, &
"orgc", 1._r4, 0._r4, orgc_i, orgc_o)
else
call mksoiltex_i_to_o(no, lookup_index, n_scid, nlay, n_mapunits, &
"orgc", 1._r4, 0._r4, orgc_i, orgc_o, &
orgc_i * bulk_i * float(100 - cfrag_i) * 0.01_r4 / 0.58_r4, &
organic_o)
end if
do l = 2,nlay
cfrag_o(no,l) = float(cfrag_i(l,1,lookup_index))
if (cfrag_o(no,l) < 0._r4) then
cfrag_o(no,l) = cfrag_o(no,l-1)
end if
end do

! Same algorithm for bulk_o as for sand_o
bulk_o(no,1) = bulk_i(1,1,lookup_index)
if (bulk_o(no,1) < 0._r4) then
do l = 2,n_scid
if (bulk_i(1,l,lookup_index) >= 0._r4) then
bulk_o(no,1) = bulk_i(1,l,lookup_index)
exit
end if
end do
end if
if (bulk_o(no,1) < 0._r4) then
if (int(bulk_o(no,1)) == -4) then ! sand dunes
bulk_o(no,:) = 1.5_r4 ! TODO Ok for sand dunes?
else
bulk_o(no,:) = 1.5_r4 ! TODO Ok for -7?
end if
end if
if (bulk_o(no,1) < 0._r4) then
write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index
call shr_sys_abort('could not find a value >= 0 for bulk_i')
end if
do l = 2,nlay
bulk_o(no,l) = bulk_i(l,1,lookup_index)
if (bulk_o(no,l) < 0._r4) then
bulk_o(no,l) = bulk_o(no,l-1)
end if
end do

! Same algorithm for phaq_o as for sand_o
phaq_o(no,1) = phaq_i(1,1,lookup_index)
if (phaq_o(no,1) < 0._r4) then
do l = 2,n_scid
if (phaq_i(1,l,lookup_index) >= 0._r4) then
phaq_o(no,1) = phaq_i(1,l,lookup_index)
exit
end if
end do
end if
if (phaq_o(no,1) < 0._r4) then
if (int(phaq_o(no,1)) == -4) then ! sand dunes
phaq_o(no,:) = 7._r4
else
phaq_o(no,:) = 7._r4
end if
end if
if (phaq_o(no,1) < 0._r4) then
write(6,*)'ERROR: at no, lookup_index = ',no,lookup_index
call shr_sys_abort('could not find a value >= 0 for phaq_i')
end if
do l = 2,nlay
phaq_o(no,l) = phaq_i(l,1,lookup_index)
if (phaq_o(no,l) < 0._r4) then
phaq_o(no,l) = phaq_o(no,l-1)
end if
end do

! ---------------------------------------------------------------
! organic_o OPTION 1, as we plan to calculate organic in the CTSM
Expand All @@ -504,11 +439,13 @@ subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o
! (calculated from orgc_i, cfrag_i, and bulk_i) to organic_o. This
! approach first calculates organic_i and then regrids to organic_o
! rather than regridding all the terms first and then calculating
! organic_o. Commented out code above belongs to that option.
do l = 1, nlay
organic_o(no,l) = orgc_o(no,l) * bulk_o(no,l) * &
(100._r4 - cfrag_o(no,l)) * 0.01_r4 / 0.58_r4
end do
! organic_o. That would be organic_o_option 2.
if (organic_o_option == 1) then
do l = 1, nlay
organic_o(no,l) = orgc_o(no,l) * bulk_o(no,l) * &
(100._r4 - cfrag_o(no,l)) * 0.01_r4 / 0.58_r4
end do
end if

end if

Expand Down
Loading