<p><b>mpetersen@lanl.gov</b> 2012-07-26 08:34:46 -0600 (Thu, 26 Jul 2012)</p><p>BRANCH COMMIT: I initialized the surface thickness (pressure) in the mpas_ocn_mpas_core.F module. This is just for initial testing.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/surface_pressure_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/surface_pressure_mrp/src/core_ocean/Registry 2012-07-26 13:12:23 UTC (rev 2057)
+++ branches/ocean_projects/surface_pressure_mrp/src/core_ocean/Registry 2012-07-26 14:34:46 UTC (rev 2058)
@@ -29,6 +29,8 @@
namelist character grid config_vert_grid_type isopycnal
namelist character grid config_pressure_type pressure
namelist real grid config_rho0 1028
+namelist real grid config_surface_rho0 1028
+namelist real grid config_max_surface_thickness 0.1
namelist logical grid config_enforce_zstar_at_restart false
namelist integer split_explicit_ts config_n_ts_iter 2
namelist integer split_explicit_ts config_n_bcl_iter_beg 2
@@ -208,9 +210,9 @@
var persistent real salinityRestore ( nCells ) 0 iro salinityRestore mesh - -
% mrp trying to figure out why these do not appear
-var persistent real windStressMonthly ( nMonths nEdges ) 0 iro windStressMonthly mesh - -
-var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 iro temperatureRestoreMonthly mesh - -
-var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 iro salinityRestoreMonthly mesh - -
+var persistent real windStressMonthly ( nMonths nEdges ) 0 - windStressMonthly mesh - -
+var persistent real temperatureRestoreMonthly ( nMonths nCells ) 0 - temperatureRestoreMonthly mesh - -
+var persistent real salinityRestoreMonthly ( nMonths nCells ) 0 - salinityRestoreMonthly mesh - -
% Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 ir u state - -
@@ -273,6 +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 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 13:12:23 UTC (rev 2057)
+++ branches/ocean_projects/surface_pressure_mrp/src/core_ocean/mpas_ocn_mpas_core.F 2012-07-26 14:34:46 UTC (rev 2058)
@@ -248,7 +248,80 @@
integer, intent(out) :: err
integer :: i, iEdge, iCell, k
integer :: err1
-
+
+! 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
+
+ ! make even layers
+ totalH = 1600.0
+ block % state % time_levs(1) % state % h % array = totalH/block % mesh % nVertLevels
+ block % mesh % hZLevel % array(:) = totalH/block % mesh % nVertLevels
+
+ ! turn off the wind.
+ !block % mesh % u_src % array = 0.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))
+
+ x_mid = 1600.0e3 ! 1600km
+ x_efold = 800.0e3 ! 800km
+ do iCell=1,block % mesh % nCells
+ xCell = block % mesh % xCell % array(iCell)
+
+! 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
+ else
+ ! Gaussian
+ !block % mesh % surfaceThickness % array(iCell) &
+ ! = config_max_surface_thickness*totalH &
+ ! * exp( -((xCell-x_mid)/x_efold)**2 )
+
+ ! top hat
+ block % mesh % surfaceThickness % array(iCell) &
+ = config_max_surface_thickness*totalH
+
+ block % state % time_levs(1) % state % h % array(:,iCell) &
+ = (1.0 -block % mesh % surfaceThickness % array(iCell)/totalH) &
+ *block % state % time_levs(1) % state % h % array(:,iCell)
+ endif
+
+ ! 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
+ ! and z-coordinates are negative below the surface.
+ k = block % mesh % nVertLevels
+ block % state % time_levs(1) % state % zMid % array(k:block % mesh % nVertLevels,iCell) = &
+ -totalH + 0.5*block % state % time_levs(1) % state % h % array(k,iCell)
+
+ do k=block % mesh % nVertLevels-1, 1, -1
+ block % state % time_levs(1) % state % zMid % array(k,iCell) = block % state % time_levs(1) % state % zMid % array(k+1,iCell) &
+ + 0.5*( block % state % time_levs(1) % state % h % array(k+1,iCell) &
+ + block % state % time_levs(1) % state % h % array(k ,iCell))
+ end do
+
+ ! constant T
+ !block % state % time_levs(1) % state % tracers % array(1,:,iCell) = 10.0
+ ! set temperature as a linear function of depth, 10C at top, 0C at bottom
+ do k=1,block % mesh % nVertLevels
+ block % state % time_levs(1) % state % tracers % array(1,k,iCell) = &
+ 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)
+
+ enddo
+
+ ! make T & S constant
+ block % state % time_levs(1) % state % tracers % array(2,:,:) = 35.0
+ block % state % time_levs(1) % state % tracers % array(3,:,:) = 1.0
+
+! mrp 120724 initialize surfaceThickness end
+
+
call ocn_initialize_advection_rk(mesh, err)
call mpas_ocn_tracer_advection_coefficients(mesh, err1)
err = ior(err, err1)
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 13:12:23 UTC (rev 2057)
+++ branches/ocean_projects/surface_pressure_mrp/src/core_ocean/mpas_ocn_tendency.F 2012-07-26 14:34:46 UTC (rev 2058)
@@ -412,7 +412,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, &
- referenceBottomDepth, ssh
+ referenceBottomDepth, ssh, surfaceThickness
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -468,6 +468,7 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelEdgeBot => grid % maxLevelEdgeBot % array
maxLevelVertexBot => grid % maxLevelVertexBot % array
+ surfaceThickness => grid % surfaceThickness % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -797,11 +798,16 @@
else
do iCell=1,nCells
- ! pressure for generalized coordinates
- ! assume atmospheric pressure at the surface is zero for now.
- pressure(1,iCell) = rho(1,iCell)*gravity &
- * 0.5*h(1,iCell)
+ ! Pressure for generalized coordinates
+ ! Start with surface pressure.
+ ! 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)
+
do k=2,maxLevelCell(iCell)
pressure(k,iCell) = pressure(k-1,iCell) &
+ 0.5*gravity*( rho(k-1,iCell)*h(k-1,iCell) &
@@ -836,7 +842,8 @@
! and z-coordinates are negative below the surface.
ssh(iCell) = -referenceBottomDepth(maxLevelCell(iCell)) &
- + sum(h(1:maxLevelCell(iCell),iCell))
+ + sum(h(1:maxLevelCell(iCell),iCell)) &
+ + surfaceThickness(iCell)
end do
</font>
</pre>