<p><b>dwj07@fsu.edu</b> 2012-07-23 15:09:33 -0600 (Mon, 23 Jul 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Merging trunk into branch to pick up recent modifications to ocean core.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_equation_of_state.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_equation_of_state.F 2012-07-23 20:57:17 UTC (rev 2045)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_equation_of_state.F 2012-07-23 21:09:33 UTC (rev 2046)
@@ -74,10 +74,13 @@
! Input: grid - grid metadata
! s - state: tracers
! k_displaced
- ! If k_displaced<=0, state % rho is returned with no displaced
- ! If k_displaced>0,the state % rhoDisplaced is returned, and is for
+ !
+ ! If k_displaced==0, state % rho is returned with no displacement
+ !
+ ! If k_displaced~=0,the state % rhoDisplaced is returned, and is for
! a parcel adiabatically displaced from its original level to level
- ! k_displaced. This does not effect the linear EOS.
+ ! k_displaced. When using the linear EOS, state % rhoDisplaced is
+ ! still filled, but depth (i.e. pressure) does not modify the output.
!
! Output: s - state: computed density
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -103,19 +106,19 @@
indexT = s % index_temperature
indexS = s % index_salinity
- if (linearEos) then
+ ! Choose to fill the array rho or rhoDisplaced
+ if (k_displaced == 0) then
rho => s % rho % array
+ else
+ rho => s % rhoDisplaced % array
+ endif
+ if (linearEos) then
+
call ocn_equation_of_state_linear_rho(grid, indexT, indexS, tracers, rho, err)
elseif (jmEos) then
- if(k_displaced == 0) then
- rho => s % rho % array
- else
- rho => s % rhoDisplaced % array
- endif
-
call ocn_equation_of_state_jm_rho(grid, k_displaced, displacement_type, indexT, indexS, tracers, rho, err)
endif
Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_mpas_core.F 2012-07-23 20:57:17 UTC (rev 2045)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_mpas_core.F 2012-07-23 21:09:33 UTC (rev 2046)
@@ -259,7 +259,7 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(1) % state, mesh)
call mpas_timer_stop("diagnostic solve", initDiagSolveTimer)
- ! Compute velocity transport, used in advection terms of h and tracer tendancy
+ ! Compute velocity transport, used in advection terms of h and tracer tendency
block % state % time_levs(1) % state % uTransport % array(:,:) &
= block % state % time_levs(1) % state % u % array(:,:) &
+ block % state % time_levs(1) % state % uBolusGM % array(:,:)
Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tendency.F 2012-07-23 20:57:17 UTC (rev 2045)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tendency.F 2012-07-23 21:09:33 UTC (rev 2046)
@@ -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
!--------------------------------------------------------------------
!
@@ -1034,9 +1036,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: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-07-23 20:57:17 UTC (rev 2045)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-07-23 21:09:33 UTC (rev 2046)
@@ -79,7 +79,7 @@
integer :: iCell, k, i, err
type (block_type), pointer :: block
- integer :: rk_step, iEdge, cell1, cell2
+ integer :: rk_step
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
@@ -161,7 +161,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, block % provis, block % diagnostics, block % mesh)
+ call ocn_filter_btr_mode_tend_u(block % tend, provis, block % mesh)
endif
call ocn_tend_scalar(block % tend, block % provis, block % diagnostics, block % mesh, dt)
@@ -212,7 +212,7 @@
call ocn_diagnostic_solve(dt, block % provis, block % mesh)
- ! Compute velocity transport, used in advection terms of h and tracer tendancy
+ ! Compute velocity transport, used in advection terms of h and tracer tendency
block % provis % uTransport % array(:,:) &
= block % provis % u % array(:,:) &
+ block % provis % uBolusGM % array(:,:)
@@ -255,54 +255,30 @@
! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
!
call mpas_timer_start("RK4-cleaup phase")
- block => domain % blocklist
- do while (associated(block))
- u => block % state % time_levs(2) % state % u % array
- tracers => block % state % time_levs(2) % state % tracers % array
- 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
+ if (config_implicit_vertical_mix) then
+ call mpas_timer_start("RK4-implicit vert mix")
+ block => domain % blocklist
+ do while(associated(block))
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ block => block % next
+ end do
- do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- tracers(:,k,iCell) = tracers(:,k,iCell) / h(k,iCell)
- end do
- end do
+ ! Update halo on u and tracers, which were 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.
+ call mpas_timer_start("RK4-implicit vert mix halos")
+ 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)
+ call mpas_timer_stop("RK4-implicit vert mix halos")
- if (config_implicit_vertical_mix) then
- call mpas_timer_start("RK4-implicit vert mix")
+ call mpas_timer_stop("RK4-implicit vert mix")
+ end if
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
+ block => domain % blocklist
+ do while (associated(block))
- !
- ! Implicit vertical solve for momentum
- !
- 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)
- !block % tend % h % array(:,:) = 0.0 ! I should not need this
- endif
-
- !
- ! Implicit vertical solve for tracers
- !
-
- call ocn_tracer_vmix_tend_implicit(block % mesh, dt, vertdifftopofcell, h, tracers, err)
- call mpas_timer_stop("RK4-implicit vert mix")
- end if
-
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
@@ -317,7 +293,7 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- ! Compute velocity transport, used in advection terms of h and tracer tendancy
+ ! Compute velocity transport, used in advection terms of h and tracer tendency
block % state % time_levs(2) % state % uTransport % array(:,:) &
= block % state % time_levs(2) % state % u % array(:,:) &
+ block % state % time_levs(2) % state % uBolusGM % array(:,:)
@@ -350,234 +326,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: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_split.F 2012-07-23 20:57:17 UTC (rev 2045)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_time_integration_split.F 2012-07-23 21:09:33 UTC (rev 2046)
@@ -293,7 +293,7 @@
! uTranport = uBcl + uBolus
! This is u used in advective terms for h and tracers
- ! in tendancy calls in stage 3.
+ ! in tendency calls in stage 3.
block % state % time_levs(2) % state % uTransport % array(k,iEdge) &
= block % mesh % edgeMask % array(k,iEdge) &
*( block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
@@ -639,7 +639,7 @@
! uTranport = uBtr + uBcl + uBolus + uCorrection
! This is u used in advective terms for h and tracers
- ! in tendancy calls in stage 3.
+ ! in tendency calls in stage 3.
block % state % time_levs(2) % state % uTransport % array(k,iEdge) &
= block % mesh % edgeMask % array(k,iEdge) &
*( block % state % time_levs(2) % state % uBtr % array( iEdge) &
@@ -825,37 +825,29 @@
! END large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- block => domain % blocklist
- do while (associated(block))
+ if (config_implicit_vertical_mix) then
+ call mpas_timer_start("se implicit vert mix")
+ block => domain % blocklist
+ do while(associated(block))
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ block => block % next
+ end do
+ ! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
+ ! this leads to lack of volume conservation. It is required because halo updates in stage 3 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.
+ call mpas_timer_start("se implicit vert mix halos")
+ 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)
+ call mpas_timer_stop("se implicit vert mix halos")
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Implicit vertical mixing, done after timestep is complete
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ call mpas_timer_stop("se implicit vert mix")
+ end if
- u => block % state % time_levs(2) % state % u % array
- tracers => block % state % time_levs(2) % state % tracers % array
- 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
+ block => domain % blocklist
+ do while (associated(block))
- if (config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
-
- ! Implicit vertical solve for momentum
- call ocn_vel_vmix_tend_implicit(block % mesh, dt, ke_edge, vertvisctopofedge, h, h_edge, u, err)
-
- ! Implicit vertical solve for tracers
- call ocn_tracer_vmix_tend_implicit(block % mesh, dt, vertdifftopofcell, h, tracers, err)
- end if
-
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
@@ -870,12 +862,18 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ ! Compute velocity transport, used in advection terms of h and tracer tendency
+ block % state % time_levs(2) % state % uTransport % array(:,:) &
+ = block % state % time_levs(2) % state % u % array(:,:) &
+ + block % state % time_levs(2) % state % uBolusGM % array(:,:)
+
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
- block % state % time_levs(2) % state % uReconstructX % array, &
- block % state % time_levs(2) % state % uReconstructY % array, &
- block % state % time_levs(2) % state % uReconstructZ % array, &
- block % state % time_levs(2) % state % uReconstructZonal % array, &
- block % state % time_levs(2) % state % uReconstructMeridional % array)
+ block % state % time_levs(2) % state % uReconstructX % array, &
+ block % state % time_levs(2) % state % uReconstructY % array, &
+ block % state % time_levs(2) % state % uReconstructZ % array, &
+ block % state % time_levs(2) % state % uReconstructZonal % array, &
+ block % state % time_levs(2) % state % uReconstructMeridional % array &
+ )
!TDR
call mpas_reconstruct(block % mesh, block % mesh % u_src % array, &
@@ -887,141 +885,15 @@
)
!TDR
-
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
-
block => block % next
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
Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-07-23 20:57:17 UTC (rev 2045)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-07-23 21:09:33 UTC (rev 2046)
@@ -77,6 +77,8 @@
real (kind=RKIND), parameter :: eps = 1.e-10
+ type (field2dReal), pointer :: high_order_horiz_flux_field
+
! Initialize pointers
dvEdge => grid % dvEdge % array
cellsOnEdge => grid % cellsOnEdge % array
@@ -100,6 +102,17 @@
nVertLevels = grid % nVertLevels
num_tracers = size(tracers,dim=1)
+ allocate(high_order_horiz_flux_field)
+ nullify(high_order_horiz_flux_field % next)
+ high_order_horiz_flux_field % block => grid % block
+ high_order_horiz_flux_field % sendList => grid % xEdge % sendList
+ high_order_horiz_flux_field % recvList => grid % xEdge % recvList
+ high_order_horiz_flux_field % copyList => grid % xEdge % copyList
+ high_order_horiz_flux_field % dimSizes(1) = nVertLevels
+ high_order_horiz_flux_field % dimSizes(2) = nEdges+1
+ allocate(high_order_horiz_flux_field % array(high_order_horiz_flux_field % dimSizes(1), high_order_horiz_flux_field % dimSizes(2)))
+ high_order_horiz_flux => high_order_horiz_flux_field % array
+
! allocate nCells arrays
allocate(tracer_new(nVertLevels, nCells+1))
@@ -112,7 +125,7 @@
allocate(flux_outgoing(nVertLevels, nCells+1))
! allocate nEdges arrays
- allocate(high_order_horiz_flux(nVertLevels, nEdges))
+! allocate(high_order_horiz_flux(nVertLevels, nEdges))
! allocate nVertLevels+1 and nCells arrays
allocate(high_order_vert_flux(nVertLevels+1, nCells))
@@ -192,6 +205,8 @@
end do ! i loop over nAdvCellsForEdge
end do ! iEdge loop
+ call mpas_dmpar_exch_halo_field(high_order_horiz_flux_field)
+
! low order upwind vertical flux (monotonic and diffused)
! Remove low order flux from the high order flux.
! Store left over high order flux in high_order_vert_flux array.
@@ -348,6 +363,7 @@
deallocate(flux_outgoing)
deallocate(high_order_horiz_flux)
deallocate(high_order_vert_flux)
+ deallocate(high_order_horiz_flux_field)
end subroutine mpas_ocn_tracer_advection_mono_tend!}}}
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_vmix.F 2012-07-23 20:57:17 UTC (rev 2045)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_vmix.F 2012-07-23 21:09:33 UTC (rev 2046)
@@ -47,7 +47,8 @@
ocn_tracer_vmix_tend_explicit, &
ocn_vel_vmix_tend_implicit, &
ocn_tracer_vmix_tend_implicit, &
- ocn_vmix_init
+ ocn_vmix_init, &
+ ocn_vmix_implicit
!--------------------------------------------------------------------
!
@@ -576,6 +577,61 @@
!***********************************************************************
!
+! routine ocn_vmix_implicit
+!
+!> \brief Driver for implicit vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine is a driver for handling implicit vertical mixing
+!> of both momentum and tracers for a block. It's intended to reduce
+!> redundant code.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_vmix_implicit(dt, grid, diagnostics, state, err)!{{{
+ real (kind=RKIND), intent(in) :: dt
+ type (mesh_type), intent(in) :: grid
+ type (diagnostics_type), intent(inout) :: diagnostics
+ type (state_type), intent(inout) :: state
+ integer, intent(out) :: err
+
+ 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
+
+ err = 0
+
+ u => state % u % array
+ tracers => state % tracers % array
+ h => state % h % array
+ h_edge => state % h_edge % array
+ ke_edge => state % ke_edge % array
+ vertViscTopOfEdge => diagnostics % vertViscTopOfEdge % array
+ vertDiffTopOfCell => diagnostics % vertDiffTopOfCell % array
+ maxLevelCell => grid % maxLevelCell % array
+
+ nCells = grid % nCells
+
+ call ocn_vmix_coefs(grid, state, diagnostics, err)
+
+ !
+ ! Implicit vertical solve for momentum
+ !
+ call ocn_vel_vmix_tend_implicit(grid, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
+
+ !
+ ! Implicit vertical solve for tracers
+ !
+
+ call ocn_tracer_vmix_tend_implicit(grid, dt, vertDiffTopOfCell, h, tracers, err)
+
+ end subroutine ocn_vmix_implicit!}}}
+
+!***********************************************************************
+!
! routine ocn_vmix_init
!
!> \brief Initializes ocean vertical mixing quantities
Modified: branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-07-23 20:57:17 UTC (rev 2045)
+++ branches/omp_blocks/multiple_blocks/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-07-23 21:09:33 UTC (rev 2046)
@@ -472,7 +472,7 @@
drhoTopOfCell = 0.0
do iCell=1,nCells
do k=2,maxLevelCell(iCell)
- drhoTopOfCell(k,iCell) = rho(k-1,iCell) - rhoDisplaced(k-1,iCell)
+ drhoTopOfCell(k,iCell) = rhoDisplaced(k-1,iCell) - rhoDisplaced(k,iCell)
end do
end do
</font>
</pre>