Skip to content

Surface thickness flux spreading in MPAS-O #323

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
maltrud wants to merge 2 commits into from
Closed

Surface thickness flux spreading in MPAS-O #323

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
06ef567
allow vertical spreading of surface thickness fluxes.
maltrud Jun 27, 2019
ab28828
fix bug in the KPP nonLocal term
maltrud Jul 10, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions src/core_ocean/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -724,6 +724,10 @@
description="The length scale of exponential decay of river runoff. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
possible_values="Any positive real number."
/>
<nml_option name="config_flux_attenuation_coefficient_prec_evap_ice" type="real" default_value="0.001" units="m"
description="The length scale of exponential decay of non-river runoff thickness fluxes. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
possible_values="Any positive real number."
/>
</nml_record>
<nml_record name="time_varying_forcing" mode="init;forward">
<nml_option name="config_use_time_varying_atmospheric_forcing" type="logical" default_value=".false." units="unitless"
Expand Down Expand Up @@ -2761,6 +2765,9 @@
<var name="surfaceFluxAttenuationCoefficientRunoff" type="real" dimensions="nCells Time" units="m"
description="The spatially-dependent length scale of exponential decay of river runoff. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
/>
<var name="surfaceFluxAttenuationCoefficientPrecEvapIce" type="real" dimensions="nCells Time" units="m"
description="The spatially-dependent length scale of exponential decay of non-river runoff thickness fluxes. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
/>
<!-- diagnostic fields for land-ice fluxes -->
<var name="landIceFrictionVelocity" type="real" dimensions="nCells Time" units="m s^{-1}"
description="The friction velocity $u_*$ under land ice"
Expand Down Expand Up @@ -2917,6 +2924,10 @@
description="Divergence of transmission through interfaces of surface fluxes below the surface layer at cell centers. These are applied only to river runoff."
packages="forwardMode;analysisMode"
/>
<var name="fractionAbsorbedPrecEvapIce" type="real" dimensions="nVertLevels nCells Time" units="fractional"
description="Divergence of transmission through interfaces of surface fluxes below the surface layer at cell centers. These are applied only to non-river runoff thickness fluxes."
packages="forwardMode;analysisMode"
/>

<!-- Input fields for coupling -->
<!-- Coupling fields associated with heat fluxes -->
Expand Down
34 changes: 26 additions & 8 deletions src/core_ocean/shared/mpas_ocn_diagnostics.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,8 +153,10 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic

real (kind=RKIND), pointer :: config_flux_attenuation_coefficient
real (kind=RKIND), pointer :: config_flux_attenuation_coefficient_runoff
real (kind=RKIND), pointer :: config_flux_attenuation_coefficient_prec_evap_ice
real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficient
real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientRunoff
real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientPrecEvapIce

real (kind=RKIND) :: areaTri1, layerThicknessVertexInv, tempRiVal, dcEdge_temp, dvEdge_temp, weightsOnEdge_temp
real (kind=RKIND), dimension(:), allocatable:: layerThicknessVertexVec
Expand All @@ -176,6 +178,8 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient', config_flux_attenuation_coefficient)
call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_runoff', &
config_flux_attenuation_coefficient_runoff)
call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_prec_evap_ice', &
config_flux_attenuation_coefficient_prec_evap_ice)
call mpas_pool_get_config(ocnConfigs, 'config_use_wetting_drying', config_use_wetting_drying)
call mpas_pool_get_config(ocnConfigs, 'config_thickness_flux_type', config_thickness_flux_type)

Expand Down Expand Up @@ -261,6 +265,8 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic

call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficient', surfaceFluxAttenuationCoefficient)
call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficientRunoff', surfaceFluxAttenuationCoefficientRunoff)
call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficientPrecEvapIce', &
surfaceFluxAttenuationCoefficientPrecEvapIce)
!
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed layerThicknessEdge value
Expand Down Expand Up @@ -852,6 +858,7 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
do iCell = 1, nCells
surfaceFluxAttenuationCoefficient(iCell) = config_flux_attenuation_coefficient
surfaceFluxAttenuationCoefficientRunoff(iCell) = config_flux_attenuation_coefficient_runoff
surfaceFluxAttenuationCoefficientPrecEvapIce(iCell) = config_flux_attenuation_coefficient_prec_evap_ice
end do
!$omp end do

