<p><b>lowrie@lanl.gov</b> 2013-03-13 11:04:50 -0600 (Wed, 13 Mar 2013)</p><p><br>
Merge with trunk.<br>
</p><hr noshade><pre><font color="gray">Index: branches/mpas_cdg_advection
===================================================================
--- branches/mpas_cdg_advection 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection 2013-03-13 17:04:50 UTC (rev 2600)
Property changes on: branches/mpas_cdg_advection
___________________________________________________________________
Modified: svn:mergeinfo
## -5,6 +5,7 ##
/branches/ocean_projects/ale_vert_coord_new:1387-1428
/branches/ocean_projects/cesm_coupling:2147-2344
/branches/ocean_projects/diagnostics_revision:2439-2462
+/branches/ocean_projects/explicit_vmix_removal:2486-2490
/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
## -15,7 +16,9 ##
/branches/ocean_projects/namelist_cleanup:2319-2414
/branches/ocean_projects/option3_b4b_test:2201-2231
/branches/ocean_projects/partial_bottom_cells:2172-2226
+/branches/ocean_projects/remove_sw_test_cases:2539-2540
/branches/ocean_projects/restart_reproducibility:2239-2272
+/branches/ocean_projects/sea_level_pressure:2488-2528
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
## -28,4 +31,4 ##
/branches/omp_blocks/multiple_blocks:1803-2084
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
-/trunk/mpas:2390-2472
+/trunk/mpas:2390-2599
\ No newline at end of property
Modified: branches/mpas_cdg_advection/namelist.input.ocean
===================================================================
--- branches/mpas_cdg_advection/namelist.input.ocean 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/namelist.input.ocean 2013-03-13 17:04:50 UTC (rev 2600)
@@ -69,7 +69,6 @@
config_Rayleigh_damping_coeff = 0.0
/
&vmix
- config_implicit_vertical_mix = .true.
config_convective_visc = 1.0
config_convective_diff = 1.0
/
@@ -140,9 +139,6 @@
config_btr_gam3_uWt2 = 1.0
config_btr_solve_SSH2 = .false.
/
-&sw_model
- config_test_case = 0
-/
&debug
config_check_zlevel_consistency = .false.
config_filter_btr_mode = .false.
Modified: branches/mpas_cdg_advection/src/core_hyd_atmos/Registry
===================================================================
--- branches/mpas_cdg_advection/src/core_hyd_atmos/Registry 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_hyd_atmos/Registry 2013-03-13 17:04:50 UTC (rev 2600)
@@ -96,7 +96,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 o edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 o localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 o cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 o cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
Modified: branches/mpas_cdg_advection/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/mpas_cdg_advection/src/core_init_nhyd_atmos/Registry 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_init_nhyd_atmos/Registry 2013-03-13 17:04:50 UTC (rev 2600)
@@ -108,7 +108,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 io edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 io localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 io cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 io cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 io cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 io verticesOnCell mesh - -
Modified: branches/mpas_cdg_advection/src/core_nhyd_atmos/Registry
===================================================================
--- branches/mpas_cdg_advection/src/core_nhyd_atmos/Registry 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_nhyd_atmos/Registry 2013-03-13 17:04:50 UTC (rev 2600)
@@ -113,7 +113,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 iro edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 iro localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 iro cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 iro cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
Index: branches/mpas_cdg_advection/src/core_ocean
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean 2013-03-13 17:04:50 UTC (rev 2600)
Property changes on: branches/mpas_cdg_advection/src/core_ocean
___________________________________________________________________
Modified: svn:mergeinfo
## -5,6 +5,7 ##
/branches/ocean_projects/ale_vert_coord_new/src/core_ocean:1387-1428
/branches/ocean_projects/cesm_coupling/src/core_ocean:2147-2344
/branches/ocean_projects/diagnostics_revision/src/core_ocean:2439-2462
+/branches/ocean_projects/explicit_vmix_removal/src/core_ocean:2486-2490
/branches/ocean_projects/gmvar/src/core_ocean:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_error/src/core_ocean:1847-1887
/branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean:754-986
## -15,7 +16,9 ##
/branches/ocean_projects/namelist_cleanup/src/core_ocean:2319-2414
/branches/ocean_projects/option3_b4b_test/src/core_ocean:2201-2231
/branches/ocean_projects/partial_bottom_cells/src/core_ocean:2172-2226
+/branches/ocean_projects/remove_sw_test_cases/src/core_ocean:2539-2540
/branches/ocean_projects/restart_reproducibility/src/core_ocean:2239-2272
+/branches/ocean_projects/sea_level_pressure/src/core_ocean:2488-2528
/branches/ocean_projects/split_explicit_mrp/src/core_ocean:1134-1138
/branches/ocean_projects/split_explicit_timestepping/src/core_ocean:1044-1097
/branches/ocean_projects/vert_adv_mrp/src/core_ocean:704-745
## -30,4 +33,4 ##
/branches/omp_blocks/openmp_test/src/core_ocean_elements:2161-2201
/branches/source_renaming/src/core_ocean:1082-1113
/branches/time_manager/src/core_ocean:924-962
-/trunk/mpas/src/core_ocean:2390-2472
+/trunk/mpas/src/core_ocean:2390-2599
\ No newline at end of property
Modified: branches/mpas_cdg_advection/src/core_ocean/Makefile
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/Makefile 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/Makefile 2013-03-13 17:04:50 UTC (rev 2600)
@@ -1,7 +1,6 @@
.SUFFIXES: .F .o
OBJS = mpas_ocn_mpas_core.o \
- mpas_ocn_test_cases.o \
mpas_ocn_advection.o \
mpas_ocn_thick_hadv.o \
mpas_ocn_thick_vadv.o \
@@ -14,7 +13,6 @@
mpas_ocn_vel_hmix_del4.o \
mpas_ocn_vel_forcing.o \
mpas_ocn_vel_forcing_windstress.o \
- mpas_ocn_vel_forcing_bottomdrag.o \
mpas_ocn_vel_forcing_rayleigh.o \
mpas_ocn_vel_pressure_grad.o \
mpas_ocn_tracer_vadv.o \
@@ -64,8 +62,6 @@
core_hyd: $(OBJS)
ar -ru libdycore.a $(OBJS)
-mpas_ocn_test_cases.o:
-
mpas_ocn_advection.o:
mpas_ocn_time_integration.o: mpas_ocn_time_integration_rk4.o mpas_ocn_time_integration_split.o
@@ -100,12 +96,10 @@
mpas_ocn_vel_hmix_del4.o:
-mpas_ocn_vel_forcing.o: mpas_ocn_vel_forcing_windstress.o mpas_ocn_vel_forcing_bottomdrag.o mpas_ocn_vel_forcing_rayleigh.o
+mpas_ocn_vel_forcing.o: mpas_ocn_vel_forcing_windstress.o mpas_ocn_vel_forcing_rayleigh.o
mpas_ocn_vel_forcing_windstress.o:
-mpas_ocn_vel_forcing_bottomdrag.o:
-
mpas_ocn_vel_forcing_rayleigh.o:
mpas_ocn_vel_coriolis.o:
@@ -176,8 +170,7 @@
mpas_ocn_monthly_forcing.o:
-mpas_ocn_mpas_core.o: mpas_ocn_test_cases.o \
- mpas_ocn_advection.o \
+mpas_ocn_mpas_core.o: mpas_ocn_advection.o \
mpas_ocn_thick_hadv.o \
mpas_ocn_gm.o \
mpas_ocn_thick_vadv.o \
@@ -189,7 +182,6 @@
mpas_ocn_vel_hmix_del4.o \
mpas_ocn_vel_forcing.o \
mpas_ocn_vel_forcing_windstress.o \
- mpas_ocn_vel_forcing_bottomdrag.o \
mpas_ocn_vel_pressure_grad.o \
mpas_ocn_tracer_vadv.o \
mpas_ocn_tracer_vadv_spline.o \
Modified: branches/mpas_cdg_advection/src/core_ocean/Registry
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/Registry 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/Registry 2013-03-13 17:04:50 UTC (rev 2600)
@@ -60,7 +60,6 @@
namelist logical Rayleigh_damping config_Rayleigh_friction .false.
namelist real Rayleigh_damping config_Rayleigh_damping_coeff 0.0
-namelist logical vmix config_implicit_vertical_mix .true.
namelist real vmix config_convective_visc 1.0
namelist real vmix config_convective_diff 1.0
@@ -121,8 +120,6 @@
namelist real split_explicit_ts config_btr_gam3_uWt2 1.0
namelist logical split_explicit_ts config_btr_solve_SSH2 .false.
-namelist integer sw_model config_test_case 0
-
namelist logical debug config_check_zlevel_consistency .false.
namelist logical debug config_filter_btr_mode .false.
namelist logical debug config_prescribe_velocity .false.
@@ -210,7 +207,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 o edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 o localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 o cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 o cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
@@ -380,3 +377,6 @@
var persistent integer edgeSignOnCell ( maxEdges nCells ) 0 - edgeSignOnCell mesh - -
var persistent integer edgeSignOnVertex ( maxEdges nVertices ) 0 - edgeSignOnVertex mesh - -
var persistent integer kiteIndexOnCell ( maxEdges nCells ) 0 - kiteIndexOnCell mesh - -
+
+% Sea surface pressure, for coupling
+var persistent real seaSurfacePressure ( nCells Time ) 0 ir seaSurfacePressure mesh - -
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_diagnostics.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_diagnostics.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_diagnostics.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -97,7 +97,7 @@
real (kind=RKIND), dimension(:), allocatable:: pTop, div_hu
real (kind=RKIND), dimension(:), pointer :: &
- bottomDepth, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, ssh
+ bottomDepth, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, ssh, seaSurfacePressure
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -158,6 +158,8 @@
maxLevelVertexBot => grid % maxLevelVertexBot % array
kiteIndexOnCell => grid % kiteIndexOnCell % array
verticesOnCell => grid % verticesOnCell % array
+
+ seaSurfacePressure => grid % seaSurfacePressure % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -422,9 +424,10 @@
else
do iCell=1,nCells
- ! pressure for generalized coordinates
- ! assume atmospheric pressure at the surface is zero for now.
- pressure(1,iCell) = rho(1,iCell)*gravity &
+ ! Pressure for generalized coordinates.
+ ! Pressure at top surface may be due to atmospheric pressure
+ ! or an ice-shelf depression.
+ pressure(1,iCell) = seaSurfacePressure(iCell) + rho(1,iCell)*gravity &
* 0.5*h(1,iCell)
do k=2,maxLevelCell(iCell)
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_mpas_core.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_mpas_core.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -7,7 +7,6 @@
use mpas_timer
use ocn_global_diagnostics
- use ocn_test_cases
use ocn_time_integration
use ocn_tendency
use ocn_diagnostics
@@ -118,15 +117,13 @@
call mpas_dmpar_abort(dminfo)
endif
- if (.not. config_do_restart) call setup_sw_test_case(domain)
-
if (config_vert_coord_movement.ne.'isopycnal') call ocn_init_vert_coord(domain)
call ocn_compute_max_level(domain)
if (.not.config_do_restart) call ocn_init_split_timestep(domain)
- write (0,'(a,a10)') ' Vertical coordinate movement is: ',config_vert_coord_movement
+ write (0,'(a,a)') ' Vertical coordinate movement is: ',trim(config_vert_coord_movement)
if (config_vert_coord_movement.ne.'isopycnal'.and. &
config_vert_coord_movement.ne.'fixed'.and. &
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tendency.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tendency.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -41,8 +41,8 @@
save
type (timer_node), pointer :: thickHadvTimer, thickVadvTimer
- type (timer_node), pointer :: velCorTimer, velVadvTimer, velPgradTimer, velHmixTimer, velForceTimer, velExpVmixTimer
- type (timer_node), pointer :: tracerHadvTimer, tracerVadvTimer, tracerHmixTimer, tracerExpVmixTimer, tracerRestoringTimer
+ type (timer_node), pointer :: velCorTimer, velVadvTimer, velPgradTimer, velHmixTimer, velForceTimer
+ type (timer_node), pointer :: tracerHadvTimer, tracerVadvTimer, tracerHmixTimer, tracerRestoringTimer
!--------------------------------------------------------------------
!
@@ -240,11 +240,6 @@
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
- if (.not.config_implicit_vertical_mix) then
- call mpas_timer_start("explicit vmix", .false., velExpVmixTimer)
- call ocn_vel_vmix_tend_explicit(grid, u, h_edge, vertvisctopofedge, tend_u, err)
- call mpas_timer_stop("explicit vmix", velExpVmixTimer)
- endif
call mpas_timer_stop("ocn_tend_u")
end subroutine ocn_tend_u!}}}
@@ -333,17 +328,7 @@
! maxval(tracers(3,1,1:nCells))
! mrp 110516 printing end
- !
- ! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
- !
- if (.not.config_implicit_vertical_mix) then
- call mpas_timer_start("explicit vmix", .false., tracerExpVmixTimer)
- call ocn_tracer_vmix_tend_explicit(grid, h, vertdifftopofcell, tracers, tend_tr, err)
-
- call mpas_timer_stop("explicit vmix", tracerExpVmixTimer)
- endif
-
! mrp 110516 printing
!print *, 'tend_tr 2',minval(tend_tr(3,1,1:nCells)),&
! maxval(tend_tr(3,1,1:nCells))
Deleted: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_test_cases.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_test_cases.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_test_cases.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -1,526 +0,0 @@
- module ocn_test_cases
-
- use mpas_grid_types
- use mpas_configure
- use mpas_constants
-
-
- contains
-
-
- subroutine setup_sw_test_case(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Configure grid metadata and model state for the shallow water test case
- ! specified in the namelist
- !
- ! Output: block - a subset (not necessarily proper) of the model domain to be
- ! initialized
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: i, iCell, iEdge, iVtx, iLevel
- type (block_type), pointer :: block_ptr
- type (dm_info) :: dminfo
-
- if (config_test_case == 0) then
- write(0,*) 'Using initial conditions supplied in input file'
-
- else if (config_test_case == 1) then
- write(0,*) ' Setting up shallow water test case 1:'
- write(0,*) ' Advection of Cosine Bell over the Pole'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 2) then
- write(0,*) ' Setup shallow water test case 2: '// &
- 'Global Steady State Nonlinear Zonal Geostrophic Flow'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_2(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 5) then
- write(0,*) ' Setup shallow water test case 5:'// &
- ' Zonal Flow over an Isolated Mountain'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_5(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 6) then
- write(0,*) ' Set up shallow water test case 6:'
- write(0,*) ' Rossby-Haurwitz Wave'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_6(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- block_ptr => block_ptr % next
- end do
-
- else
- write(0,*) 'Abort: config_test_case=',config_test_case
- write(0,*) 'Only test case 1, 2, 5, and 6 ', &
- 'are currently supported. '
- call mpas_dmpar_abort(dminfo)
- end if
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
-
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, &
- block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- end subroutine setup_sw_test_case
-
-
- subroutine sw_test_case_1(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 1: Advection of Cosine Bell over the Pole
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
- real (kind=RKIND), parameter :: h0 = 1000.0
- real (kind=RKIND), parameter :: theta_c = 0.0
- real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: alpha = pii/4.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: r, u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
- !
- ! Initialize wind field
- !
- allocate(psiVertex(grid % nVertices))
- do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * u0 * ( &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
- cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
- )
- end do
- do iEdge=1,grid % nEdges
- state % u % array(1,iEdge) = -1.0 * ( &
- psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
- psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
- ) / grid%dvEdge%array(iEdge)
- end do
- deallocate(psiVertex)
-
- !
- ! Initialize cosine bell at (theta_c, lambda_c)
- !
- do iCell=1,grid % nCells
- r = sphere_distance(theta_c, lambda_c, grid % latCell % array(iCell), grid % lonCell % array(iCell), a)
- if (r < a/3.0) then
- state % h % array(1,iCell) = (h0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
- else
- state % h % array(1,iCell) = 0.0
- end if
- end do
-
- end subroutine sw_test_case_1
-
-
- subroutine sw_test_case_2(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 2: Global Steady State Nonlinear Zonal
- ! Geostrophic Flow
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
- real (kind=RKIND), parameter :: gh0 = 29400.0
- real (kind=RKIND), parameter :: alpha = 0.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
-
- !
- ! Initialize wind field
- !
- allocate(psiVertex(grid % nVertices))
- do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * u0 * ( &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
- cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
- )
- end do
- do iEdge=1,grid % nEdges
- state % u % array(1,iEdge) = -1.0 * ( &
- psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
- psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
- ) / grid%dvEdge%array(iEdge)
- end do
- deallocate(psiVertex)
-
- !
- ! Generate rotated Coriolis field
- !
- do iEdge=1,grid % nEdges
- grid % fEdge % array(iEdge) = 2.0 * omega * &
- ( -cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(alpha) &
- )
- end do
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 2.0 * omega * &
- (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) &
- )
- end do
-
- !
- ! Initialize height field (actually, fluid thickness field)
- !
- do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
- (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
- sin(grid%latCell%array(iCell)) * cos(alpha) &
- )**2.0 &
- ) / &
- gravity
- end do
-
- end subroutine sw_test_case_2
-
-
- subroutine sw_test_case_5(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 5: Zonal Flow over an Isolated Mountain
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 20.
- real (kind=RKIND), parameter :: gh0 = 5960.0*gravity
-! real (kind=RKIND), parameter :: hs0 = 2000. original
- real (kind=RKIND), parameter :: hs0 = 250. !mrp 100204
- real (kind=RKIND), parameter :: theta_c = pii/6.0
- real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: rr = pii/9.0
- real (kind=RKIND), parameter :: alpha = 0.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: r, u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
- !
- ! Initialize wind field
- !
- allocate(psiVertex(grid % nVertices))
- do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * u0 * ( &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
- cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
- )
- end do
- do iEdge=1,grid % nEdges
- state % u % array(1,iEdge) = -1.0 * ( &
- psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
- psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
- ) / grid%dvEdge%array(iEdge)
- end do
- deallocate(psiVertex)
-
- !
- ! Generate rotated Coriolis field
- !
- do iEdge=1,grid % nEdges
- grid % fEdge % array(iEdge) = 2.0 * omega * &
- (-cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(alpha) &
- )
- end do
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 2.0 * omega * &
- (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) &
- )
- end do
-
- !
- ! Initialize mountain
- !
- do iCell=1,grid % nCells
- if (grid % lonCell % array(iCell) < 0.0) grid % lonCell % array(iCell) = grid % lonCell % array(iCell) + 2.0 * pii
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- grid % bottomDepth % array(iCell) = hs0 * (1.0 - r/rr)
- end do
-! output about mountain
-print *, 'bottomDepth',minval(grid % bottomDepth % array),sum(grid % bottomDepth % array)/grid % nCells, maxval(grid % bottomDepth % array)
-
- !
- ! Initialize tracer fields
- !
- do iCell=1,grid % nCells
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- state % tracers % array(1,1,iCell) = 1.0 - r/rr
- end do
- do iCell=1,grid % nCells
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + &
- (grid % latCell % array(iCell) - theta_c - pii/6.0)**2.0 &
- ) &
- )
- state % tracers % array(2,1,iCell) = 1.0 - r/rr
- end do
-
- !
- ! Initialize height field (actually, fluid thickness field)
- !
- do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
- (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
- sin(grid%latCell%array(iCell)) * cos(alpha) &
- )**2.0 &
- ) / &
- gravity
- state % h % array(1,iCell) = state % h % array(1,iCell) - grid % bottomDepth % array(iCell)
- end do
-
- end subroutine sw_test_case_5
-
-
- subroutine sw_test_case_6(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 6: Rossby-Haurwitz Wave
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: h0 = 8000.0
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
- !
- ! Initialize wind field
- !
- allocate(psiVertex(grid % nVertices))
- do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * a * w * sin(grid%latVertex%array(iVtx)) + &
- a *a * K * (cos(grid%latVertex%array(iVtx))**R) * &
- sin(grid%latVertex%array(iVtx)) * cos(R * grid%lonVertex%array(iVtx))
- end do
- do iEdge=1,grid % nEdges
- state % u % array(1,iEdge) = -1.0 * ( &
- psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
- psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
- ) / grid%dvEdge%array(iEdge)
- end do
- deallocate(psiVertex)
-
- !
- ! Initialize height field (actually, fluid thickness field)
- !
- do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gravity * h0 + a*a*aa(grid%latCell%array(iCell)) + &
- a*a*bb(grid%latCell%array(iCell)) * cos(R*grid%lonCell%array(iCell)) + &
- a*a*cc(grid%latCell%array(iCell)) * cos(2.0*R*grid%lonCell%array(iCell)) &
- ) / gravity
- end do
-
- end subroutine sw_test_case_6
-
-
- real function sphere_distance(lat1, lon1, lat2, lon2, radius)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
- ! sphere with given radius.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
- real (kind=RKIND) :: arg1
-
- arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
- cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
- sphere_distance = 2.*radius*asin(arg1)
-
- end function sphere_distance
-
-
- real function aa(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! A, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- real (kind=RKIND), intent(in) :: theta
-
- aa = 0.5 * w * (2.0 * omega + w) * cos(theta)**2.0 + &
- 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*R**2*cos(theta)**(-2.0))
-
- end function aa
-
-
- real function bb(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! B, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- real (kind=RKIND), intent(in) :: theta
-
- bb = (2.0*(omega + w)*K / ((R+1.0)*(R+2.0))) * cos(theta)**R * ((R**2.0 + 2.0*R + 2.0) - ((R+1.0)*cos(theta))**2.0)
-
- end function bb
-
-
- real function cc(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! C, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- real (kind=RKIND), intent(in) :: theta
-
- cc = 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 - R - 2.0)
-
- end function cc
-
-end module ocn_test_cases
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -149,11 +149,6 @@
call mpas_timer_start("RK4-tendency computations")
block => domain % blocklist
do while (associated(block))
-
- if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, err)
- end if
-
! advection of u uses u, while advection of h and tracers use uTransport.
call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
block % provis % u % array, block % provis % wTop % array, err)
@@ -201,10 +196,6 @@
end do
end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- end if
-
if (config_prescribe_velocity) then
block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
@@ -271,42 +262,35 @@
block => block % next
end do
+ call mpas_timer_start("RK4-implicit vert mix")
+ block => domain % blocklist
+ do while(associated(block))
- if (config_implicit_vertical_mix) then
- call mpas_timer_start("RK4-implicit vert mix")
- block => domain % blocklist
- do while(associated(block))
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- ! Call ocean diagnostic solve in preparation for vertical mixing. Note
- ! it is called again after vertical mixing, because u and tracers change.
- ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
- ! be computed. For kpp, more variables may be needed. Either way, this
- ! could be made more efficient by only computing what is needed for the
- ! implicit vmix routine that follows. mrp 121023.
- call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ block => block % next
+ end do
- 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 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")
- ! 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")
+ call mpas_timer_stop("RK4-implicit vert mix")
- call mpas_timer_stop("RK4-implicit vert mix")
- 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
-
if (config_prescribe_velocity) then
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -188,10 +188,6 @@
stage1_tend_time = min(split_explicit_step,2)
- if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, err)
- end if
-
! compute wTop. Use u (rather than uTransport) for momentum advection.
! Use the most recent time level available.
call ocn_wtop(block % mesh, block % state % time_levs(stage1_tend_time) % state % h % array, &
@@ -940,43 +936,37 @@
! END large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- if (config_implicit_vertical_mix) then
- call mpas_timer_start("se implicit vert mix")
- block => domain % blocklist
- do while(associated(block))
+ call mpas_timer_start("se implicit vert mix")
+ block => domain % blocklist
+ do while(associated(block))
- ! Call ocean diagnostic solve in preparation for vertical mixing. Note
- ! it is called again after vertical mixing, because u and tracers change.
- ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
- ! be computed. For kpp, more variables may be needed. Either way, this
- ! could be made more efficient by only computing what is needed for the
- ! implicit vmix routine that follows. mrp 121023.
- call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
- block => block % next
- end do
+ 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")
+ ! 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")
- call mpas_timer_stop("se implicit vert mix")
- end if
+ call mpas_timer_stop("se implicit vert mix")
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
-
if (config_prescribe_velocity) then
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -18,7 +18,6 @@
use mpas_configure
use ocn_vel_forcing_windstress
- use ocn_vel_forcing_bottomdrag
use ocn_vel_forcing_rayleigh
implicit none
@@ -115,7 +114,7 @@
!
!-----------------------------------------------------------------
- integer :: err1, err2, err3
+ integer :: err1, err2
!-----------------------------------------------------------------
!
@@ -126,11 +125,9 @@
!-----------------------------------------------------------------
call ocn_vel_forcing_windstress_tend(grid, u_src, h_edge, tend, err1)
- call ocn_vel_forcing_bottomdrag_tend(grid, u, ke_edge, h_edge, tend, err2)
- call ocn_vel_forcing_rayleigh_tend(grid, u, tend, err3)
+ call ocn_vel_forcing_rayleigh_tend(grid, u, tend, err2)
err = ior(err1, err2)
- err = ior(err, err3)
!--------------------------------------------------------------------
@@ -164,14 +161,12 @@
integer, intent(out) :: err !< Output: error flag
- integer :: err1, err2, err3
+ integer :: err1, err2
call ocn_vel_forcing_windstress_init(err1)
- call ocn_vel_forcing_bottomdrag_init(err2)
- call ocn_vel_forcing_rayleigh_init(err3)
+ call ocn_vel_forcing_rayleigh_init(err2)
err = ior(err1, err2)
- err = ior(err, err3)
!--------------------------------------------------------------------
Deleted: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -1,193 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_vel_forcing_bottomdrag
-!
-!> \brief MPAS ocean bottom drag
-!> \author Doug Jacobsen
-!> \date 16 September 2011
-!> \version SVN:$Id:$
-!> \details
-!> This module contains the routine for computing
-!> tendencies from bottom drag.
-!
-!-----------------------------------------------------------------------
-
-module ocn_vel_forcing_bottomdrag
-
- use mpas_grid_types
- use mpas_configure
-
- implicit none
- private
- save
-
- !--------------------------------------------------------------------
- !
- ! Public parameters
- !
- !--------------------------------------------------------------------
-
- !--------------------------------------------------------------------
- !
- ! Public member functions
- !
- !--------------------------------------------------------------------
-
- public :: ocn_vel_forcing_bottomdrag_tend, &
- ocn_vel_forcing_bottomdrag_init
-
- !--------------------------------------------------------------------
- !
- ! Private module variables
- !
- !--------------------------------------------------------------------
-
- logical :: bottomDragOn
- real (kind=RKIND) :: bottomDragCoef
-
-
-!***********************************************************************
-
-contains
-
-!***********************************************************************
-!
-! routine ocn_vel_forcing_bottomdrag_tend
-!
-!> \brief Computes tendency term from bottom drag
-!> \author Doug Jacobsen
-!> \date 15 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the bottom drag tendency for momentum
-!> based on current state.
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_vel_forcing_bottomdrag_tend(grid, u, ke_edge, h_edge, tend, err)!{{{
-
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- u !< Input: velocity
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- ke_edge !< Input: kinetic energy at edge
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- h_edge !< Input: thickness at edge
-
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
-
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:), intent(inout) :: &
- tend !< Input/Output: velocity tendency
-
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
-
- integer :: iEdge, nEdgesSolve, k
- integer, dimension(:), pointer :: maxLevelEdgeTop
- integer, dimension(:,:), pointer :: edgeMask
-
- !-----------------------------------------------------------------
- !
- ! call relevant routines for computing tendencies
- ! note that the user can choose multiple options and the
- ! tendencies will be added together
- !
- !-----------------------------------------------------------------
-
- err = 0
-
- if(.not.bottomDragOn) return
-
- nEdgesSolve = grid % nEdgesSolve
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- edgeMask => grid % edgeMask % array
-
- do iEdge=1,grid % nEdgesSolve
-
- k = max(maxLevelEdgeTop(iEdge), 1)
-
- ! bottom drag is the same as POP:
- ! -c |u| u where c is unitless and 1.0e-3.
- ! see POP Reference guide, section 3.4.4.
-
- tend(k,iEdge) = tend(k,iEdge)-edgeMask(k,iEdge)*(bottomDragCoef*u(k,iEdge)*sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge))
-
- enddo
-
-
-
- !--------------------------------------------------------------------
-
- end subroutine ocn_vel_forcing_bottomdrag_tend!}}}
-
-!***********************************************************************
-!
-! routine ocn_vel_forcing_bottomdrag_init
-!
-!> \brief Initializes ocean bottom drag
-!> \author Doug Jacobsen
-!> \date 16 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine initializes quantities related to bottom drag
-!> in the ocean.
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_vel_forcing_bottomdrag_init(err)!{{{
-
- !--------------------------------------------------------------------
-
- !-----------------------------------------------------------------
- !
- ! call individual init routines for each parameterization
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
-
- err = 0
-
- bottomDragOn = .false.
-
- if (.not.config_implicit_vertical_mix) then
- bottomDragOn = .true.
- bottomDragCoef = config_bottom_drag_coeff
- endif
-
- !--------------------------------------------------------------------
-
- end subroutine ocn_vel_forcing_bottomdrag_init!}}}
-
-!***********************************************************************
-
-end module ocn_vel_forcing_bottomdrag
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-! vim: foldmethod=marker
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -44,8 +44,6 @@
tridiagonal_solve_mult
public :: ocn_vmix_coefs, &
- ocn_vel_vmix_tend_explicit, &
- ocn_tracer_vmix_tend_explicit, &
ocn_vel_vmix_tend_implicit, &
ocn_tracer_vmix_tend_implicit, &
ocn_vmix_init, &
@@ -58,7 +56,6 @@
!--------------------------------------------------------------------
logical :: velVmixOn, tracerVmixOn
- logical :: explicitOn, implicitOn
!***********************************************************************
@@ -140,100 +137,6 @@
!***********************************************************************
!
-! routine ocn_vel_vmix_tendExplict
-!
-!> \brief Computes tendencies for explict momentum vertical mixing
-!> \author Doug Jacobsen
-!> \date 19 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the tendencies for explicit vertical mixing for momentum
-!> using computed coefficients.
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_vel_vmix_tend_explicit(grid, u, h_edge, vertViscTopOfEdge, tend, err)!{{{
-
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
-
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- u !< Input: velocity
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- h_edge !< Input: thickness at edge
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- vertViscTopOfEdge !< Input: vertical mixing coefficients
-
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:), intent(inout) :: &
- tend !< Input/Output: tendency information
-
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
-
- integer :: iEdge, nEdgesSolve, k, nVertLevels
-
- integer, dimension(:), pointer :: maxLevelEdgeTop
-
- real (kind=RKIND), dimension(:), allocatable :: fluxVertTop
-
- err = 0
-
- if(.not.velVmixOn) return
- if(implicitOn) return
-
- nEdgessolve = grid % nEdgesSolve
- nVertLevels = grid % nVertLevels
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
-
- allocate(fluxVertTop(nVertLevels+1))
- fluxVertTop(1) = 0.0
- do iEdge=1,nEdgesSolve
- do k=2,maxLevelEdgeTop(iEdge)
- fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &
- * ( u(k-1,iEdge) - u(k,iEdge) ) &
- * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
- enddo
- fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
-
- do k=1,maxLevelEdgeTop(iEdge)
- tend(k,iEdge) = tend(k,iEdge) &
- + (fluxVertTop(k) - fluxVertTop(k+1)) &
- / h_edge(k,iEdge)
- enddo
-
- end do
- deallocate(fluxVertTop)
- !--------------------------------------------------------------------
-
- end subroutine ocn_vel_vmix_tend_explicit!}}}
-
-!***********************************************************************
-!
! routine ocn_vel_vmix_tend_implicit
!
!> \brief Computes tendencies for implicit momentum vertical mixing
@@ -306,7 +209,6 @@
err = 0
if(.not.velVmixOn) return
- if(explicitOn) return
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
@@ -369,111 +271,6 @@
!***********************************************************************
!
-! routine ocn_tracer_vmix_tendExplict
-!
-!> \brief Computes tendencies for explict tracer vertical mixing
-!> \author Doug Jacobsen
-!> \date 19 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the tendencies for explicit vertical mixing for
-!> tracers using computed coefficients.
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_tracer_vmix_tend_explicit(grid, h, vertDiffTopOfCell, tracers, tend, err)!{{{
-
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
-
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- h !< Input: thickness at cell center
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- vertDiffTopOfCell !< Input: vertical mixing coefficients
-
- real (kind=RKIND), dimension(:,:,:), intent(in) :: &
- tracers !< Input: tracers
-
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
- tend !< Input/Output: tendency information
-
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
-
- integer :: iCell, nCellsSolve, k, iTracer, num_tracers, nVertLevels
-
- integer, dimension(:), pointer :: maxLevelCell
-
- real (kind=RKIND), dimension(:,:), allocatable :: fluxVertTop
-
- err = 0
-
- if(.not.tracerVmixOn) return
- if(implicitOn) return
-
- nCellsSolve = grid % nCellsSolve
- nVertLevels = grid % nVertLevels
- num_tracers = size(tracers, dim=1)
-
- maxLevelCell => grid % maxLevelCell % array
-
- allocate(fluxVertTop(num_tracers,nVertLevels+1))
- fluxVertTop(:,1) = 0.0
- do iCell=1,nCellsSolve
-
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! compute \kappa_v d\phi/dz
- fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
- * ( tracers(iTracer,k-1,iCell) &
- - tracers(iTracer,k ,iCell) ) &
- * 2 / (h(k-1,iCell) + h(k,iCell))
-
- enddo
- enddo
- fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
-
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! This is h d/dz( fluxVertTop) but h and dz cancel, so
- ! reduces to delta( fluxVertTop)
- tend(iTracer,k,iCell) = tend(iTracer,k,iCell) &
- + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
- enddo
- enddo
-
- enddo ! iCell loop
- deallocate(fluxVertTop)
- !--------------------------------------------------------------------
-
- end subroutine ocn_tracer_vmix_tend_explicit!}}}
-
-!***********************************************************************
-!
! routine ocn_tracer_vmix_tend_implicit
!
!> \brief Computes tendencies for implicit tracer vertical mixing
@@ -539,7 +336,6 @@
err = 0
if(.not.tracerVmixOn) return
- if(explicitOn) return
nCells = grid % nCells
nVertLevels = grid % nVertLevels
@@ -649,8 +445,7 @@
!> \version SVN:$Id$
!> \details
!> This routine initializes a variety of quantities related to
-!> vertical mixing in the ocean. This primarily determines if
-!> explicit or implicit vertical mixing is to be used.
+!> vertical mixing in the ocean.
!
!-----------------------------------------------------------------------
@@ -674,14 +469,6 @@
velVmixOn = .true.
tracerVmixOn = .true.
- explicitOn = .true.
- implicitOn = .false.
-
- if(config_implicit_vertical_mix) then
- explicitOn = .false.
- implicitOn = .true.
- end if
-
if(config_disable_u_vmix) velVmixOn = .false.
if(config_disable_tr_vmix) tracerVmixOn = .false.
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -228,26 +228,14 @@
if (RiTopOfEdge(k,iEdge)>0.0) then
vertViscTopOfEdge(k,iEdge) = vertViscTopOfEdge(k, iEdge) + config_bkrd_vert_visc &
+ config_rich_mix / (1.0 + 5.0*RiTopOfEdge(k,iEdge))**2
- ! maltrud do limiting of coefficient--should not be necessary
- ! also probably better logic could be found
+ ! maltrud do limiting of coefficient--should not be necessary
+ ! also probably better logic could be found
if (vertViscTopOfEdge(k,iEdge) > config_convective_visc) then
- if( config_implicit_vertical_mix) then
- vertViscTopOfEdge(k,iEdge) = config_convective_visc
- else
- vertViscTopOfEdge(k,iEdge) = &
- ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
- end if
+ vertViscTopOfEdge(k,iEdge) = config_convective_visc
end if
else
- ! mrp 110324 efficiency note: this if is inside iCell and k loops.
- if (config_implicit_vertical_mix) then
- ! for Ri<0 and implicit mix, use convective diffusion
- vertViscTopOfEdge(k,iEdge) = config_convective_visc
- else
- ! for Ri<0 and explicit vertical mix,
- ! use maximum diffusion allowed by CFL criterion
- vertViscTopOfEdge(k,iEdge) = vertViscTopOfEdge(k,iEdge) + ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
- end if
+ ! for Ri<0 use the convective value for the viscosity
+ vertViscTopOfEdge(k,iEdge) = config_convective_visc
end if
end do
end do
@@ -336,23 +324,11 @@
! maltrud do limiting of coefficient--should not be necessary
! also probably better logic could be found
if (vertDiffTopOfCell(k,iCell) > config_convective_diff) then
- if (config_implicit_vertical_mix) then
- vertDiffTopOfCell(k,iCell) = config_convective_diff
- else
- vertDiffTopOfCell(k,iCell) = &
- ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
- end if
+ vertDiffTopOfCell(k,iCell) = config_convective_diff
end if
else
- ! mrp 110324 efficiency note: this if is inside iCell and k loops.
- if (config_implicit_vertical_mix) then
- ! for Ri<0 and implicit mix, use convective diffusion
- vertDiffTopOfCell(k,iCell) = config_convective_diff
- else
- ! for Ri<0 and explicit vertical mix,
- ! use maximum diffusion allowed by CFL criterion
- vertDiffTopOfCell(k,iCell) = vertDiffTopOfCell(k, iCell) + ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
- end if
+ ! for Ri<0 use the convective value for the diffusion
+ vertDiffTopOfCell(k,iCell) = config_convective_diff
end if
end do
end do
Modified: branches/mpas_cdg_advection/src/core_sw/Registry
===================================================================
--- branches/mpas_cdg_advection/src/core_sw/Registry 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/core_sw/Registry 2013-03-13 17:04:50 UTC (rev 2600)
@@ -97,7 +97,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 o edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 o localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 o cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 o cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
Modified: branches/mpas_cdg_advection/src/framework/mpas_block_decomp.F
===================================================================
--- branches/mpas_cdg_advection/src/framework/mpas_block_decomp.F 2013-03-13 16:01:04 UTC (rev 2599)
+++ branches/mpas_cdg_advection/src/framework/mpas_block_decomp.F 2013-03-13 17:04:50 UTC (rev 2600)
@@ -165,7 +165,7 @@
allocate(block_id(blocks_per_proc))
allocate(block_start(blocks_per_proc))
allocate(block_count(blocks_per_proc))
-
+
do i = 1, blocks_per_proc
block_start = 0
block_count = 0
@@ -436,11 +436,11 @@
end if
else
blocks_per_proc = 0
- do i = 1, total_blocks
+ do i = 0, total_blocks-1
call mpas_get_owning_proc(dminfo, i, owning_proc)
if(owning_proc == proc_number) then
call mpas_get_local_block_id(dminfo, i, local_block_id)
- blocks_per_proc = max(blocks_per_proc, local_block_id)
+ blocks_per_proc = max(blocks_per_proc, local_block_id+1)
end if
end do
end if
</font>
</pre>