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 Nov 1, 2024
1 parent 834203f commit 40d5114
Showing 1 changed file with 145 additions and 142 deletions.
287 changes: 145 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,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
Expand All @@ -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")

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

0 comments on commit 40d5114

Please sign in to comment.