Skip to content

Commit

Permalink
Merge PR #288 'luke/ocean/gmBolus3Dvarying' into ocean/develop
Browse files Browse the repository at this point in the history 
In the GM bolus calculation, there is a specified diffusivity and phase
speed (see Ferrari et al. 2010). Currently MPAS-O assumes these values
are fixed in time and space. This PR allows separate flags to enable
2D varying (+time) phase speed (based on the first baroclinic mode phase
speed) and a separate flag for 3D (+time) varying diffusivity.
  • Loading branch information
@mark-petersen
mark-petersen committed Aug 16, 2019
2 parents c95f8f8 + 04874bf commit 0730253
Show file tree
Hide file tree
Showing 2 changed files with 115 additions and 17 deletions.
25 changes: 24 additions & 1 deletion Registry.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,6 +387,22 @@
description="If true, the standard GM for the tracer advection and mixing is turned on. This includes the redi portion which is captured in hmix."
possible_values=".true. or .false."
/>
<nml_option name="config_gm_lat_variable_c2" type="logical" default_value=".false." units="NA"
description="If true, spatially variable phase speed is used Following Ferrari et al (2010)"
possible_values=".true. or .false."
/>
<nml_option name="config_gm_kappa_lat_depth_variable" type="logical" default_value=".false." units="NA"
description="If true, spatitally and deph varying Bolus Kappa is used, following Danabasoglu et al (2007)"
possible_values=".true. or .false."
/>
<nml_option name="config_gm_min_stratification_ratio" type="real" default_value="0.1" units="NA"
description="minimum value for N2/Nmax2 ratio"
possible_values="small positive reals"
/>
<nml_option name="config_gm_min_phase_speed" type="real" default_value="0.1" units="m s^{-1}"
description="minimum allowed phase speed for spatially variable computation"
possible_values="small positive reals"
/>
<nml_option name="config_use_Redi_surface_layer_tapering" type="logical" default_value=".false." units="NA"
description="If true, the Redi K33 vertical mixing is limited in and just below the ocean boundary layer."
possible_values=".true. or .false."
Expand Down Expand Up @@ -2600,11 +2616,18 @@
description="Meridional Component of the gradient of sea surface height at cell centers."
packages="forwardMode;analysisMode"
/>

<var name="normalGMBolusVelocity" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
description="Bolus velocity in Gent-McWilliams eddy parameterization"
packages="forwardMode;analysisMode"
/>
<var name="cGMphaseSpeed" type="real" dimensions="nEdges Time" units="m s^{-1}"
description="phase speed for the bolus velocity calculation"
packages="forwardMode;analysisMode"
/>
<var name="kappaGM3D" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
description="spatially and depth varying GM kappa"
packages="forwardMode;analysisMode"
/>
<var name="GMBolusVelocityX" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1}"
description="Bolus velocity in Gent-McWilliams eddy parameterization, x-direction"
packages="forwardMode;analysisMode"
Expand Down
107 changes: 91 additions & 16 deletions shared/mpas_ocn_gm.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,10 @@ module ocn_gm
logical, pointer :: config_use_Redi_bottom_layer_tapering
real (kind=RKIND), pointer :: config_Redi_surface_layer_tapering_extent
real (kind=RKIND), pointer :: config_Redi_bottom_layer_tapering_depth

logical, pointer :: config_gm_lat_variable_c2
logical, pointer :: config_gm_kappa_lat_depth_variable
real (kind=RKIND), pointer :: config_gm_min_stratification_ratio
real (kind=RKIND), pointer :: config_gm_min_phase_speed
real (kind=RKIND), parameter :: epsGM = 1.0e-12_RKIND

!***********************************************************************
Expand Down Expand Up @@ -102,14 +105,16 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
layerThicknessEdge, gradDensityEdge, gradDensityTopOfEdge, gradDensityConstZTopOfEdge, gradZMidEdge, &
gradZMidTopOfEdge, relativeSlopeTopOfEdge, relativeSlopeTopOfCell, k33, gmStreamFuncTopOfEdge, BruntVaisalaFreqTop, &
gmStreamFuncTopOfCell, dDensityDzTopOfEdge, dDensityDzTopOfCell, relativeSlopeTapering, relativeSlopeTaperingCell, &
areaCellSum
real(kind=RKIND), dimension(:), pointer :: boundaryLayerDepth, gmBolusKappa
areaCellSum, kappaGM3D

real(kind=RKIND), dimension(:), pointer :: boundaryLayerDepth, gmBolusKappa, cGMphaseSpeed, bottomDepth
real(kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, tridiagA, tridiagB, tridiagC, rightHandSide
integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell, nEdgesOnCell
integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
integer :: i, k, iEdge, cell1, cell2, iCell, N, iter
real(kind=RKIND) :: h1, h2, areaEdge, c, BruntVaisalaFreqTopEdge, rtmp, stmp, maxSlopeK33

real(kind=RKIND) :: bottomAv, sumN2, countN2, maxN, kappaSum, ltSum

! Dimensions
integer :: nCells, nEdges
integer, pointer :: nVertLevels
Expand All @@ -124,6 +129,8 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
call mpas_pool_get_array(diagnosticsPool, 'displacedDensity', displacedDensity)
call mpas_pool_get_array(diagnosticsPool, 'zMid', zMid)

call mpas_pool_get_array(diagnosticsPool, 'cGMphaseSpeed', cGMphaseSpeed)
call mpas_pool_get_array(diagnosticsPool, 'kappaGM3D', kappaGM3D)
call mpas_pool_get_array(diagnosticsPool, 'normalGMBolusVelocity', normalGMBolusVelocity)
call mpas_pool_get_array(diagnosticsPool, 'relativeSlopeTopOfEdge', relativeSlopeTopOfEdge)
call mpas_pool_get_array(diagnosticsPool, 'relativeSlopeTopOfCell', relativeSlopeTopOfCell)
Expand All @@ -141,6 +148,7 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
if (config_use_Redi_surface_layer_tapering) call mpas_pool_get_array(diagnosticsPool, 'boundaryLayerDepth', &
boundaryLayerDepth)

call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop)
call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell)
call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
Expand Down Expand Up @@ -516,7 +524,70 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
!
!--------------------------------------------------------------------

