<p><b>mpetersen@lanl.gov</b> 2010-05-13 15:55:22 -0600 (Thu, 13 May 2010)</p><p>Clean up comments in tendency subroutines.<br>
</p><hr noshade><pre><font color="gray">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 19:32:10 UTC (rev 269)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-13 21:55:22 UTC (rev 270)
@@ -284,7 +284,6 @@
u => s % u % array
v => s % v % array
wBot => s % wBot % array
- hFluxBot => tend % wBot % array
h_edge => s % h_edge % array
circulation => s % circulation % array
vorticity => s % vorticity % array
@@ -318,6 +317,7 @@
tend_h => tend % h % array
tend_u => tend % u % array
+ hFluxBot => tend % wBot % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -327,9 +327,10 @@
u_src => grid % u_src % array
!
- ! Compute height tendency for each cell
- ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+ ! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
!
+ ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+ ! for explanation of divergence operator.
tend_h(:,:) = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
@@ -347,13 +348,15 @@
end do
end if
end do
-
do iCell=1,nCells
do k=1,nVertLevels
tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
end do
end do
+ !
+ ! height tendency: vertical advection term -d/dz(hw)
+ !
if (config_vert_grid_type.eq.'isopycnal') then
hFluxBot=0.0 ! vertical fluxes are zero
@@ -387,42 +390,23 @@
end do
endif ! coordinate type
- ! mrp 100409 end
-!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(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(hFluxBot(k,1:grid % ncellssolve)),maxval(hFluxBot(k,1:grid % ncellssolve))
-! enddo
-
!
- ! Compute u (normal) velocity tendency for each edge (cell face)
+ ! velocity tendency: vertical advection term -w du/dz
!
tend_u(:,:) = 0.0
- rho0Inv = 1.0/config_rho0
- do iEdge=1,grid % nEdgesSolve
+ if (config_vert_grid_type.eq.'zlevel') then
+ do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- ! mrp 100422 add -w*dudz term to tendancy
- ! For isopycnal mode, wBot should be zero.
- ! change nvertlevels to kmt later
do k=1,nVertLevels-1
! Average w from cell center to edge
wBotEdge = 0.5*(wBot(k,cell1)+wBot(k,cell2))
! compute dudz at vertical interface with first order derivative.
-! mrp 100511 new:
-! w_dudzBotEdge(k) = wBotEdge * (u(k,iEdge)-u(k+1,iEdge)) &
-! / (zMidZLevel(k) - zMidZLevel(k+1))
-! mrp 100511 old:
w_dudzBotEdge(k) = wBotEdge * (u(k,iEdge)-u(k+1,iEdge)) &
- * 2.0 / (h_edge(k,iEdge) + h_edge(k+1,iEdge))
+ / (zMidZLevel(k) - zMidZLevel(k+1))
end do
w_dudzBotEdge(nVertLevels) = 0.0
@@ -431,29 +415,48 @@
do k=2,nVertLevels
tend_u(k,iEdge) = - 0.5 * (w_dudzBotEdge(k-1) + w_dudzBotEdge(k))
enddo
- ! mrp 100422 end: add -w*dudz term to tendancy
+ enddo
+ endif
- ! mrp 100426: add pressure gradient
- if (config_vert_grid_type.eq.'isopycnal') then
- do k=1,nVertLevels
+ !
+ ! velocity tendency: pressure gradient
+ !
+ rho0Inv = 1.0/config_rho0
+ if (config_vert_grid_type.eq.'isopycnal') then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
tend_u(k,iEdge) = tend_u(k,iEdge) &
- (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
end do
- elseif (config_vert_grid_type.eq.'zlevel') then
- do k=1,nVertLevels
+ enddo
+ elseif (config_vert_grid_type.eq.'zlevel') then
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
tend_u(k,iEdge) = tend_u(k,iEdge) &
- rho0Inv*( pmidZLevel(k,cell2) &
- pmidZLevel(k,cell1) )/dcEdge(iEdge)
- end do
- endif
- ! mrp 100426 end: add pressure gradient
+ end do
+ enddo
+ endif
+ !
+ ! velocity tendency: -q(h u^\perp) - </font>
<font color="blue">abla K
+ ! +</font>
<font color="black">u_h(</font>
<font color="black">abla \delta + {\bf k}\times </font>
<font color="blue">abla \xi)
+ !
+ ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="red">abla vorticity
+ ! only valid for visc == constant
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
do k=1,nVertLevels
- !
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
- ! only valid for visc == constant
- !
u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
u_diffusion = visc * u_diffusion
@@ -472,6 +475,9 @@
end do
end do
+ !
+ ! velocity tendency: forcing and bottom drag
+ !
do iEdge=1,grid % nEdgesSolve
! forcing in top layer only
@@ -484,13 +490,14 @@
enddo
- ! vertical mixing
+ !
+ ! velocity tendency: vertical mixing d/dz( nu_v du/dz))
+ !
allocate(fluxVert(0:nVertLevels))
do iEdge=1,grid % nEdgesSolve
fluxVert(0) = 0.0
fluxVert(nVertLevels) = 0.0
do k=nVertLevels-1,1,-1
- ! default vertical viscosity is 2.5e-4m^2/s
fluxVert(k) = config_vert_viscosity &
* ( u(k,iEdge) - u(k+1,iEdge) ) &
/ (zMidEdge(k,iEdge) -zMidEdge(k+1,iEdge))
@@ -508,7 +515,16 @@
end do
deallocate(fluxVert)
+!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(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(hFluxBot(k,1:grid % ncellssolve)),maxval(hFluxBot(k,1:grid % ncellssolve))
+! enddo
+
+! print '(a,i5,f10.5)', &
! 'k, min(tend_u(k,:)),max(tend_u(k,:))'
! do k=1,nVertLevels
! print '(i5,10es22.12)', &
@@ -535,8 +551,7 @@
integer :: iCell, iEdge, k, iTracer, cell1, cell2, upwindCell,&
nEdges, nCells, nVertLevels
- real (kind=RKIND) :: flux, tracer_edge, r, &
- Tmid(num_tracers,grid % nVertLevels)
+ real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND) :: dist
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zsurface
@@ -546,7 +561,8 @@
tracers, tend_tr
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension(:,:), allocatable:: fluxVert
+ real (kind=RKIND), dimension(:), allocatable:: hFluxBotCol
+ real (kind=RKIND), dimension(:,:), allocatable:: fluxVert, tracerBot
real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div
u => s % u % array
@@ -568,7 +584,9 @@
nCells = grid % nCells
nVertLevels = grid % nVertLevels
- ! Horizontal advection of tracers
+ !
+ ! tracer tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( h \phi u)
+ !
tend_tr(:,:,:) = 0.0
if (config_hor_tracer_adv.eq.'centered') then
@@ -612,104 +630,74 @@
end do
endif
+ do iCell=1,grid % nCellsSolve
+ do k=1,grid % nVertLevelsSolve
+ do iTracer=1,num_tracers
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) / areaCell(iCell)
+ end do
+ end do
+ end do
- ! mrp 100409 computation of h tendency was, only had div.(huT) term
- !do iCell=1,grid % nCellsSolve
- ! do k=1,grid % nVertLevelsSolve
- ! do iTracer=1,num_tracers
- ! tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) / areaCell(iCell)
- ! end do
- ! end do
- !end do
- ! mrp 100409 replace with
- allocate(fluxVert(num_tracers,0:nVertLevels))
- Tmid=0.0
+ !
+ ! tracer tendency: vertical advection term -d/dz( h \phi w)
+ !
+ allocate(tracerBot(num_tracers,0:nVertLevels))
+ allocate(hFluxBotCol(0:nVertLevels))
+ tracerBot(:,0)=0.0
+ tracerBot(:,nVertLevels)=0.0
+ hFluxBotCol(0)=0.0
+ hFluxBotCol(nVertLevels)=0.0
do iCell=1,grid % nCellsSolve
- ! Tmid(:,k) is the tracer value at the interface between layers k and k+1
- ! Surface tracer Tmid(:,1) is not used
- ! For now, Tmid is just the average of the upper and lower points.
- ! Can use a fancier interpolation later.
- ! change nVertLevels to KMT later
if (config_vert_tracer_adv.eq.'centered') then
-
- do k=1,nVertLevels-1
- do iTracer=1,num_tracers
- Tmid(iTracer,k) = ( tracers(iTracer,k ,iCell) &
+ do k=1,nVertLevels-1
+ do iTracer=1,num_tracers
+ tracerBot(iTracer,k) = ( tracers(iTracer,k ,iCell) &
+ tracers(iTracer,k+1,iCell))/2.0
- end do
- end do
+ end do
+ end do
elseif (config_vert_tracer_adv.eq.'upwind') then
+ do k=1,nVertLevels-1
+ if (hFluxBot(k,iCell)>0.0) then
+ upwindCell = k+1
+ else
+ upwindCell = k
+ endif
+ do iTracer=1,num_tracers
+ tracerBot(iTracer,k) = tracers(iTracer,upwindCell,iCell)
+ end do
+ end do
- do k=1,nVertLevels-1
- if (hFluxBot(k,iCell)>0.0) then
- upwindCell = k+1
- else
- upwindCell = k
- endif
- do iTracer=1,num_tracers
- Tmid(iTracer,k) = tracers(iTracer,upwindCell,iCell)
- end do
- end do
-
- elseif (config_vert_tracer_adv.eq.'downwind') then
-
- do k=1,nVertLevels-1
- if (hFluxBot(k,iCell)>0.0) then
- upwindCell = k
- else
- upwindCell = k+1
- endif
- do iTracer=1,num_tracers
- Tmid(iTracer,k) = tracers(iTracer,upwindCell,iCell)
- end do
- end do
-
endif
- ! The next loop could be combined with the later two for better efficiency.
- ! Keep separate now for clarity.
- do k=1,grid % nVertLevelsSolve
+ hFluxBotCol(1:nVertLevels)=hFluxBot(1:nVertLevels,iCell)
+ do k=1,nVertLevels
do iTracer=1,num_tracers
- ! this is just div.(huT) term:
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) / areaCell(iCell)
- end do
- end do
-
- ! top layer: flux from below only
- k=1
- do iTracer=1,num_tracers
- ! this is -div.(huT) - d/dz(hwT) terms.
- ! 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) &
- - ( - hFluxBot(k,iCell)*Tmid(iTracer,k))/h(k,icell)
+ - ( hFluxBotCol(k-1)*tracerBot(iTracer,k-1) &
+ - hFluxBotCol(k )*tracerBot(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 hFluxBot(1,iCell)=0 so index can go to k=1.
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( hFluxBot(k-1,iCell)*Tmid(iTracer,k-1) &
- - hFluxBot(k ,iCell)*Tmid(iTracer,k ))/h(k,icell)
- end do
end do
- ! vertical mixing
- fluxVert(:,0) = 0.0
- fluxVert(:,nVertLevels) = 0.0
+ enddo ! iCell
+ deallocate(tracerBot,hFluxBotcol)
+
+ !
+ ! tracer tendency: vertical mixing d/dz( \kappa_v d\phi/dz))
+ !
+ allocate(fluxVert(num_tracers,0:nVertLevels))
+ fluxVert(:,0) = 0.0
+ fluxVert(:,nVertLevels) = 0.0
+ do iCell=1,grid % nCellsSolve
do k=1,nVertLevels-1
do iTracer=1,num_tracers
- ! compute kappa_v dphi/dz
- ! Default vertical diffusivity is 2.5e-5m^2/s.
- ! In POP, default vertical tracer diffusion is 0.25cm^2/s=2.5e-5m^2/s
fluxVert(iTracer,k) = config_vert_diffusion &
* (tracers(iTracer,k,iCell) - tracers(iTracer,k+1,iCell) )&
/ (zMid(k,iCell) -zMid(k+1,iCell))
enddo
enddo
+
k=1
dist = zSurface(iCell) - zBot(k,iCell)
do iTracer=1,num_tracers
@@ -727,31 +715,28 @@
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.
- ! default kappa is 2e3m^2/s, from 2e7cm^2/s in the POP Users Guide p. 39
+ !
+ ! tracer tendency: horizontal tracer diffusion
+ ! </font>
<font color="black">abla\cdot (\kappa_h </font>
<font color="blue">abla\phi )
+ !
+ ! first compute \kappa_h </font>
<font color="red">abla\phi at horizontal edges.
allocate(tr_flux(num_tracers,nVertLevels,nEdges))
tr_flux(:,:,:) = 0.0
do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells .and. cell2 <= nCells) then
-
- do k=1,nVertLevels
- do iTracer=1,num_tracers
- tr_flux(iTracer,k,iEdge) = config_hor_diffusion * &
- (Tracers(iTracer,k,cell2) - Tracers(iTracer,k,cell1))/dcEdge(iEdge)
- enddo
- enddo
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,nVertLevels
+ do iTracer=1,num_tracers
+ tr_flux(iTracer,k,iEdge) = config_hor_diffusion * &
+ (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
- !
+ ! Compute the divergence, div(\kappa_h </font>
<font color="gray">abla\phi) at cell centers
allocate(tr_div(num_tracers,nVertLevels,nCells))
tr_div(:,:,:) = 0.0
do iEdge=1,nEdges
@@ -775,8 +760,8 @@
end if
end do
- ! The term div(k grad(phi)) is then added to each cell center below.
- do iCell = 1,nCells
+ ! add </font>
<font color="black">abla\cdot (\kappa_h </font>
<font color="gray">abla\phi ) to tracer tendency
+ do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
do k = 1,nVertLevels
do iTracer=1,num_tracers
@@ -786,7 +771,6 @@
enddo
enddo
deallocate(tr_flux, tr_div)
- ! mrp 100501 end: Horizontal tracer diffusion
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
@@ -1102,18 +1086,15 @@
!
! For a isopycnal model, density should remain constant.
if (config_vert_grid_type.eq.'zlevel') then
-
- do iCell=1,nCells
- do k=1,nVertLevels
+ 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*tracers(index_temperature,k,iCell) &
+ 7.6e-4*tracers(index_salinity,k,iCell))
- end do
- end do
-
+ end do
+ end do
endif
- !mrp 100324 end
! For Isopycnal model.
</font>
</pre>