From 7ce3d1f33d137b79d4e522dca72eae0867aec8d2 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 29 Oct 2024 20:20:46 -0500 Subject: [PATCH] Fix total land-ice freshwater flux in data mode Previously, the total was only being computed when thermodynamics below ice shelves are actively computed, whereas we need to compute the total of the interface flux and the frazil flux when the interface flux comes from a data file as well. While we expect the frazil flux to be zero, these code modifications do not assume or require this to be true. --- .../shared/mpas_ocn_surface_land_ice_fluxes.F | 289 +++++++++--------- 1 file changed, 147 insertions(+), 142 deletions(-) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F b/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F index c249f202757d..94223bc4e95c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F @@ -492,168 +492,181 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & err = 0 - if (.not.landIceStandaloneOn) return + if (.not. (landIceStandaloneOn .or. landIceDataOn)) return call mpas_timer_start("land_ice_build_arrays") call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) - call mpas_pool_get_array(forcingPool, 'landIcePressure', landIcePressure) - - call mpas_pool_get_array(forcingPool, 'landIceFloatingFraction', landIceFloatingFraction) - call mpas_pool_get_array(forcingPool, 'landIceFloatingMask', landIceFloatingMask) - call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFlux) - call mpas_pool_get_array(forcingPool, 'landIceHeatFlux', landIceHeatFlux) - call mpas_pool_get_array(forcingPool, 'heatFluxToLandIce', heatFluxToLandIce) - call mpas_pool_get_array(forcingPool, 'frazilIceFreshwaterFlux', frazilIceFreshwaterFlux) call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFluxTotal', landIceFreshwaterFluxTotal) - call mpas_pool_get_array(forcingPool, 'landIceInterfaceTracers', landIceInterfaceTracers) - call mpas_pool_get_dimension(forcingPool, & - 'index_landIceInterfaceTemperature', & - indexITPtr) - call mpas_pool_get_dimension(forcingPool, & - 'index_landIceInterfaceSalinity', & - indexISPtr) - indexIT = indexITPtr - indexIS = indexISPtr - - if (useHollandJenkinsAdvDiff) then - call mpas_pool_get_array(forcingPool, 'landIceSurfaceTemperature', landIceSurfaceTemperature) - - allocate(freezeInterfaceSalinity(nCells), & - freezeInterfaceTemperature(nCells), & - freezeFreshwaterFlux(nCells), & - freezeHeatFlux(nCells), & - freezeIceHeatFlux(nCells)) - end if - nCells = nCellsArray( size(nCellsArray) ) - if (isomipOn) then !*** ISOMIP formulation + if (landIceStandaloneOn) then - !$omp parallel - !$omp do schedule(runtime) private(freshwaterFlux, heatFlux) - do iCell = 1, nCells - if (landIceFloatingMask(iCell) == 0) cycle + call mpas_pool_get_array(forcingPool, 'landIcePressure', landIcePressure) - ! linearized equaiton for the S and p dependent potential freezing temperature - landIceInterfaceTracers(indexIT,iCell) = ocn_freezing_temperature( & - salinity=landIceBoundaryLayerTracers(indexBLT,iCell), & - pressure=landIcePressure(iCell), & - inLandIceCavity=.true.) + call mpas_pool_get_array(forcingPool, 'landIceFloatingFraction', landIceFloatingFraction) + call mpas_pool_get_array(forcingPool, 'landIceFloatingMask', landIceFloatingMask) - ! using (3) and (4) from Hunter (2006) - ! or (7) from Jenkins et al. (2001) if gamma constant - ! and no heat flux into ice - ! freshwater flux = density * melt rate is in kg/m^2/s - freshwaterFlux = -rho_sw * ISOMIPgammaT * (cp_sw/latent_heat_fusion_mks) & - * (landIceInterfaceTracers(indexIT,iCell)-landIceBoundaryLayerTracers(indexBLT,iCell)) + call mpas_pool_get_array(forcingPool, 'landIceHeatFlux', landIceHeatFlux) + call mpas_pool_get_array(forcingPool, 'heatFluxToLandIce', heatFluxToLandIce) - landIceFreshwaterFlux(iCell) = landIceFloatingFraction(iCell)*freshwaterFlux - ! Using (13) from Jenkins et al. (2001) - ! heat flux is in W/s - heatFlux = cp_sw*(freshwaterFlux*landIceInterfaceTracers(indexIT,iCell) & - + rho_sw*ISOMIPgammaT & - * (landIceInterfaceTracers(indexIT,iCell)-landIceBoundaryLayerTracers(indexBLT,iCell))) - landIceHeatFlux(iCell) = landIceFloatingFraction(iCell)*heatFlux + call mpas_pool_get_array(forcingPool, 'landIceInterfaceTracers', landIceInterfaceTracers) + call mpas_pool_get_dimension(forcingPool, & + 'index_landIceInterfaceTemperature', & + indexITPtr) + call mpas_pool_get_dimension(forcingPool, & + 'index_landIceInterfaceSalinity', & + indexISPtr) + indexIT = indexITPtr + indexIS = indexISPtr - heatFluxToLandIce(iCell) = 0.0_RKIND + if (useHollandJenkinsAdvDiff) then + call mpas_pool_get_array(forcingPool, 'landIceSurfaceTemperature', landIceSurfaceTemperature) - end do - !$omp end do - !$omp end parallel - endif ! isomipOn + allocate(freezeInterfaceSalinity(nCells), & + freezeInterfaceTemperature(nCells), & + freezeFreshwaterFlux(nCells), & + freezeHeatFlux(nCells), & + freezeIceHeatFlux(nCells)) + end if - if (jenkinsOn .or. hollandJenkinsOn) then - if(useHollandJenkinsAdvDiff) then - ! melting solution - call compute_HJ99_melt_fluxes( & - landIceFloatingMask, & - landIceBoundaryLayerTracers(indexBLT,:), & - landIceBoundaryLayerTracers(indexBLS,:), & - landIceTracerTransferVelocities(indexHeatTrans,:), & - landIceTracerTransferVelocities(indexSaltTrans,:), & - landIceSurfaceTemperature, & - landIcePressure, & - landIceInterfaceTracers(indexIT,:), & - landIceInterfaceTracers(indexIS,:), & - landIceFreshwaterFlux, & - landIceHeatFlux, & - heatFluxToLandIce, & - nCells, & - err) - if(err .ne. 0) then - call mpas_log_write( & - 'compute_HJ99_melt_fluxes failed.', & - MPAS_LOG_CRIT) - end if + if (isomipOn) then !*** ISOMIP formulation - ! freezing solution - call compute_melt_fluxes( & - landIceFloatingMask, & - landIceBoundaryLayerTracers(indexBLT,:), & - landIceBoundaryLayerTracers(indexBLS,:), & - landIceTracerTransferVelocities(indexHeatTrans,:), & - landIceTracerTransferVelocities(indexSaltTrans,:), & - landIcePressure, & - freezeInterfaceTemperature, & - freezeInterfaceSalinity, & - freezeFreshwaterFlux, & - freezeHeatFlux, & - freezeIceHeatFlux, & - nCells, & - err) - if(err .ne. 0) then - call mpas_log_write( & - 'compute_melt_fluxes failed.', & - MPAS_LOG_CRIT) + !$omp parallel + !$omp do schedule(runtime) private(freshwaterFlux, heatFlux) + do iCell = 1, nCells + if (landIceFloatingMask(iCell) == 0) cycle + + ! linearized equaiton for the S and p dependent potential freezing temperature + landIceInterfaceTracers(indexIT,iCell) = ocn_freezing_temperature( & + salinity=landIceBoundaryLayerTracers(indexBLT,iCell), & + pressure=landIcePressure(iCell), & + inLandIceCavity=.true.) + + ! using (3) and (4) from Hunter (2006) + ! or (7) from Jenkins et al. (2001) if gamma constant + ! and no heat flux into ice + ! freshwater flux = density * melt rate is in kg/m^2/s + freshwaterFlux = -rho_sw * ISOMIPgammaT * (cp_sw/latent_heat_fusion_mks) & + * (landIceInterfaceTracers(indexIT,iCell)-landIceBoundaryLayerTracers(indexBLT,iCell)) + + landIceFreshwaterFlux(iCell) = landIceFloatingFraction(iCell)*freshwaterFlux + + ! Using (13) from Jenkins et al. (2001) + ! heat flux is in W/s + heatFlux = cp_sw*(freshwaterFlux*landIceInterfaceTracers(indexIT,iCell) & + + rho_sw*ISOMIPgammaT & + * (landIceInterfaceTracers(indexIT,iCell)-landIceBoundaryLayerTracers(indexBLT,iCell))) + landIceHeatFlux(iCell) = landIceFloatingFraction(iCell)*heatFlux + + heatFluxToLandIce(iCell) = 0.0_RKIND + + end do + !$omp end do + !$omp end parallel + endif ! isomipOn + + if (jenkinsOn .or. hollandJenkinsOn) then + if(useHollandJenkinsAdvDiff) then + ! melting solution + call compute_HJ99_melt_fluxes( & + landIceFloatingMask, & + landIceBoundaryLayerTracers(indexBLT,:), & + landIceBoundaryLayerTracers(indexBLS,:), & + landIceTracerTransferVelocities(indexHeatTrans,:), & + landIceTracerTransferVelocities(indexSaltTrans,:), & + landIceSurfaceTemperature, & + landIcePressure, & + landIceInterfaceTracers(indexIT,:), & + landIceInterfaceTracers(indexIS,:), & + landIceFreshwaterFlux, & + landIceHeatFlux, & + heatFluxToLandIce, & + nCells, & + err) + if(err .ne. 0) then + call mpas_log_write( & + 'compute_HJ99_melt_fluxes failed.', & + MPAS_LOG_CRIT) + end if + + ! freezing solution + call compute_melt_fluxes( & + landIceFloatingMask, & + landIceBoundaryLayerTracers(indexBLT,:), & + landIceBoundaryLayerTracers(indexBLS,:), & + landIceTracerTransferVelocities(indexHeatTrans,:), & + landIceTracerTransferVelocities(indexSaltTrans,:), & + landIcePressure, & + freezeInterfaceTemperature, & + freezeInterfaceSalinity, & + freezeFreshwaterFlux, & + freezeHeatFlux, & + freezeIceHeatFlux, & + nCells, & + err) + if(err .ne. 0) then + call mpas_log_write( & + 'compute_melt_fluxes failed.', & + MPAS_LOG_CRIT) + end if + + do iCell = 1, nCells + if ((landIceFloatingMask(iCell) == 0) .or. (landIceFreshwaterFlux(iCell) >= 0.0_RKIND)) cycle + + landIceInterfaceTracers(indexIS,iCell) = freezeInterfaceSalinity(iCell) + landIceInterfaceTracers(indexIT,iCell) = freezeInterfaceTemperature(iCell) + landIceFreshwaterFlux(iCell) = freezeFreshwaterFlux(iCell) + landIceHeatFlux(iCell) = freezeHeatFlux(iCell) + heatFluxToLandIce(iCell) = freezeIceHeatFlux(iCell) + end do + else ! not using Holland and Jenkins advection/diffusion + call compute_melt_fluxes( & + landIceFloatingMask, & + landIceBoundaryLayerTracers(indexBLT,:), & + landIceBoundaryLayerTracers(indexBLS,:), & + landIceTracerTransferVelocities(indexHeatTrans,:), & + landIceTracerTransferVelocities(indexSaltTrans,:), & + landIcePressure, & + landIceInterfaceTracers(indexIT,:), & + landIceInterfaceTracers(indexIS,:), & + landIceFreshwaterFlux, & + landIceHeatFlux, & + heatFluxToLandIce, & + nCells, & + err) + if(err .ne. 0) then + call mpas_log_write( & + 'compute_melt_fluxes failed.', & + MPAS_LOG_CRIT) + end if end if + ! modulate the fluxes by the landIceFloatingFraction do iCell = 1, nCells - if ((landIceFloatingMask(iCell) == 0) .or. (landIceFreshwaterFlux(iCell) >= 0.0_RKIND)) cycle + if (landIceFloatingMask(iCell) == 0) cycle - landIceInterfaceTracers(indexIS,iCell) = freezeInterfaceSalinity(iCell) - landIceInterfaceTracers(indexIT,iCell) = freezeInterfaceTemperature(iCell) - landIceFreshwaterFlux(iCell) = freezeFreshwaterFlux(iCell) - landIceHeatFlux(iCell) = freezeHeatFlux(iCell) - heatFluxToLandIce(iCell) = freezeIceHeatFlux(iCell) + landIceFreshwaterFlux(iCell) = landIceFloatingFraction(iCell)*landIceFreshwaterFlux(iCell) + landIceHeatFlux(iCell) = landIceFloatingFraction(iCell)*landIceHeatFlux(iCell) + heatFluxToLandIce(iCell) = landIceFloatingFraction(iCell)*heatFluxToLandIce(iCell) end do - else ! not using Holland and Jenkins advection/diffusion - call compute_melt_fluxes( & - landIceFloatingMask, & - landIceBoundaryLayerTracers(indexBLT,:), & - landIceBoundaryLayerTracers(indexBLS,:), & - landIceTracerTransferVelocities(indexHeatTrans,:), & - landIceTracerTransferVelocities(indexSaltTrans,:), & - landIcePressure, & - landIceInterfaceTracers(indexIT,:), & - landIceInterfaceTracers(indexIS,:), & - landIceFreshwaterFlux, & - landIceHeatFlux, & - heatFluxToLandIce, & - nCells, & - err) - if(err .ne. 0) then - call mpas_log_write( & - 'compute_melt_fluxes failed.', & - MPAS_LOG_CRIT) - end if - end if - ! modulate the fluxes by the landIceFloatingFraction - do iCell = 1, nCells - if (landIceFloatingMask(iCell) == 0) cycle + endif ! jenkinsOn or hollandJenkinsOn - landIceFreshwaterFlux(iCell) = landIceFloatingFraction(iCell)*landIceFreshwaterFlux(iCell) - landIceHeatFlux(iCell) = landIceFloatingFraction(iCell)*landIceHeatFlux(iCell) - heatFluxToLandIce(iCell) = landIceFloatingFraction(iCell)*heatFluxToLandIce(iCell) - end do + if(useHollandJenkinsAdvDiff) then + deallocate(freezeInterfaceSalinity, & + freezeInterfaceTemperature, & + freezeFreshwaterFlux, & + freezeHeatFlux, & + freezeIceHeatFlux) + end if - endif ! jenkinsOn or hollandJenkinsOn + endif ! landIceStandaloneOn ! Add frazil and interface melt/freeze to get total fresh water flux if ( associated(frazilIceFreshwaterFlux) ) then @@ -666,14 +679,6 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & end do end if - if(useHollandJenkinsAdvDiff) then - deallocate(freezeInterfaceSalinity, & - freezeInterfaceTemperature, & - freezeFreshwaterFlux, & - freezeHeatFlux, & - freezeIceHeatFlux) - end if - call mpas_timer_stop("land_ice_build_arrays") !--------------------------------------------------------------------