<p><b>mpetersen@lanl.gov</b> 2013-02-13 07:23:01 -0700 (Wed, 13 Feb 2013)</p><p>trunk commit, ocean core only.<br>
<br>
Merge diagnostics_revision branch to trunk. This includes the following changes:<br>
<br>
- Moved diagnostics routines from mpas_ocn_tendency.F to the new<br>
mpas_ocn_diagnostics.F<br>
- Added computation of vertical velocity and time-averaged vertical<br>
velocity. This will be used for the MOC computation in<br>
post-processing.<br>
- Added computation of Brunt-Vaisala frequency.<br>
- Changed names of hmix del2 and del4 flags. <br>
- Changed default values in Registry and namelist to be for standard<br>
120km global run.<br>
<br>
namelist.input.ocean is now created with generate_namelist.py, found<br>
in branches/tools. This revision simply has some spacing differences<br>
due to the python generation tool.<br>
<br>
On the turquiose, create a new namelist using:<br>
<br>
module load python-epd<br>
python generate_namelist.py -f src/core_ocean/Registry<br>
<br>
this makes the file namelist.input.generated.<br>
</p><hr noshade><pre><font color="gray">Index: trunk/mpas
===================================================================
--- trunk/mpas 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas 2013-02-13 14:23:01 UTC (rev 2466)
Property changes on: trunk/mpas
___________________________________________________________________
Modified: svn:mergeinfo
## -4,6 +4,7 ##
/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
\ No newline at end of property
Modified: trunk/mpas/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/namelist.input.ocean 2013-02-13 14:23:01 UTC (rev 2466)
@@ -4,19 +4,18 @@
config_stop_time = 'none'
config_run_duration = '0_06:00:00'
config_calendar_type = '360day'
- config_ncouple_per_day = 1
/
&io
config_input_name = 'grid.nc'
config_output_name = 'output.nc'
config_restart_name = 'restart.nc'
- config_restart_interval = '0_03:00:00'
+ config_restart_interval = '0_06:00:00'
config_output_interval = '0_06:00:00'
- config_stats_interval = '00_00:00:01'
- config_write_stats_on_startup = .false.
- config_write_output_on_startup = .false.
- config_frames_per_outfile = 0
- config_pio_num_iotasks = 0
+ 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
@@ -25,7 +24,6 @@
/
&grid
config_num_halos = 3
- config_enforce_grid_on_restart = .false.
config_vert_coord_movement = 'uniform_stretching'
config_alter_ICs_for_pbcs = 'zlevel_pbcs_off'
config_min_pbc_fraction = 0.10
@@ -38,23 +36,23 @@
config_proc_decomp_file_prefix = 'graph.info.part.'
/
&hmix
- config_h_ScaleWithMesh = .true.
+ config_hmix_ScaleWithMesh = .false.
config_visc_vorticity_term = .true.
config_apvm_scale_factor = 0.0
/
&hmix_del2
config_use_mom_del2 = .false.
config_use_tracer_del2 = .false.
- config_h_mom_eddy_visc2 = 0.0
- config_h_tracer_eddy_diff2 = 0.0
- config_visc_vorticity_visc2_scale = 1.0
+ config_mom_del2 = 0.0
+ config_tracer_del2 = 0.0
+ config_vorticity_del2_scale = 1.0
/
&hmix_del4
config_use_mom_del4 = .true.
config_use_tracer_del4 = .false.
- config_h_mom_eddy_visc4 = 5.0e13
- config_h_tracer_eddy_diff4 = 0.0
- config_visc_vorticity_visc4_scale = 1.0
+ config_mom_del4 = 5.0e13
+ config_tracer_del4 = 0.0
+ config_vorticity_del4_scale = 1.0
/
&hmix_Leith
config_use_Leith_del2 = .false.
@@ -122,11 +120,11 @@
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
+ 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
@@ -143,12 +141,12 @@
config_btr_solve_SSH2 = .false.
/
&sw_model
- config_test_case = 0
+ config_test_case = 0
/
&debug
config_check_zlevel_consistency = .false.
config_filter_btr_mode = .false.
- config_prescribe_velocity = .false.
+ config_prescribe_velocity = .false.
config_prescribe_thickness = .false.
config_include_KE_vertex = .false.
config_check_tracer_monotonicity = .false.
Index: trunk/mpas/src/core_ocean
===================================================================
--- trunk/mpas/src/core_ocean 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean 2013-02-13 14:23:01 UTC (rev 2466)
Property changes on: trunk/mpas/src/core_ocean
___________________________________________________________________
Modified: svn:mergeinfo
## -4,6 +4,7 ##
/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
\ No newline at end of property
Modified: trunk/mpas/src/core_ocean/Makefile
===================================================================
--- trunk/mpas/src/core_ocean/Makefile 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/Makefile 2013-02-13 14:23:01 UTC (rev 2466)
@@ -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: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/Registry 2013-02-13 14:23:01 UTC (rev 2466)
@@ -4,57 +4,55 @@
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 'none'
+namelist character time_management config_run_duration '0_06:00:00'
namelist character time_management config_calendar_type '360day'
-namelist integer time_management 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_restart_interval 'none'
-namelist character io config_output_interval '24:00:00'
-namelist character io config_stats_interval '24:00:00'
-namelist logical io config_write_stats_on_startup .false.
+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 0
+namelist integer io config_frames_per_outfile 1000
namelist integer io config_pio_num_iotasks 0
namelist integer io config_pio_stride 1
-namelist real time_integration config_dt 172.8
-namelist character time_integration config_time_integrator 'RK4'
+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 logical grid config_enforce_grid_on_restart .false.
-namelist character grid config_vert_coord_movement 'isopycnal'
-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 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 hmix config_h_ScaleWithMesh .false.
+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 .true.
-namelist logical hmix_del2 config_use_tracer_del2 .true.
-namelist real hmix_del2 config_h_mom_eddy_visc2 0.0
-namelist real hmix_del2 config_h_tracer_eddy_diff2 0.0
-namelist real hmix_del2 config_visc_vorticity_visc2_scale 1.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 .true.
-namelist real hmix_del4 config_h_mom_eddy_visc4 0.0
-namelist real hmix_del4 config_h_tracer_eddy_diff4 0.0
-namelist real hmix_del4 config_visc_vorticity_visc4_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 0.0
-namelist real hmix_Leith config_Leith_dx 0.0
-namelist real hmix_Leith config_Leith_visc2_max 1000000.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
@@ -66,8 +64,8 @@
namelist real vmix config_convective_visc 1.0
namelist real vmix config_convective_diff 1.0
-namelist logical vmix_const config_use_const_visc .true.
-namelist logical vmix_const config_use_const_diff .true.
+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
@@ -77,8 +75,8 @@
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 .true.
-namelist logical vmix_tanh config_use_tanh_diff .true.
+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
@@ -92,17 +90,17 @@
namelist real forcing config_restoreS_timescale 90.0
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_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_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 'linear'
+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
@@ -111,7 +109,7 @@
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 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
@@ -333,7 +331,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
@@ -366,6 +366,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: trunk/mpas/src/core_ocean/mpas_ocn_diagnostics.F (from rev 2462, branches/ocean_projects/diagnostics_revision/src/core_ocean/mpas_ocn_diagnostics.F)
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_diagnostics.F (rev 0)
+++ trunk/mpas/src/core_ocean/mpas_ocn_diagnostics.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -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: trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -10,6 +10,7 @@
use ocn_test_cases
use ocn_time_integration
use ocn_tendency
+ use ocn_diagnostics
use ocn_monthly_forcing
@@ -102,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)
@@ -1023,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: trunk/mpas/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -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,10 +67,6 @@
!
!--------------------------------------------------------------------
- integer :: ke_cell_flag, ke_vertex_flag
- real (kind=RKIND) :: coef_3rd_order, fCoef
-
-
!***********************************************************************
contains
@@ -238,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)
@@ -268,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
@@ -294,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
@@ -376,723 +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
- 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
- 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_coord_movement.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_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
-
- !
- ! 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, 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 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) = - 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 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.
@@ -1109,28 +383,6 @@
err = 0
- coef_3rd_order = 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_tendency_init!}}}
!***********************************************************************
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_average.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_average.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_average.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -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: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -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
@@ -166,7 +167,7 @@
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: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -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
@@ -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: trunk/mpas/src/core_ocean/mpas_ocn_tracer_advection_std.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -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: trunk/mpas/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -211,9 +211,9 @@
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.
Modified: trunk/mpas/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -233,9 +233,9 @@
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.
Modified: trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del2.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -201,13 +201,12 @@
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
Modified: trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -261,11 +261,11 @@
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
Modified: trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_leith.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2013-02-13 02:11:32 UTC (rev 2465)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2013-02-13 14:23:01 UTC (rev 2466)
@@ -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
</font>
</pre>