<p><b>mpetersen@lanl.gov</b> 2011-01-31 09:52:01 -0700 (Mon, 31 Jan 2011)</p><p>Updates to vertical advection branch. I get reasonable results for global x3 grid to 100 days using config_vert_tracer_adv: centered, linear, flux3, flux4.<br>
</p><hr noshade><pre><font color="gray">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-28 19:02:48 UTC (rev 709)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_mpas_core.F 2011-01-31 16:52:01 UTC (rev 710)
@@ -76,7 +76,7 @@
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
- integer :: i
+ integer :: i, iEdge, iCell
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
@@ -84,6 +84,29 @@
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
call reconstruct(block % state % time_levs(1) % state, mesh)
+
+
+ ! initialize velocities and tracers on land to be -1e34
+ ! The reconstructed velocity on land will have values not exactly
+ ! -1e34 due to the interpolation of reconstruction.
+
+ do iEdge=1,block % mesh % nEdges
+ ! mrp 101115 note: in order to include flux boundary conditions, the following
+ ! line will need to change. Right now, set boundary edges between land and
+ ! water to have zero velocity.
+ block % state % time_levs(1) % state % u % array( &
+ block % mesh % maxLevelEdgeTop % array(iEdge)+1 &
+ :block % mesh % maxLevelEdgeBot % array(iEdge), iEdge) = 0.0
+
+ block % state % time_levs(1) % state % u % array( &
+ block % mesh % maxLevelEdgeBot % array(iEdge)+1: &
+ block % mesh % nVertLevels,iEdge) = -1e34
+ end do
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(1) % state % tracers % array( &
+ :, block % mesh % maxLevelCell % array(iCell)+1 &
+ :block % mesh % nVertLevels,iCell) = -1e34
+ end do
if (.not. config_do_restart) then
block % state % time_levs(1) % state % xtime % scalar = 0.0
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-28 19:02:48 UTC (rev 709)
+++ branches/ocean_projects/vert_adv_mrp/src/core_ocean/module_time_integration.F 2011-01-31 16:52:01 UTC (rev 710)
@@ -71,7 +71,7 @@
type (block_type), pointer :: block
type (state_type) :: provis
- integer :: rk_step
+ integer :: rk_step, iEdge
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
@@ -265,7 +265,7 @@
type (mesh_type), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j, kmax
+ vertex1, vertex2, eoe, i, j
integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
@@ -357,27 +357,39 @@
! for z-level, only compute height tendency for top layer.
if (config_vert_grid_type.eq.'isopycnal') then
- kmax = nVertLevels
+
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ 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
+ tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ end do
+ end do
+
elseif (config_vert_grid_type.eq.'zlevel') then
- kmax = 1
- endif
-
+
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- do k=1,kmax
+ do k=1,min(1,maxLevelEdgeTop(iEdge))
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,kmax
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
- end do
+ tend_h(1,iCell) = tend_h(1,iCell) / areaCell(iCell)
end do
+ endif ! config_vert_grid_type
+
!
! height tendency: vertical advection term -d/dz(hw)
!
@@ -599,10 +611,6 @@
do iEdge=1,grid % nEdgesSolve
- ! 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
@@ -612,18 +620,17 @@
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.
+ ! forcing in top layer only
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ + u_src(1,iEdge)/rho_ref/h_edge(1,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)
+ ! 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.
- ! old bottom drag, just linear friction
- ! du/dt = u/tau where tau=100 days.
- !tend_u(k,iEdge) = tend_u(k,iEdge) &
- ! - u(k,iEdge)/(100.0*86400.0)
+ 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)
endif
@@ -882,6 +889,8 @@
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
+ if (config_vert_grid_type.eq.'zlevel') then
+
allocate(tracerTop(num_tracers,nVertLevels+1))
tracerTop(:,1)=0.0
@@ -893,39 +902,63 @@
tracerTop(iTracer,k) = &
( tracers(iTracer,k-1,iCell) &
+tracers(iTracer,k ,iCell))/2.0
- end do
- end do
- end do
+ end do
+ end do
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ 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 ) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ end do
+ end do
+ end do
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
- else
- upwindCell = k-1
- endif
- do iTracer=1,num_tracers
- tracerTop(iTracer,k) = tracers(iTracer,upwindCell,iCell)
- end do
- end do
+ do k=2,maxLevelCell(iCell)
+ if (wTop(k,iCell)>0.0) then
+ upwindCell = k
+ else
+ upwindCell = k-1
+ endif
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k) = tracers(iTracer,upwindCell,iCell)
+ end do
+ end do
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ 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 ) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ end do
+ end do
end do
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
+ 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
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ 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 ) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ end do
+ end do
+ end do
elseif (config_vert_tracer_adv.eq.'flux3') then
@@ -933,7 +966,13 @@
! namelist, if desired.
flux3Coef = 1.0
do iCell=1,grid % nCellsSolve
- do k=2,maxLevelCell(iCell)
+ k=2
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ do k=3,maxLevelCell(iCell)-1
cSignWTop = sign(flux3Coef,wTop(k,iCell))
do iTracer=1,num_tracers
tracerTop(iTracer,k) = &
@@ -944,13 +983,33 @@
)/12.
end do
end do
- end do
+ k=maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ 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 ) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ end do
+ end do
+ end do
elseif (config_vert_tracer_adv.eq.'flux4') then
do iCell=1,grid % nCellsSolve
- do k=2,maxLevelCell(iCell)
+ k=2
do iTracer=1,num_tracers
+ tracerTop(iTracer,k) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ do k=3,maxLevelCell(iCell)-1
+ do iTracer=1,num_tracers
tracerTop(iTracer,k) = &
(- tracers(iTracer,k-2,iCell) &
+7.*tracers(iTracer,k-1,iCell) &
@@ -959,23 +1018,28 @@
)/12.
end do
end do
- end do
+ k=maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ tracerTop(iTracer,k) = &
+ hRatioZLevelK(k) *tracers(iTracer,k-1,iCell) &
+ + hRatioZLevelKm1(k)*tracers(iTracer,k ,iCell)
+ end do
+ tracerTop(:,maxLevelCell(iCell)+1)=0.0
+ 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 ) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
+ 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) &
- - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
- - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
- end do
- end do
- enddo ! iCell
-
deallocate(tracerTop)
+ endif ! ZLevel
+
!
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
!
@@ -1412,17 +1476,19 @@
! Compute kinetic energy in each cell
!
ke(:,:) = 0.0
- do iCell=1,nCells
- do i=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i,iCell)
- 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,maxLevelCell(iCell)
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeBot(iEdge)
+ ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ enddo
+ end do
+ do iCell = 1,nCells
+ do k = 1,maxLevelCell(iCell)
ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- end do
- end do
+ enddo
+ enddo
!
! Compute v (tangential) velocities
@@ -1431,7 +1497,9 @@
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- do k = 1,maxLevelEdgeBot(iEdge)
+ ! mrp 101115 note: in order to include flux boundary conditions,
+ ! the following loop may need to change to maxLevelEdgeBot
+ do k = 1,maxLevelEdgeTop(iEdge)
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
end do
</font>
</pre>