<p><b>mpetersen@lanl.gov</b> 2010-05-04 14:03:20 -0600 (Tue, 04 May 2010)</p><p>Added input namelist variables for horizontal and vertical mixing coefficients.<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-04 15:59:06 UTC (rev 240)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-05-04 20:03:20 UTC (rev 241)
@@ -14,11 +14,13 @@
namelist integer restart config_restart_interval 0
namelist logical restart config_do_restart false
namelist real restart config_restart_time 172800.0
-# config_vert_grid_type: zlevel or isopycnal
namelist character grid config_vert_grid_type isopycnal
-# mrp 100426 added rho0, only used in zlevel pgrad term.
namelist real grid config_rho0 1028
+namelist real hmix config_hor_diffusion 2000.0
+namelist real vmix config_vert_viscosity 2.5e-4
+namelist real vmix config_vert_diffusion 2.5e-5
+
#
# dim type name_in_file name_in_code
#
Modified: branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F
===================================================================
--- branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F 2010-05-04 15:59:06 UTC (rev 240)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F 2010-05-04 20:03:20 UTC (rev 241)
@@ -41,7 +41,7 @@
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
real (kind=RKIND) :: localCFL, localSum
integer :: elementIndex, variableIndex, nVariables, nSums, nMaxes, nMins
- integer :: timeLevel,k
+ integer :: timeLevel,k,i
integer, parameter :: kMaxVariables = 1024 ! this must be a little more than double the number of variables to be reduced
@@ -369,6 +369,16 @@
print '(a,i5,10es25.15)', 'tracer min max',&
k,minval(tracers(k,:,1:nCellsSolve)),maxval(tracers(k,:,1:nCellsSolve))
enddo
+
+ ! print some diagnostics
+ !print '(10a)', 'itracer ilevel min tracer max tracer'
+ !do i=1,num_tracers
+ !do k = 1,nVertLevels
+ ! print '(2i5,20es12.4)', i,k, &
+ ! minval(tracers(i,k,:)), maxval(tracers(i,k,:))
+ !enddo
+ !enddo
+
! mrp 100415 temp output end
end subroutine computeGlobalDiagnostics
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-04 15:59:06 UTC (rev 240)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-05-04 20:03:20 UTC (rev 241)
@@ -21,7 +21,7 @@
type (domain_type), intent(inout) :: domain
- integer :: i, iCell, iEdge, iVtx, iLevel
+ integer :: i, iCell, iEdge, iVtx, iLevel, iTracer
type (block_type), pointer :: block_ptr
real (kind=RKIND), dimension(:), pointer :: xCell,yCell, &
hZLevel, zmidZLevel, zbotZLevel
@@ -129,9 +129,9 @@
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.and.ycell(iCell)<1000*1e3) then
! tracers(index_tracer1,iLevel,iCell) = 1.0
@@ -145,6 +145,7 @@
! 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 20 layer test
tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
@@ -212,8 +213,6 @@
print '(10a)', 'ilevel',&
' rho ',&
' min u max u ',&
- ' min T max T ',&
- ' min S max S ',&
' min u_src max u_src ', &
' min h max h ',&
' hZLevel zmidZlevel', &
@@ -221,13 +220,20 @@
do iLevel = 1,nVertLevels
print '(i5,20es12.4)', ilevel, rho(ilevel,1), &
minval(u(iLevel,:)), maxval(u(iLevel,:)), &
- minval(tracers(index_temperature,iLevel,:)), maxval(tracers(index_temperature,iLevel,:)), &
- minval(tracers(index_salinity,iLevel,:)), maxval(tracers(index_salinity,iLevel,:)), &
minval(u_src(iLevel,:)), maxval(u_src(iLevel,:)), &
minval(h(iLevel,:)), maxval(h(iLevel,:)), &
hZLevel(iLevel),zmidZlevel(iLevel),zbotZlevel(iLevel)
enddo
+ ! print some diagnostics
+ print '(10a)', 'itracer ilevel min tracer max tracer'
+ do iTracer=1,num_tracers
+ do iLevel = 1,nVertLevels
+ print '(2i5,20es12.4)', iTracer,ilevel, &
+ minval(tracers(itracer,iLevel,:)), maxval(tracers(itracer,iLevel,:))
+ enddo
+ enddo
+
block_ptr => block_ptr % next
end do
! mrp 100121 end
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-04 15:59:06 UTC (rev 240)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-04 20:03:20 UTC (rev 241)
@@ -70,7 +70,7 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
- integer :: iCell, k
+ integer :: iCell, k, i
type (block_type), pointer :: block
integer, parameter :: PROVIS = 1
@@ -161,6 +161,7 @@
if (rk_step < 4) then
block => domain % blocklist
do while (associated(block))
+
block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:) &
+ rk_substep_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
block % intermediate_step(PROVIS) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) &
@@ -523,18 +524,20 @@
fluxVert(0) = 0.0
fluxVert(nVertLevels) = 0.0
do k=nVertLevels-1,1,-1
- fluxVert(k) = ( u(k,iEdge) - u(k+1,iEdge) ) &
+ ! 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))
enddo
- ! note vertical viscosity is hardwired as 2.5e-4 here.
- fluxVert = 2.5e-4 * fluxVert
do k=1,nVertLevels
if(k.eq.1) then
dist = zSurfaceEdge(iEdge) - zBotEdge(k,iEdge)
else
dist = zBotEdge(k-1,iEdge) - zBotEdge(k,iEdge)
endif
+
tend_u(k,iEdge) = tend_u(k,iEdge) + (fluxVert(k-1) - fluxVert(k))/dist
+
enddo
end do
@@ -590,7 +593,7 @@
real (kind=RKIND) :: flux, tracer_edge, r, &
Tmid(num_tracers,grid % nVertLevels)
- real (kind=RKIND) :: visc, dist, ah
+ real (kind=RKIND) :: dist
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, zsurface
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -651,8 +654,10 @@
!end do
! mrp 100409 replace with
allocate(fluxVert(num_tracers,0:nVertLevels))
- do iCell=1,grid % nCellsSolve
+ Tmid=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.
@@ -696,21 +701,21 @@
! vertical mixing
fluxVert(:,0) = 0.0
fluxVert(:,nVertLevels) = 0.0
- do k=nVertLevels-1,1,-1
+ do k=1,nVertLevels-1
do iTracer=1,num_tracers
- ! compute dphi/dz
- fluxVert(iTracer,k) = ( tracers(iTracer,k ,iCell) - tracers(iTracer,k+1,iCell) )&
+ ! 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
- ! 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
+ + (fluxVert(iTracer,k-1) - fluxVert(iTracer,k))/dist ! orig
enddo
do k=2,nVertLevels
@@ -729,10 +734,8 @@
! 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
+ ! default kappa is 2e3m^2/s, from 2e7cm^2/s in the POP Users Guide p. 39
allocate(tr_flux(num_tracers,nVertLevels,nEdges))
-! do iEdge=1,grid % nEdgesSolve
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -740,7 +743,7 @@
do k=1,nVertLevels
do iTracer=1,num_tracers
- tr_flux(iTracer,k,iEdge) = ah * &
+ tr_flux(iTracer,k,iEdge) = config_hor_diffusion * &
(Tracers(iTracer,k,cell2) - Tracers(iTracer,k,cell1))/dcEdge(iEdge)
enddo
enddo
@@ -751,27 +754,28 @@
!
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
+ 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
+ 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.
+ ! 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
@@ -784,6 +788,15 @@
deallocate(tr_flux, tr_div)
! mrp 100501 end: Horizontal tracer diffusion
+ ! print some diagnostics
+ !do iTracer=1,num_tracers
+ !do k = 1,nVertLevels
+ ! print '(2i5,20es12.4)', iTracer,k, &
+ ! minval(tend_tr(itracer,k,:)), maxval(tend_tr(itracer,k,:))
+ !enddo
+ !enddo
+
+
end subroutine compute_scalar_tend
</font>
</pre>