<p><b>dwj07@fsu.edu</b> 2012-03-30 14:36:26 -0600 (Fri, 30 Mar 2012)</p><p><br>
-- TRUNK COMMIT --<br>
ocean core only.<br>
<br>
Merging gmvar branch to trunk.<br>
<br>
A framework for implementing GM-like parametrization is now in place. Care has been taken to make sure that the introduction of the bolus velocity works correctly with split explicit time stepping scheme. <br>
<br>
In reviewing the code, Mark Petersen has done the following tests:<br>
1. gm_var branch matches trunk bit-for-bit (11-12 digits global KE) for RK4, split explicit, and unsplit, using 120km global for several steps, for zlevel.<br>
<br>
2. convergence tests for timestepping routines on gmvar branch match trunk to at least four digits for zlevel, zstar.<br>
<br>
3. Design algorithm for ALE split explicit revised:<br>
mpas_repo/trunk/documents/ocean/current_design_doc/ale_vert_coord<br>
eqns 3.109-3.112, and 3.116-3.119<br>
<br>
4. Code review of gm implementation - everything looked fine.<br>
<br>
Qingshan Chen has done the following tests:<br>
1. GM on, isopycnal coordinates and RK4 stepping scheme, I got identical results with the before- and after-version-1699 of the gmvar branch.<br>
<br>
2. GM off and again with the isopycnal coordinate and with RK4 stepping scheme, I got identical results for the gmvar branch and for the trunk.<br>
<br>
3. GM off, with the isopycnal coordinate and with split_explicit scheme, trunk and gmvar has a difference of 10^(-9) in the thickness field on day one, which can be considered bit by bit matching.<br>
<br>
4. A series of tests with<br>
GM on, isopycnal coord, RK4, config_pressure_type = 'MontgomeryPotential',<br>
GM on, isopycnal coord, RK4, config_pressure_type = 'pressure',<br>
GM on, isopycnal coord, split explicit, config_pressure_type = 'pressure',<br>
GM on, isopycnal coord, unsplit explicit, config_pressure_type = 'pressure',<br>
the global KE on day one match by 2 digits, which indicates that the code with each of these configurations works correctly.<br>
</p><hr noshade><pre><font color="gray">
Property changes on: trunk/mpas
___________________________________________________________________
Modified: svn:mergeinfo
- /branches/cam_mpas_nh:1260-1270
/branches/ocean_projects/ale_split_exp:1437-1483
/branches/ocean_projects/ale_vert_coord:1225-1383
/branches/ocean_projects/ale_vert_coord_new:1387-1428
/branches/ocean_projects/gmvar:1214-1494
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+ /branches/cam_mpas_nh:1260-1270
/branches/ocean_projects/ale_split_exp:1437-1483
/branches/ocean_projects/ale_vert_coord:1225-1383
/branches/ocean_projects/ale_vert_coord_new:1387-1428
/branches/ocean_projects/gmvar:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
Modified: trunk/mpas/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2012-03-30 16:45:19 UTC (rev 1738)
+++ trunk/mpas/namelist.input.ocean 2012-03-30 20:36:26 UTC (rev 1739)
@@ -2,16 +2,16 @@
config_test_case = 0
config_time_integration = 'split_explicit'
config_rk_filter_btr_mode = .false.
- config_dt = 100.0
+ config_dt = 180.0
config_start_time = '0000-01-01_00:00:00'
- config_run_duration = '2000_00:00:00'
- config_stats_interval = 1920
+ config_run_duration = '1_00:00:00'
+ config_stats_interval = 480
/
&io
config_input_name = 'grid.nc'
- config_output_name = 'output.nc'
+ config_output_name = 'output..nc'
config_restart_name = 'restart.nc'
- config_output_interval = '20_00:00:00'
+ config_output_interval = '1_00:00:00'
config_frames_per_outfile = 1000000
/
&restart
@@ -19,9 +19,9 @@
config_restart_interval = '120_00:00:00'
/
&grid
- config_vert_grid_type = 'zlevel'
+ config_vert_grid_type = 'isopycnal'
config_pressure_type = 'pressure'
- config_rho0 = 1000
+ config_rho0 = 1014.65
/
&split_explicit_ts
config_n_ts_iter = 2
@@ -39,15 +39,17 @@
config_btr_solve_SSH2 = .false.
/
&hmix
- config_h_mom_eddy_visc2 = 1.0e5
+ config_h_mom_eddy_visc2 = 100.0
config_h_mom_eddy_visc4 = 0.0
+ config_h_kappa = 0.0
+ config_h_kappa_q = 0.0
config_visc_vorticity_term = .true.
config_h_tracer_eddy_diff2 = 1.0e5
config_h_tracer_eddy_diff4 = 0.0
/
&vmix
- config_vert_visc_type = 'rich'
- config_vert_diff_type = 'rich'
+ config_vert_visc_type = 'const'
+ config_vert_diff_type = 'const'
config_implicit_vertical_mix = .true.
config_convective_visc = 1.0
config_convective_diff = 1.0
@@ -71,7 +73,7 @@
config_zWidth_tanh = 100
/
&eos
- config_eos_type = 'jm'
+ config_eos_type = 'linear'
/
&advection
config_vert_tracer_adv = 'stencil'
Modified: trunk/mpas/src/core_ocean/Makefile
===================================================================
--- trunk/mpas/src/core_ocean/Makefile 2012-03-30 16:45:19 UTC (rev 1738)
+++ trunk/mpas/src/core_ocean/Makefile 2012-03-30 20:36:26 UTC (rev 1739)
@@ -5,6 +5,7 @@
mpas_ocn_advection.o \
mpas_ocn_thick_hadv.o \
mpas_ocn_thick_vadv.o \
+ mpas_ocn_gm.o \
mpas_ocn_vel_coriolis.o \
mpas_ocn_vel_vadv.o \
mpas_ocn_vel_hmix.o \
@@ -79,6 +80,8 @@
mpas_ocn_thick_vadv.o:
+mpas_ocn_gm.o:
+
mpas_ocn_vel_pressure_grad.o:
mpas_ocn_vel_vadv.o:
@@ -167,6 +170,7 @@
mpas_ocn_test_cases.o \
mpas_ocn_advection.o \
mpas_ocn_thick_hadv.o \
+ mpas_ocn_gm.o \
mpas_ocn_thick_vadv.o \
mpas_ocn_vel_coriolis.o \
mpas_ocn_vel_vadv.o \
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2012-03-30 16:45:19 UTC (rev 1738)
+++ trunk/mpas/src/core_ocean/Registry 2012-03-30 20:36:26 UTC (rev 1739)
@@ -46,6 +46,8 @@
namelist logical hmix config_include_KE_vertex false
namelist real hmix config_h_tracer_eddy_diff2 0.0
namelist real hmix config_h_tracer_eddy_diff4 0.0
+namelist real hmix config_h_kappa 0.0
+namelist real hmix config_h_kappa_q 0.0
namelist logical hmix config_rayleigh_friction false
namelist real hmix config_rayleigh_damping_coeff 0.0
namelist real hmix config_apvm_scale_factor 0.0
@@ -223,6 +225,16 @@
% Diagnostic fields: only written to output
var persistent real zMid ( nVertLevels nCells Time ) 2 io zMid state - -
var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
+var persistent real uTransport ( nVertLevels nEdges Time ) 2 - uTransport state - -
+var persistent real uBolusGM ( nVertLevels nEdges Time ) 2 - uBolusGM state - -
+var persistent real uBolusGMX ( nVertLevels nEdges Time ) 2 - uBolusGMX state - -
+var persistent real uBolusGMY ( nVertLevels nEdges Time ) 2 - uBolusGMY state - -
+var persistent real uBolusGMZ ( nVertLevels nEdges Time ) 2 - uBolusGMZ state - -
+var persistent real uBolusGMZonal ( nVertLevels nEdges Time ) 2 o uBolusGMZonal state - -
+var persistent real uBolusGMMeridional ( nVertLevels nEdges Time ) 2 o uBolusGMMeridional state - -
+var persistent real hEddyFlux ( nVertLevels nEdges Time ) 2 - hEddyFlux state - -
+var persistent real h_kappa ( nVertLevels nEdges Time ) 2 - h_kappa state - -
+var persistent real h_kappa_q ( nVertLevels nEdges Time ) 2 - h_kappa_q state - -
var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
var persistent real Vor_edge ( nVertLevels nEdges Time ) 2 - Vor_edge state - -
Copied: trunk/mpas/src/core_ocean/mpas_ocn_gm.F (from rev 1513, trunk/mpas/src/core_ocean/mpas_ocn_gm.F)
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_gm.F (rev 0)
+++ trunk/mpas/src/core_ocean/mpas_ocn_gm.F 2012-03-30 20:36:26 UTC (rev 1739)
@@ -0,0 +1,142 @@
+module ocn_gm
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_timer
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: ocn_gm_compute_uBolus
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+contains
+
+ subroutine ocn_gm_compute_uBolus(s, grid)
+ implicit none
+ type(state_type), intent(inout) :: s
+ type(mesh_type), intent(in) :: grid
+
+ real(kind=RKIND), dimension(:,:), pointer :: uBolusGM, hEddyFlux, h_edge
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer :: k, iEdge, nEdges
+
+ uBolusGM => s % uBolusGM % array
+ h_edge => s % h_edge % array
+ hEddyFlux => s % hEddyFlux % array
+
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+
+ nEdges = grid % nEdges
+
+ call ocn_gm_compute_hEddyFlux(s, grid)
+
+ if (config_vert_grid_type .EQ. 'isopycnal') then
+
+ do iEdge = 1, nEdges
+ do k = 1, maxLevelEdgeTop(iEdge)
+ uBolusGM(k,iEdge) = hEddyFlux(k,iEdge)/h_edge(k,iEdge)
+ end do
+ end do
+
+ else
+
+ ! Nothing for now for all other grid types (zlevel, zstar, ztilde)
+ uBolusGM(:,:) = 0.0
+
+ end if
+
+ end subroutine ocn_gm_compute_uBolus
+
+
+ subroutine ocn_gm_compute_hEddyFlux(s, grid)
+ implicit none
+ type(state_type), intent(inout) :: s
+ type(mesh_type), intent(in) :: grid
+
+ real(kind=RKIND), dimension(:,:), pointer :: hEddyFlux, h
+ real(kind=RKIND), dimension(:), pointer :: dcEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer :: k, cell1, cell2, iEdge, nEdges
+
+ hEddyFlux => s % hEddyFlux % array
+ h => s % h % array
+
+ dcEdge => grid % dcEdge % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+
+ nEdges = grid % nEdges
+
+ hEddyFlux(:,:) = 0.0
+
+ if (config_vert_grid_type .EQ. 'isopycnal') then
+ do iEdge = 1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ hEddyFlux(k,iEdge) = -config_h_kappa * (h(k,cell2) - h(k,cell1)) / dcEdge(iEdge)
+ end do
+ end do
+ else
+
+ !Nothing for now for all other grid types (zlevel, zstar, ztilde)
+
+ end if
+
+ end subroutine ocn_gm_compute_hEddyFlux
+
+
+
+ subroutine ocn_get_h_kappa(s, grid)
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ real(kind=RKIND), dimension(:,:), pointer :: h_kappa
+
+
+ h_kappa => s % h_kappa % array
+
+ h_kappa(:,:) = config_h_kappa
+
+
+ end subroutine ocn_get_h_kappa
+
+
+ subroutine ocn_get_h_kappa_q(s, grid)
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ real(kind=RKIND), dimension(:,:), pointer :: h_kappa_q
+
+
+ h_kappa_q => s % h_kappa_q % array
+
+ h_kappa_q(:,:) = config_h_kappa_q
+
+
+ end subroutine ocn_get_h_kappa_q
+
+end module ocn_gm
Modified: trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-30 16:45:19 UTC (rev 1738)
+++ trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-30 20:36:26 UTC (rev 1739)
@@ -19,6 +19,7 @@
use ocn_tracer_hadv
use ocn_tracer_vadv
use ocn_tracer_hmix
+ use ocn_gm
use ocn_restoring
use ocn_equation_of_state
@@ -247,6 +248,11 @@
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
+ block % state % time_levs(1) % state % uTransport % array(:,:) &
+ = block % state % time_levs(1) % state % u % array(:,:) &
+ + block % state % time_levs(1) % state % uBolusGM % array(:,:)
+
call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(1) % state, mesh)
call ocn_compute_mesh_scaling(mesh)
Modified: trunk/mpas/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-03-30 16:45:19 UTC (rev 1738)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-03-30 20:36:26 UTC (rev 1739)
@@ -25,6 +25,7 @@
use ocn_thick_hadv
use ocn_thick_vadv
+ use ocn_gm
use ocn_vel_coriolis
use ocn_vel_pressure_grad
@@ -106,13 +107,13 @@
type (state_type), intent(in) :: s !< Input: State information
type (mesh_type), intent(in) :: grid !< Input: Grid information
- real (kind=RKIND), dimension(:,:), pointer :: h_edge, u, wTop, tend_h
+ real (kind=RKIND), dimension(:,:), pointer :: h_edge, wTop, tend_h, uTransport
integer :: err
call mpas_timer_start("ocn_tend_h")
- u => s % u % array
+ uTransport => s % uTransport % array
wTop => s % wTop % array
h_edge => s % h_edge % array
@@ -129,8 +130,10 @@
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
! for explanation of divergence operator.
!
+ ! QC Comment (3/15/12): need to make sure that uTranport is the right
+ ! transport velocity here.
call mpas_timer_start("hadv", .false., thickHadvTimer)
- call ocn_thick_hadv_tend(grid, u, h_edge, tend_h, err)
+ call ocn_thick_hadv_tend(grid, uTransport, h_edge, tend_h, err)
call mpas_timer_stop("hadv", thickHadvTimer)
!
@@ -279,7 +282,7 @@
real (kind=RKIND), intent(in) :: dt !< Input: Time step
real (kind=RKIND), dimension(:,:), pointer :: &
- u, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh
+ uTransport, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
@@ -287,7 +290,7 @@
call mpas_timer_start("ocn_tend_scalar")
- u => s % u % array
+ uTransport => s % uTransport % array
h => s % h % array
wTop => s % wTop % array
tracers => s % tracers % array
@@ -298,11 +301,13 @@
tend_h => tend % h % array
allocate(uh(grid % nVertLevels, grid % nEdges+1))
-
+ !
+ ! QC Comment (3/15/12): need to make sure that uTransport is the right
+ ! transport velocity for the tracer.
do iEdge = 1, grid % nEdges
- do k = 1, grid % nVertLevels
- uh(k, iEdge) = u(k, iEdge) * h_edge(k, iEdge)
- end do
+ do k = 1, grid % nVertLevels
+ uh(k, iEdge) = uTransport(k, iEdge) * h_edge(k, iEdge)
+ end do
end do
!
@@ -323,7 +328,6 @@
call mpas_ocn_tracer_advection_tend(tracers, uh, wTop, h, h, dt, grid, tend_h, tend_tr)
call mpas_timer_stop("adv", tracerHadvTimer)
-
!
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
!
@@ -411,13 +415,15 @@
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
Vor_edge, Vor_vertex, Vor_cell, gradVor_n, gradVor_t, divergence, &
- rho, temperature, salinity, kev, kevc
+ rho, temperature, salinity, kev, kevc, uBolusGM, uTransport
real (kind=RKIND), dimension(:,:,:), pointer :: tracers, deriv_two
real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
h => s % h % array
u => s % u % array
+ uTransport => s % uTransport % array
+ uBolusGM => s % uBolusGM % array
v => s % v % array
h_edge => s % h_edge % array
circulation => s % circulation % array
@@ -469,6 +475,7 @@
boundaryCell => grid % boundaryCell % array
+
!
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
@@ -831,6 +838,16 @@
end do
+ !
+ ! Apply the GM closure as a bolus velocity
+ !
+ if (config_h_kappa .GE. epsilon(0D0)) then
+ call ocn_gm_compute_uBolus(s,grid)
+ else
+ ! mrp efficiency note: if uBolusGM is guaranteed to be zero, this can be removed.
+ uBolusGM = 0.0
+ end if
+
end subroutine ocn_diagnostic_solve!}}}
!***********************************************************************
@@ -860,7 +877,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: u,h,wTop, h_edge
+ real (kind=RKIND), dimension(:,:), pointer :: uTransport,h,wTop, h_edge
real (kind=RKIND), dimension(:,:), allocatable:: div_hu
real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col, h_weights
@@ -873,7 +890,7 @@
h => s1 % h % array
h_edge => s1 % h_edge % array
- u => s2 % u % array
+ uTransport => s2 % uTransport % array
wTop => s2 % wTop % array
areaCell => grid % areaCell % array
@@ -898,7 +915,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeBot(iEdge)
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ flux = uTransport(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
div_hu(k,cell1) = div_hu(k,cell1) + flux
div_hu(k,cell2) = div_hu(k,cell2) - flux
end do
@@ -1078,6 +1095,7 @@
end subroutine ocn_fuperp!}}}
+
!***********************************************************************
!
! routine ocn_tendency_init
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-03-30 16:45:19 UTC (rev 1738)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-03-30 20:36:26 UTC (rev 1739)
@@ -241,6 +241,11 @@
call ocn_diagnostic_solve(dt, provis, block % mesh)
+ ! Compute velocity transport, used in advection terms of h and tracer tendancy
+ provis % uTransport % array(:,:) &
+ = provis % u % array(:,:) &
+ + provis % uBolusGM % array(:,:)
+
block => block % next
end do
end if
@@ -343,6 +348,11 @@
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
+ 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, &
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-30 16:45:19 UTC (rev 1738)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-30 20:36:26 UTC (rev 1739)
@@ -312,6 +312,20 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % uBcl % array(:,:)
+ do iEdge=1,block % mesh % nEdges
+ do k=1,block % mesh % nVertLevels
+
+ ! uTranport = uBcl + uBolus
+ ! This is u used in advective terms for h and tracers
+ ! in tendancy 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) &
+ + block % state % time_levs(1) % state % uBolusGM % array(k,iEdge) )
+
+ enddo
+ end do ! iEdge
+
block => block % next
end do ! block
@@ -663,9 +677,11 @@
do iEdge=1,block % mesh % nEdges
- ! This is u^{avg}
- uTemp(:) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
- + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
+ ! velocity for uCorrection is uBtr + uBcl + uBolus
+ uTemp(:) &
+ = block % state % time_levs(2) % state % uBtr % array( iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(:,iEdge) &
+ + block % state % time_levs(1) % state % uBolusGM % array(:,iEdge)
! hSum is initialized outside the loop because on land boundaries
! maxLevelEdgeTop=0, but I want to initialize hSum with a
@@ -680,19 +696,20 @@
uCorr = ucorr_coef*(( block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum)
- ! put u^{tr}, the velocity for tracer transport, in uNew
- ! mrp 060611 not sure if boundary enforcement is needed here.
- if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
- block % state % time_levs(2) % state % u % array(:,iEdge) = 0.0
- else
- do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
- block % state % time_levs(2) % state % u % array(k,iEdge) = uTemp(k) + uCorr
- enddo
- do k=block % mesh % maxLevelEdgeTop % array(iEdge)+1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % u % array(k,iEdge) = 0.0
- end do
- endif
+ do k=1,block % mesh % nVertLevels
+ ! uTranport = uBtr + uBcl + uBolus + uCorrection
+ ! This is u used in advective terms for h and tracers
+ ! in tendancy 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) &
+ + block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
+ + block % state % time_levs(1) % state % uBolusGM % array(k,iEdge) &
+ + uCorr )
+
+ enddo
+
end do ! iEdge
deallocate(uTemp)
@@ -801,10 +818,22 @@
end do
end do ! iCell
- ! uBclNew is u'_{n+1/2}
- ! uBtrNew is {\bar u}_{avg}
- ! uNew is u^{tr}
+ do iEdge=1,block % mesh % nEdges
+ do k=1,block % mesh % nVertLevels
+
+ ! u = uBtr + uBcl
+ ! here uBcl is at time n+1/2
+ ! This is u used in next iteration or step
+ block % state % time_levs(2) % state % u % array(k,iEdge) &
+ = block % mesh % edgeMask % array(k,iEdge) &
+ *( block % state % time_levs(2) % state % uBtr % array( iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(k,iEdge) )
+
+ enddo
+
+ end do ! iEdge
+
! mrp 110512 I really only need this to compute h_edge, density, pressure, and SSH
! I can par this down later.
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
</font>
</pre>