Skip to content

Commit

Permalink
Fix total land-ice freshwater flux in data mode
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@xylar
xylar committed Oct 30, 2024
1 parent 96ef265 commit 13b1062
Showing 1 changed file with 147 additions and 142 deletions.
289 changes: 147 additions & 142 deletions components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
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)

! 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)

! 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
Expand All @@ -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")

!--------------------------------------------------------------------
Expand Down

0 comments on commit 13b1062

Please sign in to comment.