<p><b>mpetersen@lanl.gov</b> 2012-06-15 14:14:50 -0600 (Fri, 15 Jun 2012)</p><p>Merge branch vol_cons_RK_imp_mix back to trunk. A boundary update is required on the velocity field in RK4 after implicit vertical mixing. Also cleaned up the barotropic filter routines. Tested with RK4 zlevel 120km global by Mark, and isopycnal channel by Qingshan.<br>
</p><hr noshade><pre><font color="gray">
Property changes on: trunk/mpas
___________________________________________________________________
Modified: svn:mergeinfo
- /branches/atmos_physics:1672-1846
/branches/cam_mpas_nh:1260-1270
/branches/ocean_projects/ale_split_exp:1437-1483
/branches/ocean_projects/ale_vert_coord:1225-1383
/branches/ocean_projects/ale_vert_coord_new:1387-1428
/branches/ocean_projects/gmvar:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_error:1847-1887
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/monthly_forcing:1810-1867
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/ocean_projects/zstar_restart_new:1762-1770
/branches/omp_blocks/block_decomp:1374-1569
/branches/omp_blocks/ddt_reorg:1301-1414
/branches/omp_blocks/halo:1570-1638
/branches/omp_blocks/io:1639-1787
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+ /branches/atmos_physics:1672-1846
/branches/cam_mpas_nh:1260-1270
/branches/ocean_projects/ale_split_exp:1437-1483
/branches/ocean_projects/ale_vert_coord:1225-1383
/branches/ocean_projects/ale_vert_coord_new:1387-1428
/branches/ocean_projects/gmvar:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_error:1847-1887
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/monthly_forcing:1810-1867
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/ocean_projects/vol_cons_RK_imp_mix:1965-1992
/branches/ocean_projects/zstar_restart_new:1762-1770
/branches/omp_blocks/block_decomp:1374-1569
/branches/omp_blocks/ddt_reorg:1301-1414
/branches/omp_blocks/halo:1570-1638
/branches/omp_blocks/io:1639-1787
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
Modified: trunk/mpas/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-06-15 19:51:04 UTC (rev 1992)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-06-15 20:14:50 UTC (rev 1993)
@@ -70,7 +70,9 @@
ocn_diagnostic_solve, &
ocn_wtop, &
ocn_fuperp, &
- ocn_tendency_init
+ ocn_tendency_init, &
+ ocn_filter_btr_mode_u, &
+ ocn_filter_btr_mode_tend_u
!--------------------------------------------------------------------
!
@@ -1031,9 +1033,116 @@
end subroutine ocn_fuperp!}}}
+!***********************************************************************
+!
+! routine ocn_filter_btr_mode_u
+!
+!> \brief filters barotropic mode out of the velocity variable.
+!> \author Mark Petersen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine filters barotropic mode out of the velocity variable.
+!
+!-----------------------------------------------------------------------
+ subroutine ocn_filter_btr_mode_u(s, grid)!{{{
+ implicit none
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer :: iEdge, k, nEdges
+ real (kind=RKIND) :: vertSum, uhSum, hSum
+ real (kind=RKIND), dimension(:,:), pointer :: h_edge, u
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ call mpas_timer_start("ocn_filter_btr_mode_u")
+
+ u => s % u % array
+ h_edge => s % h_edge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ nEdges = grid % nEdges
+
+ do iEdge=1,nEdges
+
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
+ uhSum = h_edge(1,iEdge) * u(1,iEdge)
+ hSum = h_edge(1,iEdge)
+
+ do k=2,maxLevelEdgeTop(iEdge)
+ uhSum = uhSum + h_edge(k,iEdge) * u(k,iEdge)
+ hSum = hSum + h_edge(k,iEdge)
+ enddo
+
+ vertSum = uhSum/hSum
+ do k=1,maxLevelEdgeTop(iEdge)
+ u(k,iEdge) = u(k,iEdge) - vertSum
+ enddo
+ enddo ! iEdge
+
+ call mpas_timer_stop("ocn_filter_btr_mode_u")
+
+ end subroutine ocn_filter_btr_mode_u!}}}
+
!***********************************************************************
!
+! routine ocn_filter_btr_mode_tend_u
+!
+!> \brief ocn_filters barotropic mode out of the u tendency
+!> \author Mark Petersen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine filters barotropic mode out of the u tendency.
+!
+!-----------------------------------------------------------------------
+ subroutine ocn_filter_btr_mode_tend_u(tend, s, grid)!{{{
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer :: iEdge, k, nEdges
+ real (kind=RKIND) :: vertSum, uhSum, hSum
+ real (kind=RKIND), dimension(:,:), pointer :: h_edge, tend_u
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ call mpas_timer_start("ocn_filter_btr_mode_tend_u")
+
+ tend_u => tend % u % array
+ h_edge => s % h_edge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ nEdges = grid % nEdges
+
+ do iEdge=1,nEdges
+
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
+ uhSum = h_edge(1,iEdge) * tend_u(1,iEdge)
+ hSum = h_edge(1,iEdge)
+
+ do k=2,maxLevelEdgeTop(iEdge)
+ uhSum = uhSum + h_edge(k,iEdge) * tend_u(k,iEdge)
+ hSum = hSum + h_edge(k,iEdge)
+ enddo
+
+ vertSum = uhSum/hSum
+ do k=1,maxLevelEdgeTop(iEdge)
+ tend_u(k,iEdge) = tend_u(k,iEdge) - vertSum
+ enddo
+ enddo ! iEdge
+
+ call mpas_timer_stop("ocn_filter_btr_mode_tend_u")
+
+ end subroutine ocn_filter_btr_mode_tend_u!}}}
+
+!***********************************************************************
+!
! routine ocn_tendency_init
!
!> \brief Initializes flags used within tendency routines.
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-06-15 19:51:04 UTC (rev 1992)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-06-15 20:14:50 UTC (rev 1993)
@@ -81,21 +81,15 @@
type (state_type), target :: provis
type (state_type), pointer :: provis_ptr
- integer :: rk_step, iEdge, cell1, cell2
+ integer :: rk_step
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
- integer :: nCells, nEdges, nVertLevels, num_tracers
- real (kind=RKIND) :: coef
- real (kind=RKIND), dimension(:,:), pointer :: &
- u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell, ke_edge
+ integer :: nCells
+ real (kind=RKIND), dimension(:,:), pointer :: u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell, ke_edge
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- integer, dimension(:), pointer :: &
- maxLevelCell, maxLevelEdgeTop
- real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp
- real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
+ integer, dimension(:), pointer :: maxLevelCell
-
block => domain % blocklist
call mpas_allocate_state(block, provis, &
block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
@@ -171,7 +165,7 @@
! mrp 110718 filter btr mode out of u_tend
! still got h perturbations with just this alone. Try to set uBtr=0 after full u computation
if (config_rk_filter_btr_mode) then
- call filter_btr_mode_tend_u(block % tend, provis, block % diagnostics, block % mesh)
+ call ocn_filter_btr_mode_tend_u(block % tend, provis, block % mesh)
endif
call ocn_tend_scalar(block % tend, provis, block % diagnostics, block % mesh, dt)
@@ -276,15 +270,11 @@
h => block % state % time_levs(2) % state % h % array
h_edge => block % state % time_levs(2) % state % h_edge % array
ke_edge => block % state % time_levs(2) % state % ke_edge % array
- num_tracers = block % state % time_levs(2) % state % num_tracers
vertViscTopOfEdge => block % diagnostics % vertViscTopOfEdge % array
vertDiffTopOfCell => block % diagnostics % vertDiffTopOfCell % array
maxLevelCell => block % mesh % maxLevelCell % array
- maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
nCells = block % mesh % nCells
- nEdges = block % mesh % nEdges
- nVertLevels = block % mesh % nVertLevels
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
@@ -300,11 +290,11 @@
!
! Implicit vertical solve for momentum
!
- call ocn_vel_vmix_tend_implicit(block % mesh, dt, ke_edge, vertvisctopofedge, h, h_edge, u, err)
+ call ocn_vel_vmix_tend_implicit(block % mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
! mrp 110718 filter btr mode out of u
if (config_rk_filter_btr_mode) then
- call filter_btr_mode_u(block % state % time_levs(2) % state, block % mesh)
+ call ocn_filter_btr_mode_u(block % state % time_levs(2) % state, block % mesh)
!block % tend % h % array(:,:) = 0.0 ! I should not need this
endif
@@ -312,10 +302,26 @@
! Implicit vertical solve for tracers
!
- call ocn_tracer_vmix_tend_implicit(block % mesh, dt, vertdifftopofcell, h, tracers, err)
+ call ocn_tracer_vmix_tend_implicit(block % mesh, dt, vertDiffTopOfCell, h, tracers, err)
call mpas_timer_stop("RK4-implicit vert mix")
end if
+ block => block % next
+ end do
+
+ ! Update halo on u and tracers, which weres just updated for implicit vertical mixing. If not done,
+ ! this leads to lack of volume conservation. It is required because halo updates in RK4 are only
+ ! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
+ ! communicate the change due to implicit vertical mixing across the boundary.
+
+ if (config_implicit_vertical_mix) then
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
+ end if
+
+ block => domain % blocklist
+ do while (associated(block))
+
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
@@ -363,234 +369,6 @@
end subroutine ocn_time_integrator_rk4!}}}
- subroutine filter_btr_mode_tend_u(tend, s, d, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Filter and remove barotropic mode from the tendencies
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (diagnostics_type), intent(in) :: d
- type (mesh_type), intent(in) :: grid
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-! Some of these variables can be removed, but at a later time.
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j
-
- integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: vertSum, uhSum, hSum, sshEdge
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, dvEdge, dcEdge, areaCell, areaTriangle, &
- meshScalingDel2, meshScalingDel4
- real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
- tend_u, circulation, vorticity, ke, ke_edge, Vor_edge, &
- MontPot, wTop, divergence, vertViscTopOfEdge
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
- edgesOnEdge, edgesOnVertex
- real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
-
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-
- real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
-
- call mpas_timer_start("filter_btr_mode_tend_u")
-
- h => s % h % array
- u => s % u % array
- v => s % v % array
- wTop => s % wTop % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- ke_edge => s % ke_edge % array
- Vor_edge => s % Vor_edge % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
- vertViscTopOfEdge => d % vertViscTopOfEdge % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
-
- tend_u => tend % u % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nEdgesSolve = grid % nEdgesSolve
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- u_src => grid % u_src % array
-
- do iEdge=1,grid % nEdges
-
- uhSum = (h_edge(1,iEdge)) * tend_u(1,iEdge)
- hSum = h_edge(1,iEdge)
-
- do k=2,grid % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + h_edge(k,iEdge) * tend_u(k,iEdge)
- hSum = hSum + h_edge(k,iEdge)
- enddo
-
- vertSum = uhSum/hSum
-
- do k=1,grid % maxLevelEdgeTop % array(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) - vertSum
- enddo
-
- enddo ! iEdge
-
- call mpas_timer_stop("filter_btr_mode_tend_u")
-
- end subroutine filter_btr_mode_tend_u!}}}
-
- subroutine filter_btr_mode_u(s, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Filter and remove barotropic mode.
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-! Some of these variables can be removed, but at a later time.
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j
-
- integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: vertSum, uhSum, hSum, sshEdge
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, dvEdge, dcEdge, areaCell, areaTriangle, &
- meshScalingDel2, meshScalingDel4
- real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
- tend_u, circulation, vorticity, ke, ke_edge, Vor_edge, &
- MontPot, wTop, divergence, vertViscTopOfEdge
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
- edgesOnEdge, edgesOnVertex
- real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
-
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-
- real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
-
- call mpas_timer_start("filter_btr_mode_u")
-
- h => s % h % array
- u => s % u % array
- v => s % v % array
- wTop => s % wTop % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- ke_edge => s % ke_edge % array
- Vor_edge => s % Vor_edge % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nEdgesSolve = grid % nEdgesSolve
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
-
- u_src => grid % u_src % array
-
- do iEdge=1,grid % nEdges
-
- uhSum = (h_edge(1,iEdge)) * u(1,iEdge)
- hSum = h_edge(1,iEdge)
-
- do k=2,grid % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + h_edge(k,iEdge) * u(k,iEdge)
- hSum = hSum + h_edge(k,iEdge)
- enddo
-
- vertSum = uhSum/hSum
- do k=1,grid % maxLevelEdgeTop % array(iEdge)
- u(k,iEdge) = u(k,iEdge) - vertSum
- enddo
-
- enddo ! iEdge
-
- call mpas_timer_stop("filter_btr_mode_u")
-
- end subroutine filter_btr_mode_u!}}}
-
end module ocn_time_integration_rk4
! vim: foldmethod=marker
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-06-15 19:51:04 UTC (rev 1992)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-06-15 20:14:50 UTC (rev 1993)
@@ -895,133 +895,8 @@
end do
call mpas_timer_stop("se timestep", timer_main)
-
end subroutine ocn_time_integrator_split!}}}
- subroutine filter_btr_mode_tend_u(tend, s, d, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Filter and remove barotropic mode from the tendencies
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (diagnostics_type), intent(in) :: d
- type (mesh_type), intent(in) :: grid
-
- integer :: iEdge, k
-
- integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: vertSum, uhSum, hSum
- real (kind=RKIND), dimension(:,:), pointer :: &
- h_edge, h, u,tend_u
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: maxLevelEdgeTop
-
- call mpas_timer_start("filter_btr_mode_tend_u")
-
- h => s % h % array
- u => s % u % array
- h_edge => s % h_edge % array
-
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
-
- tend_u => tend % u % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
-
- do iEdge=1,nEdges
-
- ! hSum is initialized outside the loop because on land boundaries
- ! maxLevelEdgeTop=0, but I want to initialize hSum with a
- ! nonzero value to avoid a NaN.
- uhSum = h_edge(1,iEdge) * tend_u(1,iEdge)
- hSum = h_edge(1,iEdge)
-
- do k=2,maxLevelEdgeTop(iEdge)
- uhSum = uhSum + h_edge(k,iEdge) * tend_u(k,iEdge)
- hSum = hSum + h_edge(k,iEdge)
- enddo
-
- vertSum = uhSum/hSum
- do k=1,maxLevelEdgeTop(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) - vertSum
- enddo
- enddo ! iEdge
-
- call mpas_timer_stop("filter_btr_mode_tend_u")
-
- end subroutine filter_btr_mode_tend_u!}}}
-
- subroutine filter_btr_mode_u(s, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Filter and remove barotropic mode.
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
-
- integer :: iEdge, k
-
- integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: vertSum, uhSum, hSum
- real (kind=RKIND), dimension(:,:), pointer :: &
- h_edge, h, u
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: maxLevelEdgeTop
-
- call mpas_timer_start("filter_btr_mode_u")
-
- h => s % h % array
- u => s % u % array
- h_edge => s % h_edge % array
-
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
-
- do iEdge=1,nEdges
-
- ! hSum is initialized outside the loop because on land boundaries
- ! maxLevelEdgeTop=0, but I want to initialize hSum with a
- ! nonzero value to avoid a NaN.
- uhSum = h_edge(1,iEdge) * u(1,iEdge)
- hSum = h_edge(1,iEdge)
-
- do k=2,maxLevelEdgeTop(iEdge)
- uhSum = uhSum + h_edge(k,iEdge) * u(k,iEdge)
- hSum = hSum + h_edge(k,iEdge)
- enddo
-
- vertSum = uhSum/hSum
- do k=1,maxLevelEdgeTop(iEdge)
- u(k,iEdge) = u(k,iEdge) - vertSum
- enddo
- enddo ! iEdge
-
- call mpas_timer_stop("filter_btr_mode_u")
-
- end subroutine filter_btr_mode_u!}}}
-
end module ocn_time_integration_split
! vim: foldmethod=marker
</font>
</pre>