<p><b>mpetersen@lanl.gov</b> 2012-10-11 15:59:17 -0600 (Thu, 11 Oct 2012)</p><p>branch commit, partial_bottom_cells: Add ability to mpas code to alter initial condition to run with partial bottom cells. Also changed zstarWeights to accomodate pbcs. For zstar, zstarWeights=1 for all layers, and the h weight is included in the sum of h_tend_col in subroutine ocn_wtop.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/partial_bottom_cells/namelist.input.ocean
===================================================================
--- branches/ocean_projects/partial_bottom_cells/namelist.input.ocean 2012-10-11 17:28:10 UTC (rev 2205)
+++ branches/ocean_projects/partial_bottom_cells/namelist.input.ocean 2012-10-11 21:59:17 UTC (rev 2206)
@@ -31,6 +31,10 @@
config_pressure_type = 'pressure'
config_rho0 = 1014.65
/
+&partial_bottom_cells
+ config_alter_ICs_for_pbc = .false.
+ config_min_pbc_fraction = 0.10
+/
&split_explicit_ts
config_n_ts_iter = 2
config_n_bcl_iter_beg = 1
Modified: branches/ocean_projects/partial_bottom_cells/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/partial_bottom_cells/src/core_ocean/Registry 2012-10-11 17:28:10 UTC (rev 2205)
+++ branches/ocean_projects/partial_bottom_cells/src/core_ocean/Registry 2012-10-11 21:59:17 UTC (rev 2206)
@@ -31,7 +31,8 @@
namelist character grid config_pressure_type pressure
namelist real grid config_rho0 1028
namelist logical grid config_enforce_zstar_at_restart false
-namelist logical grid config_set_bottomDepth_to_ref false
+namelist logical partial_bottom_cells config_alter_ICs_for_pbc false
+namelist real partial_bottom_cells config_min_pbc_fraction 0.10
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
Modified: branches/ocean_projects/partial_bottom_cells/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/partial_bottom_cells/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-11 17:28:10 UTC (rev 2205)
+++ branches/ocean_projects/partial_bottom_cells/src/core_ocean/mpas_ocn_mpas_core.F 2012-10-11 21:59:17 UTC (rev 2206)
@@ -106,37 +106,37 @@
if (.not. config_do_restart) call setup_sw_test_case(domain)
+ if (config_vert_grid_type.ne.'isopycnal') call ocn_init_z_level(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
call ocn_init_split_timestep(domain)
- print *, ' Vertical grid type is: ',config_vert_grid_type
+ write (0,'(a,a10)'), ' Vertical grid type is: ',config_vert_grid_type
if (config_vert_grid_type.ne.'isopycnal'.and. &
config_vert_grid_type.ne.'zlevel'.and. &
config_vert_grid_type.ne.'zstar1'.and. &
config_vert_grid_type.ne.'zstar'.and. &
config_vert_grid_type.ne.'zstarWeights') then
- print *, ' Incorrect choice of config_vert_grid_type.'
+ write (0,*) ' Incorrect choice of config_vert_grid_type.'
call mpas_dmpar_abort(dminfo)
endif
- print *, ' Pressure type is: ',config_pressure_type
+ write (0,'(a,a10)'), ' Pressure type is: ',config_pressure_type
if (config_pressure_type.ne.'pressure'.and. &
config_pressure_type.ne.'MontgomeryPotential') then
- print *, ' Incorrect choice of config_pressure_type.'
+ write (0,*) ' Incorrect choice of config_pressure_type.'
call mpas_dmpar_abort(dminfo)
endif
if (config_filter_btr_mode.and. &
config_vert_grid_type.ne.'zlevel')then
- print *, 'filter_btr_mode has only been tested with'// &
+ write (0,*) 'filter_btr_mode has only been tested with'// &
' config_vert_grid_type=zlevel.'
call mpas_dmpar_abort(dminfo)
endif
@@ -544,14 +544,16 @@
integer :: i, iCell, iEdge, iVertex, k, nCells
type (block_type), pointer :: block
- integer :: iTracer, cell, cell1, cell2
- real (kind=RKIND) :: uhSum, hSum, hEdge1
+ integer :: iTracer, cell, cell1, cell2, indexT, indexS
+ 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.
@@ -559,6 +561,7 @@
do while (associated(block))
h => block % state % time_levs(1) % state % h % 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
@@ -585,15 +588,6 @@
refBottomDepthTopOfCell(k+1) = refBottomDepth(k)
end do
- ! Some older grids do not have the bottomDepth variable.
- ! If so, set config_set_bottomDepth_to_ref=.true. and bottomDepth will
- ! simply be initialized with refBottomDepth. This does not use pbcs.
- if (config_set_bottomDepth_to_ref) then
- do iCell = 1,nCells
- bottomDepth(iCell) = refBottomDepth(maxLevelCell(iCell))
- enddo
- endif
-
! Initialization of zstarWeights. This determines how SSH perturbations
! are distributed throughout the column.
if (config_vert_grid_type.eq.'zlevel') then
@@ -603,9 +597,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
@@ -620,16 +612,89 @@
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=.true. 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=.false. 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 (config_alter_ICs_for_pbc.and..not.config_do_restart) then
+ write (0,'(a)') ' Altering bottomDepth to avoid very thin cells.'
+ write (0,'(a)') ' Altering h, temperature, and salinity initial conditions to conform with partial bottom cells.'
+
+ allocate(minBottomDepth(nVertLevels),minBottomDepthMid(nVertLevels),zMidZLevel(nVertLevels))
+
+ indexT = block % state % time_levs(1) % state % index_temperature
+ indexS = block % state % time_levs(1) % state % index_salinity
+
+ ! 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)
+
+! original. delete soon. This version always rounds down.
+! ! For pbcs, bottom depth must be altered to prevent very thin cells.
+! bottomDepth(iCell) = max(bottomDepth(iCell),minBottomDepth(k))
+
+ 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))
+ tracers(indexT,k,iCell) = tracers(indexT,k,iCell) &
+ + (tracers(indexT,k-1,iCell) - tracers(indexT,k,iCell)) &
+ /(zMidZLevel(k-1)-zMidZLevel(k)) &
+ *(zMidPBC - zMidZLevel(k))
+ tracers(indexS,k,iCell) = tracers(indexS,k,iCell) &
+ + (tracers(indexS,k-1,iCell) - tracers(indexS,k,iCell)) &
+ /(zMidZLevel(k-1)-zMidZLevel(k)) &
+ *(zMidPBC - zMidZLevel(k))
+
+ enddo
+
+ deallocate(minBottomDepth,zMidZLevel)
+
+ elseif (.not.config_alter_ICs_for_pbc.and..not.config_do_restart) then
+
+ do iCell = 1,nCells
+ bottomDepth(iCell) = refBottomDepth(maxLevelCell(iCell))
+ enddo
+
+ endif
+
if (config_vert_grid_type.eq.'zlevel'.or. &
config_vert_grid_type.eq.'zstar'.or. &
config_vert_grid_type.eq.'zstarweights') then
do iCell = 1,nCells
- ! mrp, for debugging. Delete later:
- write (0,'(a,2i5,10f10.3)') ' iCell, maxLevelCell, bdepth, refb, reft, sumh-bottomDepth: ', &
- iCell, maxLevelCell(iCell), bottomDepth(iCell), &
- refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell)), &
- sum(h(1:maxLevelCell(iCell),iCell))-bottomDepth(iCell)
+ ! mrp, for debugging.
+ !write (0,'(a,2i5,10f10.3)') ' iCell, maxLevelCell, bdepth, refb, reft, sumh-bottomDepth: ', &
+ ! iCell, maxLevelCell(iCell), bottomDepth(iCell), &
+ ! refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell)), &
+ ! sum(h(1:maxLevelCell(iCell),iCell))-bottomDepth(iCell)
! Check if abs(ssh)>2m. If so, print warning.
if (abs(sum(h(1:maxLevelCell(iCell),iCell))-bottomDepth(iCell))>2.0) then
@@ -642,10 +707,10 @@
! Check that bottomDepth and maxLevelCell match forSome older grids do not have the bottomDepth variable.
if (bottomDepth(iCell) > refBottomDepth(maxLevelCell(iCell)).or. &
bottomDepth(iCell) < refBottomDepthTopOfCell(maxLevelCell(iCell))) then
- write (0,*) ' fatal error: bottomDepth and maxLevelCell do not match:'
- write (0,*) ' iCell, maxLevelCell(iCell), bottomDepth(iCell): ', &
+ 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,*) ' refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell)): ', &
+ write (0,'(a,10f10.2)') ' refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell)): ', &
refBottomDepth(maxLevelCell(iCell)), refBottomDepthTopOfCell(maxLevelCell(iCell))
call mpas_dmpar_finalize(domain % dminfo)
endif
@@ -835,10 +900,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 :: &
</font>
</pre>