<p><b>mpetersen@lanl.gov</b> 2010-04-22 15:22:15 -0600 (Thu, 22 Apr 2010)</p><p>Added vertical advection term to momentum equation.<br>
</p><hr noshade><pre><font color="gray">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-22 14:51:06 UTC (rev 201)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-04-22 21:22:15 UTC (rev 202)
@@ -78,8 +78,6 @@
stop
end if
-print *, 'index_temperature',index_temperature
-
! mrp 100121:
!
! Initialize u_src, rho, alpha
@@ -123,9 +121,12 @@
! tracers(2,iLevel,iCell) = 5.0 ! temperature
! tracers(3,iLevel,iCell) = 35.0 ! salinity
tracers(index_temperature,iLevel,iCell) = &
- 10.0* xCell(iCell)/2500e3 ! temperature
+ 10.0* xCell(iCell)/2500.e3
tracers(index_salinity,iLevel,iCell) = &
- 34.0 + 2* yCell(iCell)/4000e3! salinity
+ 34.0 + 2* yCell(iCell)/4000.e3
+ tracers(index_tracer1,iLevel,iCell) = 1.0
+ tracers(index_tracer2,iLevel,iCell) = &
+ (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
enddo
enddo
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-22 14:51:06 UTC (rev 201)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-04-22 21:22:15 UTC (rev 202)
@@ -114,9 +114,9 @@
rk_substep_weights(4) = 0.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do rk_step = 1, 4
! --- update halos for diagnostic variables
@@ -173,6 +173,15 @@
+ rk_substep_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell) &
) / 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(:,:)
@@ -182,6 +191,8 @@
end do
end if
+
+
!--- accumulate update (for RK4)
block => domain % blocklist
@@ -191,18 +202,6 @@
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) = &
@@ -214,11 +213,10 @@
end do
end do
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
!
! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
!
@@ -243,6 +241,22 @@
!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
@@ -266,11 +280,13 @@
type (grid_meta), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, upstream_bias
+ integer :: nCells, nEdges, nVertices, nVertLevels
- integer :: nCells, nEdges, nVertices, nVertLevels
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
+ upstream_bias, wHat, areaSumInv
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ real (kind=RKIND), dimension(grid % nVertLevels) :: w_dudz
real (kind=RKIND), dimension(:,:), pointer :: &
vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, &
tend_h, tend_u, Fv, circulation, vorticity, ke, pv_edge, &
@@ -364,7 +380,7 @@
!end do
! mrp 100409 replace with
- if (1==1) then ! isopycnal coordinates
+ if (1==2) then ! isopycnal coordinates
Fv=0.0 ! vertical fluxes are zero
do iCell=1,grid % nCellsSolve
do k=1,nVertLevels
@@ -427,7 +443,28 @@
cell2 = cellsOnEdge(2,iEdge)
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
-
+
+ ! mrp 100422 add -w*dudz term to tendancy
+ ! For isopycnal mode, Fv should be zero.
+ areaSumInv = 1.0/(areaCell(cell1)+areaCell(cell2))
+ ! change nvertlevels to kmt later
+ do k=1,nVertLevels-1
+ ! w =Fv/h. Also average w from cell center to edge
+ wHat = ( areaCell(cell1)*Fv(k,cell1)/h(k+1,cell1) &
+ + areaCell(cell2)*Fv(k,cell2)/h(k+1,cell2)) &
+ *areaSumInv
+ ! compute dudz at vertical interface with first order derivative.
+ w_dudz(k) = wHat * (u(k-1,iEdge)-u(k,iEdge)) &
+ * 2.0 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
+ end do
+ w_dudz(nVertLevels) = 0.0
+
+ tend_u(1,iEdge) = - 0.5 * w_dudz(1)
+ do k=2,nVertLevels
+ tend_u(k,iEdge) = - 0.5 * (w_dudz(k-1) + w_dudz(k))
+ enddo
+ ! mrp 100422 end: add -w*dudz term to tendancy
+
do k=1,nVertLevels
!
@@ -452,8 +489,8 @@
! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
! ) / dcEdge(iEdge)
! mrp 100112, changed to grad Montgomery potential:
- tend_u(k,iEdge) = &
- q &
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ + q &
+ u_diffusion &
- ( ke(k,cell2) - ke(k,cell1) &
+ MontPot(k,cell2) - MontPot(k,cell1) &
@@ -951,14 +988,14 @@
!
! equation of state
!
- do iCell=1,nCells
- do k=1,nVertLevels
- ! Linear equation of state, for the time being
- rho(k,iCell) = 1000.0*( 1.0 &
- - 2.5e-4*temperature(k,iCell) &
- + 7.6e-4*salinity(k,iCell))
- end do
- end do
+! do iCell=1,nCells
+! do k=1,nVertLevels
+! ! Linear equation of state, for the time being
+! rho(k,iCell) = 1000.0*( 1.0 &
+! - 2.5e-4*temperature(k,iCell) &
+! + 7.6e-4*salinity(k,iCell))
+! end do
+! end do
!mrp 100324 end
</font>
</pre>