<p><b>mpetersen@lanl.gov</b> 2010-05-21 11:38:43 -0600 (Fri, 21 May 2010)</p><p>Updated vertical advection to simply use vertical velocity w as the vertical thickness flux. This version had a successful 50 day run on the x3 grid.<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-21 17:20:37 UTC (rev 298)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-05-21 17:38:43 UTC (rev 299)
@@ -100,7 +100,10 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
h => block_ptr % time_levs(1) % state % h % array
+ u => block_ptr % time_levs(1) % state % u % array
+ rho => block_ptr % time_levs(1) % state % rho % array
tracers => block_ptr % time_levs(1) % state % tracers % array
+ u_src => block_ptr % mesh % u_src % array
xCell => block_ptr % mesh % xCell % array
yCell => block_ptr % mesh % yCell % array
@@ -132,18 +135,29 @@
call dmpar_abort(dminfo)
endif
+ ! Set tracers, if not done in grid.nc file
!tracers = 0.0
- !do iCell = 1,nCells
- ! do iLevel = 1,nVertLevels
- ! ! for 20 layer test
- ! tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
- ! tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
- ! tracers(index_tracer1,iLevel,iCell) = 1.0
- ! tracers(index_tracer2,iLevel,iCell) = &
- ! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
- ! enddo
- !enddo
+ do iCell = 1,nCells
+ do iLevel = 1,nVertLevels
+ ! for 20 layer test
+ ! tracers(index_temperature,iLevel,iCell) = 5.0 ! temperature
+ ! tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
+ ! for x3, 25 layer test
+ !tracers(index_temperature,iLevel,iCell) = 10.0 ! temperature
+ !tracers(index_salinity,iLevel,iCell) = 1.4 + iLevel*0.6 ! salinity
+
+ ! tracers(index_tracer1,iLevel,iCell) = 1.0
+ ! tracers(index_tracer2,iLevel,iCell) = &
+ ! (yCell(iCell)/4000.e3 + xCell(iCell)/2500.e3 )/2.0
+
+ rho(iLevel,iCell) = 1000.0*( 1.0 &
+ - 2.5e-4*tracers(index_temperature,iLevel,iCell) &
+ + 7.6e-4*tracers(index_salinity,iLevel,iCell))
+
+ enddo
+ enddo
+
endif
! print some diagnostics
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-21 17:20:37 UTC (rev 298)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-05-21 17:38:43 UTC (rev 299)
@@ -264,7 +264,7 @@
zMidZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wTop, &
- tend_h, tend_u, hFluxTop, circulation, vorticity, ke, pv_edge, &
+ tend_h, tend_u, circulation, vorticity, ke, pv_edge, &
divergence, MontPot, pZLevel, zMidEdge, zTopEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
@@ -315,7 +315,6 @@
tend_h => tend % h % array
tend_u => tend % u % array
- hFluxTop => tend % wTop % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -355,30 +354,21 @@
!
! height tendency: vertical advection term -d/dz(hw)
!
- if (config_vert_grid_type.eq.'isopycnal') then
- hFluxTop=0.0 ! vertical fluxes are zero
+ if (config_vert_grid_type.eq.'zlevel') then
- elseif (config_vert_grid_type.eq.'zlevel') then
-
do iCell=1,nCells
- hFluxTop(1,iCell) = 0.0
- hFluxTop(nVertLevels+1,iCell) = 0.0
- ! compute hw, the vertical height flux, using
- ! hw = -integral_zbot^ztop ( div(hu) )
- do k=nVertLevels,2,-1
- ! I use +tend_h because tend_h is -div.(hu)
- hFluxTop(k,iCell) = hFluxTop(k+1,iCell) + tend_h(k,iCell)*h(k,iCell)
- end do
+ tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
- ! add -d/dz(hw) to heigh tendency
- ! note if you substitute this into the line above this is
- ! just tend_h=0 except for the top layer.
- ! So this loop could be changed to just k=1 in the future.
- do k=1,nVertLevels
+ ! This next loop is to verify that h for levels below the first
+ ! remain constant. At a later time this could be replaced
+ ! by just tend_h(2:nVertLevels,:) = 0.0, and then there is
+ ! no need to compute the horizontal tend_h term for k=2:nVertLevels
+ ! on a zlevel grid, above.
+ do k=2,nVertLevels
tend_h(k,iCell) = tend_h(k,iCell) &
- - (hFluxTop(k,iCell) - hFluxTop(k+1,iCell))/h(k,iCell)
- end do
+ - (wTop(k,iCell) - wTop(k+1,iCell))
+ end do
end do
endif ! coordinate type
@@ -504,13 +494,9 @@
end do
deallocate(fluxVertTop)
-!print *, 'ncells,grid % nCellsSolve, size(hFluxTop)',ncells,grid % nCellsSolve, size(hFluxTop,1), size(hFluxTop,2)
-! print '(a,i5,f10.5)', &
-! 'k, min(tend_h(k,2:)),max(tend_h(k,2:)),min(hFluxTop(k,:)),max(hFluxTop(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(hFluxTop(k,1:grid % ncellssolve)),maxval(hFluxTop(k,1:grid % ncellssolve))
+! k,minval(tend_h(k,1:grid % ncellssolve)),maxval(tend_h(k,1:grid % ncellssolve))
! enddo
! print '(a,i5,f10.5)', &
@@ -545,20 +531,19 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,hFluxTop, h_edge, zMid, zTop
+ u,h,wTop, h_edge, zMid, zTop
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension(:), allocatable:: hFluxTopCol
real (kind=RKIND), dimension(:,:), allocatable:: fluxVertTop, tracerTop
real (kind=RKIND), dimension(:,:,:), allocatable::tr_flux, tr_div
u => s % u % array
h => s % h % array
+ wTop => s % wTop % array
tracers => s % tracers % array
h_edge => s % h_edge % array
- hFluxTop => tend % wTop % array
zMid => s % zMid % array
zTop => s % zTop % array
@@ -644,7 +629,7 @@
elseif (config_vert_tracer_adv.eq.'upwind') then
do k=2,nVertLevels
- if (hFluxTop(k,iCell)>0.0) then
+ if (wTop(k,iCell)>0.0) then
upwindCell = k
else
upwindCell = k-1
@@ -659,8 +644,8 @@
do k=1,nVertLevels
do iTracer=1,num_tracers
tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- - ( hFluxTop(k ,iCell)*tracerTop(iTracer,k ) &
- - hFluxTop(k+1,iCell)*tracerTop(iTracer,k+1))/h(k,icell)
+ - ( wTop(k ,iCell)*tracerTop(iTracer,k ) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1))
end do
end do
</font>
</pre>