diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_vertical_grids.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_vertical_grids.F index 5c0ac5a0e114..d2c0033c0d0b 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_vertical_grids.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_vertical_grids.F @@ -591,7 +591,7 @@ subroutine ocn_compute_z_star_layerThickness(layerThickness, restingThickness, s integer, intent(out) :: iErr integer :: k - real (kind=RKIND) :: layerStretch + real (kind=RKIND) :: layerStretch, restingThicknessSum iErr = 0 @@ -599,7 +599,10 @@ subroutine ocn_compute_z_star_layerThickness(layerThickness, restingThickness, s if (maxLevelCell<=0) return - layerStretch = (ssh + bottomDepth)/bottomDepth + do k = minLevelCell, maxLevelCell + restingThicknessSum = restingThicknessSum + restingThickness(k) + end do + layerStretch = (ssh + bottomDepth)/restingThicknessSum do k = minLevelCell, maxLevelCell layerThickness(k) = layerStretch * restingThickness(k) end do diff --git a/components/mpas-ocean/src/shared/mpas_ocn_mesh.F b/components/mpas-ocean/src/shared/mpas_ocn_mesh.F index 5e0690704e40..de21287093f8 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_mesh.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_mesh.F @@ -1568,7 +1568,7 @@ subroutine ocn_meshVertCoordInit(domain)!{{{ select case (trim(config_vert_coord_movement)) case ('fixed') - vertCoordMovementWeights = 0.0_RKIND + vertCoordMovementWeights = 1e-14_RKIND vertCoordMovementWeights(1) = 1.0_RKIND case ('uniform_stretching') @@ -1587,7 +1587,7 @@ subroutine ocn_meshVertCoordInit(domain)!{{{ (config_vert_taper_weight_depth_2 - & config_vert_taper_weight_depth_1) ! limit range to (0,1) - vertCoordMovementWeights(k) = max( 0.0_RKIND, & + vertCoordMovementWeights(k) = max( 1e-14_RKIND, & min( 1.0_RKIND, vertCoordMovementWeights(k) ) ) end do diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index a16e52f14499..89cfa07ad2bf 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -419,7 +419,7 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, & vertAleTransportTop, tendVel, err) ! Add pressure gradient - call ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & + call ocn_vel_pressure_grad_tend(ssh, gradSSH, pressure, surfacePressure, & montgomeryPotential, zMid, & density, potentialDensity, & indxTemp, indxSalt, activeTracers, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F b/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F index a4e861f6ed67..9ebe1b3645d2 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F @@ -125,10 +125,12 @@ subroutine ocn_ALE_thickness(verticalMeshPool, SSH, ALE_thickness, & nCells ! number of cells real (kind=RKIND) :: & - weightSum, &! sum of weights in vertical - thicknessSum, &! total thickness - remainder, &! track remainder in mix/max alteration - newThickness ! temp used during min/max adjustment + weightSum, &! sum of weights in vertical + restingThicknessSum, &! total resting thickness + restingThicknessDiff, &! difference from total resting thickness + weightedThicknessSum, &! total resting thickness, weighted + remainder, &! track remainder in mix/max alteration + newThickness ! temp used during min/max adjustment real (kind=RKIND), dimension(:), allocatable :: & prelim_ALE_thickness, & ! ALE thickness at new time @@ -170,33 +172,39 @@ subroutine ocn_ALE_thickness(verticalMeshPool, SSH, ALE_thickness, & #ifdef MPAS_OPENACC !xacc parallel loop & !xacc present(ALE_thickness, SSH, restingThickness, & - !xacc minLevelCell, maxLevelCell, & + !xacc minLevelCell, maxLevelCell, bottomDepth, & !xacc vertCoordMovementWeights) & - !xacc private(k, kMin, kMax, thicknessSum) + !xacc private(k, kMin, kMax, restingThicknessDiff, & + !xacc restingThicknessSum, weightedThicknessSum) #else !$omp parallel !$omp do schedule(runtime) & - !$omp private(k, kMin, kMax, thicknessSum) + !$omp private(k, kMin, kMax, restingThicknessDiff, & + !$omp restingThicknessSum, weightedThicknessSum) #endif do iCell = 1, nCells kMax = maxLevelCell(iCell) kMin = minLevelCell(iCell) - thicknessSum = 1e-14_RKIND + restingThicknessSum = 0.0_RKIND + weightedThicknessSum = 0.0_RKIND do k = kMin, kMax - thicknessSum = thicknessSum & + restingThicknessSum = restingThicknessSum & + + restingThickness(k,iCell) + weightedThicknessSum = weightedThicknessSum & + vertCoordMovementWeights(k) & * restingThickness(k,iCell) end do + restingThicknessDiff = restingThicknessSum - bottomDepth(iCell) ! Note that restingThickness is nonzero, and remaining ! terms are perturbations about zero. ! This is equation 4 and 6 in Petersen et al 2015, ! but with eqn 6 do k = kMin, kMax - ALE_thickness(k,iCell) = restingThickness(k,iCell) & - + (SSH(iCell)*vertCoordMovementWeights(k)* & - restingThickness(k,iCell) )/thicknessSum + ALE_thickness(k, iCell) = restingThickness(k, iCell) & + + ( (SSH(iCell) - restingThicknessDiff) * vertCoordMovementWeights(k) * & + restingThickness(k, iCell) ) / weightedThicknessSum end do enddo #ifndef MPAS_OPENACC diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F index b8f83cd01dce..f7e97aa13384 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_pressure_grad.F @@ -94,7 +94,7 @@ module ocn_vel_pressure_grad ! !----------------------------------------------------------------------- - subroutine ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & + subroutine ocn_vel_pressure_grad_tend(ssh, gradSSH, pressure, surfacePressure, & montgomeryPotential, zMid, & density, potentialDensity, & indxT, indxS, tracers, & @@ -113,6 +113,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & real (kind=RKIND), dimension(:), intent(in) :: & ssh, &!< [in] sea surface height + gradSSH, &!< [in] sea surface height gradient on edges surfacePressure !< [in] surface pressure real (kind=RKIND), dimension(:,:), intent(in) :: & @@ -241,7 +242,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, & do k=kMin,kMax tend(k,iEdge) = tend(k,iEdge) - & edgeMask(k,iEdge)*invdcEdge* & - ( gravity*(ssh(cell2) - ssh(cell1)) & + ( gravity*gradSSH(iEdge) & + density0Inv*( surfacePressure(cell2) & - surfacePressure(cell1)) ) end do