<p><b>dwj07@fsu.edu</b> 2013-03-11 08:55:04 -0600 (Mon, 11 Mar 2013)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Merging trunk to cvmix branch.<br>
</p><hr noshade><pre><font color="gray">Index: branches/ocean_projects/cvmix_implementation
===================================================================
--- branches/ocean_projects/cvmix_implementation 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation 2013-03-11 14:55:04 UTC (rev 2580)
Property changes on: branches/ocean_projects/cvmix_implementation
___________________________________________________________________
Modified: svn:mergeinfo
## -4,16 +4,21 ##
/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/explicit_vmix_removal:2486-2490
/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/remove_sw_test_cases:2539-2540
/branches/ocean_projects/restart_reproducibility:2239-2272
+/branches/ocean_projects/sea_level_pressure:2488-2528
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
## -26,3 +31,4 ##
/branches/omp_blocks/multiple_blocks:1803-2084
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+/trunk/mpas:2419-2579
\ No newline at end of property
Modified: branches/ocean_projects/cvmix_implementation/namelist.input.ocean
===================================================================
--- branches/ocean_projects/cvmix_implementation/namelist.input.ocean 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/namelist.input.ocean 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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.
@@ -71,7 +69,6 @@
config_Rayleigh_damping_coeff = 0.0
/
&vmix
- config_implicit_vertical_mix = .true.
config_convective_visc = 1.0
config_convective_diff = 1.0
/
@@ -86,7 +83,7 @@
config_use_rich_diff = .true.
config_bkrd_vert_visc = 1.0e-4
config_bkrd_vert_diff = 1.0e-5
- config_rich_mix = 50.0
+ config_rich_mix = 0.005
/
&vmix_tanh
config_use_tanh_visc = .false.
@@ -121,6 +118,13 @@
&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
@@ -135,13 +139,10 @@
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_velocity = .false.
config_prescribe_thickness = .false.
config_include_KE_vertex = .false.
config_check_tracer_monotonicity = .false.
Modified: branches/ocean_projects/cvmix_implementation/src/core_hyd_atmos/Registry
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_hyd_atmos/Registry 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_hyd_atmos/Registry 2013-03-11 14:55:04 UTC (rev 2580)
@@ -96,7 +96,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 o edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 o localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 o cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 o cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
Modified: branches/ocean_projects/cvmix_implementation/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_init_nhyd_atmos/Registry 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_init_nhyd_atmos/Registry 2013-03-11 14:55:04 UTC (rev 2580)
@@ -108,7 +108,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 io edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 io localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 io cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 io cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 io cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 io verticesOnCell mesh - -
Modified: branches/ocean_projects/cvmix_implementation/src/core_nhyd_atmos/Registry
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_nhyd_atmos/Registry 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_nhyd_atmos/Registry 2013-03-11 14:55:04 UTC (rev 2580)
@@ -113,7 +113,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 iro edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 iro localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 iro cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 iro cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
Index: branches/ocean_projects/cvmix_implementation/src/core_ocean
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean 2013-03-11 14:55:04 UTC (rev 2580)
Property changes on: branches/ocean_projects/cvmix_implementation/src/core_ocean
___________________________________________________________________
Modified: svn:mergeinfo
## -4,16 +4,21 ##
/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/explicit_vmix_removal/src/core_ocean:2486-2490
/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/remove_sw_test_cases/src/core_ocean:2539-2540
/branches/ocean_projects/restart_reproducibility/src/core_ocean:2239-2272
+/branches/ocean_projects/sea_level_pressure/src/core_ocean:2488-2528
/branches/ocean_projects/split_explicit_mrp/src/core_ocean:1134-1138
/branches/ocean_projects/split_explicit_timestepping/src/core_ocean:1044-1097
/branches/ocean_projects/vert_adv_mrp/src/core_ocean:704-745
## -28,3 +33,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:2419-2579
\ No newline at end of property
Modified: branches/ocean_projects/cvmix_implementation/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/Makefile 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/Makefile 2013-03-11 14:55:04 UTC (rev 2580)
@@ -1,11 +1,10 @@
.SUFFIXES: .F .o
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_gm.o \
mpas_ocn_vel_coriolis.o \
mpas_ocn_vel_vadv.o \
mpas_ocn_vel_hmix.o \
@@ -14,7 +13,6 @@
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 \
@@ -36,7 +34,7 @@
mpas_ocn_vmix_coefs_const.o \
mpas_ocn_vmix_coefs_rich.o \
mpas_ocn_vmix_coefs_tanh.o \
- mpas_ocn_vmix_cvmix.o \
+ mpas_ocn_vmix_cvmix.o \
mpas_ocn_restoring.o \
mpas_ocn_tendency.o \
mpas_ocn_tracer_advection.o \
@@ -48,15 +46,16 @@
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_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_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
all: libcvmix dycore
@@ -66,18 +65,18 @@
dycore: libcvmix $(OBJS)
ar -ru libdycore.a $(OBJS) cvmix/*.o
-mpas_ocn_test_cases.o:
-
mpas_ocn_advection.o:
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:
@@ -100,12 +99,10 @@
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_forcing.o: mpas_ocn_vel_forcing_windstress.o mpas_ocn_vel_forcing_rayleigh.o
mpas_ocn_vel_forcing_windstress.o:
-mpas_ocn_vel_forcing_bottomdrag.o:
-
mpas_ocn_vel_forcing_rayleigh.o:
mpas_ocn_vel_coriolis.o:
@@ -178,10 +175,9 @@
mpas_ocn_monthly_forcing.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_mpas_core.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 \
@@ -191,7 +187,6 @@
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 \
@@ -212,7 +207,7 @@
mpas_ocn_vmix_coefs_const.o \
mpas_ocn_vmix_coefs_rich.o \
mpas_ocn_vmix_coefs_tanh.o \
- mpas_ocn_vmix_cvmix.o \
+ mpas_ocn_vmix_cvmix.o \
mpas_ocn_restoring.o \
mpas_ocn_tracer_advection.o \
mpas_ocn_tracer_advection_std.o \
Modified: branches/ocean_projects/cvmix_implementation/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/Registry 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/Registry 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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
@@ -62,12 +60,11 @@
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 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 +74,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,20 +89,26 @@
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
+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 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
@@ -117,8 +120,6 @@
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.
@@ -204,7 +205,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 o edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 o localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 o cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 o cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
@@ -327,7 +328,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
@@ -360,8 +363,12 @@
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 - -
var persistent integer edgeSignOnVertex ( maxEdges nVertices ) 0 - edgeSignOnVertex mesh - -
var persistent integer kiteIndexOnCell ( maxEdges nCells ) 0 - kiteIndexOnCell mesh - -
+
+% Sea surface pressure, for coupling
+var persistent real seaSurfacePressure ( nCells Time ) 0 ir seaSurfacePressure mesh - -
Copied: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_diagnostics.F (from rev 2579, trunk/mpas/src/core_ocean/mpas_ocn_diagnostics.F)
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_diagnostics.F (rev 0)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_diagnostics.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -0,0 +1,848 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! 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, seaSurfacePressure
+ 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
+
+ seaSurfacePressure => grid % seaSurfacePressure % 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.
+ ! Pressure at top surface may be due to atmospheric pressure
+ ! or an ice-shelf depression.
+ pressure(1,iCell) = seaSurfacePressure(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/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_equation_of_state_linear.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_equation_of_state_linear.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_equation_of_state_linear.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_mpas_core.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_mpas_core.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -7,9 +7,9 @@
use mpas_timer
use ocn_global_diagnostics
- use ocn_test_cases
use ocn_time_integration
use ocn_tendency
+ use ocn_diagnostics
use ocn_monthly_forcing
@@ -102,6 +102,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)
@@ -115,15 +117,13 @@
call mpas_dmpar_abort(dminfo)
endif
- if (.not. config_do_restart) call setup_sw_test_case(domain)
-
if (config_vert_coord_movement.ne.'isopycnal') call ocn_init_vert_coord(domain)
call ocn_compute_max_level(domain)
if (.not.config_do_restart) call ocn_init_split_timestep(domain)
- write (0,'(a,a10)') ' Vertical coordinate movement is: ',config_vert_coord_movement
+ write (0,'(a,a)') ' Vertical coordinate movement is: ',trim(config_vert_coord_movement)
if (config_vert_coord_movement.ne.'isopycnal'.and. &
config_vert_coord_movement.ne.'fixed'.and. &
@@ -1023,7 +1023,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/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tendency.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tendency.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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,32 +23,26 @@
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
+ type (timer_node), pointer :: velCorTimer, velVadvTimer, velPgradTimer, velHmixTimer, velForceTimer
+ type (timer_node), pointer :: tracerHadvTimer, tracerVadvTimer, tracerHmixTimer, 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)
@@ -257,28 +240,23 @@
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
- if (.not.config_implicit_vertical_mix) then
- call mpas_timer_start("explicit vmix", .false., velExpVmixTimer)
- call ocn_vel_vmix_tend_explicit(grid, u, h_edge, vertvisctopofedge, tend_u, err)
- call mpas_timer_stop("explicit vmix", velExpVmixTimer)
- endif
call mpas_timer_stop("ocn_tend_u")
end subroutine ocn_tend_u!}}}
!***********************************************************************
!
-! 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 +272,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
@@ -350,17 +328,7 @@
! maxval(tracers(3,1,1:nCells))
! mrp 110516 printing end
- !
- ! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
- !
- if (.not.config_implicit_vertical_mix) then
- call mpas_timer_start("explicit vmix", .false., tracerExpVmixTimer)
- call ocn_tracer_vmix_tend_explicit(grid, h, vertdifftopofcell, tracers, tend_tr, err)
-
- call mpas_timer_stop("explicit vmix", tracerExpVmixTimer)
- endif
-
! mrp 110516 printing
!print *, 'tend_tr 2',minval(tend_tr(3,1,1:nCells)),&
! maxval(tend_tr(3,1,1:nCells))
@@ -376,723 +344,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 +368,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!}}}
!***********************************************************************
Deleted: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_test_cases.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_test_cases.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_test_cases.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -1,526 +0,0 @@
- module ocn_test_cases
-
- use mpas_grid_types
- use mpas_configure
- use mpas_constants
-
-
- contains
-
-
- subroutine setup_sw_test_case(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Configure grid metadata and model state for the shallow water test case
- ! specified in the namelist
- !
- ! Output: block - a subset (not necessarily proper) of the model domain to be
- ! initialized
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: i, iCell, iEdge, iVtx, iLevel
- type (block_type), pointer :: block_ptr
- type (dm_info) :: dminfo
-
- if (config_test_case == 0) then
- write(0,*) 'Using initial conditions supplied in input file'
-
- else if (config_test_case == 1) then
- write(0,*) ' Setting up shallow water test case 1:'
- write(0,*) ' Advection of Cosine Bell over the Pole'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 2) then
- write(0,*) ' Setup shallow water test case 2: '// &
- 'Global Steady State Nonlinear Zonal Geostrophic Flow'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_2(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 5) then
- write(0,*) ' Setup shallow water test case 5:'// &
- ' Zonal Flow over an Isolated Mountain'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_5(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- block_ptr => block_ptr % next
- end do
-
- else if (config_test_case == 6) then
- write(0,*) ' Set up shallow water test case 6:'
- write(0,*) ' Rossby-Haurwitz Wave'
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- call sw_test_case_6(block_ptr % mesh, block_ptr % state % time_levs(1) % state)
- block_ptr => block_ptr % next
- end do
-
- else
- write(0,*) 'Abort: config_test_case=',config_test_case
- write(0,*) 'Only test case 1, 2, 5, and 6 ', &
- 'are currently supported. '
- call mpas_dmpar_abort(dminfo)
- end if
-
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
-
- do i=2,nTimeLevs
- call mpas_copy_state(block_ptr % state % time_levs(i) % state, &
- block_ptr % state % time_levs(1) % state)
- end do
-
- block_ptr => block_ptr % next
- end do
-
- end subroutine setup_sw_test_case
-
-
- subroutine sw_test_case_1(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 1: Advection of Cosine Bell over the Pole
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
- real (kind=RKIND), parameter :: h0 = 1000.0
- real (kind=RKIND), parameter :: theta_c = 0.0
- real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: alpha = pii/4.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: r, u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
- !
- ! Initialize wind field
- !
- allocate(psiVertex(grid % nVertices))
- do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * u0 * ( &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
- cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
- )
- end do
- do iEdge=1,grid % nEdges
- state % u % array(1,iEdge) = -1.0 * ( &
- psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
- psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
- ) / grid%dvEdge%array(iEdge)
- end do
- deallocate(psiVertex)
-
- !
- ! Initialize cosine bell at (theta_c, lambda_c)
- !
- do iCell=1,grid % nCells
- r = sphere_distance(theta_c, lambda_c, grid % latCell % array(iCell), grid % lonCell % array(iCell), a)
- if (r < a/3.0) then
- state % h % array(1,iCell) = (h0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
- else
- state % h % array(1,iCell) = 0.0
- end if
- end do
-
- end subroutine sw_test_case_1
-
-
- subroutine sw_test_case_2(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 2: Global Steady State Nonlinear Zonal
- ! Geostrophic Flow
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 2.0 * pii * a / (12.0 * 86400.0)
- real (kind=RKIND), parameter :: gh0 = 29400.0
- real (kind=RKIND), parameter :: alpha = 0.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
-
- !
- ! Initialize wind field
- !
- allocate(psiVertex(grid % nVertices))
- do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * u0 * ( &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
- cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
- )
- end do
- do iEdge=1,grid % nEdges
- state % u % array(1,iEdge) = -1.0 * ( &
- psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
- psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
- ) / grid%dvEdge%array(iEdge)
- end do
- deallocate(psiVertex)
-
- !
- ! Generate rotated Coriolis field
- !
- do iEdge=1,grid % nEdges
- grid % fEdge % array(iEdge) = 2.0 * omega * &
- ( -cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(alpha) &
- )
- end do
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 2.0 * omega * &
- (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) &
- )
- end do
-
- !
- ! Initialize height field (actually, fluid thickness field)
- !
- do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
- (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
- sin(grid%latCell%array(iCell)) * cos(alpha) &
- )**2.0 &
- ) / &
- gravity
- end do
-
- end subroutine sw_test_case_2
-
-
- subroutine sw_test_case_5(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 5: Zonal Flow over an Isolated Mountain
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: u0 = 20.
- real (kind=RKIND), parameter :: gh0 = 5960.0*gravity
-! real (kind=RKIND), parameter :: hs0 = 2000. original
- real (kind=RKIND), parameter :: hs0 = 250. !mrp 100204
- real (kind=RKIND), parameter :: theta_c = pii/6.0
- real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: rr = pii/9.0
- real (kind=RKIND), parameter :: alpha = 0.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: r, u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
- !
- ! Initialize wind field
- !
- allocate(psiVertex(grid % nVertices))
- do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * u0 * ( &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
- cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
- )
- end do
- do iEdge=1,grid % nEdges
- state % u % array(1,iEdge) = -1.0 * ( &
- psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
- psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
- ) / grid%dvEdge%array(iEdge)
- end do
- deallocate(psiVertex)
-
- !
- ! Generate rotated Coriolis field
- !
- do iEdge=1,grid % nEdges
- grid % fEdge % array(iEdge) = 2.0 * omega * &
- (-cos(grid%lonEdge%array(iEdge)) * cos(grid%latEdge%array(iEdge)) * sin(alpha) + &
- sin(grid%latEdge%array(iEdge)) * cos(alpha) &
- )
- end do
- do iVtx=1,grid % nVertices
- grid % fVertex % array(iVtx) = 2.0 * omega * &
- (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
- sin(grid%latVertex%array(iVtx)) * cos(alpha) &
- )
- end do
-
- !
- ! Initialize mountain
- !
- do iCell=1,grid % nCells
- if (grid % lonCell % array(iCell) < 0.0) grid % lonCell % array(iCell) = grid % lonCell % array(iCell) + 2.0 * pii
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- grid % bottomDepth % array(iCell) = hs0 * (1.0 - r/rr)
- end do
-! output about mountain
-print *, 'bottomDepth',minval(grid % bottomDepth % array),sum(grid % bottomDepth % array)/grid % nCells, maxval(grid % bottomDepth % array)
-
- !
- ! Initialize tracer fields
- !
- do iCell=1,grid % nCells
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- state % tracers % array(1,1,iCell) = 1.0 - r/rr
- end do
- do iCell=1,grid % nCells
- r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + &
- (grid % latCell % array(iCell) - theta_c - pii/6.0)**2.0 &
- ) &
- )
- state % tracers % array(2,1,iCell) = 1.0 - r/rr
- end do
-
- !
- ! Initialize height field (actually, fluid thickness field)
- !
- do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
- (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
- sin(grid%latCell%array(iCell)) * cos(alpha) &
- )**2.0 &
- ) / &
- gravity
- state % h % array(1,iCell) = state % h % array(1,iCell) - grid % bottomDepth % array(iCell)
- end do
-
- end subroutine sw_test_case_5
-
-
- subroutine sw_test_case_6(grid, state)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Setup shallow water test case 6: Rossby-Haurwitz Wave
- !
- ! Reference: Williamson, D.L., et al., "A Standard Test Set for Numerical
- ! Approximations to the Shallow Water Equations in Spherical
- ! Geometry" J. of Comp. Phys., 102, pp. 211--224
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
- type (state_type), intent(inout) :: state
-
- real (kind=RKIND), parameter :: h0 = 8000.0
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- integer :: iCell, iEdge, iVtx
- real (kind=RKIND) :: u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius a
- !
- grid % xCell % array = grid % xCell % array * a
- grid % yCell % array = grid % yCell % array * a
- grid % zCell % array = grid % zCell % array * a
- grid % xVertex % array = grid % xVertex % array * a
- grid % yVertex % array = grid % yVertex % array * a
- grid % zVertex % array = grid % zVertex % array * a
- grid % xEdge % array = grid % xEdge % array * a
- grid % yEdge % array = grid % yEdge % array * a
- grid % zEdge % array = grid % zEdge % array * a
- grid % dvEdge % array = grid % dvEdge % array * a
- grid % dcEdge % array = grid % dcEdge % array * a
- grid % areaCell % array = grid % areaCell % array * a**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * a**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * a**2.0
-
- !
- ! Initialize wind field
- !
- allocate(psiVertex(grid % nVertices))
- do iVtx=1,grid % nVertices
- psiVertex(iVtx) = -a * a * w * sin(grid%latVertex%array(iVtx)) + &
- a *a * K * (cos(grid%latVertex%array(iVtx))**R) * &
- sin(grid%latVertex%array(iVtx)) * cos(R * grid%lonVertex%array(iVtx))
- end do
- do iEdge=1,grid % nEdges
- state % u % array(1,iEdge) = -1.0 * ( &
- psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
- psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
- ) / grid%dvEdge%array(iEdge)
- end do
- deallocate(psiVertex)
-
- !
- ! Initialize height field (actually, fluid thickness field)
- !
- do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gravity * h0 + a*a*aa(grid%latCell%array(iCell)) + &
- a*a*bb(grid%latCell%array(iCell)) * cos(R*grid%lonCell%array(iCell)) + &
- a*a*cc(grid%latCell%array(iCell)) * cos(2.0*R*grid%lonCell%array(iCell)) &
- ) / gravity
- end do
-
- end subroutine sw_test_case_6
-
-
- real function sphere_distance(lat1, lon1, lat2, lon2, radius)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
- ! sphere with given radius.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
- real (kind=RKIND) :: arg1
-
- arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
- cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
- sphere_distance = 2.*radius*asin(arg1)
-
- end function sphere_distance
-
-
- real function aa(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! A, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- real (kind=RKIND), intent(in) :: theta
-
- aa = 0.5 * w * (2.0 * omega + w) * cos(theta)**2.0 + &
- 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*R**2*cos(theta)**(-2.0))
-
- end function aa
-
-
- real function bb(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! B, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- real (kind=RKIND), intent(in) :: theta
-
- bb = (2.0*(omega + w)*K / ((R+1.0)*(R+2.0))) * cos(theta)**R * ((R**2.0 + 2.0*R + 2.0) - ((R+1.0)*cos(theta))**2.0)
-
- end function bb
-
-
- real function cc(theta)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! C, used in height field computation for Rossby-Haurwitz wave
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), parameter :: w = 7.848e-6
- real (kind=RKIND), parameter :: K = 7.848e-6
- real (kind=RKIND), parameter :: R = 4.0
-
- real (kind=RKIND), intent(in) :: theta
-
- cc = 0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 - R - 2.0)
-
- end function cc
-
-end module ocn_test_cases
Modified: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_average.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_average.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_average.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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
@@ -148,11 +149,6 @@
call mpas_timer_start("RK4-tendency computations")
block => domain % blocklist
do while (associated(block))
-
- if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, err)
- end if
-
! advection of u uses u, while advection of h and tracers use uTransport.
call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
block % provis % u % array, block % provis % wTop % array, err)
@@ -166,7 +162,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")
@@ -200,10 +196,6 @@
end do
end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- end if
-
if (config_prescribe_velocity) then
block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
@@ -254,7 +246,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")
@@ -270,42 +262,35 @@
block => block % next
end do
+ call mpas_timer_start("RK4-implicit vert mix")
+ block => domain % blocklist
+ do while(associated(block))
- if (config_implicit_vertical_mix) then
- call mpas_timer_start("RK4-implicit vert mix")
- block => domain % blocklist
- do while(associated(block))
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- ! Call ocean diagnostic solve in preparation for vertical mixing. Note
- ! it is called again after vertical mixing, because u and tracers change.
- ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
- ! be computed. For kpp, more variables may be needed. Either way, this
- ! could be made more efficient by only computing what is needed for the
- ! implicit vmix routine that follows. mrp 121023.
- call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ block => block % next
+ end do
- call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
- block => block % next
- end do
+ ! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
+ ! this leads to lack of volume conservation. It is required because halo updates in RK4 are only
+ ! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
+ ! communicate the change due to implicit vertical mixing across the boundary.
+ call mpas_timer_start("RK4-implicit vert mix halos")
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
+ call mpas_timer_stop("RK4-implicit vert mix halos")
- ! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
- ! this leads to lack of volume conservation. It is required because halo updates in RK4 are only
- ! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
- ! communicate the change due to implicit vertical mixing across the boundary.
- call mpas_timer_start("RK4-implicit vert mix halos")
- call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
- call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
- call mpas_timer_stop("RK4-implicit vert mix halos")
+ call mpas_timer_stop("RK4-implicit vert mix")
- call mpas_timer_stop("RK4-implicit vert mix")
- end if
-
block => domain % blocklist
do while (associated(block))
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- end if
-
if (config_prescribe_velocity) then
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
Modified: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_integration_split.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_time_integration_split.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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
@@ -187,10 +188,6 @@
stage1_tend_time = min(split_explicit_step,2)
- if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, err)
- end if
-
! compute wTop. Use u (rather than uTransport) for momentum advection.
! Use the most recent time level available.
call ocn_wtop(block % mesh, block % state % time_levs(stage1_tend_time) % state % h % array, &
@@ -802,7 +799,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
@@ -939,43 +936,37 @@
! END large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- if (config_implicit_vertical_mix) then
- call mpas_timer_start("se implicit vert mix")
- block => domain % blocklist
- do while(associated(block))
+ call mpas_timer_start("se implicit vert mix")
+ block => domain % blocklist
+ do while(associated(block))
- ! Call ocean diagnostic solve in preparation for vertical mixing. Note
- ! it is called again after vertical mixing, because u and tracers change.
- ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
- ! be computed. For kpp, more variables may be needed. Either way, this
- ! could be made more efficient by only computing what is needed for the
- ! implicit vmix routine that follows. mrp 121023.
- call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ ! Call ocean diagnostic solve in preparation for vertical mixing. Note
+ ! it is called again after vertical mixing, because u and tracers change.
+ ! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
+ ! be computed. For kpp, more variables may be needed. Either way, this
+ ! could be made more efficient by only computing what is needed for the
+ ! implicit vmix routine that follows. mrp 121023.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
- block => block % next
- end do
+ block => block % next
+ end do
- ! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
- ! this leads to lack of volume conservation. It is required because halo updates in stage 3 are only
- ! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
- ! communicate the change due to implicit vertical mixing across the boundary.
- call mpas_timer_start("se implicit vert mix halos")
- call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
- call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
- call mpas_timer_stop("se implicit vert mix halos")
+ ! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
+ ! this leads to lack of volume conservation. It is required because halo updates in stage 3 are only
+ ! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
+ ! communicate the change due to implicit vertical mixing across the boundary.
+ call mpas_timer_start("se implicit vert mix halos")
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
+ call mpas_timer_stop("se implicit vert mix halos")
- call mpas_timer_stop("se implicit vert mix")
- end if
+ call mpas_timer_stop("se implicit vert mix")
block => domain % blocklist
do while (associated(block))
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- end if
-
if (config_prescribe_velocity) then
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
Modified: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_forcing.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_forcing.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_forcing.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -18,7 +18,6 @@
use mpas_configure
use ocn_vel_forcing_windstress
- use ocn_vel_forcing_bottomdrag
use ocn_vel_forcing_rayleigh
implicit none
@@ -115,7 +114,7 @@
!
!-----------------------------------------------------------------
- integer :: err1, err2, err3
+ integer :: err1, err2
!-----------------------------------------------------------------
!
@@ -126,11 +125,9 @@
!-----------------------------------------------------------------
call ocn_vel_forcing_windstress_tend(grid, u_src, h_edge, tend, err1)
- call ocn_vel_forcing_bottomdrag_tend(grid, u, ke_edge, h_edge, tend, err2)
- call ocn_vel_forcing_rayleigh_tend(grid, u, tend, err3)
+ call ocn_vel_forcing_rayleigh_tend(grid, u, tend, err2)
err = ior(err1, err2)
- err = ior(err, err3)
!--------------------------------------------------------------------
@@ -164,14 +161,12 @@
integer, intent(out) :: err !< Output: error flag
- integer :: err1, err2, err3
+ integer :: err1, err2
call ocn_vel_forcing_windstress_init(err1)
- call ocn_vel_forcing_bottomdrag_init(err2)
- call ocn_vel_forcing_rayleigh_init(err3)
+ call ocn_vel_forcing_rayleigh_init(err2)
err = ior(err1, err2)
- err = ior(err, err3)
!--------------------------------------------------------------------
Deleted: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -1,193 +0,0 @@
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-!
-! ocn_vel_forcing_bottomdrag
-!
-!> \brief MPAS ocean bottom drag
-!> \author Doug Jacobsen
-!> \date 16 September 2011
-!> \version SVN:$Id:$
-!> \details
-!> This module contains the routine for computing
-!> tendencies from bottom drag.
-!
-!-----------------------------------------------------------------------
-
-module ocn_vel_forcing_bottomdrag
-
- use mpas_grid_types
- use mpas_configure
-
- implicit none
- private
- save
-
- !--------------------------------------------------------------------
- !
- ! Public parameters
- !
- !--------------------------------------------------------------------
-
- !--------------------------------------------------------------------
- !
- ! Public member functions
- !
- !--------------------------------------------------------------------
-
- public :: ocn_vel_forcing_bottomdrag_tend, &
- ocn_vel_forcing_bottomdrag_init
-
- !--------------------------------------------------------------------
- !
- ! Private module variables
- !
- !--------------------------------------------------------------------
-
- logical :: bottomDragOn
- real (kind=RKIND) :: bottomDragCoef
-
-
-!***********************************************************************
-
-contains
-
-!***********************************************************************
-!
-! routine ocn_vel_forcing_bottomdrag_tend
-!
-!> \brief Computes tendency term from bottom drag
-!> \author Doug Jacobsen
-!> \date 15 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the bottom drag tendency for momentum
-!> based on current state.
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_vel_forcing_bottomdrag_tend(grid, u, ke_edge, h_edge, tend, err)!{{{
-
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- u !< Input: velocity
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- ke_edge !< Input: kinetic energy at edge
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- h_edge !< Input: thickness at edge
-
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
-
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:), intent(inout) :: &
- tend !< Input/Output: velocity tendency
-
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
-
- integer :: iEdge, nEdgesSolve, k
- integer, dimension(:), pointer :: maxLevelEdgeTop
- integer, dimension(:,:), pointer :: edgeMask
-
- !-----------------------------------------------------------------
- !
- ! call relevant routines for computing tendencies
- ! note that the user can choose multiple options and the
- ! tendencies will be added together
- !
- !-----------------------------------------------------------------
-
- err = 0
-
- if(.not.bottomDragOn) return
-
- nEdgesSolve = grid % nEdgesSolve
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- edgeMask => grid % edgeMask % array
-
- do iEdge=1,grid % nEdgesSolve
-
- k = max(maxLevelEdgeTop(iEdge), 1)
-
- ! bottom drag is the same as POP:
- ! -c |u| u where c is unitless and 1.0e-3.
- ! see POP Reference guide, section 3.4.4.
-
- tend(k,iEdge) = tend(k,iEdge)-edgeMask(k,iEdge)*(bottomDragCoef*u(k,iEdge)*sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge))
-
- enddo
-
-
-
- !--------------------------------------------------------------------
-
- end subroutine ocn_vel_forcing_bottomdrag_tend!}}}
-
-!***********************************************************************
-!
-! routine ocn_vel_forcing_bottomdrag_init
-!
-!> \brief Initializes ocean bottom drag
-!> \author Doug Jacobsen
-!> \date 16 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine initializes quantities related to bottom drag
-!> in the ocean.
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_vel_forcing_bottomdrag_init(err)!{{{
-
- !--------------------------------------------------------------------
-
- !-----------------------------------------------------------------
- !
- ! call individual init routines for each parameterization
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
-
- err = 0
-
- bottomDragOn = .false.
-
- if (.not.config_implicit_vertical_mix) then
- bottomDragOn = .true.
- bottomDragCoef = config_bottom_drag_coeff
- endif
-
- !--------------------------------------------------------------------
-
- end subroutine ocn_vel_forcing_bottomdrag_init!}}}
-
-!***********************************************************************
-
-end module ocn_vel_forcing_bottomdrag
-
-!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
-! vim: foldmethod=marker
Modified: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_del2.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_leith.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vmix.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vmix.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -45,8 +45,6 @@
tridiagonal_solve_mult
public :: ocn_vmix_coefs, &
- ocn_vel_vmix_tend_explicit, &
- ocn_tracer_vmix_tend_explicit, &
ocn_vel_vmix_tend_implicit, &
ocn_tracer_vmix_tend_implicit, &
ocn_vmix_init, &
@@ -59,7 +57,6 @@
!--------------------------------------------------------------------
logical :: velVmixOn, tracerVmixOn
- logical :: explicitOn, implicitOn
!***********************************************************************
@@ -141,100 +138,6 @@
!***********************************************************************
!
-! routine ocn_vel_vmix_tendExplict
-!
-!> \brief Computes tendencies for explict momentum vertical mixing
-!> \author Doug Jacobsen
-!> \date 19 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the tendencies for explicit vertical mixing for momentum
-!> using computed coefficients.
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_vel_vmix_tend_explicit(grid, u, h_edge, vertViscTopOfEdge, tend, err)!{{{
-
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
-
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- u !< Input: velocity
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- h_edge !< Input: thickness at edge
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- vertViscTopOfEdge !< Input: vertical mixing coefficients
-
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:), intent(inout) :: &
- tend !< Input/Output: tendency information
-
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
-
- integer :: iEdge, nEdgesSolve, k, nVertLevels
-
- integer, dimension(:), pointer :: maxLevelEdgeTop
-
- real (kind=RKIND), dimension(:), allocatable :: fluxVertTop
-
- err = 0
-
- if(.not.velVmixOn) return
- if(implicitOn) return
-
- nEdgessolve = grid % nEdgesSolve
- nVertLevels = grid % nVertLevels
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
-
- allocate(fluxVertTop(nVertLevels+1))
- fluxVertTop(1) = 0.0
- do iEdge=1,nEdgesSolve
- do k=2,maxLevelEdgeTop(iEdge)
- fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &
- * ( u(k-1,iEdge) - u(k,iEdge) ) &
- * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
- enddo
- fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
-
- do k=1,maxLevelEdgeTop(iEdge)
- tend(k,iEdge) = tend(k,iEdge) &
- + (fluxVertTop(k) - fluxVertTop(k+1)) &
- / h_edge(k,iEdge)
- enddo
-
- end do
- deallocate(fluxVertTop)
- !--------------------------------------------------------------------
-
- end subroutine ocn_vel_vmix_tend_explicit!}}}
-
-!***********************************************************************
-!
! routine ocn_vel_vmix_tend_implicit
!
!> \brief Computes tendencies for implicit momentum vertical mixing
@@ -307,7 +210,6 @@
err = 0
if(.not.velVmixOn) return
- if(explicitOn) return
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
@@ -370,111 +272,6 @@
!***********************************************************************
!
-! routine ocn_tracer_vmix_tendExplict
-!
-!> \brief Computes tendencies for explict tracer vertical mixing
-!> \author Doug Jacobsen
-!> \date 19 September 2011
-!> \version SVN:$Id$
-!> \details
-!> This routine computes the tendencies for explicit vertical mixing for
-!> tracers using computed coefficients.
-!
-!-----------------------------------------------------------------------
-
- subroutine ocn_tracer_vmix_tend_explicit(grid, h, vertDiffTopOfCell, tracers, tend, err)!{{{
-
- !-----------------------------------------------------------------
- !
- ! input variables
- !
- !-----------------------------------------------------------------
-
- type (mesh_type), intent(in) :: &
- grid !< Input: grid information
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- h !< Input: thickness at cell center
-
- real (kind=RKIND), dimension(:,:), intent(in) :: &
- vertDiffTopOfCell !< Input: vertical mixing coefficients
-
- real (kind=RKIND), dimension(:,:,:), intent(in) :: &
- tracers !< Input: tracers
-
- !-----------------------------------------------------------------
- !
- ! input/output variables
- !
- !-----------------------------------------------------------------
-
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
- tend !< Input/Output: tendency information
-
- !-----------------------------------------------------------------
- !
- ! output variables
- !
- !-----------------------------------------------------------------
-
- integer, intent(out) :: err !< Output: error flag
-
- !-----------------------------------------------------------------
- !
- ! local variables
- !
- !-----------------------------------------------------------------
-
- integer :: iCell, nCellsSolve, k, iTracer, num_tracers, nVertLevels
-
- integer, dimension(:), pointer :: maxLevelCell
-
- real (kind=RKIND), dimension(:,:), allocatable :: fluxVertTop
-
- err = 0
-
- if(.not.tracerVmixOn) return
- if(implicitOn) return
-
- nCellsSolve = grid % nCellsSolve
- nVertLevels = grid % nVertLevels
- num_tracers = size(tracers, dim=1)
-
- maxLevelCell => grid % maxLevelCell % array
-
- allocate(fluxVertTop(num_tracers,nVertLevels+1))
- fluxVertTop(:,1) = 0.0
- do iCell=1,nCellsSolve
-
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! compute \kappa_v d\phi/dz
- fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
- * ( tracers(iTracer,k-1,iCell) &
- - tracers(iTracer,k ,iCell) ) &
- * 2 / (h(k-1,iCell) + h(k,iCell))
-
- enddo
- enddo
- fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
-
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! This is h d/dz( fluxVertTop) but h and dz cancel, so
- ! reduces to delta( fluxVertTop)
- tend(iTracer,k,iCell) = tend(iTracer,k,iCell) &
- + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
- enddo
- enddo
-
- enddo ! iCell loop
- deallocate(fluxVertTop)
- !--------------------------------------------------------------------
-
- end subroutine ocn_tracer_vmix_tend_explicit!}}}
-
-!***********************************************************************
-!
! routine ocn_tracer_vmix_tend_implicit
!
!> \brief Computes tendencies for implicit tracer vertical mixing
@@ -540,7 +337,6 @@
err = 0
if(.not.tracerVmixOn) return
- if(explicitOn) return
nCells = grid % nCells
nVertLevels = grid % nVertLevels
@@ -650,8 +446,7 @@
!> \version SVN:$Id$
!> \details
!> This routine initializes a variety of quantities related to
-!> vertical mixing in the ocean. This primarily determines if
-!> explicit or implicit vertical mixing is to be used.
+!> vertical mixing in the ocean.
!
!-----------------------------------------------------------------------
@@ -675,14 +470,6 @@
velVmixOn = .true.
tracerVmixOn = .true.
- explicitOn = .true.
- implicitOn = .false.
-
- if(config_implicit_vertical_mix) then
- explicitOn = .false.
- implicitOn = .true.
- end if
-
if(config_disable_u_vmix) velVmixOn = .false.
if(config_disable_tr_vmix) tracerVmixOn = .false.
Modified: branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -228,26 +228,14 @@
if (RiTopOfEdge(k,iEdge)>0.0) then
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
+ ! maltrud do limiting of coefficient--should not be necessary
+ ! also probably better logic could be found
if (vertViscTopOfEdge(k,iEdge) > config_convective_visc) then
- if( config_implicit_vertical_mix) then
- vertViscTopOfEdge(k,iEdge) = config_convective_visc
- else
- vertViscTopOfEdge(k,iEdge) = &
- ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
- end if
+ vertViscTopOfEdge(k,iEdge) = config_convective_visc
end if
else
- ! mrp 110324 efficiency note: this if is inside iCell and k loops.
- if (config_implicit_vertical_mix) then
- ! for Ri<0 and implicit mix, use convective diffusion
- vertViscTopOfEdge(k,iEdge) = config_convective_visc
- else
- ! for Ri<0 and explicit vertical mix,
- ! use maximum diffusion allowed by CFL criterion
- vertViscTopOfEdge(k,iEdge) = vertViscTopOfEdge(k,iEdge) + ((h_edge(k-1,iEdge)+h_edge(k,iEdge))/2.0)**2/config_dt/4.0
- end if
+ ! for Ri<0 use the convective value for the viscosity
+ vertViscTopOfEdge(k,iEdge) = config_convective_visc
end if
end do
end do
@@ -336,23 +324,11 @@
! maltrud do limiting of coefficient--should not be necessary
! also probably better logic could be found
if (vertDiffTopOfCell(k,iCell) > config_convective_diff) then
- if (config_implicit_vertical_mix) then
- vertDiffTopOfCell(k,iCell) = config_convective_diff
- else
- vertDiffTopOfCell(k,iCell) = &
- ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
- end if
+ vertDiffTopOfCell(k,iCell) = config_convective_diff
end if
else
- ! mrp 110324 efficiency note: this if is inside iCell and k loops.
- if (config_implicit_vertical_mix) then
- ! for Ri<0 and implicit mix, use convective diffusion
- vertDiffTopOfCell(k,iCell) = config_convective_diff
- else
- ! for Ri<0 and explicit vertical mix,
- ! use maximum diffusion allowed by CFL criterion
- vertDiffTopOfCell(k,iCell) = vertDiffTopOfCell(k, iCell) + ((h(k-1,iCell)+h(k,iCell))/2.0)**2/config_dt/4.0
- end if
+ ! for Ri<0 use the convective value for the diffusion
+ vertDiffTopOfCell(k,iCell) = config_convective_diff
end if
end do
end do
Modified: branches/ocean_projects/cvmix_implementation/src/core_sw/Registry
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/core_sw/Registry 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/core_sw/Registry 2013-03-11 14:55:04 UTC (rev 2580)
@@ -97,7 +97,7 @@
var persistent real edgeNormalVectors ( R3 nEdges ) 0 o edgeNormalVectors mesh - -
var persistent real localVerticalUnitVectors ( R3 nCells ) 0 o localVerticalUnitVectors mesh - -
-var persistent real cellTangentPlane ( R3 TWO nEdges ) 0 o cellTangentPlane mesh - -
+var persistent real cellTangentPlane ( R3 TWO nCells ) 0 o cellTangentPlane mesh - -
var persistent integer cellsOnCell ( maxEdges nCells ) 0 iro cellsOnCell mesh - -
var persistent integer verticesOnCell ( maxEdges nCells ) 0 iro verticesOnCell mesh - -
Modified: branches/ocean_projects/cvmix_implementation/src/framework/mpas_block_decomp.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/framework/mpas_block_decomp.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/framework/mpas_block_decomp.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -165,7 +165,7 @@
allocate(block_id(blocks_per_proc))
allocate(block_start(blocks_per_proc))
allocate(block_count(blocks_per_proc))
-
+
do i = 1, blocks_per_proc
block_start = 0
block_count = 0
@@ -436,11 +436,11 @@
end if
else
blocks_per_proc = 0
- do i = 1, total_blocks
+ do i = 0, total_blocks-1
call mpas_get_owning_proc(dminfo, i, owning_proc)
if(owning_proc == proc_number) then
call mpas_get_local_block_id(dminfo, i, local_block_id)
- blocks_per_proc = max(blocks_per_proc, local_block_id)
+ blocks_per_proc = max(blocks_per_proc, local_block_id+1)
end if
end do
end if
Modified: branches/ocean_projects/cvmix_implementation/src/framework/mpas_dmpar.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/framework/mpas_dmpar.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/framework/mpas_dmpar.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -25,6 +25,8 @@
module procedure mpas_dmpar_alltoall_field1d_real
module procedure mpas_dmpar_alltoall_field2d_real
module procedure mpas_dmpar_alltoall_field3d_real
+ module procedure mpas_dmpar_alltoall_field4d_real
+ module procedure mpas_dmpar_alltoall_field5d_real
end interface
private :: mpas_dmpar_alltoall_field1d_integer
@@ -32,6 +34,8 @@
private :: mpas_dmpar_alltoall_field1d_real
private :: mpas_dmpar_alltoall_field2d_real
private :: mpas_dmpar_alltoall_field3d_real
+ private :: mpas_dmpar_alltoall_field4d_real
+ private :: mpas_dmpar_alltoall_field5d_real
interface mpas_dmpar_exch_halo_field
@@ -41,6 +45,8 @@
module procedure mpas_dmpar_exch_halo_field1d_real
module procedure mpas_dmpar_exch_halo_field2d_real
module procedure mpas_dmpar_exch_halo_field3d_real
+ module procedure mpas_dmpar_exch_halo_field4d_real
+ module procedure mpas_dmpar_exch_halo_field5d_real
end interface
private :: mpas_dmpar_exch_halo_field1d_integer
@@ -49,6 +55,8 @@
private :: mpas_dmpar_exch_halo_field1d_real
private :: mpas_dmpar_exch_halo_field2d_real
private :: mpas_dmpar_exch_halo_field3d_real
+ private :: mpas_dmpar_exch_halo_field4d_real
+ private :: mpas_dmpar_exch_halo_field5d_real
interface mpas_dmpar_copy_field
module procedure mpas_dmpar_copy_field1d_integer
@@ -57,6 +65,8 @@
module procedure mpas_dmpar_copy_field1d_real
module procedure mpas_dmpar_copy_field2d_real
module procedure mpas_dmpar_copy_field3d_real
+ module procedure mpas_dmpar_copy_field4d_real
+ module procedure mpas_dmpar_copy_field5d_real
end interface
private :: mpas_dmpar_copy_field1d_integer
@@ -65,6 +75,8 @@
private :: mpas_dmpar_copy_field1d_real
private :: mpas_dmpar_copy_field2d_real
private :: mpas_dmpar_copy_field3d_real
+ private :: mpas_dmpar_copy_field4d_real
+ private :: mpas_dmpar_copy_field5d_real
contains
@@ -2810,7 +2822,610 @@
end subroutine mpas_dmpar_alltoall_field3d_real!}}}
+ subroutine mpas_dmpar_alltoall_field4d_real(fieldIn, fieldout, haloLayersIn)!{{{
+ implicit none
+
+ type (field4dReal), pointer :: fieldIn
+ type (field4dReal), pointer :: fieldOut
+ integer, dimension(:), pointer, optional :: haloLayersIn
+
+ type (field4dReal), pointer :: fieldInPtr, fieldOutPtr
+ type (mpas_exchange_list), pointer :: exchListPtr, exchListPtr2
+ type (mpas_communication_list), pointer :: sendList, recvList, commListPtr, commListPtr2
+ type (dm_info), pointer :: dminfo
+
+ logical :: comm_list_found
+
+ integer :: lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ integer :: nAdded, bufferOffset
+ integer :: mpi_ierr
+ integer :: iHalo, iBuffer, i, j, k, l
+ integer :: nHaloLayers
+ integer, dimension(:), pointer :: haloLayers
+
+ dminfo => fieldIn % block % domain % dminfo
+
+ if(present(haloLayersIn)) then
+ nHaloLayers = size(haloLayersIn)
+ allocate(haloLayers(nHaloLayers))
+ do iHalo = 1, nHaloLayers
+ haloLayers(iHalo) = haloLayersIn(iHalo)
+ end do
+ else
+ nHaloLayers = size(fieldIn % sendList % halos)
+ allocate(haloLayers(nHaloLayers))
+ do iHalo = 1, nHaloLayers
+ haloLayers(iHalo) = iHalo
+ end do
+ end if
+
+#ifdef _MPI
+ nullify(sendList)
+ nullify(recvList)
+
+ ! Setup recieve lists.
+ do iHalo = 1, nHaloLayers
+ fieldOutPtr => fieldOut
+ do while(associated(fieldOutPtr))
+ exchListPtr => fieldOutPtr % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ comm_list_found = .false.
+
+ ! Search for an already created commList to this processor.
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ if(commListPtr % procID == exchListPtr % endPointID) then
+ commListPtr % nList = commListPtr % nList + exchListPtr % nList * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3)
+ comm_list_found = .true.
+ exit
+ end if
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! If no comm list exists, create a new one.
+ if(.not. comm_list_found) then
+ if(.not.associated(recvList)) then
+ allocate(recvList)
+ nullify(recvList % next)
+ commListPtr => recvList
+ else
+ commListPtr => recvList
+ commListPtr2 => commListPtr % next
+ do while(associated(commListPtr2))
+ commListPtr => commListPtr % next
+ commListPtr2 => commListPtr % next
+ end do
+
+ allocate(commListPtr % next)
+ commListPtr => commListPtr % next
+ nullify(commListPtr % next)
+ end if
+
+ commListPtr % procID = exchListPtr % endPointID
+ commListPtr % nList = exchListPtr % nList * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3)
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldOutPtr => fieldOutPtr % next
+ end do
+ end do
+
+ ! Determine size of receive list buffers.
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldOutPtr => fieldOut
+ do while(associated(fieldOutPtr))
+ exchListPtr => fieldOutPtr % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ nAdded = max(nAdded, maxval(exchListPtr % srcList) * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3))
+ end if
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldOutPtr => fieldOutPtr % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+ commListPtr % nList = nAdded
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! Allocate buffers for recieves, and initiate mpi_irecv calls.
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ allocate(commListPtr % rbuffer(commListPtr % nList))
+ nullify(commListPtr % ibuffer)
+ call MPI_Irecv(commListPtr % rbuffer, commListPtr % nList, MPI_realKIND, commListPtr % procID, commListPtr % procID, dminfo % comm, commListPtr % reqID, mpi_ierr)
+ commListPtr => commListPtr % next
+ end do
+
+ ! Setup send lists, and determine the size of their buffers.
+ do iHalo = 1, nHaloLayers
+ fieldInPtr => fieldIn
+ do while(associated(fieldInPtr))
+ exchListPtr => fieldInPtr % sendList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ comm_list_found = .false.
+
+ ! Search for an already created commList to this processor.
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ if(commListPtr % procID == exchListPtr % endPointID) then
+ commListPtr % nList = commListPtr % nList + exchListPtr % nList * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) * fieldOutPtr % dimSizes(3)
+ comm_list_found = .true.
+ exit
+ end if
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! If no comm list exists, create a new one.
+ if(.not. comm_list_found) then
+ if(.not.associated(sendList)) then
+ allocate(sendList)
+ nullify(sendList % next)
+ commListPtr => sendList
+ else
+ commListPtr => sendList
+ commListPtr2 => commListPtr % next
+ do while(associated(commListPtr2))
+ commListPtr => commListPtr % next
+ commListPtr2 => commListPtr % next
+ end do
+
+ allocate(commListPtr % next)
+ commListPtr => commListPtr % next
+ nullify(commListPtr % next)
+ end if
+ commListPtr % procID = exchListPtr % endPointID
+ commListPtr % nList = exchListPtr % nList * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) * fieldOutPtr % dimSizes(3)
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldInPtr => fieldInPtr % next
+ end do
+ end do
+
+ ! Allocate sendLists, copy data into buffer, and initiate mpi_isends
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ allocate(commListPtr % rbuffer(commListPtr % nList))
+ nullify(commListPtr % ibuffer)
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldInPtr => fieldIn
+ do while(associated(fieldInPtr))
+ exchListPtr => fieldInPtr % sendList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ do i = 1, exchListPtr % nList
+ do j = 1, fieldInPtr % dimSizes(3)
+ do k = 1, fieldInPtr % dimSizes(2)
+ do l = 1, fieldInPtr % dimSizes(1)
+ iBuffer = (exchListPtr % destList(i)-1) * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) * fieldInPtr % dimSizes(3) &
+ + (j-1) * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) &
+ + (k-1) * fieldInPtr % dimSizes(1) + l + bufferOffset
+ commListPtr % rbuffer(iBuffer) = fieldInPtr % array(l, k, j, exchListPtr % srcList(i))
+ nAdded = nAdded + 1
+ end do
+ end do
+ end do
+ end do
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldInPtr => fieldInPtr % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+
+ call MPI_Isend(commListPtr % rbuffer, commListPtr % nlist, MPI_realKIND, &
+ commListPtr % procID, dminfo % my_proc_id, dminfo % comm, commListPtr % reqID, mpi_ierr)
+
+ commListPtr => commListPtr % next
+ end do
+
+#endif
+
+ ! Handle Local Copies. Only local copies if no MPI
+ do iHalo = 1, nHaloLayers
+ fieldInPtr => fieldIn
+ do while(associated(fieldInPtr))
+ exchListPtr => fieldInPtr % copyList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ fieldOutPtr => fieldOut
+ do while(associated(fieldOutPtr))
+ if(exchListPtr % endPointID == fieldOutPtr % block % localBlockID) then
+ do i = 1, exchListPtr % nList
+ fieldOutPtr % array(:, :, :, exchListPtr % destList(i)) = fieldInPtr % array(:, :, :, exchListPtr % srcList(i))
+ end do
+ end if
+ fieldOutPtr => fieldOutPtr % next
+ end do
+
+ exchListPtr => exchListPtr % next
+ end do
+ fieldInPtr => fieldInPtr % next
+ end do
+ end do
+
+#ifdef _MPI
+ ! Wait for MPI_Irecv's to finish, and unpack data.
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldOutPtr => fieldOut
+ do while(associated(fieldOutPtr))
+ exchListPtr => fieldOutPtr % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ do i = 1, exchListPtr % nList
+ do j = 1, fieldOutPtr % dimSizes(3)
+ do k = 1, fieldOutPtr % dimSizes(2)
+ do l = 1, fieldOutPtr % dimSizes(1)
+ iBuffer = (exchListPtr % srcList(i)-1) * fieldOutPtr % dimSizes(3) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(1) &
+ + (j-1) * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) &
+ + (k-1) * fieldOutPtr % dimSizes(1) + l + bufferOffset
+ fieldOutPtr % array(l, k, j, exchListPtr % destList(i)) = commListPtr % rbuffer(iBuffer)
+ end do
+ end do
+ end do
+ end do
+ nAdded = max(nAdded, maxval(exchListPtr % srcList) * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3))
+ end if
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldOutPtr => fieldOutPtr % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! Wait for MPI_Isend's to finish.
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ commListPtr => commListPtr % next
+ end do
+
+ ! Destroy commLists.
+ call mpas_dmpar_destroy_communication_list(sendList)
+ call mpas_dmpar_destroy_communication_list(recvList)
+#endif
+
+ deallocate(haloLayers)
+
+ end subroutine mpas_dmpar_alltoall_field4d_real!}}}
+
+ subroutine mpas_dmpar_alltoall_field5d_real(fieldIn, fieldout, haloLayersIn)!{{{
+
+ implicit none
+
+ type (field5dReal), pointer :: fieldIn
+ type (field5dReal), pointer :: fieldOut
+ integer, dimension(:), pointer, optional :: haloLayersIn
+
+ type (field5dReal), pointer :: fieldInPtr, fieldOutPtr
+ type (mpas_exchange_list), pointer :: exchListPtr, exchListPtr2
+ type (mpas_communication_list), pointer :: sendList, recvList, commListPtr, commListPtr2
+ type (dm_info), pointer :: dminfo
+
+ logical :: comm_list_found
+
+ integer :: lastPackedIdx, lastUnpackedIdx, nPacked, nUnpacked
+ integer :: nAdded, bufferOffset
+ integer :: mpi_ierr
+ integer :: iHalo, iBuffer, i, j, k, l, m
+ integer :: nHaloLayers
+ integer, dimension(:), pointer :: haloLayers
+
+ dminfo => fieldIn % block % domain % dminfo
+
+ if(present(haloLayersIn)) then
+ nHaloLayers = size(haloLayersIn)
+ allocate(haloLayers(nHaloLayers))
+ do iHalo = 1, nHaloLayers
+ haloLayers(iHalo) = haloLayersIn(iHalo)
+ end do
+ else
+ nHaloLayers = size(fieldIn % sendList % halos)
+ allocate(haloLayers(nHaloLayers))
+ do iHalo = 1, nHaloLayers
+ haloLayers(iHalo) = iHalo
+ end do
+ end if
+
+#ifdef _MPI
+ nullify(sendList)
+ nullify(recvList)
+
+ ! Setup recieve lists.
+ do iHalo = 1, nHaloLayers
+ fieldOutPtr => fieldOut
+ do while(associated(fieldOutPtr))
+ exchListPtr => fieldOutPtr % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ comm_list_found = .false.
+
+ ! Search for an already created commList to this processor.
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ if(commListPtr % procID == exchListPtr % endPointID) then
+ commListPtr % nList = commListPtr % nList + exchListPtr % nList * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3) * fieldOutPtr % dimSizes(4)
+ comm_list_found = .true.
+ exit
+ end if
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! If no comm list exists, create a new one.
+ if(.not. comm_list_found) then
+ if(.not.associated(recvList)) then
+ allocate(recvList)
+ nullify(recvList % next)
+ commListPtr => recvList
+ else
+ commListPtr => recvList
+ commListPtr2 => commListPtr % next
+ do while(associated(commListPtr2))
+ commListPtr => commListPtr % next
+ commListPtr2 => commListPtr % next
+ end do
+
+ allocate(commListPtr % next)
+ commListPtr => commListPtr % next
+ nullify(commListPtr % next)
+ end if
+
+ commListPtr % procID = exchListPtr % endPointID
+ commListPtr % nList = exchListPtr % nList * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3) * fieldOutPtr % dimSizes(4)
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldOutPtr => fieldOutPtr % next
+ end do
+ end do
+
+ ! Determine size of receive list buffers.
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldOutPtr => fieldOut
+ do while(associated(fieldOutPtr))
+ exchListPtr => fieldOutPtr % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ nAdded = max(nAdded, maxval(exchListPtr % srcList) * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3) * fieldOutPtr % dimSizes(4))
+ end if
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldOutPtr => fieldOutPtr % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+ commListPtr % nList = nAdded
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! Allocate buffers for recieves, and initiate mpi_irecv calls.
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ allocate(commListPtr % rbuffer(commListPtr % nList))
+ nullify(commListPtr % ibuffer)
+ call MPI_Irecv(commListPtr % rbuffer, commListPtr % nList, MPI_realKIND, commListPtr % procID, commListPtr % procID, dminfo % comm, commListPtr % reqID, mpi_ierr)
+ commListPtr => commListPtr % next
+ end do
+
+ ! Setup send lists, and determine the size of their buffers.
+ do iHalo = 1, nHaloLayers
+ fieldInPtr => fieldIn
+ do while(associated(fieldInPtr))
+ exchListPtr => fieldInPtr % sendList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ comm_list_found = .false.
+
+ ! Search for an already created commList to this processor.
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ if(commListPtr % procID == exchListPtr % endPointID) then
+ commListPtr % nList = commListPtr % nList + exchListPtr % nList * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) * fieldOutPtr % dimSizes(3) * fieldInPtr % dimSizes(4)
+ comm_list_found = .true.
+ exit
+ end if
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! If no comm list exists, create a new one.
+ if(.not. comm_list_found) then
+ if(.not.associated(sendList)) then
+ allocate(sendList)
+ nullify(sendList % next)
+ commListPtr => sendList
+ else
+ commListPtr => sendList
+ commListPtr2 => commListPtr % next
+ do while(associated(commListPtr2))
+ commListPtr => commListPtr % next
+ commListPtr2 => commListPtr % next
+ end do
+
+ allocate(commListPtr % next)
+ commListPtr => commListPtr % next
+ nullify(commListPtr % next)
+ end if
+ commListPtr % procID = exchListPtr % endPointID
+ commListPtr % nList = exchListPtr % nList * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) * fieldOutPtr % dimSizes(3) * fieldInPtr % dimSizes(4)
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldInPtr => fieldInPtr % next
+ end do
+ end do
+
+ ! Allocate sendLists, copy data into buffer, and initiate mpi_isends
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ allocate(commListPtr % rbuffer(commListPtr % nList))
+ nullify(commListPtr % ibuffer)
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldInPtr => fieldIn
+ do while(associated(fieldInPtr))
+ exchListPtr => fieldInPtr % sendList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ do i = 1, exchListPtr % nList
+ do j = 1, fieldInPtr % dimSizes(4)
+ do k = 1, fieldInPtr % dimSizes(3)
+ do l = 1, fieldInPtr % dimSizes(2)
+ do m = 1, fieldInPtr % dimSizes(1)
+ iBuffer = (exchListPtr % destList(i)-1) * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) * fieldInPtr % dimSizes(3) * fieldInPtr % dimSizes(4) &
+ + (j-1) * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) * fieldInPtr % dimSizes(3) &
+ + (k-1) * fieldInPtr % dimSizes(1) * fieldInPtr % dimSizes(2) &
+ + (l-1) * fieldInPtr % dimSizes(1) + m + bufferOffset
+ commListPtr % rbuffer(iBuffer) = fieldInPtr % array(m, l, k, j, exchListPtr % srcList(i))
+ nAdded = nAdded + 1
+ end do
+ end do
+ end do
+ end do
+ end do
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldInPtr => fieldInPtr % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+
+ call MPI_Isend(commListPtr % rbuffer, commListPtr % nlist, MPI_realKIND, &
+ commListPtr % procID, dminfo % my_proc_id, dminfo % comm, commListPtr % reqID, mpi_ierr)
+
+ commListPtr => commListPtr % next
+ end do
+
+#endif
+
+ ! Handle Local Copies. Only local copies if no MPI
+ do iHalo = 1, nHaloLayers
+ fieldInPtr => fieldIn
+ do while(associated(fieldInPtr))
+ exchListPtr => fieldInPtr % copyList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ fieldOutPtr => fieldOut
+ do while(associated(fieldOutPtr))
+ if(exchListPtr % endPointID == fieldOutPtr % block % localBlockID) then
+ do i = 1, exchListPtr % nList
+ fieldOutPtr % array(:, :, :, :, exchListPtr % destList(i)) = fieldInPtr % array(:, :, :, :, exchListPtr % srcList(i))
+ end do
+ end if
+ fieldOutPtr => fieldOutPtr % next
+ end do
+
+ exchListPtr => exchListPtr % next
+ end do
+ fieldInPtr => fieldInPtr % next
+ end do
+ end do
+
+#ifdef _MPI
+ ! Wait for MPI_Irecv's to finish, and unpack data.
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldOutPtr => fieldOut
+ do while(associated(fieldOutPtr))
+ exchListPtr => fieldOutPtr % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ do i = 1, exchListPtr % nList
+ do j = 1, fieldOutPtr % dimSizes(4)
+ do k = 1, fieldOutPtr % dimSizes(3)
+ do l = 1, fieldOutPtr % dimSizes(2)
+ do m = 1, fieldOutPtr % dimSizes(1)
+ iBuffer = (exchListPtr % srcList(i)-1) * fieldOutPtr % dimSizes(4) * fieldOutPtr % dimSizes(3) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(1) &
+ + (j-1) * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3) &
+ + (k-1) * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) &
+ + (l-1) * fieldOutPtr % dimSizes(1) + m + bufferOffset
+ fieldOutPtr % array(m, l, k, j, exchListPtr % destList(i)) = commListPtr % rbuffer(iBuffer)
+ end do
+ end do
+ end do
+ end do
+ end do
+ nAdded = max(nAdded, maxval(exchListPtr % srcList) * fieldOutPtr % dimSizes(1) * fieldOutPtr % dimSizes(2) * fieldOutPtr % dimSizes(3) * fieldOutPtr % dimSizes(4))
+ end if
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldOutPtr => fieldOutPtr % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! Wait for MPI_Isend's to finish.
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ commListPtr => commListPtr % next
+ end do
+
+ ! Destroy commLists.
+ call mpas_dmpar_destroy_communication_list(sendList)
+ call mpas_dmpar_destroy_communication_list(recvList)
+#endif
+
+ deallocate(haloLayers)
+
+ end subroutine mpas_dmpar_alltoall_field5d_real!}}}
+
+
+
subroutine mpas_dmpar_exch_halo_field1d_integer(field, haloLayersIn)!{{{
implicit none
@@ -4507,6 +5122,602 @@
end subroutine mpas_dmpar_exch_halo_field3d_real!}}}
+ subroutine mpas_dmpar_exch_halo_field4d_real(field, haloLayersIn)!{{{
+
+ implicit none
+
+ type (field4dReal), pointer :: field
+ integer, dimension(:), intent(in), optional :: haloLayersIn
+
+ type (dm_info), pointer :: dminfo
+ type (field4dReal), pointer :: fieldCursor, fieldCursor2
+ type (mpas_exchange_list), pointer :: exchListPtr
+ type (mpas_communication_list), pointer :: sendList, recvList, commListPtr, commListPtr2
+ integer :: mpi_ierr
+ integer :: nHaloLayers, iHalo, i, j, k, l
+ integer :: bufferOffset, nAdded
+ integer, dimension(:), pointer :: haloLayers
+
+ logical :: comm_list_found
+
+ do i = 1, 4
+ if(field % dimSizes(i) <= 0) then
+ return
+ end if
+ end do
+
+ dminfo => field % block % domain % dminfo
+
+ if(present(haloLayersIn)) then
+ nHaloLayers = size(haloLayersIn)
+ allocate(haloLayers(nHaloLayers))
+ do iHalo = 1, nHaloLayers
+ haloLayers(iHalo) = haloLayersIn(iHalo)
+ end do
+ else
+ nHaloLayers = size(field % sendList % halos)
+ allocate(haloLayers(nHaloLayers))
+ do iHalo = 1, nHaloLayers
+ haloLayers(iHalo) = iHalo
+ end do
+ end if
+
+#ifdef _MPI
+ ! Allocate communication lists, and setup dead header nodes
+ allocate(sendList)
+ nullify(sendList % next)
+ sendList % procID = -1
+ sendList % nList = 0
+
+ allocate(recvList)
+ nullify(recvList % next)
+ recvList % procID = -1
+ recvList % nList = 0
+
+ dminfo = field % block % domain % dminfo
+
+ ! Determine size of buffers for communication lists
+ fieldCursor => field
+ do while(associated(fieldCursor))
+
+ ! Need to aggregate across halo layers
+ do iHalo = 1, nHaloLayers
+
+ ! Determine size from send lists
+ exchListPtr => fieldCursor % sendList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ comm_list_found = .false.
+
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ if(commListPtr % procID == exchListPtr % endPointId) then
+ comm_list_found = .true.
+ commListPtr % nList = commListPtr % nList + exchListPtr % nList * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3)
+ exit
+ end if
+
+ commListPtr => commListPtr % next
+ end do
+
+ if(.not. comm_list_found) then
+ commListPtr => sendList
+ commListPtr2 => commListPtr % next
+ do while(associated(commListPtr2))
+ commListPtr => commListPtr % next
+ commListPtr2 => commListPtr % next
+ end do
+
+ allocate(commListPtr % next)
+ commListPtr => commListPtr % next
+ nullify(commListPtr % next)
+ commListPtr % procID = exchListPtr % endPointID
+ commListPtr % nList = exchListPtr % nList * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3)
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ ! Setup recv lists
+ exchListPtr => fieldCursor % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ comm_list_found = .false.
+
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ if(commListPtr % procID == exchListPtr % endPointId) then
+ comm_list_found = .true.
+ exit
+ end if
+
+ commListPtr => commListPtr % next
+ end do
+
+ if(.not. comm_list_found) then
+ commListPtr => recvList
+ commListPtr2 => commListPtr % next
+ do while(associated(commListPtr2))
+ commListPtr => commListPtr % next
+ commListPtr2 => commListPtr % next
+ end do
+
+ allocate(commListPtr % next)
+ commListPtr => commListPtr % next
+ nullify(commListPtr % next)
+ commListPtr % procID = exchListPtr % endPointID
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+
+ ! Remove the dead head pointer on send and recv list
+ commListPtr => sendList
+ sendList => sendList % next
+ deallocate(commListPtr)
+
+ commListPtr => recvList
+ recvList => recvList % next
+ deallocate(commListPtr)
+
+ ! Determine size of recv lists
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldCursor => field
+ do while(associated(fieldCursor))
+ exchListPtr => fieldCursor % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ nAdded = max(nAdded, maxval(exchListPtr % srcList) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3))
+ end if
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+ commListPtr % nList = bufferOffset
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! Allocate space in recv lists, and initiate mpi_irecv calls
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ allocate(commListPtr % rbuffer(commListPtr % nList))
+ nullify(commListPtr % ibuffer)
+ call MPI_Irecv(commListPtr % rbuffer, commListPtr % nList, MPI_REALKIND, commListPtr % procID, commListPtr % procID, dminfo % comm, commListPtr % reqID, mpi_ierr)
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! Allocate space in send lists, copy data into buffer, and initiate mpi_isend calls
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ allocate(commListPtr % rbuffer(commListPtr % nList))
+ nullify(commListPtr % ibuffer)
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldCursor => field
+ do while(associated(fieldCursor))
+ exchListPtr => fieldCursor % sendList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ do i = 1, exchListPtr % nList
+ do j = 1, fieldCursor % dimSizes(3)
+ do k = 1, fieldCursor % dimSizes(2)
+ do l = 1, fieldCursor % dimSizes(1)
+ commListPtr % rbuffer((exchListPtr % destList(i)-1) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3) &
+ + (j-1) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) &
+ + (k-1) * fieldCursor % dimSizes(1) + l + bufferOffset) &
+ = fieldCursor % array(l, k, j, exchListPtr % srcList(i))
+ nAdded = nAdded + 1
+ end do
+ end do
+ end do
+ end do
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+
+ call MPI_Isend(commListPtr % rbuffer, commListPtr % nList, MPI_REALKIND, commListPtr % procID, dminfo % my_proc_id, dminfo % comm, commListPtr % reqID, mpi_ierr)
+ commListPtr => commListPtr % next
+ end do
+#endif
+
+ ! Handle local copy. If MPI is off, then only local copies are performed.
+ fieldCursor => field
+ do while(associated(fieldCursor))
+ do iHalo = 1, nHaloLayers
+ exchListPtr => fieldCursor % copyList % halos(haloLayers(iHalo)) % exchList
+
+ do while(associated(exchListPtr))
+ fieldCursor2 => field
+ do while(associated(fieldCursor2))
+ if(exchListPtr % endPointID == fieldCursor2 % block % localBlockID) then
+ do i = 1, exchListPtr % nList
+ fieldCursor2 % array(:, :, :, exchListPtr % destList(i)) = fieldCursor % array(:, :, :, exchListPtr % srcList(i))
+ end do
+ end if
+
+ fieldCursor2 => fieldCursor2 % next
+ end do
+
+ exchListPtr => exchListPtr % next
+ end do
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+
+#ifdef _MPI
+
+ ! Wait for mpi_irecv to finish, and unpack data from buffer
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldCursor => field
+ do while(associated(fieldCursor))
+ exchListPtr => fieldCursor % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ do i = 1, exchListPtr % nList
+ do j = 1, fieldCursor % dimSizes(3)
+ do k = 1, fieldCursor % dimSizes(2)
+ do l = 1, fieldCursor % dimSizes(1)
+ fieldCursor % array(l, k, j, exchListPtr % destList(i)) = commListPtr % rbuffer((exchListPtr % srcList(i)-1)&
+ *fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) *fieldCursor % dimSizes(3)&
+ + (j-1)*fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) &
+ + (k-1)*fieldCursor % dimSizes(1) + l + bufferOffset)
+ end do
+ end do
+ end do
+ end do
+ nAdded = max(nAdded, maxval(exchListPtr % srcList) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3))
+ end if
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+ commListPtr => commListPtr % next
+ end do
+
+ ! wait for mpi_isend to finish.
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ commListPtr => commListPtr % next
+ end do
+
+ ! Destroy commLists.
+ call mpas_dmpar_destroy_communication_list(sendList)
+ call mpas_dmpar_destroy_communication_list(recvList)
+#endif
+
+ deallocate(haloLayers)
+
+ end subroutine mpas_dmpar_exch_halo_field4d_real!}}}
+
+ subroutine mpas_dmpar_exch_halo_field5d_real(field, haloLayersIn)!{{{
+
+ implicit none
+
+ type (field5dReal), pointer :: field
+ integer, dimension(:), intent(in), optional :: haloLayersIn
+
+ type (dm_info), pointer :: dminfo
+ type (field5dReal), pointer :: fieldCursor, fieldCursor2
+ type (mpas_exchange_list), pointer :: exchListPtr
+ type (mpas_communication_list), pointer :: sendList, recvList, commListPtr, commListPtr2
+ integer :: mpi_ierr
+ integer :: nHaloLayers, iHalo, i, j, k, l, m
+ integer :: bufferOffset, nAdded
+ integer, dimension(:), pointer :: haloLayers
+
+ logical :: comm_list_found
+
+ do i = 1, 5
+ if(field % dimSizes(i) <= 0) then
+ return
+ end if
+ end do
+
+ dminfo => field % block % domain % dminfo
+
+ if(present(haloLayersIn)) then
+ nHaloLayers = size(haloLayersIn)
+ allocate(haloLayers(nHaloLayers))
+ do iHalo = 1, nHaloLayers
+ haloLayers(iHalo) = haloLayersIn(iHalo)
+ end do
+ else
+ nHaloLayers = size(field % sendList % halos)
+ allocate(haloLayers(nHaloLayers))
+ do iHalo = 1, nHaloLayers
+ haloLayers(iHalo) = iHalo
+ end do
+ end if
+
+#ifdef _MPI
+ ! Allocate communication lists, and setup dead header nodes
+ allocate(sendList)
+ nullify(sendList % next)
+ sendList % procID = -1
+ sendList % nList = 0
+
+ allocate(recvList)
+ nullify(recvList % next)
+ recvList % procID = -1
+ recvList % nList = 0
+
+ dminfo = field % block % domain % dminfo
+
+ ! Determine size of buffers for communication lists
+ fieldCursor => field
+ do while(associated(fieldCursor))
+
+ ! Need to aggregate across halo layers
+ do iHalo = 1, nHaloLayers
+
+ ! Determine size from send lists
+ exchListPtr => fieldCursor % sendList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ comm_list_found = .false.
+
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ if(commListPtr % procID == exchListPtr % endPointId) then
+ comm_list_found = .true.
+ commListPtr % nList = commListPtr % nList + exchListPtr % nList * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3) * fieldCursor % dimSizes(4)
+ exit
+ end if
+
+ commListPtr => commListPtr % next
+ end do
+
+ if(.not. comm_list_found) then
+ commListPtr => sendList
+ commListPtr2 => commListPtr % next
+ do while(associated(commListPtr2))
+ commListPtr => commListPtr % next
+ commListPtr2 => commListPtr % next
+ end do
+
+ allocate(commListPtr % next)
+ commListPtr => commListPtr % next
+ nullify(commListPtr % next)
+ commListPtr % procID = exchListPtr % endPointID
+ commListPtr % nList = exchListPtr % nList * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3) * fieldCursor % dimSizes(4)
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ ! Setup recv lists
+ exchListPtr => fieldCursor % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ comm_list_found = .false.
+
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ if(commListPtr % procID == exchListPtr % endPointId) then
+ comm_list_found = .true.
+ exit
+ end if
+
+ commListPtr => commListPtr % next
+ end do
+
+ if(.not. comm_list_found) then
+ commListPtr => recvList
+ commListPtr2 => commListPtr % next
+ do while(associated(commListPtr2))
+ commListPtr => commListPtr % next
+ commListPtr2 => commListPtr % next
+ end do
+
+ allocate(commListPtr % next)
+ commListPtr => commListPtr % next
+ nullify(commListPtr % next)
+ commListPtr % procID = exchListPtr % endPointID
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+
+ ! Remove the dead head pointer on send and recv list
+ commListPtr => sendList
+ sendList => sendList % next
+ deallocate(commListPtr)
+
+ commListPtr => recvList
+ recvList => recvList % next
+ deallocate(commListPtr)
+
+ ! Determine size of recv lists
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldCursor => field
+ do while(associated(fieldCursor))
+ exchListPtr => fieldCursor % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ nAdded = max(nAdded, maxval(exchListPtr % srcList) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3) * fieldCursor % dimSizes(4))
+ end if
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+ commListPtr % nList = bufferOffset
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! Allocate space in recv lists, and initiate mpi_irecv calls
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ allocate(commListPtr % rbuffer(commListPtr % nList))
+ nullify(commListPtr % ibuffer)
+ call MPI_Irecv(commListPtr % rbuffer, commListPtr % nList, MPI_REALKIND, commListPtr % procID, commListPtr % procID, dminfo % comm, commListPtr % reqID, mpi_ierr)
+
+ commListPtr => commListPtr % next
+ end do
+
+ ! Allocate space in send lists, copy data into buffer, and initiate mpi_isend calls
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ allocate(commListPtr % rbuffer(commListPtr % nList))
+ nullify(commListPtr % ibuffer)
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldCursor => field
+ do while(associated(fieldCursor))
+ exchListPtr => fieldCursor % sendList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ do i = 1, exchListPtr % nList
+ do j = 1, fieldCursor % dimSizes(4)
+ do k = 1, fieldCursor % dimSizes(3)
+ do l = 1, fieldCursor % dimSizes(2)
+ do m = 1, fieldCursor % dimSizes(1)
+ commListPtr % rbuffer((exchListPtr % destList(i)-1) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3) * fieldCursor % dimSizes(4) &
+ + (j-1) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3) &
+ + (k-1) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) &
+ + (l-1) * fieldCursor % dimSizes(1) + m + bufferOffset) &
+ = fieldCursor % array(m, l, k, j, exchListPtr % srcList(i))
+ nAdded = nAdded + 1
+ end do
+ end do
+ end do
+ end do
+ end do
+ end if
+
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+
+ call MPI_Isend(commListPtr % rbuffer, commListPtr % nList, MPI_REALKIND, commListPtr % procID, dminfo % my_proc_id, dminfo % comm, commListPtr % reqID, mpi_ierr)
+ commListPtr => commListPtr % next
+ end do
+#endif
+
+ ! Handle local copy. If MPI is off, then only local copies are performed.
+ fieldCursor => field
+ do while(associated(fieldCursor))
+ do iHalo = 1, nHaloLayers
+ exchListPtr => fieldCursor % copyList % halos(haloLayers(iHalo)) % exchList
+
+ do while(associated(exchListPtr))
+ fieldCursor2 => field
+ do while(associated(fieldCursor2))
+ if(exchListPtr % endPointID == fieldCursor2 % block % localBlockID) then
+ do i = 1, exchListPtr % nList
+ fieldCursor2 % array(:, :, :, :, exchListPtr % destList(i)) = fieldCursor % array(:, :, :, :, exchListPtr % srcList(i))
+ end do
+ end if
+
+ fieldCursor2 => fieldCursor2 % next
+ end do
+
+ exchListPtr => exchListPtr % next
+ end do
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+
+#ifdef _MPI
+
+ ! Wait for mpi_irecv to finish, and unpack data from buffer
+ commListPtr => recvList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ bufferOffset = 0
+ do iHalo = 1, nHaloLayers
+ nAdded = 0
+ fieldCursor => field
+ do while(associated(fieldCursor))
+ exchListPtr => fieldCursor % recvList % halos(haloLayers(iHalo)) % exchList
+ do while(associated(exchListPtr))
+ if(exchListPtr % endPointID == commListPtr % procID) then
+ do i = 1, exchListPtr % nList
+ do j = 1, fieldCursor % dimSizes(4)
+ do k = 1, fieldCursor % dimSizes(3)
+ do l = 1, fieldCursor % dimSizes(2)
+ do m = 1, fieldCursor % dimSizes(1)
+ fieldCursor % array(m, l, k, j, exchListPtr % destList(i)) = commListPtr % rbuffer((exchListPtr % srcList(i)-1)&
+ *fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) *fieldCursor % dimSizes(3) * fieldCursor % dimSizes(4)&
+ + (j-1)*fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3) &
+ + (k-1)*fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) &
+ + (l-1)*fieldCursor % dimSizes(1) + m + bufferOffset)
+ end do
+ end do
+ end do
+ end do
+ end do
+ nAdded = max(nAdded, maxval(exchListPtr % srcList) * fieldCursor % dimSizes(1) * fieldCursor % dimSizes(2) * fieldCursor % dimSizes(3) * fieldCursor % dimSizes(4))
+ end if
+ exchListPtr => exchListPtr % next
+ end do
+
+ fieldCursor => fieldCursor % next
+ end do
+ bufferOffset = bufferOffset + nAdded
+ end do
+ commListPtr => commListPtr % next
+ end do
+
+ ! wait for mpi_isend to finish.
+ commListPtr => sendList
+ do while(associated(commListPtr))
+ call MPI_Wait(commListPtr % reqID, MPI_STATUS_IGNORE, mpi_ierr)
+ commListPtr => commListPtr % next
+ end do
+
+ ! Destroy commLists.
+ call mpas_dmpar_destroy_communication_list(sendList)
+ call mpas_dmpar_destroy_communication_list(recvList)
+#endif
+
+ deallocate(haloLayers)
+
+ end subroutine mpas_dmpar_exch_halo_field5d_real!}}}
+
subroutine mpas_dmpar_init_mulithalo_exchange_list(exchList, nHalos)!{{{
type (mpas_multihalo_exchange_list), pointer :: exchList
integer, intent(in) :: nHalos
@@ -4668,4 +5879,30 @@
end if
end subroutine mpas_dmpar_copy_field3d_real!}}}
+ subroutine mpas_dmpar_copy_field4d_real(field)!{{{
+ type (field4dReal), pointer :: field
+ type (field4dReal), pointer :: fieldCursor
+
+ if(associated(field % next)) then
+ fieldCursor => field % next
+ do while(associated(fieldCursor))
+ fieldCursor % array = field % array
+ fieldCursor => fieldCursor % next
+ end do
+ end if
+ end subroutine mpas_dmpar_copy_field4d_real!}}}
+
+ subroutine mpas_dmpar_copy_field5d_real(field)!{{{
+ type (field5dReal), pointer :: field
+ type (field5dReal), pointer :: fieldCursor
+
+ if(associated(field % next)) then
+ fieldCursor => field % next
+ do while(associated(fieldCursor))
+ fieldCursor % array = field % array
+ fieldCursor => fieldCursor % next
+ end do
+ end if
+ end subroutine mpas_dmpar_copy_field5d_real!}}}
+
end module mpas_dmpar
Modified: branches/ocean_projects/cvmix_implementation/src/framework/mpas_grid_types.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/framework/mpas_grid_types.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/framework/mpas_grid_types.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -20,8 +20,66 @@
logical :: output
end type io_info
+ ! Derived type for storing fields
+ type field5DReal
+
+ ! Back-pointer to the containing block
+ type (block_type), pointer :: block
+ ! Raw array holding field data on this block
+ real (kind=RKIND), dimension(:,:,:,:,:), pointer :: array
+
+ ! Information used by the I/O layer
+ type (io_info), pointer :: ioinfo ! to be removed later
+ character (len=StrKIND) :: fieldName
+ character (len=StrKIND), dimension(:), pointer :: constituentNames => null()
+ character (len=StrKIND), dimension(5) :: dimNames
+ integer, dimension(5) :: dimSizes
+ logical :: hasTimeDimension
+ logical :: isSuperArray
+ type (att_list_type), pointer :: attList => null()
+
+ ! Pointers to the prev and next blocks for this field on this task
+ type (field5DReal), pointer :: prev, next
+
+ ! Halo communication lists
+ type (mpas_multihalo_exchange_list), pointer :: sendList
+ type (mpas_multihalo_exchange_list), pointer :: recvList
+ type (mpas_multihalo_exchange_list), pointer :: copyList
+ end type field5DReal
+
+
! Derived type for storing fields
+ type field4DReal
+
+ ! Back-pointer to the containing block
+ type (block_type), pointer :: block
+
+ ! Raw array holding field data on this block
+ real (kind=RKIND), dimension(:,:,:,:), pointer :: array
+
+ ! Information used by the I/O layer
+ type (io_info), pointer :: ioinfo ! to be removed later
+ character (len=StrKIND) :: fieldName
+ character (len=StrKIND), dimension(:), pointer :: constituentNames => null()
+ character (len=StrKIND), dimension(4) :: dimNames
+ integer, dimension(4) :: dimSizes
+ logical :: hasTimeDimension
+ logical :: isSuperArray
+ type (att_list_type), pointer :: attList => null()
+
+ ! Pointers to the prev and next blocks for this field on this task
+ type (field4DReal), pointer :: prev, next
+
+ ! Halo communication lists
+ type (mpas_multihalo_exchange_list), pointer :: sendList
+ type (mpas_multihalo_exchange_list), pointer :: recvList
+ type (mpas_multihalo_exchange_list), pointer :: copyList
+ end type field4DReal
+
+
+
+ ! Derived type for storing fields
type field3DReal
! Back-pointer to the containing block
@@ -370,6 +428,8 @@
module procedure mpas_allocate_scratch_field1d_real
module procedure mpas_allocate_scratch_field2d_real
module procedure mpas_allocate_scratch_field3d_real
+ module procedure mpas_allocate_scratch_field4d_real
+ module procedure mpas_allocate_scratch_field5d_real
module procedure mpas_allocate_scratch_field1d_char
end interface
@@ -380,6 +440,8 @@
module procedure mpas_deallocate_scratch_field1d_real
module procedure mpas_deallocate_scratch_field2d_real
module procedure mpas_deallocate_scratch_field3d_real
+ module procedure mpas_deallocate_scratch_field4d_real
+ module procedure mpas_deallocate_scratch_field5d_real
module procedure mpas_deallocate_scratch_field1d_char
end interface
@@ -392,6 +454,8 @@
module procedure mpas_deallocate_field1d_real
module procedure mpas_deallocate_field2d_real
module procedure mpas_deallocate_field3d_real
+ module procedure mpas_deallocate_field4d_real
+ module procedure mpas_deallocate_field5d_real
module procedure mpas_deallocate_field0d_char
module procedure mpas_deallocate_field1d_char
end interface
@@ -632,6 +696,62 @@
end subroutine mpas_allocate_scratch_field3d_real!}}}
+ subroutine mpas_allocate_scratch_field4d_real(f, single_block_in)!{{{
+ type (field4dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field4dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1), f_cursor % dimSizes(2), f_cursor % dimSizes(3), f_cursor % dimSizes(4)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1), f % dimSizes(2), f % dimSizes(3), f % dimSizes(4)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field4d_real!}}}
+
+ subroutine mpas_allocate_scratch_field5d_real(f, single_block_in)!{{{
+ type (field5dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field5dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not. single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(.not.associated(f_cursor % array)) then
+ allocate(f_cursor % array(f_cursor % dimSizes(1), f_cursor % dimSizes(2), f_cursor % dimSizes(3), f_cursor % dimSizes(4), f_cursor % dimSizes(5)))
+ end if
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(.not.associated(f % array)) then
+ allocate(f % array(f % dimSizes(1), f % dimSizes(2), f % dimSizes(3), f % dimSizes(4), f % dimSizes(5)))
+ end if
+ end if
+
+ end subroutine mpas_allocate_scratch_field5d_real!}}}
+
subroutine mpas_allocate_scratch_field1d_char(f, single_block_in)!{{{
type (field1dChar), pointer :: f
logical, intent(in), optional :: single_block_in
@@ -834,6 +954,64 @@
end subroutine mpas_deallocate_scratch_field3d_real!}}}
+ subroutine mpas_deallocate_scratch_field4d_real(f, single_block_in)!{{{
+ type (field4dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field4dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field4d_real!}}}
+
+ subroutine mpas_deallocate_scratch_field5d_real(f, single_block_in)!{{{
+ type (field5dReal), pointer :: f
+ logical, intent(in), optional :: single_block_in
+ logical :: single_block
+ type (field5dReal), pointer :: f_cursor
+
+ if(present(single_block_in)) then
+ single_block = single_block_in
+ else
+ single_block = .false.
+ end if
+
+ if(.not.single_block) then
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ f_cursor => f_cursor % next
+ end do
+ else
+ if(associated(f % array)) then
+ deallocate(f % array)
+ end if
+ end if
+
+ end subroutine mpas_deallocate_scratch_field5d_real!}}}
+
subroutine mpas_deallocate_scratch_field1d_char(f, single_block_in)!{{{
type (field1dChar), pointer :: f
logical, intent(in), optional :: single_block_in
@@ -1073,6 +1251,60 @@
end subroutine mpas_deallocate_field3d_real!}}}
+ subroutine mpas_deallocate_field4d_real(f)!{{{
+ type (field4dReal), pointer :: f
+ type (field4dReal), pointer :: f_cursor
+
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f % next)) then
+ f => f % next
+ else
+ nullify(f)
+ end if
+
+ if(associated(f_cursor % ioinfo)) then
+ deallocate(f_cursor % ioinfo)
+ end if
+
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ deallocate(f_cursor)
+
+ f_cursor => f
+ end do
+
+ end subroutine mpas_deallocate_field4d_real!}}}
+
+ subroutine mpas_deallocate_field5d_real(f)!{{{
+ type (field5dReal), pointer :: f
+ type (field5dReal), pointer :: f_cursor
+
+ f_cursor => f
+ do while(associated(f_cursor))
+ if(associated(f % next)) then
+ f => f % next
+ else
+ nullify(f)
+ end if
+
+ if(associated(f_cursor % ioinfo)) then
+ deallocate(f_cursor % ioinfo)
+ end if
+
+ if(associated(f_cursor % array)) then
+ deallocate(f_cursor % array)
+ end if
+
+ deallocate(f_cursor)
+
+ f_cursor => f
+ end do
+
+ end subroutine mpas_deallocate_field5d_real!}}}
+
subroutine mpas_deallocate_field0d_char(f)!{{{
type (field0dChar), pointer :: f
type (field0dChar), pointer :: f_cursor
Modified: branches/ocean_projects/cvmix_implementation/src/framework/mpas_io.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/framework/mpas_io.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/framework/mpas_io.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -77,6 +77,7 @@
module procedure MPAS_io_get_var_real2d
module procedure MPAS_io_get_var_real3d
module procedure MPAS_io_get_var_real4d
+ module procedure MPAS_io_get_var_real5d
module procedure MPAS_io_get_var_char0d
end interface MPAS_io_get_var
@@ -91,6 +92,7 @@
module procedure MPAS_io_put_var_real2d
module procedure MPAS_io_put_var_real3d
module procedure MPAS_io_put_var_real4d
+ module procedure MPAS_io_put_var_real5d
module procedure MPAS_io_put_var_char0d
end interface MPAS_io_put_var
@@ -1146,7 +1148,7 @@
subroutine MPAS_io_get_var_generic(handle, fieldname, intVal, intArray1d, intArray2d, intArray3d, intArray4d, &
- realVal, realArray1d, realArray2d, realArray3d, realArray4d, &
+ realVal, realArray1d, realArray2d, realArray3d, realArray4d, realArray5d, &
charVal, ierr)
implicit none
@@ -1163,6 +1165,7 @@
real (kind=RKIND), dimension(:,:), intent(out), optional :: realArray2d
real (kind=RKIND), dimension(:,:,:), intent(out), optional :: realArray3d
real (kind=RKIND), dimension(:,:,:,:), intent(out), optional :: realArray4d
+ real (kind=RKIND), dimension(:,:,:,:,:), intent(out), optional :: realArray5d
character (len=*), intent(out), optional :: charVal
integer, intent(out), optional :: ierr
@@ -1263,6 +1266,10 @@
! write (0,*) ' value is real4'
call PIO_read_darray(handle % pio_file, field_cursor % fieldhandle % field_desc, field_cursor % fieldhandle % decomp % pio_iodesc, &
realArray4d, pio_ierr)
+ else if (present(realArray5d)) then
+! write (0,*) ' value is real5'
+ call PIO_read_darray(handle % pio_file, field_cursor % fieldhandle % field_desc, field_cursor % fieldhandle % decomp % pio_iodesc, &
+ realArray5d, pio_ierr)
else if (present(intArray1d)) then
! write (0,*) ' value is int1'
call PIO_read_darray(handle % pio_file, field_cursor % fieldhandle % field_desc, field_cursor % fieldhandle % decomp % pio_iodesc, &
@@ -1492,6 +1499,26 @@
end subroutine MPAS_io_get_var_real4d
+ subroutine MPAS_io_get_var_real5d(handle, fieldname, array, ierr)
+
+ implicit none
+
+ type (MPAS_IO_Handle_type), intent(inout) :: handle
+ character (len=*), intent(in) :: fieldname
+ real (kind=RKIND), dimension(:,:,:,:,:), intent(out) :: array
+ integer, intent(out), optional :: ierr
+
+ integer :: pio_ierr
+ type (fieldlist_type), pointer :: field_cursor
+
+! write(0,*) 'Called MPAS_io_get_var_real5d()'
+ if (present(ierr)) ierr = MPAS_IO_NOERR
+
+ call MPAS_io_get_var_generic(handle, fieldname, realArray5d=array, ierr=ierr)
+
+ end subroutine MPAS_io_get_var_real5d
+
+
subroutine MPAS_io_get_var_char0d(handle, fieldname, val, ierr)
implicit none
@@ -1513,7 +1540,7 @@
subroutine MPAS_io_put_var_generic(handle, fieldname, intVal, intArray1d, intArray2d, intArray3d, intArray4d, &
- realVal, realArray1d, realArray2d, realArray3d, realArray4d, &
+ realVal, realArray1d, realArray2d, realArray3d, realArray4d, realArray5d, &
charVal, ierr)
implicit none
@@ -1530,6 +1557,7 @@
real (kind=RKIND), dimension(:,:), intent(in), optional :: realArray2d
real (kind=RKIND), dimension(:,:,:), intent(in), optional :: realArray3d
real (kind=RKIND), dimension(:,:,:,:), intent(in), optional :: realArray4d
+ real (kind=RKIND), dimension(:,:,:,:,:), intent(in), optional :: realArray5d
character (len=*), intent(in), optional :: charVal
integer, intent(out), optional :: ierr
@@ -1629,6 +1657,9 @@
else if (present(realArray4d)) then
call PIO_write_darray(handle % pio_file, field_cursor % fieldhandle % field_desc, field_cursor % fieldhandle % decomp % pio_iodesc, &
realArray4d, pio_ierr)
+ else if (present(realArray5d)) then
+ call PIO_write_darray(handle % pio_file, field_cursor % fieldhandle % field_desc, field_cursor % fieldhandle % decomp % pio_iodesc, &
+ realArray5d, pio_ierr)
else if (present(intArray1d)) then
call PIO_write_darray(handle % pio_file, field_cursor % fieldhandle % field_desc, field_cursor % fieldhandle % decomp % pio_iodesc, &
intArray1d, pio_ierr)
@@ -1852,6 +1883,26 @@
end subroutine MPAS_io_put_var_real4d
+ subroutine MPAS_io_put_var_real5d(handle, fieldname, array, ierr)
+
+ implicit none
+
+ type (MPAS_IO_Handle_type), intent(inout) :: handle
+ character (len=*), intent(in) :: fieldname
+ real (kind=RKIND), dimension(:,:,:,:,:), intent(in) :: array
+ integer, intent(out), optional :: ierr
+
+ integer :: pio_ierr
+ type (fieldlist_type), pointer :: field_cursor
+
+! write(0,*) 'Called MPAS_io_put_var_real5d()'
+ if (present(ierr)) ierr = MPAS_IO_NOERR
+
+ call MPAS_io_put_var_generic(handle, fieldname, realArray5d=array, ierr=ierr)
+
+ end subroutine MPAS_io_put_var_real5d
+
+
subroutine MPAS_io_put_var_char0d(handle, fieldname, val, ierr)
implicit none
Modified: branches/ocean_projects/cvmix_implementation/src/framework/mpas_io_streams.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/framework/mpas_io_streams.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/framework/mpas_io_streams.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -19,6 +19,8 @@
type (field1dReal), pointer :: real1dField => null()
type (field2dReal), pointer :: real2dField => null()
type (field3dReal), pointer :: real3dField => null()
+ type (field4dReal), pointer :: real4dField => null()
+ type (field5dReal), pointer :: real5dField => null()
type (field0dChar), pointer :: char0dField => null()
type (field1dChar), pointer :: char1dField => null()
type (field_list_type), pointer :: next => null()
@@ -44,6 +46,8 @@
module procedure MPAS_streamAddField_1dReal
module procedure MPAS_streamAddField_2dReal
module procedure MPAS_streamAddField_3dReal
+ module procedure MPAS_streamAddField_4dReal
+ module procedure MPAS_streamAddField_5dReal
module procedure MPAS_streamAddField_0dChar
end interface MPAS_streamAddField
@@ -82,8 +86,10 @@
FIELD_1D_REAL = 6, &
FIELD_2D_REAL = 7, &
FIELD_3D_REAL = 8, &
- FIELD_0D_CHAR = 9, &
- FIELD_1D_CHAR = 10
+ FIELD_4D_REAL = 9, &
+ FIELD_5D_REAL = 10, &
+ FIELD_0D_CHAR = 11, &
+ FIELD_1D_CHAR = 12
private mergeArrays
@@ -996,6 +1002,208 @@
end subroutine MPAS_streamAddField_3dReal
+ subroutine MPAS_streamAddField_4dReal(stream, field, ierr)
+
+ implicit none
+
+ type (MPAS_Stream_type), intent(inout) :: stream
+ type (field4DReal), intent(in), target :: field
+ integer, intent(out), optional :: ierr
+
+ integer :: io_err
+ integer :: i
+ integer :: idim
+ integer :: totalDimSize, globalDimSize
+ logical :: isDecomposed
+ integer :: ndims
+ type (field4dReal), pointer :: field_ptr
+ character (len=StrKIND), dimension(5) :: dimNames
+ character (len=StrKIND), dimension(:), pointer :: dimNamesInq
+ integer, dimension(:), pointer :: indices
+ type (field_list_type), pointer :: field_list_cursor
+ type (field_list_type), pointer :: new_field_list_node
+ logical :: any_success
+ logical, dimension(:), pointer :: isAvailable
+
+ if (present(ierr)) ierr = MPAS_STREAM_NOERR
+
+ !
+ ! Sanity checks
+ !
+ if (.not. stream % isInitialized) then
+ if (present(ierr)) ierr = MPAS_STREAM_NOT_INITIALIZED
+ return
+ end if
+
+!write(0,*) '... Adding field '//trim(field % fieldName)//' to stream'
+
+ ndims = size(field % dimSizes)
+
+!write(0,*) '... field has ', ndims, ' dimensions'
+
+ !
+ ! Determine whether the field is decomposed, the indices that are owned by this task's blocks,
+ ! and the total number of outer-indices owned by this task
+ !
+#include "add_field_indices.inc"
+
+
+ any_success = .false.
+ if (field % isSuperArray) then
+!write(0,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^ we are adding a super-array'
+ allocate(isAvailable(size(field % constituentNames)))
+ isAvailable(:) = .false.
+ do i=1,size(field % constituentNames)
+ call MPAS_streamAddField_generic(stream, trim(field % constituentNames(i)), MPAS_IO_DOUBLE, field % dimNames(2:ndims), &
+ field % dimSizes(2:ndims), field % hasTimeDimension, isDecomposed, totalDimSize, globalDimSize, &
+ indices, io_err)
+ if (io_err == MPAS_STREAM_NOERR) then
+ isAvailable(i) = .true.
+ any_success = .true.
+ end if
+ end do
+ else
+ nullify(isAvailable)
+ call MPAS_streamAddField_generic(stream, trim(field % fieldName), MPAS_IO_DOUBLE, field % dimNames, field % dimSizes, &
+ field % hasTimeDimension, isDecomposed, totalDimSize, globalDimSize, indices, io_err)
+ if (io_err == MPAS_STREAM_NOERR) then
+ any_success = .true.
+ end if
+ end if
+
+ deallocate(indices)
+ if (.not. any_success) then
+ if (present(ierr)) ierr = MPAS_IO_ERR
+ return
+ end if
+
+ if (field % isSuperArray) then
+ do i=1,size(field % constituentNames)
+ call put_get_field_atts(stream % fileHandle, stream % ioDirection, trim(field % constituentNames(i)), field % attList)
+ end do
+ else
+ call put_get_field_atts(stream % fileHandle, stream % ioDirection, trim(field % fieldname), field % attList)
+ end if
+
+
+ !
+ ! Set field pointer and type in fieldList
+ !
+ new_field_list_node => stream % fieldList
+ do while (associated(new_field_list_node % next))
+ new_field_list_node => new_field_list_node % next
+ end do
+ new_field_list_node % field_type = FIELD_4D_REAL
+ new_field_list_node % real4dField => field
+ new_field_list_node % isAvailable => isAvailable
+
+!write(0,*) '... done adding field'
+!write(0,*) 'DEBUGGING : Finished adding 4d real field '//trim(field % fieldName)
+
+ end subroutine MPAS_streamAddField_4dReal
+
+
+ subroutine MPAS_streamAddField_5dReal(stream, field, ierr)
+
+ implicit none
+
+ type (MPAS_Stream_type), intent(inout) :: stream
+ type (field5DReal), intent(in), target :: field
+ integer, intent(out), optional :: ierr
+
+ integer :: io_err
+ integer :: i
+ integer :: idim
+ integer :: totalDimSize, globalDimSize
+ logical :: isDecomposed
+ integer :: ndims
+ type (field5dReal), pointer :: field_ptr
+ character (len=StrKIND), dimension(5) :: dimNames
+ character (len=StrKIND), dimension(:), pointer :: dimNamesInq
+ integer, dimension(:), pointer :: indices
+ type (field_list_type), pointer :: field_list_cursor
+ type (field_list_type), pointer :: new_field_list_node
+ logical :: any_success
+ logical, dimension(:), pointer :: isAvailable
+
+ if (present(ierr)) ierr = MPAS_STREAM_NOERR
+
+ !
+ ! Sanity checks
+ !
+ if (.not. stream % isInitialized) then
+ if (present(ierr)) ierr = MPAS_STREAM_NOT_INITIALIZED
+ return
+ end if
+
+!write(0,*) '... Adding field '//trim(field % fieldName)//' to stream'
+
+ ndims = size(field % dimSizes)
+
+!write(0,*) '... field has ', ndims, ' dimensions'
+
+ !
+ ! Determine whether the field is decomposed, the indices that are owned by this task's blocks,
+ ! and the total number of outer-indices owned by this task
+ !
+#include "add_field_indices.inc"
+
+
+ any_success = .false.
+ if (field % isSuperArray) then
+!write(0,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^ we are adding a super-array'
+ allocate(isAvailable(size(field % constituentNames)))
+ isAvailable(:) = .false.
+ do i=1,size(field % constituentNames)
+ call MPAS_streamAddField_generic(stream, trim(field % constituentNames(i)), MPAS_IO_DOUBLE, field % dimNames(2:ndims), &
+ field % dimSizes(2:ndims), field % hasTimeDimension, isDecomposed, totalDimSize, globalDimSize, &
+ indices, io_err)
+ if (io_err == MPAS_STREAM_NOERR) then
+ isAvailable(i) = .true.
+ any_success = .true.
+ end if
+ end do
+ else
+ nullify(isAvailable)
+ call MPAS_streamAddField_generic(stream, trim(field % fieldName), MPAS_IO_DOUBLE, field % dimNames, field % dimSizes, &
+ field % hasTimeDimension, isDecomposed, totalDimSize, globalDimSize, indices, io_err)
+ if (io_err == MPAS_STREAM_NOERR) then
+ any_success = .true.
+ end if
+ end if
+
+ deallocate(indices)
+ if (.not. any_success) then
+ if (present(ierr)) ierr = MPAS_IO_ERR
+ return
+ end if
+
+ if (field % isSuperArray) then
+ do i=1,size(field % constituentNames)
+ call put_get_field_atts(stream % fileHandle, stream % ioDirection, trim(field % constituentNames(i)), field % attList)
+ end do
+ else
+ call put_get_field_atts(stream % fileHandle, stream % ioDirection, trim(field % fieldname), field % attList)
+ end if
+
+
+ !
+ ! Set field pointer and type in fieldList
+ !
+ new_field_list_node => stream % fieldList
+ do while (associated(new_field_list_node % next))
+ new_field_list_node => new_field_list_node % next
+ end do
+ new_field_list_node % field_type = FIELD_5D_REAL
+ new_field_list_node % real5dField => field
+ new_field_list_node % isAvailable => isAvailable
+
+!write(0,*) '... done adding field'
+!write(0,*) 'DEBUGGING : Finished adding 3d real field '//trim(field % fieldName)
+
+ end subroutine MPAS_streamAddField_5dReal
+
+
subroutine MPAS_streamAddField_0dChar(stream, field, ierr)
implicit none
@@ -1313,6 +1521,8 @@
type (field1dReal), pointer :: field_1dreal_ptr
type (field2dReal), pointer :: field_2dreal_ptr
type (field3dReal), pointer :: field_3dreal_ptr
+ type (field4dReal), pointer :: field_4dreal_ptr
+ type (field5dReal), pointer :: field_5dreal_ptr
type (field0dChar), pointer :: field_0dchar_ptr
type (field1dChar), pointer :: field_1dchar_ptr
type (field_list_type), pointer :: field_cursor
@@ -1324,6 +1534,8 @@
real (kind=RKIND), dimension(:), pointer :: real1d_temp
real (kind=RKIND), dimension(:,:), pointer :: real2d_temp
real (kind=RKIND), dimension(:,:,:), pointer :: real3d_temp
+ real (kind=RKIND), dimension(:,:,:,:), pointer :: real4d_temp
+ real (kind=RKIND), dimension(:,:,:,:,:), pointer :: real5d_temp
if (present(ierr)) ierr = MPAS_STREAM_NOERR
@@ -1876,7 +2088,185 @@
else
deallocate(real3d_temp)
end if
+ else if (field_cursor % field_type == FIELD_4D_REAL) then
+!write(0,*) 'DEBUGGING : *************** '//trim(field_cursor % real3dField % fieldName)
+!write(0,*) 'DEBUGGING : reading a 4d real array'
+ if (field_cursor % real4dField % isSuperArray) then
+!write(0,*) 'DEBUGGING : reading a 4d real super-array'
+ ncons = size(field_cursor % real4dField % constituentNames)
+ allocate(real3d_temp(field_cursor % real4dField % dimSizes(2), &
+ field_cursor % real4dField % dimSizes(3), &
+ field_cursor % totalDimSize))
+ else
+ ncons = 1
+ allocate(real4d_temp(field_cursor % real4dField % dimSizes(1), &
+ field_cursor % real4dField % dimSizes(2), &
+ field_cursor % real4dField % dimSizes(3), &
+ field_cursor % totalDimSize))
+ end if
+
+ do j=1,ncons
+ if (field_cursor % real4dField % isSuperArray) then
+ if (.not. field_cursor % isAvailable(j)) cycle
+!write(0,*) 'DEBUGGING : calling get_var for a constitutent'
+ call MPAS_io_get_var(stream % fileHandle, field_cursor % real4dField % constituentNames(j), real3d_temp, io_err)
+ else
+ call MPAS_io_get_var(stream % fileHandle, field_cursor % real4dField % fieldName, real4d_temp, io_err)
+ end if
+ call MPAS_io_err_mesg(io_err, .false.)
+ if (io_err /= MPAS_IO_NOERR) then
+ if (present(ierr)) ierr = MPAS_IO_ERR
+ if (field_cursor % real4dField % isSuperArray) then
+ deallocate(real3d_temp)
+ else
+ deallocate(real4d_temp)
+ end if
+ return
+ end if
+
+ if (field_cursor % isDecomposed) then
+ ! Distribute field to multiple blocks
+ field_4dreal_ptr => field_cursor % real4dField
+ i = 1
+ do while (associated(field_4dreal_ptr))
+ if (trim(field_4dreal_ptr % dimNames(4)) == 'nCells') then
+ ownedSize = field_4dreal_ptr % block % mesh % nCellsSolve
+ else if (trim(field_4dreal_ptr % dimNames(4)) == 'nEdges') then
+ ownedSize = field_4dreal_ptr % block % mesh % nEdgesSolve
+ else if (trim(field_4dreal_ptr % dimNames(4)) == 'nVertices') then
+ ownedSize = field_4dreal_ptr % block % mesh % nVerticesSolve
+ else
+ ownedSize = field_4dreal_ptr % dimSizes(4)
+ end if
+
+ if (field_cursor % real4dField % isSuperArray) then
+!write(0,*) 'DEBUGGING : copying the temporary array'
+ field_4dreal_ptr % array(j, :,:,1:ownedSize) = real3d_temp(:,:,i:i+ownedSize-1)
+ else
+ field_4dreal_ptr % array(:,:,:,1:ownedSize) = real4d_temp(:,:,:,i:i+ownedSize-1)
+ end if
+ i = i + ownedSize
+ field_4dreal_ptr => field_4dreal_ptr % next
+ end do
+
+ else
+
+ if (field_cursor % real3dField % isSuperArray) then
+ call mpas_dmpar_bcast_reals(field_cursor % real4dField % block % domain % dminfo, size(real3d_temp), real3d_temp(:,1,1))
+ field_4dreal_ptr => field_cursor % real4dField
+ do while (associated(field_4dreal_ptr))
+ field_4dreal_ptr % array(j,:,:,:) = real3d_temp(:,:,:)
+ field_4dreal_ptr => field_4dreal_ptr % next
+ end do
+ else
+ call mpas_dmpar_bcast_reals(field_cursor % real4dField % block % domain % dminfo, size(real4d_temp), real4d_temp(:,1,1,1))
+ field_4dreal_ptr => field_cursor % real4dField
+ do while (associated(field_4dreal_ptr))
+ field_4dreal_ptr % array(:,:,:,:) = real4d_temp(:,:,:,:)
+ field_4dreal_ptr => field_4dreal_ptr % next
+ end do
+ end if
+ end if
+ end do
+
+ if (field_cursor % real4dField % isSuperArray) then
+ deallocate(real3d_temp)
+ else
+ deallocate(real4d_temp)
+ end if
+
+ else if (field_cursor % field_type == FIELD_5D_REAL) then
+
+!write(0,*) 'DEBUGGING : *************** '//trim(field_cursor % real3dField % fieldName)
+!write(0,*) 'DEBUGGING : reading a 4d real array'
+ if (field_cursor % real5dField % isSuperArray) then
+!write(0,*) 'DEBUGGING : reading a 4d real super-array'
+ ncons = size(field_cursor % real5dField % constituentNames)
+ allocate(real4d_temp(field_cursor % real5dField % dimSizes(2), &
+ field_cursor % real5dField % dimSizes(3), &
+ field_cursor % real5dField % dimSizes(4), &
+ field_cursor % totalDimSize))
+ else
+ ncons = 1
+ allocate(real5d_temp(field_cursor % real5dField % dimSizes(1), &
+ field_cursor % real5dField % dimSizes(2), &
+ field_cursor % real5dField % dimSizes(3), &
+ field_cursor % real5dField % dimSizes(4), &
+ field_cursor % totalDimSize))
+ end if
+
+ do j=1,ncons
+ if (field_cursor % real5dField % isSuperArray) then
+ if (.not. field_cursor % isAvailable(j)) cycle
+!write(0,*) 'DEBUGGING : calling get_var for a constitutent'
+ call MPAS_io_get_var(stream % fileHandle, field_cursor % real5dField % constituentNames(j), real4d_temp, io_err)
+ else
+ call MPAS_io_get_var(stream % fileHandle, field_cursor % real5dField % fieldName, real5d_temp, io_err)
+ end if
+ call MPAS_io_err_mesg(io_err, .false.)
+ if (io_err /= MPAS_IO_NOERR) then
+ if (present(ierr)) ierr = MPAS_IO_ERR
+ if (field_cursor % real5dField % isSuperArray) then
+ deallocate(real4d_temp)
+ else
+ deallocate(real5d_temp)
+ end if
+ return
+ end if
+
+ if (field_cursor % isDecomposed) then
+ ! Distribute field to multiple blocks
+ field_5dreal_ptr => field_cursor % real5dField
+ i = 1
+ do while (associated(field_5dreal_ptr))
+ if (trim(field_5dreal_ptr % dimNames(5)) == 'nCells') then
+ ownedSize = field_5dreal_ptr % block % mesh % nCellsSolve
+ else if (trim(field_5dreal_ptr % dimNames(5)) == 'nEdges') then
+ ownedSize = field_5dreal_ptr % block % mesh % nEdgesSolve
+ else if (trim(field_5dreal_ptr % dimNames(5)) == 'nVertices') then
+ ownedSize = field_5dreal_ptr % block % mesh % nVerticesSolve
+ else
+ ownedSize = field_5dreal_ptr % dimSizes(5)
+ end if
+
+ if (field_cursor % real5dField % isSuperArray) then
+!write(0,*) 'DEBUGGING : copying the temporary array'
+ field_5dreal_ptr % array(j,:,:,:,1:ownedSize) = real4d_temp(:,:,:,i:i+ownedSize-1)
+ else
+ field_5dreal_ptr % array(:,:,:,:,1:ownedSize) = real5d_temp(:,:,:,:,i:i+ownedSize-1)
+ end if
+ i = i + ownedSize
+ field_5dreal_ptr => field_5dreal_ptr % next
+ end do
+
+ else
+
+ if (field_cursor % real5dField % isSuperArray) then
+ call mpas_dmpar_bcast_reals(field_cursor % real5dField % block % domain % dminfo, size(real4d_temp), real4d_temp(:,1,1,1))
+ field_5dreal_ptr => field_cursor % real5dField
+ do while (associated(field_5dreal_ptr))
+ field_5dreal_ptr % array(j,:,:,:,:) = real4d_temp(:,:,:,:)
+ field_5dreal_ptr => field_5dreal_ptr % next
+ end do
+ else
+ call mpas_dmpar_bcast_reals(field_cursor % real5dField % block % domain % dminfo, size(real5d_temp), real5d_temp(:,1,1,1,1))
+ field_5dreal_ptr => field_cursor % real5dField
+ do while (associated(field_5dreal_ptr))
+ field_5dreal_ptr % array(:,:,:,:,:) = real5d_temp(:,:,:,:,:)
+ field_5dreal_ptr => field_5dreal_ptr % next
+ end do
+ end if
+ end if
+ end do
+
+ if (field_cursor % real5dField % isSuperArray) then
+ deallocate(real4d_temp)
+ else
+ deallocate(real5d_temp)
+ end if
+
+
else if (field_cursor % field_type == FIELD_0D_CHAR) then
!write(0,*) 'Reading in field '//trim(field_cursor % char0dField % fieldName)
Modified: branches/ocean_projects/cvmix_implementation/src/framework/mpas_timer.F
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/framework/mpas_timer.F 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/framework/mpas_timer.F 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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/ocean_projects/cvmix_implementation/src/registry/gen_inc.c
===================================================================
--- branches/ocean_projects/cvmix_implementation/src/registry/gen_inc.c 2013-03-11 14:47:30 UTC (rev 2579)
+++ branches/ocean_projects/cvmix_implementation/src/registry/gen_inc.c 2013-03-11 14:55:04 UTC (rev 2580)
@@ -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>