Expand Down Expand Up @@ -1367,8 +1374,9 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaCell
real (kind=RKIND), dimension(:), pointer :: penetrativeTemperatureFlux, surfaceThicknessFlux, &
surfaceBuoyancyForcing, surfaceFrictionVelocity, penetrativeTemperatureFluxOBL, &
normalVelocitySurfaceLayer, surfaceThicknessFluxRunoff
real (kind=RKIND), pointer :: config_flux_attenuation_coefficient, config_flux_attenuation_coefficient_runoff
normalVelocitySurfaceLayer, surfaceThicknessFluxRunoff, rainFlux, evaporationFlux
real (kind=RKIND), pointer :: config_flux_attenuation_coefficient, config_flux_attenuation_coefficient_runoff, &
config_flux_attenuation_coefficient_prec_evap_ice

real (kind=RKIND), dimension(:), pointer :: surfaceStress, surfaceStressMagnitude
real (kind=RKIND), dimension(:,:), pointer :: &
Expand All @@ -1388,7 +1396,8 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
! local
integer :: iCell, iEdge, i, k, err, timeLevel
integer, pointer :: indexTempFlux, indexSaltFlux
real (kind=RKIND) :: numerator, denominator, turbulentVelocitySquared, fracAbsorbed, fracAbsorbedRunoff
real (kind=RKIND) :: numerator, denominator, turbulentVelocitySquared, fracAbsorbedRunoff, &
fracAbsorbedPrecEvapIce
real (kind=RKIND) :: factor, deltaVelocitySquared, delU2, invAreaCell

type (field2DReal), pointer :: densitySurfaceDisplacedField, thermalExpansionCoeffField, salineContractionCoeffField
Expand All @@ -1406,6 +1415,8 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient', config_flux_attenuation_coefficient)
call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_runoff', &
config_flux_attenuation_coefficient_runoff)
call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_prec_evap_ice', &
config_flux_attenuation_coefficient_prec_evap_ice)
! set the parameter turbulentVelocitySquared
turbulentVelocitySquared = 0.001_RKIND

Expand Down Expand Up @@ -1441,6 +1452,8 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
call mpas_pool_get_array(forcingPool, 'penetrativeTemperatureFlux', penetrativeTemperatureFlux)
call mpas_pool_get_array(forcingPool, 'surfaceStress', surfaceStress)
call mpas_pool_get_array(forcingPool, 'surfaceStressMagnitude', surfaceStressMagnitude)
call mpas_pool_get_array(forcingPool, 'rainFlux', rainFlux)
call mpas_pool_get_array(forcingPool, 'evaporationFlux', evaporationFlux)

call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1)
call mpas_pool_get_array(tracersSurfaceFluxPool, 'activeTracersSurfaceFlux', activeTracersSurfaceFlux)
Expand Down Expand Up @@ -1489,19 +1502,24 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
!

! Compute fraction of thickness flux that is in the top model layer
fracAbsorbed = 1.0_RKIND - exp( max(-layerThickness(1, iCell) / config_flux_attenuation_coefficient, -100.0_RKIND) )
fracAbsorbedRunoff = 1.0_RKIND - exp( max(-layerThickness(1, iCell) / config_flux_attenuation_coefficient_runoff, &
-100.0_RKIND) )
fracAbsorbedPrecEvapIce = 1.0_RKIND - exp( max(-layerThickness(1, iCell) / &
config_flux_attenuation_coefficient_prec_evap_ice, -100.0_RKIND) )

Copy link
Copy Markdown
Contributor

@vanroekel vanroekel Jul 10, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we need a third fracAbsorbed for all other thickness fluxes that are not being spread, otherwise we will not have a balance of heat fluxes and the corresponding cooling from the layer expansion from the thickness flux. Similar to my comment below.

! Store the total tracer flux below in nonLocalSurfaceTemperatureFlux for use in the CVMix nonlocal
! transport code. This includes tracer forcing due to thickness

nonLocalSurfaceTracerFlux(indexTempFlux, iCell) = activeTracersSurfaceFlux(indexTempFlux,iCell) &
+ penetrativeTemperatureFlux(iCell) - penetrativeTemperatureFluxOBL(iCell) - fracAbsorbed * &
surfaceThicknessFlux(iCell) * activeTracers(indexTempFlux,1,iCell) + &
activeTracersSurfaceFluxRunoff(indexTempFlux,iCell) * fracAbsorbedRunoff
+ penetrativeTemperatureFlux(iCell) - penetrativeTemperatureFluxOBL(iCell) &
- fracAbsorbedPrecEvapIce * (rainFlux(iCell) + evaporationFlux(iCell)) * &
activeTracers(indexTempFlux,1,iCell)/rho_sw &
- fracAbsorbedRunoff * surfaceThicknessFluxRunoff(iCell) * ( &
+ min(activeTracers(indexTempFlux,1,iCell),0.0_RKIND) )/rho_sw
Copy link
Copy Markdown
Collaborator

