<p><b>mpetersen@lanl.gov</b> 2010-04-15 16:21:11 -0600 (Thu, 15 Apr 2010)</p><p>Added vertical advection to height and tracer equations.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-04-09 17:15:43 UTC (rev 188)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-04-15 22:21:11 UTC (rev 189)
@@ -92,6 +92,8 @@
var real h ( nVertLevels nCells Time ) iro h - -
var real temperature ( nVertLevels nCells Time ) iro temperature tracers dynamics
var real salinity ( nVertLevels nCells Time ) iro salinity tracers dynamics
+var real tracer1 ( nVertLevels nCells Time ) iro tracer1 tracers testing
+var real tracer2 ( nVertLevels nCells Time ) iro tracer2 tracers testing
# Diagnostic fields: only written to output
var real v ( nVertLevels nEdges Time ) o v - -
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F 2010-04-09 17:15:43 UTC (rev 188)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F 2010-04-15 22:21:11 UTC (rev 189)
@@ -36,11 +36,12 @@
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, zbot, pmid, pbot, MontPot
+ pv_cell, gradPVn, gradPVt, zmid, zbot, pmid, pbot, MontPot, w
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
real (kind=RKIND) :: localCFL, localSum
integer :: elementIndex, variableIndex, nVariables, nSums, nMaxes, nMins
- integer :: timeLevel
+ integer :: timeLevel,k
integer, parameter :: kMaxVariables = 1024 ! this must be a little more than double the number of variables to be reduced
@@ -62,7 +63,9 @@
h => state % h % array
u => state % u % array
+ tracers => state % tracers % array
v => state % v % array
+ w => state % w % array
h_edge => state % h_edge % array
circulation => state % circulation % array
vorticity => state % vorticity % array
@@ -178,6 +181,12 @@
MontPot(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
+ ! w vertical velocity
+ variableIndex = variableIndex + 1
+ call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
+ w(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ verticalSumMaxes(variableIndex))
+
nVariables = variableIndex
nSums = nVariables
nMins = nVariables
@@ -304,6 +313,10 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
+ ! w vertical velocity
+ variableIndex = variableIndex + 1
+ averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
+
! write out the data to files
if (dminfo % my_proc_id == IO_NODE) then
fileID = getFreeUnit()
@@ -341,6 +354,20 @@
state % CFLNumberGlobal % scalar = CFLNumberGlobal
deallocate(areaEdge)
+
+ ! mrp 100415 temp output
+ print '(a,10es25.15)', 'k=1 min h max h',&
+ minval(h(1,1:nCellsSolve)),maxval(h(1,1:nCellsSolve))
+ if (nVertLevels.gt.1) &
+ print '(a,10es25.15)', 'k=2+ min h max h',&
+ minval(h(2:nVertLevels,1:nCellsSolve)),maxval(h(2:nVertLevels,1:nCellsSolve))
+
+ do k=1,4
+ print '(a,i5,10es25.15)', 'tracer min max',&
+ k,minval(tracers(k,:,1:nCellsSolve)),maxval(tracers(k,:,1:nCellsSolve))
+ enddo
+ ! mrp 100415 temp output end
+
end subroutine computeGlobalDiagnostics
integer function getFreeUnit()
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-04-09 17:15:43 UTC (rev 188)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-04-15 22:21:11 UTC (rev 189)
@@ -23,7 +23,9 @@
integer :: i, iCell, iEdge, iVtx, iLevel
type (block_type), pointer :: block_ptr
+ real (kind=RKIND), dimension(:), pointer :: xCell,yCell
real (kind=RKIND), dimension(:,:), pointer :: h, u, u_src, rho
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
real (kind=RKIND) :: delta_rho
integer :: nCells, nEdges, nVertices, nVertLevels
@@ -76,6 +78,8 @@
stop
end if
+!print *, temperature_INDEX
+
! mrp 100121:
!
! Initialize u_src, rho, alpha
@@ -87,8 +91,11 @@
h => block_ptr % time_levs(1) % state % h % array
u => block_ptr % time_levs(1) % state % u % array
rho => block_ptr % time_levs(1) % state % rho % array
+ tracers => block_ptr % 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
nCells = block_ptr % mesh % nCells
nEdges = block_ptr % mesh % nEdges
@@ -106,7 +113,16 @@
end do
endif
- ! define the density for multiple layers
+ tracers = 0.0
+ do iCell = 1,nCells
+ do iLevel = 1,nVertLevels
+ tracers(3,iLevel,iCell) = xCell(iCell)
+ if (xCell(iCell)<500*1e3) tracers(4,iLevel,iCell) = 1.0
+ tracers(1,iLevel,iCell) = xCell(iCell)
+ tracers(2,iLevel,iCell) = 1.0
+ enddo
+ enddo
+
delta_rho=0.0
do iLevel = 1,nVertLevels
rho(iLevel,1) = delta_rho*(iLevel-1)
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-04-09 17:15:43 UTC (rev 188)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-04-15 22:21:11 UTC (rev 189)
@@ -190,6 +190,19 @@
+ rk_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
block % time_levs(2) % state % h % array(:,:) = block % time_levs(2) % state % h % array(:,:) &
+ rk_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
+
+ ! mrp 100415 compute vertical velocity
+ do iCell=1,block % mesh % nCellsSolve
+ do k=1,block % mesh % nVertLevels-1
+ block % time_levs(2) % state % w % array(k,iCell) = &
+ block % intermediate_step(TEND) % w % array(k,iCell) &
+ / block % time_levs(2) % state % h % array(k+1,iCell)
+ end do
+ end do
+
+ block % time_levs(2) % state % w % array(:,:) = block % intermediate_step(TEND) % w % array(:,:)
+ ! mrp 100415 compute vertical velocity end
+
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
block % time_levs(2) % state % tracers % array(:,k,iCell) = &
@@ -278,6 +291,7 @@
h => s % h % array
u => s % u % array
v => s % v % array
+ Fv => tend % w % array
h_edge => s % h_edge % array
circulation => s % circulation % array
vorticity => s % vorticity % array
@@ -309,7 +323,6 @@
tend_h => tend % h % array
tend_u => tend % u % array
- Fv => tend % w % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -343,33 +356,39 @@
end if
end do
- ! mrp100409 computation of h tendency was, only had div.(hu) term
+ ! mrp 100409 computation of h tendency was, only had div.(hu) term
!do iCell=1,grid % nCellsSolve
! do k=1,nVertLevels
! tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
! end do
!end do
- ! mrp100409 replace with
+ ! mrp 100409 replace with
+
+ if (1==1) then ! isopycnal coordinates
+ Fv=0.0 ! vertical fluxes are zero
do iCell=1,grid % nCellsSolve
do k=1,nVertLevels
+ tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ end do
+ end do
+
+ else ! z-level with vertical advection
+ do iCell=1,grid % nCellsSolve
+ do k=1,nVertLevels
! tend_h is just the -div.(hu) term at this point:
tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
end do
- k=nVertLevels
- Fv(k,iCell) = tend_h(k,iCell)*h(k,icell)
+ ! change nvertlevels to kmt later
+ ! this next line can be set permanently somewhere upon startup.
+ Fv(nVertLevels,iCell) = 0.0
- ! note if you substitute line above this is just tend_h=0.
- ! so this line can be changed to tend_h=0 in the future.
- tend_h(k,iCell) = tend_h(k,iCell) &
- - (Fv(k,iCell))/h(k,icell)
+ do k=nVertLevels,2,-1
- do k=nVertLevels-1,-1,2
-
! F^v_{k-1/2} =
! F^v_{k+1/2} - </font>
<font color="red">abla_h \cdot \left( F^h \right) h_k
! note +tend_h because tend_h is -div.(hu)
- Fv(k,iCell) = Fv(k+1,iCell) + tend_h(k,iCell)*h(k,icell)
+ Fv(k-1,iCell) = Fv(k,iCell) + tend_h(k,iCell)*h(k,iCell)
! now tend_h is div.(hu) + d/dz(hw):
! - </font>
<font color="gray">abla_h \cdot F^h_k
@@ -377,17 +396,27 @@
! note if you substitute line above this is just tend_h=0.
! so this line can be changed to tend_h=0 in the future.
tend_h(k,iCell) = tend_h(k,iCell) &
- - (Fv(k,iCell) - Fv(k+1,iCell))/h(k,icell)
+ - (Fv(k-1,iCell) - Fv(k,iCell))/h(k,iCell)
end do
k=1
- Fv(k,iCell) = 0.0
+ ! Surface tracer Fv(k,iCell) = 0.0 is not used
tend_h(k,iCell) = tend_h(k,iCell) &
- - ( - Fv(k+1,iCell))/h(k,icell)
+ - ( - Fv(k,iCell))/h(k,iCell)
end do
- ! mrp100409 end
+ endif ! coordinate type
+ ! mrp 100409 end
+!print *, 'ncells,grid % nCellsSolve, size(Fv)',ncells,grid % nCellsSolve, size(Fv,1), size(Fv,2)
+! print '(a,i5,f10.5)', &
+! 'k, min(tend_h(k,2:)),max(tend_h(k,2:)),min(Fv(k,:)),max(Fv(k,:))'
+! do k=1,nVertLevels
+! print '(i5,10es10.2)', &
+! k,minval(tend_h(k,1:grid % ncellssolve)),maxval(tend_h(k,1:grid % ncellssolve)), &
+! minval(Fv(k,1:grid % ncellssolve)),maxval(Fv(k,1:grid % ncellssolve))
+! enddo
+
#ifdef LANL_FORMULATION
!
! Compute u (normal) velocity tendency for each edge (cell face)
@@ -485,7 +514,8 @@
integer :: iCell, iEdge, k, iTracer, cell1, cell2, &
nEdges, nCells, nVertLevels
- real (kind=RKIND) :: flux, tracer_edge
+ real (kind=RKIND) :: flux, tracer_edge, &
+ Tmid(num_tracers,grid % nVertLevels)
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
@@ -497,11 +527,12 @@
integer, dimension(:,:), pointer :: cellsOnEdge
u => s % u % array
- tracers => s % tracers % array
+ h => s % h % array
+ tracers => s % tracers % array
h_edge => s % h_edge % array
+ Fv => tend % w % array
tend_tr => tend % tracers % array
- Fv => tend % w % array
areaCell => grid % areaCell % array
cellsOnEdge => grid % cellsOnEdge % array
@@ -526,13 +557,59 @@
end if
end do
+ ! mrp 100409 computation of h tendency was, only had div.(huT) term
+ !do iCell=1,grid % nCellsSolve
+ ! do k=1,grid % nVertLevelsSolve
+ ! do iTracer=1,num_tracers
+ ! tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) / areaCell(iCell)
+ ! end do
+ ! end do
+ !end do
+ ! mrp 100409 replace with
do iCell=1,grid % nCellsSolve
+
+ ! Surface tracer Tmid(:,1) is not used
+ ! For now, Tmid is just the average of the upper and lower points.
+ ! Can use a fancier interpolation later.
+ ! change nVertLevels to KMT later
+ do k=1,nVertLevels-1
+ do iTracer=1,num_tracers
+ Tmid(iTracer,k) = ( tracers(iTracer,k ,iCell) &
+ + tracers(iTracer,k+1,iCell))/2.0
+ end do
+ end do
+
+ ! The next loop could be combined with the later two for better efficiency.
+ ! Keep separate now for clarity.
do k=1,grid % nVertLevelsSolve
do iTracer=1,num_tracers
+ ! this is just div.(huT) term:
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) / areaCell(iCell)
end do
end do
+
+ ! top layer: flux from below only
+ k=1
+ do iTracer=1,num_tracers
+ ! this is -div.(huT) - d/dz(hwT) terms.
+ ! note flux at surface is Fv(1,iCell)=0 so index can go to k=1.
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ - ( - Fv(k,iCell)*Tmid(iTracer,k))/h(k,icell)
+ end do
+
+ ! change nVertLevels to KMT later
+ do k=2,nVertLevels
+ do iTracer=1,num_tracers
+ ! this is -div.(huT) - d/dz(hwT) terms.
+ ! note flux at surface is Fv(1,iCell)=0 so index can go to k=1.
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ - ( Fv(k-1,iCell)*Tmid(iTracer,k-1) &
+ - Fv(k ,iCell)*Tmid(iTracer,k ))/h(k,icell)
+ end do
+ end do
+
end do
+ ! mrp 100409 end
end subroutine compute_scalar_tend
@@ -554,19 +631,18 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
- !mrp 100112:
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ zSurface
real (kind=RKIND), dimension(:,:), pointer :: &
+ vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, &
+ tend_h, tend_u, circulation, vorticity, ke, &
+ pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
zmid, zbot, pmid, pbot, MontPot, rho
- real (kind=RKIND), dimension(:), pointer :: zSurface
- real (kind=RKIND) :: delta_p
character :: c1*6
- !mrp 100112 end
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -576,6 +652,7 @@
h => s % h % array
u => s % u % array
v => s % v % array
+ w => s % w % array
vh => s % vh % array
h_edge => s % h_edge % array
tend_h => s % h % array
@@ -701,6 +778,18 @@
end do
!
+ ! Compute vertical velocity in each cell
+ ! w(k,iCell) is w_{k-1/2}, the velocity through the interface above cell k
+ ! Upon entry, w(k,iCell) is F^v_{k-1/2}, the vertical thickness flux,
+ ! computed in compute_tend.
+ !
+! print '(a,i5,f10.5)', &
+! 'k, min(w(k,2:)),max(w(k,2:))'
+! do k=1,nVertLevels
+! print '(i5,10es10.2)', &
+! k,minval(w(k,1:grid % ncellssolve)),maxval(w(k,1:grid % ncellssolve))
+! enddo
+ !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
</font>
</pre>