<p><b>mpetersen@lanl.gov</b> 2012-10-19 13:42:08 -0600 (Fri, 19 Oct 2012)</p><p>branch commit, option3_b4b_test. Merge trunk to branch. Trunk now includes today's partial bottom cell merge.<br>
</p><hr noshade><pre><font color="gray">Index: branches/ocean_projects/option3_b4b_test
===================================================================
--- branches/ocean_projects/option3_b4b_test 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test 2012-10-19 19:42:08 UTC (rev 2229)
Property changes on: branches/ocean_projects/option3_b4b_test
___________________________________________________________________
Modified: svn:mergeinfo
## -8,6 +8,7 ##
/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/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,3 +21,4 ##
/branches/omp_blocks/multiple_blocks:1803-2084
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+/trunk/mpas:2202-2228
\ No newline at end of property
Modified: branches/ocean_projects/option3_b4b_test/namelist.input.ocean
===================================================================
--- branches/ocean_projects/option3_b4b_test/namelist.input.ocean 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/namelist.input.ocean 2012-10-19 19:42:08 UTC (rev 2229)
@@ -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/option3_b4b_test/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-10-19 19:42:08 UTC (rev 2229)
@@ -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/option3_b4b_test/src/core_ocean
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean 2012-10-19 19:42:08 UTC (rev 2229)
Property changes on: branches/ocean_projects/option3_b4b_test/src/core_ocean
___________________________________________________________________
Modified: svn:mergeinfo
## -22,3 +22,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:2202-2228
\ No newline at end of property
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/Registry 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/Registry 2012-10-19 19:42:08 UTC (rev 2229)
@@ -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
@@ -60,7 +65,7 @@
namelist real hmix config_apvm_scale_factor 0.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
@@ -162,7 +167,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 - -
@@ -192,8 +197,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 - -
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_equation_of_state_jm.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2012-10-19 19:42:08 UTC (rev 2229)
@@ -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/option3_b4b_test/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-19 19:42:08 UTC (rev 2229)
@@ -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
@@ -352,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))
@@ -391,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))
@@ -532,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
@@ -541,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.
@@ -558,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
@@ -589,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
@@ -606,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
@@ -626,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
@@ -636,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
@@ -652,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
@@ -736,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.
@@ -748,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.
@@ -763,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
@@ -788,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 :: &
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_tendency.F 2012-10-19 19:42:08 UTC (rev 2229)
@@ -411,8 +411,7 @@
real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, &
- referenceBottomDepth, ssh
+ bottomDepth, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, ssh
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -435,11 +434,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
@@ -461,9 +460,8 @@
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
@@ -859,9 +857,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)
@@ -890,10 +888,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) &
@@ -910,13 +908,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
@@ -958,7 +954,7 @@
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
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_test_cases.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_test_cases.F 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_test_cases.F 2012-10-19 19:42:08 UTC (rev 2229)
@@ -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/option3_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_time_integration_split.F 2012-10-19 19:42:08 UTC (rev 2229)
@@ -419,9 +419,20 @@
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
@@ -436,11 +447,22 @@
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
@@ -547,7 +569,19 @@
+ config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
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: 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)) &
@@ -569,9 +603,19 @@
+ 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)) &
Modified: branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/core_ocean/mpas_ocn_vmix_coefs_tanh.F 2012-10-19 19:42:08 UTC (rev 2229)
@@ -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/option3_b4b_test/src/registry/gen_inc.c
===================================================================
--- branches/ocean_projects/option3_b4b_test/src/registry/gen_inc.c 2012-10-19 17:50:51 UTC (rev 2228)
+++ branches/ocean_projects/option3_b4b_test/src/registry/gen_inc.c 2012-10-19 19:42:08 UTC (rev 2229)
@@ -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>