<p><b>mpetersen@lanl.gov</b> 2010-11-16 14:23:50 -0700 (Tue, 16 Nov 2010)</p><p>TRUNK COMMIT Ocean core: merge from branch topography2_mrp.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2010-11-16 17:49:07 UTC (rev 618)
+++ trunk/mpas/namelist.input.ocean 2010-11-16 21:23:50 UTC (rev 619)
@@ -41,6 +41,9 @@
config_vert_viscosity = 1.0e-4
config_vert_diffusion = 1.0e-4
/
+&eos
+ config_eos_type = 'linear'
+/
&advection
config_vert_tracer_adv = 'upwind'
config_tracer_adv_order = 2
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2010-11-16 17:49:07 UTC (rev 618)
+++ trunk/mpas/src/core_ocean/Registry 2010-11-16 21:23:50 UTC (rev 619)
@@ -30,6 +30,7 @@
namelist real vmix config_vmixTanhDiffMin 1.0e-5
namelist real vmix config_vmixTanhZMid -100
namelist real vmix config_vmixTanhZWidth 100
+namelist character eos config_eos_type linear
namelist character advection config_vert_tracer_adv 'centered'
namelist integer advection config_tracer_adv_order 2
namelist integer advection config_thickness_adv_order 2
@@ -106,7 +107,7 @@
var persistent real h_s ( nCells ) 0 iro h_s mesh - -
# Space needed for advection
-var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 o deriv_two mesh - -
+var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 - deriv_two mesh - -
var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
# !! NOTE: the following arrays are needed to allow the use
@@ -120,11 +121,14 @@
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
# Arrays for z-level version of mpas-ocean
-var persistent integer maxLevelsCell ( nCells ) 0 iro kmaxCell mesh - -
-var persistent integer maxLevelsEdge ( nEdges ) 0 iro kmaxEdge mesh - -
+var persistent integer maxLevelCell ( nCells ) 0 iro maxLevelCell mesh - -
+var persistent integer maxLevelEdgeTop ( nEdges ) 0 - maxLevelEdgeTop mesh - -
+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 hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
-var persistent real zMidZLevel ( nVertLevels ) 0 iro zMidZLevel mesh - -
-var persistent real zTopZLevel ( nVertLevelsP1 ) 0 iro zTopZLevel mesh - -
+var persistent real zMidZLevel ( nVertLevels ) 0 - zMidZLevel mesh - -
+var persistent real zTopZLevel ( nVertLevelsP1 ) 0 - zTopZLevel mesh - -
# Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
@@ -150,31 +154,24 @@
var persistent real tend_tracer2 ( nVertLevels nCells Time ) 1 - tracer2 tend tracers testing
# Diagnostic fields: only written to output
-var persistent real v ( nVertLevels nEdges Time ) 2 o v state - -
+var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
-var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
-var persistent real pv_edge ( nVertLevels nEdges Time ) 2 o pv_edge state - -
-var persistent real h_edge ( nVertLevels nEdges Time ) 2 o h_edge state - -
-var persistent real h_vertex ( nVertLevels nVertices Time ) 2 o h_vertex state - -
+var persistent real vorticity ( nVertLevels nVertices Time ) 2 - vorticity state - -
+var persistent real pv_edge ( nVertLevels nEdges Time ) 2 - pv_edge state - -
+var persistent real h_edge ( nVertLevels nEdges Time ) 2 - h_edge state - -
+var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
-var persistent real ke_edge ( nVertLevels nEdges Time ) 2 o ke_edge state - -
-var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 o pv_vertex state - -
+var persistent real ke_edge ( nVertLevels nEdges Time ) 2 - ke_edge state - -
+var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 - pv_vertex state - -
var persistent real pv_cell ( nVertLevels nCells Time ) 2 o pv_cell state - -
var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
var persistent real uReconstructY ( nVertLevels nCells Time ) 2 o uReconstructY state - -
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 o uReconstructZ state - -
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
-var persistent real zMid ( nVertLevels nCells Time ) 2 o zMid state - -
-var persistent real zTop ( nVertLevelsP1 nCells Time ) 2 o zTop state - -
-var persistent real zMidEdge ( nVertLevels nEdges Time ) 2 o zMidEdge state - -
-var persistent real zTopEdge ( nVertLevelsP1 nEdges Time ) 2 o zTopEdge state - -
-var persistent real p ( nVertLevels nCells Time ) 2 o p state - -
-var persistent real pTop ( nVertLevelsP1 nCells Time ) 2 o pTop state - -
-var persistent real pZLevel ( nVertLevels nCells Time ) 2 o pZLevel state - -
var persistent real MontPot ( nVertLevels nCells Time ) 2 o MontPot state - -
+var persistent real pressure ( nVertLevels nCells Time ) 2 o pressure state - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 o wTop state - -
-var persistent real ssh ( nCells Time ) 2 o ssh state - -
# Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
@@ -182,7 +179,6 @@
var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
-# xsad 10-02-05:
# Globally reduced diagnostic variables: only written to output
var persistent real areaCellGlobal ( Time ) 2 o areaCellGlobal state - -
var persistent real areaEdgeGlobal ( Time ) 2 o areaEdgeGlobal state - -
@@ -191,5 +187,4 @@
var persistent real volumeCellGlobal ( Time ) 2 o volumeCellGlobal state - -
var persistent real volumeEdgeGlobal ( Time ) 2 o volumeEdgeGlobal state - -
var persistent real CFLNumberGlobal ( Time ) 2 o CFLNumberGlobal state - -
-# xsad 10-02-05 end
Modified: trunk/mpas/src/core_ocean/module_global_diagnostics.F
===================================================================
--- trunk/mpas/src/core_ocean/module_global_diagnostics.F 2010-11-16 17:49:07 UTC (rev 618)
+++ trunk/mpas/src/core_ocean/module_global_diagnostics.F 2010-11-16 21:23:50 UTC (rev 619)
@@ -37,7 +37,7 @@
real (kind=RKIND) :: areaCellGlobal, areaEdgeGlobal, areaTriangleGlobal
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, areaEdge
real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, circulation, vorticity, ke, pv_edge, pv_vertex, &
- pv_cell, gradPVn, gradPVt, zMid, zTop, p, pTop, MontPot, wTop, rho, tracerTemp
+ pv_cell, gradPVn, gradPVt, pressure, MontPot, wTop, rho, tracerTemp
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
@@ -80,11 +80,8 @@
pv_cell => state % pv_cell % array
gradPVn => state % gradPVn % array
gradPVt => state % gradPVt % array
- zMid => state % zMid % array
- zTop => state % zTop % array
- p => state % p % array
- pTop => state % pTop % array
MontPot => state % MontPot % array
+ pressure => state % pressure % array
variableIndex = 0
! h
@@ -156,30 +153,12 @@
gradPVt(:,1:nEdgesSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! zMid
+ ! pressure
variableIndex = variableIndex + 1
- call computeFieldAreaWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), &
- zMid(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
- verticalSumMaxes(variableIndex))
-
- ! zTop
- variableIndex = variableIndex + 1
- call computeFieldAreaWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), &
- zTop(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
- verticalSumMaxes(variableIndex))
-
- ! p
- variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- p(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ pressure(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
- ! pTop
- variableIndex = variableIndex + 1
- call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
- pTop(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
- verticalSumMaxes(variableIndex))
-
! MontPot
variableIndex = variableIndex + 1
call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
@@ -309,22 +288,10 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeEdgeGlobal
- ! zMid
+ ! pressure
variableIndex = variableIndex + 1
- averages(variableIndex) = sums(variableIndex)/(areaCellGlobal*nVertLevels)
-
- ! zTop
- variableIndex = variableIndex + 1
- averages(variableIndex) = sums(variableIndex)/(areaCellGlobal*nVertLevels)
-
- ! p
- variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
- ! pTop
- variableIndex = variableIndex + 1
- averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
-
! MontPot
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
Modified: trunk/mpas/src/core_ocean/module_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/module_mpas_core.F 2010-11-16 17:49:07 UTC (rev 618)
+++ trunk/mpas/src/core_ocean/module_mpas_core.F 2010-11-16 21:23:50 UTC (rev 619)
@@ -1,6 +1,8 @@
module mpas_core
use mpas_framework
+ use dmpar
+ use test_cases
type (io_output_object) :: restart_obj
integer :: restart_frame
@@ -13,7 +15,6 @@
use configure
use grid_types
- use test_cases
implicit none
@@ -21,10 +22,24 @@
real (kind=RKIND) :: dt
type (block_type), pointer :: block
+ type (dm_info) :: dminfo
if (.not. config_do_restart) call setup_sw_test_case(domain)
+ if (config_vert_grid_type.eq.'isopycnal') then
+ print *, ' Using isopycnal coordinates'
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ print *, ' Using z-level coordinates'
+ call init_ZLevel(domain)
+ else
+ print *, ' Incorrect choice of config_vert_grid_type:',&
+ config_vert_grid_type
+ call dmpar_abort(dminfo)
+ endif
+
+ call compute_maxLevel(domain)
+
!
! Initialize core
!
@@ -208,6 +223,180 @@
end if
end subroutine mpas_timestep
+
+
+subroutine init_ZLevel(domain)
+! Initialize maxLevel and bouncary grid variables.
+
+ use grid_types
+ use configure
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: i, iCell, iEdge, iVertex, k
+ type (block_type), pointer :: block
+
+ integer :: iTracer, cell
+ real (kind=RKIND), dimension(:), pointer :: &
+ hZLevel, zMidZLevel, zTopZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: h
+ integer :: nVertLevels
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ h => block % state % time_levs(1) % state % h % array
+ hZLevel => block % mesh % hZLevel % array
+ zMidZLevel => block % mesh % zMidZLevel % array
+ zTopZLevel => block % mesh % zTopZLevel % array
+ nVertLevels = block % mesh % nVertLevels
+
+ ! These should eventually be in an input file. For now
+ ! I just read them in from h(:,1).
+ ! Upon restart, the correct hZLevel should be in restart.nc
+ if (.not. config_do_restart) hZLevel = h(:,1)
+
+ ! hZLevel should be in the grid.nc and restart.nc file,
+ ! and h for k=1 must be specified there as well.
+
+ zTopZLevel(1) = 0.0
+ do k = 1,nVertLevels
+ zMidZLevel(k) = zTopZLevel(k)-0.5*hZLevel(k)
+ zTopZLevel(k+1) = zTopZLevel(k)- hZLevel(k)
+ end do
+
+ block => block % next
+
+ end do
+
+end subroutine init_ZLevel
+
+
+subroutine compute_maxLevel(domain)
+! Initialize maxLevel and bouncary grid variables.
+
+ use grid_types
+ use configure
+ use constants
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ 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 :: &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, boundaryEdge, boundaryCell, &
+ boundaryVertex, verticesOnEdge
+
+ ! Initialize z-level grid variables from h, read in from input file.
+ block => domain % blocklist
+ do while (associated(block))
+
+ maxLevelCell => block % mesh % maxLevelCell % array
+ maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+ maxLevelEdgeBot => block % mesh % maxLevelEdgeBot % array
+ maxLevelVertexTop => block % mesh % maxLevelVertexTop % array
+ maxLevelVertexBot => block % mesh % maxLevelVertexBot % array
+ cellsOnEdge => block % mesh % cellsOnEdge % array
+ cellsOnVertex => block % mesh % cellsOnVertex % array
+ verticesOnEdge => block % mesh % verticesOnEdge % array
+ boundaryEdge => block % mesh % boundaryEdge % array
+ boundaryCell => block % mesh % boundaryCell % array
+ boundaryVertex => block % mesh % boundaryVertex % array
+
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertices = block % mesh % nVertices
+ nVertLevels = block % mesh % nVertLevels
+ vertexDegree = block % mesh % vertexDegree
+
+ ! for z-grids, maxLevelCell should be in input state
+ ! Isopycnal grid uses all vertical cells
+ if (config_vert_grid_type.eq.'isopycnal') then
+ maxLevelCell(1:nCells) = nVertLevels
+ endif
+ maxLevelCell(nCells+1) = 0
+
+ ! maxLevelEdgeTop is the minimum (shallowest) of the surrounding cells
+ do iEdge=1,nEdges
+ maxLevelEdgeTop(iEdge) = &
+ min( maxLevelCell(cellsOnEdge(1,iEdge)), &
+ maxLevelCell(cellsOnEdge(2,iEdge)) )
+ end do
+ maxLevelEdgeTop(nEdges+1) = 0
+
+ ! maxLevelEdgeBot is the maximum (deepest) of the surrounding cells
+ do iEdge=1,nEdges
+ maxLevelEdgeBot(iEdge) = &
+ max( maxLevelCell(cellsOnEdge(1,iEdge)), &
+ maxLevelCell(cellsOnEdge(2,iEdge)) )
+ end do
+ maxLevelEdgeBot(nEdges+1) = 0
+
+ ! maxLevelVertexBot is the maximum (deepest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexBot(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
+ do i=2,vertexDegree
+ maxLevelVertexBot(iVertex) = &
+ max( maxLevelVertexBot(iVertex), &
+ maxLevelCell(cellsOnVertex(i,iVertex)))
+ end do
+ end do
+ maxLevelVertexBot(nVertices+1) = 0
+
+ ! maxLevelVertexTop is the minimum (shallowest) of the surrounding cells
+ do iVertex = 1,nVertices
+ maxLevelVertexTop(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex))
+ do i=2,vertexDegree
+ maxLevelVertexTop(iVertex) = &
+ min( maxLevelVertexTop(iVertex), &
+ maxLevelCell(cellsOnVertex(i,iVertex)))
+ end do
+ end do
+ maxLevelVertexTop(nVertices+1) = 0
+
+ ! set boundary edge
+ boundaryEdge=1
+ do iEdge=1,nEdges
+ boundaryEdge(1:maxLevelEdgeTop(iEdge),iEdge)=0
+ end do
+
+ !
+ ! Find cells and vertices that have an edge on the boundary
+ !
+ boundaryCell(:,:) = 0
+ do iEdge=1,nEdges
+ do k=1,nVertLevels
+ if (boundaryEdge(k,iEdge).eq.1) then
+ boundaryCell(k,cellsOnEdge(1,iEdge)) = 1
+ boundaryCell(k,cellsOnEdge(2,iEdge)) = 1
+ boundaryVertex(k,verticesOnEdge(1,iEdge)) = 1
+ boundaryVertex(k,verticesOnEdge(2,iEdge)) = 1
+ endif
+ end do
+ end do
+
+ block => block % next
+ end do
+
+ ! Note: We do not update halos on maxLevel* variables. I want the
+ ! outside edge of a halo to be zero on each processor.
+
+end subroutine compute_maxLevel
subroutine mpas_core_finalize(domain)
Modified: trunk/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- trunk/mpas/src/core_ocean/module_test_cases.F 2010-11-16 17:49:07 UTC (rev 618)
+++ trunk/mpas/src/core_ocean/module_test_cases.F 2010-11-16 21:23:50 UTC (rev 619)
@@ -25,15 +25,6 @@
type (block_type), pointer :: block_ptr
type (dm_info) :: dminfo
- integer :: iTracer
- real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
- hZLevel, zMidZLevel, zTopZLevel, latCell,LonCell
- 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
-
if (config_test_case == 0) then
write(0,*) 'Using initial conditions supplied in input file'
@@ -95,123 +86,6 @@
block_ptr => block_ptr % next
end do
- ! Initialize z-level grid variables from h, read in from input file.
- block_ptr => domain % blocklist
- do while (associated(block_ptr))
- h => block_ptr % state % time_levs(1) % state % h % array
- u => block_ptr % state % time_levs(1) % state % u % array
- rho => block_ptr % state % time_levs(1) % state % rho % array
- tracers => block_ptr % state % time_levs(1) % state % tracers % array
- u_src => block_ptr % mesh % u_src % array
- xCell => block_ptr % mesh % xCell % array
- yCell => block_ptr % mesh % yCell % array
- latCell => block_ptr % mesh % latCell % array
- lonCell => block_ptr % mesh % lonCell % array
-
- hZLevel => block_ptr % mesh % hZLevel % array
- zMidZLevel => block_ptr % mesh % zMidZLevel % array
- zTopZLevel => block_ptr % mesh % zTopZLevel % array
-
- nCells = block_ptr % mesh % nCells
- nEdges = block_ptr % mesh % nEdges
- nVertices = block_ptr % mesh % nVertices
- nVertLevels = block_ptr % mesh % nVertLevels
-
- pi=3.1415
- ! Tracer blob in Central Pacific, away from boundaries:
- !latCenter=pi/16; lonCenter=9./8.*pi
-
- ! Tracer blob in Central Pacific, near boundaries:
- latCenter=pi*2./16; lonCenter=13./16.*pi
-
- if (config_vert_grid_type.eq.'zlevel') then
- ! These should eventually be in an input file. For now
- ! I just read them in from h(:,1).
- hZLevel = h(:,1)
- zTopZLevel(1) = 0.0
- do iLevel = 1,nVertLevels
- zMidZLevel(iLevel) = zTopZLevel(iLevel)-0.5*hZLevel(iLevel)
- zTopZLevel(iLevel+1) = zTopZLevel(iLevel)- hZLevel(iLevel)
- enddo
- if (config_vert_grid_type.eq.'isopycnal') then
- print *, ' Using isopycnal coordinates'
- elseif (config_vert_grid_type.eq.'zlevel') then
- print *, ' Using z-level coordinates'
- else
- print *, ' Incorrect choice of config_vert_grid_type:',&
- config_vert_grid_type
- call dmpar_abort(dminfo)
- endif
-
- ! Set tracers, if not done in grid.nc file
- !tracers = 0.0
- centerx = 1.0e6
- centery = 1.0e6
- dist = 2.0e5
- do iCell = 1,nCells
- dist = sqrt( (latCell(iCell)-latCenter)**2 + (lonCell(iCell)-lonCenter)**2)
- do iLevel = 1,nVertLevels
- ! for 20 layer test
- ! tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) = 5.0 ! temperature
- ! tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
-
- ! for x3, 25 layer test
- !tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) = 10.0 ! temperature
- !tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
-
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = &
- ! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
-
- !tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = block_ptr % mesh % lonCell % array(iCell)
- !tracers(block_ptr % state % time_levs(1) % state % index_tracer2,iLevel,iCell) = block_ptr % mesh % latCell % array(iCell)
-
- ! place tracer blob with radius dist at (centerx,centery)
- !if ( sqrt( (xCell(iCell)-centerx)**2 &
- ! + (yCell(iCell)-centery)**2) &
- ! .lt.dist) then
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 1.0
- !else
- ! tracers(block_ptr % state % time_levs(1) % state % index_tracer1,iLevel,iCell) = 0.0
- !endif
-
- rho(iLevel,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(block_ptr % state % time_levs(1) % state % index_temperature,iLevel,iCell) &
- + 7.6e-4*tracers(block_ptr % state % time_levs(1) % state % index_salinity,iLevel,iCell))
-
- enddo
- enddo
-
- endif
-
- ! print some diagnostics
- print '(10a)', 'ilevel',&
- ' rho ',&
- ' min u max u ',&
- ' min u_src max u_src ', &
- ' min h max h ',&
- ' hZLevel zMidZlevel', &
- ' zTopZlevel'
- do iLevel = 1,nVertLevels
- print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
- minval(u(iLevel,1:nEdges)), maxval(u(iLevel,1:nEdges)), &
- minval(u_src(iLevel,1:nEdges)), maxval(u_src(iLevel,1:nEdges)), &
- minval(h(iLevel,1:nCells)), maxval(h(iLevel,1:nCells)), &
- hZLevel(iLevel),zMidZlevel(iLevel),zTopZlevel(iLevel)
- enddo
-
- print '(10a)', 'itracer ilevel min tracer max tracer'
- do iTracer=1,block_ptr % state % time_levs(1) % state % num_tracers
- do iLevel = 1,nVertLevels
- print '(2i5,20es12.4)', iTracer,ilevel, &
- minval(tracers(itracer,iLevel,1:nCells)), maxval(tracers(itracer,iLevel,1:nCells))
- enddo
- enddo
-
- block_ptr => block_ptr % next
- end do
-
-
end subroutine setup_sw_test_case
Modified: trunk/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_ocean/module_time_integration.F 2010-11-16 17:49:07 UTC (rev 618)
+++ trunk/mpas/src/core_ocean/module_time_integration.F 2010-11-16 21:23:50 UTC (rev 619)
@@ -92,7 +92,7 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
* block % state % time_levs(1) % state % h % array(k,iCell)
end do
@@ -175,7 +175,7 @@
provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
+ rk_substep_weights(rk_step) * block % tend % h % array(:,:)
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
provis % tracers % array(:,k,iCell) = ( &
block % state % time_levs(1) % state % h % array(k,iCell) * &
block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
@@ -204,7 +204,7 @@
+ rk_weights(rk_step) * block % tend % h % array(:,:)
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
@@ -225,7 +225,7 @@
block => domain % blocklist
do while (associated(block))
do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
+ do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
/ block % state % time_levs(2) % state % h % array(k,iCell)
@@ -264,22 +264,23 @@
type (state_type), intent(in) :: s
type (mesh_type), intent(in) :: grid
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j, kmax
- integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND) :: h_mom_eddy_visc2, h_mom_eddy_visc4
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
upstream_bias, wTopEdge, rho0Inv, r
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wTop, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
tend_h, tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- divergence, MontPot, pZLevel, zMidEdge, zTopEdge
+ MontPot, wTop, divergence
type (dm_info) :: dminfo
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: &
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
edgesOnEdge, edgesOnVertex
@@ -303,12 +304,10 @@
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
- ke_edge => s % ke_edge % array
+ ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
MontPot => s % MontPot % array
- pZLevel => s % pZLevel % array
- zTopEdge => s % zTopEdge % array
- zMidEdge => s % zMidEdge % array
+ pressure => s % pressure % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -329,44 +328,52 @@
fEdge => grid % fEdge % array
zMidZLevel => grid % zMidZLevel % array
zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
tend_h => tend % h % array
tend_u => tend % u % array
nCells = grid % nCells
nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
u_src => grid % u_src % array
- h_mom_eddy_visc2 = config_h_mom_eddy_visc2
- h_mom_eddy_visc4 = config_h_mom_eddy_visc4
+ !
+ ! height tendency: start accumulating tendency terms
+ !
+ tend_h = 0.0
!
! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
!
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
! for explanation of divergence operator.
- tend_h(:,:) = 0.0
+ !
+ ! for z-level, only compute height tendency for top layer.
+
+ if (config_vert_grid_type.eq.'isopycnal') then
+ kmax = nVertLevels
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ kmax = 1
+ endif
+
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell1) = tend_h(k,cell1) - flux
- end do
- end if
- if (cell2 <= nCells) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell2) = tend_h(k,cell2) + flux
- end do
- end if
- end do
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,kmax
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell2) = tend_h(k,cell2) + flux
+ end do
+ end do
+
do iCell=1,nCells
- do k=1,nVertLevels
+ do k=1,kmax
tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
end do
end do
@@ -374,22 +381,10 @@
!
! height tendency: vertical advection term -d/dz(hw)
!
+ ! Vertical advection computed for top layer of a z grid only.
if (config_vert_grid_type.eq.'zlevel') then
-
do iCell=1,nCells
-
tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
-
- ! This next loop is to verify that h for levels below the first
- ! remain constant. At a later time this could be replaced
- ! by just tend_h(2:nVertLevels,:) = 0.0, and then there is
- ! no need to compute the horizontal tend_h term for k=2:nVertLevels
- ! on a zlevel grid, above.
- do k=2,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) &
- - (wTop(k,iCell) - wTop(k+1,iCell))
- end do
-
end do
endif ! coordinate type
@@ -401,30 +396,30 @@
!
! velocity tendency: vertical advection term -w du/dz
!
- allocate(w_dudzTopEdge(nVertLevels+1))
- w_dudzTopEdge(1) = 0.0
- w_dudzTopEdge(nVertLevels+1) = 0.0
if (config_vert_grid_type.eq.'zlevel') then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ allocate(w_dudzTopEdge(nVertLevels+1))
+ w_dudzTopEdge(1) = 0.0
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- do k=2,nVertLevels
- ! Average w from cell center to edge
- wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
+ do k=2,maxLevelEdgeTop(iEdge)
+ ! Average w from cell center to edge
+ wTopEdge = 0.5*(wTop(k,cell1)+wTop(k,cell2))
- ! compute dudz at vertical interface with first order derivative.
- w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
- / (zMidZLevel(k-1) - zMidZLevel(k))
- end do
+ ! compute dudz at vertical interface with first order derivative.
+ w_dudzTopEdge(k) = wTopEdge * (u(k-1,iEdge)-u(k,iEdge)) &
+ / (zMidZLevel(k-1) - zMidZLevel(k))
+ end do
+ w_dudzTopEdge(maxLevelEdgeTop(iEdge)+1) = 0.0
- ! Average w*du/dz from vertical interface to vertical middle of cell
- do k=1,nVertLevels
- tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
- enddo
- enddo
+ ! Average w*du/dz from vertical interface to vertical middle of cell
+ do k=1,maxLevelEdgeTop(iEdge)
+ tend_u(k,iEdge) = - 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
+ enddo
+ enddo
+ deallocate(w_dudzTopEdge)
endif
- deallocate(w_dudzTopEdge)
!
! velocity tendency: pressure gradient
@@ -434,7 +429,7 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
end do
@@ -443,10 +438,10 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - rho0Inv*( pZLevel(k,cell2) &
- - pZLevel(k,cell1) )/dcEdge(iEdge)
+ - rho0Inv*( pressure(k,cell2) &
+ - pressure(k,cell1) )/dcEdge(iEdge)
end do
enddo
endif
@@ -454,16 +449,16 @@
!
! velocity tendency: del2 dissipation, </font>
<font color="black">u_2 </font>
<font color="black">abla^2 u
! computed as </font>
<font color="black">u( </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity )
- ! strictly only valid for h_mom_eddy_visc2 == constant
+ ! strictly only valid for config_h_mom_eddy_visc2 == constant
!
- if ( h_mom_eddy_visc2 > 0.0 ) then
+ if ( config_h_mom_eddy_visc2 > 0.0 ) then
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
! Here -( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
! is - </font>
<font color="gray">abla vorticity pointing from vertex 2 to vertex 1, or equivalently
@@ -471,7 +466,7 @@
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-( vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = h_mom_eddy_visc2 * u_diffusion
+ u_diffusion = config_h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
@@ -483,9 +478,9 @@
! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="black">abla^4 u
! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity
! applied recursively.
- ! strictly only valid for h_mom_eddy_visc4 == constant
+ ! strictly only valid for config_h_mom_eddy_visc4 == constant
!
- if ( h_mom_eddy_visc4 > 0.0 ) then
+ if ( config_h_mom_eddy_visc4 > 0.0 ) then
allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_u(nVertLevels, nEdges+1))
@@ -501,12 +496,11 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
-
- delsq_u(k,iEdge) = delsq_u(k,iEdge) + u_diffusion
+ delsq_u(k,iEdge) = &
+ ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
end do
end do
@@ -516,7 +510,7 @@
do iEdge=1,nEdges
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
- dcEdge(iEdge) * delsq_u(k,iEdge)
delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
@@ -525,7 +519,7 @@
end do
do iVertex=1,nVertices
r = 1.0 / areaTriangle(iVertex)
- do k=1,nVertLevels
+ do k=1,maxLevelVertexBot(iVertex)
delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
end do
end do
@@ -535,7 +529,7 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+ delsq_u(k,iEdge)*dvEdge(iEdge)
delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
@@ -544,7 +538,7 @@
end do
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
+ do k = 1,maxLevelCell(iCell)
delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
end do
end do
@@ -557,14 +551,14 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
u_diffusion = ( delsq_divergence(k,cell2) &
- delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
-( delsq_vorticity(k,vertex2) &
- delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) - h_mom_eddy_visc4 * u_diffusion
+ tend_u(k,iEdge) = tend_u(k,iEdge) - config_h_mom_eddy_visc4 * u_diffusion
end do
end do
@@ -582,7 +576,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
@@ -600,24 +594,39 @@
!
! velocity tendency: forcing and bottom drag
!
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! know the bottom edge with nonzero velocity and place the drag there.
+
do iEdge=1,grid % nEdgesSolve
- ! forcing in top layer only
- tend_u(1,iEdge) = tend_u(1,iEdge) &
+ ! forcing in top layer only
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+ k = maxLevelEdgeTop(iEdge)
+
+ ! efficiency note: it would be nice to avoid this
+ ! if within a do. This could be done with
+ ! k = max(maxLevelEdgeTop(iEdge),1)
+ ! and then tend_u(1,iEdge) is just not used for land cells.
+
+ if (k>0) then
+
! bottom drag is the same as POP:
! -c |u| u where c is unitless and 1.0e-3.
! see POP Reference guide, section 3.4.4.
- tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- - 1.0e-3*u(nVertLevels,iEdge) &
- *sqrt(2.0*ke_edge(nVertLevels,iEdge))/h_edge(nVertLevels,iEdge)
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ - 1.0e-3*u(k,iEdge) &
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+
! old bottom drag, just linear friction
! du/dt = u/tau where tau=100 days.
- !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- ! - u(nVertLevels,iEdge)/(100.0*86400.0)
+ !tend_u(k,iEdge) = tend_u(k,iEdge) &
+ ! - u(k,iEdge)/(100.0*86400.0)
+ endif
+
enddo
!
@@ -644,20 +653,23 @@
call dmpar_abort(dminfo)
endif
- allocate(fluxVertTop(1:nVertLevels+1))
+ allocate(fluxVertTop(nVertLevels+1))
fluxVertTop(1) = 0.0
- fluxVertTop(nVertLevels+1) = 0.0
do iEdge=1,grid % nEdgesSolve
- do k=2,nVertLevels
+
+ do k=2,maxLevelEdgeTop(iEdge)
fluxVertTop(k) = vertViscTop(k) &
* ( u(k-1,iEdge) - u(k,iEdge) ) &
- / (zMidEdge(k-1,iEdge) - zMidEdge(k,iEdge))
+ * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
enddo
- do k=1,nVertLevels
+ fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
+
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
+ (fluxVertTop(k) - fluxVertTop(k+1)) &
- /(zTopEdge(k,iEdge) - zTopEdge(k+1,iEdge))
+ / h_edge(k,iEdge)
enddo
+
end do
deallocate(fluxVertTop, vertViscTop)
@@ -682,17 +694,18 @@
integer :: i, k, iCell, iEdge, iTracer, cell1, cell2, upwindCell,&
nEdges, nCells, nVertLevels, num_tracers
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
- real (kind=RKIND) :: flux, tracer_edge, r, dist
+ real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,wTop, h_edge, zMid, zTop
+ u,h,wTop, h_edge
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:,:), pointer :: boundaryEdge
type (dm_info) :: dminfo
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
real (kind=RKIND), dimension(:), pointer :: zTopZLevel
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop, boundaryMask
@@ -711,8 +724,6 @@
wTop => s % wTop % array
tracers => s % tracers % array
h_edge => s % h_edge % array
- zMid => s % zMid % array
- zTop => s % zTop % array
tend_tr => tend % tracers % array
@@ -722,6 +733,9 @@
dcEdge => grid % dcEdge % array
zTopZLevel => grid % zTopZLevel % array
boundaryEdge => grid % boundaryEdge % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
nEdges = grid % nEdges
nCells = grid % nCells
@@ -738,6 +752,11 @@
!
! tracer tendency: horizontal advection term -div( h \phi u)
!
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the compute_solve_diagnostics
+ ! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
+ ! tracer_edge at the boundary will also need to be defined for flux boundaries.
+
coef_3rd_order = 0.
if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
@@ -747,16 +766,14 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,nVertLevels
- do iTracer=1,num_tracers
- tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- end do
+ do k=1,maxLevelEdgeTop(iEdge)
+ do iTracer=1,num_tracers
+ tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
+ flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
end do
- end if
+ end do
end do
else if (config_tracer_adv_order == 3) then
@@ -765,56 +782,52 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,nVertLevels
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ do iTracer=1,num_tracers
- do iTracer=1,num_tracers
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ !-- else u <= 0:
+ else
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ end if
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- !-- else u <= 0:
- else
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- end if
-
- !-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
+ !-- update tendency
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
end do
else if (config_tracer_adv_order == 4) then
@@ -823,46 +836,42 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if an edge is not on the outer-most ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,nVertLevels
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ do iTracer=1,num_tracers
- do iTracer=1,num_tracers
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
+ 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- !-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
+ !-- update tendency
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ enddo
+ end do
end do
endif ! if (config_tracer_adv_order == 2 )
@@ -873,11 +882,10 @@
!
allocate(tracerTop(num_tracers,nVertLevels+1))
tracerTop(:,1)=0.0
- tracerTop(:,nVertLevels+1)=0.0
do iCell=1,grid % nCellsSolve
if (config_vert_tracer_adv.eq.'centered') then
- do k=2,nVertLevels
+ do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
tracerTop(iTracer,k) = ( tracers(iTracer,k-1,iCell) &
+tracers(iTracer,k ,iCell))/2.0
@@ -885,7 +893,7 @@
end do
elseif (config_vert_tracer_adv.eq.'upwind') then
- do k=2,nVertLevels
+ do k=2,maxLevelCell(iCell)
if (wTop(k,iCell)>0.0) then
upwindCell = k
else
@@ -897,8 +905,9 @@
end do
endif
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
@@ -927,7 +936,7 @@
invAreaCell1 = 1.0/areaCell(cell1)
invAreaCell2 = 1.0/areaCell(cell2)
- do k=1,grid % nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
! \kappa_2 </font>
<font color="gray">abla \phi on edge
tracer_turb_flux = config_h_tracer_eddy_diff2 &
@@ -970,7 +979,7 @@
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- do k=1,grid % nVertLevels
+ do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ dvEdge(iEdge)*h_edge(k,iEdge) &
@@ -985,9 +994,9 @@
end do
- do iCell = 1, nCells
+ do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
end do
@@ -1001,21 +1010,20 @@
invAreaCell1 = 1.0 / areaCell(cell1)
invAreaCell2 = 1.0 / areaCell(cell2)
- do k=1,grid % nVertLevels
- do iTracer=1,num_tracers
- tracer_turb_flux = config_h_tracer_eddy_diff4 &
- *( delsq_tracer(iTracer,k,cell2) &
- - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
- flux = dvEdge (iEdge) * tracer_turb_flux
+ do k=1,maxLevelEdgeTop(iEdge)
+ do iTracer=1,num_tracers
+ tracer_turb_flux = config_h_tracer_eddy_diff4 &
+ *( delsq_tracer(iTracer,k,cell2) &
+ - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
+ flux = dvEdge (iEdge) * tracer_turb_flux
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) &
- - flux * invAreaCell1 * boundaryMask(k,iEdge)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) &
- + flux * invAreaCell2 * boundaryMask(k,iEdge)
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) &
+ - flux * invAreaCell1 * boundaryMask(k,iEdge)
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) &
+ + flux * invAreaCell2 * boundaryMask(k,iEdge)
- end do
+ enddo
enddo
-
end do
deallocate(delsq_tracer)
@@ -1048,29 +1056,30 @@
allocate(fluxVertTop(num_tracers,nVertLevels+1))
fluxVertTop(:,1) = 0.0
- fluxVertTop(:,nVertLevels+1) = 0.0
do iCell=1,grid % nCellsSolve
- do k=2,nVertLevels
+
+ do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
! compute \kappa_v d\phi/dz
fluxVertTop(iTracer,k) = vertDiffTop(k) &
* (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell) )&
- / (zMid(k-1,iCell) -zMid(k,iCell))
+ * 2 / (h(k-1,iCell) + h(k,iCell))
enddo
enddo
+ fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
- do k=1,nVertLevels
- dist = zTop(k,iCell) - zTop(k+1,iCell)
+ do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
+ ! This is h d/dz( fluxVertTop) but h and dz cancel, so
+ ! reduces to delta( fluxVertTop)
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- + h(k,iCell)*(fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1))/dist
+ + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
enddo
enddo
enddo ! iCell loop
deallocate(fluxVertTop, vertDiffTop)
-
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
! 'iTracer,k, minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))'
@@ -1102,25 +1111,30 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
- integer :: nCells, nEdges, nVertices, nVertLevels
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel, ssh
+ hZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wTop, &
- circulation, vorticity, ke, ke_edge, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
+ circulation, vorticity, ke, ke_edge, MontPot, wTop, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- zMid, zTop, zMidEdge, zTopEdge, p, pTop, MontPot, rho, temperature, salinity, pZLevel
+ rho, temperature, salinity
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ real (kind=RKIND), dimension(:), allocatable:: pTop
real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
+ verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
+ boundaryEdge, boundaryCell
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot, maxLevelVertexTop
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
@@ -1144,15 +1158,8 @@
gradPVt => s % gradPVt % array
rho => s % rho % array
tracers => s % tracers % array
- zMid => s % zMid % array
- zTop => s % zTop % array
- zMidEdge => s % zMidEdge % array
- zTopEdge => s % zTopEdge % array
- p => s % p % array
- pZLevel => s % pZLevel % array
- pTop => s % pTop % array
MontPot => s % MontPot % array
- ssh => s % ssh % array
+ pressure => s % pressure % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -1173,35 +1180,27 @@
fEdge => grid % fEdge % array
hZLevel => grid % hZLevel % array
deriv_two => grid % deriv_two % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+ maxLevelVertexTop => grid % maxLevelVertexTop % array
nCells = grid % nCells
nEdges = grid % nEdges
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ vertexDegree = grid % vertexDegree
boundaryEdge => grid % boundaryEdge % array
boundaryCell => grid % boundaryCell % array
!
- ! Find those cells that have an edge on the boundary
- !
- boundaryCell(:,:) = 0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if(boundaryEdge(k,iEdge).eq.1) then
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- boundaryCell(k,cell1) = 1
- boundaryCell(k,cell2) = 1
- endif
- enddo
- enddo
-
-
- !
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
!
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
coef_3rd_order = 0.
if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
@@ -1209,141 +1208,129 @@
if (config_thickness_adv_order == 2) then
- do iEdge=1,grid % nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end if
+ do k=1,maxLevelEdgeTop(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ end do
end do
else if (config_thickness_adv_order == 3) then
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,grid % nVertLevels
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ !-- if u > 0:
+ if (u(k,iEdge) > 0) then
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ !-- else u <= 0:
+ else
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ end if
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- !-- else u <= 0:
- else
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end if
-
- end do ! do k
- end if ! if (cell1 <=
+ end do ! do k
end do ! do iEdge
else if (config_thickness_adv_order == 4) then
- do iEdge=1,grid%nEdges
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,maxLevelEdgeTop(iEdge)
- do k=1,grid % nVertLevels
+ d2fdx2_cell1 = 0.0
+ d2fdx2_cell2 = 0.0
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
+ !-- if not a boundary cell
+ if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
+ !-- all edges of cell 1
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
+ end do
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
+ !-- all edges of cell 2
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
+ end do
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
+ endif
- endif
+ h_edge(k,iEdge) = &
+ 0.5*(h(k,cell1) + h(k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- end do ! do k
- end if ! if (cell1 <=
+ end do ! do k
end do ! do iEdge
endif ! if(config_thickness_adv_order == 2)
!
- ! set the velocity in the nEdges+1 slot to zero, this is a dummy address
- ! used to when reading for edges that do not exist
+ ! set the velocity and height at dummy address
+ ! used -1e34 so error clearly occurs if these values are used.
!
- u(:,nEdges+1) = 0.0
+ u(:,nEdges+1) = -1e34
+ h(:,nCells+1) = -1e34
+ tracers(s % index_temperature,:,nCells+1) = -1e34
+ tracers(s % index_salinity,:,nCells+1) = -1e34
!
! Compute circulation and relative vorticity at each vertex
!
circulation(:,:) = 0.0
do iEdge=1,nEdges
- if (verticesOnEdge(1,iEdge) <= nVertices) then
- do k=1,nVertLevels
- circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
- if (verticesOnEdge(2,iEdge) <= nVertices) then
- do k=1,nVertLevels
- circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
- end do
- end if
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ circulation(k,vertex1) = circulation(k,vertex1) - dcEdge(iEdge) * u(k,iEdge)
+ circulation(k,vertex2) = circulation(k,vertex2) + dcEdge(iEdge) * u(k,iEdge)
+ end do
end do
do iVertex=1,nVertices
- do k=1,nVertLevels
+ do k=1,maxLevelVertexBot(iVertex)
vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
end do
end do
-
!
! Compute the divergence at each cell center
!
@@ -1351,22 +1338,16 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- enddo
- endif
- if(cell2 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- enddo
- end if
+ do k=1,maxLevelEdgeBot(iEdge)
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ enddo
end do
do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
- divergence(k,iCell) = divergence(k,iCell) * r
- enddo
+ r = 1.0 / areaCell(iCell)
+ do k = 1,maxLevelCell(iCell)
+ divergence(k,iCell) = divergence(k,iCell) * r
+ enddo
enddo
!
@@ -1376,234 +1357,190 @@
do iCell=1,nCells
do i=1,nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
end do
end do
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
end do
end do
!
- !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- if (eoe <= nEdges) then
- do k = 1,nVertLevels
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
- end do
- end if
+ do k = 1,maxLevelEdgeBot(iEdge)
+ v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
+ end do
end do
end do
!
! Compute ke on cell edges at velocity locations for quadratic bottom drag.
!
+ ! mrp 101025 efficiency note: we could get rid of ke_edge completely by
+ ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
+ ke_edge = 0.0 !mrp remove 0 for efficiency
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
- do k=1,nVertLevels
- ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
- end do
- else
- do k=1,nVertLevels
- ke_edge(k,iEdge) = 0.0
- end do
- end if
+ do k=1,maxLevelEdgeTop(iEdge)
+ ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
+ end do
end do
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
- VTX_LOOP: do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
- if (cellsOnVertex(i,iVertex) > nCells) cycle VTX_LOOP
- end do
- do k=1,nVertLevels
+ do iVertex = 1,nVertices
+ do k=1,maxLevelVertexBot(iVertex)
h_vertex = 0.0
- do i=1,grid % vertexDegree
+ do i=1,vertexDegree
h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
h_vertex = h_vertex / areaTriangle(iVertex)
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
- end do VTX_LOOP
+ end do
-
!
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
+ ! Compute pv at cell centers
+ ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
!
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
- dvEdge(iEdge)
+ pv_cell(:,:) = 0.0
+ do iVertex = 1,nVertices
+ do i=1,vertexDegree
+ iCell = cellsOnVertex(i,iVertex)
+ do k = 1,maxLevelCell(iCell)
+ pv_cell(k,iCell) = pv_cell(k,iCell) &
+ + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
+ / areaCell(iCell)
+ enddo
enddo
enddo
!
! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells and distance-1 ghost cells )
+ ! ( this computes pv_edge at all edges bounding real cells )
!
pv_edge(:,:) = 0.0
do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
- if(iEdge <= nEdges) then
- do k=1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
+ do i=1,vertexDegree
+ iEdge = edgesOnVertex(i,iVertex)
+ do k=1,maxLevelEdgeBot(iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
enddo
- endif
end do
end do
!
- ! Modify PV edge with upstream bias.
+ ! Compute gradient of PV in normal direction
+ ! ( this computes gradPVn for all edges bounding real cells )
!
+ gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
+ - pv_cell(k,cellsOnEdge(1,iEdge))) &
+ / dcEdge(iEdge)
enddo
enddo
-
!
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
+ ! Compute gradient of PV in the tangent direction
+ ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
- pv_cell(:,:) = 0.0
- do iVertex = 1, nVertices
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- if (iCell <= nCells) then
- do k = 1,nVertLevels
- pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
- enddo
- endif
- enddo
+ do iEdge = 1,nEdges
+ do k = 1,maxLevelEdgeBot(iEdge)
+ gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
+ - pv_vertex(k,verticesOnEdge(1,iEdge))) &
+ /dvEdge(iEdge)
+ enddo
enddo
!
- ! Compute gradient of PV in normal direction
- ! ( this computes gradPVn for all edges bounding real cells )
- !
- gradPVn(:,:) = 0.0
- do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
- do k = 1,nVertLevels
- gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
- dcEdge(iEdge)
- enddo
- endif
- enddo
-
-
! Modify PV edge with upstream bias.
!
do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
+ do k = 1,maxLevelEdgeBot(iEdge)
+ pv_edge(k,iEdge) = pv_edge(k,iEdge) &
+ - 0.5 * dt* ( u(k,iEdge) * gradPVn(k,iEdge) &
+ + v(k,iEdge) * gradPVt(k,iEdge) )
enddo
enddo
!
- ! Compute sea surface height
+ ! equation of state
!
- do iCell=1,nCells
- ssh(iCell) = h(1,iCell) - hZLevel(1)
- enddo
+ ! For an isopycnal model, density should remain constant.
+ if (config_vert_grid_type.eq.'zlevel') then
+ call equation_of_state(s,grid)
+ endif
!
- ! equation of state
+ ! Pressure
+ ! This section must be after computing rho
!
- ! For a isopycnal model, density should remain constant.
- if (config_vert_grid_type.eq.'zlevel') then
+ if (config_vert_grid_type.eq.'isopycnal') then
+
+ ! For Isopycnal model.
+ ! Compute pressure at top of each layer, and then
+ ! Montgomery Potential.
+ allocate(pTop(nVertLevels))
do iCell=1,nCells
- do k=1,nVertLevels
- ! Linear equation of state, for the time being
- rho(k,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*tracers(s % index_temperature,k,iCell) &
- + 7.6e-4*tracers(s % index_salinity,k,iCell))
- end do
- end do
- endif
+ ! assume atmospheric pressure at the surface is zero for now.
+ pTop(1) = 0.0
+ ! For isopycnal mode, p is the Montgomery Potential.
+ ! At top layer it is g*SSH, where SSH may be off by a
+ ! constant (ie, h_s can be relative to top or bottom)
+ MontPot(1,iCell) = gravity &
+ * (h_s(iCell) + sum(h(1:nVertLevels,iCell)))
- ! For Isopycnal model.
- ! Compute mid- and bottom-depth of each layer, from bottom up.
- ! Then compute mid- and bottom-pressure of each layer, and
- ! Montgomery Potential, from top down
- !
- do iCell=1,nCells
+ do k=2,nVertLevels
+ pTop(k) = pTop(k-1) + rho(k-1,iCell)*gravity* h(k-1,iCell)
- ! h_s is ocean topography: elevation above lowest point,
- ! and is positive. z coordinates are pos upward, and z=0
- ! is at lowest ocean point.
- zTop(nVertLevels+1,iCell) = h_s(iCell)
- do k=nVertLevels,1,-1
- zMid(k,iCell) = zTop(k+1,iCell) + 0.5*h(k,iCell)
- zTop(k,iCell) = zTop(k+1,iCell) + h(k,iCell)
- end do
+ ! from delta M = p delta / rho
+ MontPot(k,iCell) = MontPot(k-1,iCell) &
+ + pTop(k)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
+ end do
- ! assume atmospheric pressure at the surface is zero for now.
- pTop(1,iCell) = 0.0
- do k=1,nVertLevels
- delta_p = rho(k,iCell)*gravity* h(k,iCell)
- p(k ,iCell) = pTop(k,iCell) + 0.5*delta_p
- pTop(k+1,iCell) = pTop(k,iCell) + delta_p
- end do
+ end do
+ deallocate(pTop)
- MontPot(1,iCell) = gravity * zTop(1,iCell)
- do k=2,nVertLevels
- ! from delta M = p delta / rho
- MontPot(k,iCell) = MontPot(k-1,iCell) &
- + pTop(k,iCell)*(1.0/rho(k,iCell) - 1.0/rho(k-1,iCell))
- end do
+ elseif (config_vert_grid_type.eq.'zlevel') then
- end do
+ ! For z-level model.
+ ! Compute pressure at middle of each level.
+ ! At k=1, where p is pressure at a depth of hZLevel(1)/2, not
+ ! pressure at middle of layer including SSH.
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if(cell1<=nCells .and. cell2<=nCells) then
- do k=1,nVertLevels
- zTopEdge(k,iEdge) = (zTop(k,cell1)+zTop(k,cell2))/2.0
- zMidEdge(k,iEdge) = (zMid(k,cell1)+zMid(k,cell2))/2.0
- enddo
- zTopEdge(nVertLevels+1,iEdge) = ( zTop(nVertLevels+1,cell1) &
- + zTop(nVertLevels+1,cell2))/2.0
- endif
- enddo
+ do iCell=1,nCells
+ ! compute pressure for z-level coordinates
+ ! assume atmospheric pressure at the surface is zero for now.
+ pressure(1,iCell) = rho(1,iCell)*gravity &
+ * (h(1,iCell)-0.5*hZLevel(1))
- ! For z-level model.
- ! Compute pressure at middle of each level.
- ! pZLevel and p should only differ at k=1, where p is
- ! pressure at middle of layer including SSH, and pZLevel is
- ! pressure at a depth of hZLevel(1)/2.
- !
- do iCell=1,nCells
- ! compute pressure for z-level coordinates
- ! assume atmospheric pressure at the surface is zero for now.
- pZLevel(1,iCell) = rho(1,iCell)*gravity &
- * (h(1,iCell)-0.5*hZLevel(1))
- do k=2,nVertLevels
- pZLevel(k,iCell) = pZLevel(k-1,iCell) &
- + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
- + rho(k ,iCell)*hZLevel(k ))
- end do
+ do k=2,maxLevelCell(iCell)
+ pressure(k,iCell) = pressure(k-1,iCell) &
+ + 0.5*gravity*( rho(k-1,iCell)*hZLevel(k-1) &
+ + rho(k ,iCell)*hZLevel(k ))
+ end do
- end do
+ end do
- ! compute vertical velocity
+ endif
+
+ !
+ ! vertical velocity
+ !
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
wTop=0.0
@@ -1614,34 +1551,26 @@
! Compute div(u) for each cell
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
!
- allocate(div_u(nVertLevels,nCells))
+ allocate(div_u(nVertLevels,nCells+1))
div_u(:,:) = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge)
- div_u(k,cell1) = div_u(k,cell1) + flux
- end do
- end if
- if (cell2 <= nCells) then
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge)
- div_u(k,cell2) = div_u(k,cell2) - flux
- end do
- end if
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ flux = u(k,iEdge) * dvEdge(iEdge)
+ div_u(k,cell1) = div_u(k,cell1) + flux
+ div_u(k,cell2) = div_u(k,cell2) - flux
+ end do
end do
do iCell=1,nCells
- do k=1,nVertLevels
+ do k=1,maxLevelCell(iCell)
div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
end do
! Vertical velocity at bottom is zero.
- ! this next line can be set permanently somewhere upon startup.
- wTop(nVertLevels+1,iCell) = 0.0
- do k=nVertLevels,1,-1
+ wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+ do k=maxLevelCell(iCell),1,-1
wTop(k,iCell) = wTop(k+1,iCell) - div_u(k,iCell)*h(k,iCell)
end do
@@ -1650,7 +1579,6 @@
endif
-
end subroutine compute_solve_diagnostics
@@ -1696,4 +1624,261 @@
end subroutine enforce_boundaryEdge
+
+ subroutine equation_of_state(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ integer, dimension(:), pointer :: maxLevelCell
+ real (kind=RKIND), dimension(:,:), pointer :: rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer :: nCells, iCell, k
+ type (dm_info) :: dminfo
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+
+ maxLevelCell => grid % maxLevelCell % array
+ nCells = grid % nCells
+
+ if (config_eos_type.eq.'linear') then
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+ ! Linear equation of state
+ rho(k,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(s % index_temperature,k,iCell) &
+ + 7.6e-4*tracers(s % index_salinity,k,iCell))
+ end do
+ end do
+
+ elseif (config_eos_type.eq.'jm') then
+
+ call equation_of_state_jm(s, grid)
+
+ else
+ print *, ' Incorrect choice of config_eos_type:',&
+ config_eos_type
+ call dmpar_abort(dminfo)
+ endif
+
+ end subroutine equation_of_state
+
+
+ subroutine equation_of_state_jm(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! This module contains routines necessary for computing the density
+ ! from model temperature and salinity using an equation of state.
+ !
+ ! The UNESCO equation of state computed using the
+ ! potential-temperature-based bulk modulus from Jackett and
+ ! McDougall, JTECH, Vol.12, pp 381-389, April, 1995.
+ !
+ ! Input: grid - grid metadata
+ ! s - state: tracers
+ !
+ ! Output: s - state: computed density
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
+
+
+ type (dm_info) :: dminfo
+ integer :: iEdge, iCell, iVertex, k
+
+ integer :: nCells, nEdges, nVertices, nVertLevels
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ zMidZLevel, pRefEOS
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND) :: &
+ TQ,SQ, &! adjusted T,S
+ BULK_MOD, &! Bulk modulus
+ RHO_S, &! density at the surface
+ DRDT0, &! d(density)/d(temperature), for surface
+ DRDS0, &! d(density)/d(salinity ), for surface
+ DKDT, &! d(bulk modulus)/d(pot. temp.)
+ DKDS, &! d(bulk modulus)/d(salinity )
+ SQR,DENOMK, &! work arrays
+ WORK1, WORK2, WORK3, WORK4, T2, depth
+
+ real (kind=RKIND) :: &
+ tmin, tmax, &! valid temperature range for level k
+ smin, smax ! valid salinity range for level k
+
+ real (kind=RKIND) :: p, p2 ! temporary pressure scalars
+
+!-----------------------------------------------------------------------
+!
+! UNESCO EOS constants and JMcD bulk modulus constants
+!
+!-----------------------------------------------------------------------
+
+ !*** for density of fresh water (standard UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ unt0 = 999.842594, &
+ unt1 = 6.793952e-2, &
+ unt2 = -9.095290e-3, &
+ unt3 = 1.001685e-4, &
+ unt4 = -1.120083e-6, &
+ unt5 = 6.536332e-9
+
+ !*** for dependence of surface density on salinity (UNESCO)
+
+ real (kind=RKIND), parameter :: &
+ uns1t0 = 0.824493 , &
+ uns1t1 = -4.0899e-3, &
+ uns1t2 = 7.6438e-5, &
+ uns1t3 = -8.2467e-7, &
+ uns1t4 = 5.3875e-9, &
+ unsqt0 = -5.72466e-3, &
+ unsqt1 = 1.0227e-4, &
+ unsqt2 = -1.6546e-6, &
+ uns2t0 = 4.8314e-4
+
+ !*** from Table A1 of Jackett and McDougall
+
+ real (kind=RKIND), parameter :: &
+ bup0s0t0 = 1.965933e+4, &
+ bup0s0t1 = 1.444304e+2, &
+ bup0s0t2 = -1.706103 , &
+ bup0s0t3 = 9.648704e-3, &
+ bup0s0t4 = -4.190253e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0s1t0 = 5.284855e+1, &
+ bup0s1t1 = -3.101089e-1, &
+ bup0s1t2 = 6.283263e-3, &
+ bup0s1t3 = -5.084188e-5
+
+ real (kind=RKIND), parameter :: &
+ bup0sqt0 = 3.886640e-1, &
+ bup0sqt1 = 9.085835e-3, &
+ bup0sqt2 = -4.619924e-4
+
+ real (kind=RKIND), parameter :: &
+ bup1s0t0 = 3.186519 , &
+ bup1s0t1 = 2.212276e-2, &
+ bup1s0t2 = -2.984642e-4, &
+ bup1s0t3 = 1.956415e-6
+
+ real (kind=RKIND), parameter :: &
+ bup1s1t0 = 6.704388e-3, &
+ bup1s1t1 = -1.847318e-4, &
+ bup1s1t2 = 2.059331e-7, &
+ bup1sqt0 = 1.480266e-4
+
+ real (kind=RKIND), parameter :: &
+ bup2s0t0 = 2.102898e-4, &
+ bup2s0t1 = -1.202016e-5, &
+ bup2s0t2 = 1.394680e-7, &
+ bup2s1t0 = -2.040237e-6, &
+ bup2s1t1 = 6.128773e-8, &
+ bup2s1t2 = 6.207323e-10
+
+ rho => s % rho % array
+ tracers => s % tracers % array
+
+ nCells = grid % nCells
+ maxLevelCell => grid % maxLevelCell % array
+ nVertLevels = grid % nVertLevels
+ zMidZLevel => grid % zMidZLevel % array
+
+
+! Jackett and McDougall
+ tmin = -2.0 ! valid pot. temp. range
+ tmax = 40.0
+ smin = 0.0 ! valid salinity, in psu
+ smax = 42.0
+
+ ! This could be put in a startup routine.
+ ! Note I am using zMidZlevel, so pressure on top level does
+ ! not include SSH contribution. I am not sure if that matters.
+
+! This function computes pressure in bars from depth in meters
+! using a mean density derived from depth-dependent global
+! average temperatures and salinities from Levitus 1994, and
+! integrating using hydrostatic balance.
+
+ allocate(pRefEOS(nVertLevels))
+ do k = 1,nVertLevels
+ depth = -zMidZLevel(k)
+ pRefEOS(k) = 0.059808*(exp(-0.025*depth) - 1.0) &
+ + 0.100766*depth + 2.28405e-7*depth**2
+ enddo
+
+ do iCell=1,nCells
+ do k=1,maxLevelCell(iCell)
+
+ p = pRefEOS(k)
+ p2 = pRefEOS(k)*pRefEOS(k)
+
+ SQ = max(min(tracers(s%index_salinity,k,iCell),smax),smin)
+ TQ = max(min(tracers(s%index_temperature,k,iCell),tmax),tmin)
+
+ SQR = sqrt(SQ)
+ T2 = TQ*TQ
+
+ !***
+ !*** first calculate surface (p=0) values from UNESCO eqns.
+ !***
+
+ WORK1 = uns1t0 + uns1t1*TQ + &
+ (uns1t2 + uns1t3*TQ + uns1t4*T2)*T2
+ WORK2 = SQR*(unsqt0 + unsqt1*TQ + unsqt2*T2)
+
+ RHO_S = unt1*TQ + (unt2 + unt3*TQ + (unt4 + unt5*TQ)*T2)*T2 &
+ + (uns2t0*SQ + WORK1 + WORK2)*SQ
+
+ !***
+ !*** now calculate bulk modulus at pressure p from
+ !*** Jackett and McDougall formula
+ !***
+
+ WORK3 = bup0s1t0 + bup0s1t1*TQ + &
+ (bup0s1t2 + bup0s1t3*TQ)*T2 + &
+ p *(bup1s1t0 + bup1s1t1*TQ + bup1s1t2*T2) + &
+ p2*(bup2s1t0 + bup2s1t1*TQ + bup2s1t2*T2)
+ WORK4 = SQR*(bup0sqt0 + bup0sqt1*TQ + bup0sqt2*T2 + &
+ bup1sqt0*p)
+
+ BULK_MOD = bup0s0t0 + bup0s0t1*TQ + &
+ (bup0s0t2 + bup0s0t3*TQ + bup0s0t4*T2)*T2 + &
+ p *(bup1s0t0 + bup1s0t1*TQ + &
+ (bup1s0t2 + bup1s0t3*TQ)*T2) + &
+ p2*(bup2s0t0 + bup2s0t1*TQ + bup2s0t2*T2) + &
+ SQ*(WORK3 + WORK4)
+
+ DENOMK = 1.0/(BULK_MOD - p)
+
+ rho(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
+
+ end do
+ end do
+
+ deallocate(pRefEOS)
+
+ end subroutine equation_of_state_jm
+
end module time_integration
</font>
</pre>