From df704c755a6701af6115bd27bc813eaf3ccc4ae6 Mon Sep 17 00:00:00 2001 From: Darin Scott Comeau Date: Tue, 8 May 2018 14:57:24 -0600 Subject: [PATCH] Adding iceberg terms to sea ice momentum equation. --- .../shared/mpas_seaice_berg_velocity_solver.F | 4 ++++ .../shared/mpas_seaice_velocity_solver.F | 15 +++++++++++---- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/src/core_seaice/shared/mpas_seaice_berg_velocity_solver.F b/src/core_seaice/shared/mpas_seaice_berg_velocity_solver.F index 34d64ba871..1fbda0a977 100644 --- a/src/core_seaice/shared/mpas_seaice_berg_velocity_solver.F +++ b/src/core_seaice/shared/mpas_seaice_berg_velocity_solver.F @@ -1846,6 +1846,10 @@ subroutine seaice_berg_forcing_for_ice(domain) call MPAS_pool_get_array(bergVelocitySolverPool, "bergStressU", bergStressU) call MPAS_pool_get_array(bergVelocitySolverPool, "bergStressV", bergStressV) + ! initialize + bergStressU(:) = 0.0_RKIND + bergStressV(:) = 0.0_RKIND + ! interpolate cell to vertex call seaice_interpolate_cell_to_vertex(& meshPool, & diff --git a/src/core_seaice/shared/mpas_seaice_velocity_solver.F b/src/core_seaice/shared/mpas_seaice_velocity_solver.F index e960bc82e7..fed451410a 100644 --- a/src/core_seaice/shared/mpas_seaice_velocity_solver.F +++ b/src/core_seaice/shared/mpas_seaice_velocity_solver.F @@ -2731,7 +2731,8 @@ subroutine solve_velocity_revised(domain) type(MPAS_pool_type), pointer :: & velocitySolverPool, & - icestatePool + icestatePool, & + bergVelocitySolverPool integer, dimension(:), pointer :: & solveVelocity @@ -2751,7 +2752,9 @@ subroutine solve_velocity_revised(domain) surfaceTiltForceV, & oceanStressU, & oceanStressV, & - oceanStressCoeff + oceanStressCoeff, & + bergStressU, & + bergStressV real(kind=RKIND), pointer :: & dynamicsTimeStep @@ -2778,6 +2781,7 @@ subroutine solve_velocity_revised(domain) call MPAS_pool_get_subpool(block % structs, "velocity_solver", velocitySolverPool) call MPAS_pool_get_subpool(block % structs, "icestate", icestatePool) + call MPAS_pool_get_subpool(block % structs, "berg_velocity_solver", bergVelocitySolverPool) call MPAS_pool_get_array(icestatePool, "totalMassVertex", totalMassVertex) @@ -2798,6 +2802,9 @@ subroutine solve_velocity_revised(domain) call MPAS_pool_get_array(velocitySolverPool, "oceanStressCoeff", oceanStressCoeff) call MPAS_pool_get_array(velocitySolverPool, "dynamicsTimeStep", dynamicsTimeStep) + call MPAS_pool_get_array(bergVelocitySolverPool, "bergStressU", bergStressU) + call MPAS_pool_get_array(bergVelocitySolverPool, "bergStressV", bergStressV) + do iVertex = 1, nVerticesSolve if (solveVelocity(iVertex) == 1) then @@ -2816,12 +2823,12 @@ subroutine solve_velocity_revised(domain) ! right hand side of matrix solve rightHandSide(1) = stressDivergenceU(iVertex) + airStressVertexU(iVertex) + surfaceTiltForceU(iVertex) + & - oceanStressCoeff(iVertex) * oceanStressU(iVertex) + & + oceanStressCoeff(iVertex) * oceanStressU(iVertex) + bergStressU(iVertex) + & (totalMassVertex(iVertex) * (numericalInertiaCoefficient * uVelocity(iVertex) + & uVelocityInitial(iVertex))) / dynamicsTimeStep rightHandSide(2) = stressDivergenceV(iVertex) + airStressVertexV(iVertex) + surfaceTiltForceV(iVertex) + & - oceanStressCoeff(iVertex) * oceanStressV(iVertex) + & + oceanStressCoeff(iVertex) * oceanStressV(iVertex) + bergStressV(iVertex) + & (totalMassVertex(iVertex) * (numericalInertiaCoefficient * vVelocity(iVertex) + & vVelocityInitial(iVertex))) / dynamicsTimeStep