<p><b>dwj07@fsu.edu</b> 2012-11-19 13:17:58 -0700 (Mon, 19 Nov 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Merging trunk with cesm coupling branch.<br>
</p><hr noshade><pre><font color="gray">Index: branches/ocean_projects/cesm_coupling/src
===================================================================
--- branches/ocean_projects/cesm_coupling/src 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src 2012-11-19 20:17:58 UTC (rev 2315)
Property changes on: branches/ocean_projects/cesm_coupling/src
___________________________________________________________________
Added: svn:mergeinfo
## -0,0 +1,23 ##
+/branches/atmos_physics/src:1672-1846
+/branches/cam_mpas_nh/src:1260-1270
+/branches/ocean_projects/ale_split_exp/src:1437-1483
+/branches/ocean_projects/ale_vert_coord/src:1225-1383
+/branches/ocean_projects/ale_vert_coord_new/src:1387-1428
+/branches/ocean_projects/gmvar/src:1214-1514,1517-1738
+/branches/ocean_projects/imp_vert_mix_error/src:1847-1887
+/branches/ocean_projects/imp_vert_mix_mrp/src:754-986
+/branches/ocean_projects/monotonic_advection/src:1499-1640
+/branches/ocean_projects/monthly_forcing/src:1810-1867
+/branches/ocean_projects/split_explicit_mrp/src:1134-1138
+/branches/ocean_projects/split_explicit_timestepping/src:1044-1097
+/branches/ocean_projects/vert_adv_mrp/src:704-745
+/branches/ocean_projects/vol_cons_RK_imp_mix/src:1965-1992
+/branches/ocean_projects/zstar_restart_new/src:1762-1770
+/branches/omp_blocks/block_decomp/src:1374-1569
+/branches/omp_blocks/ddt_reorg/src:1301-1414
+/branches/omp_blocks/halo/src:1570-1638
+/branches/omp_blocks/io/src:1639-1787
+/branches/omp_blocks/multiple_blocks/src:1803-2084
+/branches/source_renaming/src:1082-1113
+/branches/time_manager/src:924-962
+/trunk/mpas/src:2147-2314
\ No newline at end of property
Modified: branches/ocean_projects/cesm_coupling/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -6229,10 +6229,12 @@
nearest_distance = current_distance
do i = 1, nEdgesOnCell(current_cell)
iCell = cellsOnCell(i,current_cell)
- d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0_RKIND)
- if (d < nearest_distance) then
- nearest_cell = iCell
- nearest_distance = d
+ if (iCell <= nCells) then
+ d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0_RKIND)
+ if (d < nearest_distance) then
+ nearest_cell = iCell
+ nearest_distance = d
+ end if
end if
end do
end do
@@ -6281,10 +6283,12 @@
end if
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
- d = sphere_distance(latEdge(iEdge), lonEdge(iEdge), target_lat, target_lon, 1.0_RKIND)
- if (d < nearest_distance) then
- nearest_edge = iEdge
- nearest_distance = d
+ if (iEdge <= nEdges) then
+ d = sphere_distance(latEdge(iEdge), lonEdge(iEdge), target_lat, target_lon, 1.0_RKIND)
+ if (d < nearest_distance) then
+ nearest_edge = iEdge
+ nearest_distance = d
+ end if
end if
end do
end do
Index: branches/ocean_projects/cesm_coupling/src/core_ocean
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean 2012-11-19 20:17:58 UTC (rev 2315)
Property changes on: branches/ocean_projects/cesm_coupling/src/core_ocean
___________________________________________________________________
Added: svn:mergeinfo
## -0,0 +1,29 ##
+/branches/atmos_physics/src/core_ocean:1672-1846
+/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-1514,1517-1738
+/branches/ocean_projects/imp_vert_mix_error/src/core_ocean:1847-1887
+/branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean:754-986
+/branches/ocean_projects/leith_mrp/src/core_ocean:2182-2241
+/branches/ocean_projects/monotonic_advection/src/core_ocean:1499-1640
+/branches/ocean_projects/monthly_forcing/src/core_ocean:1810-1867
+/branches/ocean_projects/option3_b4b_test/src/core_ocean:2201-2231
+/branches/ocean_projects/partial_bottom_cells/src/core_ocean:2172-2226
+/branches/ocean_projects/restart_reproducibility/src/core_ocean:2239-2272
+/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/ocean_projects/vol_cons_RK_imp_mix/src/core_ocean:1965-1992
+/branches/ocean_projects/zstar_restart_new/src/core_ocean:1762-1770
+/branches/omp_blocks/block_decomp/src/core_ocean:1374-1569
+/branches/omp_blocks/ddt_reorg/src/core_ocean:1301-1414
+/branches/omp_blocks/halo/src/core_ocean:1570-1638
+/branches/omp_blocks/io/src/core_ocean:1639-1787
+/branches/omp_blocks/multiple_blocks/src/core_ocean:1803-2084
+/branches/omp_blocks/openmp_test/src/core_ocean:2107-2144
+/branches/omp_blocks/openmp_test/src/core_ocean_elements:2161-2201
+/branches/source_renaming/src/core_ocean:1082-1113
+/branches/time_manager/src/core_ocean:924-962
+/trunk/mpas/src/core_ocean:2147-2314
\ No newline at end of property
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/Makefile 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/Makefile 2012-11-19 20:17:58 UTC (rev 2315)
@@ -10,6 +10,7 @@
mpas_ocn_vel_vadv.o \
mpas_ocn_vel_hmix.o \
mpas_ocn_vel_hmix_del2.o \
+ mpas_ocn_vel_hmix_leith.o \
mpas_ocn_vel_hmix_del4.o \
mpas_ocn_vel_forcing.o \
mpas_ocn_vel_forcing_windstress.o \
@@ -87,10 +88,12 @@
mpas_ocn_vel_vadv.o:
-mpas_ocn_vel_hmix.o: mpas_ocn_vel_hmix_del2.o mpas_ocn_vel_hmix_del4.o
+mpas_ocn_vel_hmix.o: mpas_ocn_vel_hmix_del2.o mpas_ocn_vel_hmix_leith.o mpas_ocn_vel_hmix_del4.o
mpas_ocn_vel_hmix_del2.o:
+mpas_ocn_vel_hmix_leith.o:
+
mpas_ocn_vel_hmix_del4.o:
mpas_ocn_vel_forcing.o: mpas_ocn_vel_forcing_windstress.o mpas_ocn_vel_forcing_bottomdrag.o mpas_ocn_vel_forcing_rayleigh.o
@@ -179,6 +182,7 @@
mpas_ocn_vel_vadv.o \
mpas_ocn_vel_hmix.o \
mpas_ocn_vel_hmix_del2.o \
+ mpas_ocn_vel_hmix_leith.o \
mpas_ocn_vel_hmix_del4.o \
mpas_ocn_vel_forcing.o \
mpas_ocn_vel_forcing_windstress.o \
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/Registry 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/Registry 2012-11-19 20:17:58 UTC (rev 2315)
@@ -20,11 +20,12 @@
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 integer io config_pio_num_iotasks 0
-namelist integer io config_pio_stride 1
+namelist integer io config_pio_num_iotasks 0
+namelist integer io config_pio_stride 1
+namelist logical io config_write_output_on_startup true
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 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
@@ -32,6 +33,10 @@
namelist character grid config_pressure_type pressure
namelist real grid config_rho0 1028
namelist logical grid config_enforce_zstar_at_restart false
+namelist character partial_bottom_cells config_alter_ICs_for_pbcs zlevel_pbcs_off
+namelist real partial_bottom_cells config_min_pbc_fraction 0.10
+namelist logical partial_bottom_cells config_check_ssh_consistency true
+namelist logical partial_bottom_cells config_check_zlevel_consistency 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
@@ -59,9 +64,13 @@
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
+namelist logical hmix_leith config_use_leith_del2 false
+namelist real hmix_leith config_leith_parameter 0.0
+namelist real hmix_leith config_leith_dx 0.0
+namelist real hmix_leith config_leith_visc2_max 1000000.0
namelist character vmix config_vert_visc_type const
namelist character vmix config_vert_diff_type const
-namelist logical vmix config_implicit_vertical_mix .true.
+namelist logical vmix config_implicit_vertical_mix true
namelist real vmix config_convective_visc 1.0
namelist real vmix config_convective_diff 1.0
namelist real vmix config_bottom_drag_coeff 1.0e-3
@@ -105,7 +114,7 @@
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
dim nVertLevelsP1 nVertLevels+1
-dim nMonths 1
+dim nMonths nMonths
%
% var persistence type name_in_file ( dims ) time_levs iro- name_in_code struct super-array array_class
@@ -163,7 +172,7 @@
var persistent real kiteAreasOnVertex ( vertexDegree nVertices ) 0 iro kiteAreasOnVertex mesh - -
var persistent real fEdge ( nEdges ) 0 iro fEdge mesh - -
var persistent real fVertex ( nVertices ) 0 iro fVertex mesh - -
-var persistent real h_s ( nCells ) 0 iro h_s mesh - -
+var persistent real bottomDepth ( nCells ) 0 iro bottomDepth mesh - -
% Space needed for advection
var persistent real deriv_two ( maxEdges2 TWO nEdges ) 0 - deriv_two mesh - -
@@ -193,26 +202,26 @@
var persistent integer maxLevelEdgeBot ( nEdges ) 0 - maxLevelEdgeBot mesh - -
var persistent integer maxLevelVertexTop ( nVertices ) 0 - maxLevelVertexTop mesh - -
var persistent integer maxLevelVertexBot ( nVertices ) 0 - maxLevelVertexBot mesh - -
-var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
-var persistent real referenceBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - referenceBottomDepthTopOfCell mesh - -
+var persistent real refBottomDepth ( nVertLevels ) 0 iro refBottomDepth mesh - -
+var persistent real refBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - refBottomDepthTopOfCell 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 - -
-var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
+% Boundary conditions and masks
+var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 - boundaryEdge mesh - -
+var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 - boundaryVertex mesh - -
var persistent integer boundaryCell ( nVertLevels nCells ) 0 - boundaryCell mesh - -
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 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 - windStressMonthly mesh - -
-%var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 - temperatureRestoreMonthly mesh - -
-%var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 - salinityRestoreMonthly mesh - -
+% Forcing variables.
+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 windStressMonthly ( nMonths nEdges ) 0 ir windStressMonthly mesh - -
+var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 ir temperatureRestoreMonthly mesh - -
+var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 ir 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 - -
@@ -231,7 +240,7 @@
var persistent real tend_tracer1 ( nVertLevels nCells Time ) 1 - tracer1 tend tracers testing
% state variables for Split Explicit timesplitting
-var persistent real uBtr ( nEdges Time ) 2 - uBtr state - -
+var persistent real uBtr ( nEdges Time ) 2 r uBtr state - -
var persistent real ssh ( nCells Time ) 2 o ssh state - -
var persistent real uBtrSubcycle ( nEdges Time ) 2 - uBtrSubcycle state - -
var persistent real sshSubcycle ( nCells Time ) 2 - sshSubcycle state - -
@@ -277,6 +286,7 @@
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
+var persistent real viscosity ( nVertLevels nEdges Time ) 2 o viscosity state - -
% Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
@@ -308,3 +318,8 @@
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 - -
+
+% Sign fields, for openmp and bit reproducibility without branching statements.
+var persistent integer edgeSignOnCell ( maxEdges nCells ) 0 - edgeSignOnCell mesh - -
+var persistent integer edgeSignOnVertex ( maxEdges nVertices ) 0 - edgeSignOnVertex mesh - -
+var persistent integer kiteIndexOnCell ( maxEdges nCells ) 0 - kiteIndexOnCell mesh - -
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_equation_of_state_jm.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -96,7 +96,7 @@
real (kind=RKIND), dimension(:), pointer :: &
- referenceBottomDepth, pRefEOS
+ refBottomDepth, pRefEOS
real (kind=RKIND), dimension(:,:), intent(inout) :: &
rho
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
@@ -197,7 +197,7 @@
nCells = grid % nCells
maxLevelCell => grid % maxLevelCell % array
nVertLevels = grid % nVertLevels
- referenceBottomDepth => grid % referenceBottomDepth % array
+ refBottomDepth => grid % refBottomDepth % array
! Jackett and McDougall
@@ -214,14 +214,14 @@
allocate(pRefEOS(nVertLevels),p(nVertLevels),p2(nVertLevels))
! This could be put in the init routine.
- ! Note I am using referenceBottomDepth, so pressure on top level does
+ ! Note I am using refBottomDepth, so pressure on top level does
! not include SSH contribution. I am not sure if that matters, but
! POP does it the same way.
- depth = 0.5*referenceBottomDepth(1)
+ depth = 0.5*refBottomDepth(1)
pRefEOS(1) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ 0.100766*depth + 2.28405e-7*depth**2
do k = 2,nVertLevels
- depth = 0.5*(referenceBottomDepth(k)+referenceBottomDepth(k-1))
+ depth = 0.5*(refBottomDepth(k)+refBottomDepth(k-1))
pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ 0.100766*depth + 2.28405e-7*depth**2
enddo
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_monthly_forcing.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_monthly_forcing.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_monthly_forcing.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -176,6 +176,8 @@
if(config_use_monthly_forcing) then
monthlyForcingOn = .true.
+
+ write (0,'(a)') " Monthly forcing is on. Make sure monthly forcing variables include iro in Registry, and are in your initial condition or restart file."
end if
!--------------------------------------------------------------------
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_mpas_core.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_mpas_core.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -106,37 +106,37 @@
if (.not. config_do_restart) call setup_sw_test_case(domain)
+ if (config_vert_grid_type.ne.'isopycnal') call ocn_init_vert_coord(domain)
+
call ocn_compute_max_level(domain)
- 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)
+ if (.not.config_do_restart) call ocn_init_split_timestep(domain)
- print *, ' Vertical grid type is: ',config_vert_grid_type
+ write (0,'(a,a10)') ' Vertical grid type is: ',config_vert_grid_type
if (config_vert_grid_type.ne.'isopycnal'.and. &
config_vert_grid_type.ne.'zlevel'.and. &
config_vert_grid_type.ne.'zstar1'.and. &
config_vert_grid_type.ne.'zstar'.and. &
config_vert_grid_type.ne.'zstarWeights') then
- print *, ' Incorrect choice of config_vert_grid_type.'
+ write (0,*) ' Incorrect choice of config_vert_grid_type.'
call mpas_dmpar_abort(dminfo)
endif
- print *, ' Pressure type is: ',config_pressure_type
+ write (0,'(a,a10)') ' Pressure type is: ',config_pressure_type
if (config_pressure_type.ne.'pressure'.and. &
config_pressure_type.ne.'MontgomeryPotential') then
- print *, ' Incorrect choice of config_pressure_type.'
+ write (0,*) ' Incorrect choice of config_pressure_type.'
call mpas_dmpar_abort(dminfo)
endif
if (config_filter_btr_mode.and. &
config_vert_grid_type.ne.'zlevel')then
- print *, 'filter_btr_mode has only been tested with'// &
+ write (0,*) 'filter_btr_mode has only been tested with'// &
' config_vert_grid_type=zlevel.'
call mpas_dmpar_abort(dminfo)
endif
@@ -249,6 +249,7 @@
integer :: i, iEdge, iCell, k
integer :: err1
+ call ocn_setup_sign_and_index_fields(mesh)
call ocn_initialize_advection_rk(mesh, err)
call mpas_ocn_tracer_advection_coefficients(mesh, err1)
err = ior(err, err1)
@@ -264,8 +265,6 @@
= 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)
call mpas_rbf_interp_initialize(mesh)
@@ -351,7 +350,9 @@
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
write(0,*) 'Initial time ', trim(timeStamp)
- call ocn_write_output_frame(output_obj, output_frame, domain)
+ if (config_write_output_on_startup) then
+ call ocn_write_output_frame(output_obj, output_frame, domain)
+ endif
block_ptr => domain % blocklist
do while(associated(block_ptr))
@@ -390,7 +391,8 @@
if (mpas_is_alarm_ringing(clock, outputAlarmID, ierr=ierr)) then
call mpas_reset_clock_alarm(clock, outputAlarmID, ierr=ierr)
- ! output_frame will always be > 1 here unless it was reset after the maximum number of frames per outfile was reached
+ ! output_frame will always be > 1 here unless it was reset after the
+ ! maximum number of frames per outfile was reached.
if(output_frame == 1) then
call mpas_output_state_finalize(output_obj, domain % dminfo)
call mpas_output_state_init(output_obj, domain, "OUTPUT", trim(timeStamp))
@@ -531,8 +533,9 @@
end subroutine mpas_timestep!}}}
- subroutine ocn_init_z_level(domain)!{{{
- ! Initialize zlevel-type variables
+ subroutine ocn_init_vert_coord(domain)!{{{
+ ! Initialize zlevel-type variables and adjust initial conditions for
+ ! partial bottom cells.
use mpas_grid_types
use mpas_configure
@@ -540,43 +543,55 @@
implicit none
type (domain_type), intent(inout) :: domain
+ type (dm_info) :: dminfo
- integer :: i, iCell, iEdge, iVertex, k
+ integer :: i, iCell, iEdge, iVertex, k, nCells, num_tracers
type (block_type), pointer :: block
integer :: iTracer, cell, cell1, cell2
- real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, &
- referenceBottomDepthTopOfCell, zstarWeight, hZLevel
+ real (kind=RKIND) :: uhSum, hSum, hEdge1, zMidPBC
+
+ integer, dimension(:), pointer :: maxLevelCell
+ real (kind=RKIND), dimension(:), pointer :: refBottomDepth, &
+ refBottomDepthTopOfCell, zstarWeight, hZLevel, bottomDepth
+ real (kind=RKIND), dimension(:), allocatable :: minBottomDepth, minBottomDepthMid, zMidZLevel
real (kind=RKIND), dimension(:,:), pointer :: h
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer :: nVertLevels
+ logical :: consistentSSH
! 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
- referenceBottomDepthTopOfCell => block % mesh % referenceBottomDepthTopOfCell % array
+ tracers => block % state % time_levs(1) % state % tracers % array
+ refBottomDepth => block % mesh % refBottomDepth % array
+ refBottomDepthTopOfCell => block % mesh % refBottomDepthTopOfCell % array
+ bottomDepth => block % mesh % bottomDepth % array
zstarWeight => block % mesh % zstarWeight % array
hZLevel => block % mesh % hZLevel % array
+ maxLevelCell => block % mesh % maxLevelCell % array
+
+ nCells = block % mesh % nCells
nVertLevels = block % mesh % nVertLevels
+ num_tracers = size(tracers, dim=1)
! mrp 120208 right now hZLevel is in the grid.nc file.
- ! We would like to transition to using referenceBottomDepth
+ ! We would like to transition to using refBottomDepth
! 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) = hZLevel(1)
+ refBottomDepth(1) = hZLevel(1)
do k = 2,nVertLevels
- referenceBottomDepth(k) = referenceBottomDepth(k-1) + hZLevel(k)
+ refBottomDepth(k) = refBottomDepth(k-1) + hZLevel(k)
end do
! TopOfCell needed where zero depth for the very top may be referenced.
- referenceBottomDepthTopOfCell(1) = 0.0
+ refBottomDepthTopOfCell(1) = 0.0
do k = 1,nVertLevels
- referenceBottomDepthTopOfCell(k+1) = referenceBottomDepth(k)
+ refBottomDepthTopOfCell(k+1) = refBottomDepth(k)
end do
! Initialization of zstarWeights. This determines how SSH perturbations
@@ -588,9 +603,7 @@
elseif (config_vert_grid_type.eq.'zstar') then
- do k = 1,nVertLevels
- zstarWeight(k) = hZLevel(k)
- enddo
+ zstarWeight = 1.0
elseif (config_vert_grid_type.eq.'zstarWeights') then
@@ -605,10 +618,125 @@
endif
+ ! Initial condition files (ocean.nc, produced by basin) include a realistic
+ ! bottomDepth variable and h,T,S variables for full thickness cells.
+ ! If running with pbcs, set config_alter_ICs_for_pbc='zlevel_pbcs_on'. Then thin pbc cells
+ ! will be changed, and h,T,S will be altered to match the pbcs.
+ ! If running without pbcs, set config_alter_ICs_for_pbc='zlevel_pbcs_off'. Then
+ ! bottomDepth will be altered so it is full cells everywhere.
+ ! If your input file does not include bottomDepth, the false option will
+ ! initialize bottomDepth correctly for a non-pbc run.
+
+
+ if (.not.config_do_restart) then
+
+ if (config_alter_ICs_for_pbcs.eq.'zlevel_pbcs_on') then
+
+ write (0,'(a)') ' Altering bottomDepth to avoid very thin cells.'
+ write (0,'(a)') ' Altering h and tracer initial conditions to conform with partial bottom cells.'
+
+ allocate(minBottomDepth(nVertLevels),minBottomDepthMid(nVertLevels),zMidZLevel(nVertLevels))
+
+ ! min_pbc_fraction restricts pbcs from being too small.
+ ! A typical value is 10%, so pbcs must occupy at least 10% of the cell thickness.
+ ! If min_pbc_fraction = 0.0, bottomDepth gives the actual depth for that cell.
+ ! If min_pbc_fraction = 1.0, bottomDepth reverts to discrete z-level depths, same
+ ! as partial_bottom_cells = .false.
+
+ do k=1,nVertLevels
+ minBottomDepth(k) = refBottomDepth(k) - (1.0-config_min_pbc_fraction)*hZLevel(k)
+ minBottomDepthMid(k) = 0.5*(minBottomDepth(k) + refBottomDepthTopOfCell(k))
+ zMidZLevel(k) = - 0.5*(refBottomDepth(k) + refBottomDepthTopOfCell(k))
+ enddo
+
+ do iCell=1,nCells
+ k = maxLevelCell(iCell)
+
+ if (bottomDepth(iCell).lt.minBottomDepthMid(k)) then
+ ! Round up to cell above
+ maxLevelCell(iCell) = maxLevelCell(iCell) - 1
+ bottomDepth(iCell) = refBottomDepth(maxLevelCell(iCell))
+ elseif (bottomDepth(iCell).lt.minBottomDepth(k)) then
+ ! Round down cell to the min_pbc_fraction.
+ bottomDepth(iCell) = minBottomDepth(k)
+ endif
+ k = maxLevelCell(iCell)
+
+ ! Alter thickness of bottom level to account for PBC
+ h(k,iCell) = bottomDepth(iCell) - refBottomDepthTopOfCell(k)
+
+ ! Linearly interpolate the initial T&S for new location of bottom cell for PBCs
+ zMidPBC = -0.5*(bottomDepth(iCell) + refBottomDepthTopOfCell(k))
+
+ do iTracer=1,num_tracers
+ tracers(iTracer,k,iCell) = tracers(iTracer,k,iCell) &
+ + (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell)) &
+ /(zMidZLevel(k-1)-zMidZLevel(k)) &
+ *(zMidPBC - zMidZLevel(k))
+ enddo
+
+ enddo
+
+ deallocate(minBottomDepth,zMidZLevel)
+
+ elseif (config_alter_ICs_for_pbcs.eq.'zlevel_pbcs_off') then
+
+ do iCell = 1,nCells
+ bottomDepth(iCell) = refBottomDepth(maxLevelCell(iCell))
+ enddo
+
+ elseif (config_alter_ICs_for_pbcs.eq.'off') then
+ ! No action taken. This is for isopycnal or sigma coordinates,
+ ! or if ICs were already altered upon start-up.
+
+ else
+
+ write (0,*) ' Incorrect choice of config_alter_ICs_for_pbcs.'
+ call mpas_dmpar_abort(dminfo)
+
+ endif
+ endif
+
+ if (config_check_ssh_consistency) then
+ consistentSSH = .true.
+ do iCell = 1,nCells
+ ! Check if abs(ssh)>2m. If so, print warning.
+ if (abs(sum(h(1:maxLevelCell(iCell),iCell))-bottomDepth(iCell))>2.0) then
+ consistentSSH = .false.
+#ifdef MPAS_DEBUG
+ write (0,'(a)') ' Warning: abs(sum(h)-bottomDepth)>2m. Most likely, initial h does not match bottomDepth.'
+ write (0,*) ' iCell, K=maxLevelCell(iCell), bottomDepth(iCell),sum(h),bottomDepth,hZLevel(K),h(K): ', &
+ iCell, maxLevelCell(iCell), bottomDepth(iCell),sum(h(1:maxLevelCell(iCell),iCell)),bottomDepth(iCell), &
+ hZLevel(maxLevelCell(iCell)), h(maxLevelCell(iCell),iCell)
+#endif
+ endif
+ enddo
+
+ if (.not. consistentSSH) then
+ write(0,*) 'Warning: SSH is not consistent. Most likely, initial h does not match bottomDepth.'
+ end if
+ endif
+
+ if (config_check_zlevel_consistency) then
+ do iCell = 1,nCells
+ ! Check that bottomDepth and maxLevelCell match. Some older grids do not have the bottomDepth variable.
+ if (bottomDepth(iCell) > refBottomDepth(maxLevelCell(iCell)).or. &
+ bottomDepth(iCell) < refBottomDepthTopOfCell(maxLevelCell(iCell))) then
+ write (0,'(a)') ' fatal error: bottomDepth and maxLevelCell do not match:'
+ write (0,'(a,2i5,10f10.2)') ' iCell, maxLevelCell(iCell), bottomDepth(iCell): ', &
+ iCell, maxLevelCell(iCell), bottomDepth(iCell)
+ write (0,'(a,10f10.2)') ' refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell)): ', &
+ refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell))
+ call mpas_dmpar_abort(dminfo)
+ endif
+
+ enddo
+ endif
+
block => block % next
end do
- end subroutine ocn_init_z_level!}}}
+ end subroutine ocn_init_vert_coord!}}}
subroutine ocn_init_split_timestep(domain)!{{{
! Initialize splitting variables
@@ -625,7 +753,7 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:), pointer :: refBottomDepth
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -635,7 +763,7 @@
do while (associated(block))
h => block % state % time_levs(1) % state % h % array
- referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ refBottomDepth => block % mesh % refBottomDepth % array
nVertLevels = block % mesh % nVertLevels
! Compute barotropic velocity at first timestep
@@ -651,7 +779,7 @@
if (config_filter_btr_mode) then
do iCell=1,block % mesh % nCells
block % state % time_levs(1) % state % h % array(1,iCell) &
- = block % mesh % referenceBottomDepth % array(1)
+ = block % mesh % refBottomDepth % array(1)
enddo
endif
@@ -735,7 +863,7 @@
real (kind=RKIND) :: hSum, sumZstarWeights
real (kind=RKIND), dimension(:), pointer :: hZLevel, zstarWeight, &
- referenceBottomDepth
+ refBottomDepth
real (kind=RKIND), dimension(:,:), pointer :: h
! Initialize z-level grid variables from h, read in from input file.
@@ -747,7 +875,7 @@
hZLevel => block % mesh % hZLevel % array
maxLevelCell => block % mesh % maxLevelCell % array
zstarWeight => block % mesh % zstarWeight % array
- referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ refBottomDepth => block % mesh % refBottomDepth % array
do iCell=1,block % mesh % nCells
! Compute the total column thickness, hSum, and the sum of zstar weights.
@@ -762,7 +890,7 @@
! where zeta is SSH and W_k are weights
do k = 1,maxLevelCell(iCell)
h(k,iCell) = hZLevel(k) &
- + (hSum - referenceBottomDepth(maxLevelCell(iCell))) &
+ + (hSum - refBottomDepth(maxLevelCell(iCell))) &
* zstarWeight(k)/sumZstarWeights
enddo
@@ -787,10 +915,6 @@
integer :: i, iCell, iEdge, iVertex, k
type (block_type), pointer :: block
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
- real (kind=RKIND) :: centerx, centery
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
integer, dimension(:), pointer :: &
@@ -960,6 +1084,72 @@
end subroutine ocn_compute_mesh_scaling!}}}
+ subroutine ocn_setup_sign_and_index_fields(mesh)!{{{
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnCell, edgesOnVertex, cellsOnVertex, cellsOnEdge, verticesOnCell, verticesOnEdge
+ integer, dimension(:,:), pointer :: edgeSignOnCell, edgeSignOnVertex, kiteIndexOnCell
+
+ integer :: nCells, nEdges, nVertices, vertexDegree
+ integer :: iCell, iEdge, iVertex, i, j, k
+
+ nCells = mesh % nCells
+ nEdges = mesh % nEdges
+ nVertices = mesh % nVertices
+ vertexDegree = mesh % vertexDegree
+
+ nEdgesOnCell => mesh % nEdgesOnCell % array
+ edgesOnCell => mesh % edgeSOnCell % array
+ edgesOnVertex => mesh % edgesOnVertex % array
+ cellsOnVertex => mesh % cellsOnVertex % array
+ cellsOnEdge => mesh % cellsOnEdge % array
+ verticesOnCell => mesh % verticesOnCell % array
+ verticesOnEdge => mesh % verticesOnEdge % array
+ edgeSignOnCell => mesh % edgeSignOnCell % array
+ edgeSignOnVertex => mesh % edgeSignOnVertex % array
+ kiteIndexOnCell => mesh % kiteIndexOnCell % array
+
+ edgeSignOnCell = 0.0_RKIND
+ edgeSignOnVertex = 0.0_RKIND
+ kiteIndexOnCell = 0.0_RKIND
+
+ do iCell = 1, nCells
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+
+ ! Vector points from cell 1 to cell 2
+ if(iCell == cellsOnEdge(1, iEdge)) then
+ edgeSignOnCell(i, iCell) = -1
+ else
+ edgeSignOnCell(i, iCell) = 1
+ end if
+
+ do j = 1, vertexDegree
+ if(cellsOnVertex(j, iVertex) == iCell) then
+ kiteIndexOnCell(i, iCell) = j
+ end if
+ end do
+ end do
+ end do
+
+ do iVertex = 1, nVertices
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+
+ ! Vector points from vertex 1 to vertex 2
+ if(iVertex == verticesOnEdge(1, iEdge)) then
+ edgeSignOnVertex(i, iVertex) = -1
+ else
+ edgeSignOnVertex(i, iVertex) = 1
+ end if
+ end do
+ end do
+
+ end subroutine ocn_setup_sign_and_index_fields!}}}
+
end module mpas_core
! vim: foldmethod=marker
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tendency.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tendency.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -172,7 +172,7 @@
real (kind=RKIND), dimension(:,:), pointer :: &
h_edge, h, u, rho, zMid, pressure, &
- tend_u, circulation, vorticity, ke, ke_edge, Vor_edge, &
+ tend_u, circulation, vorticity, viscosity, ke, ke_edge, Vor_edge, &
MontPot, wTop, divergence, vertViscTopOfEdge
real (kind=RKIND), dimension(:,:), pointer :: u_src
@@ -186,6 +186,7 @@
wTop => s % wTop % array
zMid => s % zMid % array
h_edge => s % h_edge % array
+ viscosity => s % viscosity % array
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
@@ -237,7 +238,7 @@
! strictly only valid for config_h_mom_eddy_visc2 == constant
!
call mpas_timer_start("hmix", .false., velHmixTimer)
- call ocn_vel_hmix_tend(grid, divergence, vorticity, tend_u, err)
+ call ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend_u, err)
call mpas_timer_stop("hmix", velHmixTimer)
!
@@ -404,15 +405,14 @@
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- verticesOnEdge, edgesOnEdge, edgesOnVertex,boundaryCell
+ verticesOnEdge, edgesOnEdge, edgesOnVertex,boundaryCell, kiteIndexOnCell, verticesOnCell, edgeSignOnVertex, edgeSignOnCell, edgesOnCell
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2, coef_3rd_order, r_tmp, invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, invLength, h_vertex
real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, &
- referenceBottomDepth, ssh
+ bottomDepth, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, ssh
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -435,11 +435,11 @@
kev => s % kev % array
kevc => s % kevc % array
ke_edge => s % ke_edge % array
- Vor_edge => s % Vor_edge % array
- Vor_vertex => s % Vor_vertex % array
- Vor_cell => s % Vor_cell % array
- gradVor_n => s % gradVor_n % array
- gradVor_t => s % gradVor_t % array
+ Vor_edge => s % Vor_edge % array
+ Vor_vertex => s % Vor_vertex % array
+ Vor_cell => s % Vor_cell % array
+ gradVor_n => s % gradVor_n % array
+ gradVor_t => s % gradVor_t % array
rho => s % rho % array
MontPot => s % MontPot % array
pressure => s % pressure % array
@@ -454,20 +454,22 @@
verticesOnEdge => grid % verticesOnEdge % array
nEdgesOnCell => grid % nEdgesOnCell % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
edgesOnEdge => grid % edgesOnEdge % array
edgesOnVertex => grid % edgesOnVertex % array
dcEdge => grid % dcEdge % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
+ bottomDepth => grid % bottomDepth % array
fVertex => grid % fVertex % array
- referenceBottomDepth => grid % referenceBottomDepth % array
deriv_two => grid % deriv_two % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
maxLevelVertexBot => grid % maxLevelVertexBot % array
+ kiteIndexOnCell => grid % kiteIndexOnCell % array
+ verticesOnCell => grid % verticesOnCell % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -477,7 +479,10 @@
boundaryCell => grid % boundaryCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
!
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
@@ -578,40 +583,33 @@
divergence(:,:) = 0.0
ke(:,:) = 0.0
v(:,:) = 0.0
- do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+ do iVertex = 1, nVertices
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ do k = 1, maxLevelVertexBot(iVertex)
+ r_tmp = dcEdge(iEdge) * u(k, iEdge)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ circulation(k, iVertex) = circulation(k, iVertex) + edgeSignOnVertex(i, iVertex) * r_tmp
+ vorticity(k, iVertex) = vorticity(k, iVertex) + edgeSignOnVertex(i, iVertex) * r_tmp * invAreaTri1
+ end do
+ end do
+ end do
- invAreaTri1 = 1.0 / areaTriangle(vertex1)
- invAreaTri2 = 1.0 / areaTriangle(vertex2)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ r_tmp = dvEdge(iEdge) * u(k, iEdge) * invAreaCell1
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(cell1), 1.0_RKIND)
- invAreaCell2 = 1.0 / max(areaCell(cell2), 1.0_RKIND)
-
- do k=1,maxLevelEdgeBot(iEdge)
- ! Compute circulation and relative vorticity at each vertex
- r_tmp = dcEdge(iEdge) * u(k,iEdge)
- circulation(k,vertex1) = circulation(k,vertex1) - r_tmp
- circulation(k,vertex2) = circulation(k,vertex2) + r_tmp
-
- vorticity(k, vertex1) = vorticity(k, vertex1) - r_tmp * invAreaTri1
- vorticity(k, vertex2) = vorticity(k, vertex2) + r_tmp * invAreaTri2
-
- ! Compute the divergence at each cell center
- r_tmp = dvEdge(iEdge) * u(k, iEdge)
- divergence(k,cell1) = divergence(k,cell1) + r_tmp * invAreaCell1
- divergence(k,cell2) = divergence(k,cell2) - r_tmp * invAreaCell2
-
- ! Compute kinetic energy in each cell
- r_tmp = r_tmp * dcEdge(iEdge) * u(k,iEdge)
- ke(k,cell1) = ke(k,cell1) + 0.25 * r_tmp * invAreaCell1
- ke(k,cell2) = ke(k,cell2) + 0.25 * r_tmp * invAreaCell2
+ divergence(k, iCell) = divergence(k, iCell) - edgeSignOnCell(i, iCell) * r_tmp
+ ke(k, iCell) = ke(k, iCell) + 0.25 * r_tmp * dcEdge(iEdge) * u(k,iEdge)
+ end do
end do
+ end do
+ do iEdge=1,nEdges
! Compute v (tangential) velocities
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
@@ -628,29 +626,27 @@
! Compute kinetic energy in each vertex
!
kev(:,:) = 0.0; kevc(:,:) = 0.0
- do iEdge=1,nEdges*ke_vertex_flag
- do k=1,nVertLevels
- r_tmp = dcEdge(iEdge) * dvEdge(iEdge) * u(k, iEdge)**2
- kev(k,verticesOnEdge(1,iEdge)) = kev(k,verticesOnEdge(1,iEdge)) + r_tmp
- kev(k,verticesOnEdge(2,iEdge)) = kev(k,verticesOnEdge(2,iEdge)) + r_tmp
- end do
+ do iVertex = 1, nVertices*ke_vertex_flag
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ r_tmp = dcEdge(iEdge) * dvEdge(iEdge) * 0.25 / areaTriangle(iVertex)
+ do k = 1, nVertLevels
+ kev(k, iVertex) = kev(k, iVertex) + r_tmp * u(k, iEdge)**2
+ end do
+ end do
end do
- do iVertex = 1,nVertices*ke_vertex_flag
- do k=1,nVertLevels
- kev(k,iVertex) = kev(k,iVertex) / areaTriangle(iVertex) * 0.25
- enddo
- enddo
- do iVertex = 1, nVertices*ke_vertex_flag
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0_RKIND)
- do k=1,nVertLevels
- kevc(k,iCell) = kevc(k,iCell) + kiteAreasOnVertex(i, iVertex) * kev(k, iVertex) * invAreaCell1
- enddo
- enddo
- enddo
+ do iCell = 1, nCells*ke_vertex_flag
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ j = kiteIndexOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+ do k = 1, nVertLevels
+ kevc(k, iCell) = kevc(k, iCell) + kiteAreasOnVertex(j, iVertex) * kev(k, iVertex) * invAreaCell1
+ end do
+ end do
+ end do
+
!
! Compute kinetic energy in each cell by blending ke and kevc
!
@@ -692,30 +688,26 @@
Vor_cell(:,:) = 0.0
Vor_edge(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- iEdge = edgesOnVertex(i,iVertex)
+ do iEdge = 1, nEdges
+ vertex1 = verticesOnEdge(1, iEdge)
+ vertex2 = verticesOnEdge(2, iEdge)
+ do k = 1, maxLevelEdgeBot(iEdge)
+ Vor_edge(k, iEdge) = 0.5 * (Vor_vertex(k, vertex1) + Vor_vertex(k, vertex2))
+ end do
+ end do
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0_RKIND)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
- ! Compute pv at cell centers
- ! ( this computes Vor_cell for all real cells and distance-1 ghost cells )
- do k = 1,maxLevelCell(iCell)
- Vor_cell(k,iCell) = Vor_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * Vor_vertex(k, iVertex) * invAreaCell1
- enddo
+ do i = 1, nEdgesOnCell(iCell)
+ j = kiteIndexOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ Vor_cell(k, iCell) = Vor_cell(k, iCell) + kiteAreasOnVertex(j, iVertex) * Vor_vertex(k, iVertex) * invAreaCell1
+ end do
+ end do
+ end do
- ! Compute pv at the edges
- ! ( this computes Vor_edge at all edges bounding real cells )
- do k=1,maxLevelEdgeBot(iEdge)
- Vor_edge(k,iEdge) = Vor_edge(k,iEdge) + 0.5 * Vor_vertex(k,iVertex)
- enddo
- enddo
- enddo
-
-! gradVor_n(:,:) = 0.0
-! gradVor_t(:,:) = 0.0
do iEdge = 1,nEdges
cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)
@@ -779,9 +771,9 @@
pTop(1) = 0.0
! For isopycnal mode, p is the Montgomery Potential.
! At top layer it is g*SSH, where SSH may be off by a
- ! constant (ie, h_s can be relative to top or bottom)
+ ! constant (ie, bottomDepth can be relative to top or bottom)
MontPot(1,iCell) = gravity &
- * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
+ * (bottomDepth(iCell) + sum(h(1:nVertLevels,iCell)))
do k=2,nVertLevels
pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
@@ -810,10 +802,10 @@
! Compute zMid, the z-coordinate of the middle of the layer.
! This is used for the rho g grad z momentum term.
- ! Note the negative sign, since referenceBottomDepth is positive
+ ! Note the negative sign, since bottomDepth is positive
! and z-coordinates are negative below the surface.
k = maxLevelCell(iCell)
- zMid(k:nVertLevels,iCell) = -referenceBottomDepth(k) + 0.5*h(k,iCell)
+ zMid(k:nVertLevels,iCell) = -bottomDepth(iCell) + 0.5*h(k,iCell)
do k=maxLevelCell(iCell)-1, 1, -1
zMid(k,iCell) = zMid(k+1,iCell) &
@@ -830,13 +822,11 @@
!
do iCell=1,nCells
! Start at the bottom where we know the depth, and go up.
- ! The bottom depth for this cell is
- ! referenceBottomDepth(maxLevelCell(iCell)).
- ! Note the negative sign, since referenceBottomDepth is positive
+ ! The bottom depth for this cell is bottomDepth(iCell).
+ ! Note the negative sign, since bottomDepth is positive
! and z-coordinates are negative below the surface.
- ssh(iCell) = -referenceBottomDepth(maxLevelCell(iCell)) &
- + sum(h(1:maxLevelCell(iCell),iCell))
+ ssh(iCell) = - bottomDepth(iCell) + sum(h(1:maxLevelCell(iCell),iCell))
end do
@@ -864,39 +854,68 @@
!> This routine computes the vertical velocity in the top layer for the ocean
!
!-----------------------------------------------------------------------
- subroutine ocn_wtop(s1,s2, grid)!{{{
- implicit none
+ subroutine ocn_wtop(grid,h,h_edge,u,wTop, err)!{{{
- type (state_type), intent(inout) :: s1 !< Input/Output: State 1 information
- type (state_type), intent(inout) :: s2 !< Input/Output: State 2 information
- type (mesh_type), intent(in) :: grid !< Input: Grid information
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h_edge !< Input: h interpolated to an edge
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u !< Input: velocity
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: &
+ wTop !< Output: vertical transport at top edge
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum, invAreaCell
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
real (kind=RKIND), dimension(:), pointer :: &
- 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
+ dvEdge, areaCell, zstarWeight
+ real (kind=RKIND), dimension(:), allocatable:: div_hu, h_tend_col
+ real (kind=RKIND) :: div_hu_btr
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
- boundaryEdge, boundaryCell
+ boundaryEdge, boundaryCell, edgeSignOnCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot, maxLevelVertexTop
- h => s1 % h % array
- h_edge => s1 % h_edge % array
- uTransport => s2 % uTransport % array
- wTop => s2 % wTop % array
+ err = 0
+ nEdgesOnCell => grid % nEdgesOnCell % array
areaCell => grid % areaCell % array
cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
dvEdge => grid % dvEdge % array
@@ -906,67 +925,57 @@
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
- allocate(div_hu(nVertLevels,nCells+1), div_hu_btr(nCells+1), &
- h_tend_col(nVertLevels))
+ if (config_vert_grid_type.eq.'isopycnal') then
+ ! set vertical velocity to zero in isopycnal case
+ wTop=0.0_RKIND
+ return
+ end if
+
+ allocate(div_hu(nVertLevels), h_tend_col(nVertLevels))
+
!
! Compute div(h^{edge} u) for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
!
- div_hu(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(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
- end do
do iCell=1,nCells
- div_hu_btr(iCell) = 0.0
- do k=1,maxLevelCell(iCell)
- div_hu(k,iCell) = div_hu(k,iCell) / areaCell(iCell)
- div_hu_btr(iCell) = div_hu_btr(iCell) + div_hu(k,iCell)
- end do
- end do
+ div_hu(:) = 0.0_RKIND
+ div_hu_btr = 0.0_RKIND
+ hSum = 0.0_RKIND
+ invAreaCell = 1.0_RKIND / areaCell(iCell)
- !
- ! vertical velocity through layer interface
- !
- !dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
- if (config_vert_grid_type.eq.'isopycnal') then
- ! set vertical velocity to zero in isopycnal case
- wTop=0.0
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
- else ! zlevel or zstar type vertical grid
+ do k = 1, maxLevelEdgeBot(iEdge)
+ flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
+ flux = edgeSignOnCell(i, iCell) * flux * invAreaCell
+ div_hu(k) = div_hu(k) - flux
+ div_hu_btr = div_hu_btr - flux
+ end do
+ end do
- do iCell=1,nCells
+ do k = 1, maxLevelCell(iCell)
+ h_tend_col(k) = - zstarWeight(k) * h(k, iCell) * div_hu_btr
+ hSum = hSum + zstarWeight(k) * h(k, iCell)
+ end do
- hSum = 0.0
- do k=1,maxLevelCell(iCell)
- h_tend_col(k) = - zstarWeight(k)*h(k,iCell)*div_hu_btr(iCell)
- hSum = hSum + zstarWeight(k)*h(k,iCell)
- end do
- if(hSum > 0.0) then
- h_tend_col = h_tend_col / hSum
- else
- end if
+ if(hSum > 0.0) then
+ h_tend_col = h_tend_col / hSum
+ end if
- ! 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
+ ! Vertical velocity through layer interface at top and
+ ! bottom is zero.
+ wTop(1,iCell) = 0.0_RKIND
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND
+ do k=maxLevelCell(iCell),2,-1
+ wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k) - h_tend_col(k)
end do
+ end do
- endif
+ deallocate(div_hu, h_tend_col)
- deallocate(div_hu, div_hu_btr, h_tend_col)
-
end subroutine ocn_wtop!}}}
!***********************************************************************
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_test_cases.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_test_cases.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_test_cases.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -342,10 +342,10 @@
do iCell=1,grid % nCells
if (grid % lonCell % array(iCell) < 0.0) grid % lonCell % array(iCell) = grid % lonCell % array(iCell) + 2.0 * pii
r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- grid % h_s % array(iCell) = hs0 * (1.0 - r/rr)
+ grid % bottomDepth % array(iCell) = hs0 * (1.0 - r/rr)
end do
! output about mountain
-print *, 'h_s',minval(grid % h_s % array),sum(grid % h_s % array)/grid % nCells, maxval(grid % h_s % array)
+print *, 'bottomDepth',minval(grid % bottomDepth % array),sum(grid % bottomDepth % array)/grid % nCells, maxval(grid % bottomDepth % array)
!
! Initialize tracer fields
@@ -372,7 +372,7 @@
)**2.0 &
) / &
gravity
- state % h % array(1,iCell) = state % h % array(1,iCell) - grid % h_s % array(iCell)
+ state % h % array(1,iCell) = state % h % array(1,iCell) - grid % bottomDepth % array(iCell)
end do
end subroutine sw_test_case_5
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_thick_hadv.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_thick_hadv.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_thick_hadv.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -101,13 +101,13 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k
+ integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k, i
integer :: iCell, nCells
- integer, dimension(:), pointer :: maxLevelEdgeBot, MaxLevelCell
- integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelEdgeBot, MaxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgeSignOnCell
- real (kind=RKIND) :: flux, invAreaCell1, invAreaCell2
+ real (kind=RKIND) :: flux, invAreaCell, invAreaCell1, invAreaCell2
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
!-----------------------------------------------------------------
@@ -130,20 +130,20 @@
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend(k,cell1) = tend(k,cell1) - flux
- tend(k,cell2) = tend(k,cell2) + flux
- end do
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
+ do iCell = 1, nCells
+ invAreaCell = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelEdgeBot(iEdge)
+ flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
+ tend(k, iCell) = tend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell
+ end do
+ end do
end do
- do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- tend(k,iCell) = tend(k,iCell) / areaCell(iCell)
- end do
- end do
!--------------------------------------------------------------------
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -149,17 +149,19 @@
block => domain % blocklist
do while (associated(block))
- ! mrp 111206 put ocn_wtop call at top for ALE
- call ocn_wtop(block % provis, block % provis, block % mesh)
-
if (.not.config_implicit_vertical_mix) then
call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, err)
end if
- call ocn_tend_h(block % tend, block % provis, block % mesh)
+
+ ! advection of u uses u, while advection of h and tracers use uTransport.
+ call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
+ block % provis % u % array, block % provis % wTop % array, err)
call ocn_tend_u(block % tend, block % provis, block % diagnostics, block % mesh)
- ! mrp 110718 filter btr mode out of u_tend
- ! still got h perturbations with just this alone. Try to set uBtr=0 after full u computation
+ call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
+ block % provis % uTransport % array, block % provis % wTop % array, err)
+ call ocn_tend_h(block % tend, block % provis, block % mesh)
+
if (config_rk_filter_btr_mode) then
call ocn_filter_btr_mode_tend_u(block % tend, block % provis, block % mesh)
endif
@@ -236,9 +238,8 @@
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- ( block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell) )&
- / block % state % time_levs(2) % state % h % array(k, iCell)
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
end do
end do
@@ -256,10 +257,33 @@
! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
!
call mpas_timer_start("RK4-cleaup phase")
+
+ ! Rescale tracers
+ block => domain % blocklist
+ do while(associated(block))
+ do iCell = 1, block % mesh % nCells
+ do k = 1, block % mesh % maxLevelCell % array(iCell)
+ block % state % time_levs(2) % state % tracers % array(:, k, iCell) = block % state % time_levs(2) % state % tracers % array(:, k, iCell) &
+ / block % state % time_levs(2) % state % h % array(k, iCell)
+ end do
+ end do
+ block => block % next
+ end do
+
+
if (config_implicit_vertical_mix) then
call mpas_timer_start("RK4-implicit vert mix")
block => domain % blocklist
do while(associated(block))
+
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
block => block % next
end do
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_time_integration_split.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_time_integration_split.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -86,9 +86,9 @@
type (dm_info) :: dminfo
integer :: iCell, i,k,j, iEdge, cell1, cell2, split_explicit_step, split, &
eoe, oldBtrSubcycleTime, newBtrSubcycleTime, uPerpTime, BtrCorIter, &
- n_bcl_iter(config_n_ts_iter)
+ n_bcl_iter(config_n_ts_iter), stage1_tend_time
type (block_type), pointer :: block
- real (kind=RKIND) :: uhSum, hSum, flux, sshEdge, &
+ real (kind=RKIND) :: uhSum, hSum, flux, sshEdge, hEdge1, &
CoriolisTerm, uCorr, temp, temp_h, coef, FBtr_coeff, sshCell1, sshCell2
integer :: num_tracers, ucorr_coef, err
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -117,6 +117,9 @@
! The baroclinic velocity needs be recomputed at the beginning of a
! timestep because the implicit vertical mixing is conducted on the
! total u. We keep uBtr from the previous timestep.
+ ! Note that uBcl may now include a barotropic component, because the
+ ! weights h have changed. That is OK, because the GBtrForcing variable
+ ! subtracts out the barotropic component from the baroclinic.
block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
= block % state % time_levs(1) % state % u % array(k,iEdge) &
- block % state % time_levs(1) % state % uBtr % array( iEdge)
@@ -181,10 +184,22 @@
block => domain % blocklist
do while (associated(block))
+
+ stage1_tend_time = min(split_explicit_step,2)
+
if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
+ call ocn_vmix_coefs(block % mesh, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, err)
end if
- call ocn_tend_u(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+
+ ! compute wTop. Use u (rather than uTransport) for momentum advection.
+ ! Use the most recent time level available.
+ call ocn_wtop(block % mesh, block % state % time_levs(stage1_tend_time) % state % h % array, &
+ block % state % time_levs(stage1_tend_time) % state % h_edge % array, &
+ block % state % time_levs(stage1_tend_time) % state % u % array, &
+ block % state % time_levs(stage1_tend_time) % state % wTop % array, err)
+
+ call ocn_tend_u(block % tend, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, block % mesh)
+
block => block % next
end do
@@ -246,6 +261,7 @@
= 0.5*( &
block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
+ uTemp(k) - dt * block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
+
enddo
enddo ! iEdge
@@ -408,21 +424,63 @@
! config_btr_gam1_uWt1=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam1_uWt1= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
+
+ do iCell = 1, block % mesh % nCells
+ do i = 1, block % mesh % nEdgesOnCell % array(iCell)
+ iEdge = block % mesh % edgesOnCell % array(i, iCell)
+
+ cell1 = block % mesh % cellsOnEdge % array(1, iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2, iEdge)
+
+ sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
+
+ ! method 0: orig, works only without pbc:
+ !hSum = sshEdge + block % mesh % refBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ ! method 1, matches method 0 without pbcs, works with pbcs.
+ hSum = sshEdge + min(block % mesh % bottomDepth % array(cell1), &
+ block % mesh % bottomDepth % array(cell2))
+
+ ! method 2: may be better than method 1.
+ ! Take average of full thickness at two neighboring cells.
+ !hSum = sshEdge + 0.5 *( block % mesh % bottomDepth % array(cell1) &
+ ! + block % mesh % bottomDepth % array(cell2) )
+
+
+ flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * hSum
+
+ block % tend % ssh % array(iCell) = block % tend % ssh % array(iCell) + block % mesh % edgeSignOncell % array(i, iCell) * flux &
+ * block % mesh % dvEdge % array(iEdge)
+
+ end do
+ end do
+
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
-
+
sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
- hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
-
+
+ ! method 0: orig, works only without pbc:
+ !hSum = sshEdge + block % mesh % refBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ ! method 1, matches method 0 without pbcs, works with pbcs.
+ hSum = sshEdge + min(block % mesh % bottomDepth % array(cell1), &
+ block % mesh % bottomDepth % array(cell2))
+
+ ! method 2: may be better than method 1.
+ ! take average of full thickness at two neighboring cells
+ !hSum = sshEdge + 0.5 *( block % mesh % bottomDepth % array(cell1) &
+ ! + block % mesh % bottomDepth % array(cell2) )
+
flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
* hSum
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
-
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ FBtr_coeff*flux
end do
@@ -452,6 +510,8 @@
block => domain % blocklist
do while (associated(block))
+ allocate(utemp(block % mesh % nEdges+1))
+ uTemp(:) = block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -461,7 +521,8 @@
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
eoe = block % mesh % edgesOnEdge % array(i,iEdge)
CoriolisTerm = CoriolisTerm + block % mesh % weightsOnEdge % array(i,iEdge) &
- * block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ !* block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ * uTemp(eoe) &
* block % mesh % fEdge % array(eoe)
end do
@@ -478,6 +539,7 @@
+ dt/config_n_btr_subcycles *(CoriolisTerm - gravity *(sshCell2 - sshCell1) /block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1,iEdge)
end do
+ deallocate(uTemp)
block => block % next
end do ! block
@@ -502,6 +564,45 @@
! config_btr_gam3_uWt2=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
+
+ do iCell = 1, block % mesh % nCells
+ do i = 1, block % mesh % nEdgesOnCell % array(iCell)
+ iEdge = block % mesh % edgesOnCell % array(i, iCell)
+
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ ! SSH is a linear combination of SSHold and SSHnew.
+ sshCell1 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
+ sshCell2 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+
+ sshEdge = 0.5 * (sshCell1 + sshCell2)
+
+ ! method 0: orig, works only without pbc:
+ !hSum = sshEdge + block % mesh % refBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ ! method 1, matches method 0 without pbcs, works with pbcs.
+ hSum = sshEdge + min(block % mesh % bottomDepth % array(cell1), &
+ block % mesh % bottomDepth % array(cell2))
+
+ ! method 2: may be better than method 1.
+ ! take average of full thickness at two neighboring cells
+ !hSum = sshEdge + 0.5 *( block % mesh % bottomDepth % array(cell1) &
+ ! + block % mesh % bottomDepth % array(cell2) )
+
+
+ flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * hSum
+
+ block % tend % ssh % array(iCell) = block % tend % ssh % array(iCell) + block % mesh % edgeSignOnCell % array(i, iCell) * flux &
+ * block % mesh % dvEdge % array(iEdge)
+
+ end do
+ end do
+
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -511,17 +612,24 @@
+ config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
sshCell2 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+ sshEdge = 0.5 * (sshCell1 + sshCell2)
- sshEdge = 0.5 * (sshCell1 + sshCell2)
- hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+ ! method 0: orig, works only without pbc:
+ !hSum = sshEdge + block % mesh % refBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ ! method 1, matches method 0 without pbcs, works with pbcs.
+ hSum = sshEdge + min(block % mesh % bottomDepth % array(cell1), &
+ block % mesh % bottomDepth % array(cell2))
+
+ ! method 2, better, I think.
+ ! take average of full thickness at two neighboring cells
+ !hSum = sshEdge + 0.5 *( block % mesh % bottomDepth % array(cell1) &
+ ! + block % mesh % bottomDepth % array(cell2) )
flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
* hSum
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
-
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) + flux
end do
@@ -675,8 +783,14 @@
! dwj: 02/22/12 splitting thickness and tracer tendency computations and halo updates to allow monotonic advection.
block => domain % blocklist
do while (associated(block))
- call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
+ ! compute wTop. Use uTransport for advection of h and tracers.
+ ! Use time level 1 values of h and h_edge because h has not yet been computed for time level 2.
+ call ocn_wtop(block % mesh, block % state % time_levs(1) % state % h % array, &
+ block % state % time_levs(1) % state % h_edge % array, &
+ block % state % time_levs(2) % state % uTransport % array, &
+ block % state % time_levs(2) % state % wTop % array, err)
+
call ocn_tend_h(block % tend, block % state % time_levs(2) % state, block % mesh)
block => block % next
end do
@@ -829,7 +943,17 @@
call mpas_timer_start("se implicit vert mix")
block => domain % blocklist
do while(associated(block))
+
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+
block => block % next
end do
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -64,7 +64,7 @@
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
integer, dimension(:), pointer :: nAdvCellsForEdge, maxLevelCell, maxLevelEdgeTop, nEdgesOnCell
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask, edgeSignOnCell, edgesOnCell
real (kind=RKIND) :: flux_upwind, tracer_min_new, tracer_max_new, tracer_upwind_new, scale_factor
real (kind=RKIND) :: flux, tracer_weight, invDvEdge, invAreaCell1, invAreaCell2
@@ -95,6 +95,8 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
highOrderAdvectionMask => grid % highOrderAdvectionMask % array
lowOrderAdvectionMask => grid % lowOrderAdvectionMask % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
nCells = grid % nCells
nCellsSolve = grid % nCellsSolve
@@ -205,8 +207,6 @@
end do ! i loop over nAdvCellsForEdge
end do ! iEdge loop
- call mpas_dmpar_exch_halo_field(high_order_horiz_flux_field)
-
! low order upwind vertical flux (monotonic and diffused)
! Remove low order flux from the high order flux.
! Store left over high order flux in high_order_vert_flux array.
@@ -214,7 +214,7 @@
do iCell = 1, nCells
do k = 2, maxLevelCell(iCell)
! dwj 02/03/12 Ocean and Atmosphere are different in vertical
-! flux_upwind = max(0.,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.,w(k,iCell))*tracer_cur(k,iCell)
+! flux_upwind = max(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.0_RKIND,w(k,iCell))*tracer_cur(k,iCell)
flux_upwind = min(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.0_RKIND,w(k,iCell))*tracer_cur(k,iCell)
upwind_tendency(k-1,iCell) = upwind_tendency(k-1,iCell) + flux_upwind
upwind_tendency(k ,iCell) = upwind_tendency(k ,iCell) - flux_upwind
@@ -227,8 +227,8 @@
! it is negative
do k = 1, maxLevelCell(iCell)
! dwj 02/03/12 Ocean and Atmosphere are different in vertical
-! flux_incoming (k,iCell) = -(min(0.,high_order_vert_flux(k+1,iCell))-max(0.,high_order_vert_flux(k,iCell)))
-! flux_outgoing(k,iCell) = -(max(0.,high_order_vert_flux(k+1,iCell))-min(0.,high_order_vert_flux(k,iCell)))
+! flux_incoming (k,iCell) = -(min(0.0_RKIND,high_order_vert_flux(k+1,iCell))-max(0.0_RKIND,high_order_vert_flux(k,iCell)))
+! flux_outgoing(k,iCell) = -(max(0.0_RKIND,high_order_vert_flux(k+1,iCell))-min(0.0_RKIND,high_order_vert_flux(k,iCell)))
flux_incoming (k, iCell) = max(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - min(0.0_RKIND, high_order_vert_flux(k, iCell))
flux_outgoing(k, iCell) = min(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - max(0.0_RKIND, high_order_vert_flux(k, iCell))
@@ -249,18 +249,27 @@
do k = 1, maxLevelEdgeTop(iEdge)
flux_upwind = dvEdge(iEdge) * (max(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell2))
high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) - flux_upwind
-
- upwind_tendency(k,cell1) = upwind_tendency(k,cell1) - flux_upwind * invAreaCell1
- upwind_tendency(k,cell2) = upwind_tendency(k,cell2) + flux_upwind * invAreaCell2
-
- ! Accumulate remaining high order fluxes
- flux_outgoing(k,cell1) = flux_outgoing(k,cell1) - max(0.0_RKIND,high_order_horiz_flux(k,iEdge)) * invAreaCell1
- flux_incoming (k,cell1) = flux_incoming (k,cell1) - min(0.0_RKIND,high_order_horiz_flux(k,iEdge)) * invAreaCell1
- flux_outgoing(k,cell2) = flux_outgoing(k,cell2) + min(0.0_RKIND,high_order_horiz_flux(k,iEdge)) * invAreaCell2
- flux_incoming (k,cell2) = flux_incoming (k,cell2) + max(0.0_RKIND,high_order_horiz_flux(k,iEdge)) * invAreaCell2
end do ! k loop
end do ! iEdge loop
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ flux_upwind = dvEdge(iEdge) * (max(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell2))
+
+ upwind_tendency(k,iCell) = upwind_tendency(k,iCell) + edgeSignOncell(i, iCell) * flux_upwind * invAreaCell1
+
+ ! Accumulate remaining high order fluxes
+ flux_outgoing(k,iCell) = flux_outgoing(k,iCell) + min(0.0_RKIND, edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge)) * invAreaCell1
+ flux_incoming(k,iCell) = flux_incoming(k,iCell) + max(0.0_RKIND, edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge)) * invAreaCell1
+ end do
+ end do
+ end do
+
! Build the factors for the FCT
! Computed using the bounds that were computed previously, and the bounds on the newly updated value
! Factors are placed in the flux_incoming and flux_outgoing arrays
@@ -295,8 +304,8 @@
do k = 2, maxLevelCell(iCell)
flux = high_order_vert_flux(k,iCell)
! dwj 02/03/12 Ocean and Atmosphere are different in vertical.
-! flux = max(0.,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell)) &
-! + min(0.,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell))
+! flux = max(0.0_RKIND,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell)) &
+! + min(0.0_RKIND,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell))
flux = max(0.0_RKIND,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell)) &
+ min(0.0_RKIND,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell))
high_order_vert_flux(k,iCell) = flux
@@ -304,24 +313,20 @@
end do ! iCell loop
! Accumulate the scaled high order horizontal tendencies
- do iEdge = 1, nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
- do k=1,maxLevelEdgeTop(iEdge)
- tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+ if(config_check_monotonicity) then
+ tracer_new(k, iCell) = tracer_new(k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ end if
+ end do
+ end do
+ end do
- if (config_check_monotonicity) then
- !tracer_new holds a tendency for now.
- tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
- end if
- end do ! k loop
- end do ! iEdge loop
-
! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
do iCell = 1, nCellsSolve
do k = 1,maxLevelCell(iCell)
@@ -388,7 +393,7 @@
if ( config_horiz_tracer_adv_order == 3) then
coef_3rd_order = config_coef_3rd_order
else if(config_horiz_tracer_adv_order == 2 .or. config_horiz_tracer_adv_order == 4) then
- coef_3rd_order = 0.0_RKIND
+ coef_3rd_order = 0.0
end if
end subroutine mpas_ocn_tracer_advection_mono_init!}}}
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -107,13 +107,13 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, nVertLevels, cell1, cell2
- integer :: k, iTracer, num_tracers
+ integer :: iCell, iEdge, nCells, nEdges, nVertLevels, cell1, cell2
+ integer :: i, k, iTracer, num_tracers
integer, dimension(:,:), allocatable :: boundaryMask
- integer, dimension(:), pointer :: maxLevelEdgeTop
- integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask
+ integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask, edgesOnCell, edgeSignOnCell
real (kind=RKIND) :: invAreaCell1, invAreaCell2
real (kind=RKIND) :: tracer_turb_flux, flux, r_tmp
@@ -134,6 +134,7 @@
if (.not.del2On) return
nEdges = grid % nEdges
+ nCells = grid % nCells
nVertLevels = grid % nVertLevels
num_tracers = size(tracers, dim=1)
@@ -145,31 +146,37 @@
dcEdge => grid % dcEdge % array
meshScalingDel2 => grid % meshScalingDel2 % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0/areaCell(cell1)
- invAreaCell2 = 1.0/areaCell(cell2)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOncell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers
+ r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
+
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers
! \kappa_2 </font>
<font color="red">abla \phi on edge
- tracer_turb_flux = tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)
+ tracer_turb_flux = tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1)
! div(h \kappa_2 </font>
<font color="gray">abla \phi) at cell center
- flux = h_edge(k,iEdge) * tracer_turb_flux * edgeMask(k, iEdge) * r_tmp
+ flux = h_edge(k, iEdge) * tracer_turb_flux * edgeMask(k, iEdge) * r_tmp
- tend(iTracer,k,cell1) = tend(iTracer,k,cell1) + flux * invAreaCell1
- tend(iTracer,k,cell2) = tend(iTracer,k,cell2) - flux * invAreaCell2
- end do
- end do
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) - edgeSignOnCell(i, iCell) * flux * invAreaCell1
+ end do
+ end do
+ end do
end do
+
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del2_tend!}}}
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -108,10 +108,10 @@
!-----------------------------------------------------------------
integer :: iEdge, nEdges, num_tracers, nVertLevels, nCells
- integer :: iTracer, k, iCell, cell1, cell2
+ integer :: iTracer, k, iCell, cell1, cell2, i
- integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
- integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge, edgesOnCell, edgeSignOnCell
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux, flux, invdcEdge, r_tmp1, r_tmp2
@@ -148,56 +148,55 @@
edgeMask => grid % edgeMask % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
delsq_tracer(:,:,:) = 0.0
! first del2: div(h </font>
<font color="red">abla \phi) at cell center
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ invdcEdge = dvEdge(iEdge) / dcEdge(iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- invdcEdge = 1.0 / dcEdge(iEdge)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers * edgeMask(k, iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+ r_tmp1 = invdcEdge * h_edge(k, iEdge) * tracers(iTracer, k, cell1)
+ r_tmp2 = invdcEdge * h_edge(k, iEdge) * tracers(iTracer, k, cell2)
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers * edgeMask(k, iEdge)
-
- r_tmp1 = dvEdge(iEdge) * h_edge(k,iEdge) * invdcEdge
-
- r_tmp2 = r_tmp1 * tracers(iTracer,k,cell2)
- r_tmp1 = r_tmp1 * tracers(iTracer,k,cell1)
-
- delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) + (r_tmp2 - r_tmp1) * invAreaCell1
- delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) - (r_tmp2 - r_tmp1) * invAreaCell2
- end do
- end do
+ delsq_tracer(iTracer, k, iCell) = delsq_tracer(iTracer, k, iCell) - edgeSignOnCell(i, iCell) * (r_tmp2 - r_tmp1) * invAreaCell1
+ end do
+ end do
+ end do
end do
! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ cell1 = cellsOnEdge(1, iEdge)
+ cell2 = cellsOnedge(2, iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+ invdcEdge = meshScalingDel4(iEdge) * dvEdge(iEdge) * eddyDiff4 / dcEdge(iEdge)
- invdcEdge = 1.0 / dcEdge(iEdge)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers * edgeMask(k, iEdge)
+ tracer_turb_flux = (delsq_tracer(iTracer, k, cell2) - delsq_tracer(iTracer, k, cell1))
+
+ flux = tracer_turb_flux * invdcEdge
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers * edgeMask(k,iEdge)
- tracer_turb_flux = meshScalingDel4(iEdge) * eddyDiff4 &
- * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) &
- * invdcEdge
-
- flux = dvEdge (iEdge) * tracer_turb_flux
-
- tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux * invAreaCell1
- tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux * invAreaCell2
- enddo
- enddo
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell1
+ end do
+ end do
+ end do
end do
deallocate(delsq_tracer)
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -20,6 +20,7 @@
use mpas_configure
use mpas_timer
use ocn_vel_hmix_del2
+ use ocn_vel_hmix_leith
use ocn_vel_hmix_del4
implicit none
@@ -47,7 +48,7 @@
!
!--------------------------------------------------------------------
- type (timer_node), pointer :: del2Timer, del4Timer
+ type (timer_node), pointer :: del2Timer, leithTimer, del4Timer
!***********************************************************************
@@ -72,7 +73,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vel_hmix_tend(grid, divergence, vorticity, tend, err)!{{{
+ subroutine ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -98,6 +99,9 @@
real (kind=RKIND), dimension(:,:), intent(inout) :: &
tend !< Input/Output: velocity tendency
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ viscosity !< Input: viscosity
+
!-----------------------------------------------------------------
!
! output variables
@@ -112,7 +116,7 @@
!
!-----------------------------------------------------------------
- integer :: err1, err2
+ integer :: err1, err2, err3
!-----------------------------------------------------------------
!
@@ -122,14 +126,21 @@
!
!-----------------------------------------------------------------
+ viscosity = 0.0
+
call mpas_timer_start("del2", .false., del2Timer)
- call ocn_vel_hmix_del2_tend(grid, divergence, vorticity, tend, err1)
+ call ocn_vel_hmix_del2_tend(grid, divergence, vorticity, viscosity, tend, err1)
call mpas_timer_stop("del2", del2Timer)
+
+ call mpas_timer_start("leith", .false., leithTimer)
+ call ocn_vel_hmix_leith_tend(grid, divergence, vorticity, viscosity, tend, err2)
+ call mpas_timer_stop("leith", leithTimer)
+
call mpas_timer_start("del4", .false., del4Timer)
- call ocn_vel_hmix_del4_tend(grid, divergence, vorticity, tend, err2)
+ call ocn_vel_hmix_del4_tend(grid, divergence, vorticity, tend, err3)
call mpas_timer_stop("del4", del4Timer)
- err = ior(err1, err2)
+ err = ior(ior(err1, err2),err3)
!--------------------------------------------------------------------
@@ -163,12 +174,13 @@
integer, intent(out) :: err !< Output: error flag
- integer :: err1, err2
+ integer :: err1, err2, err3
call ocn_vel_hmix_del2_init(err1)
- call ocn_vel_hmix_del4_init(err2)
+ call ocn_vel_hmix_leith_init(err2)
+ call ocn_vel_hmix_del4_init(err3)
- err = ior(err1, err2)
+ err = ior(ior(err1, err2),err3)
!--------------------------------------------------------------------
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_del2.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -70,7 +70,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vel_hmix_del2_tend(grid, divergence, vorticity, tend, err)!{{{
+ subroutine ocn_vel_hmix_del2_tend(grid, divergence, vorticity, viscosity, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -96,6 +96,8 @@
real (kind=RKIND), dimension(:,:), intent(inout) :: &
tend !< Input/Output: velocity tendency
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ viscosity !< Input: viscosity
!-----------------------------------------------------------------
!
@@ -111,12 +113,11 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdgesSolve, cell1, cell2, vertex1, vertex2
- integer :: k
+ integer :: iEdge, nEdgesSolve, cell1, cell2, vertex1, vertex2, k
integer, dimension(:), pointer :: maxLevelEdgeTop
integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
- real (kind=RKIND) :: u_diffusion, invLength1, invLength2
+ real (kind=RKIND) :: u_diffusion, invLength1, invLength2, visc2
real (kind=RKIND), dimension(:), pointer :: meshScalingDel2, &
dcEdge, dvEdge
@@ -158,10 +159,12 @@
-viscVortCoef &
*( vorticity(k,vertex2) - vorticity(k,vertex1) ) * invLength2
- u_diffusion = meshScalingDel2(iEdge) * eddyVisc2 * u_diffusion
+ visc2 = meshScalingDel2(iEdge) * eddyVisc2
- tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * u_diffusion
+ tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * visc2 * u_diffusion
+ viscosity(k,iEdge) = viscosity(k,iEdge) + visc2
+
end do
end do
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -112,22 +112,22 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, cell1, cell2, vertex1, vertex2, k
+ integer :: iEdge, cell1, cell2, vertex1, vertex2, k, i
integer :: iCell, iVertex
- integer :: nVertices, nVertLevels, nCells, nEdges, nEdgesSolve
+ integer :: nVertices, nVertLevels, nCells, nEdges, nEdgesSolve, vertexDegree
- integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelVertexBot, &
- maxLevelCell
- integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
+ integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelVertexTop, &
+ maxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask, edgesOnVertex, edgesOnCell, edgeSignOnVertex, edgeSignOnCell
real (kind=RKIND) :: u_diffusion, invAreaCell1, invAreaCell2, invAreaTri1, &
- invAreaTri2, invDcEdge, invDvEdge, r_tmp, delsq_u
+ invAreaTri2, invDcEdge, invDvEdge, r_tmp
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaTriangle, &
meshScalingDel4, areaCell
real (kind=RKIND), dimension(:,:), allocatable :: delsq_divergence, &
- delsq_circulation, delsq_vorticity
+ delsq_circulation, delsq_vorticity, delsq_u
err = 0
@@ -138,8 +138,10 @@
nEdgesSolve = grid % nEdgessolve
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
+
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
+ maxLevelVertexTop => grid % maxLevelVertexTop % array
maxLevelCell => grid % maxLevelCell % array
cellsOnEdge => grid % cellsOnEdge % array
verticesOnEdge => grid % verticesOnEdge % array
@@ -149,43 +151,57 @@
areaCell => grid % areaCell % array
meshScalingDel4 => grid % meshScalingDel4 % array
edgeMask => grid % edgeMask % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+ allocate(delsq_u(nVertLEvels, nEdges+1))
allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_vorticity(nVertLevels, nVertices+1))
+ delsq_u(:,:) = 0.0
delsq_vorticity(:,:) = 0.0
delsq_divergence(:,:) = 0.0
- do iEdge=1,nEdges
+ !Compute delsq_u
+ do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- invAreaTri1 = 1.0 / areaTriangle(vertex1)
- invAreaTri2 = 1.0 / areaTriangle(vertex2)
-
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
-
invDcEdge = 1.0 / dcEdge(iEdge)
invDvEdge = 1.0 / dvEdge(iEdge)
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 &
+ delsq_u(k, iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
-viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDcEdge * sqrt(3.0) ! TDR
+ end do
+ end do
- ! vorticity using </font>
<font color="red">abla^2 u
- r_tmp = dcEdge(iEdge) * delsq_u
- delsq_vorticity(k,vertex1) = delsq_vorticity(k,vertex1) - r_tmp * invAreaTri1
- delsq_vorticity(k,vertex2) = delsq_vorticity(k,vertex2) + r_tmp * invAreaTri2
+ ! Compute delsq_vorticity
+ do iVertex = 1, nVertices
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ do k = 1, maxLevelVertexTop(iVertex)
+ delsq_vorticity(k, iVertex) = delsq_vorticity(k, iVertex) + edgeSignOnVertex(i, iVertex) * dcEdge(iEdge) * delsq_u(k, iEdge) * invAreaTri1
+ end do
+ end do
+ end do
- ! Divergence using </font>
<font color="gray">abla^2 u
- r_tmp = dvEdge(iEdge) * delsq_u
- delsq_divergence(k, cell1) = delsq_divergence(k,cell1) + r_tmp * invAreaCell1
- delsq_divergence(k, cell2) = delsq_divergence(k,cell2) - r_tmp * invAreaCell2
+ ! Compute delsq_divergence
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ delsq_divergence(k, iCell) = delsq_divergence(k, iCell) - edgeSignOnCell(i, iCell) * dvEdge(iEdge) * delsq_u(k, iEdge) * invAreaCell1
+ end do
end do
end do
@@ -209,6 +225,7 @@
end do
end do
+ deallocate(delsq_u)
deallocate(delsq_divergence)
deallocate(delsq_vorticity)
Copied: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_leith.F (from rev 2314, trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_leith.F)
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_leith.F (rev 0)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -0,0 +1,235 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! ocn_vel_hmix_leith
+!
+!> \brief Ocean horizontal mixing - Leith parameterization
+!> \author Mark Petersen
+!> \date 22 October 2012
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal mixing
+!> tendencies using the Leith parameterization.
+!
+!-----------------------------------------------------------------------
+
+module ocn_vel_hmix_leith
+
+ use mpas_grid_types
+ use mpas_configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: ocn_vel_hmix_leith_tend, &
+ ocn_vel_hmix_leith_init
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: hmixLeithOn !< integer flag to determine whether leith chosen
+
+ real (kind=RKIND) :: &
+ viscVortCoef
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine ocn_vel_hmix_leith_tend
+!
+!> \brief Computes tendency term for horizontal momentum mixing with Leith parameterization
+!> \author Mark Petersen, Todd Ringler
+!> \date 22 October 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for momentum
+!> based on the Leith closure. The Leith closure is the
+!> enstrophy-cascade analogy to the Smagorinsky (1963) energy-cascade
+!> closure, i.e. Leith (1996) assumes an inertial range of enstrophy flux
+!> moving toward the grid scale. The assumption of an enstrophy cascade
+!> and dimensional analysis produces right-hand-side dissipation,
+!> $\bf{D}$, of velocity of the form
+!> $ {\bf D} = </font>
<font color="black">abla \cdot \left( </font>
<font color="black">u_\ast </font>
<font color="blue">abla {\bf u} \right)
+!> = </font>
<font color="black">abla \cdot \left( \gamma \left| </font>
<font color="blue">abla \omega \right|
+!> \left( \Delta x \right)^3 </font>
<font color="blue">abla \bf{u} \right)
+!> where $\omega$ is the relative vorticity and $\gamma$ is a non-dimensional,
+!> $O(1)$ parameter. We set $\gamma=1$.
+
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_vel_hmix_leith_tend(grid, divergence, vorticity, viscosity, tend, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ divergence !< Input: velocity divergence
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ vorticity !< Input: vorticity
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ viscosity !< Input: viscosity
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdgesSolve, cell1, cell2, vertex1, vertex2, k
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
+
+ real (kind=RKIND) :: u_diffusion, invLength1, invLength2, visc2
+ real (kind=RKIND), dimension(:), pointer :: meshScaling, &
+ dcEdge, dvEdge
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this mixing is not selected
+ !
+ !-----------------------------------------------------------------
+
+ err = 0
+
+ if(.not.hmixLeithOn) return
+
+ nEdgesSolve = grid % nEdgesSolve
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ meshScaling => grid % meshScaling % array
+ edgeMask => grid % edgeMask % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+
+ do iEdge=1,nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ invLength1 = 1.0 / dcEdge(iEdge)
+ invLength2 = 1.0 / dvEdge(iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ ! Here -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ ! is - </font>
<font color="blue">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
+ ! + k \times </font>
<font color="blue">abla vorticity pointing from cell1 to cell2.
+
+ u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) * invLength1 &
+ -viscVortCoef &
+ *( vorticity(k,vertex2) - vorticity(k,vertex1) ) * invLength2
+
+ ! Here the first line is (\delta x)^3
+ ! the second line is |</font>
<font color="blue">abla \omega|
+ ! and u_diffusion is </font>
<font color="gray">abla^2 u (see formula for $\bf{D}$ above).
+ visc2 = ( config_leith_parameter * config_leith_dx * meshScaling(iEdge) / 3.14)**3 &
+ * abs( vorticity(k,vertex2) - vorticity(k,vertex1) ) * invLength1 * sqrt(3.0)
+ visc2 = min(visc2, config_leith_visc2_max)
+
+ tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * visc2 * u_diffusion
+
+ viscosity(k,iEdge) = viscosity(k,iEdge) + visc2
+
+ end do
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_vel_hmix_leith_tend!}}}
+
+!***********************************************************************
+!
+! routine ocn_vel_hmix_leith_init
+!
+!> \brief Initializes ocean momentum horizontal mixing with Leith parameterization
+!> \author Mark Petersen
+!> \date 22 October 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> Leith parameterization for horizontal momentum mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_vel_hmix_leith_init(err)!{{{
+
+
+ integer, intent(out) :: err !< Output: error flag
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ err = 0
+
+ hmixLeithOn = .false.
+
+ if (config_use_leith_del2) then
+ hmixLeithOn = .true.
+
+ if (config_visc_vorticity_term) then
+ viscVortCoef = config_visc_vorticity_visc2_scale
+ else
+ viscVortCoef = 0.0
+ endif
+
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_vel_hmix_leith_init!}}}
+
+!***********************************************************************
+
+end module ocn_vel_hmix_leith
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -427,13 +427,13 @@
!
!-----------------------------------------------------------------
- integer :: nVertLevels, nCells, nEdges, iCell, iEdge, k
+ integer :: nVertLevels, nCells, nEdges, iCell, iEdge, k, i
integer :: cell1, cell2
- integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot
- integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOncell, edgeSignOnCell
- real (kind=RKIND) :: coef
+ real (kind=RKIND) :: coef, invAreaCell
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaCell
real (kind=RKIND), dimension(:,:), allocatable :: drhoTopOfCell, du2TopOfCell, &
drhoTopOfEdge, du2TopOfEdge
@@ -453,6 +453,9 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
areaCell => grid % areaCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
allocate( &
drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges), &
@@ -498,21 +501,16 @@
! interpolate du2TopOfEdge to du2TopOfCell
du2TopOfCell = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeBot(iEdge)
- du2TopOfCell(k,cell1) = du2TopOfCell(k,cell1) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
- du2TopOfCell(k,cell2) = du2TopOfCell(k,cell2) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
- end do
+ do iCell = 1, nCells
+ invAreaCell = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+
+ do k = 2, maxLevelEdgeBot(iEdge)
+ du2TopOfCell(k, iCell) = du2TopOfCell(k, iCell) + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k, iEdge) * invAreaCell
+ end do
+ end do
end do
- do iCell = 1,nCells
- do k = 2,maxLevelCell(iCell)
- du2TopOfCell(k,iCell) = du2TopOfCell(k,iCell) / areaCell(iCell)
- end do
- end do
! compute RiTopOfEdge using drhoTopOfEdge and du2TopOfEdge
! coef = -g/rho_0/2
Modified: branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -177,22 +177,22 @@
integer :: k, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:), pointer :: refBottomDepth
err = 0
if(.not.tanhViscOn) return
nVertLevels = grid % nVertLevels
- referenceBottomDepth => grid % referenceBottomDepth % array
+ refBottomDepth => grid % refBottomDepth % array
- ! referenceBottomDepth is used here for simplicity. Using zMid and h, which
+ ! refBottomDepth is used here for simplicity. Using zMid and h, which
! vary in time, would give the exact location of the top, but it
! would only change the diffusion value very slightly.
vertViscTopOfEdge = 0.0
do k=2,nVertLevels
vertViscTopOfEdge(k,:) = -(config_max_visc_tanh-config_min_visc_tanh)/2.0 &
- *tanh((referenceBottomDepth(k-1)+config_ZMid_tanh) &
+ *tanh((refBottomDepth(k-1)+config_ZMid_tanh) &
/config_zWidth_tanh) &
+ (config_max_visc_tanh+config_min_visc_tanh)/2
end do
@@ -250,22 +250,22 @@
integer :: k, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:), pointer :: refBottomDepth
err = 0
if(.not.tanhDiffOn) return
nVertLevels = grid % nVertLevels
- referenceBottomDepth => grid % referenceBottomDepth % array
+ refBottomDepth => grid % refBottomDepth % array
- ! referenceBottomDepth is used here for simplicity. Using zMid and h, which
+ ! refBottomDepth is used here for simplicity. Using zMid and h, which
! vary in time, would give the exact location of the top, but it
! would only change the diffusion value very slightly.
vertDiffTopOfCell = 0.0
do k=2,nVertLevels
vertDiffTopOfCell(k,:) = -(config_max_diff_tanh-config_min_diff_tanh)/2.0 &
- *tanh((referenceBottomDepth(k-1)+config_ZMid_tanh) &
+ *tanh((refBottomDepth(k-1)+config_ZMid_tanh) &
/config_zWidth_tanh) &
+ (config_max_diff_tanh+config_min_diff_tanh)/2
end do
Modified: branches/ocean_projects/cesm_coupling/src/framework/mpas_dmpar.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/framework/mpas_dmpar.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/framework/mpas_dmpar.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -2831,6 +2831,12 @@
logical :: comm_list_found
+ do i = 1, 1
+ if(field % dimSizes(i) <= 0) then
+ return
+ end if
+ end do
+
dminfo => field % block % domain % dminfo
if(present(haloLayersIn)) then
@@ -3105,6 +3111,12 @@
logical :: comm_list_found
+ do i = 1, 2
+ if(field % dimSizes(i) <= 0) then
+ return
+ end if
+ end do
+
dminfo => field % block % domain % dminfo
if(present(haloLayersIn)) then
@@ -3380,6 +3392,12 @@
logical :: comm_list_found
+ do i = 1, 3
+ if(field % dimSizes(i) <= 0) then
+ return
+ end if
+ end do
+
dminfo => field % block % domain % dminfo
if(present(haloLayersIn)) then
@@ -3661,6 +3679,12 @@
logical :: comm_list_found
+ do i = 1, 1
+ if(field % dimSizes(i) <= 0) then
+ return
+ end if
+ end do
+
dminfo => field % block % domain % dminfo
if(present(haloLayersIn)) then
@@ -3932,7 +3956,13 @@
integer, dimension(:), pointer :: haloLayers
logical :: comm_list_found
-
+
+ do i = 1, 2
+ if(field % dimSizes(i) <= 0) then
+ return
+ end if
+ end do
+
dminfo => field % block % domain % dminfo
if(present(haloLayersIn)) then
@@ -4210,6 +4240,12 @@
logical :: comm_list_found
+ do i = 1, 3
+ if(field % dimSizes(i) <= 0) then
+ return
+ end if
+ end do
+
dminfo => field % block % domain % dminfo
if(present(haloLayersIn)) then
Modified: branches/ocean_projects/cesm_coupling/src/framework/mpas_grid_types.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/framework/mpas_grid_types.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/framework/mpas_grid_types.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -363,6 +363,26 @@
type (dm_info), pointer :: dminfo
end type domain_type
+ interface mpas_allocate_scratch_field
+ module procedure mpas_allocate_scratch_field1d_integer
+ module procedure mpas_allocate_scratch_field2d_integer
+ module procedure mpas_allocate_scratch_field3d_integer
+ module procedure mpas_allocate_scratch_field1d_real
+ module procedure mpas_allocate_scratch_field2d_real
+ module procedure mpas_allocate_scratch_field3d_real
+ module procedure mpas_allocate_scratch_field1d_char
+ end interface
+
+ interface mpas_deallocate_scratch_field
+ module procedure mpas_deallocate_scratch_field1d_integer
+ module procedure mpas_deallocate_scratch_field2d_integer
+ module procedure mpas_deallocate_scratch_field3d_integer
+ module procedure mpas_deallocate_scratch_field1d_real
+ module procedure mpas_deallocate_scratch_field2d_real
+ module procedure mpas_deallocate_scratch_field3d_real
+ module procedure mpas_deallocate_scratch_field1d_char
+ end interface
+
interface mpas_deallocate_field
module procedure mpas_deallocate_field0d_integer
module procedure mpas_deallocate_field1d_integer
@@ -444,6 +464,406 @@
end subroutine mpas_deallocate_domain!}}}
+ subroutine mpas_allocate_scratch_field1d_integer(f, single_block_in)!{{{
+ type (field1dInteger), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field1dInteger), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field1d_integer!}}}
+
+ subroutine mpas_allocate_scratch_field2d_integer(f, single_block_in)!{{{
+ type (field2dInteger), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field2dInteger), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1), f_cursor % dimSizes(2)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1), f % dimSizes(2)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field2d_integer!}}}
+
+ subroutine mpas_allocate_scratch_field3d_integer(f, single_block_in)!{{{
+ type (field3dInteger), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field3dInteger), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1), f_cursor % dimSizes(2), f_cursor % dimSizes(3)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1), f % dimSizes(2), f % dimSizes(3)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field3d_integer!}}}
+
+ subroutine mpas_allocate_scratch_field1d_real(f, single_block_in)!{{{
+ type (field1dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field1dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field1d_real!}}}
+
+ subroutine mpas_allocate_scratch_field2d_real(f, single_block_in)!{{{
+ type (field2dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field2dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1), f_cursor % dimSizes(2)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1), f % dimSizes(2)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field2d_real!}}}
+
+ subroutine mpas_allocate_scratch_field3d_real(f, single_block_in)!{{{
+ type (field3dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field3dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1), f_cursor % dimSizes(2), f_cursor % dimSizes(3)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1), f % dimSizes(2), f % dimSizes(3)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field3d_real!}}}
+
+ subroutine mpas_allocate_scratch_field1d_char(f, single_block_in)!{{{
+ type (field1dChar), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field1dChar), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field1d_char!}}}
+
+ subroutine mpas_deallocate_scratch_field1d_integer(f, single_block_in)!{{{
+ type (field1dInteger), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field1dInteger), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field1d_integer!}}}
+
+ subroutine mpas_deallocate_scratch_field2d_integer(f, single_block_in)!{{{
+ type (field2dInteger), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field2dInteger), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field2d_integer!}}}
+
+ subroutine mpas_deallocate_scratch_field3d_integer(f, single_block_in)!{{{
+ type (field3dInteger), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field3dInteger), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field3d_integer!}}}
+
+ subroutine mpas_deallocate_scratch_field1d_real(f, single_block_in)!{{{
+ type (field1dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field1dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field1d_real!}}}
+
+ subroutine mpas_deallocate_scratch_field2d_real(f, single_block_in)!{{{
+ type (field2dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field2dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field2d_real!}}}
+
+ subroutine mpas_deallocate_scratch_field3d_real(f, single_block_in)!{{{
+ type (field3dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field3dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field3d_real!}}}
+
+ subroutine mpas_deallocate_scratch_field1d_char(f, single_block_in)!{{{
+ type (field1dChar), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field1dChar), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field1d_char!}}}
+
+
subroutine mpas_deallocate_field0d_integer(f)!{{{
type (field0dInteger), pointer :: f
type (field0dInteger), pointer :: f_cursor
Modified: branches/ocean_projects/cesm_coupling/src/framework/mpas_io.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/framework/mpas_io.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/framework/mpas_io.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -375,6 +375,7 @@
if (present(ierr)) ierr = MPAS_IO_ERR_PIO
deallocate(new_dimlist_node % dimhandle)
deallocate(new_dimlist_node)
+ write(0,*) 'WARNING: Dimension ', trim(dimname), ' not in input file.'
dimsize = -1
return
end if
@@ -559,6 +560,7 @@
if (present(ierr)) ierr = MPAS_IO_ERR_PIO
deallocate(new_fieldlist_node % fieldhandle)
deallocate(new_fieldlist_node)
+ write(0,*) 'WARNING: Variable ', trim(fieldname), ' not in input file.'
return
end if
!write(0,*) 'Inquired about variable ID', new_fieldlist_node % fieldhandle % fieldid
Modified: branches/ocean_projects/cesm_coupling/src/framework/mpas_io_input.F
===================================================================
--- branches/ocean_projects/cesm_coupling/src/framework/mpas_io_input.F 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/framework/mpas_io_input.F 2012-11-19 20:17:58 UTC (rev 2315)
@@ -262,7 +262,7 @@
call mpas_dmpar_abort(domain % dminfo)
end if
!write(0,*) 'MGD DEBUGGING time = ', input_obj % time
- write(0,*) 'Restarting model from time ', timeStamp
+ write(0,*) 'Restarting model from time ', trim(timeStamp)
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -516,7 +516,7 @@
integer :: ierr
- call MPAS_readStream(input_obj % io_stream, 1, ierr)
+ call MPAS_readStream(input_obj % io_stream, input_obj % time, ierr)
end subroutine mpas_read_and_distribute_fields!}}}
Modified: branches/ocean_projects/cesm_coupling/src/registry/gen_inc.c
===================================================================
--- branches/ocean_projects/cesm_coupling/src/registry/gen_inc.c 2012-11-19 17:38:23 UTC (rev 2314)
+++ branches/ocean_projects/cesm_coupling/src/registry/gen_inc.c 2012-11-19 20:17:58 UTC (rev 2315)
@@ -143,8 +143,8 @@
fortprintf(fd, " call mpas_dmpar_abort(dminfo)</font>
<font color="black">");
fortprintf(fd, " else if (ierr < 0) then</font>
<font color="black">");
fortprintf(fd, " write(0,*) \'Namelist record &%s not found; using default values for this namelist\'\'s variables\'</font>
<font color="red">",nls_ptr->record);
- fortprintf(fd, " rewind(funit)</font>
<font color="black">");
fortprintf(fd, " end if</font>
<font color="blue">");
+ fortprintf(fd, " rewind(funit)</font>
<font color="gray">");
dict_insert(dictionary, nls_ptr->record);
}
@@ -813,6 +813,10 @@
fortprintf(fd, " %s %% %s %% fieldName = \'%s\'</font>
<font color="black">", group_ptr->name, var_ptr->name_in_code, var_ptr->name_in_file);
fortprintf(fd, " %s %% %s %% isSuperArray = .false.</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code);
if (var_ptr->ndims > 0) {
+ if(var_ptr->persistence == SCRATCH){
+ fortprintf(fd, " ! SCRATCH VARIABLE</font>
<font color="blue">");
+ fortprintf(fd, " nullify(%s %% %s %% array)</font>
<font color="gray">", group_ptr->name, var_ptr->name_in_code);
+ } else if(var_ptr->persistence == PERSISTENT){
fortprintf(fd, " allocate(%s %% %s %% array(", group_ptr->name, var_ptr->name_in_code);
dimlist_ptr = var_ptr->dimlist;
if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
@@ -843,6 +847,7 @@
else if (var_ptr->vtype == CHARACTER)
fortprintf(fd, " %s %% %s %% array = \'\'</font>
<font color="gray">", group_ptr->name, var_ptr->name_in_code ); /* initialize field to zero */
+ }
dimlist_ptr = var_ptr->dimlist;
i = 1;
while (dimlist_ptr) {
@@ -869,7 +874,7 @@
i++;
dimlist_ptr = dimlist_ptr->next;
}
- }
+ }
if (var_ptr->timedim) fortprintf(fd, " %s %% %s %% hasTimeDimension = .true.</font>
<font color="black">", group_ptr->name, var_ptr->name_in_code);
else fortprintf(fd, " %s %% %s %% hasTimeDimension = .false.</font>
<font color="gray">", group_ptr->name, var_ptr->name_in_code);
@@ -934,14 +939,18 @@
var_list_ptr2 = var_list_ptr;
var_list_ptr = var_list_ptr->next;
}
- fortprintf(fd, " deallocate(%s %% %s %% array)</font>
<font color="blue">", group_ptr->name, var_list_ptr2->var->super_array);
+ fortprintf(fd, " if(associated(%s %% %s %% array)) then</font>
<font color="blue">", group_ptr->name, var_list_ptr2->var->super_array);
+ fortprintf(fd, " deallocate(%s %% %s %% array)</font>
<font color="blue">", group_ptr->name, var_list_ptr2->var->super_array);
+ fortprintf(fd, " end if</font>
<font color="black">");
fortprintf(fd, " deallocate(%s %% %s %% ioinfo)</font>
<font color="black">", group_ptr->name, var_list_ptr2->var->super_array);
fortprintf(fd, " call mpas_deallocate_attlist(%s %% %s %% attList)</font>
<font color="black">", group_ptr->name, var_list_ptr2->var->super_array);
fortprintf(fd, " deallocate(%s %% %s)</font>
<font color="black"></font>
<font color="red">", group_ptr->name, var_list_ptr2->var->super_array);
}
else {
if (var_ptr->ndims > 0) {
- fortprintf(fd, " deallocate(%s %% %s %% array)</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code);
+ fortprintf(fd, " if(associated(%s %% %s %% array)) then</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code);
+ fortprintf(fd, " deallocate(%s %% %s %% array)</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code);
+ fortprintf(fd, " end if</font>
<font color="black">");
fortprintf(fd, " deallocate(%s %% %s %% ioinfo)</font>
<font color="black">", group_ptr->name, var_ptr->name_in_code);
fortprintf(fd, " call mpas_deallocate_attlist(%s %% %s %% attList)</font>
<font color="black">", group_ptr->name, var_ptr->name_in_code);
fortprintf(fd, " deallocate(%s %% %s)</font>
<font color="black"></font>
<font color="gray">", group_ptr->name, var_ptr->name_in_code);
@@ -1110,8 +1119,10 @@
{
var_list_ptr = group_ptr->vlist;
var_list_ptr = var_list_ptr->next;
+
+ if (!var_list_ptr) break;
+
var_ptr = var_list_ptr->var;
-
int ntime_levs = 1;
@@ -2126,6 +2137,7 @@
dimlist_ptr = var_ptr->dimlist;
i = 1;
+ if(var_ptr->persistence == PERSISTENT){
while (dimlist_ptr) {
if (i == var_ptr->ndims) {
@@ -2172,6 +2184,7 @@
i++;
dimlist_ptr = dimlist_ptr -> next;
}
+ }
if (var_list_ptr) var_list_ptr = var_list_ptr->next;
}
</font>
</pre>