From 0d7e8e09ac63c5e008b7de529a872c81a888f148 Mon Sep 17 00:00:00 2001 From: Courtney Shafer Date: Mon, 26 Aug 2024 14:54:04 -0600 Subject: [PATCH 1/8] Modified ocean extrapolation routine to take in a 3d ocean cavity mask --- .../mpas-albany-landice/src/Registry.xml | 8 +++++- .../src/mode_forward/mpas_li_ocean_extrap.F | 25 ++++++++++++++++--- 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 9ce37b56fdfc..4dac552d19b4 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -810,6 +810,7 @@ + - + + @@ -1689,6 +1691,10 @@ is the value of that variable from the *previous* time level! + + 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 16bc4eddd94b..1ac3b8cb63f1 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 @@ -91,6 +91,7 @@ subroutine li_ocean_extrap_solve(domain, err) real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing, ismip6shelfMelt_zBndsOcean real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography integer, dimension(:), pointer :: origOceanMaskHoriz + integer, dimension(:,:), pointer :: orig3dOceanCavityMask ! CAS 8/16/2024 integer, dimension(:,:), pointer :: validOceanMask, validOceanMaskOrig, availOceanMask !masks to pass to flood-fill routine integer, dimension(:), pointer :: seedOceanMaskHoriz, growOceanMaskHoriz, seedOceanMaskHorizInit integer, dimension(:), allocatable :: seedOceanMaskHorizOld @@ -127,6 +128,7 @@ subroutine li_ocean_extrap_solve(domain, err) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing) call mpas_pool_get_array(extrapOceanDataPool, 'origOceanMaskHoriz', origOceanMaskHoriz) + call mpas_pool_get_array(extrapOceanDataPool, 'orig3dOceanCavityMask', orig3dOceanCavityMask) ! CAS 8/16/2024 call mpas_pool_get_array(extrapOceanDataPool, 'validOceanMask', validOceanMask) call mpas_pool_get_array(extrapOceanDataPool, 'validOceanMaskOrig', validOceanMaskOrig) call mpas_pool_get_array(extrapOceanDataPool, 'availOceanMask', availOceanMask) @@ -189,11 +191,25 @@ subroutine li_ocean_extrap_solve(domain, err) if ( (seedOceanMaskHoriz(iCell) == 1) .and. (bedTopography(iCell) < layerTop) ) then availOceanMask(iLayer,iCell) = 1 endif - if ( (origOceanMaskHoriz(iCell) == 1) .and. (bedTopography(iCell) < layerTop) ) then - validOceanMask(iLayer,iCell) = 1 - endif enddo enddo + + ! CAS 8/16/2024 Hijacking Holly's code to test using the 3D SORRM cavity TF field. We set 3d valid ocean mask to original 3d ocean cavity mask. + ! We don't need to loop through the layers if we're using a 3d mask already + validOceanMask(:,:) = orig3dOceanCavityMask(:,:) + + ! if ( (origOceanMaskHoriz(iCell) == 1) .and. (bedTopography(iCell) < layerTop) ) then + ! validOceanMask(iLayer,iCell) = 1 + ! endif + ! enddo + !enddo + + call mpas_log_write('==HH==: updating halos for the avail/valid ocean masks') + ! Update halos + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'availOceanMask') + call mpas_dmpar_field_halo_exch(domain, 'validOceanMask') + call mpas_timer_stop("halo updates") ! save the initial validOceanMask validOceanMaskOrig(:,:) = validOceanMask(:,:) @@ -201,9 +217,10 @@ subroutine li_ocean_extrap_solve(domain, err) ! initialize the TF field TFocean(:,:) = ismip6shelfMelt_3dThermalForcing(:,:) * validOceanMask(:,:) ! initialize the invalid data locations with fill value + ! CAS 8/19/2024 Changed availOceanMask in the if check to validOceanMask. We want to make sure all invalid values outside the validOceanMask are set to the invalid_value_TF do iCell = 1, nCellsSolve do iLayer = 1, nISMIP6OceanLayers - if ( availOceanMask(iLayer,iCell) == 0 ) then + if ( validOceanMask(iLayer,iCell) == 0 ) then TFocean(iLayer,iCell) = invalid_value_TF endif enddo From 93c7ad7843d6339c484a4e8c09dc40f1782fef83 Mon Sep 17 00:00:00 2001 From: Courtney Shafer Date: Mon, 26 Aug 2024 15:02:08 -0600 Subject: [PATCH 2/8] Modified which ocean layers are picked for creating TFocean and added check for invalid TF values --- .../src/mode_forward/mpas_li_iceshelf_melt.F | 33 +++++++++++++++---- 1 file changed, 27 insertions(+), 6 deletions(-) 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..b24683b9aad0 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 @@ -1386,6 +1386,7 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & !----------------------------------------------------------------- integer :: iCell, iNeighbor, iEdge, nEmptyNeighbors real (kind=RKIND), pointer :: rhoi + real (kind=RKIND), pointer :: invalid_value_TF ! CAS 8/20/2024 integer, pointer :: nCells, nISMIP6OceanLayers integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell real (kind=RKIND) :: waterDepth @@ -1402,7 +1403,7 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & real (kind=RKIND), dimension(:), pointer :: ismip6shelfMelt_deltaT real (kind=RKIND), dimension(:), pointer :: zOcean real (kind=RKIND), dimension(:), pointer :: areaCell - integer, dimension(:), pointer :: cellMask, edgeMask, nEdgesOnCell + integer, dimension(:), pointer :: cellMask, edgeMask, nEdgesOnCell, indexToCellID ! CAS 8/22/2024 real (kind=RKIND), pointer :: aSubglacial ! param A real (kind=RKIND), pointer :: alphaSubglacial ! param alpha real (kind=RKIND), pointer :: B ! param B @@ -1425,6 +1426,7 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & ! Get sea level, bedTopography, ice density call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi) call mpas_pool_get_config(liConfigs, 'config_sea_level', seaLevel) + call mpas_pool_get_config(liConfigs, 'config_invalid_value_TF', invalid_value_TF) ! CAS 8/20/2024 ! Get melt parameters call mpas_pool_get_config(liConfigs, 'config_beta_ocean_thermal_forcing', betaTF) @@ -1452,6 +1454,7 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) + call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID) ! CAS 8/22/2024 call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) call mpas_pool_get_array(geometryPool, 'faceMeltingThickness', faceMeltingThickness) @@ -1461,6 +1464,8 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & call mpas_pool_get_array(geometryPool, 'dtFaceMeltingCFLratio', dtFaceMeltingCFLratio) if ( config_use_3d_thermal_forcing_for_face_melt ) then + call mpas_log_write("config_use_3d_thermal_forcing_for_face_melt is .true.") + call mpas_pool_get_dimension(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zOcean', zOcean) @@ -1477,14 +1482,30 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & enddo ! If bed is shallower than first layer, use TF from the first layer. ! If bed is deeper than the bottomg ocean layer, use TF from the bottom layer. - if ( (kk == 1) .or. ( (kk == nISMIP6OceanLayers) .and. (zOcean(kk) >= bedTopography(iCell)) ) ) then - TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk, iCell) + ismip6shelfMelt_deltaT(iCell) + !if ( (kk == 1) .or. ( (kk == nISMIP6OceanLayers) .and. (zOcean(kk) >= bedTopography(iCell)) ) ) then + ! TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk, iCell) + ismip6shelfMelt_deltaT(iCell) ! For all other bed depths, interpolate linearly between layers above and below bed depth. + !else + ! TFocean(iCell) = ( (zOcean(kk-1) - bedTopography(iCell)) * ismip6shelfMelt_3dThermalForcing(kk, iCell) + & + ! (bedTopography(iCell) - zOcean(kk)) * ismip6shelfMelt_3dThermalForcing(kk-1, iCell) ) / & + ! (zOcean(kk-1) - zOcean(kk)) + ismip6shelfMelt_deltaT(iCell) + !endif + + if ( kk == 1 ) then + TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk, iCell) + ismip6shelfMelt_deltaT(iCell) else - TFocean(iCell) = ( (zOcean(kk-1) - bedTopography(iCell)) * ismip6shelfMelt_3dThermalForcing(kk, iCell) + & - (bedTopography(iCell) - zOcean(kk)) * ismip6shelfMelt_3dThermalForcing(kk-1, iCell) ) / & - (zOcean(kk-1) - zOcean(kk)) + ismip6shelfMelt_deltaT(iCell) + TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk-1, iCell) + ismip6shelfMelt_deltaT(iCell) endif + + ! check if any invalid TFocean value is used for calculating face melting speed + if ( (TFocean(iCell) == invalid_value_TF) .or. (TFocean(iCell) > 1.0e1) ) then + call mpas_log_write("grounded_face_melt_ismip6: Invalid value for TFocean. " // & + "TFocean(iCell)=$r zOcean(kk)=$r bedTopography(iCell)=$r kk=$i indexToCellID=$i", MPAS_LOG_ERR, & + intArgs=(/kk, indexToCellID(iCell)/), & + realArgs=(/TFocean(iCell), zOcean(kk), bedTopography(iCell)/) ) + err = ior(err, 1) + endif + endif end do endif From 932736e57dd0179afe997cc46b3d682100f6921b Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 28 Aug 2024 10:51:32 -0600 Subject: [PATCH 3/8] Revert "Modified which ocean layers are picked for creating TFocean and added check for invalid TF values" This reverts commit 93c7ad7843d6339c484a4e8c09dc40f1782fef83. --- .../src/mode_forward/mpas_li_iceshelf_melt.F | 33 ++++--------------- 1 file changed, 6 insertions(+), 27 deletions(-) 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 b24683b9aad0..37d6578dd331 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 @@ -1386,7 +1386,6 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & !----------------------------------------------------------------- integer :: iCell, iNeighbor, iEdge, nEmptyNeighbors real (kind=RKIND), pointer :: rhoi - real (kind=RKIND), pointer :: invalid_value_TF ! CAS 8/20/2024 integer, pointer :: nCells, nISMIP6OceanLayers integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell real (kind=RKIND) :: waterDepth @@ -1403,7 +1402,7 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & real (kind=RKIND), dimension(:), pointer :: ismip6shelfMelt_deltaT real (kind=RKIND), dimension(:), pointer :: zOcean real (kind=RKIND), dimension(:), pointer :: areaCell - integer, dimension(:), pointer :: cellMask, edgeMask, nEdgesOnCell, indexToCellID ! CAS 8/22/2024 + integer, dimension(:), pointer :: cellMask, edgeMask, nEdgesOnCell real (kind=RKIND), pointer :: aSubglacial ! param A real (kind=RKIND), pointer :: alphaSubglacial ! param alpha real (kind=RKIND), pointer :: B ! param B @@ -1426,7 +1425,6 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & ! Get sea level, bedTopography, ice density call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi) call mpas_pool_get_config(liConfigs, 'config_sea_level', seaLevel) - call mpas_pool_get_config(liConfigs, 'config_invalid_value_TF', invalid_value_TF) ! CAS 8/20/2024 ! Get melt parameters call mpas_pool_get_config(liConfigs, 'config_beta_ocean_thermal_forcing', betaTF) @@ -1454,7 +1452,6 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) - call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID) ! CAS 8/22/2024 call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) call mpas_pool_get_array(geometryPool, 'faceMeltingThickness', faceMeltingThickness) @@ -1464,8 +1461,6 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & call mpas_pool_get_array(geometryPool, 'dtFaceMeltingCFLratio', dtFaceMeltingCFLratio) if ( config_use_3d_thermal_forcing_for_face_melt ) then - call mpas_log_write("config_use_3d_thermal_forcing_for_face_melt is .true.") - call mpas_pool_get_dimension(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zOcean', zOcean) @@ -1482,30 +1477,14 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & enddo ! If bed is shallower than first layer, use TF from the first layer. ! If bed is deeper than the bottomg ocean layer, use TF from the bottom layer. - !if ( (kk == 1) .or. ( (kk == nISMIP6OceanLayers) .and. (zOcean(kk) >= bedTopography(iCell)) ) ) then - ! TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk, iCell) + ismip6shelfMelt_deltaT(iCell) - ! For all other bed depths, interpolate linearly between layers above and below bed depth. - !else - ! TFocean(iCell) = ( (zOcean(kk-1) - bedTopography(iCell)) * ismip6shelfMelt_3dThermalForcing(kk, iCell) + & - ! (bedTopography(iCell) - zOcean(kk)) * ismip6shelfMelt_3dThermalForcing(kk-1, iCell) ) / & - ! (zOcean(kk-1) - zOcean(kk)) + ismip6shelfMelt_deltaT(iCell) - !endif - - if ( kk == 1 ) then + if ( (kk == 1) .or. ( (kk == nISMIP6OceanLayers) .and. (zOcean(kk) >= bedTopography(iCell)) ) ) then TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk, iCell) + ismip6shelfMelt_deltaT(iCell) + ! For all other bed depths, interpolate linearly between layers above and below bed depth. else - TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk-1, iCell) + ismip6shelfMelt_deltaT(iCell) + TFocean(iCell) = ( (zOcean(kk-1) - bedTopography(iCell)) * ismip6shelfMelt_3dThermalForcing(kk, iCell) + & + (bedTopography(iCell) - zOcean(kk)) * ismip6shelfMelt_3dThermalForcing(kk-1, iCell) ) / & + (zOcean(kk-1) - zOcean(kk)) + ismip6shelfMelt_deltaT(iCell) endif - - ! check if any invalid TFocean value is used for calculating face melting speed - if ( (TFocean(iCell) == invalid_value_TF) .or. (TFocean(iCell) > 1.0e1) ) then - call mpas_log_write("grounded_face_melt_ismip6: Invalid value for TFocean. " // & - "TFocean(iCell)=$r zOcean(kk)=$r bedTopography(iCell)=$r kk=$i indexToCellID=$i", MPAS_LOG_ERR, & - intArgs=(/kk, indexToCellID(iCell)/), & - realArgs=(/TFocean(iCell), zOcean(kk), bedTopography(iCell)/) ) - err = ior(err, 1) - endif - endif end do endif From 5aa744cb63cc22b4f0ce9fc8be8f0bc3513eea33 Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 28 Aug 2024 21:22:02 -0600 Subject: [PATCH 4/8] Handle seafloor and lakes in extrapolation module This commit makes adjustments to how the seafloor and inland seas (locations below sea level not connected to the global ocean) are handled for extraplation. The adjustments are applied primarily to availOceanMask, which is the mask of to where ocean data should be extrapolated. Changes made: * Adjust availOceanMask to extend one layer below the seafloor (needed so that facemelting has a valid value one level below the seafloor for interpolation to the seafloor elevation) * to support this, generate error if ismip6shelfMelt_zBndsOcean is not populated, because that field is needed for the seafloor detection * Adjust availOceanMask to ignore inland seas - we will not attempt to flood fill into areas below sea level that are not connected to the open ocean. * to support this, create mask of marine locations connected to global open ocean * Adjustment to where invalid values are assigned to avoid inserting them in inland seas. This will give inland seas either the value from the ocean data if it exists, or else TF=0 --- .../src/mode_forward/Makefile | 2 +- .../src/mode_forward/mpas_li_ocean_extrap.F | 97 +++++++++++++++---- 2 files changed, 78 insertions(+), 21 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/Makefile b/components/mpas-albany-landice/src/mode_forward/Makefile index 2ca4ce74b36d..3c9e0061ec91 100644 --- a/components/mpas-albany-landice/src/mode_forward/Makefile +++ b/components/mpas-albany-landice/src/mode_forward/Makefile @@ -82,7 +82,7 @@ mpas_li_velocity_external.o: Interface_velocity_solver.o mpas_li_bedtopo.o: mpas_li_advection.o -mpas_li_ocean_extrap.o: +mpas_li_ocean_extrap.o: mpas_li_calving.o Interface_velocity_solver.o: 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 1ac3b8cb63f1..17e8e57688ec 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 @@ -64,6 +64,7 @@ module li_ocean_extrap !----------------------------------------------------------------------- subroutine li_ocean_extrap_solve(domain, err) + use li_calving, only: li_flood_fill !----------------------------------------------------------------- ! input variables @@ -89,6 +90,7 @@ subroutine li_ocean_extrap_solve(domain, err) real (kind=RKIND) :: layerTop real (kind=RKIND), dimension(:,:), pointer :: TFocean, TFoceanOld real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing, ismip6shelfMelt_zBndsOcean + real (kind=RKIND), dimension(:), pointer :: ismip6shelfMelt_zOcean real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography integer, dimension(:), pointer :: origOceanMaskHoriz integer, dimension(:,:), pointer :: orig3dOceanCavityMask ! CAS 8/16/2024 @@ -100,6 +102,9 @@ subroutine li_ocean_extrap_solve(domain, err) integer, dimension(:,:), pointer :: cellsOnCell integer :: iCell, jCell, iLayer, iNeighbor, iter, err_tmp integer :: GlobalLoopCount, newMaskCountGlobal + type (field1dInteger), pointer :: seedMaskField + type (field1dInteger), pointer :: growMaskField + integer, dimension(:), pointer :: connectedMarineMask, growMask !masks to pass to flood-fill routine ! No init is needed. err = 0 @@ -180,36 +185,79 @@ subroutine li_ocean_extrap_solve(domain, err) enddo deallocate(seedOceanMaskHorizOld) + ! Calculate mask of connected ocean + call mpas_pool_get_subpool(domain % blocklist % structs, 'scratch', scratchPool) + call mpas_pool_get_field(scratchPool, 'seedMask', seedMaskField) + call mpas_allocate_scratch_field(seedMaskField, single_block_in = .true.) + connectedMarineMask => seedMaskField % array + connectedmarineMask(:) = 0 + call mpas_pool_get_field(scratchPool, 'growMask', growMaskField) + call mpas_allocate_scratch_field(growMaskField, single_block_in = .true.) + growMask => growMaskField % array + growMask(:) = 0 + + do iCell = 1, nCells + ! seedMask = open ocean cells in contact with the domain boundary + if ((bedTopography(iCell) < config_sea_level) .and. (thickness(iCell) == 0.0_RKIND)) then + do iNeighbor = 1, nEdgesOnCell(iCell) + if (cellsOnCell(iNeighbor, iCell) == nCells + 1) then + connectedMarineMask(iCell) = 1 + exit ! no need to keep checking neighbors + endif + enddo + endif + ! growMask - all marine cells + if (bedTopography(iCell) < config_sea_level) then + growMask(iCell) = 1 + endif + enddo + ! now create mask of all marine locations connected to open ocean - to be used below to screen out lakes + call li_flood_fill(connectedMarineMask, growMask, domain) + ! make it a 3D mask based on the topography (loop through nISMIP6OceanLayers) call mpas_pool_get_dimension(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers) + call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zOcean', ismip6shelfMelt_zOcean) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zBndsOcean', ismip6shelfMelt_zBndsOcean) + ! check for valid data + do iLayer = 1, nISMIP6OceanLayers + if (ismip6shelfMelt_zOcean(iLayer) >= 0.0_RKIND) then + call mpas_log_write("ismip6shelfMelt_zOcean has invalid value of $r in layer $i", MPAS_LOG_ERR, & + realArgs=(/ismip6shelfMelt_zOcean(iLayer)/), intArgs=(/iLayer/)) + err = ior(err, 1) + endif + if ((ismip6shelfMelt_zBndsOcean(1,iLayer) > 0.0_RKIND) .or. & + (ismip6shelfMelt_zBndsOcean(1,iLayer) < ismip6shelfMelt_zOcean(iLayer))) then + call mpas_log_write("ismip6shelfMelt_zBndsOcean(1,:) has invalid value of $r in layer $i", MPAS_LOG_ERR, & + realArgs=(/ismip6shelfMelt_zBndsOcean(1,iLayer)/), intArgs=(/iLayer/)) + err = ior(err, 1) + endif + if ((ismip6shelfMelt_zBndsOcean(2,iLayer) >= 0.0_RKIND) .or. & + (ismip6shelfMelt_zBndsOcean(2,iLayer) > ismip6shelfMelt_zOcean(iLayer))) then + call mpas_log_write("ismip6shelfMelt_zBndsOcean(2,:) has invalid value of $r in layer $i", MPAS_LOG_ERR, & + realArgs=(/ismip6shelfMelt_zBndsOcean(2,iLayer)/), intArgs=(/iLayer/)) + err = ior(err, 1) + endif + enddo availOceanMask(:,:) = 0 validOceanMask(:,:) = 0 do iCell = 1, nCells do iLayer = 1, nISMIP6OceanLayers layerTop = ismip6shelfMelt_zBndsOcean(1, iLayer) - if ( (seedOceanMaskHoriz(iCell) == 1) .and. (bedTopography(iCell) < layerTop) ) then - availOceanMask(iLayer,iCell) = 1 + if ( (seedOceanMaskHoriz(iCell) == 1) .and. (connectedMarineMask(iCell) == 1)) then + if (bedTopography(iCell) < layerTop) then + availOceanMask(iLayer,iCell) = 1 + else + ! keep the first layer below the seafloor in the region to be filled + ! this ensures linear interpolation from above and below the seafloor is possible + availOceanMask(iLayer,iCell) = 1 + exit ! stop looping over levels after we've included the first level below the seafloor + endif endif enddo enddo - - ! CAS 8/16/2024 Hijacking Holly's code to test using the 3D SORRM cavity TF field. We set 3d valid ocean mask to original 3d ocean cavity mask. - ! We don't need to loop through the layers if we're using a 3d mask already - validOceanMask(:,:) = orig3dOceanCavityMask(:,:) - ! if ( (origOceanMaskHoriz(iCell) == 1) .and. (bedTopography(iCell) < layerTop) ) then - ! validOceanMask(iLayer,iCell) = 1 - ! endif - ! enddo - !enddo - - call mpas_log_write('==HH==: updating halos for the avail/valid ocean masks') - ! Update halos - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'availOceanMask') - call mpas_dmpar_field_halo_exch(domain, 'validOceanMask') - call mpas_timer_stop("halo updates") + ! CAS 8/16/2024 Hijacking Holly's code to test using the 3D SORRM cavity TF field. We set 3d valid ocean mask to original 3d ocean cavity mask. + validOceanMask(:,:) = orig3dOceanCavityMask(:,:) ! save the initial validOceanMask validOceanMaskOrig(:,:) = validOceanMask(:,:) @@ -217,15 +265,24 @@ subroutine li_ocean_extrap_solve(domain, err) ! initialize the TF field TFocean(:,:) = ismip6shelfMelt_3dThermalForcing(:,:) * validOceanMask(:,:) ! initialize the invalid data locations with fill value - ! CAS 8/19/2024 Changed availOceanMask in the if check to validOceanMask. We want to make sure all invalid values outside the validOceanMask are set to the invalid_value_TF do iCell = 1, nCellsSolve do iLayer = 1, nISMIP6OceanLayers - if ( validOceanMask(iLayer,iCell) == 0 ) then + if ((connectedMarineMask(iCell) == 0) .and. (bedTopography(iCell) < config_sea_level)) then + ! Don't assign invalid value to lakes/inland seas disconnected from global ocean + ! Let them retain the existing value: + ! This will take on the valid ocean data value where it exists or + ! zero where valid ocean data does not exist + elseif (validOceanMask(iLayer,iCell) == 0) then + ! everywhere else where valid ocean data does not exist, insert invalid value outside of validOceanMask TFocean(iLayer,iCell) = invalid_value_TF endif enddo enddo + ! deallocate scratch fields used for flood fill + call mpas_deallocate_scratch_field(seedMaskField, single_block_in=.true.) + call mpas_deallocate_scratch_field(growMaskField, single_block_in=.true.) + ! flood-fill the valid ocean mask and TF field through ! horizontal and vertial extrapolation ! get initial 3D valid data based on the original ISMIP6 field From 75c194029391fd607699a96510b81b68a2832b8b Mon Sep 17 00:00:00 2001 From: Matthew Hoffman Date: Wed, 28 Aug 2024 21:33:03 -0600 Subject: [PATCH 5/8] Replace verbose per-iteration log message with a summary Eliminates 100s of repetitive log messages per timestep --- .../src/mode_forward/mpas_li_ocean_extrap.F | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) 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 17e8e57688ec..0677b9c63c67 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 @@ -384,6 +384,7 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali integer :: iCell, jCell, iLayer, iNeighbor, iter integer :: localLoopCount integer :: nValidNeighb, newValidCount, newMaskCountLocalAccum, newMaskCountGlobal + integer :: newMaskCountTotal err = 0 @@ -413,6 +414,7 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali ! initialize the local loop and count for validOceanMask localLoopCount = 0 + newMaskCountTotal = 0 newMaskCountGlobal = 1 call mpas_log_write('Weight given to the cell with valid data from extrapolation: $r', realArgs=(/weightCell/)) do while ( newMaskCountGlobal > 0 ) @@ -474,16 +476,17 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali ! update count of cells added to mask globally call mpas_dmpar_sum_int(domain % dminfo, newMaskCountLocalAccum, newMaskCountGlobal) - call mpas_log_write('Horizontal extrap: Added total $i new cells to validOceanMask', intArgs=(/newMaskCountGlobal/)) + newMaskCountTotal = newMaskCountTotal + newMaskCountGlobal + !call mpas_log_write('Horizontal extrap: Added total $i new cells to validOceanMask', intArgs=(/newMaskCountGlobal/)) enddo - call mpas_log_write('Horizontal extrapolation done after $i loops', intArgs=(/localLoopCount/)) + call mpas_log_write('Horizontal extrapolation done after $i iterations. Added total of $i cells across all processors', & + intArgs=(/localLoopCount, newMaskCountTotal/)) deallocate(validOceanMaskOld) - - end subroutine horizontal_extrapolation !----------------------------------------------------------------------- + !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ! ! routine vertical_extrapolation From bba3a427e65e4e8df3de7807bee04dc3e23b1823 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 5 Sep 2024 17:42:51 -0600 Subject: [PATCH 6/8] Add salt + temp to ocean extrapolation Makes changes to mpas_li_ocean_extrap.F to enable the choice to either extrapolate thermal forcing (original implementation) or extrapolate both ocean salinity and temperature. In the latter case, thermal forcing will eventually be recalculated after the extrapolation process. The option to extrapolate salt/temp is activated by setting config_recalculate_thermal_forcing to true. --- .../mpas-albany-landice/src/Registry.xml | 24 +- .../src/mode_forward/mpas_li_ocean_extrap.F | 217 +++++++++++------- 2 files changed, 144 insertions(+), 97 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 4dac552d19b4..ac90339fdad2 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -454,16 +454,16 @@ possible_values="any non-negative value" /> - + + - - - @@ -1713,11 +1711,11 @@ is the value of that variable from the *previous* time level! - - 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 0677b9c63c67..8e742aaeee8b 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 @@ -11,11 +11,11 @@ ! li_ocean_extrap ! !> \MPAS land-ice ocean-data extrapolation driver -!> \author Holly Han -!> \date January 2022 +!> \author Holly Han, modified by Alex Hager +!> \date January 2022 (modified September 2024) !> \details !> This module contains the routines for extrapolating -!> ocean data (e.g., temperature, salinity, thermal forcing) +!> ocean data (e.g., temperature and salinity, or thermal forcing) !> into ice draft ! !----------------------------------------------------------------------- @@ -45,6 +45,8 @@ module li_ocean_extrap !-------------------------------------------------------------------- ! Private module variables !-------------------------------------------------------------------- + + real (kind=RKIND), parameter :: invalid_value = 1.0e36_RKIND !*********************************************************************** @@ -84,11 +86,10 @@ subroutine li_ocean_extrap_solve(domain, err) ! local variables !----------------------------------------------------------------- logical, pointer :: config_ocean_data_extrapolation - real(kind=RKIND), pointer :: config_sea_level, invalid_value_TF + real(kind=RKIND), pointer :: config_sea_level type (block_type), pointer :: block type (mpas_pool_type), pointer :: scratchPool, geometryPool, meshPool, extrapOceanDataPool real (kind=RKIND) :: layerTop - real (kind=RKIND), dimension(:,:), pointer :: TFocean, TFoceanOld real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing, ismip6shelfMelt_zBndsOcean real (kind=RKIND), dimension(:), pointer :: ismip6shelfMelt_zOcean real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography @@ -105,6 +106,7 @@ subroutine li_ocean_extrap_solve(domain, err) type (field1dInteger), pointer :: seedMaskField type (field1dInteger), pointer :: growMaskField integer, dimension(:), pointer :: connectedMarineMask, growMask !masks to pass to flood-fill routine + real (kind=RKIND), dimension(:,:,:), allocatable :: oceanField ! No init is needed. err = 0 @@ -120,6 +122,7 @@ subroutine li_ocean_extrap_solve(domain, err) ! initialize the ocean data and mask fields call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) call mpas_pool_get_config(liConfigs, 'config_ocean_data_extrap_ncells_extra', nCellsExtra) + call mpas_pool_get_config(liConfigs, 'config_recalculate_thermal_forcing', config_recalcuate_thermal_forcing) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'extrapOceanData', extrapOceanDataPool) @@ -140,9 +143,8 @@ subroutine li_ocean_extrap_solve(domain, err) call mpas_pool_get_array(extrapOceanDataPool, 'seedOceanMaskHoriz', seedOceanMaskHoriz) call mpas_pool_get_array(extrapOceanDataPool, 'growOceanMaskHoriz', growOceanMaskHoriz) call mpas_pool_get_array(extrapOceanDataPool, 'seedOceanMaskHorizInit', seedOceanMaskHorizInit) - call mpas_pool_get_array(extrapOceanDataPool, 'TFoceanOld', TFoceanOld) - call mpas_pool_get_array(extrapOceanDataPool, 'TFocean', TFocean) - call mpas_pool_get_config(liConfigs, 'config_invalid_value_TF', invalid_value_TF) + + allocate(oceanField(nISMIP6OceanLayers,nCells+1,2)) ! create a 2D mask based on the open ocean + floating ice + grounded ice mask to be used for defining the location of the n-extra cells into the grounded ice allocate(seedOceanMaskHorizOld(nCells+1)) @@ -165,7 +167,7 @@ subroutine li_ocean_extrap_solve(domain, err) do iter = 1, nCellsExtra do iCell = 1, nCellsSolve if ( growOceanMaskHoriz(iCell) == 1 .and. seedOceanMaskHorizOld(iCell) == 0 ) then - do iNeighbor = 1, nEdgesOnCell(iCell) + do iNeighbor /TFoce= 1, nEdgesOnCell(iCell) jCell = cellsOnCell(iNeighbor, iCell) if ( seedOceanMaskHorizOld(jCell) == 1 ) then seedOceanMaskHoriz(iCell) = 1 @@ -262,8 +264,22 @@ subroutine li_ocean_extrap_solve(domain, err) ! save the initial validOceanMask validOceanMaskOrig(:,:) = validOceanMask(:,:) - ! initialize the TF field - TFocean(:,:) = ismip6shelfMelt_3dThermalForcing(:,:) * validOceanMask(:,:) + ! initialize the extrapolated fields + if (config_recalculate_thermal_forcing) then + + real (kind=RKIND), dimension(:,:), pointer :: MPAS_3dOceanTemperature, MPAS_3dOceanSalinity + call mpas_pool_get_array(extrapOceanDataPool, 'MPAS_3dOceanTemperature', oceanTemp) + call mpas_pool_get_array(extrapOceanDataPool, 'MPAS_3dOceanSalinity', oceanSalt) + + !combine fields into one variable to expedite extrapolation + oceanField(:,:,1) = MPAS_3dOceanTemperature(:,:) * validOceanMask(:,:) + oceanField(:,:,2) = MPAS_3dOceanSalinity(:,:) * validOceanMask(:,:) + else + + oceanField(:,:,1) = ismip6shelfMelt_3dThermalForcing * validOceanMask(:,:)i + oceanField(:,:,2) = invalid_value + endif + ! initialize the invalid data locations with fill value do iCell = 1, nCellsSolve do iLayer = 1, nISMIP6OceanLayers @@ -274,7 +290,7 @@ subroutine li_ocean_extrap_solve(domain, err) ! zero where valid ocean data does not exist elseif (validOceanMask(iLayer,iCell) == 0) then ! everywhere else where valid ocean data does not exist, insert invalid value outside of validOceanMask - TFocean(iLayer,iCell) = invalid_value_TF + oceanField(iLayer,iCell,:) = invalid_value endif enddo enddo @@ -300,23 +316,31 @@ subroutine li_ocean_extrap_solve(domain, err) endif ! call the horizontal extrapolation routine call mpas_timer_start("horizontal scheme") - call horizontal_extrapolation(domain, availOceanMask, validOceanMask, validOceanMaskOrig, TFocean, err_tmp) + call horizontal_extrapolation(domain, availOceanMask, validOceanMask, validOceanMaskOrig, oceanField, err_tmp) err = ior(err, err_tmp) call mpas_timer_stop("horizontal scheme") ! call the vertical extrapolation routine call mpas_timer_start("vertical scheme") - call vertical_extrapolation(domain, availOceanMask, validOceanMask, newMaskCountGlobal, TFocean, err_tmp) + call vertical_extrapolation(domain, availOceanMask, validOceanMask, newMaskCountGlobal, oceanField, err_tmp) err = ior(err, err_tmp) 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 ! Reassign extrapolated TF back to primary TF field - ismip6shelfMelt_3dThermalForcing(:,:) = TFocean(:,:) + if (config_recalculate_thermal_forcing) then + MPAS_3dOceanTemperature(:,:) = oceanField(:,:,1) + MPAS_3dOceanSalinity(:,:) = oceanField(:,:,2) + else + ismip6shelfMelt_3dThermalForcing(:,:) = oceanField(:,:,1) + endif + + deallocate(oceanField) + endif !-------------------------------------------------------------------- @@ -337,8 +361,8 @@ end subroutine li_ocean_extrap_solve ! ! routine horizontal_extrapolation ! !> \brief Extrapolate validOceanMask horizontally -!> \author Holly Kyeore Han -!> \date November 2023 +!> \author Holly Kyeore Han, modified by Alex Hager +!> \date November 2023 (modified September 2024) !> \details !> This routine extrapolates takes the initialized availOceanMask !> and validOceanMask and extrapolates validOceanMask in horizontal @@ -348,7 +372,7 @@ end subroutine li_ocean_extrap_solve !----------------------------------------------------------------------- - subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, validOceanMaskOrig, TFocean, err) + subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, validOceanMaskOrig, oceanField, err) !----------------------------------------------------------------- ! input variables @@ -361,7 +385,7 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali type (domain_type), intent(inout) :: domain !< Input/Output: domain object integer, dimension(:,:), pointer, intent(inout) :: validOceanMask integer, intent(inout) :: err !< Output: error flag - real (kind=RKIND), dimension(:,:), pointer, intent(inout) :: TFocean + real (kind=RKIND), dimension(:,:), pointer, intent(inout) :: oceanField !----------------------------------------------------------------- ! output variables @@ -372,19 +396,21 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali !----------------------------------------------------------------- type (block_type), pointer :: block type (mpas_pool_type), pointer :: scratchPool, geometryPool, meshPool, extrapOceanDataPool - real (kind=RKIND) :: layerTop, TFsum, areaSum, weightCellLocal + real (kind=RKIND) :: layerTop, fieldSum, areaSum, weightCellLocal real (kind=RKIND), pointer :: weightCell integer, dimension(:,:), allocatable :: validOceanMaskOld real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing - real (kind=RKIND), dimension(:,:), pointer :: TFoceanOld + real (kind=RKIND), dimension(:,:,:), pointer :: oceanFieldOld real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography, areaCell integer, pointer :: nCells, nCellsSolve, nISMIP6OceanLayers, nCellsExtra integer, dimension(:), pointer :: cellMask, nEdgesOnCell integer, dimension(:,:), pointer :: cellsOnCell - integer :: iCell, jCell, iLayer, iNeighbor, iter + integer :: iCell, jCell, iLayer, iNeighbor, iter, iField integer :: localLoopCount integer :: nValidNeighb, newValidCount, newMaskCountLocalAccum, newMaskCountGlobal integer :: newMaskCountTotal + logical, pointer :: config_recalculate_thermal_forcing + real (kind=RKIND), dimension(:,:,:), allocatable :: oceanFieldOld err = 0 @@ -392,6 +418,7 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali block => domain % blocklist call mpas_pool_get_config(liConfigs, 'config_ocean_data_extrap_ncells_extra', nCellsExtra) call mpas_pool_get_config(liConfigs,'config_weight_value_cell', weightCell) + call mpas_pool_get_config(liConfigs, 'config_recalculate_thermal_forcing', config_recalculate_thermal_forcing) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'extrapOceanData', extrapOceanDataPool) @@ -404,13 +431,13 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(extrapOceanDataPool, 'TFoceanOld', TFoceanOld) !TFoceanOld(:,:) = 1.0 ! for debugging, set TF to ones to make it easy to verify the horizonal/vertical averaging ! perform horizontal extrapolation until the validOceanMask is unchanged allocate(validOceanMaskOld(nISMIP6OceanLayers,nCells+1)) + allocate(oceanFieldOld(nISMIP6OceanLayers,nCells+1,2)) validOceanMaskOld(:,:) = validOceanMask(:,:) - TFoceanOld(:,:) = TFocean(:,:) + oceanFieldOld(:,:,:) = oceanField(:,:,:) ! initialize the local loop and count for validOceanMask localLoopCount = 0 @@ -422,46 +449,52 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali newMaskCountLocalAccum = 0 do iCell = 1, nCellsSolve do iLayer = 1, nISMIP6OceanLayers - if ( (availOceanMask(iLayer,iCell) == 1) .and. (validOceanMaskOrig(iLayer,iCell) == 0) ) then - TFsum = 0.0 - areaSum = 0.0 - newValidCount = 0 - nValidNeighb = 0 - do iNeighbor = 1, nEdgesOnCell(iCell) - jCell = cellsOnCell(iNeighbor, iCell) - 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) - err = ior(err,1) - else - TFsum = TFsum + (TFoceanOld(iLayer,jCell) * areaCell(jCell)) + do iField = 1, 2 + if ( (availOceanMask(iLayer,iCell) == 1) .and. (validOceanMaskOrig(iLayer,iCell) == 0) ) then + fieldSum = 0.0 + areaSum = 0.0 + newValidCount = 0 + nValidNeighb = 0 + do iNeighbor = 1, nEdgesOnCell(iCell) + jCell = cellsOnCell(iNeighbor, iCell) + if ( validOceanMaskOld(iLayer,jCell) == 1 ) then + if ( oceanFieldOld(iLayer,jCell,iField) > 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) + err = ior(err,1) + else + fieldSum = fieldSum + (oceanFieldOld(iLayer,jCell,iField) * areaCell(jCell)) + endif + areaSum = areaSum + areaCell(jCell) + nValidNeighb = nValidNeighb + 1 endif - areaSum = areaSum + areaCell(jCell) - nValidNeighb = nValidNeighb + 1 + enddo + if ( validOceanMaskOld(iLayer,iCell) == 0 .and. nValidNeighb > 0 ) then + ! if current cell is not valid, set its weight to zero + weightCellLocal = 0.0_RKIND + validOceanMask(iLayer,iCell) = 1 + newValidCount = 1 + else + weightCellLocal = weightCell endif - enddo - if ( validOceanMaskOld(iLayer,iCell) == 0 .and. nValidNeighb > 0 ) then - ! if current cell is not valid, set its weight to zero - weightCellLocal = 0.0_RKIND - validOceanMask(iLayer,iCell) = 1 - newValidCount = 1 - else - weightCellLocal = weightCell + ! perform area-weighted averaging of ocean fields + if ( nValidNeighb == 0 ) then + oceanField(iLayer,iCell,iField) = oceanFieldOld(iLayer,iCell,iField) + else + oceanField(iLayer,iCell,iField) = ( weightCellLocal * oceanFieldOld(iLayer,iCell,iField) * areaCell(iCell) + & + & ((1.0_RKIND - weightCellLocal) * (fieldSum / nValidNeighb)) ) / & + & ( weightCellLocal * areaCell(iCell) + & + & (1.0_RKIND - weightCellLocal) * (areaSum / nValidNeighb) ) + endif + ! Accumulate cells added locally until we do the next global reduce + newMaskCountLocalAccum = newMaskCountLocalAccum + newValidCount endif - ! perform area-weighted averaging of the thermal forcing field - if ( nValidNeighb == 0 ) then - TFocean(iLayer,iCell) = TFoceanOld(iLayer,iCell) - else - TFocean(iLayer,iCell) = ( weightCellLocal * TFoceanOld(iLayer,iCell) * areaCell(iCell) + & - & ((1.0_RKIND - weightCellLocal) * (TFsum / nValidNeighb)) ) / & - & ( weightCellLocal * areaCell(iCell) + & - & (1.0_RKIND - weightCellLocal) * (areaSum / nValidNeighb) ) + if (.not. config_recalculate_thermal_forcing) then + ! Don't repeat loop if only extrapolating thermal forcing + exit endif - ! Accumulate cells added locally until we do the next global reduce - newMaskCountLocalAccum = newMaskCountLocalAccum + newValidCount - endif + enddo enddo enddo @@ -482,6 +515,7 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali call mpas_log_write('Horizontal extrapolation done after $i iterations. Added total of $i cells across all processors', & intArgs=(/localLoopCount, newMaskCountTotal/)) deallocate(validOceanMaskOld) + deallocate(oceanFieldOld) end subroutine horizontal_extrapolation !----------------------------------------------------------------------- @@ -492,18 +526,19 @@ end subroutine horizontal_extrapolation ! ! routine vertical_extrapolation ! !> \brief Extrapolate validOceanMask vertically -!> \author Holly Kyeore Han -!> \date November 2023 +!> \author Holly Kyeore Han, modified by Alex Hager +!> \date November 2023 (modified September 2024) !> \details !> This routine extrapolates the horizontally extrapolated !> validOceanMask through the vertical layers of the ocean. !> The vertical extrapolation is completed once local new mask count !> stops updating. The output of the routine is an updated -!> validOceanMask, thermal forcing fields and newMaskCountLocal. +!> validOceanMask, ocean forcing fields (temperature and salinity, or +!> thermal forcing) and newMaskCountLocal. !----------------------------------------------------------------------- - subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMaskCountGlobal, TFocean, err) + subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMaskCountGlobal, oceanField, err) use li_constants, only: oceanFreezingTempDepthDependence @@ -518,7 +553,7 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas type (domain_type), intent(inout) :: domain !< Input/Output: domain object integer, dimension(:,:), pointer, intent(inout) :: validOceanMask integer, intent(inout) :: err !< Output: error flag - real (kind=RKIND), dimension(:,:), pointer, intent(inout) :: TFocean + real (kind=RKIND), dimension(:,:), pointer, intent(inout) :: oceanField !----------------------------------------------------------------- ! output variables @@ -530,19 +565,21 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas !----------------------------------------------------------------- type (block_type), pointer :: block type (mpas_pool_type), pointer :: scratchPool, geometryPool, meshPool, extrapOceanDataPool - real (kind=RKIND) :: layerTop, TFsum, areaSum + real (kind=RKIND) :: layerTop, fieldSum, areaSum real (kind=RKIND), dimension(:), pointer :: ismip6shelfMelt_zOcean real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography, areaCell integer, pointer :: nCells, nCellsSolve, nISMIP6OceanLayers integer, dimension(:), pointer :: cellMask, nEdgesOnCell integer, dimension(:,:), pointer :: cellsOnCell - integer :: iCell, jCell, iLayer, iNeighbor, iter + integer :: iCell, jCell, iLayer, iNeighbor, iter, iField integer :: localLoopCount, newMaskCountLocalAccum + logical, pointer :: config_recalculate_thermal_forcing err = 0 ! initialize the ocean data and mask fields block => domain % blocklist + call mpas_pool_get_config(liConfigs, 'config_recalculate_thermal_forcing', config_recalculate_thermal_forcing) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'extrapOceanData', extrapOceanDataPool) @@ -562,22 +599,34 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas newMaskCountLocalAccum = 0 do iCell = 1, nCellsSolve do iLayer = 2, nISMIP6OceanLayers - if ( (availOceanMask(iLayer,iCell) == 1) .and. (validOceanMask(iLayer,iCell) == 0) ) then - 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) - err = ior(err,1) - else - TFocean(iLayer,iCell) = TFocean(iLayer-1,iCell) - & - (ismip6shelfMelt_zOcean(iLayer-1) - ismip6shelfMelt_zOcean(iLayer)) * & - oceanFreezingTempDepthDependence - endif - validOceanMask(iLayer,iCell) = 1 - newMaskCountLocalAccum = newMaskCountLocalAccum + 1 - endif - endif + do iField = 1, 2 + if ( (availOceanMask(iLayer,iCell) == 1) .and. (validOceanMask(iLayer,iCell) == 0) ) then + if ( validOceanMask(iLayer-1,iCell) == 1 ) then + if (oceanField(iLayer-1,iCell,iField) > 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) + err = ior(err,1) + else + ! << TO DO >>: This will have to be done differently for temp/salt + if (config_recalculate_thermal_forcing) then + !if extrapolating temp/salt, assign value equal to overlying cell + oceanField(iLayer,iCell,iField) = oceanField(iLayer-1,iCell,iField) + else + !if extrapolating thermal forcing, need to account for depth-dependent difference in freezing temp + oceanField(iLayer,iCell,iField) = oceanField(iLayer-1,iCell,iField) - & + (ismip6shelfMelt_zOcean(iLayer-1) - ismip6shelfMelt_zOcean(iLayer)) * & + oceanFreezingTempDepthDependence + endif + validOceanMask(iLayer,iCell) = 1 + newMaskCountLocalAccum = newMaskCountLocalAccum + 1 + endif + endif + + if (.not. config_recalculate_thermal_forcing) then + exit + endif + enddo enddo enddo From 735add86d4d279a3e08da445887339c85a4d4ea1 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Mon, 16 Sep 2024 09:32:32 -0600 Subject: [PATCH 7/8] Add thermal forcing calculation routine Add recalculate_thermal_forcing subroutine to ocean extrapolation. By setting config_recalculate_thermal_forcing to true, the user can now choose to calculate thermal forcing at the glacier face using any of the methods outlined in Hager et al. (2024). If config_recalculate_thermal_forcing is false then MALI will default to using the extrapolated version of ismip6shelfMelt_3dThermalForcing from the input file. --- .../mpas-albany-landice/src/Registry.xml | 18 +- .../src/mode_forward/mpas_li_ocean_extrap.F | 301 ++++++++++++++++-- 2 files changed, 293 insertions(+), 26 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index ac90339fdad2..6162f4249c9f 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -457,11 +457,19 @@ description="If true, extrapolate ocean data (temperature and salinity, or thermal forcing) from external source into underneath the ice draft." possible_values=".true. or .false." /> - - + + + >: If using "retreat" athermal forcing paramterization, then need to translate T/S instead of extrap routine + ! call the extrapolation scheme call mpas_log_write('ocean data will be extrapolated into the MALI ice draft') block => domain % blocklist @@ -122,7 +134,8 @@ subroutine li_ocean_extrap_solve(domain, err) ! initialize the ocean data and mask fields call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) call mpas_pool_get_config(liConfigs, 'config_ocean_data_extrap_ncells_extra', nCellsExtra) - call mpas_pool_get_config(liConfigs, 'config_recalculate_thermal_forcing', config_recalcuate_thermal_forcing) + call mpas_pool_get_config(liConfigs, 'config_recalculate_thermal_forcing', config_recalculate_thermal_forcing) + call mpas_pool_get_config(liConfigs, 'config_TF_param_type', config_TF_param_type) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'extrapOceanData', extrapOceanDataPool) @@ -145,7 +158,7 @@ subroutine li_ocean_extrap_solve(domain, err) call mpas_pool_get_array(extrapOceanDataPool, 'seedOceanMaskHorizInit', seedOceanMaskHorizInit) allocate(oceanField(nISMIP6OceanLayers,nCells+1,2)) - + ! create a 2D mask based on the open ocean + floating ice + grounded ice mask to be used for defining the location of the n-extra cells into the grounded ice allocate(seedOceanMaskHorizOld(nCells+1)) seedOceanMaskHorizOld(:) = 0 @@ -167,7 +180,7 @@ subroutine li_ocean_extrap_solve(domain, err) do iter = 1, nCellsExtra do iCell = 1, nCellsSolve if ( growOceanMaskHoriz(iCell) == 1 .and. seedOceanMaskHorizOld(iCell) == 0 ) then - do iNeighbor /TFoce= 1, nEdgesOnCell(iCell) + do iNeighbor = 1, nEdgesOnCell(iCell) jCell = cellsOnCell(iNeighbor, iCell) if ( seedOceanMaskHorizOld(jCell) == 1 ) then seedOceanMaskHoriz(iCell) = 1 @@ -249,10 +262,18 @@ subroutine li_ocean_extrap_solve(domain, err) if (bedTopography(iCell) < layerTop) then availOceanMask(iLayer,iCell) = 1 else - ! keep the first layer below the seafloor in the region to be filled - ! this ensures linear interpolation from above and below the seafloor is possible - availOceanMask(iLayer,iCell) = 1 - exit ! stop looping over levels after we've included the first level below the seafloor + ! If recalculating thermal forcing with ISMIP6retreat or AMretreat methods, then ignore + ! topography and allow for translation of water properties directly to glacier face. Otherwise, + ! cut off availOceanMask at first cell below seaflooor; this cell is necessary to keep in + ! order to ensure interpolation between cells directly above and below sea level is possible + + if ((config_recalculate_thermal_forcing) .and. & + ( (config_TF_param_type == 'ISMIP6retreat') .or. (config_TF_param_type == 'AMretreat'))) then + availOceanMask(iLayer,iCell) = 1 + else + availOceanMask(iLayer,iCell) = 1 + exit + endif endif endif enddo @@ -267,16 +288,15 @@ subroutine li_ocean_extrap_solve(domain, err) ! initialize the extrapolated fields if (config_recalculate_thermal_forcing) then - real (kind=RKIND), dimension(:,:), pointer :: MPAS_3dOceanTemperature, MPAS_3dOceanSalinity - call mpas_pool_get_array(extrapOceanDataPool, 'MPAS_3dOceanTemperature', oceanTemp) - call mpas_pool_get_array(extrapOceanDataPool, 'MPAS_3dOceanSalinity', oceanSalt) + call mpas_pool_get_array(extrapOceanDataPool, 'MPAS_3dOceanTemperature', MPAS_3dOceanTemperature) + call mpas_pool_get_array(extrapOceanDataPool, 'MPAS_3dOceanSalinity', MPAS_3dOceanSalinity) !combine fields into one variable to expedite extrapolation oceanField(:,:,1) = MPAS_3dOceanTemperature(:,:) * validOceanMask(:,:) oceanField(:,:,2) = MPAS_3dOceanSalinity(:,:) * validOceanMask(:,:) else - oceanField(:,:,1) = ismip6shelfMelt_3dThermalForcing * validOceanMask(:,:)i + oceanField(:,:,1) = ismip6shelfMelt_3dThermalForcing * validOceanMask(:,:) oceanField(:,:,2) = invalid_value endif @@ -385,7 +405,7 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali type (domain_type), intent(inout) :: domain !< Input/Output: domain object integer, dimension(:,:), pointer, intent(inout) :: validOceanMask integer, intent(inout) :: err !< Output: error flag - real (kind=RKIND), dimension(:,:), pointer, intent(inout) :: oceanField + real (kind=RKIND), dimension(:,:,:), allocatable, intent(inout) :: oceanField !----------------------------------------------------------------- ! output variables @@ -400,7 +420,6 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali real (kind=RKIND), pointer :: weightCell integer, dimension(:,:), allocatable :: validOceanMaskOld real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing - real (kind=RKIND), dimension(:,:,:), pointer :: oceanFieldOld real (kind=RKIND), dimension(:), pointer :: thickness, bedTopography, areaCell integer, pointer :: nCells, nCellsSolve, nISMIP6OceanLayers, nCellsExtra integer, dimension(:), pointer :: cellMask, nEdgesOnCell @@ -432,7 +451,6 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - !TFoceanOld(:,:) = 1.0 ! for debugging, set TF to ones to make it easy to verify the horizonal/vertical averaging ! perform horizontal extrapolation until the validOceanMask is unchanged allocate(validOceanMaskOld(nISMIP6OceanLayers,nCells+1)) allocate(oceanFieldOld(nISMIP6OceanLayers,nCells+1,2)) @@ -498,14 +516,13 @@ subroutine horizontal_extrapolation(domain, availOceanMask, validOceanMask, vali enddo enddo - ! update halo for validOceanMask and ocean data + ! update halo for validOceanMask call mpas_timer_start("halo updates") call mpas_dmpar_field_halo_exch(domain, 'validOceanMask') - call mpas_dmpar_field_halo_exch(domain, 'TFocean') call mpas_timer_stop("halo updates") validOceanMaskOld(:,:) = validOceanMask(:,:) - TFoceanOld(:,:) = TFocean(:,:) + oceanFieldOld(:,:,:) = oceanField(:,:,:) ! update count of cells added to mask globally call mpas_dmpar_sum_int(domain % dminfo, newMaskCountLocalAccum, newMaskCountGlobal) @@ -553,7 +570,7 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas type (domain_type), intent(inout) :: domain !< Input/Output: domain object integer, dimension(:,:), pointer, intent(inout) :: validOceanMask integer, intent(inout) :: err !< Output: error flag - real (kind=RKIND), dimension(:,:), pointer, intent(inout) :: oceanField + real (kind=RKIND), dimension(:,:,:), allocatable, intent(inout) :: oceanField !----------------------------------------------------------------- ! output variables @@ -608,7 +625,6 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas MPAS_LOG_ERR) err = ior(err,1) else - ! << TO DO >>: This will have to be done differently for temp/salt if (config_recalculate_thermal_forcing) then !if extrapolating temp/salt, assign value equal to overlying cell oceanField(iLayer,iCell,iField) = oceanField(iLayer-1,iCell,iField) @@ -617,6 +633,7 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas oceanField(iLayer,iCell,iField) = oceanField(iLayer-1,iCell,iField) - & (ismip6shelfMelt_zOcean(iLayer-1) - ismip6shelfMelt_zOcean(iLayer)) * & oceanFreezingTempDepthDependence + endif endif validOceanMask(iLayer,iCell) = 1 newMaskCountLocalAccum = newMaskCountLocalAccum + 1 @@ -626,6 +643,7 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas if (.not. config_recalculate_thermal_forcing) then exit endif + enddo enddo enddo @@ -633,7 +651,6 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas ! update halo for validOceanMask and ocean data call mpas_timer_start("halo updates") call mpas_dmpar_field_halo_exch(domain, 'validOceanMask') - call mpas_dmpar_field_halo_exch(domain, 'TFocean') call mpas_timer_stop("halo updates") ! update count of cells added to mask globally call mpas_dmpar_sum_int(domain % dminfo, newMaskCountLocalAccum, newMaskCountGlobal) @@ -641,6 +658,248 @@ subroutine vertical_extrapolation(domain, availOceanMask, validOceanMask, newMas end subroutine vertical_extrapolation !----------------------------------------------------------------------- + + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ! routine recalculate_thermal_forcing +! +!> \brief recalculate thermal forcing from extrapolated temp/salt +!> \author Alex Hager +!> \date September 2024 +!> \details +!> Use extrapolated temperature and salt to recalculate thermal forcing +!> throughout the extrapolated region. User specifies which thermal forcing +!> parameterization from Hager et al. (2023) to use with +!> config_TF_param_type. Output thermal forcing fields are 2D. +!----------------------------------------------------------------------- + + subroutine recalculate_thermal_forcing(domain, err) + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + type (domain_type), intent(inout) :: domain !< Input/Output: domain object + integer, intent(inout) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + type (block_type), pointer :: block + type (mpas_pool_type), pointer :: geometryPool, meshPool, extrapOceanDataPool + character (len=StrKIND), pointer :: config_TF_param_type + real (kind=RKIND), pointer :: config_TF_iceberg_depth + real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing + real (kind=RKIND), dimension(:), pointer :: ismip6_2dThermalForcing + real (kind=RKIND), dimension(:,:), pointer :: MPAS_3dOceanTemperature + real (kind=RKIND), dimension(:,:), pointer :: MPAS_3dOceanSalinity + real (kind=RKIND), dimension(:), pointer :: ismip6shelfMelt_zOcean + real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_zBndsOcean + real (kind=RKIND), dimension(:), pointer :: icebergFjordMask + real (kind=RKIND), dimension(:), pointer :: bedTopography + real (kind=RKIND), pointer :: dzAccumulated, dz, total, gadeSlope + integer, pointer :: nISMIP6OceanLayers, nCells, iCell, iLayer, counter + integer, dimension(:), allocatable :: icebergCellVertID + integer, dimension(:), allocatable :: oceanCellVertID + real (kind=RKIND), dimension(:), pointer :: freezingTemperature + allocate(icebergCellVertID(nISMIP6OceanLayers)) + allocate(oceanCellVertID(nISMIP6OceanLayers)) + allocate(freezingTemperature(nISMIP6OceanLayers)) + + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) + call mpas_pool_get_config(liConfigs, 'config_TF_param_type', config_TF_param_type) + call mpas_pool_get_config(liConfigs, 'config_TF_iceberg_depth', config_TF_iceberg_depth) + call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing) + call mpas_pool_get_array(geometryPool, 'ismip6_2dThermalForcing', ismip6_2dThermalForcing) + call mpas_pool_get_array(geometryPool, 'MPAS_3dOceanTemperature', MPAS_3dOceanTemperature) + call mpas_pool_get_array(geometryPool, 'MPAS_3dOceanSalinity', MPAS_3dOceanSalinity) + call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zOcean', ismip6shelfMelt_zOcean) + call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zBndsOcean', ismip6shelfMelt_zBndsOcean) + call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) + call mpas_pool_get_array(geometryPool, 'icebergFjordMask', icebergFjordMask) + call mpas_pool_get_array(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers) + call mpas_pool_get_array(meshPool, 'nCells', nCells) + + if (config_TF_param_type == 'AMretreat') then + !depth-bin-weighted average of profiles + do iCell = 1, nCells + ! From Jenkins 2011 and Slater 2020 + ismip6shelfMelt_3dThermalForcing(:,iCell) = MPAS_3dOceanTemperature(:,iCell) - (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) + + ismip6_2dThermalForcing = 0.0_RKIND + dzAccumulated = 0.0_RKIND + do iLayer = 1, nISMIP6OceanLayers + if (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer)) then + dz = abs(ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer)) + ismip6_2dThermalForcing(iCell) = (ismip6_2dThermalForcing(iCell) * dzAccumulated + & + ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) / (dzAccumulated + dz) + dzAccumulated = dzAccumulated + dz + endif + enddo + enddo + + elseif (config_TF_param_type == 'ISMIP6retreat') then + ! simple average between 200-500m + do iCell = 1, nCells + total = 0.0_RKIND + counter = 0 + do iLayer = 1, nISMIP6OceanLayers + if ( (ismip6shelfMelt_zOcean(iLayer) >= 200.0_RKIND) .and. ( ismip6shelfMelt_zOcean(iLayer) <= 500.0_RKIND) & + .and. (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer))) then + total = ismip6shelfMelt_3dThermalForcing(iLayer,iCell) + total + counter = counter + 1 + endif + enddo + ismip6_2dThermalForcing(iCell) = total/counter + enddo + + elseif (config_TF_param_type == 'AMmelt') then + ! depth-bin-weighted average of profiles + do iCell = 1, nCells + ! << NOTE>> : Using TF equations from Jenkins 2011 for now but Slater 2020 uses gsw toolbox version. May need to change this eventually + ismip6shelfMelt_3dThermalForcing(:,iCell) = MPAS_3dOceanTemperature(:,iCell) - (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) + + ismip6_2dThermalForcing = 0.0_RKIND + dzAccumulated = 0.0_RKIND + do iLayer = 1, nISMIP6OceanLayers + if (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer)) then + dz = ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer) + ismip6_2dThermalForcing(iCell) = (ismip6_2dThermalForcing(iCell) * dzAccumulated & + + ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) / (dzAccumulated + dz) + dzAccumulated = dzAccumulated + dz + endif + enddo + enddo + + elseif (config_TF_param_type == 'ISMIP6melt') then + ! define TF at grounding line + do iCell = 1, nCells + oceanCellVertID(:) = 0 + do iLayer = 1, nISMIP6OceanLayers + if (ismip6shelfMelt_zOcean(iLayer) >= bedTopography(iCell)) then + oceanCellVertID(iLayer) = iLayer + endif + enddo + ismip6_2dThermalForcing(iCell) = ismip6shelfMelt_3dThermalForcing(maxval(oceanCellVertID),iCell) + enddo + + elseif (config_TF_param_type == 'AMberg') then + ! Iceberg Parameteriation from Hager 2023 + + do iCell = 1, nCells + ! find cell just above config_TF_iceberg_depth + icebergCellVertID(:) = 0 + do iLayer = 1, nISMIP6OceanLayers + if ( (ismip6shelfMelt_zOcean(iLayer) <= config_TF_iceberg_depth) & + .and. (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer)) )then + icebergCellVertID(iLayer) = iLayer + endif + enddo + + ! calculate freezing temperature. + ! << NOTE >>: Again we are using the simplified equation from Jenkins 2011 to do this. May eventually want to + ! switch to method in gsw toolbox to do this to be consistent with Hager 2023 + freezingTemperature = (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) + + ! Calculate Gade slope (slope of submarine melt mixing line) using cell just above iceberg depth as reference point. Assumes ice at pressure melting temperature + gadeSlope = (1/MPAS_3dOceanSalinity(maxval(icebergCellVertID),iCell)) * (L/cp - (freezingTemperature(maxval(icebergCellVertID)) & + - MPAS_3dOceanTemperature(maxval(icebergCellVertID),iCell))) + + ! Now that we have the slope of the melt line, alter ocean temperatures along submarine melt line + do iLayer = 1, nISMIP6OceanLayers + if (icebergCellVertID(iLayer) /= 0) then + MPAS_3dOceanTemperature(iLayer, iCell) = gadeSlope * (MPAS_3dOceanSalinity(iLayer, iCell) - MPAS_3dOceanSalinity(maxval(icebergCellVertID),iCell)) & + + MPAS_3dOceanTemperature(iLayer, iCell) + + ! Ensure no artificially adjusted temperatures fall below the freezing temperature + if (MPAS_3dOceanTemperature(iLayer, iCell) < freezingTemperature(iLayer)) then + MPAS_3dOceanTemperature(iLayer, iCell) = freezingTemperature(iLayer) + endif + + ismip6shelfMelt_3dThermalForcing(iLayer, iCell) = MPAS_3dOceanTemperature(iLayer, iCell) - freezingTemperature(iLayer) + endif + enddo + + ismip6_2dThermalForcing = 0.0_RKIND + dzAccumulated = 0.0_RKIND + do iLayer = 1, nISMIP6OceanLayers + dz = ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer) + ismip6_2dThermalForcing(iCell) = (ismip6_2dThermalForcing(iCell) * dzAccumulated & + + ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) / (dzAccumulated + dz) + dzAccumulated = dzAccumulated + dz + enddo + enddo + + elseif (config_TF_param_type == 'AMfit') then + + do iCell = 1, nCells + if (icebergFjordMask(iCell) == 1) then + + ! find cell just above config_TF_iceberg_depth + icebergCellVertID(:) = 0 + do iLayer = 1, nISMIP6OceanLayers + if (ismip6shelfMelt_zOcean(iLayer) <= config_TF_iceberg_depth) then + icebergCellVertID(iLayer) = iLayer + endif + enddo + + ! calculate freezing temperature. + ! << NOTE >>: Again we are using the simplified equation from Jenkins 2011 to do this. May eventually want to + ! switch to method in gsw toolbox to do this to be consistent with Hager 2023 + freezingTemperature = (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) + + ! Calculate Gade slope (slope of submarine melt mixing line) using cell just above iceberg depth as reference point. Assumes ice at pressure melting temperature + gadeSlope = (1/MPAS_3dOceanSalinity(maxval(icebergCellVertID),iCell)) * (L/cp - (freezingTemperature(maxval(icebergCellVertID)) & + - MPAS_3dOceanTemperature(maxval(icebergCellVertID),iCell))) + + ! Now that we have the slope of the melt line, alter ocean temperatures along submarine melt line + do iLayer = 1, nISMIP6OceanLayers + if (icebergCellVertID(iLayer) /= 0) then + MPAS_3dOceanTemperature(iLayer, iCell) = gadeSlope * (MPAS_3dOceanSalinity(iLayer, iCell) - MPAS_3dOceanSalinity(maxval(icebergCellVertID),iCell)) & + + MPAS_3dOceanTemperature(iLayer, iCell) + + ! Ensure no artificially adjusted temperatures fall below the freezing temperature + if (MPAS_3dOceanTemperature(iLayer, iCell) < freezingTemperature(iLayer)) then + MPAS_3dOceanTemperature(iLayer, iCell) = freezingTemperature(iLayer) + endif + + ismip6shelfMelt_3dThermalForcing(iLayer, iCell) = MPAS_3dOceanTemperature(iLayer, iCell) - freezingTemperature(iLayer) + endif + enddo + endif + + ! << NOTE>> : Using TF equations from Jenkins 2011 for now but Slater 2020 uses gsw toolbox version. May need to change this eventually + ismip6shelfMelt_3dThermalForcing(:,iCell) = MPAS_3dOceanTemperature(:,iCell) - (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) + + ismip6_2dThermalForcing = 0.0_RKIND + dzAccumulated = 0.0_RKIND + do iLayer = 1, nISMIP6OceanLayers + dz = ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer) + ismip6_2dThermalForcing(iCell) = (ismip6_2dThermalForcing(iCell) * dzAccumulated & + + ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) / (dzAccumulated + dz) + dzAccumulated = dzAccumulated + dz + enddo + enddo + endif + + deallocate(icebergCellVertID) + deallocate(freezingTemperature) + + block => block % next + enddo + end subroutine recalculate_thermal_forcing end module li_ocean_extrap + !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| From 8edf1784dfb22547ae49673b6adcf69a919d2b12 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Mon, 4 Nov 2024 13:01:14 -0700 Subject: [PATCH 8/8] 2d from 3d fields moved to mpas_li_iceshelf_melt.F Moves consolidation of 3d thermal forcing profiles into 2d scalars, previously done in mpas_li_ocean_extrap.F, to mpas_li_iceshelf_melt.F. This is done to be consistent with MALI configurations that do not use the ocean extrapolation framework. --- .../src/mode_forward/mpas_li_iceshelf_melt.F | 68 +++++++-- .../src/mode_forward/mpas_li_ocean_extrap.F | 137 ++---------------- 2 files changed, 63 insertions(+), 142 deletions(-) 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..2e0824a0989a 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 @@ -1466,27 +1466,63 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zOcean', zOcean) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_deltaT', ismip6shelfMelt_deltaT) - ! Calculate thermal forcing based on bedTopography - TFocean(:) = 0.0_RKIND - do iCell = 1, nCells - if (bedTopography(iCell) < 0.0_RKIND) then - kk = 1 - ! Find deepest zOcean above ocean floor - do while ( (zOcean(kk) >= bedTopography(iCell)) .and. (kk < nISMIP6OceanLayers) ) - kk = kk + 1 + !Consolidate 3d profile to 2d depending on parameterization type + if ( (config_recalculate_thermal_forcing) .and. ( (config_TF_param_type == 'AMretreat') & + .or. (config_TF_param_type == 'AMmelt') .or. (config_TF_param_type == 'AMberg') & + .or. (config_TF_param_type == 'AMfit') ) ) then + ! Depth weighted average + TFocean(:) = 0.0_RKIND + do iCell = 1, nCells + dzAccumulated = 0.0_RKIND + do iLayer = 1, nISMIP6OceanLayers + if (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer)) then + dz = abs(ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer)) + TFocean(iCell) = (TFocean(iCell) * dzAccumulated + ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) & + / (dzAccumulated + dz) + ismip6shelfMelt_deltaT(iCell) + dzAccumulated = dzAccumulated + dz + endif enddo - ! If bed is shallower than first layer, use TF from the first layer. - ! If bed is deeper than the bottomg ocean layer, use TF from the bottom layer. - if ( (kk == 1) .or. ( (kk == nISMIP6OceanLayers) .and. (zOcean(kk) >= bedTopography(iCell)) ) ) then - TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk, iCell) + ismip6shelfMelt_deltaT(iCell) - ! For all other bed depths, interpolate linearly between layers above and below bed depth. - else - TFocean(iCell) = ( (zOcean(kk-1) - bedTopography(iCell)) * ismip6shelfMelt_3dThermalForcing(kk, iCell) + & + enddo + + elseif ( (config_recalculate_thermal_forcing) .and. (config_TF_param_type == 'ISMIP6retreat') ) then + ! simple average between 200-500m + TFocean(:) = 0.0_RKIND + do iCell = 1, nCells + total = 0.0_RKIND + counter = 0 + do iLayer = 1, nISMIP6OceanLayers + if ( (ismip6shelfMelt_zOcean(iLayer) >= 200.0_RKIND) .and. ( ismip6shelfMelt_zOcean(iLayer) <= 500.0_RKIND) & + .and. (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer))) then + total = ismip6shelfMelt_3dThermalForcing(iLayer,iCell) + total + counter = counter + 1 + endif + enddo + TFocean(iCell) = total/counter + ismip6shelfMelt_deltaT(iCell) + enddo + + elseif ( (config_TF_param_type == 'ISMIP6melt') .or. (.not. config_recalculate_thermal_forcing) ) then + ! Calculate thermal forcing based solely on bedTopography. Default ISMIP6 method and ISMIP6melt are equivalent + TFocean(:) = 0.0_RKIND + do iCell = 1, nCells + if (bedTopography(iCell) < 0.0_RKIND) then + kk = 1 + ! Find deepest zOcean above ocean floor + do while ( (zOcean(kk) >= bedTopography(iCell)) .and. (kk < nISMIP6OceanLayers) ) + kk = kk + 1 + enddo + ! If bed is shallower than first layer, use TF from the first layer. + ! If bed is deeper than the bottomg ocean layer, use TF from the bottom layer. + if ( (kk == 1) .or. ( (kk == nISMIP6OceanLayers) .and. (zOcean(kk) >= bedTopography(iCell)) ) ) then + TFocean(iCell) = ismip6shelfMelt_3dThermalForcing(kk, iCell) + ismip6shelfMelt_deltaT(iCell) + ! For all other bed depths, interpolate linearly between layers above and below bed depth. + else + TFocean(iCell) = ( (zOcean(kk-1) - bedTopography(iCell)) * ismip6shelfMelt_3dThermalForcing(kk, iCell) + & (bedTopography(iCell) - zOcean(kk)) * ismip6shelfMelt_3dThermalForcing(kk-1, iCell) ) / & (zOcean(kk-1) - zOcean(kk)) + ismip6shelfMelt_deltaT(iCell) + endif endif + enddo endif - end do endif faceMeltSpeed(:) = 0.0_RKIND 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 cd7dc6f91432..a1831ca9e869 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 @@ -731,73 +731,18 @@ subroutine recalculate_thermal_forcing(domain, err) call mpas_pool_get_array(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers) call mpas_pool_get_array(meshPool, 'nCells', nCells) - if (config_TF_param_type == 'AMretreat') then - !depth-bin-weighted average of profiles - do iCell = 1, nCells + do iCell = 1, nCells + !Main difference between these parameterization is how we collapse into 2d, which happens in mpas_li_ice_shelt_melt.F + !Difference between '...melt' and '...retreat' parameterizations happens in main body of mpas_li_ocean_extrap.F + if ( (config_TF_param_type == 'AMretreat') .or. (config_TF_param_type == 'ISMIP6retreat') & + (config_TF_param_type == 'AMmelt') .or. (config_TF_param_type == 'ISMIP6melt') & + ((config_TF_param_type == 'AMfit') .and. (icebergFjordMask(iCell))) ) then + ! From Jenkins 2011 and Slater 2020 ismip6shelfMelt_3dThermalForcing(:,iCell) = MPAS_3dOceanTemperature(:,iCell) - (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) - - ismip6_2dThermalForcing = 0.0_RKIND - dzAccumulated = 0.0_RKIND - do iLayer = 1, nISMIP6OceanLayers - if (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer)) then - dz = abs(ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer)) - ismip6_2dThermalForcing(iCell) = (ismip6_2dThermalForcing(iCell) * dzAccumulated + & - ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) / (dzAccumulated + dz) - dzAccumulated = dzAccumulated + dz - endif - enddo - enddo - - elseif (config_TF_param_type == 'ISMIP6retreat') then - ! simple average between 200-500m - do iCell = 1, nCells - total = 0.0_RKIND - counter = 0 - do iLayer = 1, nISMIP6OceanLayers - if ( (ismip6shelfMelt_zOcean(iLayer) >= 200.0_RKIND) .and. ( ismip6shelfMelt_zOcean(iLayer) <= 500.0_RKIND) & - .and. (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer))) then - total = ismip6shelfMelt_3dThermalForcing(iLayer,iCell) + total - counter = counter + 1 - endif - enddo - ismip6_2dThermalForcing(iCell) = total/counter - enddo - - elseif (config_TF_param_type == 'AMmelt') then - ! depth-bin-weighted average of profiles - do iCell = 1, nCells - ! << NOTE>> : Using TF equations from Jenkins 2011 for now but Slater 2020 uses gsw toolbox version. May need to change this eventually - ismip6shelfMelt_3dThermalForcing(:,iCell) = MPAS_3dOceanTemperature(:,iCell) - (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) - - ismip6_2dThermalForcing = 0.0_RKIND - dzAccumulated = 0.0_RKIND - do iLayer = 1, nISMIP6OceanLayers - if (bedTopography(iCell) >= ismip6shelfMelt_zOcean(iLayer)) then - dz = ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer) - ismip6_2dThermalForcing(iCell) = (ismip6_2dThermalForcing(iCell) * dzAccumulated & - + ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) / (dzAccumulated + dz) - dzAccumulated = dzAccumulated + dz - endif - enddo - enddo - elseif (config_TF_param_type == 'ISMIP6melt') then - ! define TF at grounding line - do iCell = 1, nCells - oceanCellVertID(:) = 0 - do iLayer = 1, nISMIP6OceanLayers - if (ismip6shelfMelt_zOcean(iLayer) >= bedTopography(iCell)) then - oceanCellVertID(iLayer) = iLayer - endif - enddo - ismip6_2dThermalForcing(iCell) = ismip6shelfMelt_3dThermalForcing(maxval(oceanCellVertID),iCell) - enddo - - elseif (config_TF_param_type == 'AMberg') then - ! Iceberg Parameteriation from Hager 2023 - - do iCell = 1, nCells + elseif ( (config_TF_param_type == 'AMberg') .or. ( (config_TF_param_type == 'AMfit') .and. (icebergFjordMask(iCell))) ) then + ! Alter 3d profile following the iceberg Parameterization from Hager 2023 for fjord where icebergs are prevalent ! find cell just above config_TF_iceberg_depth icebergCellVertID(:) = 0 do iLayer = 1, nISMIP6OceanLayers @@ -830,68 +775,8 @@ subroutine recalculate_thermal_forcing(domain, err) ismip6shelfMelt_3dThermalForcing(iLayer, iCell) = MPAS_3dOceanTemperature(iLayer, iCell) - freezingTemperature(iLayer) endif enddo - - ismip6_2dThermalForcing = 0.0_RKIND - dzAccumulated = 0.0_RKIND - do iLayer = 1, nISMIP6OceanLayers - dz = ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer) - ismip6_2dThermalForcing(iCell) = (ismip6_2dThermalForcing(iCell) * dzAccumulated & - + ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) / (dzAccumulated + dz) - dzAccumulated = dzAccumulated + dz - enddo - enddo - - elseif (config_TF_param_type == 'AMfit') then - - do iCell = 1, nCells - if (icebergFjordMask(iCell) == 1) then - - ! find cell just above config_TF_iceberg_depth - icebergCellVertID(:) = 0 - do iLayer = 1, nISMIP6OceanLayers - if (ismip6shelfMelt_zOcean(iLayer) <= config_TF_iceberg_depth) then - icebergCellVertID(iLayer) = iLayer - endif - enddo - - ! calculate freezing temperature. - ! << NOTE >>: Again we are using the simplified equation from Jenkins 2011 to do this. May eventually want to - ! switch to method in gsw toolbox to do this to be consistent with Hager 2023 - freezingTemperature = (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) - - ! Calculate Gade slope (slope of submarine melt mixing line) using cell just above iceberg depth as reference point. Assumes ice at pressure melting temperature - gadeSlope = (1/MPAS_3dOceanSalinity(maxval(icebergCellVertID),iCell)) * (L/cp - (freezingTemperature(maxval(icebergCellVertID)) & - - MPAS_3dOceanTemperature(maxval(icebergCellVertID),iCell))) - - ! Now that we have the slope of the melt line, alter ocean temperatures along submarine melt line - do iLayer = 1, nISMIP6OceanLayers - if (icebergCellVertID(iLayer) /= 0) then - MPAS_3dOceanTemperature(iLayer, iCell) = gadeSlope * (MPAS_3dOceanSalinity(iLayer, iCell) - MPAS_3dOceanSalinity(maxval(icebergCellVertID),iCell)) & - + MPAS_3dOceanTemperature(iLayer, iCell) - - ! Ensure no artificially adjusted temperatures fall below the freezing temperature - if (MPAS_3dOceanTemperature(iLayer, iCell) < freezingTemperature(iLayer)) then - MPAS_3dOceanTemperature(iLayer, iCell) = freezingTemperature(iLayer) - endif - - ismip6shelfMelt_3dThermalForcing(iLayer, iCell) = MPAS_3dOceanTemperature(iLayer, iCell) - freezingTemperature(iLayer) - endif - enddo - endif - - ! << NOTE>> : Using TF equations from Jenkins 2011 for now but Slater 2020 uses gsw toolbox version. May need to change this eventually - ismip6shelfMelt_3dThermalForcing(:,iCell) = MPAS_3dOceanTemperature(:,iCell) - (gamma1 * MPAS_3dOceanSalinity(:,iCell) + gamma2 + gamma3 * ismip6shelfMelt_zOcean) - - ismip6_2dThermalForcing = 0.0_RKIND - dzAccumulated = 0.0_RKIND - do iLayer = 1, nISMIP6OceanLayers - dz = ismip6shelfMelt_zBndsOcean(2,iLayer) - ismip6shelfMelt_zBndsOcean(1,iLayer) - ismip6_2dThermalForcing(iCell) = (ismip6_2dThermalForcing(iCell) * dzAccumulated & - + ismip6shelfMelt_3dThermalForcing(iLayer,iCell) * dz) / (dzAccumulated + dz) - dzAccumulated = dzAccumulated + dz - enddo - enddo - endif + endif + enddo deallocate(icebergCellVertID) deallocate(freezingTemperature)