<p><b>lowrie@lanl.gov</b> 2013-02-13 14:29:30 -0700 (Wed, 13 Feb 2013)</p><p>Merged lastest trunk changes to mpas_cdg_advection branch.<br>
</p><hr noshade><pre><font color="gray">Index: branches/mpas_cdg_advection
===================================================================
--- branches/mpas_cdg_advection 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection 2013-02-13 21:29:30 UTC (rev 2473)
Property changes on: branches/mpas_cdg_advection
___________________________________________________________________
Modified: svn:mergeinfo
## -4,12 +4,15 ##
/branches/ocean_projects/ale_vert_coord:1225-1383
/branches/ocean_projects/ale_vert_coord_new:1387-1428
/branches/ocean_projects/cesm_coupling:2147-2344
+/branches/ocean_projects/diagnostics_revision:2439-2462
/branches/ocean_projects/gmvar:1214-1514,1517-1738
/branches/ocean_projects/imp_vert_mix_error:1847-1887
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/leith_mrp:2182-2241
+/branches/ocean_projects/linear_eos:2435-2437
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/monthly_forcing:1810-1867
+/branches/ocean_projects/namelist_cleanup:2319-2414
/branches/ocean_projects/option3_b4b_test:2201-2231
/branches/ocean_projects/partial_bottom_cells:2172-2226
/branches/ocean_projects/restart_reproducibility:2239-2272
## -25,3 +28,4 ##
/branches/omp_blocks/multiple_blocks:1803-2084
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+/trunk/mpas:2390-2472
\ No newline at end of property
Modified: branches/mpas_cdg_advection/namelist.input.ocean
===================================================================
--- branches/mpas_cdg_advection/namelist.input.ocean 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/namelist.input.ocean 2013-02-13 21:29:30 UTC (rev 2473)
@@ -1,110 +1,167 @@
-&sw_model
- config_test_case = 0
- config_time_integration = 'split_explicit'
- config_rk_filter_btr_mode = .false.
- config_dt = 180.0
- config_start_time = '0000-01-01_00:00:00'
- config_run_duration = '1_00:00:00'
- config_stats_interval = 480
+&time_management
+ config_do_restart = .false.
+ config_start_time = '0000-01-01_00:00:00'
+ config_stop_time = 'none'
+ config_run_duration = '0_06:00:00'
+ config_calendar_type = '360day'
/
&io
- config_input_name = 'grid.nc'
- config_output_name = 'output.nc'
- config_restart_name = 'restart.nc'
- config_output_interval = '1_00:00:00'
- config_frames_per_outfile = 1000000
- config_pio_num_iotasks = 0
- config_pio_stride = 1
- config_write_output_on_startup = .true.
+ config_input_name = 'grid.nc'
+ config_output_name = 'output.nc'
+ config_restart_name = 'restart.nc'
+ config_restart_interval = '0_06:00:00'
+ config_output_interval = '0_06:00:00'
+ config_stats_interval = '0_01:00:00'
+ config_write_stats_on_startup = .true.
+ config_write_output_on_startup = .true.
+ config_frames_per_outfile = 1000
+ config_pio_num_iotasks = 0
+ config_pio_stride = 1
/
+&time_integration
+ config_dt = 3000.0
+ config_time_integrator = 'split_explicit'
+/
+&grid
+ config_num_halos = 3
+ config_vert_coord_movement = 'uniform_stretching'
+ config_alter_ICs_for_pbcs = 'zlevel_pbcs_off'
+ config_min_pbc_fraction = 0.10
+ config_check_ssh_consistency = .true.
+/
&decomposition
- config_number_of_blocks = 0
- config_block_decomp_file_prefix = 'graph.info.part.'
- config_explicit_proc_decomp = .false.
- config_proc_decomp_file_prefix = 'graph.info.part.'
+ config_block_decomp_file_prefix = 'graph.info.part.'
+ config_number_of_blocks = 0
+ config_explicit_proc_decomp = .false.
+ config_proc_decomp_file_prefix = 'graph.info.part.'
/
-&restart
- config_do_restart = .false.
- config_restart_interval = '120_00:00:00'
+&hmix
+ config_hmix_ScaleWithMesh = .false.
+ config_visc_vorticity_term = .true.
+ config_apvm_scale_factor = 0.0
/
-&grid
- config_vert_grid_type = 'isopycnal'
- config_pressure_type = 'pressure'
- config_rho0 = 1014.65
+&hmix_del2
+ config_use_mom_del2 = .false.
+ config_use_tracer_del2 = .false.
+ config_mom_del2 = 0.0
+ config_tracer_del2 = 0.0
+ config_vorticity_del2_scale = 1.0
/
-&partial_bottom_cells
- config_alter_ICs_for_pbcs = 'off'
- config_min_pbc_fraction = 0.10
- config_check_ssh_consistency = .true.
- config_check_zlevel_consistency = .false.
+&hmix_del4
+ config_use_mom_del4 = .true.
+ config_use_tracer_del4 = .false.
+ config_mom_del4 = 5.0e13
+ config_tracer_del4 = 0.0
+ config_vorticity_del4_scale = 1.0
/
-&split_explicit_ts
- config_n_ts_iter = 2
- config_n_bcl_iter_beg = 1
- config_n_bcl_iter_mid = 2
- config_n_bcl_iter_end = 2
- config_n_btr_subcycles = 10
- config_n_btr_cor_iter = 2
- config_u_correction = .true.
- config_filter_btr_mode = .false.
- config_btr_subcycle_loop_factor = 2
- config_btr_gam1_uWt1 = 0.5
- config_btr_gam2_SSHWt1 = 1.0
- config_btr_gam3_uWt2 = 1.0
- config_btr_solve_SSH2 = .false.
+&hmix_Leith
+ config_use_Leith_del2 = .false.
+ config_Leith_parameter = 1.0
+ config_Leith_dx = 15000.0
+ config_Leith_visc2_max = 2.5e3
/
-&hmix
- config_h_mom_eddy_visc2 = 100.0
- config_h_mom_eddy_visc4 = 0.0
- config_h_kappa = 0.0
- config_h_kappa_q = 0.0
- config_visc_vorticity_term = .true.
- config_h_tracer_eddy_diff2 = 1.0e5
- config_h_tracer_eddy_diff4 = 0.0
+&standard_GM
+ config_h_kappa = 0.0
+ config_h_kappa_q = 0.0
/
-&hmix_leith
- config_use_leith_del2 = .false.
- config_leith_parameter = 1.0
- config_leith_dx = 15000.0
- config_leith_visc2_max = 2.5e3
+&Rayleigh_damping
+ config_Rayleigh_friction = .false.
+ config_Rayleigh_damping_coeff = 0.0
/
&vmix
- config_vert_visc_type = 'const'
- config_vert_diff_type = 'const'
- config_implicit_vertical_mix = .true.
- config_convective_visc = 1.0
- config_convective_diff = 1.0
- config_bottom_drag_coeff = 1.0e-3
+ config_implicit_vertical_mix = .true.
+ config_convective_visc = 1.0
+ config_convective_diff = 1.0
/
&vmix_const
- config_vert_visc = 1.0e-5
- config_vert_diff = 1.0e-5
+ config_use_const_visc = .false.
+ config_use_const_diff = .false.
+ config_vert_visc = 2.5e-4
+ config_vert_diff = 2.5e-5
/
&vmix_rich
- config_bkrd_vert_visc = 1.0e-5
- config_bkrd_vert_diff = 1.0e-5
- config_rich_mix = 0.005
+ config_use_rich_visc = .true.
+ config_use_rich_diff = .true.
+ config_bkrd_vert_visc = 1.0e-4
+ config_bkrd_vert_diff = 1.0e-5
+ config_rich_mix = 0.005
/
&vmix_tanh
- config_max_visc_tanh = 2.5e-1
- config_min_visc_tanh = 1.0e-4
- config_max_diff_tanh = 2.5e-2
- config_min_diff_tanh = 1.0e-5
- config_zMid_tanh = -100
- config_zWidth_tanh = 100
+ config_use_tanh_visc = .false.
+ config_use_tanh_diff = .false.
+ config_max_visc_tanh = 2.5e-1
+ config_min_visc_tanh = 1.0e-4
+ config_max_diff_tanh = 2.5e-2
+ config_min_diff_tanh = 1.0e-5
+ config_zMid_tanh = -100
+ config_zWidth_tanh = 100
/
-&eos
- config_eos_type = 'linear'
+&forcing
+ config_use_monthly_forcing = .false.
+ config_restoreTS = .false.
+ config_restoreT_timescale = 90.0
+ config_restoreS_timescale = 90.0
/
&advection
- config_vert_tracer_adv = 'stencil'
- config_vert_tracer_adv_order = 2
- config_horiz_tracer_adv_order = 2
- config_thickness_adv_order = 2
- config_monotonic = .false.
+ config_vert_tracer_adv = 'stencil'
+ config_vert_tracer_adv_order = 3
+ config_horiz_tracer_adv_order = 3
+ config_coef_3rd_order = 0.25
+ config_monotonic = .true.
/
-&restore
- config_restoreTS = .false.
- config_restoreT_timescale = 90.0
- config_restoreS_timescale = 90.0
+&bottom_drag
+ config_bottom_drag_coeff = 1.0e-3
/
+&pressure_gradient
+ config_pressure_gradient_type = 'pressure_and_zmid'
+ config_rho0 = 1014.65
+/
+&eos
+ config_eos_type = 'jm'
+/
+&eos_linear
+ config_eos_linear_alpha = 2.55e-1
+ config_eos_linear_beta = 7.64e-1
+ config_eos_linear_Tref = 19.0
+ config_eos_linear_Sref = 35.0
+ config_eos_linear_rhoref = 1025.022
+/
+&split_explicit_ts
+ config_n_ts_iter = 2
+ config_n_bcl_iter_beg = 1
+ config_n_bcl_iter_mid = 2
+ config_n_bcl_iter_end = 2
+ config_n_btr_subcycles = 20
+ config_n_btr_cor_iter = 2
+ config_u_correction = .true.
+ config_btr_subcycle_loop_factor = 2
+ config_btr_gam1_uWt1 = 0.5
+ config_btr_gam2_SSHWt1 = 1.0
+ config_btr_gam3_uWt2 = 1.0
+ config_btr_solve_SSH2 = .false.
+/
+&sw_model
+ config_test_case = 0
+/
+&debug
+ config_check_zlevel_consistency = .false.
+ config_filter_btr_mode = .false.
+ config_prescribe_velocity = .false.
+ config_prescribe_thickness = .false.
+ config_include_KE_vertex = .false.
+ config_check_tracer_monotonicity = .false.
+ config_disable_h_all_tend = .false.
+ config_disable_h_hadv = .false.
+ config_disable_h_vadv = .false.
+ config_disable_u_all_tend = .false.
+ config_disable_u_coriolis = .false.
+ config_disable_u_pgrad = .false.
+ config_disable_u_hmix = .false.
+ config_disable_u_windstress = .false.
+ config_disable_u_vmix = .false.
+ config_disable_u_vadv = .false.
+ config_disable_tr_all_tend = .false.
+ config_disable_tr_adv = .false.
+ config_disable_tr_hmix = .false.
+ config_disable_tr_vmix = .false.
+/
Index: branches/mpas_cdg_advection/src/core_ocean
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean 2013-02-13 21:29:30 UTC (rev 2473)
Property changes on: branches/mpas_cdg_advection/src/core_ocean
___________________________________________________________________
Modified: svn:mergeinfo
## -4,12 +4,15 ##
/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/cesm_coupling/src/core_ocean:2147-2344
+/branches/ocean_projects/diagnostics_revision/src/core_ocean:2439-2462
/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/linear_eos/src/core_ocean:2435-2437
/branches/ocean_projects/monotonic_advection/src/core_ocean:1499-1640
/branches/ocean_projects/monthly_forcing/src/core_ocean:1810-1867
+/branches/ocean_projects/namelist_cleanup/src/core_ocean:2319-2414
/branches/ocean_projects/option3_b4b_test/src/core_ocean:2201-2231
/branches/ocean_projects/partial_bottom_cells/src/core_ocean:2172-2226
/branches/ocean_projects/restart_reproducibility/src/core_ocean:2239-2272
## -27,3 +30,4 ##
/branches/omp_blocks/openmp_test/src/core_ocean_elements:2161-2201
/branches/source_renaming/src/core_ocean:1082-1113
/branches/time_manager/src/core_ocean:924-962
+/trunk/mpas/src/core_ocean:2390-2472
\ No newline at end of property
Modified: branches/mpas_cdg_advection/src/core_ocean/Makefile
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/Makefile 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/Makefile 2013-02-13 21:29:30 UTC (rev 2473)
@@ -3,59 +3,61 @@
OBJS = mpas_ocn_mpas_core.o \
mpas_ocn_test_cases.o \
mpas_ocn_advection.o \
- mpas_ocn_thick_hadv.o \
- mpas_ocn_thick_vadv.o \
- mpas_ocn_gm.o \
- mpas_ocn_vel_coriolis.o \
- 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 \
- mpas_ocn_vel_forcing_bottomdrag.o \
- mpas_ocn_vel_forcing_rayleigh.o \
- mpas_ocn_vel_pressure_grad.o \
- mpas_ocn_tracer_vadv.o \
- mpas_ocn_tracer_vadv_spline.o \
- mpas_ocn_tracer_vadv_spline2.o \
- mpas_ocn_tracer_vadv_spline3.o \
- mpas_ocn_tracer_vadv_stencil.o \
- mpas_ocn_tracer_vadv_stencil2.o \
- mpas_ocn_tracer_vadv_stencil3.o \
- mpas_ocn_tracer_vadv_stencil4.o \
- mpas_ocn_tracer_hadv.o \
- mpas_ocn_tracer_hadv2.o \
- mpas_ocn_tracer_hadv3.o \
- mpas_ocn_tracer_hadv4.o \
- mpas_ocn_tracer_hmix.o \
- mpas_ocn_tracer_hmix_del2.o \
- mpas_ocn_tracer_hmix_del4.o \
- mpas_ocn_vmix.o \
- mpas_ocn_vmix_coefs_const.o \
- mpas_ocn_vmix_coefs_rich.o \
- mpas_ocn_vmix_coefs_tanh.o \
- mpas_ocn_restoring.o \
- mpas_ocn_tendency.o \
- mpas_ocn_tracer_advection.o \
- mpas_ocn_tracer_advection_std.o \
- mpas_ocn_tracer_advection_std_hadv.o \
- mpas_ocn_tracer_advection_std_vadv.o \
- mpas_ocn_tracer_advection_std_vadv2.o \
- mpas_ocn_tracer_advection_std_vadv3.o \
- mpas_ocn_tracer_advection_std_vadv4.o \
- mpas_ocn_tracer_advection_mono.o \
- mpas_ocn_tracer_advection_helpers.o \
+ mpas_ocn_thick_hadv.o \
+ mpas_ocn_thick_vadv.o \
+ mpas_ocn_gm.o \
+ mpas_ocn_vel_coriolis.o \
+ mpas_ocn_vel_vadv.o \
+ mpas_ocn_vel_hmix.o \
+ 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 \
+ mpas_ocn_vel_pressure_grad.o \
+ mpas_ocn_tracer_vadv.o \
+ mpas_ocn_tracer_vadv_spline.o \
+ mpas_ocn_tracer_vadv_spline2.o \
+ mpas_ocn_tracer_vadv_spline3.o \
+ mpas_ocn_tracer_vadv_stencil.o \
+ mpas_ocn_tracer_vadv_stencil2.o \
+ mpas_ocn_tracer_vadv_stencil3.o \
+ mpas_ocn_tracer_vadv_stencil4.o \
+ mpas_ocn_tracer_hadv.o \
+ mpas_ocn_tracer_hadv2.o \
+ mpas_ocn_tracer_hadv3.o \
+ mpas_ocn_tracer_hadv4.o \
+ mpas_ocn_tracer_hmix.o \
+ mpas_ocn_tracer_hmix_del2.o \
+ mpas_ocn_tracer_hmix_del4.o \
+ mpas_ocn_vmix.o \
+ mpas_ocn_vmix_coefs_const.o \
+ mpas_ocn_vmix_coefs_rich.o \
+ mpas_ocn_vmix_coefs_tanh.o \
+ mpas_ocn_restoring.o \
+ mpas_ocn_tendency.o \
+ mpas_ocn_diagnostics.o \
+ mpas_ocn_tracer_advection.o \
+ mpas_ocn_tracer_advection_std.o \
+ mpas_ocn_tracer_advection_std_hadv.o \
+ mpas_ocn_tracer_advection_std_vadv.o \
+ mpas_ocn_tracer_advection_std_vadv2.o \
+ mpas_ocn_tracer_advection_std_vadv3.o \
+ mpas_ocn_tracer_advection_std_vadv4.o \
+ mpas_ocn_tracer_advection_mono.o \
+ mpas_ocn_tracer_advection_helpers.o \
mpas_ocn_time_integration.o \
mpas_ocn_time_integration_rk4.o \
mpas_ocn_time_integration_split.o \
- mpas_ocn_equation_of_state.o \
- mpas_ocn_equation_of_state_jm.o \
- mpas_ocn_equation_of_state_linear.o \
+ mpas_ocn_equation_of_state.o \
+ mpas_ocn_equation_of_state_jm.o \
+ mpas_ocn_equation_of_state_linear.o \
+ mpas_ocn_diagnostics.o \
mpas_ocn_global_diagnostics.o \
- mpas_ocn_time_average.o \
- mpas_ocn_monthly_forcing.o
+ mpas_ocn_time_average.o \
+ mpas_ocn_monthly_forcing.o
all: core_hyd
@@ -68,12 +70,14 @@
mpas_ocn_time_integration.o: mpas_ocn_time_integration_rk4.o mpas_ocn_time_integration_split.o
-mpas_ocn_time_integration_rk4.o:
+mpas_ocn_time_integration_rk4.o: mpas_ocn_tendency.o mpas_ocn_diagnostics.o
-mpas_ocn_time_integration_split.o:
+mpas_ocn_time_integration_split.o: mpas_ocn_tendency.o mpas_ocn_diagnostics.o
mpas_ocn_tendency.o: mpas_ocn_time_average.o
+mpas_ocn_diagnostics.o: mpas_ocn_time_average.o
+
mpas_ocn_global_diagnostics.o:
mpas_ocn_time_average.o:
@@ -172,61 +176,61 @@
mpas_ocn_monthly_forcing.o:
-mpas_ocn_mpas_core.o: mpas_ocn_mpas_core.o \
- mpas_ocn_test_cases.o \
- mpas_ocn_advection.o \
- mpas_ocn_thick_hadv.o \
+mpas_ocn_mpas_core.o: mpas_ocn_test_cases.o \
+ mpas_ocn_advection.o \
+ mpas_ocn_thick_hadv.o \
mpas_ocn_gm.o \
- mpas_ocn_thick_vadv.o \
- mpas_ocn_vel_coriolis.o \
- mpas_ocn_vel_vadv.o \
- 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 \
- mpas_ocn_vel_forcing_bottomdrag.o \
- mpas_ocn_vel_pressure_grad.o \
- mpas_ocn_tracer_vadv.o \
- mpas_ocn_tracer_vadv_spline.o \
- mpas_ocn_tracer_vadv_spline2.o \
- mpas_ocn_tracer_vadv_spline3.o \
- mpas_ocn_tracer_vadv_stencil.o \
- mpas_ocn_tracer_vadv_stencil2.o \
- mpas_ocn_tracer_vadv_stencil3.o \
- mpas_ocn_tracer_vadv_stencil4.o \
- mpas_ocn_tracer_hadv.o \
- mpas_ocn_tracer_hadv2.o \
- mpas_ocn_tracer_hadv3.o \
- mpas_ocn_tracer_hadv4.o \
- mpas_ocn_tracer_hmix.o \
- mpas_ocn_tracer_hmix_del2.o \
- mpas_ocn_tracer_hmix_del4.o \
- mpas_ocn_vmix.o \
- mpas_ocn_vmix_coefs_const.o \
- mpas_ocn_vmix_coefs_rich.o \
- mpas_ocn_vmix_coefs_tanh.o \
- mpas_ocn_restoring.o \
- mpas_ocn_tracer_advection.o \
- mpas_ocn_tracer_advection_std.o \
- mpas_ocn_tracer_advection_std_hadv.o \
- mpas_ocn_tracer_advection_std_vadv.o \
- mpas_ocn_tracer_advection_std_vadv2.o \
- mpas_ocn_tracer_advection_std_vadv3.o \
- mpas_ocn_tracer_advection_std_vadv4.o \
- mpas_ocn_tracer_advection_mono.o \
- mpas_ocn_tracer_advection_helpers.o \
- mpas_ocn_tendency.o \
- mpas_ocn_time_integration.o \
- mpas_ocn_time_integration_rk4.o \
- mpas_ocn_time_integration_split.o \
- mpas_ocn_equation_of_state.o \
- mpas_ocn_equation_of_state_jm.o \
- mpas_ocn_equation_of_state_linear.o \
- mpas_ocn_global_diagnostics.o \
- mpas_ocn_time_average.o \
- mpas_ocn_monthly_forcing.o
+ mpas_ocn_thick_vadv.o \
+ mpas_ocn_vel_coriolis.o \
+ 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 \
+ mpas_ocn_vel_forcing_bottomdrag.o \
+ mpas_ocn_vel_pressure_grad.o \
+ mpas_ocn_tracer_vadv.o \
+ mpas_ocn_tracer_vadv_spline.o \
+ mpas_ocn_tracer_vadv_spline2.o \
+ mpas_ocn_tracer_vadv_spline3.o \
+ mpas_ocn_tracer_vadv_stencil.o \
+ mpas_ocn_tracer_vadv_stencil2.o \
+ mpas_ocn_tracer_vadv_stencil3.o \
+ mpas_ocn_tracer_vadv_stencil4.o \
+ mpas_ocn_tracer_hadv.o \
+ mpas_ocn_tracer_hadv2.o \
+ mpas_ocn_tracer_hadv3.o \
+ mpas_ocn_tracer_hadv4.o \
+ mpas_ocn_tracer_hmix.o \
+ mpas_ocn_tracer_hmix_del2.o \
+ mpas_ocn_tracer_hmix_del4.o \
+ mpas_ocn_vmix.o \
+ mpas_ocn_vmix_coefs_const.o \
+ mpas_ocn_vmix_coefs_rich.o \
+ mpas_ocn_vmix_coefs_tanh.o \
+ mpas_ocn_restoring.o \
+ mpas_ocn_tracer_advection.o \
+ mpas_ocn_tracer_advection_std.o \
+ mpas_ocn_tracer_advection_std_hadv.o \
+ mpas_ocn_tracer_advection_std_vadv.o \
+ mpas_ocn_tracer_advection_std_vadv2.o \
+ mpas_ocn_tracer_advection_std_vadv3.o \
+ mpas_ocn_tracer_advection_std_vadv4.o \
+ mpas_ocn_tracer_advection_mono.o \
+ mpas_ocn_tracer_advection_helpers.o \
+ mpas_ocn_tendency.o \
+ mpas_ocn_diagnostics.o \
+ mpas_ocn_time_integration.o \
+ mpas_ocn_time_integration_rk4.o \
+ mpas_ocn_time_integration_split.o \
+ mpas_ocn_equation_of_state.o \
+ mpas_ocn_equation_of_state_jm.o \
+ mpas_ocn_equation_of_state_linear.o \
+ mpas_ocn_global_diagnostics.o \
+ mpas_ocn_time_average.o \
+ mpas_ocn_monthly_forcing.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/mpas_cdg_advection/src/core_ocean/Registry
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/Registry 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/Registry 2013-02-13 21:29:30 UTC (rev 2473)
@@ -1,103 +1,149 @@
%
% namelist type namelist_record name default_value
%
-namelist integer sw_model config_test_case 5
-namelist character sw_model config_time_integration RK4
-namelist logical sw_model config_rk_filter_btr_mode false
-namelist real sw_model config_dt 172.8
-namelist character sw_model config_calendar_type 360day
-namelist character sw_model config_start_time 0000-01-01_00:00:00
-namelist character sw_model config_stop_time none
-namelist character sw_model config_run_duration none
-namelist integer sw_model config_stats_interval 100
-namelist logical sw_model config_initial_stats false
-namelist logical sw_model config_prescribe_velocity false
-namelist logical sw_model config_prescribe_thickness false
-namelist integer sw_model config_num_halos 3
-namelist integer sw_model config_ncouple_per_day 1
-namelist character io config_input_name grid.nc
-namelist character io config_output_name output.nc
-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 logical time_management config_do_restart .false.
+namelist character time_management config_start_time '0000-01-01_00:00:00'
+namelist character time_management config_stop_time 'none'
+namelist character time_management config_run_duration '0_06:00:00'
+namelist character time_management config_calendar_type '360day'
+
+namelist character io config_input_name 'grid.nc'
+namelist character io config_output_name 'output.nc'
+namelist character io config_restart_name 'restart.nc'
+namelist character io config_restart_interval '0_06:00:00'
+namelist character io config_output_interval '0_06:00:00'
+namelist character io config_stats_interval '0_01:00:00'
+namelist logical io config_write_stats_on_startup .true.
+namelist logical io config_write_output_on_startup .true.
+namelist integer io config_frames_per_outfile 1000
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 real time_integration config_dt 3000.0
+namelist character time_integration config_time_integrator 'split_explicit'
+
+namelist integer grid config_num_halos 3
+namelist character grid config_vert_coord_movement 'uniform_stretching'
+namelist character grid config_alter_ICs_for_pbcs 'zlevel_pbcs_off'
+namelist real grid config_min_pbc_fraction 0.10
+namelist logical grid config_check_ssh_consistency .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 character decomposition config_proc_decomp_file_prefix graph.info.part.
-namelist logical restart config_do_restart false
-namelist character restart config_restart_interval none
-namelist character grid config_vert_grid_type isopycnal
-namelist character grid config_pressure_type pressure
-namelist real grid config_rho0 1028
-namelist logical grid config_enforce_zstar_at_restart false
-namelist 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
-namelist integer split_explicit_ts config_n_bcl_iter_end 2
-namelist integer split_explicit_ts config_n_btr_subcycles 20
-namelist integer split_explicit_ts config_n_btr_cor_iter 2
-namelist logical split_explicit_ts config_u_correction true
-namelist logical split_explicit_ts config_filter_btr_mode false
-namelist integer split_explicit_ts config_btr_subcycle_loop_factor 2
-namelist real split_explicit_ts config_btr_gam1_uWt1 0.5
-namelist real split_explicit_ts config_btr_gam2_SSHWt1 1.0
-namelist real split_explicit_ts config_btr_gam3_uWt2 1.0
-namelist logical split_explicit_ts config_btr_solve_SSH2 false
-namelist logical hmix config_h_ScaleWithMesh false
-namelist real hmix config_h_mom_eddy_visc2 0.0
-namelist real hmix config_h_mom_eddy_visc4 0.0
-namelist logical hmix config_visc_vorticity_term true
-namelist real hmix config_visc_vorticity_visc2_scale 1.0
-namelist real hmix config_visc_vorticity_visc4_scale 1.0
-namelist logical hmix config_include_KE_vertex false
-namelist real hmix config_h_tracer_eddy_diff2 0.0
-namelist real hmix config_h_tracer_eddy_diff4 0.0
-namelist real hmix config_h_kappa 0.0
-namelist real hmix config_h_kappa_q 0.0
-namelist logical hmix config_rayleigh_friction false
-namelist real hmix config_rayleigh_damping_coeff 0.0
-namelist real hmix config_apvm_scale_factor 0.0
-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 decomposition config_explicit_proc_decomp .false.
+namelist character decomposition config_proc_decomp_file_prefix 'graph.info.part.'
+
+namelist logical hmix config_hmix_ScaleWithMesh .false.
+namelist logical hmix config_visc_vorticity_term .true.
+namelist real hmix config_apvm_scale_factor 0.0
+
+namelist logical hmix_del2 config_use_mom_del2 .false.
+namelist logical hmix_del2 config_use_tracer_del2 .false.
+namelist real hmix_del2 config_mom_del2 0.0
+namelist real hmix_del2 config_tracer_del2 0.0
+namelist real hmix_del2 config_vorticity_del2_scale 1.0
+
+namelist logical hmix_del4 config_use_mom_del4 .true.
+namelist logical hmix_del4 config_use_tracer_del4 .false.
+namelist real hmix_del4 config_mom_del4 5.0e13
+namelist real hmix_del4 config_tracer_del4 0.0
+namelist real hmix_del4 config_vorticity_del4_scale 1.0
+
+namelist logical hmix_Leith config_use_Leith_del2 .false.
+namelist real hmix_Leith config_Leith_parameter 1.0
+namelist real hmix_Leith config_Leith_dx 15000.0
+namelist real hmix_Leith config_Leith_visc2_max 2.5e3
+
+namelist real standard_GM config_h_kappa 0.0
+namelist real standard_GM config_h_kappa_q 0.0
+
+namelist logical Rayleigh_damping config_Rayleigh_friction .false.
+namelist real Rayleigh_damping config_Rayleigh_damping_coeff 0.0
+
+namelist logical vmix config_implicit_vertical_mix .true.
namelist real vmix config_convective_visc 1.0
namelist real vmix config_convective_diff 1.0
-namelist real vmix config_bottom_drag_coeff 1.0e-3
+
+namelist logical vmix_const config_use_const_visc .false.
+namelist logical vmix_const config_use_const_diff .false.
namelist real vmix_const config_vert_visc 2.5e-4
namelist real vmix_const config_vert_diff 2.5e-5
+
+namelist logical vmix_rich config_use_rich_visc .true.
+namelist logical vmix_rich config_use_rich_diff .true.
namelist real vmix_rich config_bkrd_vert_visc 1.0e-4
namelist real vmix_rich config_bkrd_vert_diff 1.0e-5
namelist real vmix_rich config_rich_mix 0.005
+
+namelist logical vmix_tanh config_use_tanh_visc .false.
+namelist logical vmix_tanh config_use_tanh_diff .false.
namelist real vmix_tanh config_max_visc_tanh 2.5e-1
namelist real vmix_tanh config_min_visc_tanh 1.0e-4
namelist real vmix_tanh config_max_diff_tanh 2.5e-2
namelist real vmix_tanh config_min_diff_tanh 1.0e-5
namelist real vmix_tanh config_zMid_tanh -100
namelist real vmix_tanh config_zWidth_tanh 100
-namelist character eos config_eos_type linear
-namelist character advection config_vert_tracer_adv stencil
-namelist integer advection config_vert_tracer_adv_order 4
-namelist integer advection config_horiz_tracer_adv_order 2
-namelist integer advection config_thickness_adv_order 2
+
+namelist logical forcing config_use_monthly_forcing .false.
+namelist logical forcing config_restoreTS .false.
+namelist real forcing config_restoreT_timescale 90.0
+namelist real forcing config_restoreS_timescale 90.0
+
+namelist character advection config_vert_tracer_adv 'stencil'
+namelist integer advection config_vert_tracer_adv_order 3
+namelist integer advection config_horiz_tracer_adv_order 3
namelist real advection config_coef_3rd_order 0.25
-namelist logical advection config_monotonic false
-namelist logical advection config_check_monotonicity false
-namelist logical restore config_restoreTS false
-namelist real restore config_restoreT_timescale 90.0
-namelist real restore config_restoreS_timescale 90.0
-namelist logical restore config_use_monthly_forcing false
+namelist logical advection config_monotonic .true.
+namelist real bottom_drag config_bottom_drag_coeff 1.0e-3
+
+namelist character pressure_gradient config_pressure_gradient_type 'pressure_and_zmid'
+namelist real pressure_gradient config_rho0 1014.65
+
+namelist character eos config_eos_type 'jm'
+
+namelist real eos_linear config_eos_linear_alpha 2.55e-1
+namelist real eos_linear config_eos_linear_beta 7.64e-1
+namelist real eos_linear config_eos_linear_Tref 19.0
+namelist real eos_linear config_eos_linear_Sref 35.0
+namelist real eos_linear config_eos_linear_rhoref 1025.022
+
+namelist integer split_explicit_ts config_n_ts_iter 2
+namelist integer split_explicit_ts config_n_bcl_iter_beg 1
+namelist integer split_explicit_ts config_n_bcl_iter_mid 2
+namelist integer split_explicit_ts config_n_bcl_iter_end 2
+namelist integer split_explicit_ts config_n_btr_subcycles 20
+namelist integer split_explicit_ts config_n_btr_cor_iter 2
+namelist logical split_explicit_ts config_u_correction .true.
+namelist integer split_explicit_ts config_btr_subcycle_loop_factor 2
+namelist real split_explicit_ts config_btr_gam1_uWt1 0.5
+namelist real split_explicit_ts config_btr_gam2_SSHWt1 1.0
+namelist real split_explicit_ts config_btr_gam3_uWt2 1.0
+namelist logical split_explicit_ts config_btr_solve_SSH2 .false.
+
+namelist integer sw_model config_test_case 0
+
+namelist logical debug config_check_zlevel_consistency .false.
+namelist logical debug config_filter_btr_mode .false.
+namelist logical debug config_prescribe_velocity .false.
+namelist logical debug config_prescribe_thickness .false.
+namelist logical debug config_include_KE_vertex .false.
+namelist logical debug config_check_tracer_monotonicity .false.
+namelist logical debug config_disable_h_all_tend .false.
+namelist logical debug config_disable_h_hadv .false.
+namelist logical debug config_disable_h_vadv .false.
+namelist logical debug config_disable_u_all_tend .false.
+namelist logical debug config_disable_u_coriolis .false.
+namelist logical debug config_disable_u_pgrad .false.
+namelist logical debug config_disable_u_hmix .false.
+namelist logical debug config_disable_u_windstress .false.
+namelist logical debug config_disable_u_vmix .false.
+namelist logical debug config_disable_u_vadv .false.
+namelist logical debug config_disable_tr_all_tend .false.
+namelist logical debug config_disable_tr_adv .false.
+namelist logical debug config_disable_tr_hmix .false.
+namelist logical debug config_disable_tr_vmix .false.
+
%
% dim type name_in_file name_in_code
%
@@ -198,7 +244,7 @@
% Arrays required for reconstruction of velocity field
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
-% Arrays for z-level version of mpas-ocean
+% Arrays for non-isopycnal version of mpas-ocean
var persistent integer maxLevelCell ( nCells ) 0 iro maxLevelCell mesh - -
var persistent integer maxLevelEdgeTop ( nEdges ) 0 - maxLevelEdgeTop mesh - -
var persistent integer maxLevelEdgeBot ( nEdges ) 0 - maxLevelEdgeBot mesh - -
@@ -207,7 +253,7 @@
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 - -
+var persistent real vertCoordMovementWeights ( nVertLevels ) 0 iro vertCoordMovementWeights mesh - -
% Boundary conditions and masks
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 - boundaryEdge mesh - -
@@ -248,7 +294,7 @@
var persistent real tend_tracer1Moms ( nVertBasis nHorizBasis nVertLevels nCells Time ) 1 - tracer1Moms tend tracerMoms testingMoms
% state variables for Split Explicit timesplitting
-var persistent real uBtr ( nEdges Time ) 2 r 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 - -
@@ -293,7 +339,9 @@
var persistent real MontPot ( nVertLevels nCells Time ) 2 - MontPot state - -
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
+var persistent real vertVelocityTop ( nVertLevelsP1 nCells Time ) 2 - vertVelocityTop state - -
var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
+var persistent real BruntVaisalaFreqTop ( nVertLevels nCells Time ) 2 o BruntVaisalaFreqTop state - -
var persistent real viscosity ( nVertLevels nEdges Time ) 2 o viscosity state - -
% Other diagnostic variables: neither read nor written to any files
@@ -326,6 +374,7 @@
var persistent real acc_uReconstructMeridionalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridionalVar state - -
var persistent real acc_u ( nVertLevels nEdges Time ) 2 o acc_u state - -
var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 o acc_uVar state - -
+var persistent real acc_vertVelocityTop ( nVertLevelsP1 nCells Time ) 2 o acc_vertVelocityTop state - -
% Sign fields, for openmp and bit reproducibility without branching statements.
var persistent integer edgeSignOnCell ( maxEdges nCells ) 0 - edgeSignOnCell mesh - -
Copied: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_diagnostics.F (from rev 2472, trunk/mpas/src/core_ocean/mpas_ocn_diagnostics.F)
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_diagnostics.F (rev 0)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_diagnostics.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -0,0 +1,845 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! ocn_diagnostics
+!
+!> \brief MPAS ocean diagnostics driver
+!> \author Mark Petersen
+!> \date 23 September 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routines for computing
+!> diagnostic variables, and other quantities such as wTop.
+!
+!-----------------------------------------------------------------------
+
+module ocn_diagnostics
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_constants
+ use mpas_timer
+
+ use ocn_gm
+ use ocn_equation_of_state
+
+ implicit none
+ private
+ save
+
+ type (timer_node), pointer :: diagEOSTimer
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: ocn_diagnostic_solve, &
+ ocn_wtop, &
+ ocn_fuperp, &
+ ocn_filter_btr_mode_u, &
+ ocn_filter_btr_mode_tend_u, &
+ ocn_diagnostics_init
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ integer :: ke_cell_flag, ke_vertex_flag
+ real (kind=RKIND) :: coef_3rd_order, fCoef
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine ocn_diagnostic_solve
+!
+!> \brief Computes diagnostic variables
+!> \author Mark Petersen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the diagnostic variables for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_diagnostic_solve(dt, s, grid)!{{{
+ implicit none
+
+ real (kind=RKIND), intent(in) :: dt !< Input: Time step
+ type (state_type), intent(inout) :: s !< Input/Output: State information
+ type (mesh_type), intent(in) :: grid !< Input: Grid information
+
+
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
+ integer :: boundaryMask, velMask, nCells, nEdges, nVertices, nVertLevels, vertexDegree, err
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
+ 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, coef
+
+ real (kind=RKIND), dimension(:), allocatable:: pTop, div_hu
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ 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, &
+ Vor_edge, Vor_vertex, Vor_cell, gradVor_n, gradVor_t, divergence, &
+ rho, rhoDisplaced, temperature, salinity, kev, kevc, uBolusGM, uTransport, &
+ vertVelocityTop, BruntVaisalaFreqTop
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers, deriv_two
+ character :: c1*6
+
+ h => s % h % array
+ u => s % u % array
+ uTransport => s % uTransport % array
+ uBolusGM => s % uBolusGM % array
+ v => s % v % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ 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
+ rho => s % rho % array
+ rhoDisplaced=> s % rhoDisplaced % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+ zMid => s % zMid % array
+ ssh => s % ssh % array
+ tracers => s % tracers % array
+ vertVelocityTop => s % vertVelocityTop % array
+ BruntVaisalaFreqTop => s % BruntVaisalaFreqTop % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ 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
+ bottomDepth => grid % bottomDepth % array
+ fVertex => grid % fVertex % 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
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
+
+ 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
+ !
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
+
+ ! initialize h_edge to avoid divide by zero and NaN problems.
+ h_edge = -1.0e34
+ coef_3rd_order = config_coef_3rd_order
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
+ end do
+
+ !
+ ! set the velocity and height at dummy address
+ ! used -1e34 so error clearly occurs if these values are used.
+ !
+ u(:,nEdges+1) = -1e34
+ h(:,nCells+1) = -1e34
+ tracers(s % index_temperature,:,nCells+1) = -1e34
+ tracers(s % index_salinity,:,nCells+1) = -1e34
+
+ circulation(:,:) = 0.0
+ vorticity(:,:) = 0.0
+ divergence(:,:) = 0.0
+ vertVelocityTop(:,:)=0.0
+ ke(:,:) = 0.0
+ v(:,:) = 0.0
+ 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)
+
+ 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
+
+ allocate(div_hu(nVertLevels))
+ do iCell = 1, nCells
+ div_hu(:) = 0.0
+ 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
+
+ divergence(k, iCell) = divergence(k, iCell) - edgeSignOnCell(i, iCell) * r_tmp
+ div_hu(k) = div_hu(k) - h_edge(k, iEdge) * edgeSignOnCell(i, iCell) * r_tmp
+ ke(k, iCell) = ke(k, iCell) + 0.25 * r_tmp * dcEdge(iEdge) * u(k,iEdge)
+ end do
+ end do
+ ! Vertical velocity at bottom (maxLevelCell(iCell)+1) is zero, initialized above.
+ do k=maxLevelCell(iCell),1,-1
+ vertVelocityTop(k,iCell) = vertVelocityTop(k+1,iCell) - div_hu(k)
+ end do
+ end do
+ deallocate(div_hu)
+
+ do iEdge=1,nEdges
+ ! Compute v (tangential) velocities
+ do i=1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(i,iEdge)
+ ! mrp 101115 note: in order to include flux boundary conditions,
+ ! the following loop may need to change to maxLevelEdgeBot
+ do k = 1,maxLevelEdgeTop(iEdge)
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
+ end do
+
+ end do
+
+ !
+ ! Compute kinetic energy in each vertex
+ !
+ kev(:,:) = 0.0; kevc(:,:) = 0.0
+ 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 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
+ !
+ do iCell=1,nCells*ke_vertex_flag
+ do k=1,nVertLevels
+ ke(k,iCell) = 5.0/8.0*ke(k,iCell) + 3.0/8.0*kevc(k,iCell)
+ end do
+ end do
+
+ !
+ ! Compute ke on cell edges at velocity locations for quadratic bottom drag.
+ !
+ ! mrp 101025 efficiency note: we could get rid of ke_edge completely by
+ ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
+ end do
+ end do
+
+ !
+ ! Compute height at vertices, pv at vertices, and average pv to edge locations
+ ! ( this computes Vor_vertex at all vertices bounding real cells and distance-1 ghost cells )
+ !
+ do iVertex = 1,nVertices
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
+ do k=1,maxLevelVertexBot(iVertex)
+ h_vertex = 0.0
+ do i=1,vertexDegree
+ h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
+ end do
+ h_vertex = h_vertex * invAreaTri1
+
+ Vor_vertex(k,iVertex) = (fCoef*fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
+ end do
+ end do
+
+ Vor_cell(:,:) = 0.0
+ Vor_edge(:,:) = 0.0
+ 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
+
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+
+ 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
+
+ do iEdge = 1,nEdges
+ cell1 = cellsOnEdge(1, iEdge)
+ cell2 = cellsOnEdge(2, iEdge)
+ vertex1 = verticesOnedge(1, iEdge)
+ vertex2 = verticesOnedge(2, iEdge)
+
+ invLength = 1.0 / dcEdge(iEdge)
+ ! Compute gradient of PV in normal direction
+ ! ( this computes gradVor_n for all edges bounding real cells )
+ do k=1,maxLevelEdgeTop(iEdge)
+ gradVor_n(k,iEdge) = (Vor_cell(k,cell2) - Vor_cell(k,cell1)) * invLength
+ enddo
+
+ invLength = 1.0 / dvEdge(iEdge)
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradVor_t at all edges bounding real cells and distance-1 ghost cells )
+ do k = 1,maxLevelEdgeBot(iEdge)
+ gradVor_t(k,iEdge) = (Vor_vertex(k,vertex2) - Vor_vertex(k,vertex1)) * invLength
+ enddo
+
+ enddo
+
+ !
+ ! Modify PV edge with upstream bias.
+ !
+ do iEdge = 1,nEdges
+ do k = 1,maxLevelEdgeBot(iEdge)
+ Vor_edge(k,iEdge) = Vor_edge(k,iEdge) &
+ - config_apvm_scale_factor * dt* ( u(k,iEdge) * gradVor_n(k,iEdge) &
+ + v(k,iEdge) * gradVor_t(k,iEdge) )
+ enddo
+ enddo
+
+ !
+ ! equation of state
+ !
+ ! For an isopycnal model, density should remain constant.
+ ! For zlevel, calculate in-situ density
+ if (config_vert_coord_movement.ne.'isopycnal') then
+ call mpas_timer_start("equation of state", .false., diagEOSTimer)
+
+ ! compute in-place density
+ call ocn_equation_of_state_rho(s, grid, 0, 'relative', err)
+
+ ! compute rhoDisplaced, the potential density referenced to the top layer
+ call ocn_equation_of_state_rho(s, grid, 1, 'relative', err)
+
+ call mpas_timer_stop("equation of state", diagEOSTimer)
+ endif
+
+ !
+ ! Pressure
+ ! This section must be after computing rho
+ !
+ ! dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
+ if (config_pressure_gradient_type.eq.'MontgomeryPotential') then
+
+ ! For Isopycnal model.
+ ! Compute pressure at top of each layer, and then
+ ! Montgomery Potential.
+ allocate(pTop(nVertLevels))
+ do iCell=1,nCells
+
+ ! assume atmospheric pressure at the surface is zero for now.
+ 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, bottomDepth can be relative to top or bottom)
+ MontPot(1,iCell) = gravity &
+ * (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)
+
+ ! from delta M = p delta / rho
+ MontPot(k,iCell) = MontPot(k-1,iCell) &
+ + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
+ end do
+
+ end do
+ deallocate(pTop)
+
+ else
+
+ do iCell=1,nCells
+ ! pressure for generalized coordinates
+ ! assume atmospheric pressure at the surface is zero for now.
+ pressure(1,iCell) = rho(1,iCell)*gravity &
+ * 0.5*h(1,iCell)
+
+ do k=2,maxLevelCell(iCell)
+ pressure(k,iCell) = pressure(k-1,iCell) &
+ + 0.5*gravity*( rho(k-1,iCell)*h(k-1,iCell) &
+ + rho(k ,iCell)*h(k ,iCell))
+ end do
+
+ ! 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 bottomDepth is positive
+ ! and z-coordinates are negative below the surface.
+ k = maxLevelCell(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) &
+ + 0.5*( h(k+1,iCell) &
+ + h(k ,iCell))
+ end do
+
+ end do
+
+ endif
+
+ !
+ ! Brunt-Vaisala frequency
+ !
+ coef = -gravity/config_rho0
+ do iCell=1,nCells
+ BruntVaisalaFreqTop(1,iCell) = 0.0
+ do k=2,maxLevelCell(iCell)
+ BruntVaisalaFreqTop(k,iCell) = coef * (rhoDisplaced(k-1,iCell) - rhoDisplaced(k,iCell)) &
+ / (zMid(k-1,iCell) - zMid(k,iCell))
+ end do
+ end do
+
+ !
+ ! Sea Surface Height
+ !
+ do iCell=1,nCells
+ ! Start at the bottom where we know the depth, and go up.
+ ! 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) = - bottomDepth(iCell) + sum(h(1:maxLevelCell(iCell),iCell))
+
+ end do
+
+ !
+ ! Apply the GM closure as a bolus velocity
+ !
+ if (config_h_kappa .GE. epsilon(0D0)) then
+ call ocn_gm_compute_uBolus(s,grid)
+ else
+ ! mrp efficiency note: if uBolusGM is guaranteed to be zero, this can be removed.
+ uBolusGM = 0.0
+ end if
+
+ end subroutine ocn_diagnostic_solve!}}}
+
+!***********************************************************************
+!
+! routine ocn_wtop
+!
+!> \brief Computes vertical transport
+!> \author Mark Petersen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the vertical transport through the top of each
+!> cell.
+!
+!-----------------------------------------------------------------------
+ subroutine ocn_wtop(grid,h,h_edge,u,wTop, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! 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: transport
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(out) :: &
+ wTop !< Output: vertical transport at top of cell
+
+ 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, invAreaCell
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ dvEdge, areaCell, vertCoordMovementWeights
+ 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, edgeSignOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot, maxLevelVertexTop
+
+ 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
+ vertCoordMovementWeights => grid % vertCoordMovementWeights % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+
+
+ if (config_vert_coord_movement.eq.'isopycnal') then
+ ! set vertical transport to zero in isopycnal case
+ wTop=0.0
+ 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.
+ !
+
+ do iCell=1,nCells
+ div_hu(:) = 0.0
+ div_hu_btr = 0.0
+ hSum = 0.0
+ invAreaCell = 1.0 / areaCell(iCell)
+
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+
+ do k = 1, maxLevelEdgeBot(iEdge)
+ flux = h_edge(k, iEdge) * u(k, iEdge) * dvEdge(iEdge) * edgeSignOnCell(i, iCell) * invAreaCell
+ div_hu(k) = div_hu(k) - flux
+ div_hu_btr = div_hu_btr - flux
+ end do
+ end do
+
+ do k = 1, maxLevelCell(iCell)
+ h_tend_col(k) = - vertCoordMovementWeights(k) * h(k, iCell) * div_hu_btr
+ hSum = hSum + vertCoordMovementWeights(k) * h(k, iCell)
+ end do
+
+ if(hSum > 0.0) then
+ h_tend_col = h_tend_col / hSum
+ end if
+
+ ! Vertical transport 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) - h_tend_col(k)
+ end do
+ end do
+
+ deallocate(div_hu, h_tend_col)
+
+ end subroutine ocn_wtop!}}}
+
+!***********************************************************************
+!
+! routine ocn_fuperp
+!
+!> \brief Computes f u_perp
+!> \author Mark Petersen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes f u_perp for the ocean
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_fuperp(s, grid)!{{{
+ implicit none
+
+ type (state_type), intent(inout) :: s !< Input/Output: State information
+ type (mesh_type), intent(in) :: grid !< Input: Grid information
+
+! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
+! Some of these variables can be removed, but at a later time.
+ integer :: iEdge, cell1, cell2, eoe, i, j, k
+
+ integer :: nEdgesSolve
+ real (kind=RKIND), dimension(:), pointer :: fEdge
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, u, uBcl
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnEdge
+
+ call mpas_timer_start("ocn_fuperp")
+
+ u => s % u % array
+ uBcl => s % uBcl % array
+ weightsOnEdge => grid % weightsOnEdge % array
+ fEdge => grid % fEdge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+
+ fEdge => grid % fEdge % array
+
+ nEdgesSolve = grid % nEdgesSolve
+
+ !
+ ! Put f*uBcl^{perp} in u as a work variable
+ !
+ do iEdge=1,nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ u(k,iEdge) = 0.0
+ do j = 1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ u(k,iEdge) = u(k,iEdge) + weightsOnEdge(j,iEdge) * uBcl(k,eoe) * fEdge(eoe)
+ end do
+ end do
+ end do
+
+ call mpas_timer_stop("ocn_fuperp")
+
+ end subroutine ocn_fuperp!}}}
+
+!***********************************************************************
+!
+! routine ocn_filter_btr_mode_u
+!
+!> \brief filters barotropic mode out of the velocity variable.
+!> \author Mark Petersen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine filters barotropic mode out of the velocity variable.
+!
+!-----------------------------------------------------------------------
+ subroutine ocn_filter_btr_mode_u(s, grid)!{{{
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer :: iEdge, k, nEdges
+ real (kind=RKIND) :: vertSum, uhSum, hSum
+ real (kind=RKIND), dimension(:,:), pointer :: h_edge, u
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ call mpas_timer_start("ocn_filter_btr_mode_u")
+
+ u => s % u % array
+ h_edge => s % h_edge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ nEdges = grid % nEdges
+
+ do iEdge=1,nEdges
+
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
+ uhSum = h_edge(1,iEdge) * u(1,iEdge)
+ hSum = h_edge(1,iEdge)
+
+ do k=2,maxLevelEdgeTop(iEdge)
+ uhSum = uhSum + h_edge(k,iEdge) * u(k,iEdge)
+ hSum = hSum + h_edge(k,iEdge)
+ enddo
+
+ vertSum = uhSum/hSum
+ do k=1,maxLevelEdgeTop(iEdge)
+ u(k,iEdge) = u(k,iEdge) - vertSum
+ enddo
+ enddo ! iEdge
+
+ call mpas_timer_stop("ocn_filter_btr_mode_u")
+
+ end subroutine ocn_filter_btr_mode_u!}}}
+
+!***********************************************************************
+!
+! routine ocn_filter_btr_mode_tend_u
+!
+!> \brief ocn_filters barotropic mode out of the u tendency
+!> \author Mark Petersen
+!> \date 23 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine filters barotropic mode out of the u tendency.
+!
+!-----------------------------------------------------------------------
+ subroutine ocn_filter_btr_mode_tend_u(tend, s, grid)!{{{
+ implicit none
+
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer :: iEdge, k, nEdges
+ real (kind=RKIND) :: vertSum, uhSum, hSum
+ real (kind=RKIND), dimension(:,:), pointer :: h_edge, tend_u
+
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ call mpas_timer_start("ocn_filter_btr_mode_tend_u")
+
+ tend_u => tend % u % array
+ h_edge => s % h_edge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ nEdges = grid % nEdges
+
+ do iEdge=1,nEdges
+
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
+ uhSum = h_edge(1,iEdge) * tend_u(1,iEdge)
+ hSum = h_edge(1,iEdge)
+
+ do k=2,maxLevelEdgeTop(iEdge)
+ uhSum = uhSum + h_edge(k,iEdge) * tend_u(k,iEdge)
+ hSum = hSum + h_edge(k,iEdge)
+ enddo
+
+ vertSum = uhSum/hSum
+ do k=1,maxLevelEdgeTop(iEdge)
+ tend_u(k,iEdge) = tend_u(k,iEdge) - vertSum
+ enddo
+ enddo ! iEdge
+
+ call mpas_timer_stop("ocn_filter_btr_mode_tend_u")
+
+ end subroutine ocn_filter_btr_mode_tend_u!}}}
+
+!***********************************************************************
+!
+! routine ocn_diagnostics_init
+!
+!> \brief Initializes flags used within diagnostics routines.
+!> \author Mark Petersen
+!> \date 4 November 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes flags related to quantities computed within
+!> other diagnostics routines.
+!
+!-----------------------------------------------------------------------
+ subroutine ocn_diagnostics_init(err)!{{{
+ integer, intent(out) :: err !< Output: Error flag
+
+ err = 0
+
+ if(config_include_KE_vertex) then
+ ke_vertex_flag = 1
+ ke_cell_flag = 0
+ else
+ ke_vertex_flag = 0
+ ke_cell_flag = 1
+ endif
+
+ if (trim(config_time_integrator) == 'RK4') then
+ ! for RK4, PV is really PV = (eta+f)/h
+ fCoef = 1
+ elseif (trim(config_time_integrator) == 'split_explicit' &
+ .or.trim(config_time_integrator) == 'unsplit_explicit') then
+ ! for split explicit, PV is eta/h because f is added separately to the momentum forcing.
+ ! mrp temp, new should be:
+ fCoef = 0
+ ! old, for testing:
+ ! fCoef = 1
+ end if
+
+ end subroutine ocn_diagnostics_init!}}}
+
+!***********************************************************************
+
+end module ocn_diagnostics
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_equation_of_state.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_equation_of_state.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_equation_of_state.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -158,7 +158,7 @@
linearEos = .false.
jmEos = .false.
- if(config_vert_grid_type.ne.'isopycnal') then
+ if(config_vert_coord_movement.ne.'isopycnal') then
eosON = .true.
if (config_eos_type.eq.'linear') then
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_equation_of_state_linear.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_equation_of_state_linear.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_equation_of_state_linear.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -82,27 +82,21 @@
integer, intent(in) :: indexT, indexS
integer, intent(out) :: err
- real (kind=RKIND), parameter :: rho_ref = 1025.022 ! kg / m^3
- real (kind=RKIND), parameter :: alpha = 2.55e-1 ! kg / m^3 / K (dT/dRho)
- real (kind=RKIND), parameter :: beta = 7.64e-1 ! kg / m^3 / psu (dS/dRho)
- real (kind=RKIND), parameter :: T_ref = 19.0 ! K
- real (kind=RKIND), parameter :: S_ref = 35.0 ! psu
- real (kind=RKIND), parameter :: rho_prime_ref = rho_ref + alpha * T_ref - beta * S_ref
-
integer, dimension(:), pointer :: maxLevelCell
integer :: nCells, iCell, k
type (dm_info) :: dminfo
- maxLevelCell => grid % maxLevelCell % array
- nCells = grid % nCells
+ maxLevelCell => grid % maxLevelCell % array
+ nCells = grid % nCells
err = 0
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
! Linear equation of state
- ! rho = rho_ref - alpha * (T - T_ref) + beta * (S - S_ref)
- rho(k,iCell) = rho_prime_ref - alpha*tracers(indexT,k,iCell) + beta*tracers(indexS,k,iCell)
+ rho(k,iCell) = config_eos_linear_rhoref &
+ - config_eos_linear_alpha * (tracers(indexT,k,iCell)-config_eos_linear_Tref) &
+ + config_eos_linear_beta * (tracers(indexS,k,iCell)-config_eos_linear_Sref)
end do
end do
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_gm.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_gm.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_gm.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -30,7 +30,7 @@
contains
- subroutine ocn_gm_compute_uBolus(s, grid)
+ subroutine ocn_gm_compute_uBolus(s, grid)!{{{
implicit none
type(state_type), intent(inout) :: s
type(mesh_type), intent(in) :: grid
@@ -50,7 +50,7 @@
call ocn_gm_compute_hEddyFlux(s, grid)
- if (config_vert_grid_type .EQ. 'isopycnal') then
+ if (config_vert_coord_movement .EQ. 'isopycnal') then
do iEdge = 1, nEdges
do k = 1, maxLevelEdgeTop(iEdge)
@@ -65,10 +65,9 @@
end if
- end subroutine ocn_gm_compute_uBolus
+ end subroutine ocn_gm_compute_uBolus!}}}
-
- subroutine ocn_gm_compute_hEddyFlux(s, grid)
+ subroutine ocn_gm_compute_hEddyFlux(s, grid)!{{{
implicit none
type(state_type), intent(inout) :: s
type(mesh_type), intent(in) :: grid
@@ -90,7 +89,7 @@
hEddyFlux(:,:) = 0.0
- if (config_vert_grid_type .EQ. 'isopycnal') then
+ if (config_vert_coord_movement .EQ. 'isopycnal') then
do iEdge = 1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -104,12 +103,10 @@
end if
- end subroutine ocn_gm_compute_hEddyFlux
+ end subroutine ocn_gm_compute_hEddyFlux!}}}
+ subroutine ocn_get_h_kappa(s, grid)!{{{
-
- subroutine ocn_get_h_kappa(s, grid)
-
type (state_type), intent(inout) :: s
type (mesh_type), intent(in) :: grid
@@ -121,11 +118,10 @@
h_kappa(:,:) = config_h_kappa
- end subroutine ocn_get_h_kappa
+ end subroutine ocn_get_h_kappa!}}}
+ subroutine ocn_get_h_kappa_q(s, grid)!{{{
- subroutine ocn_get_h_kappa_q(s, grid)
-
type (state_type), intent(inout) :: s
type (mesh_type), intent(in) :: grid
@@ -137,6 +133,6 @@
h_kappa_q(:,:) = config_h_kappa_q
- end subroutine ocn_get_h_kappa_q
+ end subroutine ocn_get_h_kappa_q!}}}
end module ocn_gm
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_mpas_core.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_mpas_core.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -10,13 +10,18 @@
use ocn_test_cases
use ocn_time_integration
use ocn_tendency
+ use ocn_diagnostics
use ocn_monthly_forcing
+ use ocn_thick_hadv
+ use ocn_thick_vadv
+
use ocn_vel_pressure_grad
use ocn_vel_vadv
use ocn_vel_hmix
use ocn_vel_forcing
+ use ocn_vel_coriolis
use ocn_tracer_hadv
use ocn_tracer_vadv
@@ -64,9 +69,16 @@
! Initialize submodules before initializing blocks.
call ocn_timestep_init(err)
+ call ocn_thick_hadv_init(err_tmp)
+ err = ior(err, err_tmp)
+
+ call ocn_thick_vadv_init(err_tmp)
+ err = ior(err, err_tmp)
+
+ call ocn_vel_coriolis_init(err_tmp)
+ err = ior(err, err_tmp)
call ocn_vel_pressure_grad_init(err_tmp)
err = ior(err, err_tmp)
-
call ocn_vel_vadv_init(err_tmp)
err = ior(err, err_tmp)
call ocn_vel_hmix_init(err_tmp)
@@ -91,6 +103,8 @@
call ocn_tendency_init(err_tmp)
err = ior(err,err_tmp)
+ call ocn_diagnostics_init(err_tmp)
+ err = ior(err,err_tmp)
call mpas_ocn_tracer_advection_init(err_tmp)
err = ior(err,err_tmp)
@@ -106,38 +120,38 @@
if (.not. config_do_restart) call setup_sw_test_case(domain)
- if (config_vert_grid_type.ne.'isopycnal') call ocn_init_vert_coord(domain)
+ if (config_vert_coord_movement.ne.'isopycnal') call ocn_init_vert_coord(domain)
call ocn_compute_max_level(domain)
- if (config_enforce_zstar_at_restart) then
- call ocn_init_h_zstar(domain)
- endif
-
if (.not.config_do_restart) call ocn_init_split_timestep(domain)
- write (0,'(a,a10)') ' Vertical grid type is: ',config_vert_grid_type
+ write (0,'(a,a10)') ' Vertical coordinate movement is: ',config_vert_coord_movement
- 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
- write (0,*) ' Incorrect choice of config_vert_grid_type.'
+ if (config_vert_coord_movement.ne.'isopycnal'.and. &
+ config_vert_coord_movement.ne.'fixed'.and. &
+ config_vert_coord_movement.ne.'uniform_stretching'.and. &
+ config_vert_coord_movement.ne.'user_specified') then
+ write (0,*) ' Incorrect choice of config_vert_coord_movement.'
call mpas_dmpar_abort(dminfo)
endif
- write (0,'(a,a10)') ' Pressure type is: ',config_pressure_type
- if (config_pressure_type.ne.'pressure'.and. &
- config_pressure_type.ne.'MontgomeryPotential') then
- write (0,*) ' Incorrect choice of config_pressure_type.'
+ write (0,'(a,a)') ' Pressure type is: ',trim(config_pressure_gradient_type)
+ if (config_pressure_gradient_type.ne.'pressure_and_zmid'.and. &
+ config_pressure_gradient_type.ne.'MontgomeryPotential') then
+ write (0,*) ' Incorrect choice of config_pressure_gradient_type.'
call mpas_dmpar_abort(dminfo)
endif
+ if(config_vert_coord_movement .ne. 'isopycnal' .and. config_pressure_gradient_type .eq. 'MontgomeryPotential') then
+ write (0,*) ' Incorrect combination of config_vert_coord_movement and config_pressure_gradient_type'
+ call mpas_dmpar_abort(dminfo)
+ end if
+
if (config_filter_btr_mode.and. &
- config_vert_grid_type.ne.'zlevel')then
+ config_vert_coord_movement.ne.'fixed')then
write (0,*) 'filter_btr_mode has only been tested with'// &
- ' config_vert_grid_type=zlevel.'
+ ' config_vert_coord_movement=fixed.'
call mpas_dmpar_abort(dminfo)
endif
@@ -162,7 +176,7 @@
! input arguement into mpas_init. Ask about that later. For now, there will be
! no initial statistics write.
- if (config_initial_stats) then
+ if (config_write_stats_on_startup) then
call mpas_timer_start("global diagnostics", .false., globalDiagTimer)
call ocn_compute_global_diagnostics(domain, 1 , 0, dt)
call mpas_timer_stop("global diagnostics", globalDiagTimer)
@@ -230,11 +244,11 @@
!TODO: use this code if we desire to convert config_stats_interval to alarms
!(must also change config_stats_interval type to character)
! set stats alarm, if necessary
- !if (trim(config_stats_interval) /= "none") then
- ! call mpas_set_timeInterval(alarmTimeStep, timeString=config_stats_interval, ierr=ierr)
- ! alarmStartTime = startTime + alarmTimeStep
- ! call mpas_add_clock_alarm(clock, statsAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
- !end if
+ if (trim(config_stats_interval) /= "none") then
+ call mpas_set_timeInterval(alarmTimeStep, timeString=config_stats_interval, ierr=ierr)
+ alarmStartTime = startTime + alarmTimeStep
+ call mpas_add_clock_alarm(clock, statsAlarmID, alarmStartTime, alarmTimeStep, ierr=ierr)
+ end if
call mpas_get_time(curr_time=startTime, dateTimeString=startTimeStamp, ierr=ierr)
@@ -514,30 +528,31 @@
call ocn_timestep(domain, dt, timeStamp)
- if (config_stats_interval > 0) then
- if (mod(itimestep, config_stats_interval) == 0) then
- call mpas_timer_start("global diagnostics", .false., globalDiagTimer)
- call ocn_compute_global_diagnostics(domain, 2, itimestep, dt);
- call mpas_timer_stop("global diagnostics", globalDiagTimer)
- end if
- end if
+ !if (config_stats_interval > 0) then
+ ! if (mod(itimestep, config_stats_interval) == 0) then
+ ! call mpas_timer_start("global diagnostics", .false., globalDiagTimer)
+ ! call ocn_compute_global_diagnostics(domain, 2, itimestep, dt);
+ ! call mpas_timer_stop("global diagnostics", globalDiagTimer)
+ ! end if
+ !end if
!TODO: replace the above code block with this if we desire to convert config_stats_interval to use alarms
- !if (mpas_is_alarm_ringing(clock, statsAlarmID, ierr=ierr)) then
- ! call mpas_reset_clock_alarm(clock, statsAlarmID, ierr=ierr)
+ if (mpas_is_alarm_ringing(clock, statsAlarmID, ierr=ierr)) then
+ call mpas_reset_clock_alarm(clock, statsAlarmID, ierr=ierr)
- ! block_ptr => domain % blocklist
- ! if (associated(block_ptr % next)) then
- ! write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
- ! 'that there is only one block per processor.'
- ! end if
+! block_ptr => domain % blocklist
+! if (associated(block_ptr % next)) then
+! write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
+! 'that there is only one block per processor.'
+! end if
- ! call mpas_timer_start("global diagnostics")
- ! call ocn_compute_global_diagnostics(domain % dminfo, &
- ! block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
- ! timeStamp, dt)
- ! call mpas_timer_stop("global diagnostics")
- !end if
+ call mpas_timer_start("global diagnostics")
+ call ocn_compute_global_diagnostics(domain, 2, itimestep, dt);
+ ! call ocn_compute_global_diagnostics(domain % dminfo, &
+ ! block_ptr % state % time_levs(2) % state, block_ptr % mesh, &
+ ! timeStamp, dt)
+ call mpas_timer_stop("global diagnostics")
+ end if
end subroutine mpas_timestep!}}}
@@ -561,7 +576,7 @@
integer, dimension(:), pointer :: maxLevelCell
real (kind=RKIND), dimension(:), pointer :: refBottomDepth, &
- refBottomDepthTopOfCell, zstarWeight, hZLevel, bottomDepth
+ refBottomDepthTopOfCell, vertCoordMovementWeights, hZLevel, bottomDepth
real (kind=RKIND), dimension(:), allocatable :: minBottomDepth, minBottomDepthMid, zMidZLevel
real (kind=RKIND), dimension(:,:), pointer :: h
@@ -578,7 +593,7 @@
refBottomDepth => block % mesh % refBottomDepth % array
refBottomDepthTopOfCell => block % mesh % refBottomDepthTopOfCell % array
bottomDepth => block % mesh % bottomDepth % array
- zstarWeight => block % mesh % zstarWeight % array
+ vertCoordMovementWeights => block % mesh % vertCoordMovementWeights % array
hZLevel => block % mesh % hZLevel % array
maxLevelCell => block % mesh % maxLevelCell % array
@@ -602,28 +617,17 @@
refBottomDepthTopOfCell(k+1) = refBottomDepth(k)
end do
- ! Initialization of zstarWeights. This determines how SSH perturbations
+ ! Initialization of vertCoordMovementWeights. This determines how SSH perturbations
! are distributed throughout the column.
- if (config_vert_grid_type.eq.'zlevel') then
+ if (config_vert_coord_movement.eq.'fixed') then
- zstarWeight = 0.0
- zstarWeight(1) = 1.0
+ vertCoordMovementWeights = 0.0
+ vertCoordMovementWeights(1) = 1.0
- elseif (config_vert_grid_type.eq.'zstar') then
+ elseif (config_vert_coord_movement.eq.'uniform_stretching') then
- zstarWeight = 1.0
+ vertCoordMovementWeights = 1.0
- elseif (config_vert_grid_type.eq.'zstarWeights') then
-
- ! This is a test with other weights, just to make sure zstar functions
- ! using variable weights.
-
- zstarWeight = 0.0
- zstarWeight(1:5) = 1.0
- do k=1,10
- zstarWeight(5+k) = 1.0-k*0.1
- end do
-
endif
! Initial condition files (ocean.nc, produced by basin) include a realistic
@@ -776,13 +780,13 @@
! Compute barotropic velocity at first timestep
! This is only done upon start-up.
- if (trim(config_time_integration) == 'unsplit_explicit') then
+ if (trim(config_time_integrator) == 'unsplit_explicit') then
block % state % time_levs(1) % state % uBtr % array(:) = 0.0
block % state % time_levs(1) % state % uBcl % array(:,:) &
= block % state % time_levs(1) % state % u % array(:,:)
- elseif (trim(config_time_integration) == 'split_explicit') then
+ elseif (trim(config_time_integrator) == 'split_explicit') then
if (config_filter_btr_mode) then
do iCell=1,block % mesh % nCells
@@ -851,64 +855,6 @@
end subroutine ocn_init_split_timestep!}}}
- subroutine ocn_init_h_zstar(domain)!{{{
- ! If changing from zlevel to zstar, compute h based on zstar weights,
- ! where SSH is distributed through the layers. We only change h.
- ! We do not remap the tracer variables, so this breaks total global
- ! conservation.
-
- use mpas_grid_types
- use mpas_configure
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- type (block_type), pointer :: block
-
- integer :: i, iCell, iEdge, iVertex, k, nVertLevels
- integer, dimension(:), pointer :: maxLevelCell
-
- real (kind=RKIND) :: hSum, sumZstarWeights
- real (kind=RKIND), dimension(:), pointer :: hZLevel, zstarWeight, &
- refBottomDepth
- real (kind=RKIND), dimension(:,:), pointer :: h
-
- ! Initialize z-level grid variables from h, read in from input file.
- block => domain % blocklist
- do while (associated(block))
-
- h => block % state % time_levs(1) % state % h % array
- nVertLevels = block % mesh % nVertLevels
- hZLevel => block % mesh % hZLevel % array
- maxLevelCell => block % mesh % maxLevelCell % array
- zstarWeight => block % mesh % zstarWeight % array
- refBottomDepth => block % mesh % refBottomDepth % array
-
- do iCell=1,block % mesh % nCells
- ! Compute the total column thickness, hSum, and the sum of zstar weights.
- hSum = 0.0
- sumZstarWeights = 0.0
- do k = 1,maxLevelCell(iCell)
- hSum = hSum + h(k,iCell)
- sumZstarWeights = sumZstarWeights + zstarWeight(k)
- enddo
-
- ! h_k = h_k^{zlevel} + zeta * W_k/sum(W_k)
- ! where zeta is SSH and W_k are weights
- do k = 1,maxLevelCell(iCell)
- h(k,iCell) = hZLevel(k) &
- + (hSum - refBottomDepth(maxLevelCell(iCell))) &
- * zstarWeight(k)/sumZstarWeights
- enddo
-
- enddo
-
- block => block % next
- end do
-
- end subroutine ocn_init_h_zstar!}}}
-
subroutine ocn_compute_max_level(domain)!{{{
! Initialize maxLevel and bouncary grid variables.
@@ -959,7 +905,7 @@
! for z-grids, maxLevelCell should be in input state
! Isopycnal grid uses all vertical cells
- if (config_vert_grid_type.eq.'isopycnal') then
+ if (config_vert_coord_movement.eq.'isopycnal') then
maxLevelCell(1:nCells) = nVertLevels
endif
maxLevelCell(nCells+1) = 0
@@ -1080,7 +1026,7 @@
meshScalingDel2(:) = 1.0
meshScalingDel4(:) = 1.0
meshScaling(:) = 1.0
- if (config_h_ScaleWithMesh) then
+ if (config_hmix_ScaleWithMesh) then
do iEdge=1,mesh%nEdges
cell1 = mesh % cellsOnEdge % array(1,iEdge)
cell2 = mesh % cellsOnEdge % array(2,iEdge)
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tendency.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tendency.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -8,9 +8,7 @@
!> \version SVN:$Id:$
!> \details
!> This module contains the routines for computing
-!> various tendencies for the ocean. As well as routines
-!> for computing diagnostic variables, and other quantities
-!> such as wTop.
+!> tendency terms for the ocean primitive equations.
!
!-----------------------------------------------------------------------
@@ -25,29 +23,23 @@
use ocn_thick_hadv
use ocn_thick_vadv
- use ocn_gm
use ocn_vel_coriolis
use ocn_vel_pressure_grad
use ocn_vel_vadv
use ocn_vel_hmix
use ocn_vel_forcing
+ use ocn_vmix
use ocn_tracer_hadv
use ocn_tracer_vadv
use ocn_tracer_hmix
use ocn_restoring
- use ocn_equation_of_state
- use ocn_vmix
-
- use ocn_time_average
-
implicit none
private
save
- type (timer_node), pointer :: diagEOSTimer
type (timer_node), pointer :: thickHadvTimer, thickVadvTimer
type (timer_node), pointer :: velCorTimer, velVadvTimer, velPgradTimer, velHmixTimer, velForceTimer, velExpVmixTimer
type (timer_node), pointer :: tracerHadvTimer, tracerVadvTimer, tracerHmixTimer, tracerExpVmixTimer, tracerRestoringTimer
@@ -66,13 +58,8 @@
public :: ocn_tend_h, &
ocn_tend_u, &
- ocn_tend_scalar, &
- ocn_diagnostic_solve, &
- ocn_wtop, &
- ocn_fuperp, &
- ocn_tendency_init, &
- ocn_filter_btr_mode_u, &
- ocn_filter_btr_mode_tend_u
+ ocn_tend_tracer, &
+ ocn_tendency_init
!--------------------------------------------------------------------
!
@@ -80,11 +67,6 @@
!
!--------------------------------------------------------------------
- integer :: hadv2nd, hadv3rd, hadv4th
- integer :: ke_cell_flag, ke_vertex_flag
- real (kind=RKIND) :: coef_3rd_order, fCoef
-
-
!***********************************************************************
contains
@@ -126,6 +108,8 @@
!
tend_h = 0.0
+ if(config_disable_h_all_tend) return
+
!
! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
!
@@ -206,6 +190,8 @@
! mrp 110516 efficiency: could remove next line and have first tend_u operation not be additive
tend_u(:,:) = 0.0
+ if(config_disable_u_all_tend) return
+
!
! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
!
@@ -225,7 +211,7 @@
! velocity tendency: pressure gradient
!
call mpas_timer_start("pressure grad", .false., velPgradTimer)
- if (config_pressure_type.eq.'MontgomeryPotential') then
+ if (config_pressure_gradient_type.eq.'MontgomeryPotential') then
call ocn_vel_pressure_grad_tend(grid, MontPot, zMid, rho, tend_u, err)
else
call ocn_vel_pressure_grad_tend(grid, pressure, zMid, rho, tend_u, err)
@@ -235,7 +221,7 @@
!
! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="black">abla^2 u
! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity )
- ! strictly only valid for config_h_mom_eddy_visc2 == constant
+ ! strictly only valid for config_mom_del2 == constant
!
call mpas_timer_start("hmix", .false., velHmixTimer)
call ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend_u, err)
@@ -265,17 +251,17 @@
!***********************************************************************
!
-! routine ocn_tendSalar
+! routine ocn_tend_tracer
!
-!> \brief Computes scalar tendency
+!> \brief Computes tracer tendency
!> \author Doug Jacobsen
!> \date 23 September 2011
!> \version SVN:$Id$
!> \details
-!> This routine computes the scalar (tracer) tendency for the ocean
+!> This routine computes tracer tendencies for the ocean
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_scalar(tend, s, d, grid, dt)!{{{
+ subroutine ocn_tend_tracer(tend, s, d, grid, dt)!{{{
implicit none
type (tend_type), intent(inout) :: tend !< Input/Output: Tendency structure
@@ -291,7 +277,7 @@
integer :: err, iEdge, k
- call mpas_timer_start("ocn_tend_scalar")
+ call mpas_timer_start("ocn_tend_tracer")
uTransport => s % uTransport % array
h => s % h % array
@@ -303,6 +289,13 @@
tend_tr => tend % tracers % array
tend_h => tend % h % array
+ !
+ ! initialize tracer tendency (RHS of tracer equation) to zero.
+ !
+ tend_tr(:,:,:) = 0.0
+
+ if(config_disable_tr_all_tend) return
+
allocate(uh(grid % nVertLevels, grid % nEdges+1))
!
! QC Comment (3/15/12): need to make sure that uTransport is the right
@@ -314,11 +307,6 @@
end do
!
- ! initialize tracer tendency (RHS of tracer equation) to zero.
- !
- tend_tr(:,:,:) = 0.0
-
- !
! tracer tendency: horizontal advection term -div( h \phi u)
!
! mrp 101115 note: in order to include flux boundary conditions, we will need to
@@ -371,790 +359,14 @@
call mpas_timer_stop("restoring", tracerRestoringTimer)
10 format(2i8,10e20.10)
- call mpas_timer_stop("ocn_tend_scalar")
+ call mpas_timer_stop("ocn_tend_tracer")
deallocate(uh)
- end subroutine ocn_tend_scalar!}}}
+ end subroutine ocn_tend_tracer!}}}
!***********************************************************************
!
-! routine ocn_diagnostic_solve
-!
-!> \brief Computes diagnostic variables
-!> \author Doug Jacobsen
-!> \date 23 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the diagnostic variables for the ocean
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_diagnostic_solve(dt, s, grid)!{{{
- implicit none
-
- real (kind=RKIND), intent(in) :: dt !< Input: Time step
- type (state_type), intent(inout) :: s !< Input/Output: State information
- type (mesh_type), intent(in) :: grid !< Input: Grid information
-
-
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- integer :: boundaryMask, velMask, nCells, nEdges, nVertices, nVertLevels, vertexDegree, err
-
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
- maxLevelVertexBot
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- 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 :: &
- 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, &
- Vor_edge, Vor_vertex, Vor_cell, gradVor_n, gradVor_t, divergence, &
- rho, temperature, salinity, kev, kevc, uBolusGM, uTransport
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers, deriv_two
- real (kind=RKIND), dimension(:,:), allocatable:: div_u
- character :: c1*6
-
- h => s % h % array
- u => s % u % array
- uTransport => s % uTransport % array
- uBolusGM => s % uBolusGM % array
- v => s % v % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- 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
- rho => s % rho % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
- zMid => s % zMid % array
- ssh => s % ssh % array
- tracers => s % tracers % array
-
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- 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
- bottomDepth => grid % bottomDepth % array
- fVertex => grid % fVertex % 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
- nVertices = grid % nVertices
- nVertLevels = grid % nVertLevels
- vertexDegree = grid % vertexDegree
-
- 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
- !
- ! mrp 101115 note: in order to include flux boundary conditions, we will need to
- ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
-
- ! initialize h_edge to avoid divide by zero and NaN problems.
- h_edge = -1.0e34
- coef_3rd_order = config_coef_3rd_order
-
- do iEdge=1,nEdges*hadv2nd
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end do
-
- do iEdge=1,nEdges*hadv3rd
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- boundaryMask = abs(transfer(.not.(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0), boundaryMask))
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1) * boundaryMask
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2) * boundaryMask
-
- !-- all edges of cell 1
- do i=1, nEdgesOnCell(cell1) * boundaryMask
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, nEdgesOnCell(cell2) * boundaryMask
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- velMask = 2*(abs(transfer(u(k,iEdge) <= 0, velMask))) - 1
-
- h_edge(k,iEdge) = 0.5*(h(k,cell1) + h(k,cell2)) - (dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- + velMask * (dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
-
- end do ! do k
- end do ! do iEdge
-
- do iEdge=1,nEdges*hadv4th
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- boundaryMask = abs(transfer(.not.(boundaryCell(k,cell1) == 0 .and. boundaryCell(k,cell2) == 0), boundaryMask))
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1) * boundaryMask
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2) * boundaryMask
-
- !-- all edges of cell 1
- do i=1, nEdgesOnCell(cell1) * boundaryMask
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, nEdgesOnCell(cell2) * boundaryMask
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- end do ! do k
- end do ! do iEdge
-
- !
- ! set the velocity and height at dummy address
- ! used -1e34 so error clearly occurs if these values are used.
- !
- u(:,nEdges+1) = -1e34
- h(:,nCells+1) = -1e34
- tracers(s % index_temperature,:,nCells+1) = -1e34
- tracers(s % index_salinity,:,nCells+1) = -1e34
-
- circulation(:,:) = 0.0
- vorticity(:,:) = 0.0
- divergence(:,:) = 0.0
- ke(:,:) = 0.0
- v(:,:) = 0.0
- 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)
-
- 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
-
- 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
-
- 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)
- ! mrp 101115 note: in order to include flux boundary conditions,
- ! the following loop may need to change to maxLevelEdgeBot
- do k = 1,maxLevelEdgeTop(iEdge)
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
- end do
- end do
-
- end do
-
- !
- ! Compute kinetic energy in each vertex
- !
- kev(:,:) = 0.0; kevc(:,:) = 0.0
- 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 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
- !
- do iCell=1,nCells*ke_vertex_flag
- do k=1,nVertLevels
- ke(k,iCell) = 5.0/8.0*ke(k,iCell) + 3.0/8.0*kevc(k,iCell)
- end do
- end do
-
- !
- ! Compute ke on cell edges at velocity locations for quadratic bottom drag.
- !
- ! mrp 101025 efficiency note: we could get rid of ke_edge completely by
- ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
- end do
- end do
-
- !
- ! Compute height at vertices, pv at vertices, and average pv to edge locations
- ! ( this computes Vor_vertex at all vertices bounding real cells and distance-1 ghost cells )
- !
- do iVertex = 1,nVertices
- invAreaTri1 = 1.0 / areaTriangle(iVertex)
- do k=1,maxLevelVertexBot(iVertex)
- h_vertex = 0.0
- do i=1,vertexDegree
- h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
- end do
- h_vertex = h_vertex * invAreaTri1
-
- Vor_vertex(k,iVertex) = (fCoef*fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
- end do
- end do
-
- Vor_cell(:,:) = 0.0
- Vor_edge(:,:) = 0.0
- 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
-
- do iCell = 1, nCells
- invAreaCell1 = 1.0 / areaCell(iCell)
-
- 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
-
- do iEdge = 1,nEdges
- cell1 = cellsOnEdge(1, iEdge)
- cell2 = cellsOnEdge(2, iEdge)
- vertex1 = verticesOnedge(1, iEdge)
- vertex2 = verticesOnedge(2, iEdge)
-
- invLength = 1.0 / dcEdge(iEdge)
- ! Compute gradient of PV in normal direction
- ! ( this computes gradVor_n for all edges bounding real cells )
- do k=1,maxLevelEdgeTop(iEdge)
- gradVor_n(k,iEdge) = (Vor_cell(k,cell2) - Vor_cell(k,cell1)) * invLength
- enddo
-
- invLength = 1.0 / dvEdge(iEdge)
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradVor_t at all edges bounding real cells and distance-1 ghost cells )
- do k = 1,maxLevelEdgeBot(iEdge)
- gradVor_t(k,iEdge) = (Vor_vertex(k,vertex2) - Vor_vertex(k,vertex1)) * invLength
- enddo
-
- enddo
-
- !
- ! Modify PV edge with upstream bias.
- !
- do iEdge = 1,nEdges
- do k = 1,maxLevelEdgeBot(iEdge)
- Vor_edge(k,iEdge) = Vor_edge(k,iEdge) &
- - config_apvm_scale_factor * dt* ( u(k,iEdge) * gradVor_n(k,iEdge) &
- + v(k,iEdge) * gradVor_t(k,iEdge) )
- enddo
- enddo
-
- !
- ! equation of state
- !
- ! For an isopycnal model, density should remain constant.
- ! For zlevel, calculate in-situ density
- if (config_vert_grid_type.ne.'isopycnal') then
- call mpas_timer_start("equation of state", .false., diagEOSTimer)
- call ocn_equation_of_state_rho(s, grid, 0, 'relative', err)
- ! mrp 110324 In order to visualize rhoDisplaced, include the following
- call ocn_equation_of_state_rho(s, grid, 1, 'relative', err)
- call mpas_timer_stop("equation of state", diagEOSTimer)
- endif
-
- !
- ! Pressure
- ! This section must be after computing rho
- !
- ! dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
- if (config_pressure_type.eq.'MontgomeryPotential') then
-
- ! For Isopycnal model.
- ! Compute pressure at top of each layer, and then
- ! Montgomery Potential.
- allocate(pTop(nVertLevels))
- do iCell=1,nCells
-
- ! assume atmospheric pressure at the surface is zero for now.
- 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, bottomDepth can be relative to top or bottom)
- MontPot(1,iCell) = gravity &
- * (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)
-
- ! from delta M = p delta / rho
- MontPot(k,iCell) = MontPot(k-1,iCell) &
- + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
- end do
-
- end do
- deallocate(pTop)
-
- else
-
- do iCell=1,nCells
- ! pressure for generalized coordinates
- ! assume atmospheric pressure at the surface is zero for now.
- pressure(1,iCell) = rho(1,iCell)*gravity &
- * 0.5*h(1,iCell)
-
- do k=2,maxLevelCell(iCell)
- pressure(k,iCell) = pressure(k-1,iCell) &
- + 0.5*gravity*( rho(k-1,iCell)*h(k-1,iCell) &
- + rho(k ,iCell)*h(k ,iCell))
- end do
-
- ! 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 bottomDepth is positive
- ! and z-coordinates are negative below the surface.
- k = maxLevelCell(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) &
- + 0.5*( h(k+1,iCell) &
- + h(k ,iCell))
- end do
-
- end do
-
- endif
-
- !
- ! Sea Surface Height
- !
- do iCell=1,nCells
- ! Start at the bottom where we know the depth, and go up.
- ! 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) = - bottomDepth(iCell) + sum(h(1:maxLevelCell(iCell),iCell))
-
- end do
-
- !
- ! Apply the GM closure as a bolus velocity
- !
- if (config_h_kappa .GE. epsilon(0D0)) then
- call ocn_gm_compute_uBolus(s,grid)
- else
- ! mrp efficiency note: if uBolusGM is guaranteed to be zero, this can be removed.
- uBolusGM = 0.0
- end if
-
- end subroutine ocn_diagnostic_solve!}}}
-
-!***********************************************************************
-!
-! routine ocn_wtop
-!
-!> \brief Computes vertical velocity
-!> \author Doug Jacobsen
-!> \date 23 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the vertical velocity in the top layer for the ocean
-!
-!-----------------------------------------------------------------------
- subroutine ocn_wtop(grid,h,h_edge,u,wTop, err)!{{{
-
- !-----------------------------------------------------------------
- !
- ! 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, invAreaCell
-
- integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
-
-
- real (kind=RKIND), dimension(:), pointer :: &
- 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, edgeSignOnCell
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
- maxLevelVertexBot, maxLevelVertexTop
-
- 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
- zstarWeight => grid % zstarWeight % array
-
- nCells = grid % nCells
- nEdges = grid % nEdges
- nVertLevels = grid % 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.
- !
-
- do iCell=1,nCells
- div_hu(:) = 0.0_RKIND
- div_hu_btr = 0.0_RKIND
- hSum = 0.0_RKIND
- invAreaCell = 1.0_RKIND / 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)
- 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 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
-
- 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_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
-
- deallocate(div_hu, h_tend_col)
-
- end subroutine ocn_wtop!}}}
-
-!***********************************************************************
-!
-! routine ocn_fuperp
-!
-!> \brief Computes f u_perp
-!> \author Doug Jacobsen
-!> \date 23 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes f u_perp for the ocean
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_fuperp(s, grid)!{{{
- implicit none
-
- type (state_type), intent(inout) :: s !< Input/Output: State information
- type (mesh_type), intent(in) :: grid !< Input: Grid information
-
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-! Some of these variables can be removed, but at a later time.
- integer :: iEdge, cell1, cell2, eoe, i, j, k
-
- integer :: nEdgesSolve
- real (kind=RKIND), dimension(:), pointer :: fEdge
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, u, uBcl
- type (dm_info) :: dminfo
-
- integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnEdge
- integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnEdge
-
- call mpas_timer_start("ocn_fuperp")
-
- u => s % u % array
- uBcl => s % uBcl % array
- weightsOnEdge => grid % weightsOnEdge % array
- fEdge => grid % fEdge % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- cellsOnEdge => grid % cellsOnEdge % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
-
- fEdge => grid % fEdge % array
-
- nEdgesSolve = grid % nEdgesSolve
-
- !
- ! Put f*uBcl^{perp} in u as a work variable
- !
- do iEdge=1,nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
-
- u(k,iEdge) = 0.0
- do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- u(k,iEdge) = u(k,iEdge) + weightsOnEdge(j,iEdge) * uBcl(k,eoe) * fEdge(eoe)
- end do
- end do
- end do
-
- call mpas_timer_stop("ocn_fuperp")
-
- end subroutine ocn_fuperp!}}}
-
-!***********************************************************************
-!
-! routine ocn_filter_btr_mode_u
-!
-!> \brief filters barotropic mode out of the velocity variable.
-!> \author Mark Petersen
-!> \date 23 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine filters barotropic mode out of the velocity variable.
-!
-!-----------------------------------------------------------------------
- subroutine ocn_filter_btr_mode_u(s, grid)!{{{
- implicit none
-
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
-
- integer :: iEdge, k, nEdges
- real (kind=RKIND) :: vertSum, uhSum, hSum
- real (kind=RKIND), dimension(:,:), pointer :: h_edge, u
- integer, dimension(:), pointer :: maxLevelEdgeTop
-
- call mpas_timer_start("ocn_filter_btr_mode_u")
-
- u => s % u % array
- h_edge => s % h_edge % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- nEdges = grid % nEdges
-
- do iEdge=1,nEdges
-
- ! hSum is initialized outside the loop because on land boundaries
- ! maxLevelEdgeTop=0, but I want to initialize hSum with a
- ! nonzero value to avoid a NaN.
- uhSum = h_edge(1,iEdge) * u(1,iEdge)
- hSum = h_edge(1,iEdge)
-
- do k=2,maxLevelEdgeTop(iEdge)
- uhSum = uhSum + h_edge(k,iEdge) * u(k,iEdge)
- hSum = hSum + h_edge(k,iEdge)
- enddo
-
- vertSum = uhSum/hSum
- do k=1,maxLevelEdgeTop(iEdge)
- u(k,iEdge) = u(k,iEdge) - vertSum
- enddo
- enddo ! iEdge
-
- call mpas_timer_stop("ocn_filter_btr_mode_u")
-
- end subroutine ocn_filter_btr_mode_u!}}}
-
-!***********************************************************************
-!
-! routine ocn_filter_btr_mode_tend_u
-!
-!> \brief ocn_filters barotropic mode out of the u tendency
-!> \author Mark Petersen
-!> \date 23 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine filters barotropic mode out of the u tendency.
-!
-!-----------------------------------------------------------------------
- subroutine ocn_filter_btr_mode_tend_u(tend, s, grid)!{{{
- implicit none
-
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (mesh_type), intent(in) :: grid
-
- integer :: iEdge, k, nEdges
- real (kind=RKIND) :: vertSum, uhSum, hSum
- real (kind=RKIND), dimension(:,:), pointer :: h_edge, tend_u
-
- integer, dimension(:), pointer :: maxLevelEdgeTop
-
- call mpas_timer_start("ocn_filter_btr_mode_tend_u")
-
- tend_u => tend % u % array
- h_edge => s % h_edge % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- nEdges = grid % nEdges
-
- do iEdge=1,nEdges
-
- ! hSum is initialized outside the loop because on land boundaries
- ! maxLevelEdgeTop=0, but I want to initialize hSum with a
- ! nonzero value to avoid a NaN.
- uhSum = h_edge(1,iEdge) * tend_u(1,iEdge)
- hSum = h_edge(1,iEdge)
-
- do k=2,maxLevelEdgeTop(iEdge)
- uhSum = uhSum + h_edge(k,iEdge) * tend_u(k,iEdge)
- hSum = hSum + h_edge(k,iEdge)
- enddo
-
- vertSum = uhSum/hSum
- do k=1,maxLevelEdgeTop(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) - vertSum
- enddo
- enddo ! iEdge
-
- call mpas_timer_stop("ocn_filter_btr_mode_tend_u")
-
- end subroutine ocn_filter_btr_mode_tend_u!}}}
-
-!***********************************************************************
-!
! routine ocn_tendency_init
!
!> \brief Initializes flags used within tendency routines.
@@ -1171,49 +383,6 @@
err = 0
- coef_3rd_order = 0.
-
- if (config_thickness_adv_order == 2) then
- hadv2nd = 1
- hadv3rd = 0
- hadv4th = 0
- else if (config_thickness_adv_order == 3) then
- hadv2nd = 0
- hadv3rd = 1
- hadv4th = 0
-
- if(config_monotonic) then
- coef_3rd_order = 0.25
- else
- coef_3rd_order = 1.0
- endif
- else if (config_thickness_adv_order == 4) then
- hadv2nd = 0
- hadv3rd = 0
- hadv4th = 1
- end if
-
-
- if(config_include_KE_vertex) then
- ke_vertex_flag = 1
- ke_cell_flag = 0
- else
- ke_vertex_flag = 0
- ke_cell_flag = 1
- endif
-
- if (trim(config_time_integration) == 'RK4') then
- ! for RK4, PV is really PV = (eta+f)/h
- fCoef = 1
- elseif (trim(config_time_integration) == 'split_explicit' &
- .or.trim(config_time_integration) == 'unsplit_explicit') then
- ! for split explicit, PV is eta/h because f is added separately to the momentum forcing.
- ! mrp temp, new should be:
- fCoef = 0
- ! old, for testing:
- ! fCoef = 1
- end if
-
end subroutine ocn_tendency_init!}}}
!***********************************************************************
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_thick_hadv.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_thick_hadv.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_thick_hadv.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -41,8 +41,9 @@
! Private module variables
!
!--------------------------------------------------------------------
+
+ logical :: thickHadvOn
-
!***********************************************************************
contains
@@ -120,6 +121,8 @@
err = 0
+ if(.not.thickHadvOn) return
+
nEdges = grid % nEdges
nCells = grid % nCells
nVertLevels = grid % nVertLevels
@@ -175,6 +178,10 @@
integer, intent(out) :: err !< Output: error flag
+ thickHadvOn = .true.
+
+ if(config_disable_h_hadv) thickHadvOn = .false.
+
err = 0
!--------------------------------------------------------------------
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_thick_vadv.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_thick_vadv.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_thick_vadv.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -42,6 +42,7 @@
!
!--------------------------------------------------------------------
+ logical :: thickVadvOn
!***********************************************************************
@@ -111,6 +112,8 @@
err = 0
+ if(.not.thickVadvOn) return
+
maxLevelCell => grid % maxLevelCell % array
nCells = grid % nCells
@@ -151,6 +154,10 @@
!-----------------------------------------------------------------
integer, intent(out) :: err !< Output: error flag
+
+ thickVadvOn = .true.
+
+ if(config_disable_h_vadv) thickVadvOn = .false.
err = 0
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_average.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_average.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_average.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -15,7 +15,7 @@
real (kind=RKIND), dimension(:), pointer :: acc_ssh, acc_sshVar
real (kind=RKIND), dimension(:,:), pointer :: acc_uReconstructZonal, acc_uReconstructMeridional, acc_uReconstructZonalVar, acc_uReconstructMeridionalVar
- real (kind=RKIND), dimension(:,:), pointer :: acc_u, acc_uVar
+ real (kind=RKIND), dimension(:,:), pointer :: acc_u, acc_uVar, acc_vertVelocityTop
nAccumulate => state % nAccumulate % scalar
@@ -27,6 +27,7 @@
acc_uReconstructMeridionalVar => state % acc_uReconstructMeridionalVar % array
acc_u => state % acc_u % array
acc_uVar => state % acc_uVar % array
+ acc_vertVelocityTop => state % acc_vertVelocityTop % array
nAccumulate = 0
@@ -38,6 +39,7 @@
acc_uReconstructMeridionalVar = 0.0
acc_u = 0.0
acc_uVar = 0.0
+ acc_vertVelocityTop = 0.0
end subroutine ocn_time_average_init!}}}
@@ -48,13 +50,13 @@
real (kind=RKIND), pointer :: nAccumulate, old_nAccumulate
real (kind=RKIND), dimension(:), pointer :: ssh
- real (kind=RKIND), dimension(:,:), pointer :: uReconstructZonal, uReconstructMeridional, u
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructZonal, uReconstructMeridional, u, vertVelocityTop
- real (kind=RKIND), dimension(:,:), pointer :: acc_u, acc_uVar
+ real (kind=RKIND), dimension(:,:), pointer :: acc_u, acc_uVar, acc_vertVelocityTop
real (kind=RKIND), dimension(:,:), pointer :: acc_uReconstructZonal, acc_uReconstructMeridional, acc_uReconstructZonalVar, acc_uReconstructMeridionalVar
real (kind=RKIND), dimension(:), pointer :: acc_ssh, acc_sshVar
- real (kind=RKIND), dimension(:,:), pointer :: old_acc_u, old_acc_uVar
+ real (kind=RKIND), dimension(:,:), pointer :: old_acc_u, old_acc_uVar, old_acc_vertVelocityTop
real (kind=RKIND), dimension(:,:), pointer :: old_acc_uReconstructZonal, old_acc_uReconstructMeridional, old_acc_uReconstructZonalVar, old_acc_uReconstructMeridionalVar
real (kind=RKIND), dimension(:), pointer :: old_acc_ssh, old_acc_sshVar
@@ -65,6 +67,7 @@
uReconstructZonal => state % uReconstructZonal % array
uReconstructMeridional => state % uReconstructMeridional % array
u => state % u % array
+ vertVelocityTop => state % vertVelocityTop % array
acc_ssh => state % acc_ssh % array
acc_sshVar => state % acc_sshVar % array
@@ -74,6 +77,7 @@
acc_uReconstructMeridionalVar => state % acc_uReconstructMeridionalVar % array
acc_u => state % acc_u % array
acc_uVar => state % acc_uVar % array
+ acc_vertVelocityTop => state % acc_vertVelocityTop % array
old_acc_ssh => old_state % acc_ssh % array
old_acc_sshVar => old_state % acc_sshVar % array
@@ -83,6 +87,7 @@
old_acc_uReconstructMeridionalVar => old_state % acc_uReconstructMeridionalVar % array
old_acc_u => old_state % acc_u % array
old_acc_uVar => old_state % acc_uVar % array
+ old_acc_vertVelocityTop => old_state % acc_vertVelocityTop % array
acc_ssh = old_acc_ssh + ssh
acc_sshVar = old_acc_sshVar + ssh**2
@@ -92,6 +97,7 @@
acc_uReconstructMeridionalVar = old_acc_uReconstructMeridionalVar + uReconstructMeridional**2
acc_u = old_acc_u + u
acc_uVar = old_acc_uVar + u**2
+ acc_vertVelocityTop = old_acc_vertVelocityTop + vertVelocityTop
nAccumulate = old_nAccumulate + 1
end subroutine ocn_time_average_accumulate!}}}
@@ -103,7 +109,7 @@
real (kind=RKIND), dimension(:), pointer :: acc_ssh, acc_sshVar
real (kind=RKIND), dimension(:,:), pointer :: acc_uReconstructZonal, acc_uReconstructMeridional, acc_uReconstructZonalVar, acc_uReconstructMeridionalVar
- real (kind=RKIND), dimension(:,:), pointer :: acc_u, acc_uVar
+ real (kind=RKIND), dimension(:,:), pointer :: acc_u, acc_uVar, acc_vertVelocityTop
nAccumulate => state % nAccumulate % scalar
@@ -115,6 +121,7 @@
acc_uReconstructMeridionalVar => state % acc_uReconstructMeridionalVar % array
acc_u => state % acc_u % array
acc_uVar => state % acc_uVar % array
+ acc_vertVelocityTop => state % acc_vertVelocityTop % array
if(nAccumulate > 0) then
acc_ssh = acc_ssh / nAccumulate
@@ -125,6 +132,7 @@
acc_uReconstructMeridionalVar = acc_uReconstructMeridionalVar / nAccumulate
acc_u = acc_u / nAccumulate
acc_uVar = acc_uVar / nAccumulate
+ acc_vertVelocityTop = acc_vertVelocityTop / nAccumulate
end if
end subroutine ocn_time_average_normalize!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -121,14 +121,14 @@
rk4On = .false.
splitOn = .false.
- if (trim(config_time_integration) == 'RK4') then
+ if (trim(config_time_integrator) == 'RK4') then
rk4On = .true.
- elseif (trim(config_time_integration) == 'split_explicit' &
- .or.trim(config_time_integration) == 'unsplit_explicit') then
+ elseif (trim(config_time_integrator) == 'split_explicit' &
+ .or.trim(config_time_integrator) == 'unsplit_explicit') then
splitOn = .true.
else
err = 1
- write(*,*) 'Incorrect choice for config_time_integration:', trim(config_time_integration)
+ write(*,*) 'Incorrect choice for config_time_integrator:', trim(config_time_integrator)
write(*,*) ' choices are: RK4, split_explicit, unsplit_explicit'
endif
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -22,6 +22,7 @@
use mpas_timer
use ocn_tendency
+ use ocn_diagnostics
use ocn_equation_of_state
use ocn_vmix
@@ -137,7 +138,7 @@
call mpas_timer_start("RK4-diagnostic halo update")
call mpas_dmpar_exch_halo_field(domain % blocklist % provis % Vor_edge)
- if (config_h_mom_eddy_visc4 > 0.0) then
+ if (config_mom_del4 > 0.0) then
call mpas_dmpar_exch_halo_field(domain % blocklist % provis % divergence)
call mpas_dmpar_exch_halo_field(domain % blocklist % provis % vorticity)
end if
@@ -162,11 +163,11 @@
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
+ if (config_filter_btr_mode) then
call ocn_filter_btr_mode_tend_u(block % tend, block % provis, block % mesh)
endif
- call ocn_tend_scalar(block % tend, block % provis, block % diagnostics, block % mesh, dt)
+ call ocn_tend_tracer(block % tend, block % provis, block % diagnostics, block % mesh, dt)
block => block % next
end do
call mpas_timer_stop("RK4-tendency computations")
@@ -254,7 +255,7 @@
call mpas_timer_stop("RK4-main loop")
!
- ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
+ ! A little clean up at the end: rescale tracer fields and compute diagnostics for new state
!
call mpas_timer_start("RK4-cleaup phase")
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -24,6 +24,7 @@
use mpas_timer
use ocn_tendency
+ use ocn_diagnostics
use ocn_equation_of_state
use ocn_vmix
@@ -167,7 +168,7 @@
call mpas_timer_start("se halo diag", .false., timer_halo_diagnostic)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % Vor_edge)
- if (config_h_mom_eddy_visc4 > 0.0) then
+ if (config_mom_del4 > 0.0) then
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % divergence)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % vorticity)
end if
@@ -209,9 +210,9 @@
do j=1,n_bcl_iter(split_explicit_step)
! Use this G coefficient to avoid an if statement within the iEdge loop.
- if (trim(config_time_integration) == 'unsplit_explicit') then
+ if (trim(config_time_integrator) == 'unsplit_explicit') then
split = 0
- elseif (trim(config_time_integration) == 'split_explicit') then
+ elseif (trim(config_time_integrator) == 'split_explicit') then
split = 1
endif
@@ -294,7 +295,7 @@
oldBtrSubcycleTime = 1
newBtrSubcycleTime = 2
- if (trim(config_time_integration) == 'unsplit_explicit') then
+ if (trim(config_time_integrator) == 'unsplit_explicit') then
block => domain % blocklist
do while (associated(block))
@@ -321,7 +322,7 @@
block => block % next
end do ! block
- elseif (trim(config_time_integration) == 'split_explicit') then
+ elseif (trim(config_time_integrator) == 'split_explicit') then
! Initialize variables for barotropic subcycling
block => domain % blocklist
@@ -802,7 +803,7 @@
block => domain % blocklist
do while (associated(block))
- call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh, dt)
+ call ocn_tend_tracer(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh, dt)
block => block % next
end do
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -32,6 +32,7 @@
mpas_ocn_tracer_advection_coefficients, &
mpas_ocn_tracer_advection_tend
+ logical :: tracerAdvOn
logical :: monotonicOn
contains
@@ -251,6 +252,8 @@
type (mesh_type), intent(in) :: grid !< Input: grid information
real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !< Input: Thickness tendency information
+ if(.not. tracerAdvOn) return
+
if(monotonicOn) then
call mpas_ocn_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)
else
@@ -279,6 +282,10 @@
err = 0
+ tracerAdvOn = .true.
+
+ if(config_disable_tr_adv) tracerAdvOn = .false.
+
call mpas_ocn_tracer_advection_std_init(err_tmp)
call mpas_ocn_tracer_advection_mono_init(err_tmp)
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -147,7 +147,7 @@
upwind_tendency(k, iCell) = 0.0
!tracer_new is supposed to be the "new" tracer state. This allows bounds checks.
- if (config_check_monotonicity) then
+ if (config_check_tracer_monotonicity) then
tracer_new(k,iCell) = 0.0
end if
end do ! k loop
@@ -320,7 +320,7 @@
do k = 1, maxLevelEdgeTop(iEdge)
tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
- if(config_check_monotonicity) then
+ if(config_check_tracer_monotonicity) then
tracer_new(k, iCell) = tracer_new(k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
end if
end do
@@ -332,7 +332,7 @@
do k = 1,maxLevelCell(iCell)
tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
- if (config_check_monotonicity) then
+ if (config_check_tracer_monotonicity) then
!tracer_new holds a tendency for now. Only for a check on monotonicity
tracer_new(k, iCell) = tracer_new(k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
@@ -342,7 +342,7 @@
end do ! k loop
end do ! iCell loop
- if (config_check_monotonicity) then
+ if (config_check_tracer_monotonicity) then
!build min and max bounds on old and new tracer for check on monotonicity.
do iCell = 1, nCellsSolve
do k = 1, maxLevelCell(iCell)
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection_std.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -51,7 +51,7 @@
! Input: s - current model state
! grid - grid metadata
!
- ! Output: tend - computed scalar tendencies
+ ! Output: tend - computed tracer tendencies
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -47,6 +47,7 @@
!
!--------------------------------------------------------------------
+ logical :: tracerHmixOn
type (timer_node), pointer :: del2Timer, del4Timer
@@ -122,6 +123,8 @@
!
!-----------------------------------------------------------------
+ if(.not.tracerHmixOn) return
+
call mpas_timer_start("del2", .false., del2Timer)
call ocn_tracer_hmix_del2_tend(grid, h_edge, tracers, tend, err1)
call mpas_timer_stop("del2", del2Timer)
@@ -165,6 +168,9 @@
integer :: err1, err2
+ tracerHmixOn = .true.
+
+ if(config_disable_tr_hmix) tracerHmixOn = .false.
call ocn_tracer_hmix_del2_init(err1)
call ocn_tracer_hmix_del4_init(err2)
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -211,11 +211,13 @@
del2on = .false.
- if ( config_h_tracer_eddy_diff2 > 0.0 ) then
+ if ( config_tracer_del2 > 0.0 ) then
del2On = .true.
- eddyDiff2 = config_h_tracer_eddy_diff2
+ eddyDiff2 = config_tracer_del2
endif
+ if(.not.config_use_tracer_del2) del2on = .false.
+
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del2_init!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -233,11 +233,13 @@
err = 0
del4on = .false.
- if ( config_h_tracer_eddy_diff4 > 0.0 ) then
+ if ( config_tracer_del4 > 0.0 ) then
del4On = .true.
- eddyDiff4 = config_h_tracer_eddy_diff4
+ eddyDiff4 = config_tracer_del4
endif
+ if(.not.config_use_tracer_del4) del4on = .false.
+
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del4_init!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_vadv.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_vadv.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_tracer_vadv.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -172,7 +172,7 @@
err = 0
vadvOn = .false.
- if (config_vert_grid_type.ne.'isopycnal') then
+ if (config_vert_coord_movement.ne.'isopycnal') then
vadvOn = .true.
call ocn_tracer_vadv_stencil_init(err1)
call ocn_tracer_vadv_spline_init(err2)
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_coriolis.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_coriolis.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_coriolis.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -43,6 +43,7 @@
!
!--------------------------------------------------------------------
+ logical :: coriolisOn
!***********************************************************************
@@ -116,6 +117,8 @@
err = 0
+ if(.not.coriolisOn) return
+
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
cellsOnEdge => grid % cellsOnEdge % array
@@ -181,6 +184,10 @@
err = 0
+ coriolisOn = .true.
+
+ if(config_disable_u_coriolis) coriolisOn = .false.
+
!--------------------------------------------------------------------
end subroutine ocn_vel_coriolis_init!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing_windstress.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing_windstress.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_forcing_windstress.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -168,6 +168,8 @@
windStressOn = .true.
+ if(config_disable_u_windstress) windStressOn = .false.
+
err = 0
!--------------------------------------------------------------------
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -48,6 +48,7 @@
!
!--------------------------------------------------------------------
+ logical :: hmixOn
type (timer_node), pointer :: del2Timer, leithTimer, del4Timer
@@ -126,6 +127,8 @@
!
!-----------------------------------------------------------------
+ if(.not.hmixOn) return
+
viscosity = 0.0
call mpas_timer_start("del2", .false., del2Timer)
@@ -176,12 +179,16 @@
integer :: err1, err2, err3
+ hmixOn = .true.
+
call ocn_vel_hmix_del2_init(err1)
call ocn_vel_hmix_leith_init(err2)
call ocn_vel_hmix_del4_init(err3)
err = ior(ior(err1, err2),err3)
+ if(config_disable_u_hmix) hmixOn = .false.
+
!--------------------------------------------------------------------
end subroutine ocn_vel_hmix_init!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_del2.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -201,19 +201,20 @@
hmixDel2On = .false.
- if ( config_h_mom_eddy_visc2 > 0.0 ) then
+ if ( config_mom_del2 > 0.0 ) then
hmixDel2On = .true.
- eddyVisc2 = config_h_mom_eddy_visc2
+ eddyVisc2 = config_mom_del2
-
if (config_visc_vorticity_term) then
- viscVortCoef = config_visc_vorticity_visc2_scale
+ viscVortCoef = config_vorticity_del2_scale
else
viscVortCoef = 0.0
endif
endif
+ if(.not.config_use_mom_del2) hmixDel2On = .false.
+
!--------------------------------------------------------------------
end subroutine ocn_vel_hmix_del2_init!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -261,17 +261,19 @@
hmixDel4On = .false.
- if ( config_h_mom_eddy_visc4 > 0.0 ) then
+ if ( config_mom_del4 > 0.0 ) then
hmixDel4On = .true.
- eddyVisc4 = config_h_mom_eddy_visc4
+ eddyVisc4 = config_mom_del4
if (config_visc_vorticity_term) then
- viscVortCoef = config_visc_vorticity_visc4_scale
+ viscVortCoef = config_vorticity_del4_scale
else
viscVortCoef = 0.0
endif
endif
+ if(.not.config_use_mom_del4) hmixDel4On = .false.
+
!--------------------------------------------------------------------
end subroutine ocn_vel_hmix_del4_init!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_leith.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -216,7 +216,7 @@
hmixLeithOn = .true.
if (config_visc_vorticity_term) then
- viscVortCoef = config_visc_vorticity_visc2_scale
+ viscVortCoef = config_vorticity_del2_scale
else
viscVortCoef = 0.0
endif
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_pressure_grad.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -44,6 +44,7 @@
!
!--------------------------------------------------------------------
+ logical :: pgradOn
real (kind=RKIND) :: rho0Inv, grho0Inv
@@ -113,6 +114,8 @@
err = 0
+ if(.not.pgradOn) return
+
nEdgesSolve = grid % nEdgesSolve
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
cellsOnEdge => grid % cellsOnEdge % array
@@ -186,7 +189,9 @@
err = 0
- if (config_pressure_type.eq.'MontgomeryPotential') then
+ pgradOn = .true.
+
+ if (config_pressure_gradient_type.eq.'MontgomeryPotential') then
rho0Inv = 1.0
grho0Inv = 0.0
else
@@ -194,8 +199,10 @@
grho0Inv = gravity/config_rho0
end if
+ if(config_disable_u_pgrad) pgradOn = .false.
+
!--------------------------------------------------------------------
end subroutine ocn_vel_pressure_grad_init!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_vadv.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_vadv.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vel_vadv.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -183,10 +183,12 @@
err = 0
velVadvOn = .false.
- if (config_vert_grid_type.ne.'isopycnal') then
+ if (config_vert_coord_movement.ne.'isopycnal') then
velVadvOn = .true.
end if
+ if(config_disable_u_vadv) velVadvOn = .false.
+
!--------------------------------------------------------------------
end subroutine ocn_vel_vadv_init!}}}
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -15,6 +15,7 @@
module ocn_vmix
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
use mpas_timer
@@ -56,6 +57,7 @@
!
!--------------------------------------------------------------------
+ logical :: velVmixOn, tracerVmixOn
logical :: explicitOn, implicitOn
!***********************************************************************
@@ -123,6 +125,9 @@
err = 0
+ d % vertViscTopOfEdge % array = 0.0_RKIND
+ d % vertDiffTopOfCell % array = 0.0_RKIND
+
call ocn_vmix_coefs_const_build(grid, s, d, err1)
call ocn_vmix_coefs_tanh_build(grid, s, d, err2)
call ocn_vmix_coefs_rich_build(grid, s, d, err3)
@@ -198,6 +203,7 @@
err = 0
+ if(.not.velVmixOn) return
if(implicitOn) return
nEdgessolve = grid % nEdgesSolve
@@ -299,6 +305,7 @@
err = 0
+ if(.not.velVmixOn) return
if(explicitOn) return
nEdges = grid % nEdges
@@ -425,6 +432,7 @@
err = 0
+ if(.not.tracerVmixOn) return
if(implicitOn) return
nCellsSolve = grid % nCellsSolve
@@ -530,6 +538,7 @@
err = 0
+ if(.not.tracerVmixOn) return
if(explicitOn) return
nCells = grid % nCells
@@ -662,6 +671,9 @@
err = 0
+ velVmixOn = .true.
+ tracerVmixOn = .true.
+
explicitOn = .true.
implicitOn = .false.
@@ -670,6 +682,9 @@
implicitOn = .true.
end if
+ if(config_disable_u_vmix) velVmixOn = .false.
+ if(config_disable_tr_vmix) tracerVmixOn = .false.
+
call ocn_vmix_coefs_const_init(err1)
call ocn_vmix_coefs_tanh_init(err2)
call ocn_vmix_coefs_rich_init(err3)
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_const.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_const.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_const.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -121,7 +121,6 @@
!-----------------------------------------------------------------
err = 0
- if((.not.constViscOn) .and. (.not.constDiffOn)) return
vertViscTopOfEdge => d % vertViscTopOfEdge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
@@ -185,7 +184,7 @@
if(.not.constViscOn) return
- vertViscTopOfEdge = constVisc
+ vertViscTopOfEdge = vertViscTopOfEdge + constVisc
!--------------------------------------------------------------------
@@ -241,7 +240,7 @@
if(.not.constDiffOn) return
- vertDiffTopOfCell = constDiff
+ vertDiffTopOfCell = vertDiffTopOfCell + constDiff
!--------------------------------------------------------------------
@@ -279,18 +278,20 @@
err = 0
- constViscOn = .false.
- constDiffOn = .false.
+ constViscOn = config_use_const_visc
+ constDiffOn = config_use_const_diff
+ constVisc = config_vert_visc
+ constDiff = config_vert_diff
- if (config_vert_visc_type.eq.'const') then
- constViscOn = .true.
- constVisc = config_vert_visc
- endif
+! if (config_vert_visc_type.eq.'const') then
+! constViscOn = .true.
+! constVisc = config_vert_visc
+! endif
- if (config_vert_diff_type.eq.'const') then
- constDiffOn = .true.
- constDiff = config_vert_diff
- endif
+! if (config_vert_diff_type.eq.'const') then
+! constDiffOn = .true.
+! constDiff = config_vert_diff
+! endif
!--------------------------------------------------------------------
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -124,7 +124,6 @@
!-----------------------------------------------------------------
err = 0
- if((.not.richViscOn) .and. (.not.richDiffOn)) return
indexT = s%index_temperature
indexS = s%index_salinity
@@ -222,13 +221,12 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- vertViscTopOfEdge = 0.0
do iEdge = 1,nEdges
do k = 2,maxLevelEdgeTop(iEdge)
! mrp 110324 efficiency note: this if is inside iEdge and k loops.
! Perhaps there is a more efficient way to do this.
if (RiTopOfEdge(k,iEdge)>0.0) then
- vertViscTopOfEdge(k,iEdge) = config_bkrd_vert_visc &
+ vertViscTopOfEdge(k,iEdge) = vertViscTopOfEdge(k, iEdge) + config_bkrd_vert_visc &
+ config_rich_mix / (1.0 + 5.0*RiTopOfEdge(k,iEdge))**2
! maltrud do limiting of coefficient--should not be necessary
! also probably better logic could be found
@@ -248,8 +246,7 @@
else
! for Ri<0 and explicit vertical mix,
! use maximum diffusion allowed by CFL criterion
- vertViscTopOfEdge(k,iEdge) = &
- ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
+ vertViscTopOfEdge(k,iEdge) = vertViscTopOfEdge(k,iEdge) + ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
end if
end if
end do
@@ -326,14 +323,13 @@
maxLevelCell => grid % maxLevelCell % array
- vertDiffTopOfCell = 0.0
coef = -gravity/config_rho0/2.0
do iCell = 1,nCells
do k = 2,maxLevelCell(iCell)
! mrp 110324 efficiency note: this if is inside iCell and k loops.
! Perhaps there is a more efficient way to do this.
if (RiTopOfCell(k,iCell)>0.0) then
- vertDiffTopOfCell(k,iCell) = config_bkrd_vert_diff &
+ vertDiffTopOfCell(k,iCell) = vertDiffTopOfCell(k, iCell) + config_bkrd_vert_diff &
+ (config_bkrd_vert_visc &
+ config_rich_mix / (1.0 + 5.0*RiTopOfCell(k,iCell))**2) &
/ (1.0 + 5.0*RiTopOfCell(k,iCell))
@@ -355,8 +351,7 @@
else
! for Ri<0 and explicit vertical mix,
! use maximum diffusion allowed by CFL criterion
- vertDiffTopOfCell(k,iCell) = &
- ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
+ vertDiffTopOfCell(k,iCell) = vertDiffTopOfCell(k, iCell) + ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
end if
end if
end do
@@ -573,16 +568,16 @@
err = 0
- richViscOn = .false.
- richDiffOn = .false.
+ richViscOn = config_use_rich_visc
+ richDiffOn = config_use_rich_diff
- if (config_vert_visc_type.eq.'rich') then
- richViscOn = .true.
- endif
+! if (config_vert_visc_type.eq.'rich') then
+! richViscOn = .true.
+! endif
- if (config_vert_diff_type.eq.'rich') then
- richDiffOn = .true.
- endif
+! if (config_vert_diff_type.eq.'rich') then
+! richDiffOn = .true.
+! endif
!--------------------------------------------------------------------
Modified: branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F
===================================================================
--- branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -115,7 +115,6 @@
!-----------------------------------------------------------------
err = 0
- if((.not.tanhViscOn) .and. (.not.tanhDiffOn)) return
vertViscTopOfEdge => d % vertViscTopOfEdge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
@@ -189,9 +188,8 @@
! 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 &
+ vertViscTopOfEdge(k,:) = vertViscTopOfEdge(k,:)-(config_max_visc_tanh-config_min_visc_tanh)/2.0 &
*tanh((refBottomDepth(k-1)+config_ZMid_tanh) &
/config_zWidth_tanh) &
+ (config_max_visc_tanh+config_min_visc_tanh)/2
@@ -262,9 +260,8 @@
! 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 &
+ vertDiffTopOfCell(k,:) = vertDiffTopOfCell(k,:)-(config_max_diff_tanh-config_min_diff_tanh)/2.0 &
*tanh((refBottomDepth(k-1)+config_ZMid_tanh) &
/config_zWidth_tanh) &
+ (config_max_diff_tanh+config_min_diff_tanh)/2
@@ -305,16 +302,16 @@
err = 0
- tanhViscOn = .false.
- tanhDiffOn = .false.
+ tanhViscOn = config_use_tanh_visc
+ tanhDiffOn = config_use_tanh_diff
- if (config_vert_visc_type.eq.'tanh') then
- tanhViscOn = .true.
- endif
+! if (config_vert_visc_type.eq.'tanh') then
+! tanhViscOn = .true.
+! endif
- if (config_vert_diff_type.eq.'tanh') then
- tanhDiffOn = .true.
- endif
+! if (config_vert_diff_type.eq.'tanh') then
+! tanhDiffOn = .true.
+! endif
!--------------------------------------------------------------------
Modified: branches/mpas_cdg_advection/src/framework/mpas_timer.F
===================================================================
--- branches/mpas_cdg_advection/src/framework/mpas_timer.F 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/framework/mpas_timer.F 2013-02-13 21:29:30 UTC (rev 2473)
@@ -269,7 +269,7 @@
timer_ptr%avg_time = 0.0d0
percent = 0.0d0
else
- timer_ptr%avg_time = timer_ptr%avg_time/timer_ptr%calls
+ timer_ptr%avg_time = timer_ptr%total_time/timer_ptr%calls
percent = timer_ptr%total_time/total_ptr%total_time
endif
Modified: branches/mpas_cdg_advection/src/registry/gen_inc.c
===================================================================
--- branches/mpas_cdg_advection/src/registry/gen_inc.c 2013-02-13 18:29:41 UTC (rev 2472)
+++ branches/mpas_cdg_advection/src/registry/gen_inc.c 2013-02-13 21:29:30 UTC (rev 2473)
@@ -718,35 +718,37 @@
var_list_ptr3 = var_list_ptr3->next;
}
- fortprintf(fd, " allocate(%s %% %s %% array(%i, ", group_ptr->name, var_ptr2->super_array, i);
- dimlist_ptr = var_ptr2->dimlist;
- if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
- !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
- !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
- if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
- else fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_file);
- else
- if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s", dimlist_ptr->dim->name_in_file);
- else fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
- dimlist_ptr = dimlist_ptr->next;
- while (dimlist_ptr) {
+ if(var_ptr2->persistence == PERSISTENT){
+ fortprintf(fd, " allocate(%s %% %s %% array(%i, ", group_ptr->name, var_ptr2->super_array, i);
+ dimlist_ptr = var_ptr2->dimlist;
if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
!strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
- if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
- else fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_file);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_file);
else
- if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_file);
- else fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
- }
- fortprintf(fd, "))</font>
<font color="red">");
- if (var_ptr->vtype == INTEGER)
- fortprintf(fd, " %s %% %s %% array = 0</font>
<font color="red">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
- else if (var_ptr->vtype == REAL)
- fortprintf(fd, " %s %% %s %% array = 0.0</font>
<font color="red">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
- else if (var_ptr->vtype == CHARACTER)
- fortprintf(fd, " %s %% %s %% array = \'\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
+ while (dimlist_ptr) {
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_file);
+ else
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ dimlist_ptr = dimlist_ptr->next;
+ }
+ fortprintf(fd, "))</font>
<font color="blue">");
+ if (var_ptr->vtype == INTEGER)
+ fortprintf(fd, " %s %% %s %% array = 0</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
+ else if (var_ptr->vtype == REAL)
+ fortprintf(fd, " %s %% %s %% array = 0.0</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
+ else if (var_ptr->vtype == CHARACTER)
+ fortprintf(fd, " %s %% %s %% array = \'\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
+ }
fortprintf(fd, " %s %% %s %% dimSizes(1) = %i</font>
<font color="black">", group_ptr->name, var_ptr2->super_array, i);
fortprintf(fd, " %s %% %s %% dimNames(1) = \'num_%s\'</font>
<font color="gray">", group_ptr->name, var_ptr2->super_array, var_ptr2->super_array);
@@ -757,8 +759,14 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
if (!dimlist_ptr->dim->namelist_defined) {
- fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="red">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_code);
- fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_file);
+ if (var_ptr2->persistence == PERSISTENT){
+ fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_code);
+ fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_file);
+ }
+ else {
+ fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s+1</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_code);
+ fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_file);
+ }
}
else {
fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="gray">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_file);
@@ -814,7 +822,6 @@
fortprintf(fd, " %s %% %s %% isSuperArray = .false.</font>
<font color="red">", 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="black">");
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);
@@ -855,8 +862,14 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
if (!dimlist_ptr->dim->namelist_defined) {
- fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="red">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_code);
- fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_file);
+ if(var_ptr->persistence == PERSISTENT){
+ fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_code);
+ fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_file);
+ }
+ else {
+ fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s+1</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_code);
+ fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_file);
+ }
}
else {
fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="gray">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_file);
@@ -999,10 +1012,12 @@
fortprintf(fd, " dest %% %s %% scalar = src %% %s %% scalar</font>
<font color="blue">", var_ptr2->super_array, var_ptr2->super_array);
}
else {
+ if (var_ptr->persistence == PERSISTENT){
if (var_ptr->ndims > 0)
fortprintf(fd, " dest %% %s %% array = src %% %s %% array</font>
<font color="black">", var_ptr->name_in_code, var_ptr->name_in_code);
else
fortprintf(fd, " dest %% %s %% scalar = src %% %s %% scalar</font>
<font color="blue">", var_ptr->name_in_code, var_ptr->name_in_code);
+ }
var_list_ptr = var_list_ptr->next;
}
}
</font>
</pre>