@xylar xylar Jul 10, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't freshwater fluxes from land ice also be included here?

Copy link
Copy Markdown
Collaborator

@xylar xylar Jul 10, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess sea-ice fluxes aren't included so maybe not, but what's special about rain and evap?

Copy link
Copy Markdown
Collaborator

@xylar xylar Jul 10, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm pretty sure we want to include all these other surfaceThicknessFlux terms that this change removes:
https://github.com/MPAS-Dev/MPAS-Model/blob/ocean/develop/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F#L299-L300
https://github.com/MPAS-Dev/MPAS-Model/blob/ocean/develop/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F#L269

Copy link
Copy Markdown
Collaborator

@xylar xylar Jul 10, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@vanroekel?

Copy link
Copy Markdown
Contributor

@vanroekel vanroekel Jul 10, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In looking at this a bit, I agree with @xylar here. I think the problem is that @maltrud is talking about the thickness flux spreading, but Xylar raises an important point that there is still a bug in the non local flux for KPP, which is more appropriate to #305 . Landice and seaice freshwater do get added to the heat flux and thickness flux such that the heat content doesn't change, but we need to correct for that additional flux to retain a net zero flux for KPP. As the code is right now a seaice or landice freshwater flux will create an artificially strong surface buoyancy flux which would reduce mixing.


nonLocalSurfaceTracerFlux(indexSaltFlux,iCell) = activeTracersSurfaceFlux(indexSaltFlux,iCell) &
- fracAbsorbed * surfaceThicknessFlux(iCell) * activeTracers(indexSaltFlux,1,iCell)
- fracAbsorbedPrecEvapIce * surfaceThicknessFlux(iCell) * activeTracers(indexSaltFlux,1,iCell) &
- fracAbsorbedRunoff * surfaceThicknessFluxRunoff(iCell) * activeTracers(indexSaltFlux,1,iCell)

surfaceBuoyancyForcing(iCell) = thermalExpansionCoeff (1,iCell) &
* nonLocalSurfaceTracerFlux(indexTempFlux,iCell) &
Expand Down
26 changes: 25 additions & 1 deletion src/core_ocean/shared/mpas_ocn_forcing.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,9 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, diagno

real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficient
real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientRunoff
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, fractionAbsorbed, fractionAbsorbedRunoff
real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientPrecEvapIce
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, fractionAbsorbed, fractionAbsorbedRunoff, &
fractionAbsorbedPrecEvapIce

integer :: iCell, k, timeLevel, nCells

Expand All @@ -137,11 +139,17 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, diagno
call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficientRunoff', &
surfaceFluxAttenuationCoefficientRunoff)

call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficientPrecEvapIce', &
surfaceFluxAttenuationCoefficientPrecEvapIce)

call mpas_pool_get_array(forcingPool, 'fractionAbsorbed', fractionAbsorbed)
call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', fractionAbsorbedRunoff)
call mpas_pool_get_array(forcingPool, 'fractionAbsorbedPrecEvapIce', fractionAbsorbedPrecEvapIce)

nCells = nCellsArray( 2 )

! first do surface fluxes that are not associated with thickness fluxes (radiative, sensible, etc)

do iCell = 1, nCells
zTop = 0.0_RKIND
transmissionCoeffTop = ocn_forcing_transmission(zTop, surfaceFluxAttenuationCoefficient(iCell))
Expand Down Expand Up @@ -172,6 +180,22 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, diagno
end do
end do

! now do non-runoff thickness fluxes

do iCell = 1, nCells
zTop = 0.0_RKIND
transmissionCoeffTop = ocn_forcing_transmission(zTop, surfaceFluxAttenuationCoefficientPrecEvapIce(iCell))
do k = 1, maxLevelCell(iCell)
zBot = zTop - layerThickness(k,iCell)
transmissionCoeffBot = ocn_forcing_transmission(zBot, surfaceFluxAttenuationCoefficientPrecEvapIce(iCell))

fractionAbsorbedPrecEvapIce(k, iCell) = transmissionCoeffTop - transmissionCoeffBot

zTop = zBot
transmissionCoeffTop = transmissionCoeffBot
end do
end do

end subroutine ocn_forcing_build_fraction_absorbed_array!}}}

