From 40d5114cc7ebae278cf299f852d2580a8889b4c7 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 | 287 +++++++++--------- 1 file changed, 145 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..fe489029a139 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,179 @@ 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 + call mpas_pool_get_array(forcingPool, 'landIceFloatingMask', landIceFloatingMask) 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) + call mpas_pool_get_array(forcingPool, 'landIceFloatingFraction', landIceFloatingFraction) + call mpas_pool_get_array(forcingPool, 'landIceHeatFlux', landIceHeatFlux) + call mpas_pool_get_array(forcingPool, 'heatFluxToLandIce', heatFluxToLandIce) + + + 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 - ! 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.) + if (isomipOn) then !*** ISOMIP formulation - ! 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)) + !$omp parallel + !$omp do schedule(runtime) private(freshwaterFlux, heatFlux) + do iCell = 1, nCells + if (landIceFloatingMask(iCell) == 0) cycle - landIceFreshwaterFlux(iCell) = landIceFloatingFraction(iCell)*freshwaterFlux + ! 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 (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 + ! 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)) - heatFluxToLandIce(iCell) = 0.0_RKIND + landIceFreshwaterFlux(iCell) = landIceFloatingFraction(iCell)*freshwaterFlux - end do - !$omp end do - !$omp end parallel - endif ! isomipOn + ! 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 - 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 + heatFluxToLandIce(iCell) = 0.0_RKIND - ! 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 do + !$omp end do + !$omp end parallel + end if ! 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 + end if ! 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 + end if ! landIceStandaloneOn ! Add frazil and interface melt/freeze to get total fresh water flux if ( associated(frazilIceFreshwaterFlux) ) then @@ -666,14 +677,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") !--------------------------------------------------------------------