<p><b>mpetersen@lanl.gov</b> 2011-12-06 15:52:00 -0700 (Tue, 06 Dec 2011)</p><p>Remove zlevel variables whereever possible, replace with full h variables. In particular, remove the hZLevel ratios in the vadv schemes. I still need to remove zlevel references from spline3.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_mpas_core.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_mpas_core.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -472,7 +472,7 @@
type (block_type), pointer :: block
integer :: iTracer, cell, cell1, cell2
- real (kind=RKIND) :: uhSum, hSum, sshEdge
+ real (kind=RKIND) :: uhSum, hSum, hEdge1
real (kind=RKIND), dimension(:), pointer :: &
hZLevel, zMidZLevel, zTopZLevel, &
hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
@@ -502,21 +502,15 @@
! hZLevel should be in the grid.nc and restart.nc file,
! and h for k=1 must be specified there as well.
+!!!!!!!! mrp 111206
+!!!! I still need to change spline3 from hZLevel to h variables
+!!!! then I can remove these variables.
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
- 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
-
! mrp 110601 For now, h is the variable saved in the restart file
! I am computing SSH here. In the future, could make smaller
! restart files for z-Level runs by saving SSH only.
@@ -550,21 +544,20 @@
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- sshEdge = 0.5*( &
- block % state % time_levs(1) % state % ssh % array(cell1) &
- + block % state % time_levs(1) % state % ssh % array(cell2) )
-
! uBtr = sum(u)/sum(h) on each column
- uhSum = (sshEdge + block % mesh % hZLevel % array(1)) &
- * block % state % time_levs(1) % state % u % array(1,iEdge)
- hSum = sshEdge + block % mesh % hZLevel % array(1)
+ uhSum = 0.0
+ hSum = 0.0
- do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ ! ocn_diagnostic_solve has not yet been called, so compute hEdge
+ ! just for this edge.
+ hEdge1 = 0.5*( &
+ block % state % time_levs(1) % state % h % array(k,cell1) &
+ + block % state % time_levs(1) % state % h % array(k,cell2) )
+
uhSum = uhSum &
- + block % mesh % hZLevel % array(k) &
- *block % state % time_levs(1) % state % u % array(k,iEdge)
- hSum = hSum &
- + block % mesh % hZLevel % array(k)
+ + hEdge1*block % state % time_levs(1) % state % u % array(k,iEdge)
+ hSum = hSum + hEdge1
enddo
block % state % time_levs(1) % state % uBtr % array(iEdge) = uhSum/hsum
@@ -592,32 +585,6 @@
endif
-!print *, '11 u ',minval(domain % blocklist % state % time_levs(1) % state % u % array(:,1:domain % blocklist % mesh % nEdgesSolve)), &
-! maxval(domain % blocklist % state % time_levs(1) % state % u % array(:,1:domain % blocklist % mesh % nEdgesSolve))
-!print *, '11 uBtr ',minval(domain % blocklist % state % time_levs(1) % state % uBtr % array(1:domain % blocklist % mesh % nEdgesSolve)), &
-! maxval(domain % blocklist % state % time_levs(1) % state % uBtr % array(1:domain % blocklist % mesh % nEdgesSolve))
-!print *, '11 uBcl ',minval(domain % blocklist % state % time_levs(1) % state % uBcl % array(:,1:domain % blocklist % mesh % nEdgesSolve)), &
-! maxval(domain % blocklist % state % time_levs(1) % state % uBcl % array(:,1:domain % blocklist % mesh % nEdgesSolve))
-
-
-! mrp temp testing - is uBcl vert sum zero?
-! do iEdge=1,block % mesh % nEdges
-! uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * block % state % time_levs(1) % state % uBcl % array(1,iEdge)
-! hSum = sshEdge + block % mesh % hZLevel % array(1)
-
-! do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
-! uhSum = uhSum + block % mesh % hZLevel % array(k) * block % state % time_levs(1) % state % uBcl % array(k,iEdge)
-! hSum = hSum + block % mesh % hZLevel % array(k)
-! enddo
-! block % state % time_levs(1) % state % FBtr % array(iEdge) = uhSum/hSum
-
-! enddo ! iEdge
-
-!print *, 'uBcl vert sum IC',minval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve)), &
-! maxval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve))
-
-! mrp temp testing - is uBcl vert sum zero? end
-
block => block % next
end do
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tendency.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tendency.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -453,7 +453,7 @@
maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
real (kind=RKIND), dimension(:), pointer :: zTopZLevel,zMidZLevel, &
- hRatioZLevelK, hRatioZLevelKm1, meshScalingDel2, meshScalingDel4
+ meshScalingDel2, meshScalingDel4
real (kind=RKIND), dimension(:), allocatable:: tracer2ndDer, tracersIn, tracersOut, posZMidZLevel, &
posZTopZLevel
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, boundaryMask
@@ -486,8 +486,6 @@
dcEdge => grid % dcEdge % array
zTopZLevel => grid % zTopZLevel % array
zMidZLevel => grid % zMidZLevel % array
- hRatioZLevelK => grid % hRatioZLevelK % array
- hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
boundaryEdge => grid % boundaryEdge % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
@@ -531,7 +529,7 @@
call mpas_timer_start("ocn_tend_scalar-vert adv")
- call ocn_tracer_vadv_tend(grid, wtop, tracers, tend_tr, err)
+ call ocn_tracer_vadv_tend(grid, h, wtop, tracers, tend_tr, err)
call mpas_timer_stop("ocn_tend_scalar-vert adv")
@@ -618,7 +616,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel, zTopZLevel
+ zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -682,7 +680,6 @@
h_s => grid % h_s % array
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
- hZLevel => grid % hZLevel % array
zTopZLevel => grid % zTopZLevel % array
deriv_two => grid % deriv_two % array
maxLevelCell => grid % maxLevelCell % array
@@ -707,9 +704,6 @@
! mrp 101115 note: in order to include flux boundary conditions, we will need to
! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
- ! mrp 110516 efficiency note: For z-level, only do this on level 1. h_edge for all
- ! lower levels is defined by hZlevel.
-
call mpas_timer_start("ocn_diagnostic_solve-hEdge")
coef_3rd_order = 0.
@@ -1086,27 +1080,11 @@
end do
deallocate(pTop)
- elseif (config_vert_grid_type.eq.'zlevel') then
+ else
- ! 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 iCell=1,nCells
- ! compute pressure for z-level coordinates
+ ! pressure for generalized coordinates
! assume atmospheric pressure at the surface is zero for now.
-
-! mrp 111206 pressure for z-level, old
-! pressure(1,iCell) = rho(1,iCell)*gravity &
-! * (h(1,iCell)-0.5*hZLevel(1))
-! 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
-
- ! pressure for generalized coordinates
pressure(1,iCell) = rho(1,iCell)*gravity &
* 0.5*h(1,iCell)
@@ -1117,6 +1095,7 @@
end do
! Compute zMid, the z-coordinate of the middle of the layer.
+ ! This is used for the rho g grad z momentum term.
k = maxLevelCell(iCell)
zMid(k:nVertLevels,iCell) = zTopZLevel(k+1) + 0.5*h(k,iCell)
@@ -1173,8 +1152,7 @@
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: u,wTop, h_edge
real (kind=RKIND), dimension(:,:), allocatable:: div_u
@@ -1193,7 +1171,6 @@
areaCell => grid % areaCell % array
cellsOnEdge => grid % cellsOnEdge % array
- hZLevel => grid % hZLevel % array
maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
dvEdge => grid % dvEdge % array
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_time_integration_rk4.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_time_integration_rk4.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -469,15 +469,12 @@
do iEdge=1,grid % nEdges
- ! I am using hZLevel here. This assumes that SSH is zero everywhere already,
- ! which should be the case if the barotropic mode is filtered.
- ! The more general case is to use sshEdge or h_edge.
- uhSum = (grid % hZLevel % array(1)) * tend_u(1,iEdge)
- hSum = grid % hZLevel % array(1)
+ uhSum = (h_edge(1,iEdge)) * tend_u(1,iEdge)
+ hSum = h_edge(1,iEdge)
do k=2,grid % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + grid % hZLevel % array(k) * tend_u(k,iEdge)
- hSum = hSum + grid % hZLevel % array(k)
+ uhSum = uhSum + h_edge(k,iEdge) * tend_u(k,iEdge)
+ hSum = hSum + h_edge(k,iEdge)
enddo
vertSum = uhSum/hSum
@@ -589,15 +586,12 @@
do iEdge=1,grid % nEdges
- ! I am using hZLevel here. This assumes that SSH is zero everywhere already,
- ! which should be the case if the barotropic mode is filtered.
- ! The more general case is to use sshedge or h_edge.
- uhSum = (grid % hZLevel % array(1)) * u(1,iEdge)
- hSum = grid % hZLevel % array(1)
+ uhSum = (h_edge(1,iEdge)) * u(1,iEdge)
+ hSum = h_edge(1,iEdge)
do k=2,grid % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + grid % hZLevel % array(k) * u(k,iEdge)
- hSum = hSum + grid % hZLevel % array(k)
+ uhSum = uhSum + h_edge(k,iEdge) * u(k,iEdge)
+ hSum = hSum + h_edge(k,iEdge)
enddo
vertSum = uhSum/hSum
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -70,7 +70,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_vadv_tend(grid, wTop, tracers, tend, err)!{{{
+ subroutine ocn_tracer_vadv_tend(grid, h, wTop, tracers, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -79,6 +79,7 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h, & !< Input: layer thickness
wTop !< Input: vertical tracer in top layer
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
@@ -124,8 +125,8 @@
if(.not.vadvOn) return
- call ocn_tracer_vadv_stencil_tend(grid, wTop, tracers, tend, err1)
- call ocn_tracer_vadv_spline_tend(grid, wTop, tracers, tend, err2)
+ call ocn_tracer_vadv_stencil_tend(grid, h, wTop, tracers, tend, err1)
+ call ocn_tracer_vadv_spline_tend(grid, h, wTop, tracers, tend, err2)
err = ior(err1, err2)
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -70,7 +70,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_vadv_spline_tend(grid, wTop, tracers, tend, err)!{{{
+ subroutine ocn_tracer_vadv_spline_tend(grid, h, wTop, tracers, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -79,6 +79,7 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h, & !< Input: layer thickness
wTop !< Input: vertical tracer in top layer
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
@@ -124,8 +125,8 @@
if(.not.splineOn) return
- call ocn_tracer_vadv_spline2_tend(grid, wTop, tracers, tend, err1)
- call ocn_tracer_vadv_spline3_tend(grid, wTop, tracers, tend, err2)
+ call ocn_tracer_vadv_spline2_tend(grid, h, wTop, tracers, tend, err1)
+ call ocn_tracer_vadv_spline3_tend(grid, h, wTop, tracers, tend, err2)
err = ior(err1, err2)
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline2.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline2.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline2.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -64,7 +64,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_vadv_spline2_tend(grid, wTop, tracers, tend, err)!{{{
+ subroutine ocn_tracer_vadv_spline2_tend(grid, h, wTop, tracers, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -73,6 +73,7 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h, & !< Input: layer thickness
wTop !< Input: vertical tracer in top layer
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
@@ -108,8 +109,6 @@
integer, dimension(:), pointer :: maxLevelCell
- real (kind=RKIND), dimension(:), pointer :: hRatioZLevelK, hRatioZLevelKm1
-
real (kind=RKIND), dimension(:,:,:), allocatable :: tracerTop
!-----------------------------------------------------------------
@@ -133,20 +132,18 @@
num_tracers = size(tracers, 1)
maxLevelCell => grid % maxLevelCell % array
- hRatioZLevelK => grid % hRatioZLevelK % array
- hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
-
allocate(tracerTop(num_tracers,nVertLevels+1,nCells))
do iCell=1,nCellsSolve
tracerTop(:,1,iCell) = tracers(:,1,iCell)
do k=2,maxLevelCell(iCell)
do iTracer=1,num_tracers
- ! 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.
+ ! Note h on the k side is multiplied by tracer at k-1
+ ! and h on the Km1 (k-1) side is mult. by tracer at k.
tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ ( h(k ,iCell)*tracers(iTracer,k-1,iCell) &
+ + h(k-1,iCell)*tracers(iTracer,k ,iCell) ) &
+ / (h(k-1,iCell) + h(k,iCell))
end do
end do
tracerTop(:,maxLevelCell(iCell)+1,iCell) = tracers(:,maxLevelCell(iCell),iCell)
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -65,7 +65,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_vadv_spline3_tend(grid, wTop, tracers, tend, err)!{{{
+ subroutine ocn_tracer_vadv_spline3_tend(grid, h, wTop, tracers, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -74,6 +74,7 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h, & !< Input: layer thickness
wTop !< Input: vertical tracer in top layer
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
@@ -109,8 +110,7 @@
integer, dimension(:), pointer :: maxLevelCell
- real (kind=RKIND), dimension(:), pointer :: hRatioZLevelK, &
- hRatioZLevelKm1, zTopZLevel, zMidZLevel
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel, zMidZLevel
real (kind=RKIND), dimension(:), allocatable :: tracer2ndDer, &
tracersIn, tracersOut, posZMidZLevel, posZTopZLevel
@@ -137,8 +137,8 @@
num_tracers = size(tracers, 1)
maxLevelCell => grid % maxLevelCell % array
- hRatioZLevelK => grid % hRatioZLevelK % array
- hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
+!!!!!!!! mrp 111206
+!!!! I still need to change spline3 from hZLevel to h variables
zMidZLevel => grid % zMidZLevel % array
zTopZLevel => grid % zTopZLevel % array
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -71,7 +71,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_vadv_stencil_tend(grid, wTop, tracers, tend, err)!{{{
+ subroutine ocn_tracer_vadv_stencil_tend(grid, h, wTop, tracers, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -80,6 +80,7 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h, & !< Input: layer thickness
wTop !< Input: vertical tracer in top layer
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
@@ -126,8 +127,8 @@
if(.not. stencilOn) return
call ocn_tracer_vadv_stencil2_tend(grid, wTop, tracers, tend, err1)
- call ocn_tracer_vadv_stencil3_tend(grid, wTop, tracers, tend, err2)
- call ocn_tracer_vadv_stencil4_tend(grid, wTop, tracers, tend, err3)
+ call ocn_tracer_vadv_stencil3_tend(grid, h, wTop, tracers, tend, err2)
+ call ocn_tracer_vadv_stencil4_tend(grid, h, wTop, tracers, tend, err3)
err = ior(err1, ior(err2, err3))
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil3.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil3.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil3.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -64,7 +64,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_vadv_stencil3_tend(grid, wTop, tracers, tend, err)!{{{
+ subroutine ocn_tracer_vadv_stencil3_tend(grid, h, wTop, tracers, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -73,6 +73,7 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h, & !< Input: layer thickness
wTop !< Input: vertical tracer in top layer
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
@@ -110,7 +111,6 @@
integer, dimension(:), pointer :: maxLevelCell
real (kind=RKIND) :: cSignWTop, flux3Coef
- real (kind=RKIND), dimension(:), pointer :: hRatioZLevelK, hRatioZLevelKm1
real (kind=RKIND), dimension(:,:,:), allocatable :: tracerTop
@@ -131,8 +131,6 @@
num_tracers = size(tracers, 1)
nVertLevels = grid % nVertLevels
maxLevelCell => grid % maxLevelCell % array
- hRatioZLevelK => grid % hRatioZLevelK % array
- hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
call mpas_timer_start("compute_scalar_tend-vert adv stencil 3")
@@ -148,9 +146,10 @@
tracerTop(:,1,iCell) = tracers(:,1,iCell)
k=2
do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ tracerTop(iTracer,k,iCell) = &
+ ( h(k,iCell)*tracers(iTracer,k-1,iCell) &
+ + h(k-1,iCell)*tracers(iTracer,k ,iCell) ) &
+ / (h(k-1,iCell) + h(k,iCell))
end do
do k=3,maxLevelCell(iCell)-1
cSignWTop = sign(flux3Coef,wTop(k,iCell))
@@ -165,9 +164,10 @@
end do
k=maxLevelCell(iCell)
do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ tracerTop(iTracer,k,iCell) = &
+ ( h(k,iCell)*tracers(iTracer,k-1,iCell) &
+ + h(k-1,iCell)*tracers(iTracer,k ,iCell) ) &
+ / (h(k-1,iCell) + h(k,iCell))
end do
tracerTop(:,maxLevelCell(iCell)+1,iCell) = tracers(:,maxLevelCell(iCell),iCell)
end do
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil4.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil4.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_tracer_vadv_stencil4.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -64,7 +64,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_vadv_stencil4_tend(grid, wTop, tracers, tend, err)!{{{
+ subroutine ocn_tracer_vadv_stencil4_tend(grid, h, wTop, tracers, tend, err)!{{{
!-----------------------------------------------------------------
!
@@ -73,6 +73,7 @@
!-----------------------------------------------------------------
real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h, & !< Input: layer thickness
wTop !< Input: vertical tracer in top layer
real (kind=RKIND), dimension(:,:,:), intent(in) :: &
@@ -110,7 +111,6 @@
integer, dimension(:), pointer :: maxLevelCell
real (kind=RKIND) :: cSingWTop, flux3Coef
- real (kind=RKIND), dimension(:), pointer :: hRatioZLevelK, hRatioZLevelKm1
real (kind=RKIND), dimension(:,:,:), allocatable :: tracerTop
@@ -133,8 +133,6 @@
num_tracers = size(tracers, 1)
nVertLevels = grid % nVertLevels
maxLevelCell => grid % maxLevelCell % array
- hRatioZLevelK => grid % hRatioZLevelK % array
- hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
allocate(tracerTop(num_tracers,nVertLevels+1,nCells))
@@ -144,9 +142,10 @@
tracerTop(:,1,iCell) = tracers(:,1,iCell)
k=2
do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ tracerTop(iTracer,k,iCell) = &
+ ( h(k ,iCell)*tracers(iTracer,k-1,iCell) &
+ + h(k-1,iCell)*tracers(iTracer,k ,iCell) ) &
+ / (h(k-1,iCell) + h(k,iCell))
end do
do k=3,maxLevelCell(iCell)-1
do iTracer=1,num_tracers
@@ -160,9 +159,10 @@
end do
k=maxLevelCell(iCell)
do iTracer=1,num_tracers
- tracerTop(iTracer,k,iCell) = &
- hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
- + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ tracerTop(iTracer,k,iCell) = &
+ ( h(k ,iCell)*tracers(iTracer,k-1,iCell) &
+ + h(k-1,iCell)*tracers(iTracer,k ,iCell) ) &
+ / (h(k-1,iCell) + h(k,iCell))
end do
tracerTop(:,maxLevelCell(iCell)+1,iCell) = tracers(:,maxLevelCell(iCell),iCell)
end do
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_pressure_grad.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -127,16 +127,6 @@
end do
enddo
else
-! mrp 111206 pressure for z-level, old, delete before trunk merge.
-! do iEdge=1,nEdgesSolve
-! cell1 = cellsOnEdge(1,iEdge)
-! cell2 = cellsOnEdge(2,iEdge)
-! do k=1,maxLevelEdgeTop(iEdge)
-! tend(k,iEdge) = tend(k,iEdge) &
-! - rho0Inv*( pressure(k,cell2) &
-! - pressure(k,cell1) )/dcEdge(iEdge)
-! end do
-! enddo
! pressure for generalized coordinates
! -1/rho_0 (grad p_k + rho g grad z_k^{mid})
Modified: branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vmix.F 2011-12-06 21:35:39 UTC (rev 1240)
+++ branches/ocean_projects/ale_vert_coord/src/core_ocean/mpas_ocn_vmix.F 2011-12-06 22:52:00 UTC (rev 1241)
@@ -319,7 +319,6 @@
! mrp 110315 efficiency note: for z-level, could precompute
! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
! h_edge is computed in compute_solve_diag, and is not available yet.
- ! This could be removed if hZLevel used instead.
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
</font>
</pre>