<p><b>mpetersen@lanl.gov</b> 2010-06-01 11:17:03 -0600 (Tue, 01 Jun 2010)</p><p>Changes to compare with POP run: add quadratic bottom drag, SSH variable, and some clean up.<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-06-01 15:18:45 UTC (rev 318)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/Registry 2010-06-01 17:17:03 UTC (rev 319)
@@ -125,6 +125,7 @@
var real pv_edge ( nVertLevels nEdges Time ) o pv_edge - -
var real h_edge ( nVertLevels nEdges Time ) o h_edge - -
var real ke ( nVertLevels nCells Time ) o ke - -
+var real ke_edge ( nVertLevels nEdges Time ) o ke_edge - -
var real pv_vertex ( nVertLevels nVertices Time ) o pv_vertex - -
var real pv_cell ( nVertLevels nCells Time ) o pv_cell - -
var real uReconstructX ( nVertLevels nCells Time ) o uReconstructX - -
@@ -139,6 +140,7 @@
var real pZLevel ( nVertLevels nCells Time ) o pZLevel - -
var real MontPot ( nVertLevels nCells Time ) o MontPot - -
var real wTop ( nVertLevelsP1 nCells Time ) o wTop - -
+var real ssh ( nCells Time ) o ssh - -
# Other diagnostic variables: neither read nor written to any files
var real vh ( nVertLevels nEdges Time ) - vh - -
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-06-01 15:18:45 UTC (rev 318)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_global_diagnostics.F 2010-06-01 17:17:03 UTC (rev 319)
@@ -32,12 +32,14 @@
integer, intent(in) :: timeIndex
real (kind=RKIND), intent(in) :: dt
- integer :: nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, nCellsGlobal, nEdgesGlobal, nVerticesGlobal
+ integer :: nVertLevels, nCellsSolve, nEdgesSolve, nVerticesSolve, nCellsGlobal, nEdgesGlobal, nVerticesGlobal, iTracer
+
real (kind=RKIND) :: areaCellGlobal, areaEdgeGlobal, areaTriangleGlobal
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, areaEdge
real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, circulation, vorticity, ke, pv_edge, pv_vertex, &
- pv_cell, gradPVn, gradPVt, zMid, zTop, p, pTop, MontPot, wTop, rho
+ pv_cell, gradPVn, gradPVt, zMid, zTop, p, pTop, MontPot, wTop, rho, tracerTemp
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+
real (kind=RKIND) :: volumeCellGlobal, volumeEdgeGlobal, CFLNumberGlobal
real (kind=RKIND) :: localCFL, localSum
integer :: elementIndex, variableIndex, nVariables, nSums, nMaxes, nMins
@@ -188,6 +190,17 @@
wTop(:,1:nCellsSolve), sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
verticalSumMaxes(variableIndex))
+ ! Tracers
+ allocate(tracerTemp(nVertLevels,nCellsSolve))
+ do iTracer=1,num_tracers
+ variableIndex = variableIndex + 1
+ tracerTemp = Tracers(iTracer,:,1:nCellsSolve)
+ call computeFieldVolumeWeightedLocalStats(dminfo, nVertLevels, nCellsSolve, areaCell(1:nCellsSolve), h(:,1:nCellsSolve), &
+ tracerTemp, sums(variableIndex), mins(variableIndex), maxes(variableIndex), verticalSumMins(variableIndex), &
+ verticalSumMaxes(variableIndex))
+ enddo
+ deallocate(tracerTemp)
+
nVariables = variableIndex
nSums = nVariables
nMins = nVariables
@@ -318,6 +331,12 @@
variableIndex = variableIndex + 1
averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
+ ! Tracers
+ do iTracer=1,num_tracers
+ variableIndex = variableIndex + 1
+ averages(variableIndex) = sums(variableIndex)/volumeCellGlobal
+ enddo
+
! write out the data to files
if (dminfo % my_proc_id == IO_NODE) then
fileID = getFreeUnit()
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-06-01 15:18:45 UTC (rev 318)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_test_cases.F 2010-06-01 17:17:03 UTC (rev 319)
@@ -144,12 +144,12 @@
! 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_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
+ 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) &
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-06-01 15:18:45 UTC (rev 318)
+++ branches/ocean_projects/z_level_mrp/mpas/src/core_ocean/module_time_integration.F 2010-06-01 17:17:03 UTC (rev 319)
@@ -261,10 +261,10 @@
upstream_bias, wTopEdge, rho0Inv
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- zMidZLevel
+ zMidZLevel, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, wTop, &
- tend_h, tend_u, circulation, vorticity, ke, pv_edge, &
+ tend_h, tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
divergence, MontPot, pZLevel, zMidEdge, zTopEdge
type (dm_info) :: dminfo
@@ -273,7 +273,7 @@
cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
edgesOnEdge, edgesOnVertex
real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge, vertVisc
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge, vertViscTop
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
@@ -287,6 +287,7 @@
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
+ ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
MontPot => s % MontPot % array
pZLevel => s % pZLevel % array
@@ -311,6 +312,7 @@
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
tend_h => tend % h % array
tend_u => tend % u % array
@@ -464,21 +466,31 @@
do iEdge=1,grid % nEdgesSolve
! forcing in top layer only
- tend_u(1,iEdge) = tend_u(1,iEdge) + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
+ tend_u(1,iEdge) = tend_u(1,iEdge) &
+ + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
- ! bottom drag - change from nVertLevels to KMT later.
- ! du/dt = u/tau where tau=100 days.
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
- - u(nVertLevels,iEdge)/(100.0*86400.0)
+ - 1.0e-3*u(nVertLevels,iEdge) &
+ *sqrt(2.0*ke_edge(nVertLevels,iEdge))
+! I do not trust the v yet, need to test this one:
+! *sqrt(u(nVertLevels,iEdge)**2 + v(nVertLevels,iEdge)**2)
+ ! old bottom drag, just linear friction
+ ! du/dt = u/tau where tau=100 days.
+ !tend_u(nVertLevels,iEdge) = tend_u(nVertLevels,iEdge) &
+ ! - u(nVertLevels,iEdge)/(100.0*86400.0)
+
enddo
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
- allocate(vertVisc(nVertLevels))
+ allocate(vertViscTop(nVertLevels+1))
if (config_vert_visc_type.eq.'const') then
- vertVisc = config_vert_viscosity
+ vertViscTop = config_vert_viscosity
elseif (config_vert_visc_type.eq.'tanh') then
if (config_vert_grid_type.ne.'zlevel') then
write(0,*) 'Abort: config_vert_visc_type.eq.tanh may only', &
@@ -486,9 +498,9 @@
call dmpar_abort(dminfo)
endif
- do k=1,nVertLevels
- vertVisc(k) = -(config_vmixTanhViscMax-config_vmixTanhViscMin)/2.0 &
- *tanh(-(zMidZLevel(k)-config_vmixTanhZMid) &
+ do k=1,nVertLevels+1
+ vertViscTop(k) = -(config_vmixTanhViscMax-config_vmixTanhViscMin)/2.0 &
+ *tanh(-(zTopZLevel(k)-config_vmixTanhZMid) &
/config_vmixTanhZWidth) &
+ (config_vmixTanhViscMax+config_vmixTanhViscMin)/2
enddo
@@ -502,7 +514,7 @@
fluxVertTop(nVertLevels+1) = 0.0
do iEdge=1,grid % nEdgesSolve
do k=2,nVertLevels
- fluxVertTop(k) = vertVisc(k) &
+ fluxVertTop(k) = vertViscTop(k) &
* ( u(k-1,iEdge) - u(k,iEdge) ) &
/ (zMidEdge(k-1,iEdge) - zMidEdge(k,iEdge))
enddo
@@ -512,10 +524,16 @@
/(zTopEdge(k,iEdge) - zTopEdge(k+1,iEdge))
enddo
end do
- deallocate(fluxVertTop, vertVisc)
+ deallocate(fluxVertTop, vertViscTop)
+!print *, 'u_src'
! do k=1,nVertLevels
! print '(i5,10es10.2)', &
+! k,minval(u_src(k,:)),maxval(u_src(k,:))
+! enddo
+
+! do k=1,nVertLevels
+! print '(i5,10es10.2)', &
! k,minval(tend_h(k,1:grid % ncellssolve)),maxval(tend_h(k,1:grid % ncellssolve))
! enddo
@@ -739,7 +757,7 @@
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
- allocate(vertDiffTop(nVertLevels))
+ allocate(vertDiffTop(nVertLevels+1))
if (config_vert_diff_type.eq.'const') then
vertDiffTop = config_vert_diffusion
elseif (config_vert_diff_type.eq.'tanh') then
@@ -749,7 +767,7 @@
call dmpar_abort(dminfo)
endif
- do k=1,nVertLevels
+ do k=1,nVertLevels+1
vertDiffTop(k) = -(config_vmixTanhDiffMax-config_vmixTanhDiffMin)/2.0 &
*tanh(-(zTopZLevel(k)-config_vmixTanhZMid) &
/config_vmixTanhZWidth) &
@@ -816,17 +834,17 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p, pTop_temp
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, delta_p
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
- hZLevel
+ hZLevel, ssh
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, w, wTop, &
- circulation, vorticity, ke, &
+ circulation, vorticity, ke, ke_edge, &
pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
zMid, zTop, zMidEdge, zTopEdge, p, pTop, MontPot, rho, temperature, salinity, pZLevel
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -847,6 +865,7 @@
vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
+ ke_edge => s % ke_edge % array
pv_edge => s % pv_edge % array
pv_vertex => s % pv_vertex % array
pv_cell => s % pv_cell % array
@@ -862,6 +881,7 @@
pZLevel => s % pZLevel % array
pTop => s % pTop % array
MontPot => s % MontPot % array
+ ssh => s % ssh % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -982,6 +1002,7 @@
end do
!
+ !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
@@ -997,6 +1018,34 @@
end do
!
+ ! Compute ke on cell edges at velocity locations for quadratic bottom drag.
+ ! I did not trust this method, so I averaged from centers instead.
+ ! See next section.
+ !
+ !do iEdge=1,nEdges
+ ! do k=1,nVertLevels
+ ! ke_edge(k,iEdge) = 0.5*(u(k,iEdge)**2 + v(k,iEdge)**2 )
+ ! end do
+ !end do
+
+ !
+ ! Compute ke on cell edges at velocity locations for quadratic bottom drag.
+ !
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 <= nCells .and. cell2 <= nCells) then
+ do k=1,nVertLevels
+ ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
+ end do
+ else
+ do k=1,nVertLevels
+ ke_edge(k,iEdge) = 0.0
+ end do
+ end if
+ end do
+
+ !
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
@@ -1093,6 +1142,13 @@
enddo
!
+ ! Compute sea surface height
+ !
+ do iCell=1,nCells
+ ssh(iCell) = h(1,iCell) - hZLevel(1)
+ enddo
+
+ !
! equation of state
!
! For a isopycnal model, density should remain constant.
@@ -1166,11 +1222,9 @@
! assume atmospheric pressure at the surface is zero for now.
pZLevel(1,iCell) = rho(1,iCell)*gravity &
* (h(1,iCell)-0.5*hZLevel(1))
- pTop_temp = rho(1,iCell)*gravity*h(1,iCell)
do k=2,nVertLevels
delta_p = rho(k,iCell)*gravity*hZLevel(k)
pZLevel(k,iCell) = pZLevel(k-1,iCell) + 0.5*delta_p
- pTop_temp = pTop_temp + delta_p
end do
end do
</font>
</pre>