<p><b>mhoffman@lanl.gov</b> 2012-05-15 11:21:01 -0600 (Tue, 15 May 2012)</p><p>BRANCH COMMIT -- land_ice<br>
<br>
Merging the trunk version of the ocean core to the land ice implement_core branch to allow access to the most up-to-date version of the advection routines there.<br>
r1652 through r1912<br>
</p><hr noshade><pre><font color="gray">
Property changes on: branches/land_ice_projects/implement_core/src/core_ocean
___________________________________________________________________
Added: svn:mergeinfo
+ /branches/cam_mpas_nh/src/core_ocean:1260-1270
/branches/ocean_projects/ale_split_exp/src/core_ocean:1437-1483
/branches/ocean_projects/ale_vert_coord/src/core_ocean:1225-1383
/branches/ocean_projects/ale_vert_coord_new/src/core_ocean:1387-1428
/branches/ocean_projects/gmvar/src/core_ocean:1214-1494
/branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean:754-986
/branches/ocean_projects/monotonic_advection/src/core_ocean:1499-1640
/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
/branches/source_renaming/src/core_ocean:1082-1113
/branches/time_manager/src/core_ocean:924-962
/trunk/mpas/src/core_ocean:1561-1912
Modified: branches/land_ice_projects/implement_core/src/core_ocean/Makefile
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/Makefile 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/Makefile 2012-05-15 17:21:01 UTC (rev 1913)
@@ -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 \
@@ -52,7 +53,8 @@
mpas_ocn_equation_of_state_jm.o \
mpas_ocn_equation_of_state_linear.o \
mpas_ocn_global_diagnostics.o \
- mpas_ocn_time_average.o
+ mpas_ocn_time_average.o \
+ mpas_ocn_monthly_forcing.o
all: core_hyd
@@ -79,6 +81,8 @@
mpas_ocn_thick_vadv.o:
+mpas_ocn_gm.o:
+
mpas_ocn_vel_pressure_grad.o:
mpas_ocn_vel_vadv.o:
@@ -163,10 +167,13 @@
mpas_ocn_equation_of_state_linear.o:
+mpas_ocn_monthly_forcing.o:
+
mpas_ocn_mpas_core.o: mpas_ocn_mpas_core.o \
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 \
@@ -214,7 +221,8 @@
mpas_ocn_equation_of_state_jm.o \
mpas_ocn_equation_of_state_linear.o \
mpas_ocn_global_diagnostics.o \
- mpas_ocn_time_average.o
+ mpas_ocn_time_average.o \
+ mpas_ocn_monthly_forcing.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/land_ice_projects/implement_core/src/core_ocean/Registry
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/Registry 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/Registry 2012-05-15 17:21:01 UTC (rev 1913)
@@ -18,12 +18,18 @@
namelist character io config_restart_name restart.nc
namelist character io config_output_interval 24:00:00
namelist integer io config_frames_per_outfile 0
-namelist character io config_decomp_file_prefix graph.info.part.
+namelist integer io config_pio_num_iotasks 0
+namelist integer io config_pio_stride 1
+namelist character decomposition config_block_decomp_file_prefix graph.info.part.
+namelist integer decomposition config_number_of_blocks 0
+namelist logical decomposition config_explicit_proc_decomp .false.
+namelist character decomposition config_proc_decomp_file_prefix graph.info.part.
namelist logical restart config_do_restart false
namelist character restart config_restart_interval none
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 +52,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
@@ -77,6 +85,7 @@
namelist logical restore config_restoreTS false
namelist real restore config_restoreT_timescale 90.0
namelist real restore config_restoreS_timescale 90.0
+namelist logical restore config_use_monthly_forcing false
%
% dim type name_in_file name_in_code
@@ -94,6 +103,7 @@
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
dim nVertLevelsP1 nVertLevels+1
+dim nMonths nMonths
%
% var persistence type name_in_file ( dims ) time_levs iro- name_in_code struct super-array array_class
@@ -124,6 +134,7 @@
var persistent real meshDensity ( nCells ) 0 iro meshDensity mesh - -
var persistent real meshScalingDel2 ( nEdges ) 0 ro meshScalingDel2 mesh - -
var persistent real meshScalingDel4 ( nEdges ) 0 ro meshScalingDel4 mesh - -
+var persistent real meshScaling ( nEdges ) 0 ro meshScaling mesh - -
var persistent integer cellsOnEdge ( TWO nEdges ) 0 iro cellsOnEdge mesh - -
var persistent integer nEdgesOnCell ( nCells ) 0 iro nEdgesOnCell mesh - -
@@ -183,6 +194,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 - -
@@ -191,10 +203,15 @@
var persistent integer edgeMask ( nVertLevels nEdges ) 0 o edgeMask mesh - -
var persistent integer vertexMask ( nVertLevels nVertices ) 0 o vertexMask mesh - -
var persistent integer cellMask ( nVertLevels nCells ) 0 o cellMask mesh - -
-var persistent real u_src ( nVertLevels nEdges ) 0 ir u_src mesh - -
-var persistent real temperatureRestore ( nCells ) 0 ir temperatureRestore mesh - -
-var persistent real salinityRestore ( nCells ) 0 ir salinityRestore mesh - -
+var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
+var persistent real temperatureRestore ( nCells ) 0 iro temperatureRestore mesh - -
+var persistent real salinityRestore ( nCells ) 0 iro salinityRestore mesh - -
+% mrp trying to figure out why these do not appear
+var persistent real windStressMonthly ( nMonths nEdges ) 0 iro windStressMonthly mesh - -
+var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 iro temperatureRestoreMonthly mesh - -
+var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 iro salinityRestoreMonthly mesh - -
+
% Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 ir u state - -
var persistent real h ( nVertLevels nCells Time ) 2 iro h state - -
@@ -223,6 +240,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 - -
@@ -239,6 +266,11 @@
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 - uReconstructZ state - -
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
+var persistent real uSrcReconstructX ( nVertLevels nCells Time ) 2 - uSrcReconstructX state - -
+var persistent real uSrcReconstructY ( nVertLevels nCells Time ) 2 - uSrcReconstructY state - -
+var persistent real uSrcReconstructZ ( nVertLevels nCells Time ) 2 - uSrcReconstructZ state - -
+var persistent real uSrcReconstructZonal ( nVertLevels nCells Time ) 2 o uSrcReconstructZonal state - -
+var persistent real uSrcReconstructMeridional ( nVertLevels nCells Time ) 2 o uSrcReconstructMeridional state - -
var persistent real MontPot ( nVertLevels nCells Time ) 2 - MontPot state - -
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
@@ -272,5 +304,5 @@
var persistent real acc_uReconstructMeridional ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridional state - -
var persistent real acc_uReconstructZonalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructZonalVar state - -
var persistent real acc_uReconstructMeridionalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridionalVar state - -
-var persistent real acc_u ( nVertLevels nEdges Time ) 2 o acc_u state - -
-var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 o acc_uVar state - -
+var persistent real acc_u ( nVertLevels nEdges Time ) 2 o acc_u state - -
+var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 o acc_uVar state - -
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_advection.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_advection.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_advection.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -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, &
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_equation_of_state.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_equation_of_state.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_equation_of_state.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -14,6 +14,7 @@
module ocn_equation_of_state
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
use ocn_equation_of_state_linear
@@ -86,7 +87,7 @@
type (mesh_type), intent(in) :: grid
integer, intent(out) :: err
integer :: k_displaced
- character(len=8), intent(in) :: displacement_type
+ character(len=*), intent(in) :: displacement_type
integer, dimension(:), pointer :: maxLevelCell
real (kind=RKIND), dimension(:,:), pointer :: rho
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_equation_of_state_jm.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -14,6 +14,7 @@
module ocn_equation_of_state_jm
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
@@ -85,7 +86,7 @@
type (mesh_type), intent(in) :: grid
integer :: k_displaced, indexT, indexS
- character(len=8), intent(in) :: displacement_type
+ character(len=*), intent(in) :: displacement_type
integer, intent(out) :: err
type (dm_info) :: dminfo
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_global_diagnostics.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_global_diagnostics.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_global_diagnostics.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -46,13 +46,15 @@
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal, localCFL, localSum, areaCellGlobal, areaEdgeGlobal, areaTriangleGlobal
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, areaEdge
- real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, circulation, vorticity, ke, Vor_edge, Vor_vertex, &
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, vorticity, ke, Vor_edge, Vor_vertex, &
Vor_cell, gradVor_n, gradVor_t, pressure, MontPot, wTop, rho, tracerTemp
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND), dimension(kMaxVariables) :: sums, mins, maxes, averages, verticalSumMins, verticalSumMaxes, reductions
real (kind=RKIND), dimension(kMaxVariables) :: sums_tmp, mins_tmp, maxes_tmp, averages_tmp, verticalSumMins_tmp, verticalSumMaxes_tmp
+ real (kind=RKIND), dimension(:,:), allocatable :: enstrophy
+
block => domain % blocklist
dminfo => domain % dminfo
@@ -90,7 +92,6 @@
v => state % v % array
wTop => state % wTop % array
h_edge => state % h_edge % array
- circulation => state % circulation % array
vorticity => state % vorticity % array
ke => state % ke % array
Vor_edge => state % Vor_edge % array
@@ -144,9 +145,9 @@
verticalSumMins(variableIndex) = min(verticalSumMins(variableIndex), verticalSumMins_tmp(variableIndex))
verticalSumMaxes(variableIndex) = max(verticalSumMaxes(variableIndex), verticalSumMaxes_tmp(variableIndex))
- ! circulation
+ ! vorticity
variableIndex = variableIndex + 1
- call ocn_compute_field_local_stats(dminfo, nVertLevels, nVerticesSolve, circulation(:,1:nVerticesSolve), &
+ call ocn_compute_field_local_stats(dminfo, nVertLevels, nVerticesSolve, vorticity(:,1:nVerticesSolve), &
sums_tmp(variableIndex), mins_tmp(variableIndex), maxes_tmp(variableIndex), verticalSumMins_tmp(variableIndex), verticalSumMaxes_tmp(variableIndex))
sums(variableIndex) = sums(variableIndex) + sums_tmp(variableIndex)
mins(variableIndex) = min(mins(variableIndex), mins_tmp(variableIndex))
@@ -154,11 +155,14 @@
verticalSumMins(variableIndex) = min(verticalSumMins(variableIndex), verticalSumMins_tmp(variableIndex))
verticalSumMaxes(variableIndex) = max(verticalSumMaxes(variableIndex), verticalSumMaxes_tmp(variableIndex))
- ! vorticity
+ ! vorticity**2
+ allocate(enstrophy(nVertLevels,nVerticesSolve))
+ enstrophy(:,:)=vorticity(:,1:nVerticesSolve)**2
variableIndex = variableIndex + 1
call ocn_compute_field_area_weighted_local_stats(dminfo, nVertLevels, nVerticesSolve, areaTriangle(1:nVerticesSolve), &
- vorticity(:,1:nVerticesSolve), sums_tmp(variableIndex), mins_tmp(variableIndex), maxes_tmp(variableIndex), &
+ enstrophy(:,:), sums_tmp(variableIndex), mins_tmp(variableIndex), maxes_tmp(variableIndex), &
verticalSumMins_tmp(variableIndex), verticalSumMaxes_tmp(variableIndex))
+ deallocate(enstrophy)
sums(variableIndex) = sums(variableIndex) + sums_tmp(variableIndex)
mins(variableIndex) = min(mins(variableIndex), mins_tmp(variableIndex))
maxes(variableIndex) = max(maxes(variableIndex), maxes_tmp(variableIndex))
@@ -362,7 +366,7 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/(areaEdgeGlobal*nVertLevels)
- ! circulation
+ ! vorticity
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/(nVerticesGlobal*nVertLevels)
Copied: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_gm.F (from rev 1912, trunk/mpas/src/core_ocean/mpas_ocn_gm.F)
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_gm.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_gm.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -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
Copied: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_monthly_forcing.F (from rev 1912, trunk/mpas/src/core_ocean/mpas_ocn_monthly_forcing.F)
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_monthly_forcing.F (rev 0)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_monthly_forcing.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -0,0 +1,190 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! ocn_monthly_forcing
+!
+!> \brief MPAS ocean monthly forcing
+!> \author Doug Jacobsen
+!> \date 04/25/12
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for building the forcing arrays,
+!> if monthly forcing is used.
+!
+!-----------------------------------------------------------------------
+
+module ocn_monthly_forcing
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_timekeeping
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: ocn_build_forcing_arrays, &
+ ocn_monthly_forcing_init
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: monthlyForcingOn !< Flag to turn on/off resotring
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine ocn_build_forcing_arrays
+!
+!> \brief Determines the forcing array used for the monthly forcing.
+!> \author Doug Jacobsen
+!> \date 04/25/12
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the forcing arrays used later in MPAS.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_build_forcing_arrays(timeStamp, grid, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ type(MPAS_Time_type), intent(in) :: timeStamp
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (mesh_type), intent(inout) :: &
+ grid !< Input: grid information
+
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: Error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+
+ real (kind=RKIND), dimension(:,:), pointer :: temperatureRestoreMonthly
+ real (kind=RKIND), dimension(:,:), pointer :: salinityRestoreMonthly
+ real (kind=RKIND), dimension(:,:), pointer :: windStressMonthly
+ real (kind=RKIND), dimension(:), pointer :: temperatureRestore
+ real (kind=RKIND), dimension(:), pointer :: salinityRestore
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ integer :: iCell, iEdge, nCells, nEdges, nMonths, k
+ integer :: iMonth, iMonthP1, iDayInMonth, ierr
+ real (kind=RKIND) :: data, dataP1, weight, weightP1
+
+ err = 0
+
+ if(.not.monthlyForcingOn) return
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nMonths = grid % nMonths
+
+ temperatureRestore => grid % temperatureRestore % array
+ salinityRestore => grid % salinityRestore % array
+ u_src => grid % u_src % array
+
+ temperatureRestoreMonthly => grid % temperatureRestoreMonthly % array
+ salinityRestoreMonthly => grid % salinityRestoreMonthly % array
+ windStressMonthly => grid % windStressMonthly % array
+
+ call mpas_get_time(timeStamp, MM = iMonth, DD = iDayInMonth, ierr = ierr)
+
+ err = ierr
+
+ iMonthP1 = mod(iMonth, nMonths) + 1
+
+ weight = 1.0 - (iDayInMonth-1) / 30.0
+ weightP1 = 1.0 - weight
+
+ do iCell=1,nCells
+ ! Interpolate between iMonth and iMonthP1 records, using iDayInMonth
+ data = temperatureRestoreMonthly(iMonth,iCell)
+ dataP1 = temperatureRestoreMonthly(iMonthP1,iCell)
+ temperatureRestore(iCell) = data * weight + dataP1 * weightP1
+ data = salinityRestoreMonthly(iMonth,iCell)
+ dataP1 = salinityRestoreMonthly(iMonthP1,iCell)
+ salinityRestore(iCell) = data * weight + dataP1 * weightP1
+ end do
+
+ do iEdge=1,nEdges
+ ! Interpolate between iMonth and iMonthP1 records, using iDayInMonth
+ data = windStressMonthly(iMonth,iEdge)
+ dataP1 = windStressMonthly(iMonthP1,iEdge)
+ u_src(1,iEdge) = data * weight + dataP1 * weightP1
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_build_forcing_arrays!}}}
+
+!***********************************************************************
+!
+! routine ocn_monthly_forcing_init
+!
+!> \brief Initializes monthly forcing module
+!> \author Doug Jacobsen
+!> \date 04/25/12
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes the monthly forcing module.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_monthly_forcing_init(err)!{{{
+
+ integer, intent(out) :: err !< Output: error flag
+
+ err = 0
+
+ monthlyForcingOn = .false.
+
+ if(config_use_monthly_forcing) then
+ monthlyForcingOn = .true.
+ end if
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_monthly_forcing_init!}}}
+
+!***********************************************************************
+
+end module ocn_monthly_forcing
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_mpas_core.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_mpas_core.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -11,6 +11,8 @@
use ocn_time_integration
use ocn_tendency
+ use ocn_monthly_forcing
+
use ocn_vel_pressure_grad
use ocn_vel_vadv
use ocn_vel_hmix
@@ -19,6 +21,7 @@
use ocn_tracer_hadv
use ocn_tracer_vadv
use ocn_tracer_hmix
+ use ocn_gm
use ocn_restoring
use ocn_equation_of_state
@@ -27,7 +30,7 @@
use ocn_time_average
- type (io_output_object) :: restart_obj
+ type (io_output_object), save :: restart_obj
integer :: current_outfile_frames
@@ -92,6 +95,9 @@
call mpas_ocn_tracer_advection_init(err_tmp)
err = ior(err,err_tmp)
+ call ocn_monthly_forcing_init(err_tmp)
+ err = ior(err, err_tmp)
+
call mpas_timer_init(domain)
if(err.eq.1) then
@@ -104,6 +110,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 +259,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)
@@ -261,6 +278,16 @@
block % state % time_levs(1) % state % uReconstructMeridional % array &
)
+!TDR
+ call mpas_reconstruct(mesh, mesh % u_src % array, &
+ block % state % time_levs(1) % state % uSrcReconstructX % array, &
+ block % state % time_levs(1) % state % uSrcReconstructY % array, &
+ block % state % time_levs(1) % state % uSrcReconstructZ % array, &
+ block % state % time_levs(1) % state % uSrcReconstructZonal % array, &
+ block % state % time_levs(1) % state % uSrcReconstructMeridional % array &
+ )
+!TDR
+
! initialize velocities and tracers on land to be -1e34
! The reconstructed velocity on land will have values not exactly
! -1e34 due to the interpolation of reconstruction.
@@ -298,6 +325,7 @@
subroutine mpas_core_run(domain, output_obj, output_frame)!{{{
+ use mpas_kind_types
use mpas_grid_types
use mpas_io_output
use mpas_timer
@@ -313,7 +341,7 @@
type (block_type), pointer :: block_ptr
type (MPAS_Time_Type) :: currTime
- character(len=32) :: timeStamp
+ character(len=StrKIND) :: timeStamp
integer :: ierr
! Eventually, dt should be domain specific
@@ -321,7 +349,7 @@
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
- write(0,*) 'Initial time ', timeStamp
+ write(0,*) 'Initial time ', trim(timeStamp)
call ocn_write_output_frame(output_obj, output_frame, domain)
block_ptr => domain % blocklist
@@ -341,7 +369,13 @@
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
- write(0,*) 'Doing timestep ', timeStamp
+ write(0,*) 'Doing timestep ', trim(timeStamp)
+
+ block_ptr => domain % blocklist
+ do while(associated(block_ptr))
+ call ocn_build_forcing_arrays(currTime, block_ptr % mesh, ierr)
+ block_ptr => block_ptr % next
+ end do
call mpas_timer_start("time integration", .false., timeIntTimer)
call mpas_timestep(domain, itimestep, dt, timeStamp)
@@ -451,6 +485,7 @@
subroutine mpas_timestep(domain, itimestep, dt, timeStamp)!{{{
+ use mpas_kind_types
use mpas_grid_types
implicit none
@@ -493,7 +528,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 +542,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 +555,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 +564,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 +575,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 +709,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.
@@ -776,23 +931,26 @@
type (mesh_type), intent(inout) :: mesh
integer :: iEdge, cell1, cell2
- real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4
+ real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4, meshScaling
meshDensity => mesh % meshDensity % array
meshScalingDel2 => mesh % meshScalingDel2 % array
meshScalingDel4 => mesh % meshScalingDel4 % array
+ meshScaling => mesh % meshScaling % array
!
! Compute the scaling factors to be used in the del2 and del4 dissipation
!
meshScalingDel2(:) = 1.0
meshScalingDel4(:) = 1.0
+ meshScaling(:) = 1.0
if (config_h_ScaleWithMesh) then
do iEdge=1,mesh%nEdges
cell1 = mesh % cellsOnEdge % array(1,iEdge)
cell2 = mesh % cellsOnEdge % array(2,iEdge)
- meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/12.0)
- meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/6.0)
+ meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(3.0/4.0) ! goes as dc**3
+ meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(3.0/4.0) ! goes as dc**3
+ meshScaling(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(1.0/4.0)
end do
end if
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_tendency.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_tendency.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -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!}}}
!***********************************************************************
@@ -859,10 +876,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, &
@@ -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
@@ -881,13 +898,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
@@ -898,7 +916,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
@@ -920,51 +938,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
@@ -977,37 +958,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!}}}
@@ -1078,6 +1031,7 @@
end subroutine ocn_fuperp!}}}
+
!***********************************************************************
!
! routine ocn_tendency_init
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -78,7 +78,8 @@
integer :: iCell, k, i, err
type (block_type), pointer :: block
- type (state_type) :: provis
+ type (state_type), target :: provis
+ type (state_type), pointer :: provis_ptr
integer :: rk_step, iEdge, cell1, cell2
@@ -96,10 +97,13 @@
block => domain % blocklist
- call mpas_allocate_state(provis, &
+ call mpas_allocate_state(block, provis, &
block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
- block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels )
+ block % mesh % nVertices, block % mesh % vertexDegree, block % mesh % nVertLevels, block % mesh % nMonths )
+ provis_ptr => provis
+ call mpas_create_state_links(provis_ptr)
+
!
! Initialize time_levs(2) with state at current time
! Initialize first RK state
@@ -142,23 +146,11 @@
! --- update halos for diagnostic variables
call mpas_timer_start("RK4-diagnostic halo update")
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % Vor_edge % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
- if (config_h_mom_eddy_visc4 > 0.0) then
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % divergence % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % vorticity % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nVertices, &
- block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- end if
-
- block => block % next
- end do
+ call mpas_dmpar_exch_halo_field(provis % Vor_edge)
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ call mpas_dmpar_exch_halo_field(provis % divergence)
+ call mpas_dmpar_exch_halo_field(provis % vorticity)
+ end if
call mpas_timer_stop("RK4-diagnostic halo update")
! --- compute tendencies
@@ -190,19 +182,9 @@
! --- update halos for prognostic variables
call mpas_timer_start("RK4-pronostic halo update")
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call mpas_dmpar_exch_halo_field3d_real(domain % dminfo, block % tend % tracers % array(:,:,:), &
- block % tend % num_tracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
+ call mpas_dmpar_exch_halo_field(domain % blocklist % tend % u)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % tend % h)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % tend % tracers)
call mpas_timer_stop("RK4-pronostic halo update")
! --- compute next substep state
@@ -241,6 +223,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 +330,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, &
@@ -351,6 +343,16 @@
block % state % time_levs(2) % state % uReconstructMeridional % array &
)
+!TDR
+ call mpas_reconstruct(block % mesh, block % mesh % u_src % array, &
+ block % state % time_levs(2) % state % uSrcReconstructX % array, &
+ block % state % time_levs(2) % state % uSrcReconstructY % array, &
+ block % state % time_levs(2) % state % uSrcReconstructZ % array, &
+ block % state % time_levs(2) % state % uSrcReconstructZonal % array, &
+ block % state % time_levs(2) % state % uSrcReconstructMeridional % array &
+ )
+!TDR
+
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
block => block % next
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_time_integration_split.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_time_integration_split.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -163,27 +163,11 @@
! --- update halos for diagnostic variables
call mpas_timer_start("se halo diag", .false., timer_halo_diagnostic)
- block => domain % blocklist
- do while (associated(block))
-
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, &
- block % state % time_levs(2) % state % Vor_edge % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
- if (config_h_mom_eddy_visc4 > 0.0) then
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, &
- block % state % time_levs(2) % state % divergence % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, &
- block % state % time_levs(2) % state % vorticity % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nVertices, &
- block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- end if
-
- block => block % next
- end do
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % Vor_edge)
+ if (config_h_mom_eddy_visc4 > 0.0) then
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % divergence)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % vorticity)
+ end if
call mpas_timer_stop("se halo diag", timer_halo_diagnostic)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -272,15 +256,7 @@
end do
call mpas_timer_start("se halo ubcl", .false., timer_halo_ubcl)
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, &
- block % state % time_levs(2) % state % uBcl % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
- block => block % next
- end do
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % uBcl)
call mpas_timer_stop("se halo ubcl", timer_halo_ubcl)
end do ! do j=1,config_n_bcl_iter
@@ -312,6 +288,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
@@ -392,15 +382,7 @@
! boundary update on uBtrNew
call mpas_timer_start("se halo ubtr", .false., timer_halo_ubtr)
- block => domain % blocklist
- do while (associated(block))
-
- call mpas_dmpar_exch_halo_field1d_real(domain % dminfo, &
- block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:), &
- block % mesh % nEdges, block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
- block => block % next
- end do ! block
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle)
call mpas_timer_stop("se halo ubtr", timer_halo_ubtr)
endif ! config_btr_gam1_uWt1>1.0e-12
@@ -459,15 +441,7 @@
! boundary update on SSHnew
call mpas_timer_start("se halo ssh", .false., timer_halo_ssh)
- block => domain % blocklist
- do while (associated(block))
-
- call mpas_dmpar_exch_halo_field1d_real(domain % dminfo, &
- block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
- block % mesh % nCells, block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-
- block => block % next
- end do ! block
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle)
call mpas_timer_stop("se halo ssh", timer_halo_ssh)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -510,14 +484,7 @@
! boundary update on uBtrNew
call mpas_timer_start("se halo ubtr", .false., timer_halo_ubtr)
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field1d_real(domain % dminfo, &
- block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:), &
- block % mesh % nEdges, block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
- block => block % next
- end do ! block
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle)
call mpas_timer_stop("se halo ubtr", timer_halo_ubtr)
end do !do BtrCorIter=1,config_n_btr_cor_iter
@@ -570,15 +537,8 @@
! boundary update on SSHnew
call mpas_timer_start("se halo ssh", .false., timer_halo_ssh)
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field1d_real(domain % dminfo, &
- block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
- block % mesh % nCells, block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-
- block => block % next
- end do ! block
- call mpas_timer_stop("se halo ssh", timer_halo_ssh)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle)
+ call mpas_timer_stop("se halo ssh", timer_halo_ssh)
endif ! config_btr_solve_SSH2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -629,14 +589,7 @@
! boundary update on F
call mpas_timer_start("se halo F", .false., timer_halo_f)
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field1d_real(domain % dminfo, &
- block % state % time_levs(1) % state % FBtr % array(:), &
- block % mesh % nEdges, block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
- block => block % next
- end do ! block
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(1) % state % FBtr)
call mpas_timer_stop("se halo F", timer_halo_f)
@@ -663,9 +616,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 +635,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)
@@ -727,13 +683,7 @@
! update halo for thickness and tracer tendencies
call mpas_timer_start("se halo h", .false., timer_halo_h)
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
+ call mpas_dmpar_exch_halo_field(domain % blocklist % tend % h)
call mpas_timer_stop("se halo h", timer_halo_h)
block => domain % blocklist
@@ -745,13 +695,7 @@
! update halo for thickness and tracer tendencies
call mpas_timer_start("se halo tracers", .false., timer_halo_tracers)
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field3d_real(domain % dminfo, block % tend % tracers % array(:,:,:), &
- block % tend % num_tracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
+ call mpas_dmpar_exch_halo_field(domain % blocklist % tend % tracers)
call mpas_timer_stop("se halo tracers", timer_halo_tracers)
block => domain % blocklist
@@ -801,10 +745,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)
@@ -913,6 +869,7 @@
end if
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+
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, &
@@ -920,6 +877,17 @@
block % state % time_levs(2) % state % uReconstructZonal % array, &
block % state % time_levs(2) % state % uReconstructMeridional % array)
+!TDR
+ call mpas_reconstruct(block % mesh, block % mesh % u_src % array, &
+ block % state % time_levs(2) % state % uSrcReconstructX % array, &
+ block % state % time_levs(2) % state % uSrcReconstructY % array, &
+ block % state % time_levs(2) % state % uSrcReconstructZ % array, &
+ block % state % time_levs(2) % state % uSrcReconstructZonal % array, &
+ block % state % time_levs(2) % state % uSrcReconstructMeridional % array &
+ )
+!TDR
+
+
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -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
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -175,7 +175,7 @@
do k=1,maxLevelEdgeTop(iEdge)
! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
delsq_u = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
- -viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDvEdge
+ -viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDcEdge * sqrt(3.0) ! TDR
! vorticity using </font>
<font color="gray">abla^2 u
r_tmp = dcEdge(iEdge) * delsq_u
@@ -203,7 +203,7 @@
do k=1,maxLevelEdgeTop(iEdge)
u_diffusion = (delsq_divergence(k,cell2) - delsq_divergence(k,cell1)) * invDcEdge &
- -viscVortCoef * (delsq_vorticity(k,vertex2) - delsq_vorticity(k,vertex1) ) * invDvEdge
+ -viscVortCoef * (delsq_vorticity(k,vertex2) - delsq_vorticity(k,vertex1) ) * invDcEdge * sqrt(3.0) ! TDR
tend(k,iEdge) = tend(k,iEdge) - edgeMask(k, iEdge) * u_diffusion * r_tmp
end do
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vmix.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vmix.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -288,13 +288,13 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, k, cell1, cell2, nVertLevels
+ integer :: iEdge, nEdges, k, cell1, cell2, nVertLevels, N
integer, dimension(:), pointer :: maxLevelEdgeTop
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension(:), allocatable :: A, C, uTemp
+ real (kind=RKIND), dimension(:), allocatable :: A, B, C, uTemp
err = 0
@@ -305,43 +305,55 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
cellsOnEdge => grid % cellsOnEdge % array
- allocate(A(nVertLevels),C(nVertLevels),uTemp(nVertLevels))
+ allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),uTemp(nVertLevels))
+ A(1)=0
do iEdge=1,nEdges
- if (maxLevelEdgeTop(iEdge).gt.0) then
+ N=maxLevelEdgeTop(iEdge)
+ if (N.gt.0) then
- ! Compute A(k), C(k) for momentum
- ! mrp 110315 efficiency note: for z-level, could precompute
- ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
- ! h_edge is computed in compute_solve_diag, and is not available yet.
+ ! Compute A(k), B(k), C(k)
+ ! h_edge is computed in compute_solve_diag, and is not available yet,
+ ! so recompute h_edge here.
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
+ do k=1,N
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
end do
- do k=1,maxLevelEdgeTop(iEdge)-1
- A(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &
+ ! A is lower diagonal term
+ do k=2,N
+ A(k) = -2.0*dt*vertViscTopOfEdge(k,iEdge) &
+ / (h_edge(k-1,iEdge) + h_edge(k,iEdge)) &
+ / h_edge(k,iEdge)
+ enddo
+
+ ! C is upper diagonal term
+ do k=1,N-1
+ C(k) = -2.0*dt*vertViscTopOfEdge(k+1,iEdge) &
/ (h_edge(k,iEdge) + h_edge(k+1,iEdge)) &
/ h_edge(k,iEdge)
enddo
- A(maxLevelEdgeTop(iEdge)) = -dt*config_bottom_drag_coeff &
- *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
- C(1) = 1 - A(1)
- do k=2,maxLevelEdgeTop(iEdge)
- C(k) = 1 - A(k) - A(k-1)
+ ! B is diagonal term
+ B(1) = 1 - C(1)
+ do k=2,N-1
+ B(k) = 1 - A(k) - C(k)
enddo
- call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
+ ! Apply bottom drag boundary condition on the viscous term
+ B(N) = 1 - A(N) + dt*config_bottom_drag_coeff &
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
- u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
- u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
+ call tridiagonal_solve(A(2:N),B,C(1:N-1),u(:,iEdge),uTemp,N)
+ u(1:N,iEdge) = uTemp(1:N)
+ u(N+1:nVertLevels,iEdge) = 0.0
+
end if
end do
- deallocate(A,C,uTemp)
+ deallocate(A,B,C,uTemp)
!--------------------------------------------------------------------
@@ -444,8 +456,7 @@
+ fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
enddo
enddo
-!print '(a,50e12.2)', 'fluxVertTop',fluxVertTop(3,1:maxLevelCell(iCell)+1)
-!print '(a,50e12.2)', 'tend_tr ',tend_tr(3,1,1:maxLevelCell(iCell))
+
enddo ! iCell loop
deallocate(fluxVertTop)
!--------------------------------------------------------------------
@@ -509,11 +520,11 @@
!
!-----------------------------------------------------------------
- integer :: iCell, nCells, k, nVertLevels, num_tracers
+ integer :: iCell, nCells, k, nVertLevels, num_tracers, N
integer, dimension(:), pointer :: maxLevelCell
- real (kind=RKIND), dimension(:), allocatable :: A, C
+ real (kind=RKIND), dimension(:), allocatable :: A,B,C
real (kind=RKIND), dimension(:,:), allocatable :: tracersTemp
err = 0
@@ -525,32 +536,39 @@
num_tracers = size(tracers, dim=1)
maxLevelCell => grid % maxLevelCell % array
- allocate(A(nVertLevels),C(nVertLevels), tracersTemp(num_tracers,nVertLevels))
+ allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),tracersTemp(num_tracers,nVertLevels))
do iCell=1,nCells
- ! Compute A(k), C(k) for tracers
- ! mrp 110315 efficiency note: for z-level, could precompute
- ! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
- do k=1,maxLevelCell(iCell)-1
- A(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &
+ ! Compute A(k), B(k), C(k) for tracers
+ N = maxLevelCell(iCell)
+
+ ! A is lower diagonal term
+ A(1)=0
+ do k=2,N
+ A(k) = -2.0*dt*vertDiffTopOfCell(k,iCell) &
+ / (h(k-1,iCell) + h(k,iCell)) / h(k,iCell)
+ enddo
+
+ ! C is upper diagonal term
+ do k=1,N-1
+ C(k) = -2.0*dt*vertDiffTopOfCell(k+1,iCell) &
/ (h(k,iCell) + h(k+1,iCell)) / h(k,iCell)
enddo
+ C(N) = 0.0
- A(maxLevelCell(iCell)) = 0.0
-
- C(1) = 1 - A(1)
- do k=2,maxLevelCell(iCell)
- C(k) = 1 - A(k) - A(k-1)
+ ! B is diagonal term
+ do k=1,N
+ B(k) = 1 - A(k) - C(k)
enddo
- call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &
- tracersTemp, maxLevelCell(iCell), nVertLevels,num_tracers)
+ call tridiagonal_solve_mult(A(2:N),B,C(1:N-1),tracers(:,:,iCell), &
+ tracersTemp, N, nVertLevels,num_tracers)
- tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
- tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
+ tracers(:,1:N,iCell) = tracersTemp(:,1:N)
+ tracers(:,N+1:nVertLevels,iCell) = -1e34
end do
- deallocate(A,C,tracersTemp)
+ deallocate(A,B,C,tracersTemp)
!--------------------------------------------------------------------
Modified: branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-05-15 16:09:17 UTC (rev 1912)
+++ branches/land_ice_projects/implement_core/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-05-15 17:21:01 UTC (rev 1913)
@@ -142,8 +142,8 @@
tracers => s % tracers % array
call mpas_timer_start("eos rich", .false., richEOSTimer)
- call ocn_equation_of_state_rho(s, grid, 0, 'relative', err)
- call ocn_equation_of_state_rho(s, grid, 1, 'relative', err)
+ call ocn_equation_of_state_rho(s, grid, 0,'relative', err)
+ call ocn_equation_of_state_rho(s, grid, 1,'relative', err)
call mpas_timer_stop("eos rich", richEOSTimer)
call ocn_vmix_get_rich_numbers(grid, indexT, indexS, u, h, h_edge, &
</font>
</pre>