<p><b>mpetersen@lanl.gov</b> 2010-05-03 14:35:52 -0600 (Mon, 03 May 2010)</p><p>Added horizontal and vertical tracer diffusion to zlevel branch.<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-05-03 20:22:06 UTC (rev 235)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-05-03 20:35:52 UTC (rev 236)
@@ -129,17 +129,24 @@
tracers = 0.0
do iCell = 1,nCells
! for three layer test
- tracers(index_salinity,1,iCell) = 2.0 ! salinity
- tracers(index_salinity,2,iCell) = 6.0 ! salinity
- tracers(index_salinity,3,iCell) = 13.0 ! salinity
+ !tracers(index_salinity,1,iCell) = 2.0 ! salinity
+ !tracers(index_salinity,2,iCell) = 6.0 ! salinity
+ !tracers(index_salinity,3,iCell) = 13.0 ! salinity
do iLevel = 1,nVertLevels
- !if (xCell(iCell)<500*1e3) tracers(4,iLevel,iCell) = 1.0
+! if (xCell(iCell)<500*1e3.and.ycell(iCell)<1000*1e3) then
+! tracers(index_tracer1,iLevel,iCell) = 1.0
+! tracers(index_tracer2,iLevel,iCell) = 1.0e4
+! endif
+! if (ycell(iCell)>1400*1e3.and.ycell(iCell)<1500*1e3) then
+! tracers(index_tracer1,iLevel,iCell) = 1.0
+! tracers(index_tracer2,iLevel,iCell) = 1.0e4
+! endif
! tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
! tracers(index_temperature,iLevel,iCell) = 6.0-iLevel ! temperature
! tracers(index_temperature,iLevel,iCell) = 6.0-iLevel+yCell(iCell)/4000.e3 ! temperature
! tracers(index_salinity,iLevel,iCell) = 35.0 ! salinity
-
- ! for three layer test
+ ! for 20 layer test
+ tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
! tracers(index_temperature,iLevel,iCell) = &
@@ -149,6 +156,7 @@
tracers(index_tracer1,iLevel,iCell) = 1.0
tracers(index_tracer2,iLevel,iCell) = &
(yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
+! tracers(index_tracer1,iLevel,iCell) = mod(ilevel,2)
rho(iLevel,iCell) = 1000.0*( 1.0 &
- 2.5e-4*tracers(index_temperature,iLevel,iCell) &
+ 7.6e-4*tracers(index_salinity,iLevel,iCell))
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-03 20:22:06 UTC (rev 235)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-03 20:35:52 UTC (rev 236)
@@ -534,7 +534,7 @@
else
dist = zBotEdge(k-1,iEdge) - zBotEdge(k,iEdge)
endif
- tend_u(k,iEdge) = tend_u(k,iEdge) - (fluxVert(k-1) - fluxVert(k))/dist
+ tend_u(k,iEdge) = tend_u(k,iEdge) + (fluxVert(k-1) - fluxVert(k))/dist
enddo
end do
@@ -587,33 +587,41 @@
integer :: iCell, iEdge, k, iTracer, cell1, cell2, &
nEdges, nCells, nVertLevels
- real (kind=RKIND) :: flux, tracer_edge, &
+ real (kind=RKIND) :: flux, tracer_edge, r, &
Tmid(num_tracers,grid % nVertLevels)
+ real (kind=RKIND) :: visc, dist, ah
real (kind=RKIND), dimension(:), pointer :: &
- h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zsurface
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,Fv, h_edge
+ u,h,Fv, h_edge, zmid, zbot
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
integer, dimension(:,:), pointer :: cellsOnEdge
+ real (kind=RKIND), dimension(:,:), allocatable:: fluxVert
+ real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div
u => s % u % array
h => s % h % array
tracers => s % tracers % array
h_edge => s % h_edge % array
Fv => tend % w % array
+ zmid => s % zmid % array
+ zbot => s % zbot % array
+ zsurface => s % zsurface % array
tend_tr => tend % tracers % array
areaCell => grid % areaCell % array
cellsOnEdge => grid % cellsOnEdge % array
dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
nEdges = grid % nEdges
nCells = grid % nCells
nVertLevels = grid % nVertLevels
+ ! Horizontal advection of tracers
tend_tr(:,:,:) = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -632,6 +640,7 @@
end if
end do
+
! mrp 100409 computation of h tendency was, only had div.(huT) term
!do iCell=1,grid % nCellsSolve
! do k=1,grid % nVertLevelsSolve
@@ -641,6 +650,7 @@
! end do
!end do
! mrp 100409 replace with
+ allocate(fluxVert(num_tracers,0:nVertLevels))
do iCell=1,grid % nCellsSolve
! Surface tracer Tmid(:,1) is not used
@@ -683,9 +693,97 @@
end do
end do
- end do
+ ! vertical mixing
+ fluxVert(:,0) = 0.0
+ fluxVert(:,nVertLevels) = 0.0
+ do k=nVertLevels-1,1,-1
+ do iTracer=1,num_tracers
+ ! compute dphi/dz
+ fluxVert(iTracer,k) = ( tracers(iTracer,k ,iCell) - tracers(iTracer,k+1,iCell) )&
+ / (zMid(k,iCell) -zMid(k+1,iCell))
+ enddo
+ enddo
+ ! note vertical diffusivity is hardwired as 2.5e-4 here.
+ ! In POP, default vertical tracer diffusion is 0.25cm^2/s=2.5e-5m^2/s
+ fluxVert = 2.5e-4 * fluxVert
+ k=1
+ dist = zSurface(iCell) - zBot(k,iCell)
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ + (fluxVert(iTracer,k-1) - fluxVert(iTracer,k))/dist
+ enddo
+
+ do k=2,nVertLevels
+ dist = zBot(k-1,iCell) - zBot(k,iCell)
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ + (fluxVert(iTracer,k-1) - fluxVert(iTracer,k))/dist
+ enddo
+ enddo
+
+ enddo ! iCell loop
+ deallocate(fluxVert)
! mrp 100409 end
+ ! mrp 100501: Horizontal tracer diffusion
+ ! first compute kappa*grad(phi) at edges.
+ ! note this could be combined with the above loops for tracer flux.
+ ! leave it separate for now for clarity.
+ ! Hardwire ah as 2e3m^2/s, from 2e7cm^2/s in the POP Users Guide p. 39
+ ah = 2.0e3 ! original
+ allocate(tr_flux(num_tracers,nVertLevels,nEdges))
+! do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+
+ do k=1,nVertLevels
+ do iTracer=1,num_tracers
+ tr_flux(iTracer,k,iEdge) = ah * &
+ (Tracers(iTracer,k,cell2) - Tracers(iTracer,k,cell1))/dcEdge(iEdge)
+ enddo
+ enddo
+ endif
+ enddo
+
+ ! Compute the divergence, div(k grad(phi)) at each cell center
+ !
+ allocate(tr_div(num_tracers,nVertLevels,nCells))
+ tr_div(:,:,:) = 0.0
+! do iEdge=1,grid % nEdgesSolve
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0) then
+ do k=1,nVertLevels
+ do iTracer=1,num_tracers
+ tr_div(iTracer,k,cell1) = tr_div(iTracer,k,cell1) + tr_flux(iTracer,k,iEdge)*dvEdge(iEdge)
+ enddo
+ enddo
+ endif
+ if(cell2 > 0) then
+ do k=1,nVertLevels
+ do iTracer=1,num_tracers
+ tr_div(iTracer,k,cell2) = tr_div(iTracer,k,cell2) - tr_flux(iTracer,k,iEdge)*dvEdge(iEdge)
+ enddo
+ enddo
+ end if
+ end do
+
+ ! The term div(k grad(phi)) is then added to each cell center below.
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,nVertLevels
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ + tr_div(iTracer,k,iCell)*r
+ enddo
+ enddo
+ enddo
+ deallocate(tr_flux, tr_div)
+ ! mrp 100501 end: Horizontal tracer diffusion
+
end subroutine compute_scalar_tend
@@ -1017,7 +1115,7 @@
do k = 1,nVertLevels
if(uBC(k,iEdge).eq.1) then
h1 = 0.0
- h2 = 0.0
+ h2 = 0.0
if(cell1.gt.0) h1 = h(k,cell1)
if(cell2.gt.0) h2 = h(k,cell2)
pv_edge(k,iEdge) = fEdge(iEdge) / ( max(h1,h2) )
</font>
</pre>