<p><b>mpetersen@lanl.gov</b> 2011-01-26 10:53:01 -0700 (Wed, 26 Jan 2011)</p><p>BRANCH COMMIT: Added high order vertical advection.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/vert_adv_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_ocean/Registry 2011-01-26 00:17:05 UTC (rev 705)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/Registry 2011-01-26 17:53:01 UTC (rev 706)
@@ -130,6 +130,9 @@
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
var persistent real zMidZLevel ( nVertLevels ) 0 - zMidZLevel mesh - -
var persistent real zTopZLevel ( nVertLevelsP1 ) 0 - zTopZLevel mesh - -
+var persistent real hMeanTopZLevel ( nVertLevels ) 0 - hMeanTopZLevel mesh - -
+var persistent real hRatioZLevelK ( nVertLevels ) 0 - hRatioZLevelK mesh - -
+var persistent real hRatioZLevelKm1 ( nVertLevels ) 0 - hRatioZLevelKm1 mesh - -
# Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
Modified: branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_mpas_core.F 2011-01-26 00:17:05 UTC (rev 705)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_mpas_core.F 2011-01-26 17:53:01 UTC (rev 706)
@@ -224,7 +224,7 @@
block_ptr => domain % blocklist
if (associated(block_ptr % next)) then
write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
- 'that there is only one block per processor.'
+ 'that there is only one block per processor.'
end if
call timer_start("global diagnostics")
@@ -253,7 +253,8 @@
integer :: iTracer, cell
real (kind=RKIND), dimension(:), pointer :: &
- hZLevel, zMidZLevel, zTopZLevel
+ hZLevel, zMidZLevel, zTopZLevel, &
+ hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -266,6 +267,9 @@
zMidZLevel => block % mesh % zMidZLevel % array
zTopZLevel => block % mesh % zTopZLevel % array
nVertLevels = block % mesh % nVertLevels
+ hMeanTopZLevel => block % mesh % hMeanTopZLevel % array
+ hRatioZLevelK => block % mesh % hRatioZLevelK % array
+ hRatioZLevelKm1 => block % mesh % hRatioZLevelKm1 % array
! These should eventually be in an input file. For now
! I just read them in from h(:,1).
@@ -281,6 +285,15 @@
zTopZLevel(k+1) = zTopZLevel(k)- hZLevel(k)
end do
+ hMeanTopZLevel(1) = 0.0
+ hRatioZLevelK(1) = 0.0
+ hRatioZLevelKm1(1) = 0.0
+ do k = 2,nVertLevels
+ hMeanTopZLevel(k) = 0.5*(hZLevel(k-1) + hZLevel(k))
+ hRatioZLevelK(k) = 0.5*hZLevel(k)/hMeanTopZLevel(k)
+ hRatioZLevelKm1(k) = 0.5*hZLevel(k-1)/hMeanTopZLevel(k)
+ end do
+
block => block % next
end do
Modified: branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-01-26 00:17:05 UTC (rev 705)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-01-26 17:53:01 UTC (rev 706)
@@ -707,7 +707,8 @@
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
- real (kind=RKIND), dimension(:), pointer :: zTopZLevel
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel, &
+ hRatioZLevelK, hRatioZLevelKm1
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop, boundaryMask
real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div, delsq_tracer
@@ -715,9 +716,8 @@
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
+ real (kind=RKIND) :: coef_3rd_order, flux3Coef, cSignWTop
-
u => s % u % array
h => s % h % array
boundaryCell=> grid % boundaryCell % array
@@ -732,6 +732,8 @@
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
zTopZLevel => grid % zTopZLevel % array
+ hRatioZLevelK => grid % hRatioZLevelK % array
+ hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
boundaryEdge => grid % boundaryEdge % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
@@ -882,17 +884,22 @@
!
allocate(tracerTop(num_tracers,nVertLevels+1))
tracerTop(:,1)=0.0
- do iCell=1,grid % nCellsSolve
- if (config_vert_tracer_adv.eq.'centered') then
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tracerTop(iTracer,k) = ( tracers(iTracer,k-1,iCell) &
- +tracers(iTracer,k ,iCell))/2.0
+ if (config_vert_tracer_adv.eq.'centered') then
+
+ do iCell=1,grid % nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k) = &
+ ( tracers(iTracer,k-1,iCell) &
+ +tracers(iTracer,k ,iCell))/2.0
+ end do
end do
- end do
+ end do
- elseif (config_vert_tracer_adv.eq.'upwind') then
+ elseif (config_vert_tracer_adv.eq.'upwind') then
+
+ do iCell=1,grid % nCellsSolve
do k=2,maxLevelCell(iCell)
if (wTop(k,iCell)>0.0) then
upwindCell = k
@@ -903,10 +910,61 @@
tracerTop(iTracer,k) = tracers(iTracer,upwindCell,iCell)
end do
end do
+ end do
- endif
+ elseif (config_vert_tracer_adv.eq.'linear') then
+
+ do iCell=1,grid % nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ ! Linear interpolation of tracer value at Top interface.
+ ! Note hRatio on the k side is multiplied by tracer at k-1
+ ! and hRatio on the Km1 (k-1) side is mult. by tracer at k.
+ tracerTop(iTracer,k) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ end do
+ end do
+
+ elseif (config_vert_tracer_adv.eq.'flux3') then
+
+ ! Hardwire flux3Coeff at 1.0 for now. Could add this to the
+ ! namelist, if desired.
+ flux3Coef = 1.0
+ do iCell=1,grid % nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ cSignWTop = sign(flux3Coef,wTop(k,iCell))
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k) = &
+ ( (-1.+ cSignWTop)*tracers(iTracer,k-2,iCell) &
+ +( 7.-3.*cSignWTop)*tracers(iTracer,k-1,iCell) &
+ +( 7.+3.*cSignWTop)*tracers(iTracer,k ,iCell) &
+ +(-1.- cSignWTop)*tracers(iTracer,k+1,iCell) &
+ )/12.
+ end do
+ end do
+ end do
+
+ elseif (config_vert_tracer_adv.eq.'flux4') then
+
+ do iCell=1,grid % nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k) = &
+ (- tracers(iTracer,k-2,iCell) &
+ +7.*tracers(iTracer,k-1,iCell) &
+ +7.*tracers(iTracer,k ,iCell) &
+ - tracers(iTracer,k+1,iCell) &
+ )/12.
+ end do
+ end do
+ end do
+
+ endif
tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ do iCell=1,grid % nCellsSolve
do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
@@ -914,8 +972,8 @@
- wTop(k+1,iCell)*tracerTop(iTracer,k+1))
end do
end do
+ enddo ! iCell
- enddo ! iCell
deallocate(tracerTop)
!
</font>
</pre>