diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F index 37d6578dd331..daa1d5cca87d 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_iceshelf_melt.F @@ -1162,6 +1162,8 @@ end subroutine calc_iceshelf_draft_info subroutine iceshelf_melt_ismip6(domain, err) + use li_constants, only: oceanFreezingTempDepthDependence + !----------------------------------------------------------------- ! input/output variables !----------------------------------------------------------------- @@ -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) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_ocean_extrap.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_ocean_extrap.F index 024c0df4c5f5..de501d629605 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_ocean_extrap.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_ocean_extrap.F @@ -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 @@ -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 @@ -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) @@ -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)) @@ -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 @@ -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) @@ -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) - &