Skip to content

Commit

Permalink
Merge branch 'matthewhoffman/mali/tf-vert-interp-fix' into MALI-Dev/d…
Browse files Browse the repository at this point in the history 
…evelop

The ISMIP6 code for interpolating TF vertically does not give correct
values for depths below the deepest data z-level. This commit handles
that case by using the deepest TF value corrected by the freezing
temperature depth-dependence. It also handles the less critical
situation of depths above the shallowest data level by using the
shallowest value unmodified.

* origin/matthewhoffman/mali/tf-vert-interp-fix:
  Update error message
  Update zOcean error check
  Improve error checking in ocean extrap
  Fix new indexing error introduced trying to fix indexing errors
  Fix indexing errors
  Modify TFocean vertical interpolation for depths outside of data range
  • Loading branch information
@matthewhoffman
matthewhoffman committed Jan 7, 2025
2 parents 8fa61a6 + 2f684b3 commit 37848c5
Show file tree
Hide file tree
Showing 2 changed files with 57 additions and 22 deletions.
63 changes: 46 additions & 17 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F
View file
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 >= 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 non-positive 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
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

0 comments on commit 37848c5

Please sign in to comment.