c = config_gravWaveSpeed_trunc**2

cGMphaseSpeed(:) = config_gravWaveSpeed_trunc
if (config_gm_lat_variable_c2) then
nEdges = nEdgesArray( 3 )
!$omp do schedule(runtime) private(cell1, cell2, k, ltSum, BruntVaisalaFreqTopEdge, sumN2, countN2, bottomAv)
do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
sumN2 = 0.0
countN2 = 0
ltSum = 0.0

do k=2,maxLevelEdgeTop(iEdge)

BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)

sumN2 = sumN2 + BruntVaisalaFreqTopEdge*layerThicknessEdge(k,iEdge)
ltSum = ltSum + layerThicknessEdge(k,iEdge)
countN2 = countN2 +1

enddo

if(countN2 > 0) cGMphaseSpeed(iEdge) = max(config_gm_min_phase_speed ,sqrt(sumN2/ltSum)*ltSum / 3.141592_RKIND)

enddo
!$omp end do

else
cGMphaseSpeed(:) = config_gravWaveSpeed_trunc
endif

!$omp do schedule(runtime)
do iEdge=1,nEdges
kappaGM3D(:,iEdge) = gmBolusKappa(iEdge)
enddo
!$omp end do

if (config_gm_kappa_lat_depth_variable) then

!$omp do schedule(runtime) private(cell1, cell2, k, BruntVaisalaFreqTopEdge, maxN)
do iEdge = 1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)

maxN = -1.0_RKIND
do k=2,maxLevelEdgeTop(iEdge)
BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)

maxN = max(maxN,BruntVaisalaFreqTopEdge)

enddo

do k=2,maxLevelEdgeTop(iEdge)
BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)

kappaGM3D(k,iEdge) = gmBolusKappa(iEdge)*max(config_gm_min_stratification_ratio, &
BruntVaisalaFreqTopEdge / (maxN + 1.0E-10_RKIND))
enddo
enddo
!$omp end do
endif

nEdges = nEdgesArray( 3 )

Expand All @@ -533,34 +604,34 @@ subroutine ocn_gm_compute_Bolus_velocity(diagnosticsPool, meshPool, scratchPool)
k = 2
BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1,iEdge) &
tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1,iEdge) &
* layerThicknessEdge(k,iEdge)) - BruntVaisalaFreqTopEdge
tridiagC(k-1) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k, iEdge) &
tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) &
/ (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge))
rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)

! Second to next to the last rows
do k = 3, maxLevelEdgeTop(iEdge)-1
BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
tridiagA(k-2) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k-1, iEdge) &
tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1, iEdge) &
/ (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge))
tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1, iEdge) &
tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) &
* layerThicknessEdge(k, iEdge) ) - BruntVaisalaFreqTopEdge
tridiagC(k-1) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k, iEdge) &
tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) &
/ (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge))
rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
end do

! Last row
k = maxLevelEdgeTop(iEdge)
BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND)
tridiagA(k-2) = 2.0_RKIND * config_gravWaveSpeed_trunc**2 / layerThicknessEdge(k-1,iEdge) &
tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1,iEdge) &
/ (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge))
tridiagB(k-1) = - 2.0_RKIND * config_gravWaveSpeed_trunc**2 / (layerThicknessEdge(k-1, iEdge) &
tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) &
* layerThicknessEdge(k, iEdge)) - BruntVaisalaFreqTopEdge
rightHandSide(k-1) = gmBolusKappa(iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)
rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge)

! Total number of rows
N = maxLevelEdgeTop(iEdge) - 1
Expand Down Expand Up @@ -747,7 +818,10 @@ subroutine ocn_gm_init(domain, err)!{{{
call mpas_pool_get_config(ocnConfigs, 'config_use_Redi_bottom_layer_tapering',config_use_Redi_bottom_layer_tapering)
call mpas_pool_get_config(ocnConfigs, 'config_Redi_surface_layer_tapering_extent',config_Redi_surface_layer_tapering_extent)
call mpas_pool_get_config(ocnConfigs, 'config_Redi_bottom_layer_tapering_depth',config_Redi_bottom_layer_tapering_depth)

call mpas_pool_get_config(ocnConfigs, 'config_gm_lat_variable_c2',config_gm_lat_variable_c2)
call mpas_pool_get_config(ocnConfigs, 'config_gm_kappa_lat_depth_variable', config_gm_kappa_lat_depth_variable)
call mpas_pool_get_config(ocnConfigs, 'config_gm_min_stratification_ratio', config_gm_min_stratification_ratio)
call mpas_pool_get_config(ocnConfigs, 'config_gm_min_phase_speed', config_gm_min_phase_speed)

block => domain % blocklist
do while (associated(block))
Expand Down Expand Up @@ -788,6 +862,7 @@ subroutine ocn_gm_init(domain, err)!{{{
call mpas_dmpar_finalize(domain % dminfo)
end if


block => block % next
end do
end subroutine ocn_gm_init!}}}
Expand Down

0 comments on commit 0730253

Please sign in to comment.