<p><b>mpetersen@lanl.gov</b> 2012-10-22 09:37:15 -0600 (Mon, 22 Oct 2012)</p><p>branch commit, leith_mrp: Merge trunk to branch. This includes pbc and bit-for-bit updates.<br>
</p><hr noshade><pre><font color="gray">Index: branches/ocean_projects/leith_mrp
===================================================================
--- branches/ocean_projects/leith_mrp 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp 2012-10-22 15:37:15 UTC (rev 2236)
Property changes on: branches/ocean_projects/leith_mrp
___________________________________________________________________
Modified: svn:mergeinfo
## -8,6 +8,8 ##
/branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/monotonic_advection:1499-1640
/branches/ocean_projects/monthly_forcing:1810-1867
+/branches/ocean_projects/option3_b4b_test:2201-2231
+/branches/ocean_projects/partial_bottom_cells:2172-2226
/branches/ocean_projects/split_explicit_mrp:1134-1138
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
## -20,4 +22,4 ##
/branches/omp_blocks/multiple_blocks:1803-2084
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
-/trunk/mpas:2182-2183
+/trunk/mpas:2182-2235
\ No newline at end of property
Modified: branches/ocean_projects/leith_mrp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/leith_mrp/namelist.input.ocean 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/namelist.input.ocean 2012-10-22 15:37:15 UTC (rev 2236)
@@ -9,12 +9,13 @@
/
&io
config_input_name = 'grid.nc'
- config_output_name = 'output..nc'
+ config_output_name = 'output.nc'
config_restart_name = 'restart.nc'
config_output_interval = '1_00:00:00'
config_frames_per_outfile = 1000000
config_pio_num_iotasks = 0
config_pio_stride = 1
+ config_write_output_on_startup = .true.
/
&decomposition
config_number_of_blocks = 0
@@ -31,6 +32,12 @@
config_pressure_type = 'pressure'
config_rho0 = 1014.65
/
+&partial_bottom_cells
+ config_alter_ICs_for_pbcs = 'off'
+ config_min_pbc_fraction = 0.10
+ config_check_ssh_consistency = .true.
+ config_check_zlevel_consistency = .false.
+/
&split_explicit_ts
config_n_ts_iter = 2
config_n_bcl_iter_beg = 1
Modified: branches/ocean_projects/leith_mrp/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -6229,10 +6229,12 @@
nearest_distance = current_distance
do i = 1, nEdgesOnCell(current_cell)
iCell = cellsOnCell(i,current_cell)
- d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0_RKIND)
- if (d < nearest_distance) then
- nearest_cell = iCell
- nearest_distance = d
+ if (iCell <= nCells) then
+ d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0_RKIND)
+ if (d < nearest_distance) then
+ nearest_cell = iCell
+ nearest_distance = d
+ end if
end if
end do
end do
@@ -6281,10 +6283,12 @@
end if
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
- d = sphere_distance(latEdge(iEdge), lonEdge(iEdge), target_lat, target_lon, 1.0_RKIND)
- if (d < nearest_distance) then
- nearest_edge = iEdge
- nearest_distance = d
+ if (iEdge <= nEdges) then
+ d = sphere_distance(latEdge(iEdge), lonEdge(iEdge), target_lat, target_lon, 1.0_RKIND)
+ if (d < nearest_distance) then
+ nearest_edge = iEdge
+ nearest_distance = d
+ end if
end if
end do
end do
Index: branches/ocean_projects/leith_mrp/src/core_ocean
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean 2012-10-22 15:37:15 UTC (rev 2236)
Property changes on: branches/ocean_projects/leith_mrp/src/core_ocean
___________________________________________________________________
Added: svn:mergeinfo
## -0,0 +1,27 ##
+/branches/atmos_physics/src/core_ocean:1672-1846
+/branches/cam_mpas_nh/src/core_ocean:1260-1270
+/branches/ocean_projects/ale_split_exp/src/core_ocean:1437-1483
+/branches/ocean_projects/ale_vert_coord/src/core_ocean:1225-1383
+/branches/ocean_projects/ale_vert_coord_new/src/core_ocean:1387-1428
+/branches/ocean_projects/gmvar/src/core_ocean:1214-1514,1517-1738
+/branches/ocean_projects/imp_vert_mix_error/src/core_ocean:1847-1887
+/branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean:754-986
+/branches/ocean_projects/monotonic_advection/src/core_ocean:1499-1640
+/branches/ocean_projects/monthly_forcing/src/core_ocean:1810-1867
+/branches/ocean_projects/option3_b4b_test/src/core_ocean:2201-2231
+/branches/ocean_projects/partial_bottom_cells/src/core_ocean:2172-2226
+/branches/ocean_projects/split_explicit_mrp/src/core_ocean:1134-1138
+/branches/ocean_projects/split_explicit_timestepping/src/core_ocean:1044-1097
+/branches/ocean_projects/vert_adv_mrp/src/core_ocean:704-745
+/branches/ocean_projects/vol_cons_RK_imp_mix/src/core_ocean:1965-1992
+/branches/ocean_projects/zstar_restart_new/src/core_ocean:1762-1770
+/branches/omp_blocks/block_decomp/src/core_ocean:1374-1569
+/branches/omp_blocks/ddt_reorg/src/core_ocean:1301-1414
+/branches/omp_blocks/halo/src/core_ocean:1570-1638
+/branches/omp_blocks/io/src/core_ocean:1639-1787
+/branches/omp_blocks/multiple_blocks/src/core_ocean:1803-2084
+/branches/omp_blocks/openmp_test/src/core_ocean:2107-2144
+/branches/omp_blocks/openmp_test/src/core_ocean_elements:2161-2201
+/branches/source_renaming/src/core_ocean:1082-1113
+/branches/time_manager/src/core_ocean:924-962
+/trunk/mpas/src/core_ocean:2182-2235
\ No newline at end of property
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/Registry 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/Registry 2012-10-22 15:37:15 UTC (rev 2236)
@@ -19,11 +19,12 @@
namelist character io config_restart_name restart.nc
namelist character io config_output_interval 24:00:00
namelist integer io config_frames_per_outfile 0
-namelist integer io config_pio_num_iotasks 0
-namelist integer io config_pio_stride 1
+namelist integer io config_pio_num_iotasks 0
+namelist integer io config_pio_stride 1
+namelist logical io config_write_output_on_startup true
namelist character decomposition config_block_decomp_file_prefix graph.info.part.
namelist integer decomposition config_number_of_blocks 0
-namelist logical decomposition config_explicit_proc_decomp .false.
+namelist logical decomposition config_explicit_proc_decomp false
namelist character decomposition config_proc_decomp_file_prefix graph.info.part.
namelist logical restart config_do_restart false
namelist character restart config_restart_interval none
@@ -31,6 +32,10 @@
namelist character grid config_pressure_type pressure
namelist real grid config_rho0 1028
namelist logical grid config_enforce_zstar_at_restart false
+namelist character partial_bottom_cells config_alter_ICs_for_pbcs zlevel_pbcs_off
+namelist real partial_bottom_cells config_min_pbc_fraction 0.10
+namelist logical partial_bottom_cells config_check_ssh_consistency true
+namelist logical partial_bottom_cells config_check_zlevel_consistency false
namelist integer split_explicit_ts config_n_ts_iter 2
namelist integer split_explicit_ts config_n_bcl_iter_beg 2
namelist integer split_explicit_ts config_n_bcl_iter_mid 2
@@ -64,7 +69,7 @@
namelist real hmix_leith config_leith_visc2_max 1000000.0
namelist character vmix config_vert_visc_type const
namelist character vmix config_vert_diff_type const
-namelist logical vmix config_implicit_vertical_mix .true.
+namelist logical vmix config_implicit_vertical_mix true
namelist real vmix config_convective_visc 1.0
namelist real vmix config_convective_diff 1.0
namelist real vmix config_bottom_drag_coeff 1.0e-3
@@ -166,7 +171,7 @@
var persistent real kiteAreasOnVertex ( vertexDegree nVertices ) 0 iro kiteAreasOnVertex mesh - -
var persistent real fEdge ( nEdges ) 0 iro fEdge mesh - -
var persistent real fVertex ( nVertices ) 0 iro fVertex mesh - -
-var persistent real h_s ( nCells ) 0 iro h_s mesh - -
+var persistent real bottomDepth ( nCells ) 0 iro bottomDepth mesh - -
% Space needed for advection
var persistent real deriv_two ( maxEdges2 TWO nEdges ) 0 - deriv_two mesh - -
@@ -196,8 +201,8 @@
var persistent integer maxLevelEdgeBot ( nEdges ) 0 - maxLevelEdgeBot mesh - -
var persistent integer maxLevelVertexTop ( nVertices ) 0 - maxLevelVertexTop mesh - -
var persistent integer maxLevelVertexBot ( nVertices ) 0 - maxLevelVertexBot mesh - -
-var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
-var persistent real referenceBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - referenceBottomDepthTopOfCell mesh - -
+var persistent real refBottomDepth ( nVertLevels ) 0 iro refBottomDepth mesh - -
+var persistent real refBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - refBottomDepthTopOfCell mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
var persistent real zstarWeight ( nVertLevels ) 0 - zstarWeight mesh - -
@@ -280,7 +285,7 @@
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
-var persistent real viscosity ( nVertLevels nCells Time ) 2 o viscosity state - -
+var persistent real viscosity ( nVertLevels nEdges Time ) 2 o viscosity state - -
% Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
@@ -312,3 +317,8 @@
var persistent real acc_uReconstructMeridionalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridionalVar state - -
var persistent real acc_u ( nVertLevels nEdges Time ) 2 o acc_u state - -
var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 o acc_uVar state - -
+
+% Sign fields, for openmp and bit reproducibility without branching statements.
+var persistent integer edgeSignOnCell ( maxEdges nCells ) 0 - edgeSignOnCell mesh - -
+var persistent integer edgeSignOnVertex ( maxEdges nVertices ) 0 - edgeSignOnVertex mesh - -
+var persistent integer kiteIndexOnCell ( maxEdges nCells ) 0 - kiteIndexOnCell mesh - -
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_equation_of_state_jm.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -96,7 +96,7 @@
real (kind=RKIND), dimension(:), pointer :: &
- referenceBottomDepth, pRefEOS
+ refBottomDepth, pRefEOS
real (kind=RKIND), dimension(:,:), intent(inout) :: &
rho
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
@@ -197,7 +197,7 @@
nCells = grid % nCells
maxLevelCell => grid % maxLevelCell % array
nVertLevels = grid % nVertLevels
- referenceBottomDepth => grid % referenceBottomDepth % array
+ refBottomDepth => grid % refBottomDepth % array
! Jackett and McDougall
@@ -214,14 +214,14 @@
allocate(pRefEOS(nVertLevels),p(nVertLevels),p2(nVertLevels))
! This could be put in the init routine.
- ! Note I am using referenceBottomDepth, so pressure on top level does
+ ! Note I am using refBottomDepth, so pressure on top level does
! not include SSH contribution. I am not sure if that matters, but
! POP does it the same way.
- depth = 0.5*referenceBottomDepth(1)
+ depth = 0.5*refBottomDepth(1)
pRefEOS(1) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ 0.100766*depth + 2.28405e-7*depth**2
do k = 2,nVertLevels
- depth = 0.5*(referenceBottomDepth(k)+referenceBottomDepth(k-1))
+ depth = 0.5*(refBottomDepth(k)+refBottomDepth(k-1))
pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ 0.100766*depth + 2.28405e-7*depth**2
enddo
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -106,10 +106,10 @@
if (.not. config_do_restart) call setup_sw_test_case(domain)
+ if (config_vert_grid_type.ne.'isopycnal') call ocn_init_vert_coord(domain)
+
call ocn_compute_max_level(domain)
- call ocn_init_z_level(domain)
-
if (config_enforce_zstar_at_restart) then
call ocn_init_h_zstar(domain)
endif
@@ -249,6 +249,7 @@
integer :: i, iEdge, iCell, k
integer :: err1
+ call ocn_setup_sign_and_index_fields(mesh)
call ocn_initialize_advection_rk(mesh, err)
call mpas_ocn_tracer_advection_coefficients(mesh, err1)
err = ior(err, err1)
@@ -351,7 +352,9 @@
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
write(0,*) 'Initial time ', trim(timeStamp)
- call ocn_write_output_frame(output_obj, output_frame, domain)
+ if (config_write_output_on_startup) then
+ call ocn_write_output_frame(output_obj, output_frame, domain)
+ endif
block_ptr => domain % blocklist
do while(associated(block_ptr))
@@ -390,7 +393,8 @@
if (mpas_is_alarm_ringing(clock, outputAlarmID, ierr=ierr)) then
call mpas_reset_clock_alarm(clock, outputAlarmID, ierr=ierr)
- ! output_frame will always be > 1 here unless it was reset after the maximum number of frames per outfile was reached
+ ! output_frame will always be > 1 here unless it was reset after the
+ ! maximum number of frames per outfile was reached.
if(output_frame == 1) then
call mpas_output_state_finalize(output_obj, domain % dminfo)
call mpas_output_state_init(output_obj, domain, "OUTPUT", trim(timeStamp))
@@ -531,8 +535,9 @@
end subroutine mpas_timestep!}}}
- subroutine ocn_init_z_level(domain)!{{{
- ! Initialize zlevel-type variables
+ subroutine ocn_init_vert_coord(domain)!{{{
+ ! Initialize zlevel-type variables and adjust initial conditions for
+ ! partial bottom cells.
use mpas_grid_types
use mpas_configure
@@ -540,16 +545,21 @@
implicit none
type (domain_type), intent(inout) :: domain
+ type (dm_info) :: dminfo
- integer :: i, iCell, iEdge, iVertex, k
+ integer :: i, iCell, iEdge, iVertex, k, nCells, num_tracers
type (block_type), pointer :: block
integer :: iTracer, cell, cell1, cell2
- real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, &
- referenceBottomDepthTopOfCell, zstarWeight, hZLevel
+ real (kind=RKIND) :: uhSum, hSum, hEdge1, zMidPBC
+
+ integer, dimension(:), pointer :: maxLevelCell
+ real (kind=RKIND), dimension(:), pointer :: refBottomDepth, &
+ refBottomDepthTopOfCell, zstarWeight, hZLevel, bottomDepth
+ real (kind=RKIND), dimension(:), allocatable :: minBottomDepth, minBottomDepthMid, zMidZLevel
real (kind=RKIND), dimension(:,:), pointer :: h
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer :: nVertLevels
! Initialize z-level grid variables from h, read in from input file.
@@ -557,26 +567,32 @@
do while (associated(block))
h => block % state % time_levs(1) % state % h % array
- referenceBottomDepth => block % mesh % referenceBottomDepth % array
- referenceBottomDepthTopOfCell => block % mesh % referenceBottomDepthTopOfCell % array
+ tracers => block % state % time_levs(1) % state % tracers % array
+ refBottomDepth => block % mesh % refBottomDepth % array
+ refBottomDepthTopOfCell => block % mesh % refBottomDepthTopOfCell % array
+ bottomDepth => block % mesh % bottomDepth % array
zstarWeight => block % mesh % zstarWeight % array
hZLevel => block % mesh % hZLevel % array
+ maxLevelCell => block % mesh % maxLevelCell % array
+
+ nCells = block % mesh % nCells
nVertLevels = block % mesh % nVertLevels
+ num_tracers = size(tracers, dim=1)
! mrp 120208 right now hZLevel is in the grid.nc file.
- ! We would like to transition to using referenceBottomDepth
+ ! We would like to transition to using refBottomDepth
! as the defining variable instead, and will transition soon.
! When the transition is done, hZLevel can be removed from
! registry and the following four lines deleted.
- referenceBottomDepth(1) = hZLevel(1)
+ refBottomDepth(1) = hZLevel(1)
do k = 2,nVertLevels
- referenceBottomDepth(k) = referenceBottomDepth(k-1) + hZLevel(k)
+ refBottomDepth(k) = refBottomDepth(k-1) + hZLevel(k)
end do
! TopOfCell needed where zero depth for the very top may be referenced.
- referenceBottomDepthTopOfCell(1) = 0.0
+ refBottomDepthTopOfCell(1) = 0.0
do k = 1,nVertLevels
- referenceBottomDepthTopOfCell(k+1) = referenceBottomDepth(k)
+ refBottomDepthTopOfCell(k+1) = refBottomDepth(k)
end do
! Initialization of zstarWeights. This determines how SSH perturbations
@@ -588,9 +604,7 @@
elseif (config_vert_grid_type.eq.'zstar') then
- do k = 1,nVertLevels
- zstarWeight(k) = hZLevel(k)
- enddo
+ zstarWeight = 1.0
elseif (config_vert_grid_type.eq.'zstarWeights') then
@@ -605,10 +619,117 @@
endif
+ ! Initial condition files (ocean.nc, produced by basin) include a realistic
+ ! bottomDepth variable and h,T,S variables for full thickness cells.
+ ! If running with pbcs, set config_alter_ICs_for_pbc='zlevel_pbcs_on'. Then thin pbc cells
+ ! will be changed, and h,T,S will be altered to match the pbcs.
+ ! If running without pbcs, set config_alter_ICs_for_pbc='zlevel_pbcs_off'. Then
+ ! bottomDepth will be altered so it is full cells everywhere.
+ ! If your input file does not include bottomDepth, the false option will
+ ! initialize bottomDepth correctly for a non-pbc run.
+
+
+ if (.not.config_do_restart) then
+
+ if (config_alter_ICs_for_pbcs.eq.'zlevel_pbcs_on') then
+
+ write (0,'(a)') ' Altering bottomDepth to avoid very thin cells.'
+ write (0,'(a)') ' Altering h and tracer initial conditions to conform with partial bottom cells.'
+
+ allocate(minBottomDepth(nVertLevels),minBottomDepthMid(nVertLevels),zMidZLevel(nVertLevels))
+
+ ! min_pbc_fraction restricts pbcs from being too small.
+ ! A typical value is 10%, so pbcs must occupy at least 10% of the cell thickness.
+ ! If min_pbc_fraction = 0.0, bottomDepth gives the actual depth for that cell.
+ ! If min_pbc_fraction = 1.0, bottomDepth reverts to discrete z-level depths, same
+ ! as partial_bottom_cells = .false.
+
+ do k=1,nVertLevels
+ minBottomDepth(k) = refBottomDepth(k) - (1.0-config_min_pbc_fraction)*hZLevel(k)
+ minBottomDepthMid(k) = 0.5*(minBottomDepth(k) + refBottomDepthTopOfCell(k))
+ zMidZLevel(k) = - 0.5*(refBottomDepth(k) + refBottomDepthTopOfCell(k))
+ enddo
+
+ do iCell=1,nCells
+ k = maxLevelCell(iCell)
+
+ if (bottomDepth(iCell).lt.minBottomDepthMid(k)) then
+ ! Round up to cell above
+ maxLevelCell(iCell) = maxLevelCell(iCell) - 1
+ bottomDepth(iCell) = refBottomDepth(maxLevelCell(iCell))
+ elseif (bottomDepth(iCell).lt.minBottomDepth(k)) then
+ ! Round down cell to the min_pbc_fraction.
+ bottomDepth(iCell) = minBottomDepth(k)
+ endif
+ k = maxLevelCell(iCell)
+
+ ! Alter thickness of bottom level to account for PBC
+ h(k,iCell) = bottomDepth(iCell) - refBottomDepthTopOfCell(k)
+
+ ! Linearly interpolate the initial T&S for new location of bottom cell for PBCs
+ zMidPBC = -0.5*(bottomDepth(iCell) + refBottomDepthTopOfCell(k))
+
+ do iTracer=1,num_tracers
+ tracers(iTracer,k,iCell) = tracers(iTracer,k,iCell) &
+ + (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell)) &
+ /(zMidZLevel(k-1)-zMidZLevel(k)) &
+ *(zMidPBC - zMidZLevel(k))
+ enddo
+
+ enddo
+
+ deallocate(minBottomDepth,zMidZLevel)
+
+ elseif (config_alter_ICs_for_pbcs.eq.'zlevel_pbcs_off') then
+
+ do iCell = 1,nCells
+ bottomDepth(iCell) = refBottomDepth(maxLevelCell(iCell))
+ enddo
+
+ elseif (config_alter_ICs_for_pbcs.eq.'off') then
+ ! No action taken. This is for isopycnal or sigma coordinates,
+ ! or if ICs were already altered upon start-up.
+
+ else
+
+ write (0,*) ' Incorrect choice of config_alter_ICs_for_pbcs.'
+ call mpas_dmpar_abort(dminfo)
+
+ endif
+ endif
+
+ if (config_check_ssh_consistency) then
+ do iCell = 1,nCells
+ ! Check if abs(ssh)>2m. If so, print warning.
+ if (abs(sum(h(1:maxLevelCell(iCell),iCell))-bottomDepth(iCell))>2.0) then
+ write (0,'(a)') ' Warning: abs(sum(h)-bottomDepth)>2m. Most likely, initial h does not match bottomDepth.'
+ write (0,*) ' iCell, K=maxLevelCell(iCell), bottomDepth(iCell),sum(h),bottomDepth,hZLevel(K),h(K): ', &
+ iCell, maxLevelCell(iCell), bottomDepth(iCell),sum(h(1:maxLevelCell(iCell),iCell)),bottomDepth(iCell), &
+ hZLevel(maxLevelCell(iCell)), h(maxLevelCell(iCell),iCell)
+ endif
+ enddo
+ endif
+
+ if (config_check_zlevel_consistency) then
+ do iCell = 1,nCells
+ ! Check that bottomDepth and maxLevelCell match. Some older grids do not have the bottomDepth variable.
+ if (bottomDepth(iCell) > refBottomDepth(maxLevelCell(iCell)).or. &
+ bottomDepth(iCell) < refBottomDepthTopOfCell(maxLevelCell(iCell))) then
+ write (0,'(a)') ' fatal error: bottomDepth and maxLevelCell do not match:'
+ write (0,'(a,2i5,10f10.2)') ' iCell, maxLevelCell(iCell), bottomDepth(iCell): ', &
+ iCell, maxLevelCell(iCell), bottomDepth(iCell)
+ write (0,'(a,10f10.2)') ' refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell)): ', &
+ refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell))
+ call mpas_dmpar_abort(dminfo)
+ endif
+
+ enddo
+ endif
+
block => block % next
end do
- end subroutine ocn_init_z_level!}}}
+ end subroutine ocn_init_vert_coord!}}}
subroutine ocn_init_split_timestep(domain)!{{{
! Initialize splitting variables
@@ -625,7 +746,7 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:), pointer :: refBottomDepth
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -635,7 +756,7 @@
do while (associated(block))
h => block % state % time_levs(1) % state % h % array
- referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ refBottomDepth => block % mesh % refBottomDepth % array
nVertLevels = block % mesh % nVertLevels
! Compute barotropic velocity at first timestep
@@ -651,7 +772,7 @@
if (config_filter_btr_mode) then
do iCell=1,block % mesh % nCells
block % state % time_levs(1) % state % h % array(1,iCell) &
- = block % mesh % referenceBottomDepth % array(1)
+ = block % mesh % refBottomDepth % array(1)
enddo
endif
@@ -735,7 +856,7 @@
real (kind=RKIND) :: hSum, sumZstarWeights
real (kind=RKIND), dimension(:), pointer :: hZLevel, zstarWeight, &
- referenceBottomDepth
+ refBottomDepth
real (kind=RKIND), dimension(:,:), pointer :: h
! Initialize z-level grid variables from h, read in from input file.
@@ -747,7 +868,7 @@
hZLevel => block % mesh % hZLevel % array
maxLevelCell => block % mesh % maxLevelCell % array
zstarWeight => block % mesh % zstarWeight % array
- referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ refBottomDepth => block % mesh % refBottomDepth % array
do iCell=1,block % mesh % nCells
! Compute the total column thickness, hSum, and the sum of zstar weights.
@@ -762,7 +883,7 @@
! where zeta is SSH and W_k are weights
do k = 1,maxLevelCell(iCell)
h(k,iCell) = hZLevel(k) &
- + (hSum - referenceBottomDepth(maxLevelCell(iCell))) &
+ + (hSum - refBottomDepth(maxLevelCell(iCell))) &
* zstarWeight(k)/sumZstarWeights
enddo
@@ -787,10 +908,6 @@
integer :: i, iCell, iEdge, iVertex, k
type (block_type), pointer :: block
- real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND) :: delta_rho, pi, latCenter, lonCenter, dist
- real (kind=RKIND) :: centerx, centery
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
integer, dimension(:), pointer :: &
@@ -960,6 +1077,72 @@
end subroutine ocn_compute_mesh_scaling!}}}
+ subroutine ocn_setup_sign_and_index_fields(mesh)!{{{
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnCell, edgesOnVertex, cellsOnVertex, cellsOnEdge, verticesOnCell, verticesOnEdge
+ integer, dimension(:,:), pointer :: edgeSignOnCell, edgeSignOnVertex, kiteIndexOnCell
+
+ integer :: nCells, nEdges, nVertices, vertexDegree
+ integer :: iCell, iEdge, iVertex, i, j, k
+
+ nCells = mesh % nCells
+ nEdges = mesh % nEdges
+ nVertices = mesh % nVertices
+ vertexDegree = mesh % vertexDegree
+
+ nEdgesOnCell => mesh % nEdgesOnCell % array
+ edgesOnCell => mesh % edgeSOnCell % array
+ edgesOnVertex => mesh % edgesOnVertex % array
+ cellsOnVertex => mesh % cellsOnVertex % array
+ cellsOnEdge => mesh % cellsOnEdge % array
+ verticesOnCell => mesh % verticesOnCell % array
+ verticesOnEdge => mesh % verticesOnEdge % array
+ edgeSignOnCell => mesh % edgeSignOnCell % array
+ edgeSignOnVertex => mesh % edgeSignOnVertex % array
+ kiteIndexOnCell => mesh % kiteIndexOnCell % array
+
+ edgeSignOnCell = 0.0_RKIND
+ edgeSignOnVertex = 0.0_RKIND
+ kiteIndexOnCell = 0.0_RKIND
+
+ do iCell = 1, nCells
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+
+ ! Vector points from cell 1 to cell 2
+ if(iCell == cellsOnEdge(1, iEdge)) then
+ edgeSignOnCell(i, iCell) = -1
+ else
+ edgeSignOnCell(i, iCell) = 1
+ end if
+
+ do j = 1, vertexDegree
+ if(cellsOnVertex(j, iVertex) == iCell) then
+ kiteIndexOnCell(i, iCell) = j
+ end if
+ end do
+ end do
+ end do
+
+ do iVertex = 1, nVertices
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+
+ ! Vector points from vertex 1 to vertex 2
+ if(iVertex == verticesOnEdge(1, iEdge)) then
+ edgeSignOnVertex(i, iVertex) = -1
+ else
+ edgeSignOnVertex(i, iVertex) = 1
+ end if
+ end do
+ end do
+
+ end subroutine ocn_setup_sign_and_index_fields!}}}
+
end module mpas_core
! vim: foldmethod=marker
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tendency.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tendency.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -405,15 +405,14 @@
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- verticesOnEdge, edgesOnEdge, edgesOnVertex,boundaryCell
+ verticesOnEdge, edgesOnEdge, edgesOnVertex,boundaryCell, kiteIndexOnCell, verticesOnCell, edgeSignOnVertex, edgeSignOnCell, edgesOnCell
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2, coef_3rd_order, r_tmp, invAreaCell1, invAreaCell2, invAreaTri1, invAreaTri2, invLength, h_vertex
real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, &
- referenceBottomDepth, ssh
+ bottomDepth, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, ssh
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -436,11 +435,11 @@
kev => s % kev % array
kevc => s % kevc % array
ke_edge => s % ke_edge % array
- Vor_edge => s % Vor_edge % array
- Vor_vertex => s % Vor_vertex % array
- Vor_cell => s % Vor_cell % array
- gradVor_n => s % gradVor_n % array
- gradVor_t => s % gradVor_t % array
+ Vor_edge => s % Vor_edge % array
+ Vor_vertex => s % Vor_vertex % array
+ Vor_cell => s % Vor_cell % array
+ gradVor_n => s % gradVor_n % array
+ gradVor_t => s % gradVor_t % array
rho => s % rho % array
MontPot => s % MontPot % array
pressure => s % pressure % array
@@ -455,20 +454,22 @@
verticesOnEdge => grid % verticesOnEdge % array
nEdgesOnCell => grid % nEdgesOnCell % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
edgesOnEdge => grid % edgesOnEdge % array
edgesOnVertex => grid % edgesOnVertex % array
dcEdge => grid % dcEdge % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
+ bottomDepth => grid % bottomDepth % array
fVertex => grid % fVertex % array
- referenceBottomDepth => grid % referenceBottomDepth % array
deriv_two => grid % deriv_two % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
maxLevelVertexBot => grid % maxLevelVertexBot % array
+ kiteIndexOnCell => grid % kiteIndexOnCell % array
+ verticesOnCell => grid % verticesOnCell % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -478,7 +479,10 @@
boundaryCell => grid % boundaryCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
!
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
@@ -579,40 +583,33 @@
divergence(:,:) = 0.0
ke(:,:) = 0.0
v(:,:) = 0.0
- do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+ do iVertex = 1, nVertices
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ do k = 1, maxLevelVertexBot(iVertex)
+ r_tmp = dcEdge(iEdge) * u(k, iEdge)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ circulation(k, iVertex) = circulation(k, iVertex) + edgeSignOnVertex(i, iVertex) * r_tmp
+ vorticity(k, iVertex) = vorticity(k, iVertex) + edgeSignOnVertex(i, iVertex) * r_tmp * invAreaTri1
+ end do
+ end do
+ end do
- invAreaTri1 = 1.0 / areaTriangle(vertex1)
- invAreaTri2 = 1.0 / areaTriangle(vertex2)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ r_tmp = dvEdge(iEdge) * u(k, iEdge) * invAreaCell1
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(cell1), 1.0)
- invAreaCell2 = 1.0 / max(areaCell(cell2), 1.0)
-
- do k=1,maxLevelEdgeBot(iEdge)
- ! Compute circulation and relative vorticity at each vertex
- r_tmp = dcEdge(iEdge) * u(k,iEdge)
- circulation(k,vertex1) = circulation(k,vertex1) - r_tmp
- circulation(k,vertex2) = circulation(k,vertex2) + r_tmp
-
- vorticity(k, vertex1) = vorticity(k, vertex1) - r_tmp * invAreaTri1
- vorticity(k, vertex2) = vorticity(k, vertex2) + r_tmp * invAreaTri2
-
- ! Compute the divergence at each cell center
- r_tmp = dvEdge(iEdge) * u(k, iEdge)
- divergence(k,cell1) = divergence(k,cell1) + r_tmp * invAreaCell1
- divergence(k,cell2) = divergence(k,cell2) - r_tmp * invAreaCell2
-
- ! Compute kinetic energy in each cell
- r_tmp = r_tmp * dcEdge(iEdge) * u(k,iEdge)
- ke(k,cell1) = ke(k,cell1) + 0.25 * r_tmp * invAreaCell1
- ke(k,cell2) = ke(k,cell2) + 0.25 * r_tmp * invAreaCell2
+ divergence(k, iCell) = divergence(k, iCell) - edgeSignOnCell(i, iCell) * r_tmp
+ ke(k, iCell) = ke(k, iCell) + 0.25 * r_tmp * dcEdge(iEdge) * u(k,iEdge)
+ end do
end do
+ end do
+ do iEdge=1,nEdges
! Compute v (tangential) velocities
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
@@ -629,29 +626,27 @@
! Compute kinetic energy in each vertex
!
kev(:,:) = 0.0; kevc(:,:) = 0.0
- do iEdge=1,nEdges*ke_vertex_flag
- do k=1,nVertLevels
- r_tmp = dcEdge(iEdge) * dvEdge(iEdge) * u(k, iEdge)**2
- kev(k,verticesOnEdge(1,iEdge)) = kev(k,verticesOnEdge(1,iEdge)) + r_tmp
- kev(k,verticesOnEdge(2,iEdge)) = kev(k,verticesOnEdge(2,iEdge)) + r_tmp
- end do
+ do iVertex = 1, nVertices*ke_vertex_flag
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ r_tmp = dcEdge(iEdge) * dvEdge(iEdge) * 0.25 / areaTriangle(iVertex)
+ do k = 1, nVertLevels
+ kev(k, iVertex) = kev(k, iVertex) + r_tmp * u(k, iEdge)**2
+ end do
+ end do
end do
- do iVertex = 1,nVertices*ke_vertex_flag
- do k=1,nVertLevels
- kev(k,iVertex) = kev(k,iVertex) / areaTriangle(iVertex) * 0.25
- enddo
- enddo
- do iVertex = 1, nVertices*ke_vertex_flag
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0)
- do k=1,nVertLevels
- kevc(k,iCell) = kevc(k,iCell) + kiteAreasOnVertex(i, iVertex) * kev(k, iVertex) * invAreaCell1
- enddo
- enddo
- enddo
+ do iCell = 1, nCells*ke_vertex_flag
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ j = kiteIndexOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+ do k = 1, nVertLevels
+ kevc(k, iCell) = kevc(k, iCell) + kiteAreasOnVertex(j, iVertex) * kev(k, iVertex) * invAreaCell1
+ end do
+ end do
+ end do
+
!
! Compute kinetic energy in each cell by blending ke and kevc
!
@@ -693,30 +688,26 @@
Vor_cell(:,:) = 0.0
Vor_edge(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- iEdge = edgesOnVertex(i,iVertex)
+ do iEdge = 1, nEdges
+ vertex1 = verticesOnEdge(1, iEdge)
+ vertex2 = verticesOnEdge(2, iEdge)
+ do k = 1, maxLevelEdgeBot(iEdge)
+ Vor_edge(k, iEdge) = 0.5 * (Vor_vertex(k, vertex1) + Vor_vertex(k, vertex2))
+ end do
+ end do
- !dwj: 02/23/12 arraCell(nCells+1) is still 0, this is a temporary fix
- invAreaCell1 = 1.0 / max(areaCell(iCell), 1.0)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
- ! Compute pv at cell centers
- ! ( this computes Vor_cell for all real cells and distance-1 ghost cells )
- do k = 1,maxLevelCell(iCell)
- Vor_cell(k,iCell) = Vor_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * Vor_vertex(k, iVertex) * invAreaCell1
- enddo
+ do i = 1, nEdgesOnCell(iCell)
+ j = kiteIndexOnCell(i, iCell)
+ iVertex = verticesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ Vor_cell(k, iCell) = Vor_cell(k, iCell) + kiteAreasOnVertex(j, iVertex) * Vor_vertex(k, iVertex) * invAreaCell1
+ end do
+ end do
+ end do
- ! Compute pv at the edges
- ! ( this computes Vor_edge at all edges bounding real cells )
- do k=1,maxLevelEdgeBot(iEdge)
- Vor_edge(k,iEdge) = Vor_edge(k,iEdge) + 0.5 * Vor_vertex(k,iVertex)
- enddo
- enddo
- enddo
-
-! gradVor_n(:,:) = 0.0
-! gradVor_t(:,:) = 0.0
do iEdge = 1,nEdges
cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)
@@ -780,9 +771,9 @@
pTop(1) = 0.0
! For isopycnal mode, p is the Montgomery Potential.
! At top layer it is g*SSH, where SSH may be off by a
- ! constant (ie, h_s can be relative to top or bottom)
+ ! constant (ie, bottomDepth can be relative to top or bottom)
MontPot(1,iCell) = gravity &
- * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
+ * (bottomDepth(iCell) + sum(h(1:nVertLevels,iCell)))
do k=2,nVertLevels
pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
@@ -811,10 +802,10 @@
! Compute zMid, the z-coordinate of the middle of the layer.
! This is used for the rho g grad z momentum term.
- ! Note the negative sign, since referenceBottomDepth is positive
+ ! Note the negative sign, since bottomDepth is positive
! and z-coordinates are negative below the surface.
k = maxLevelCell(iCell)
- zMid(k:nVertLevels,iCell) = -referenceBottomDepth(k) + 0.5*h(k,iCell)
+ zMid(k:nVertLevels,iCell) = -bottomDepth(iCell) + 0.5*h(k,iCell)
do k=maxLevelCell(iCell)-1, 1, -1
zMid(k,iCell) = zMid(k+1,iCell) &
@@ -831,13 +822,11 @@
!
do iCell=1,nCells
! Start at the bottom where we know the depth, and go up.
- ! The bottom depth for this cell is
- ! referenceBottomDepth(maxLevelCell(iCell)).
- ! Note the negative sign, since referenceBottomDepth is positive
+ ! The bottom depth for this cell is bottomDepth(iCell).
+ ! Note the negative sign, since bottomDepth is positive
! and z-coordinates are negative below the surface.
- ssh(iCell) = -referenceBottomDepth(maxLevelCell(iCell)) &
- + sum(h(1:maxLevelCell(iCell),iCell))
+ ssh(iCell) = - bottomDepth(iCell) + sum(h(1:maxLevelCell(iCell),iCell))
end do
@@ -873,20 +862,20 @@
type (mesh_type), intent(in) :: grid !< Input: Grid information
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum, invAreaCell
integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zstarWeight
+ dvEdge, areaCell, zstarWeight
real (kind=RKIND), dimension(:,:), pointer :: uTransport,h,wTop, h_edge
- real (kind=RKIND), dimension(:,:), allocatable:: div_hu
- real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col
+ real (kind=RKIND), dimension(:), allocatable:: div_hu, h_tend_col
+ real (kind=RKIND) :: div_hu_btr
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
- boundaryEdge, boundaryCell
+ boundaryEdge, boundaryCell, edgeSignOnCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot, maxLevelVertexTop
@@ -896,8 +885,11 @@
uTransport => s2 % uTransport % array
wTop => s2 % wTop % array
+ 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
@@ -907,67 +899,57 @@
nEdges = grid % nEdges
nVertLevels = grid % nVertLevels
- allocate(div_hu(nVertLevels,nCells+1), div_hu_btr(nCells+1), &
- h_tend_col(nVertLevels))
+ if (config_vert_grid_type.eq.'isopycnal') then
+ ! set vertical velocity to zero in isopycnal case
+ wTop=0.0_RKIND
+ return
+ end if
+
+ allocate(div_hu(nVertLevels), h_tend_col(nVertLevels))
+
!
! Compute div(h^{edge} u) for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
!
- div_hu(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- flux = uTransport(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- div_hu(k,cell1) = div_hu(k,cell1) + flux
- div_hu(k,cell2) = div_hu(k,cell2) - flux
- end do
- end do
do iCell=1,nCells
- div_hu_btr(iCell) = 0.0
- do k=1,maxLevelCell(iCell)
- div_hu(k,iCell) = div_hu(k,iCell) / areaCell(iCell)
- div_hu_btr(iCell) = div_hu_btr(iCell) + div_hu(k,iCell)
- end do
- end do
+ div_hu(:) = 0.0_RKIND
+ div_hu_btr = 0.0_RKIND
+ hSum = 0.0_RKIND
+ invAreaCell = 1.0_RKIND / areaCell(iCell)
- !
- ! vertical velocity through layer interface
- !
- !dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
- if (config_vert_grid_type.eq.'isopycnal') then
- ! set vertical velocity to zero in isopycnal case
- wTop=0.0
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
- else ! zlevel or zstar type vertical grid
+ do k = 1, maxLevelEdgeBot(iEdge)
+ flux = uTransport(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
+ flux = edgeSignOnCell(i, iCell) * flux * invAreaCell
+ div_hu(k) = div_hu(k) - flux
+ div_hu_btr = div_hu_btr - flux
+ end do
+ end do
- do iCell=1,nCells
+ do k = 1, maxLevelCell(iCell)
+ h_tend_col(k) = - zstarWeight(k) * h(k, iCell) * div_hu_btr
+ hSum = hSum + zstarWeight(k) * h(k, iCell)
+ end do
- hSum = 0.0
- do k=1,maxLevelCell(iCell)
- h_tend_col(k) = - zstarWeight(k)*h(k,iCell)*div_hu_btr(iCell)
- hSum = hSum + zstarWeight(k)*h(k,iCell)
- end do
- if(hSum > 0.0) then
- h_tend_col = h_tend_col / hSum
- else
- end if
+ if(hSum > 0.0) then
+ h_tend_col = h_tend_col / hSum
+ end if
- ! Vertical velocity through layer interface at top and
- ! bottom is zero.
- wTop(1,iCell) = 0.0
- wTop(maxLevelCell(iCell)+1,iCell) = 0.0
- do k=maxLevelCell(iCell),2,-1
- wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k,iCell) - h_tend_col(k)
- end do
+ ! Vertical velocity through layer interface at top and
+ ! bottom is zero.
+ wTop(1,iCell) = 0.0_RKIND
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND
+ do k=maxLevelCell(iCell),2,-1
+ wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k) - h_tend_col(k)
end do
+ end do
- endif
+ deallocate(div_hu, h_tend_col)
- deallocate(div_hu, div_hu_btr, h_tend_col)
-
end subroutine ocn_wtop!}}}
!***********************************************************************
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_test_cases.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_test_cases.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_test_cases.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -342,10 +342,10 @@
do iCell=1,grid % nCells
if (grid % lonCell % array(iCell) < 0.0) grid % lonCell % array(iCell) = grid % lonCell % array(iCell) + 2.0 * pii
r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- grid % h_s % array(iCell) = hs0 * (1.0 - r/rr)
+ grid % bottomDepth % array(iCell) = hs0 * (1.0 - r/rr)
end do
! output about mountain
-print *, 'h_s',minval(grid % h_s % array),sum(grid % h_s % array)/grid % nCells, maxval(grid % h_s % array)
+print *, 'bottomDepth',minval(grid % bottomDepth % array),sum(grid % bottomDepth % array)/grid % nCells, maxval(grid % bottomDepth % array)
!
! Initialize tracer fields
@@ -372,7 +372,7 @@
)**2.0 &
) / &
gravity
- state % h % array(1,iCell) = state % h % array(1,iCell) - grid % h_s % array(iCell)
+ state % h % array(1,iCell) = state % h % array(1,iCell) - grid % bottomDepth % array(iCell)
end do
end subroutine sw_test_case_5
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_thick_hadv.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_thick_hadv.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_thick_hadv.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -101,13 +101,13 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k
+ integer :: iEdge, nEdges, cell1, cell2, nVertLevels, k, i
integer :: iCell, nCells
- integer, dimension(:), pointer :: maxLevelEdgeBot, MaxLevelCell
- integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelEdgeBot, MaxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgeSignOnCell
- real (kind=RKIND) :: flux, invAreaCell1, invAreaCell2
+ real (kind=RKIND) :: flux, invAreaCell, invAreaCell1, invAreaCell2
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
!-----------------------------------------------------------------
@@ -130,20 +130,20 @@
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,maxLevelEdgeBot(iEdge)
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend(k,cell1) = tend(k,cell1) - flux
- tend(k,cell2) = tend(k,cell2) + flux
- end do
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
+ do iCell = 1, nCells
+ invAreaCell = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelEdgeBot(iEdge)
+ flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k, iEdge)
+ tend(k, iCell) = tend(k, iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell
+ end do
+ end do
end do
- do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- tend(k,iCell) = tend(k,iCell) / areaCell(iCell)
- end do
- end do
!--------------------------------------------------------------------
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -408,21 +408,63 @@
! config_btr_gam1_uWt1=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam1_uWt1= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
+
+ do iCell = 1, block % mesh % nCells
+ do i = 1, block % mesh % nEdgesOnCell % array(iCell)
+ iEdge = block % mesh % edgesOnCell % array(i, iCell)
+
+ cell1 = block % mesh % cellsOnEdge % array(1, iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2, iEdge)
+
+ sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
+
+ ! method 0: orig, works only without pbc:
+ !hSum = sshEdge + block % mesh % refBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ ! method 1, matches method 0 without pbcs, works with pbcs.
+ hSum = sshEdge + min(block % mesh % bottomDepth % array(cell1), &
+ block % mesh % bottomDepth % array(cell2))
+
+ ! method 2: may be better than method 1.
+ ! Take average of full thickness at two neighboring cells.
+ !hSum = sshEdge + 0.5 *( block % mesh % bottomDepth % array(cell1) &
+ ! + block % mesh % bottomDepth % array(cell2) )
+
+
+ flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * hSum
+
+ block % tend % ssh % array(iCell) = block % tend % ssh % array(iCell) + block % mesh % edgeSignOncell % array(i, iCell) * flux &
+ * block % mesh % dvEdge % array(iEdge)
+
+ end do
+ end do
+
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
-
+
sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
- hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
-
+
+ ! method 0: orig, works only without pbc:
+ !hSum = sshEdge + block % mesh % refBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ ! method 1, matches method 0 without pbcs, works with pbcs.
+ hSum = sshEdge + min(block % mesh % bottomDepth % array(cell1), &
+ block % mesh % bottomDepth % array(cell2))
+
+ ! method 2: may be better than method 1.
+ ! take average of full thickness at two neighboring cells
+ !hSum = sshEdge + 0.5 *( block % mesh % bottomDepth % array(cell1) &
+ ! + block % mesh % bottomDepth % array(cell2) )
+
flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
* hSum
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
-
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ FBtr_coeff*flux
end do
@@ -452,6 +494,8 @@
block => domain % blocklist
do while (associated(block))
+ allocate(utemp(block % mesh % nEdges+1))
+ uTemp(:) = block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -461,7 +505,8 @@
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
eoe = block % mesh % edgesOnEdge % array(i,iEdge)
CoriolisTerm = CoriolisTerm + block % mesh % weightsOnEdge % array(i,iEdge) &
- * block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ !* block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ * uTemp(eoe) &
* block % mesh % fEdge % array(eoe)
end do
@@ -478,6 +523,7 @@
+ dt/config_n_btr_subcycles *(CoriolisTerm - gravity *(sshCell2 - sshCell1) /block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1,iEdge)
end do
+ deallocate(uTemp)
block => block % next
end do ! block
@@ -502,6 +548,45 @@
! config_btr_gam3_uWt2=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
+
+ do iCell = 1, block % mesh % nCells
+ do i = 1, block % mesh % nEdgesOnCell % array(iCell)
+ iEdge = block % mesh % edgesOnCell % array(i, iCell)
+
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ ! SSH is a linear combination of SSHold and SSHnew.
+ sshCell1 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
+ sshCell2 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+
+ sshEdge = 0.5 * (sshCell1 + sshCell2)
+
+ ! method 0: orig, works only without pbc:
+ !hSum = sshEdge + block % mesh % refBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ ! method 1, matches method 0 without pbcs, works with pbcs.
+ hSum = sshEdge + min(block % mesh % bottomDepth % array(cell1), &
+ block % mesh % bottomDepth % array(cell2))
+
+ ! method 2: may be better than method 1.
+ ! take average of full thickness at two neighboring cells
+ !hSum = sshEdge + 0.5 *( block % mesh % bottomDepth % array(cell1) &
+ ! + block % mesh % bottomDepth % array(cell2) )
+
+
+ flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * hSum
+
+ block % tend % ssh % array(iCell) = block % tend % ssh % array(iCell) + block % mesh % edgeSignOnCell % array(i, iCell) * flux &
+ * block % mesh % dvEdge % array(iEdge)
+
+ end do
+ end do
+
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -511,17 +596,24 @@
+ config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
sshCell2 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+ sshEdge = 0.5 * (sshCell1 + sshCell2)
- sshEdge = 0.5 * (sshCell1 + sshCell2)
- hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+ ! method 0: orig, works only without pbc:
+ !hSum = sshEdge + block % mesh % refBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
+
+ ! method 1, matches method 0 without pbcs, works with pbcs.
+ hSum = sshEdge + min(block % mesh % bottomDepth % array(cell1), &
+ block % mesh % bottomDepth % array(cell2))
+
+ ! method 2, better, I think.
+ ! take average of full thickness at two neighboring cells
+ !hSum = sshEdge + 0.5 *( block % mesh % bottomDepth % array(cell1) &
+ ! + block % mesh % bottomDepth % array(cell2) )
flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
* hSum
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
-
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) + flux
end do
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -64,7 +64,7 @@
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
integer, dimension(:), pointer :: nAdvCellsForEdge, maxLevelCell, maxLevelEdgeTop, nEdgesOnCell
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask, edgeSignOnCell, edgesOnCell
real (kind=RKIND) :: flux_upwind, tracer_min_new, tracer_max_new, tracer_upwind_new, scale_factor
real (kind=RKIND) :: flux, tracer_weight, invDvEdge, invAreaCell1, invAreaCell2
@@ -95,6 +95,8 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
highOrderAdvectionMask => grid % highOrderAdvectionMask % array
lowOrderAdvectionMask => grid % lowOrderAdvectionMask % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
nCells = grid % nCells
nCellsSolve = grid % nCellsSolve
@@ -230,8 +232,8 @@
! flux_incoming (k,iCell) = -(min(0.,high_order_vert_flux(k+1,iCell))-max(0.,high_order_vert_flux(k,iCell)))
! flux_outgoing(k,iCell) = -(max(0.,high_order_vert_flux(k+1,iCell))-min(0.,high_order_vert_flux(k,iCell)))
- flux_incoming (k, iCell) = max(0.0, high_order_vert_flux(k+1, iCell)) - min(0.0, high_order_vert_flux(k, iCell))
- flux_outgoing(k, iCell) = min(0.0, high_order_vert_flux(k+1, iCell)) - max(0.0, high_order_vert_flux(k, iCell))
+ flux_incoming (k, iCell) = max(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - min(0.0_RKIND, high_order_vert_flux(k, iCell))
+ flux_outgoing(k, iCell) = min(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - max(0.0_RKIND, high_order_vert_flux(k, iCell))
end do ! k Loop
end do ! iCell Loop
@@ -249,18 +251,27 @@
do k = 1, maxLevelEdgeTop(iEdge)
flux_upwind = dvEdge(iEdge) * (max(0.,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.,uh(k,iEdge))*tracer_cur(k,cell2))
high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) - flux_upwind
-
- upwind_tendency(k,cell1) = upwind_tendency(k,cell1) - flux_upwind * invAreaCell1
- upwind_tendency(k,cell2) = upwind_tendency(k,cell2) + flux_upwind * invAreaCell2
-
- ! Accumulate remaining high order fluxes
- flux_outgoing(k,cell1) = flux_outgoing(k,cell1) - max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
- flux_incoming (k,cell1) = flux_incoming (k,cell1) - min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
- flux_outgoing(k,cell2) = flux_outgoing(k,cell2) + min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
- flux_incoming (k,cell2) = flux_incoming (k,cell2) + max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
end do ! k loop
end do ! iEdge loop
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ flux_upwind = dvEdge(iEdge) * (max(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.0_RKIND,uh(k,iEdge))*tracer_cur(k,cell2))
+
+ upwind_tendency(k,iCell) = upwind_tendency(k,iCell) + edgeSignOncell(i, iCell) * flux_upwind * invAreaCell1
+
+ ! Accumulate remaining high order fluxes
+ flux_outgoing(k,iCell) = flux_outgoing(k,iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ flux_incoming(k,iCell) = flux_incoming(k,iCell) - edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ end do
+ end do
+ end do
+
! Build the factors for the FCT
! Computed using the bounds that were computed previously, and the bounds on the newly updated value
! Factors are placed in the flux_incoming and flux_outgoing arrays
@@ -304,24 +315,20 @@
end do ! iCell loop
! Accumulate the scaled high order horizontal tendencies
- do iEdge = 1, nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
- do k=1,maxLevelEdgeTop(iEdge)
- tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+ if(config_check_monotonicity) then
+ tracer_new(k, iCell) = tracer_new(k, iCell) + edgeSignOnCell(i, iCell) * high_order_horiz_flux(k, iEdge) * invAreaCell1
+ end if
+ end do
+ end do
+ end do
- if (config_check_monotonicity) then
- !tracer_new holds a tendency for now.
- tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
- end if
- end do ! k loop
- end do ! iEdge loop
-
! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
do iCell = 1, nCellsSolve
do k = 1,maxLevelCell(iCell)
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -107,13 +107,13 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, nEdges, nVertLevels, cell1, cell2
- integer :: k, iTracer, num_tracers
+ integer :: iCell, iEdge, nCells, nEdges, nVertLevels, cell1, cell2
+ integer :: i, k, iTracer, num_tracers
integer, dimension(:,:), allocatable :: boundaryMask
- integer, dimension(:), pointer :: maxLevelEdgeTop
- integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask
+ integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask, edgesOnCell, edgeSignOnCell
real (kind=RKIND) :: invAreaCell1, invAreaCell2
real (kind=RKIND) :: tracer_turb_flux, flux, r_tmp
@@ -134,6 +134,7 @@
if (.not.del2On) return
nEdges = grid % nEdges
+ nCells = grid % nCells
nVertLevels = grid % nVertLevels
num_tracers = size(tracers, dim=1)
@@ -145,31 +146,37 @@
dcEdge => grid % dcEdge % array
meshScalingDel2 => grid % meshScalingDel2 % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0/areaCell(cell1)
- invAreaCell2 = 1.0/areaCell(cell2)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOncell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
-
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers
+ r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
+
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers
! \kappa_2 </font>
<font color="red">abla \phi on edge
- tracer_turb_flux = tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)
+ tracer_turb_flux = tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1)
! div(h \kappa_2 </font>
<font color="gray">abla \phi) at cell center
- flux = h_edge(k,iEdge) * tracer_turb_flux * edgeMask(k, iEdge) * r_tmp
+ flux = h_edge(k, iEdge) * tracer_turb_flux * edgeMask(k, iEdge) * r_tmp
- tend(iTracer,k,cell1) = tend(iTracer,k,cell1) + flux * invAreaCell1
- tend(iTracer,k,cell2) = tend(iTracer,k,cell2) - flux * invAreaCell2
- end do
- end do
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) - edgeSignOnCell(i, iCell) * flux * invAreaCell1
+ end do
+ end do
+ end do
end do
+
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del2_tend!}}}
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -108,10 +108,10 @@
!-----------------------------------------------------------------
integer :: iEdge, nEdges, num_tracers, nVertLevels, nCells
- integer :: iTracer, k, iCell, cell1, cell2
+ integer :: iTracer, k, iCell, cell1, cell2, i
- integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
- integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge, edgesOnCell, edgeSignOnCell
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux, flux, invdcEdge, r_tmp1, r_tmp2
@@ -148,56 +148,55 @@
edgeMask => grid % edgeMask % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+
allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
delsq_tracer(:,:,:) = 0.0
! first del2: div(h </font>
<font color="red">abla \phi) at cell center
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ invdcEdge = dvEdge(iEdge) / dcEdge(iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- invdcEdge = 1.0 / dcEdge(iEdge)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers * edgeMask(k, iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+ r_tmp1 = invdcEdge * h_edge(k, iEdge) * tracers(iTracer, k, cell1)
+ r_tmp2 = invdcEdge * h_edge(k, iEdge) * tracers(iTracer, k, cell2)
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers * edgeMask(k, iEdge)
-
- r_tmp1 = dvEdge(iEdge) * h_edge(k,iEdge) * invdcEdge
-
- r_tmp2 = r_tmp1 * tracers(iTracer,k,cell2)
- r_tmp1 = r_tmp1 * tracers(iTracer,k,cell1)
-
- delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) + (r_tmp2 - r_tmp1) * invAreaCell1
- delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) - (r_tmp2 - r_tmp1) * invAreaCell2
- end do
- end do
+ delsq_tracer(iTracer, k, iCell) = delsq_tracer(iTracer, k, iCell) - edgeSignOnCell(i, iCell) * (r_tmp2 - r_tmp1) * invAreaCell1
+ end do
+ end do
+ end do
end do
! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ cell1 = cellsOnEdge(1, iEdge)
+ cell2 = cellsOnedge(2, iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+ invdcEdge = meshScalingDel4(iEdge) * dvEdge(iEdge) * eddyDiff4 / dcEdge(iEdge)
- invdcEdge = 1.0 / dcEdge(iEdge)
+ do k = 1, maxLevelEdgeTop(iEdge)
+ do iTracer = 1, num_tracers * edgeMask(k, iEdge)
+ tracer_turb_flux = (delsq_tracer(iTracer, k, cell2) - delsq_tracer(iTracer, k, cell1))
+
+ flux = tracer_turb_flux * invdcEdge
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,num_tracers * edgeMask(k,iEdge)
- tracer_turb_flux = meshScalingDel4(iEdge) * eddyDiff4 &
- * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) &
- * invdcEdge
-
- flux = dvEdge (iEdge) * tracer_turb_flux
-
- tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux * invAreaCell1
- tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux * invAreaCell2
- enddo
- enddo
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * invAreaCell1
+ end do
+ end do
+ end do
end do
deallocate(delsq_tracer)
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_del2.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -165,8 +165,7 @@
tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * visc2 * u_diffusion
- viscosity(k,cell1) = viscosity(k,cell1) + visc2/nEdgesOnCell(cell1)
- viscosity(k,cell2) = viscosity(k,cell2) + visc2/nEdgesOnCell(cell2)
+ viscosity(k,iEdge) = viscosity(k,iEdge) + visc2
end do
end do
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -112,22 +112,22 @@
!
!-----------------------------------------------------------------
- integer :: iEdge, cell1, cell2, vertex1, vertex2, k
+ integer :: iEdge, cell1, cell2, vertex1, vertex2, k, i
integer :: iCell, iVertex
- integer :: nVertices, nVertLevels, nCells, nEdges, nEdgesSolve
+ integer :: nVertices, nVertLevels, nCells, nEdges, nEdgesSolve, vertexDegree
integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelVertexBot, &
- maxLevelCell
- integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
+ maxLevelCell, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask, edgesOnVertex, edgesOnCell, edgeSignOnVertex, edgeSignOnCell
real (kind=RKIND) :: u_diffusion, invAreaCell1, invAreaCell2, invAreaTri1, &
- invAreaTri2, invDcEdge, invDvEdge, r_tmp, delsq_u
+ invAreaTri2, invDcEdge, invDvEdge, r_tmp
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaTriangle, &
meshScalingDel4, areaCell
real (kind=RKIND), dimension(:,:), allocatable :: delsq_divergence, &
- delsq_circulation, delsq_vorticity
+ delsq_circulation, delsq_vorticity, delsq_u
err = 0
@@ -138,6 +138,8 @@
nEdgesSolve = grid % nEdgessolve
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
+
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelVertexBot => grid % maxLevelVertexBot % array
maxLevelCell => grid % maxLevelCell % array
@@ -149,43 +151,57 @@
areaCell => grid % areaCell % array
meshScalingDel4 => grid % meshScalingDel4 % array
edgeMask => grid % edgeMask % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+ allocate(delsq_u(nVertLEvels, nEdges+1))
allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_vorticity(nVertLevels, nVertices+1))
+ delsq_u(:,:) = 0.0
delsq_vorticity(:,:) = 0.0
delsq_divergence(:,:) = 0.0
- do iEdge=1,nEdges
+ !Compute delsq_u
+ do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- invAreaTri1 = 1.0 / areaTriangle(vertex1)
- invAreaTri2 = 1.0 / areaTriangle(vertex2)
-
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
-
invDcEdge = 1.0 / dcEdge(iEdge)
invDvEdge = 1.0 / dvEdge(iEdge)
do k=1,maxLevelEdgeTop(iEdge)
! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- delsq_u = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
+ delsq_u(k, iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
-viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDcEdge * sqrt(3.0) ! TDR
+ end do
+ end do
- ! vorticity using </font>
<font color="red">abla^2 u
- r_tmp = dcEdge(iEdge) * delsq_u
- delsq_vorticity(k,vertex1) = delsq_vorticity(k,vertex1) - r_tmp * invAreaTri1
- delsq_vorticity(k,vertex2) = delsq_vorticity(k,vertex2) + r_tmp * invAreaTri2
+ ! Compute delsq_vorticity
+ do iVertex = 1, nVertices
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
+ do i = 1, vertexDegree
+ iEdge = edgesOnVertex(i, iVertex)
+ do k = 1, maxLevelVertexBot(iVertex)
+ delsq_vorticity(k, iVertex) = delsq_vorticity(k, iVertex) + edgeSignOnVertex(i, iVertex) * dcEdge(iEdge) * delsq_u(k, iEdge) * invAreaTri1
+ end do
+ end do
+ end do
- ! Divergence using </font>
<font color="gray">abla^2 u
- r_tmp = dvEdge(iEdge) * delsq_u
- delsq_divergence(k, cell1) = delsq_divergence(k,cell1) + r_tmp * invAreaCell1
- delsq_divergence(k, cell2) = delsq_divergence(k,cell2) - r_tmp * invAreaCell2
+ ! Compute delsq_divergence
+ do iCell = 1, nCells
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+ do k = 1, maxLevelCell(iCell)
+ delsq_divergence(k, iCell) = delsq_divergence(k, iCell) - edgeSignOnCell(i, iCell) * dvEdge(iEdge) * delsq_u(k, iEdge) * invAreaCell1
+ end do
end do
end do
@@ -209,6 +225,7 @@
end do
end do
+ deallocate(delsq_u)
deallocate(delsq_divergence)
deallocate(delsq_vorticity)
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_leith.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vel_hmix_leith.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -167,8 +167,7 @@
tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * visc2 * u_diffusion
- viscosity(k,cell1) = viscosity(k,cell1) + visc2/nEdgesOnCell(cell1)
- viscosity(k,cell2) = viscosity(k,cell2) + visc2/nEdgesOnCell(cell2)
+ viscosity(k,iEdge) = viscosity(k,iEdge) + visc2
end do
end do
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -427,13 +427,13 @@
!
!-----------------------------------------------------------------
- integer :: nVertLevels, nCells, nEdges, iCell, iEdge, k
+ integer :: nVertLevels, nCells, nEdges, iCell, iEdge, k, i
integer :: cell1, cell2
- integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot
- integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOncell, edgeSignOnCell
- real (kind=RKIND) :: coef
+ real (kind=RKIND) :: coef, invAreaCell
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaCell
real (kind=RKIND), dimension(:,:), allocatable :: drhoTopOfCell, du2TopOfCell, &
drhoTopOfEdge, du2TopOfEdge
@@ -453,6 +453,9 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
areaCell => grid % areaCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
allocate( &
drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges), &
@@ -498,21 +501,16 @@
! interpolate du2TopOfEdge to du2TopOfCell
du2TopOfCell = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=2,maxLevelEdgeBot(iEdge)
- du2TopOfCell(k,cell1) = du2TopOfCell(k,cell1) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
- du2TopOfCell(k,cell2) = du2TopOfCell(k,cell2) &
- + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k,iEdge)
- end do
+ do iCell = 1, nCells
+ invAreaCell = 1.0 / areaCell(iCell)
+ do i = 1, nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i, iCell)
+
+ do k = 2, maxLevelEdgeBot(iEdge)
+ du2TopOfCell(k, iCell) = du2TopOfCell(k, iCell) + 0.5 * dcEdge(iEdge) * dvEdge(iEdge) * du2TopOfEdge(k, iEdge) * invAreaCell
+ end do
+ end do
end do
- do iCell = 1,nCells
- do k = 2,maxLevelCell(iCell)
- du2TopOfCell(k,iCell) = du2TopOfCell(k,iCell) / areaCell(iCell)
- end do
- end do
! compute RiTopOfEdge using drhoTopOfEdge and du2TopOfEdge
! coef = -g/rho_0/2
Modified: branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F
===================================================================
--- branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2012-10-22 15:37:15 UTC (rev 2236)
@@ -177,22 +177,22 @@
integer :: k, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:), pointer :: refBottomDepth
err = 0
if(.not.tanhViscOn) return
nVertLevels = grid % nVertLevels
- referenceBottomDepth => grid % referenceBottomDepth % array
+ refBottomDepth => grid % refBottomDepth % array
- ! referenceBottomDepth is used here for simplicity. Using zMid and h, which
+ ! refBottomDepth is used here for simplicity. Using zMid and h, which
! vary in time, would give the exact location of the top, but it
! would only change the diffusion value very slightly.
vertViscTopOfEdge = 0.0
do k=2,nVertLevels
vertViscTopOfEdge(k,:) = -(config_max_visc_tanh-config_min_visc_tanh)/2.0 &
- *tanh((referenceBottomDepth(k-1)+config_ZMid_tanh) &
+ *tanh((refBottomDepth(k-1)+config_ZMid_tanh) &
/config_zWidth_tanh) &
+ (config_max_visc_tanh+config_min_visc_tanh)/2
end do
@@ -250,22 +250,22 @@
integer :: k, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth
+ real (kind=RKIND), dimension(:), pointer :: refBottomDepth
err = 0
if(.not.tanhDiffOn) return
nVertLevels = grid % nVertLevels
- referenceBottomDepth => grid % referenceBottomDepth % array
+ refBottomDepth => grid % refBottomDepth % array
- ! referenceBottomDepth is used here for simplicity. Using zMid and h, which
+ ! refBottomDepth is used here for simplicity. Using zMid and h, which
! vary in time, would give the exact location of the top, but it
! would only change the diffusion value very slightly.
vertDiffTopOfCell = 0.0
do k=2,nVertLevels
vertDiffTopOfCell(k,:) = -(config_max_diff_tanh-config_min_diff_tanh)/2.0 &
- *tanh((referenceBottomDepth(k-1)+config_ZMid_tanh) &
+ *tanh((refBottomDepth(k-1)+config_ZMid_tanh) &
/config_zWidth_tanh) &
+ (config_max_diff_tanh+config_min_diff_tanh)/2
end do
Modified: branches/ocean_projects/leith_mrp/src/registry/gen_inc.c
===================================================================
--- branches/ocean_projects/leith_mrp/src/registry/gen_inc.c 2012-10-19 22:24:29 UTC (rev 2235)
+++ branches/ocean_projects/leith_mrp/src/registry/gen_inc.c 2012-10-22 15:37:15 UTC (rev 2236)
@@ -143,8 +143,8 @@
fortprintf(fd, " call mpas_dmpar_abort(dminfo)</font>
<font color="black">");
fortprintf(fd, " else if (ierr < 0) then</font>
<font color="black">");
fortprintf(fd, " write(0,*) \'Namelist record &%s not found; using default values for this namelist\'\'s variables\'</font>
<font color="red">",nls_ptr->record);
- fortprintf(fd, " rewind(funit)</font>
<font color="black">");
fortprintf(fd, " end if</font>
<font color="blue">");
+ fortprintf(fd, " rewind(funit)</font>
<font color="black">");
dict_insert(dictionary, nls_ptr->record);
}
</font>
</pre>