diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 9ce37b56fdfc..6162f4249c9f 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -454,16 +454,24 @@ possible_values="any non-negative value" /> + + + - + - + + - - @@ -1689,6 +1697,10 @@ is the value of that variable from the *previous* time level! + + @@ -1707,11 +1719,11 @@ is the value of that variable from the *previous* time level! - - 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_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 16bc4eddd94b..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 @@ -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,7 +45,16 @@ module li_ocean_extrap !-------------------------------------------------------------------- ! Private module variables !-------------------------------------------------------------------- - + + real (kind=RKIND), parameter :: invalid_value = 1.0e36_RKIND + + ! Coefficients in Jenkins 2011 shelf melt parameterizations + real (kind=RKIND), parameter :: gamma1 = -5.7E-2 + real (kind=RKIND), parameter :: gamma2 = 8.8E-2 + real (kind=RKIND), parameter :: gamma3 = 7.61E-4 + + real (kind=RKIND), parameter :: cp = 3974 !specific heat (J/kg), using value from Cowton 2015 + real (kind=RKIND), parameter :: L = 3.34E5 !latent heat of melting !*********************************************************************** contains @@ -64,6 +73,7 @@ module li_ocean_extrap !----------------------------------------------------------------------- subroutine li_ocean_extrap_solve(domain, err) + use li_calving, only: li_flood_fill !----------------------------------------------------------------- ! input variables @@ -83,23 +93,31 @@ 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 + logical, pointer :: config_recalculate_thermal_forcing + character (len=StrKIND), pointer :: config_TF_param_type + 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 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 + real (kind=RKIND), dimension(:,:), pointer :: MPAS_3dOceanSalinity, MPAS_3dOceanTemperature integer, dimension(:), allocatable :: seedOceanMaskHorizOld integer, pointer :: nCells, nCellsSolve, nISMIP6OceanLayers, nCellsExtra integer, dimension(:), pointer :: cellMask, nEdgesOnCell 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 + real (kind=RKIND), dimension(:,:,:), allocatable :: oceanField + ! No init is needed. err = 0 err_tmp = 0 @@ -107,6 +125,8 @@ subroutine li_ocean_extrap_solve(domain, err) call mpas_pool_get_config(liConfigs, 'config_ocean_data_extrapolation', config_ocean_data_extrapolation) if ( config_ocean_data_extrapolation ) then + ! << TO DO >>: 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 @@ -114,6 +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_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) @@ -127,16 +149,16 @@ 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) 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)) seedOceanMaskHorizOld(:) = 0 @@ -178,37 +200,125 @@ 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 - endif - if ( (origOceanMaskHoriz(iCell) == 1) .and. (bedTopography(iCell) < layerTop) ) then - validOceanMask(iLayer,iCell) = 1 + if ( (seedOceanMaskHoriz(iCell) == 1) .and. (connectedMarineMask(iCell) == 1)) then + if (bedTopography(iCell) < layerTop) then + availOceanMask(iLayer,iCell) = 1 + else + ! 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 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. + validOceanMask(:,:) = orig3dOceanCavityMask(:,:) + ! save the initial validOceanMask validOceanMaskOrig(:,:) = validOceanMask(:,:) - ! initialize the TF field - TFocean(:,:) = ismip6shelfMelt_3dThermalForcing(:,:) * validOceanMask(:,:) + ! initialize the extrapolated fields + if (config_recalculate_thermal_forcing) then + + 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(:,:) + oceanField(:,:,2) = invalid_value + endif + ! initialize the invalid data locations with fill value do iCell = 1, nCellsSolve do iLayer = 1, nISMIP6OceanLayers - if ( availOceanMask(iLayer,iCell) == 0 ) then - TFocean(iLayer,iCell) = invalid_value_TF + 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 + oceanField(iLayer,iCell,:) = invalid_value 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 @@ -226,23 +336,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 !-------------------------------------------------------------------- @@ -263,8 +381,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 @@ -274,7 +392,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 @@ -287,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) :: TFocean + real (kind=RKIND), dimension(:,:,:), allocatable, intent(inout) :: oceanField !----------------------------------------------------------------- ! output variables @@ -298,18 +416,20 @@ 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 :: 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 @@ -317,6 +437,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) @@ -329,16 +450,16 @@ 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 + 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 ) @@ -346,87 +467,95 @@ 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 - ! 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) - 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) - - + deallocate(oceanFieldOld) 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 @@ -441,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) :: TFocean + real (kind=RKIND), dimension(:,:,:), allocatable, intent(inout) :: oceanField !----------------------------------------------------------------- ! output variables @@ -453,19 +582,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) @@ -485,29 +616,41 @@ 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 + 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 + endif + validOceanMask(iLayer,iCell) = 1 + newMaskCountLocalAccum = newMaskCountLocalAccum + 1 + endif + endif + + if (.not. config_recalculate_thermal_forcing) then + exit + endif + + enddo enddo enddo ! 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) @@ -515,6 +658,133 @@ 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) + + 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) + + 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 + 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 + endif + enddo + + deallocate(icebergCellVertID) + deallocate(freezingTemperature) + + block => block % next + enddo + end subroutine recalculate_thermal_forcing end module li_ocean_extrap + !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||