<p><b>mpetersen@lanl.gov</b> 2012-07-26 14:02:51 -0600 (Thu, 26 Jul 2012)</p><p>BRANCH COMMIT surface_pressure_mrp<br>
<br>
This branch can run successfully with an imposed pressure at the<br>
top surface. The critical parts to add to the algorithm are:<br>
<br>
1. Add the variable surfacePressure to the Registry:<br>
var persistent real surfacePressure ( nCells ) 1 iro surfacePressure mesh - -<br>
<br>
2. Add surfacePressure to the top layer when computing pressure<br>
pressure(1,iCell) = surfacePressure(iCell) + rho(1,iCell)*gravity*0.5*h(1,iCell)<br>
<br>
Nothing else is required. The SSH variable does not include any<br>
height contribution from the surfacePressure. Split explicit still<br>
works because 1/rho0 grad(surfacePressure) ends up in the Gbar term,<br>
so it is added to the barotropic momentum equation to counter the<br>
depressed sea surface height.<br>
<br>
The branch was tested as follows:<br>
<br>
1. Zero initial velocity, constant initial T&S, no wind stress, in a<br>
periodic hex channel 40x20x40 80km gridcells, 1600m deep, with a<br>
tophat 1200km wide in the x-direction in surfacePressure. Tests<br>
were made with peak surfacePressure equivalent to 10%, 50%, 90%<br>
of the total depth. Rho is step 2 above is same as rest of domain.<br>
Use a linear EOS, as JM is depth-dependant. One expects the fluid<br>
to remain at rest here, because the grad(p) and grad(zMid) terms<br>
should balance. Global KE remains below 1e-22 to 90 days using<br>
split explicit for all three cases, and below 1e-19 for RK4.<br>
<br>
2. Same as #1 but with stratification and "wind" forcing on top layer.<br>
Here a velocity develops because the rho in the surface thickness<br>
is not the same as rho of the displaced water column, and surface<br>
"wind" stress. This is simply a stability test. All cases run to<br>
90 days stably for both RK4 and split explicit, except 90% for<br>
split explicit hits a nan with 90%. Split explicit works stably<br>
with 70% or less.<br>
<br>
Because the changes required are so minimal (two lines) I will not<br>
merge this to the trunk, because it would then require a<br>
surfacePressure field to be added to the input files. This branch<br>
also contains the tophat and Gaussian initialization of<br>
surfacePressure in mpas_ocn_mpas_core.F.<br>
</p><hr noshade><pre><font color="gray">Index: branches/ocean_projects/surface_pressure_mrp
===================================================================
--- branches/ocean_projects/surface_pressure_mrp 2012-07-26 18:23:23 UTC (rev 2064)
+++ branches/ocean_projects/surface_pressure_mrp 2012-07-26 20:02:51 UTC (rev 2065)
Property changes on: branches/ocean_projects/surface_pressure_mrp
___________________________________________________________________
Modified: svn:mergeinfo
## -19,3 +19,4 ##
/branches/omp_blocks/io:1639-1787
/branches/source_renaming:1082-1113
/branches/time_manager:924-962
+/trunk/mpas:2049-2055
\ No newline at end of property
Modified: branches/ocean_projects/surface_pressure_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/surface_pressure_mrp/src/core_ocean/Registry 2012-07-26 18:23:23 UTC (rev 2064)
+++ branches/ocean_projects/surface_pressure_mrp/src/core_ocean/Registry 2012-07-26 20:02:51 UTC (rev 2065)
@@ -275,7 +275,7 @@
var persistent real uSrcReconstructMeridional ( nVertLevels nCells Time ) 2 o uSrcReconstructMeridional state - -
var persistent real MontPot ( nVertLevels nCells Time ) 2 - MontPot state - -
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
-var persistent real surfaceThickness ( nCells ) 1 iro surfaceThickness mesh - -
+var persistent real surfacePressure ( nCells ) 1 iro surfacePressure mesh - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
Modified: branches/ocean_projects/surface_pressure_mrp/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/surface_pressure_mrp/src/core_ocean/mpas_ocn_mpas_core.F 2012-07-26 18:23:23 UTC (rev 2064)
+++ branches/ocean_projects/surface_pressure_mrp/src/core_ocean/mpas_ocn_mpas_core.F 2012-07-26 20:02:51 UTC (rev 2065)
@@ -252,7 +252,7 @@
! mrp 120724 initialize surfaceThickness and other variables for surface_pressure branch
! This is only for initial testing. Later on, it should go in the grid file.
- real (kind=RKIND) :: x_mid, x_efold, xCell, totalH
+ real (kind=RKIND) :: x_mid, x_efold, xCell, totalH, surfaceThickness, surfaceRho
! make even layers
totalH = 1600.0
@@ -262,6 +262,12 @@
! turn off the wind.
!block % mesh % u_src % array = 0.0
+ ! surfacePressure is the variable used by the code. Here I
+ ! specify a surface thickness, and then back compute pressure
+ ! with p=rho g h. The rho from the linear equation of state,
+ ! for T=10 and S=35, is:
+ surfaceRho = 1000*(1-2.5e-4*10 + 7.6e-4*35.0)
+
print *, 'config_max_surface_thickness',config_max_surface_thickness
print *, 'config_max_surface_thickness*H',config_max_surface_thickness &
*sum(block % state % time_levs(1) % state % h % array(:,1))
@@ -273,22 +279,23 @@
! if (abs(xCell-x_mid)>3*x_efold) then ! Gaussian
if (abs(xCell-x_mid)>x_efold) then ! tophat
- block % mesh % surfaceThickness % array(iCell) = 0.0
+ surfaceThickness = 0.0
else
! Gaussian
- !block % mesh % surfaceThickness % array(iCell) &
- ! = config_max_surface_thickness*totalH &
+ !surfaceThickness = config_max_surface_thickness*totalH &
! * exp( -((xCell-x_mid)/x_efold)**2 )
! top hat
- block % mesh % surfaceThickness % array(iCell) &
- = config_max_surface_thickness*totalH
+ surfaceThickness = config_max_surface_thickness*totalH
block % state % time_levs(1) % state % h % array(:,iCell) &
- = (1.0 -block % mesh % surfaceThickness % array(iCell)/totalH) &
+ = (1.0 -surfaceThickness/totalH) &
*block % state % time_levs(1) % state % h % array(:,iCell)
endif
+ block % mesh % surfacePressure % array(iCell) = surfaceRho*gravity*surfaceThickness
+
+
! Compute zMid, the z-coordinate of the middle of the layer.
! This is used for the rho g grad z momentum term.
! Note the negative sign, since referenceBottomDepth is positive
@@ -311,7 +318,7 @@
10.0 + block % state % time_levs(1) % state % zMid % array(k,iCell)/totalH*10.0
end do
- print *, 'xCell/1e3,surfaceThickness',xCell/1.0e3,block % mesh % surfaceThickness % array(iCell)
+ print *, 'xCell/1e3,surfaceThickness',xCell/1.0e3,surfaceThickness
enddo
Modified: branches/ocean_projects/surface_pressure_mrp/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/surface_pressure_mrp/src/core_ocean/mpas_ocn_tendency.F 2012-07-26 18:23:23 UTC (rev 2064)
+++ branches/ocean_projects/surface_pressure_mrp/src/core_ocean/mpas_ocn_tendency.F 2012-07-26 20:02:51 UTC (rev 2065)
@@ -412,7 +412,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, &
- referenceBottomDepth, ssh, surfaceThickness
+ referenceBottomDepth, ssh, surfacePressure
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -468,7 +468,7 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
maxLevelVertexBot => grid % maxLevelVertexBot % array
- surfaceThickness => grid % surfaceThickness % array
+ surfacePressure => grid % surfacePressure % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -803,10 +803,7 @@
! mrp 120724 use same rho as top cell right now.
! Later this would be a rho that represents the land ice.
-! pressure(1,iCell) = config_surface_rho0*gravity*surfaceThickness(iCell) &
-! + rho(1,iCell)*gravity*0.5*h(1,iCell)
- pressure(1,iCell) = rho(1,iCell)*gravity*surfaceThickness(iCell) &
- + rho(1,iCell)*gravity*0.5*h(1,iCell)
+ pressure(1,iCell) = surfacePressure(iCell) + rho(1,iCell)*gravity*0.5*h(1,iCell)
do k=2,maxLevelCell(iCell)
pressure(k,iCell) = pressure(k-1,iCell) &
@@ -842,8 +839,7 @@
! and z-coordinates are negative below the surface.
ssh(iCell) = -referenceBottomDepth(maxLevelCell(iCell)) &
- + sum(h(1:maxLevelCell(iCell),iCell)) &
- + surfaceThickness(iCell)
+ + sum(h(1:maxLevelCell(iCell),iCell))
end do
</font>
</pre>