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

Modify TFocean vertical interpolation for depths outside of data range #134

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

Modify TFocean vertical interpolation for depths outside of data range #134

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3570441
Modify TFocean vertical interpolation for depths outside of data range
matthewhoffman Nov 26, 2024
e2055b7
Fix indexing errors
matthewhoffman Nov 28, 2024
a48ed3b
Fix new indexing error introduced trying to fix indexing errors
matthewhoffman Nov 30, 2024
cbb5ee3
Improve error checking in ocean extrap
matthewhoffman Dec 1, 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
63 changes: 46 additions & 17 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1162,6 +1162,8 @@ end subroutine calc_iceshelf_draft_info

subroutine iceshelf_melt_ismip6(domain, err)

use li_constants, only: oceanFreezingTempDepthDependence

!-----------------------------------------------------------------
! input/output variables
!-----------------------------------------------------------------
Expand Down Expand Up @@ -1249,38 +1251,65 @@ subroutine iceshelf_melt_ismip6(domain, err)
mean_TF(:) = 0.d0
IS_area(:) = 0.d0

! Check zOcean for valid values
if (minval(zOcean) == 0.0_RKIND) then
call mpas_log_write("Invalid value for zOcean. It should have negative values but min value of 0.0 was found", &
MPAS_LOG_ERR)
err = ior(err, 1)
endif
if (maxval(zOcean) > 0.0_RKIND) then
call mpas_log_write("Invalid value for zOcean. It should have negative values but max value greater than 0.0 was found", &
MPAS_LOG_ERR)
err = ior(err, 1)
endif

do iCell = 1, nCellsSolve

if ( li_mask_is_floating_ice(cellMask(iCell)) ) then
! 1 - Linear interpolation of the thermal forcing on the ice draft depth :
ksup=1
do kk=2,nOceanLayers-1
ksup=0
do kk=1,nOceanLayers
if ( zOcean(kk) >= lowerSurface(iCell) ) ksup = kk
enddo
kinf = ksup + 1
if ((zOcean(ksup)-zOcean(kinf)) == 0) then
call mpas_log_write("iceshelf_melt_ismip6: Invalid value for zOcean. " // &
"ksup=$i kinf=$i zOcean(ksup)=$r zOcean(kinf)=$r indexToCellID=$i lowerSurface=$r", MPAS_LOG_ERR, &
intArgs=(/ksup, kinf, indexToCellID(iCell)/), &
realArgs=(/zOcean(ksup), zOcean(kinf), lowerSurface(iCell) /) )
err = ior(err, 1)
endif

!call mpas_log_write("kinf=$i, zOcean(kinf)=$r, TFocean=$r",realArgs=(/zOcean(kinf),TFocean(kinf,iCell)/), &
! intArgs=(/kinf/))
!call mpas_log_write("ksup=$i, zOcean(ksup)=$r, TFocean=$r",realArgs=(/zOcean(ksup),TFocean(ksup,iCell)/), &
! intArgs=(/ksup/))

! check if any invalid TFocean value is used for calculating TF at the draft
if ( (TFocean(kinf, iCell) == invalid_value_TF .or. TFocean(ksup, iCell) == invalid_value_TF) ) then
call mpas_log_write("iceshelf_melt_ismip6: Invalid value for TFocean. " // &
"ksup=$i kinf=$i TFocean(ksup, iCell)=$r TFocean(kinf,iCell)=$r indexToCellID=$i", MPAS_LOG_ERR, &
intArgs=(/ksup, kinf, indexToCellID(iCell)/), &
realArgs=(/TFocean(ksup, iCell), TFocean(kinf, iCell) /) )
err = ior(err, 1)
if (ksup >= 1) then
if (TFocean(ksup, iCell) == invalid_value_TF) then
call mpas_log_write("iceshelf_melt_ismip6: Invalid value for TFocean. " // &
"ksup=$i TFocean(ksup, iCell)=$r indexToCellID=$i", MPAS_LOG_ERR, &
intArgs=(/ksup, indexToCellID(iCell)/), &
realArgs=(/TFocean(ksup, iCell) /) )
err = ior(err, 1)
endif
endif
if (kinf <= nOceanLayers) then
if (TFocean(kinf, iCell) == invalid_value_TF) then
call mpas_log_write("iceshelf_melt_ismip6: Invalid value for TFocean. " // &
"kinf=$i TFocean(kinf,iCell)=$r indexToCellID=$i", MPAS_LOG_ERR, &
intArgs=(/kinf, indexToCellID(iCell)/), &
realArgs=(/TFocean(kinf, iCell) /) )
err = ior(err, 1)
endif
endif

