<p><b>mpetersen@lanl.gov</b> 2010-05-13 13:32:10 -0600 (Thu, 13 May 2010)</p><p>Remove NCAR formulation from zlevel ocean branch, which computes workpv*vh rather than q+u_diffusion in tend_u. Change variable name w to wBot, Fv to hFluxBot.<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-05-13 12:33:31 UTC (rev 268)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-05-13 19:32:10 UTC (rev 269)
@@ -131,7 +131,6 @@
var real pBot ( nVertLevels nCells Time ) o pBot - -
var real pMidZLevel ( nVertLevels nCells Time ) o pmidZLevel - -
var real MontPot ( nVertLevels nCells Time ) o MontPot - -
-var real w ( nVertLevels nCells Time ) o w - -
var real wBot ( nVertLevels nCells Time ) o wBot - -
# Other diagnostic variables: neither read nor written to any files
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-05-13 12:33:31 UTC (rev 268)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-13 19:32:10 UTC (rev 269)
@@ -8,10 +8,8 @@
use vector_reconstruction
! xsad 10-02-05 end
-
contains
-
subroutine timestep(domain, dt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step
@@ -174,14 +172,6 @@
) / block % intermediate_step(PROVIS) % h % array(k,iCell)
end do
- ! mrp 100415 compute vertical velocity
- do k=1,block % mesh % nVertLevels-1
- block % time_levs(PROVIS) % state % w % array(k,iCell) = &
- block % intermediate_step(TEND) % w % array(k,iCell) &
- / block % time_levs(PROVIS) % state % h % array(k+1,iCell)
- end do
- ! mrp 100415 compute vertical velocity end
-
end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
@@ -242,22 +232,6 @@
!call reconstruct(block % time_levs(2) % state, block % mesh)
! xsad 10-02-09 end
- ! mrp 100415 compute vertical velocity
- ! I am cheating here a little bit. Because I only compute the
- ! vertical flux within the compute_tend subroutine, it only shows
- ! up in the intermediate_step(TEND) % w variable. So this is
- ! w at the last stage of RK4. To fix it, I would need to compute
- ! the vertical flux within diagnostics to obtain w there using
- ! time_levs(2) for the horizontal 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(PROVIS) % state % h % array(k+1,iCell)
- end do
- end do
- ! mrp 100415 compute vertical velocity end
-
block => block % next
end do
@@ -290,8 +264,8 @@
zMidZLevel
real (kind=RKIND), dimension(grid % nVertLevels) :: w_dudzBotEdge
real (kind=RKIND), dimension(:,:), pointer :: &
- vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wBot, &
- tend_h, tend_u, Fv, circulation, vorticity, ke, pv_edge, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wBot, &
+ tend_h, tend_u, hFluxBot, circulation, vorticity, ke, pv_edge, &
divergence, MontPot, pmidZLevel, zMidEdge, zBotEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -310,14 +284,13 @@
u => s % u % array
v => s % v % array
wBot => s % wBot % array
- Fv => tend % w % array
+ hFluxBot => tend % wBot % array
h_edge => s % h_edge % array
circulation => s % circulation % array
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
pv_edge => s % pv_edge % array
- vh => s % vh % array
MontPot => s % MontPot % array
pmidZLevel => s % pmidZLevel % array
zBotEdge => s % zBotEdge % array
@@ -382,21 +355,21 @@
end do
if (config_vert_grid_type.eq.'isopycnal') then
- Fv=0.0 ! vertical fluxes are zero
+ hFluxBot=0.0 ! vertical fluxes are zero
elseif (config_vert_grid_type.eq.'zlevel') then
do iCell=1,nCells
! change nvertlevels to kmt later
! this next line can be set permanently somewhere upon startup.
- Fv(nVertLevels,iCell) = 0.0
+ hFluxBot(nVertLevels,iCell) = 0.0
do k=nVertLevels,2,-1
! 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-1,iCell) = Fv(k,iCell) + tend_h(k,iCell)*h(k,iCell)
+ hFluxBot(k-1,iCell) = hFluxBot(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
@@ -404,28 +377,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-1,iCell) - Fv(k,iCell))/h(k,iCell)
+ - (hFluxBot(k-1,iCell) - hFluxBot(k,iCell))/h(k,iCell)
end do
k=1
- ! Surface tracer Fv(k,iCell) = 0.0 is not used
+ ! Surface tracer hFluxBot(k,iCell) = 0.0 is not used
tend_h(k,iCell) = tend_h(k,iCell) &
- - ( - Fv(k,iCell))/h(k,iCell)
+ - ( - hFluxBot(k,iCell))/h(k,iCell)
end do
endif ! coordinate type
! mrp 100409 end
-!print *, 'ncells,grid % nCellsSolve, size(Fv)',ncells,grid % nCellsSolve, size(Fv,1), size(Fv,2)
+!print *, 'ncells,grid % nCellsSolve, size(hFluxBot)',ncells,grid % nCellsSolve, size(hFluxBot,1), size(hFluxBot,2)
! print '(a,i5,f10.5)', &
-! 'k, min(tend_h(k,2:)),max(tend_h(k,2:)),min(Fv(k,:)),max(Fv(k,:))'
+! 'k, min(tend_h(k,2:)),max(tend_h(k,2:)),min(hFluxBot(k,:)),max(hFluxBot(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))
+! minval(hFluxBot(k,1:grid % ncellssolve)),maxval(hFluxBot(k,1:grid % ncellssolve))
! enddo
-#ifdef LANL_FORMULATION
!
! Compute u (normal) velocity tendency for each edge (cell face)
!
@@ -536,8 +508,6 @@
end do
deallocate(fluxVert)
-#endif
-
! print '(a,i5,f10.5)', &
! 'k, min(tend_u(k,:)),max(tend_u(k,:))'
! do k=1,nVertLevels
@@ -545,32 +515,6 @@
! k,minval(tend_u(k,1:grid % nedgessolve)),maxval(tend_u(k,1:grid % nedgessolve))
! enddo
-#ifdef NCAR_FORMULATION
- !
- ! Compute u (normal) velocity tendency for each edge (cell face)
- !
- tend_u(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,nVertLevels
- vorticity_abs = fEdge(iEdge) + (circulation(k,vertex1) + circulation(k,vertex2)) / &
- (areaTriangle(vertex1) + areaTriangle(vertex2))
-
- workpv = 2.0 * vorticity_abs / (h(k,cell1) + h(k,cell2))
-
- tend_u(k,iEdge) = workpv * vh(k,iEdge) - &
- (ke(k,cell2) - ke(k,cell1) + &
- gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ) / &
- dcEdge(iEdge)
- end do
- end do
-#endif
-
end subroutine compute_tend
@@ -597,7 +541,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zsurface
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,Fv, h_edge, zmid, zbot
+ u,h,hFluxBot, h_edge, zmid, zbot
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -609,7 +553,7 @@
h => s % h % array
tracers => s % tracers % array
h_edge => s % h_edge % array
- Fv => tend % w % array
+ hFluxBot => tend % wBot % array
zmid => s % zmid % array
zbot => s % zbot % array
zsurface => s % zsurface % array
@@ -699,7 +643,7 @@
elseif (config_vert_tracer_adv.eq.'upwind') then
do k=1,nVertLevels-1
- if (Fv(k,iCell)>0.0) then
+ if (hFluxBot(k,iCell)>0.0) then
upwindCell = k+1
else
upwindCell = k
@@ -712,7 +656,7 @@
elseif (config_vert_tracer_adv.eq.'downwind') then
do k=1,nVertLevels-1
- if (Fv(k,iCell)>0.0) then
+ if (hFluxBot(k,iCell)>0.0) then
upwindCell = k
else
upwindCell = k+1
@@ -737,19 +681,19 @@
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.
+ ! note flux at surface is hFluxBot(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)
+ - ( - hFluxBot(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.
+ ! note flux at surface is hFluxBot(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)
+ - ( hFluxBot(k-1,iCell)*Tmid(iTracer,k-1) &
+ - hFluxBot(k ,iCell)*Tmid(iTracer,k ))/h(k,icell)
end do
end do
@@ -884,7 +828,7 @@
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
zSurface, hZLevel, zSurfaceEdge
real (kind=RKIND), dimension(:,:), pointer :: &
- vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wBot, &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wBot, &
circulation, vorticity, ke, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
zmid, zbot, zmidEdge, zbotEdge, pmid, pbot, MontPot, rho, temperature, salinity, pmidZLevel
@@ -900,9 +844,7 @@
h => s % h % array
u => s % u % array
v => s % v % array
- w => s % w % array
wBot => s % wBot % array
- vh => s % vh % array
h_edge => s % h_edge % array
circulation => s % circulation % array
vorticity => s % vorticity % array
@@ -1059,23 +1001,7 @@
end do
end do
-#ifdef NCAR_FORMULATION
!
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
- vh(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- do j=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- do k=1,nVertLevels
- vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * u(k,eoe) * h_edge(k,eoe)
- end do
- end do
- end do
-#endif
-
-
- !
! 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 )
!
@@ -1265,39 +1191,33 @@
elseif (config_vert_grid_type.eq.'zlevel') then
- !
- ! Compute div(u) for each cell
- ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
- !
- allocate(div_u(nVertLevels,nCells))
- div_u(:,:) = 0.0
- do iEdge=1,nEdges
+ !
+ ! Compute div(u) for each cell
+ ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+ !
+ allocate(div_u(nVertLevels,nCells))
+ div_u(:,:) = 0.0
+ do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells) then
do k=1,nVertLevels
-! should be this
flux = u(k,iEdge) * dvEdge(iEdge)
-! testing with this
-! flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
div_u(k,cell1) = div_u(k,cell1) + flux
end do
end if
if (cell2 <= nCells) then
do k=1,nVertLevels
-! should be this
flux = u(k,iEdge) * dvEdge(iEdge)
-! testing with this
-! flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
div_u(k,cell2) = div_u(k,cell2) - flux
end do
end if
- end do
+ end do
- do iCell=1,nCells
- do k=1,nVertLevels
- div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
- end do
+ do iCell=1,nCells
+ do k=1,nVertLevels
+ div_u(k,iCell) = div_u(k,iCell) / areaCell(iCell)
+ end do
! change nvertlevels to kmt later
! this next line can be set permanently somewhere upon startup.
@@ -1308,9 +1228,9 @@
end do
end do
+ deallocate(div_u)
endif
- deallocate(div_u)
end subroutine compute_solve_diagnostics
</font>
</pre>