!***********************************************************************
Expand Down
6 changes: 4 additions & 2 deletions src/core_ocean/shared/mpas_ocn_tendency.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,8 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{
real (kind=RKIND), dimension(:), pointer :: surfaceThicknessFlux
real (kind=RKIND), dimension(:), pointer :: surfaceThicknessFluxRunoff
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge, &
vertAleTransportTop, tend_layerThickness, normalTransportVelocity, fractionAbsorbed, fractionAbsorbedRunoff
vertAleTransportTop, tend_layerThickness, normalTransportVelocity, fractionAbsorbed, fractionAbsorbedRunoff, &
fractionAbsorbedPrecEvapIce

integer, pointer :: nCells
integer :: err, iCell
Expand All @@ -138,6 +139,7 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{
call mpas_pool_get_array(forcingPool, 'surfaceThicknessFluxRunoff', surfaceThicknessFluxRunoff)
call mpas_pool_get_array(forcingPool, 'fractionAbsorbed', fractionAbsorbed)
call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', fractionAbsorbedRunoff)
call mpas_pool_get_array(forcingPool, 'fractionAbsorbedPrecEvapIce', fractionAbsorbedPrecEvapIce)

!
! height tendency: start accumulating tendency terms
Expand Down Expand Up @@ -179,7 +181,7 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{
! surface flux tendency
!

call ocn_thick_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbsorbedRunoff, layerThickness, &
call ocn_thick_surface_flux_tend(meshPool, fractionAbsorbedPrecEvapIce, fractionAbsorbedRunoff, layerThickness, &
surfaceThicknessFlux, surfaceThicknessFluxRunoff, tend_layerThickness, err)

!
Expand Down
19 changes: 10 additions & 9 deletions src/core_ocean/shared/mpas_ocn_thick_surface_flux.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ module ocn_thick_surface_flux
!
!-----------------------------------------------------------------------

subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, transmissionCoefficientsRunoff, &
subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficientsPrecEvapIce, transmissionCoefficientsRunoff, &
layerThickness, surfaceThicknessFlux, surfaceThicknessFluxRunoff, tend, err)!{{{
!-----------------------------------------------------------------
!
Expand All @@ -83,7 +83,7 @@ subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, trans
meshPool !< Input: mesh information

real (kind=RKIND), dimension(:,:), intent(in) :: &
transmissionCoefficients, &!< Input: Coefficients for the transmission of surface fluxes
transmissionCoefficientsPrecEvapIce, &!< Input: Coefficients for transmission of non-runoff surface thickness fluxes
transmissionCoefficientsRunoff !< Input: Coefficients for the transmission of surface fluxes due to river runoff

real (kind=RKIND), dimension(:,:), intent(in) :: &
Expand Down Expand Up @@ -122,7 +122,7 @@ subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, trans
integer, dimension(:), pointer :: nCellsArray
integer, dimension(:), pointer :: maxLevelCell

real (kind=RKIND) :: remainingFlux, remainingFluxRunoff
real (kind=RKIND) :: remainingFluxPrecEvapIce, remainingFluxRunoff

err = 0

Expand All @@ -136,20 +136,21 @@ subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, trans

nCells = nCellsArray( 1 )

!$omp do schedule(runtime) private(remainingFlux, remainingFluxRunoff, k)
!$omp do schedule(runtime) private(remainingFluxPrecEvapIce, remainingFluxRunoff, k)
do iCell = 1, nCells
remainingFlux = 1.0_RKIND
remainingFluxPrecEvapIce = 1.0_RKIND
remainingFluxRunoff = 1.0_RKIND
do k = 1, maxLevelCell(iCell)
remainingFlux = remainingFlux - transmissionCoefficients(k, iCell)
remainingFluxPrecEvapIce = remainingFluxPrecEvapIce - transmissionCoefficientsPrecEvapIce(k, iCell)
remainingFluxRunoff = remainingFluxRunoff - transmissionCoefficientsRunoff(k, iCell)

tend(k, iCell) = tend(k, iCell) + surfaceThicknessFlux(iCell) * transmissionCoefficients(k, iCell) &
tend(k, iCell) = tend(k, iCell) + surfaceThicknessFlux(iCell) * transmissionCoefficientsPrecEvapIce(k, iCell) &
+ surfaceThicknessFluxRunoff(iCell) * transmissionCoefficientsRunoff(k, iCell)
end do

if(maxLevelCell(iCell) > 0 .and. remainingFlux > 0.0_RKIND) then
tend(maxLevelCell(iCell), iCell) = tend(maxLevelCell(iCell), iCell) + remainingFlux * surfaceThicknessFlux(iCell)
if(maxLevelCell(iCell) > 0 .and. remainingFluxPrecEvapIce > 0.0_RKIND) then
tend(maxLevelCell(iCell), iCell) = tend(maxLevelCell(iCell), iCell) &
+ remainingFluxPrecEvapIce * surfaceThicknessFlux(iCell)
end if

if(maxLevelCell(iCell) > 0 .and. remainingFluxRunoff > 0.0_RKIND) then
Expand Down