TFdraft(iCell) = ( (zOcean(ksup)-lowerSurface(iCell)) * TFocean(kinf, iCell) &
+ (lowerSurface(iCell)-zOcean(kinf)) * TFocean(ksup, iCell) ) / (zOcean(ksup)-zOcean(kinf))
if (ksup == 0) then
! For depths shallower than shallowest layer center, use shallowest layer
TFdraft(iCell) = TFocean(1, iCell)
elseif (ksup == nOceanLayers) then
! for depths below the deepest layer center, use deepest layer corrected for Tfreeze
TFdraft(iCell) = TFocean(nOceanLayers, iCell) - &
(zOcean(nOceanLayers) - lowerSurface(iCell)) * oceanFreezingTempDepthDependence
else
! for depths between the first and last layer centers, linearly interpolate
TFdraft(iCell) = ( (zOcean(ksup)-lowerSurface(iCell)) * TFocean(kinf, iCell) &
+ (lowerSurface(iCell)-zOcean(kinf)) * TFocean(ksup, iCell) ) / (zOcean(ksup)-zOcean(kinf))
endif

! 2 - Mean Thermal forcing in individual basins (NB: fortran norm while basins start at zero):
mean_TF(basinNumber(iCell)+1) = mean_TF(basinNumber(iCell)+1) + areaCell(iCell) * TFdraft(iCell)
Expand Down
16 changes: 11 additions & 5 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_ocean_extrap.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -308,7 +308,7 @@ subroutine li_ocean_extrap_solve(domain, err)
call mpas_timer_stop("vertical scheme")

if (err > 0) then
call mpas_log_write("Ocean extraolation main iteration loop has encountered an error", MPAS_LOG_ERR)
call mpas_log_write("Ocean extrapolation main iteration loop has encountered an error", MPAS_LOG_ERR)
return
endif
enddo
Expand Down Expand Up @@ -378,6 +378,7 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali
real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography, areaCell
integer, pointer :: nCells, nCellsSolve, nISMIP6OceanLayers, nCellsExtra
integer, dimension(:), pointer :: cellMask, nEdgesOnCell
integer, dimension(:), pointer :: indexToCellID
integer, dimension(:,:), pointer :: cellsOnCell
integer :: iCell, jCell, iLayer, iNeighbor, iter
integer :: localLoopCount
Expand All @@ -399,6 +400,7 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
Expand Down Expand Up @@ -430,8 +432,9 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali
if ( validOceanMaskOld(iLayer,jCell) == 1 ) then
if ( TFoceanOld(iLayer,jCell) > 1.0e6_RKIND) then
! raise error if an invalid ocean data value is used
call mpas_log_write("ocean data value used for extrapolation is invalid", &
MPAS_LOG_ERR)
call mpas_log_write("ocean data value used for extrapolation is invalid " // &
"in horizontal_extrapolation: cell id=$i, iLayer=$i, TF=$r", &
MPAS_LOG_ERR, intArgs=(/indexToCellID(jCell), iLayer/), realArgs=(/TFoceanOld(iLayer,jCell)/))
err = ior(err,1)
else
TFsum = TFsum + (TFoceanOld(iLayer,jCell) * areaCell(jCell))
Expand Down Expand Up @@ -533,6 +536,7 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas
real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography, areaCell
integer, pointer :: nCells, nCellsSolve, nISMIP6OceanLayers
integer, dimension(:), pointer :: cellMask, nEdgesOnCell
integer, dimension(:), pointer :: indexToCellID
integer, dimension(:,:), pointer :: cellsOnCell
integer :: iCell, jCell, iLayer, iNeighbor, iter
integer :: localLoopCount, newMaskCountLocalAccum
Expand All @@ -550,6 +554,7 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas
call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell)
call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
Expand All @@ -564,8 +569,9 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas
if ( validOceanMask(iLayer-1,iCell) == 1 ) then
if (TFocean(iLayer-1,iCell) > 1.0e6_RKIND) then
! raise error if an invalid ocean data value is used
call mpas_log_write("ocean data value used for extrapolation is invalid", &
MPAS_LOG_ERR)
call mpas_log_write("ocean data value used for extrapolation is invalid " // &
"in vertical_extrapolation: cell id=$i, iLayer=$i, TF=$r", &
MPAS_LOG_ERR, intArgs=(/indexToCellID(iCell), iLayer-1/), realArgs=(/TFocean(iLayer-1,iCell)/))
err = ior(err,1)
else
TFocean(iLayer,iCell) = TFocean(iLayer-1,iCell) - &
Expand Down