<p><b>qchen3@fsu.edu</b> 2012-04-12 14:38:58 -0600 (Thu, 12 Apr 2012)</p><p>BRANCH COMMIT<br>
<br>
Merged updates from trunk.<br>
</p><hr noshade><pre><font color="gray">
Property changes on: branches/ocean_projects/accgm
___________________________________________________________________
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
/trunk/mpas:1618-1679
+ /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/ocean_projects/zstar_restart_new:1762-1770
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
/trunk/mpas:1618-1778
Modified: branches/ocean_projects/accgm/namelist.input.ocean
===================================================================
--- branches/ocean_projects/accgm/namelist.input.ocean 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/namelist.input.ocean 2012-04-12 20:38:58 UTC (rev 1779)
@@ -9,7 +9,7 @@
/
&io
config_input_name = 'grid.nc'
- config_output_name = 'output.nc'
+ config_output_name = 'output..nc'
config_restart_name = 'restart.nc'
config_output_interval = '1_00:00:00'
config_frames_per_outfile = 1000000
@@ -41,6 +41,8 @@
&hmix
config_h_mom_eddy_visc2 = 0.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
@@ -71,7 +73,7 @@
config_zWidth_tanh = 100
/
&eos
- config_eos_type = 'jm'
+ config_eos_type = 'linear'
/
&advection
config_vert_tracer_adv = 'stencil'
Modified: branches/ocean_projects/accgm/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/Makefile 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/src/core_ocean/Makefile 2012-04-12 20:38:58 UTC (rev 1779)
@@ -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: branches/ocean_projects/accgm/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/Registry 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/src/core_ocean/Registry 2012-04-12 20:38:58 UTC (rev 1779)
@@ -24,6 +24,7 @@
namelist character grid config_vert_grid_type isopycnal
namelist character grid config_pressure_type pressure
namelist real grid config_rho0 1028
+namelist logical grid config_enforce_zstar_at_restart false
namelist integer split_explicit_ts config_n_ts_iter 2
namelist integer split_explicit_ts config_n_bcl_iter_beg 2
namelist integer split_explicit_ts config_n_bcl_iter_mid 2
@@ -46,6 +47,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
@@ -183,6 +186,7 @@
var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
var persistent real referenceBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - referenceBottomDepthTopOfCell mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
+var persistent real zstarWeight ( nVertLevels ) 0 - zstarWeight mesh - -
% Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
@@ -224,6 +228,16 @@
var persistent real z ( nVertLevels nCells Time ) 2 o z state - -
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 - -
Modified: branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_advection.F
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_advection.F 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_advection.F 2012-04-12 20:38:58 UTC (rev 1779)
@@ -123,9 +123,9 @@
do i=1,n
advCells(i+1) = cell_list(i)
- xc(i) = grid % xCell % array(advCells(i+1))/a
- yc(i) = grid % yCell % array(advCells(i+1))/a
- zc(i) = grid % zCell % array(advCells(i+1))/a
+ xc(i) = grid % xCell % array(advCells(i+1))/grid % sphere_radius
+ yc(i) = grid % yCell % array(advCells(i+1))/grid % sphere_radius
+ zc(i) = grid % zCell % array(advCells(i+1))/grid % sphere_radius
end do
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
@@ -143,7 +143,7 @@
xc(i+1), yc(i+1), zc(i+1), &
xc(ip2), yc(ip2), zc(ip2) )
- dl_sphere(i) = a*arc_length( xc(1), yc(1), zc(1), &
+ dl_sphere(i) = grid % sphere_radius*arc_length( xc(1), yc(1), zc(1), &
xc(i+1), yc(i+1), zc(i+1) )
end do
@@ -172,8 +172,8 @@
if ( iCell .ne. grid % CellsOnEdge % array(1,iEdge)) &
angle_2d(i) = angle_2d(i) - pii
-! xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
-! yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+! xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
+! yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
xp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * cos(angle_2d(i))
yp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * sin(angle_2d(i))
@@ -274,12 +274,12 @@
if (ip1 > n-1) ip1 = 1
iEdge = grid % EdgesOnCell % array (i,iCell)
- xv1 = grid % xVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- yv1 = grid % yVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- zv1 = grid % zVertex % array(grid % verticesOnEdge % array (1,iedge))/a
- xv2 = grid % xVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/a
- zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/a
+ xv1 = grid % xVertex % array(grid % verticesOnEdge % array (1,iedge))/grid % sphere_radius
+ yv1 = grid % yVertex % array(grid % verticesOnEdge % array (1,iedge))/grid % sphere_radius
+ zv1 = grid % zVertex % array(grid % verticesOnEdge % array (1,iedge))/grid % sphere_radius
+ xv2 = grid % xVertex % array(grid % verticesOnEdge % array (2,iedge))/grid % sphere_radius
+ yv2 = grid % yVertex % array(grid % verticesOnEdge % array (2,iedge))/grid % sphere_radius
+ zv2 = grid % zVertex % array(grid % verticesOnEdge % array (2,iedge))/grid % sphere_radius
if ( grid % on_a_sphere ) then
call ocn_arc_bisect( xv1, yv1, zv1, &
Copied: branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_gm.F (from rev 1778, trunk/mpas/src/core_ocean/mpas_ocn_gm.F)
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_gm.F (rev 0)
+++ branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_gm.F 2012-04-12 20:38:58 UTC (rev 1779)
@@ -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: branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_mpas_core.F 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_mpas_core.F 2012-04-12 20:38:58 UTC (rev 1779)
@@ -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
@@ -104,6 +105,12 @@
call ocn_init_z_level(domain)
+ if (config_enforce_zstar_at_restart) then
+ call ocn_init_h_zstar(domain)
+ endif
+
+ call ocn_init_split_timestep(domain)
+
print *, ' Vertical grid type is: ',config_vert_grid_type
if (config_vert_grid_type.ne.'isopycnal'.and. &
@@ -247,6 +254,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)
@@ -493,7 +505,7 @@
end subroutine mpas_timestep!}}}
subroutine ocn_init_z_level(domain)!{{{
- ! Initialize maxLevel and bouncary grid variables.
+ ! Initialize zlevel-type variables
use mpas_grid_types
use mpas_configure
@@ -507,7 +519,8 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, referenceBottomDepthTopOfCell
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, &
+ referenceBottomDepthTopOfCell, zstarWeight, hZLevel
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -519,6 +532,8 @@
h => block % state % time_levs(1) % state % h % array
referenceBottomDepth => block % mesh % referenceBottomDepth % array
referenceBottomDepthTopOfCell => block % mesh % referenceBottomDepthTopOfCell % array
+ zstarWeight => block % mesh % zstarWeight % array
+ hZLevel => block % mesh % hZLevel % array
nVertLevels = block % mesh % nVertLevels
! mrp 120208 right now hZLevel is in the grid.nc file.
@@ -526,9 +541,9 @@
! as the defining variable instead, and will transition soon.
! When the transition is done, hZLevel can be removed from
! registry and the following four lines deleted.
- referenceBottomDepth(1) = block % mesh % hZLevel % array(1)
+ referenceBottomDepth(1) = hZLevel(1)
do k = 2,nVertLevels
- referenceBottomDepth(k) = referenceBottomDepth(k-1) + block % mesh % hZLevel % array(k)
+ referenceBottomDepth(k) = referenceBottomDepth(k-1) + hZLevel(k)
end do
! TopOfCell needed where zero depth for the very top may be referenced.
@@ -537,6 +552,65 @@
referenceBottomDepthTopOfCell(k+1) = referenceBottomDepth(k)
end do
+ ! Initialization of zstarWeights. This determines how SSH perturbations
+ ! are distributed throughout the column.
+ if (config_vert_grid_type.eq.'zlevel') then
+
+ zstarWeight = 0.0
+ zstarWeight(1) = 1.0
+
+ elseif (config_vert_grid_type.eq.'zstar') then
+
+ do k = 1,nVertLevels
+ zstarWeight(k) = hZLevel(k)
+ enddo
+
+ elseif (config_vert_grid_type.eq.'zstarWeights') then
+
+ ! This is a test with other weights, just to make sure zstar functions
+ ! using variable weights.
+
+ zstarWeight = 0.0
+ zstarWeight(1:5) = 1.0
+ do k=1,10
+ zstarWeight(5+k) = 1.0-k*0.1
+ end do
+
+ endif
+
+ block => block % next
+ end do
+
+ end subroutine ocn_init_z_level!}}}
+
+ subroutine ocn_init_split_timestep(domain)!{{{
+ ! Initialize splitting variables
+
+ use mpas_grid_types
+ use mpas_configure
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ integer :: iTracer, cell, cell1, cell2
+ real (kind=RKIND) :: uhSum, hSum, hEdge1
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+
+ real (kind=RKIND), dimension(:,:), pointer :: h
+ integer :: nVertLevels
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ h => block % state % time_levs(1) % state % h % array
+ referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ nVertLevels = block % mesh % nVertLevels
+
! Compute barotropic velocity at first timestep
! This is only done upon start-up.
if (trim(config_time_integration) == 'unsplit_explicit') then
@@ -612,8 +686,66 @@
block => block % next
end do
- end subroutine ocn_init_z_level!}}}
+ end subroutine ocn_init_split_timestep!}}}
+ subroutine ocn_init_h_zstar(domain)!{{{
+ ! If changing from zlevel to zstar, compute h based on zstar weights,
+ ! where SSH is distributed through the layers. We only change h.
+ ! We do not remap the tracer variables, so this breaks total global
+ ! conservation.
+
+ use mpas_grid_types
+ use mpas_configure
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ type (block_type), pointer :: block
+
+ integer :: i, iCell, iEdge, iVertex, k, nVertLevels
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: hSum, sumZstarWeights
+ real (kind=RKIND), dimension(:), pointer :: hZLevel, zstarWeight, &
+ referenceBottomDepth
+ real (kind=RKIND), dimension(:,:), pointer :: h
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ h => block % state % time_levs(1) % state % h % array
+ nVertLevels = block % mesh % nVertLevels
+ hZLevel => block % mesh % hZLevel % array
+ maxLevelCell => block % mesh % maxLevelCell % array
+ zstarWeight => block % mesh % zstarWeight % array
+ referenceBottomDepth => block % mesh % referenceBottomDepth % array
+
+ do iCell=1,block % mesh % nCells
+ ! Compute the total column thickness, hSum, and the sum of zstar weights.
+ hSum = 0.0
+ sumZstarWeights = 0.0
+ do k = 1,maxLevelCell(iCell)
+ hSum = hSum + h(k,iCell)
+ sumZstarWeights = sumZstarWeights + zstarWeight(k)
+ enddo
+
+ ! h_k = h_k^{zlevel} + zeta * W_k/sum(W_k)
+ ! where zeta is SSH and W_k are weights
+ do k = 1,maxLevelCell(iCell)
+ h(k,iCell) = hZLevel(k) &
+ + (hSum - referenceBottomDepth(maxLevelCell(iCell))) &
+ * zstarWeight(k)/sumZstarWeights
+ enddo
+
+ enddo
+
+ block => block % next
+ end do
+
+ end subroutine ocn_init_h_zstar!}}}
+
subroutine ocn_compute_max_level(domain)!{{{
! Initialize maxLevel and bouncary grid variables.
Modified: branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_tendency.F 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_tendency.F 2012-04-12 20:38:58 UTC (rev 1779)
@@ -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, z
+ rho, temperature, salinity, kev, kevc, z, 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
z => s % z % array
h_edge => s % h_edge % array
@@ -471,6 +477,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
@@ -853,6 +860,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!}}}
!***********************************************************************
@@ -881,10 +898,10 @@
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: u,h,wTop, h_edge
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zstarWeight
+ 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
+ real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
@@ -895,7 +912,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
@@ -903,13 +920,14 @@
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
dvEdge => grid % dvEdge % array
+ zstarWeight => grid % zstarWeight % array
nCells = grid % nCells
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
allocate(div_hu(nVertLevels,nCells+1), div_hu_btr(nCells+1), &
- h_tend_col(nVertLevels), h_weights(nVertLevels))
+ h_tend_col(nVertLevels))
!
! Compute div(h^{edge} u) for each cell
@@ -920,7 +938,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
@@ -942,51 +960,14 @@
! set vertical velocity to zero in isopycnal case
wTop=0.0
- elseif (config_vert_grid_type.eq.'zlevel') then
+ else ! zlevel or zstar type vertical grid
do iCell=1,nCells
- ! Vertical velocity through layer interface at top and
- ! bottom is zero.
- wTop(1,iCell) = 0.0
- wTop(maxLevelCell(iCell)+1,iCell) = 0.0
- do k=maxLevelCell(iCell),2,-1
- wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell)
- end do
- end do
- elseif (config_vert_grid_type.eq.'zstar1') then
-
- ! This is a testing setting. The computation is similar to zstar,
- ! but the weights are all in the top layer, so is a bit-for-bit
- ! match with zlevel.
-
- do iCell=1,nCells
-
- h_tend_col = 0.0
- h_tend_col(1) = - div_hu_btr(iCell)
-
- ! Vertical velocity through layer interface at top and
- ! bottom is zero.
- wTop(1,iCell) = 0.0
- wTop(maxLevelCell(iCell)+1,iCell) = 0.0
-
- do k=maxLevelCell(iCell),2,-1
- wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell) - h_tend_col(k)
- end do
- end do
-
- elseif (config_vert_grid_type.eq.'zstar') then
-
- ! Distribute the change in total column height due to the external
- ! mode, div_hu_btr, among all the layers. Distribute in proportion
- ! to the layer thickness.
-
- do iCell=1,nCells
-
hSum = 0.0
do k=1,maxLevelCell(iCell)
- h_tend_col(k) = - h(k,iCell)*div_hu_btr(iCell)
- hSum = hSum + h(k,iCell)
+ h_tend_col(k) = - zstarWeight(k)*h(k,iCell)*div_hu_btr(iCell)
+ hSum = hSum + zstarWeight(k)*h(k,iCell)
end do
h_tend_col = h_tend_col / hSum
@@ -999,37 +980,9 @@
end do
end do
- elseif (config_vert_grid_type.eq.'zstarWeights') then
-
- ! This is a test with other weights, not meant to be permanent.
-
- h_weights = 0.0
- h_weights(1:5) = 1.0
- do k=1,10
- h_weights(5+k) = 1.0-k*0.1
- end do
-
- do iCell=1,nCells
-
- hSum = 0.0
- do k=1,maxLevelCell(iCell)
- h_tend_col(k) = - h_weights(k)*h(k,iCell)*div_hu_btr(iCell)
- hSum = hSum + h_weights(k)*h(k,iCell)
- end do
- h_tend_col = h_tend_col / hSum
-
- ! Vertical velocity through layer interface at top and
- ! bottom is zero.
- wTop(1,iCell) = 0.0
- wTop(maxLevelCell(iCell)+1,iCell) = 0.0
- do k=maxLevelCell(iCell),2,-1
- wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell) - h_tend_col(k)
- end do
- end do
-
endif
- deallocate(div_hu, div_hu_btr, h_tend_col, h_weights)
+ deallocate(div_hu, div_hu_btr, h_tend_col)
end subroutine ocn_wtop!}}}
@@ -1100,6 +1053,7 @@
end subroutine ocn_fuperp!}}}
+
!***********************************************************************
!
! routine ocn_tendency_init
Modified: branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-04-12 20:38:58 UTC (rev 1779)
@@ -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: branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_time_integration_split.F 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_time_integration_split.F 2012-04-12 20:38:58 UTC (rev 1779)
@@ -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)
Modified: branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-04-12 18:44:46 UTC (rev 1778)
+++ branches/ocean_projects/accgm/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-04-12 20:38:58 UTC (rev 1779)
@@ -324,19 +324,17 @@
if (config_check_monotonicity) then
!build min and max bounds on old and new tracer for check on monotonicity.
- cur_min = minval(tracer_cur(:,1:nCellsSolve))
- cur_max = maxval(tracer_cur(:,1:nCellsSolve))
- new_min = minval(tracer_new(:,1:nCellsSolve))
- new_max = maxval(tracer_new(:,1:nCellsSolve))
+ do iCell = 1, nCellsSolve
+ do k = 1, maxLevelCell(iCell)
+ if(tracer_new(k,iCell) < tracer_min(k, iCell)-eps) then
+ write(*,*) 'Minimum out of bounds on tracer ', iTracer, tracer_min(k, iCell), tracer_new(k,iCell)
+ end if
- !check monotonicity
- if(new_min < cur_min-eps) then
- write(*,*) 'Minimum out of bounds on tracer ', iTracer, cur_min, new_min
- end if
-
- if(new_max > cur_max+eps) then
- write(*,*) 'Maximum out of bounds on tracer ', iTracer, cur_max, new_max
- end if
+ if(tracer_new(k,iCell) > tracer_max(k,iCell)+eps) then
+ write(*,*) 'Maximum out of bounds on tracer ', iTracer, tracer_max(k, iCell), tracer_new(k,iCell)
+ end if
+ end do
+ end do
end if
end do ! iTracer loop
</font>
</pre>