<p><b>mpetersen@lanl.gov</b> 2010-02-08 10:44:38 -0700 (Mon, 08 Feb 2010)</p><p>Added the following:<br>
1. Computation of zmid, zbot, pmid, pbot, and Montgomery potential in <br>
subroutine compute_solve_diagnostics in module_time_integration.F<br>
2. Use of Montgomery potential in<br>
subroutine compute_tend in module_time_integration.F<br>
3. Compute statistics, for one processor only, in <br>
subroutine write_stats in module_sw_solver.F<br>
4. Initialization of h, rho, alpha (1/rho), and u_src<br>
module_test_cases.F<br>
5. Addition of new variables in Registry.<br>
<br>
Note that to run a multilayer case, use<br>
EXPAND_LEVELS = -DEXPAND_LEVELS=4<br>
in Makefile.<br>
<br>
Also, statistics are placed in a stats subdirectory that must be<br>
created before the run.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean/Registry/Registry
===================================================================
--- branches/ocean/Registry/Registry 2010-02-05 23:51:55 UTC (rev 114)
+++ branches/ocean/Registry/Registry 2010-02-08 17:44:38 UTC (rev 115)
@@ -6,6 +6,9 @@
namelist real sw_model config_dt 172.8
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
+# mrp 100120:
+namelist integer sw_model config_stats_interval 100
+# mrp 100120 end
namelist real sw_model config_visc 0.0
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
@@ -76,6 +79,10 @@
var real fEdge ( nEdges ) iro fEdge - -
var real fVertex ( nVertices ) iro fVertex - -
var real h_s ( nCells ) iro h_s - -
+# mrp 100112:
+var real rho ( nVertLevels ) iro rho - -
+var real alpha ( nVertLevels ) iro alpha - -
+# mrp 100112 end
# Arrays required for reconstruction of velocity field
var real matrix_reconstruct ( maxEdges maxEdges nCells ) - matrix_reconstruct - -
@@ -104,7 +111,16 @@
var real uReconstructX ( nVertLevels nCells Time ) o uReconstructX - -
var real uReconstructY ( nVertLevels nCells Time ) o uReconstructY - -
var real uReconstructZ ( nVertLevels nCells Time ) o uReconstructZ - -
+# mrp 100112:
+var real zmid ( nVertLevels nCells Time ) o zmid - -
+var real zbot ( nVertLevels nCells Time ) o zbot - -
+var real zSurface ( nCells Time ) o zSurface - -
+var real pmid ( nVertLevels nCells Time ) o pmid - -
+var real pbot ( nVertLevels nCells Time ) o pbot - -
+var real MontPot ( nVertLevels nCells Time ) o MontPot - -
+# mrp 100112 end
+
# Other diagnostic variables: neither read nor written to any files
var real vh ( nVertLevels nEdges Time ) - vh - -
var real circulation ( nVertLevels nVertices Time ) - circulation - -
Modified: branches/ocean/namelist.input
===================================================================
--- branches/ocean/namelist.input 2010-02-05 23:51:55 UTC (rev 114)
+++ branches/ocean/namelist.input 2010-02-08 17:44:38 UTC (rev 115)
@@ -1,12 +1,19 @@
&sw_model
- config_test_case = 0
+ config_test_case = 5
config_time_integration = 'RK4'
config_dt = 300.0
- config_ntimesteps = 288000
- config_output_interval = 288
+ config_ntimesteps = 3000
+ config_output_interval = 300
+ config_stats_interval = 10
config_visc = 4.0
/
+&io
+ config_input_name = grid.nc
+ config_output_name = output.nc
+ config_restart_name = restart.nc
+/
+
&restart
config_restart_interval = 864000
config_do_restart = .false.
Modified: branches/ocean/src/module_sw_solver.F
===================================================================
--- branches/ocean/src/module_sw_solver.F 2010-02-05 23:51:55 UTC (rev 114)
+++ branches/ocean/src/module_sw_solver.F 2010-02-08 17:44:38 UTC (rev 115)
@@ -51,6 +51,9 @@
! Before integrating, write out the initial state
output_frame = 1
restart_frame = 1
+ ! mrp 100120:
+ call write_stats(domain, 0, dt)
+ ! mrp 100120 end
call output_state_init(output_obj, domain, "OUTPUT")
call write_output_frame(output_obj, domain)
@@ -58,7 +61,7 @@
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
do itimestep = 1,ntimesteps
- write(0,*) 'Doing timestep ', itimestep
+ !write(0,*) 'Doing timestep ', itimestep
call timer_start("time integration")
call timestep(domain, dt)
call timer_stop("time integration")
@@ -66,6 +69,12 @@
! Move time level 2 fields back into time level 1 for next time step
call shift_time_levels(domain)
+ ! mrp 100120:
+ if (mod(itimestep, config_stats_interval) == 0) then
+ call write_stats(domain, itimestep, dt)
+ end if
+ ! mrp 100120 end
+
if (mod(itimestep, config_output_interval) == 0) then
call write_output_frame(output_obj, domain)
end if
@@ -110,7 +119,130 @@
end subroutine write_output_frame
+! mrp 100120:
+ subroutine write_stats(domain, itimestep, dt)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute statistics
+ !
+ ! Input/Output: domain - contains model state; diagnostic field are computed
+ ! before returning
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! mrp 100120
+! This is a first cut at global variable output. It has many restrictions
+! and items to correct:
+! 1. This works on a single processors only. We need global min,max,sums
+! across processors before the write (c1,'(i6.5)') 100000 line.
+! 2. There should be a variable for the number of blocks owned by this processor,
+! but I could not find it. That should replace my funny n variable.
+! 3. This nblocks variable should replace the first dim of _stat variables
+! 4. The avg is not cell-area weighted.
+! 5. Boundary cells are included in the sum and avg, and should not.
+! 6. The assignments like
+! case(1); var => block % time_levs(1) % state % h % array
+! are long and prone to errors. This should be automated with a registry
+! entry for statistics writes, which is a different frequency from binary output.
+! 7. The nvars is hardcoded, but this should really come from the registry.
+! 8. I only did 3D variables. 2D variables will need a var2D(:) rather than
+! var(:,:) for the 3D - here var(nCells,nLevels), etc.
+
+ implicit none
+
+ type (domain_type), intent(in) :: domain
+ integer, intent(in) :: itimestep
+ real (kind=RKIND), intent(in) :: dt
+
+ integer :: i, j, k, n, nvars, num
+ type (block_type), pointer :: block
+ character :: c1*6
+ real (kind=RKIND), dimension(20,20) :: &
+ min_stat, max_stat, sum_stat, avg_stat, colmin_stat, colmax_stat
+ real (kind=RKIND), dimension(20) :: &
+ min_stat_all, max_stat_all, sum_stat_all, avg_stat_all, &
+ colmin_stat_all, colmax_stat_all
+ real (kind=RKIND), dimension(:,:), pointer :: var
+ real (kind=RKIND) :: current_time
+
+ nvars=17
+
+ n=0
+ block => domain % blocklist
+
+ do while (associated(block))
+ n=n+1
+ do j=1,nvars
+ num=block % mesh % nCells
+ select case(j)
+ case(1); var => block % time_levs(1) % state % h % array
+ case(2); var => block % time_levs(1) % state % u % array
+ case(3); var => block % time_levs(1) % state % v % array
+ case(4); var => block % time_levs(1) % state % h_edge % array
+ num=block % mesh % nEdges
+ case(5); var => block % time_levs(1) % state % circulation % array
+ case(6); var => block % time_levs(1) % state % vorticity % array
+ case(7); var => block % time_levs(1) % state % ke % array
+ case(8); var => block % time_levs(1) % state % pv_edge % array
+ num=block % mesh % nEdges
+ case(9); var => block % time_levs(1) % state % pv_vertex % array
+ num=block % mesh % nVertices
+ case(10); var => block % time_levs(1) % state % pv_cell % array
+ case(11); var => block % time_levs(1) % state % gradPVn % array
+ num=block % mesh % nEdges
+ case(12); var => block % time_levs(1) % state % gradPVt % array
+ num=block % mesh % nEdges
+ case(13); var => block % time_levs(1) % state % zmid % array
+ case(14); var => block % time_levs(1) % state % zbot % array
+ case(15); var => block % time_levs(1) % state % pmid % array
+ case(16); var => block % time_levs(1) % state % pbot % array
+ case(17); var => block % time_levs(1) % state % MontPot % array
+ end select
+ colmin_stat(n,j) = minval(sum(var,1))
+ colmax_stat(n,j) = maxval(sum(var,1))
+ min_stat(n,j) = minval(var)
+ max_stat(n,j) = maxval(var)
+ sum_stat(n,j) = sum(var)
+ avg_stat(n,j) = sum(var)/num/block % mesh % nVertLevels
+ enddo
+ block => block % next
+ end do
+
+ do j=1,nvars
+ min_stat_all(j) = minval(min_stat(1:n,j))
+ max_stat_all(j) = maxval(max_stat(1:n,j))
+ sum_stat_all(j) = sum(sum_stat(1:n,j))
+ avg_stat_all(j) = sum(avg_stat(1:n,j))/n
+ colmin_stat_all(j) = minval(colmin_stat(1:n,j))
+ colmax_stat_all(j) = maxval(colmax_stat(1:n,j))
+ enddo
+
+ write (c1,'(i6.5)') 100000
+
+ open(51,file='stats/min.'//c1(2:6)//'.txt',ACCESS='append')
+ open(52,file='stats/max.'//c1(2:6)//'.txt',ACCESS='append')
+ open(53,file='stats/sum.'//c1(2:6)//'.txt',ACCESS='append')
+ open(54,file='stats/avg.'//c1(2:6)//'.txt',ACCESS='append')
+ open(55,file='stats/time.'//c1(2:6)//'.txt',ACCESS='append')
+ open(56,file='stats/colmin.'//c1(2:6)//'.txt',ACCESS='append')
+ open(57,file='stats/colmax.'//c1(2:6)//'.txt',ACCESS='append')
+ write (51,'(100es24.16)') min_stat_all(1:nvars)
+ write (52,'(100es24.16)') max_stat_all(1:nvars)
+ write (53,'(100es24.16)') sum_stat_all(1:nvars)
+ write (54,'(100es24.16)') avg_stat_all(1:nvars)
+ write (55,'(i5,100es24.16)') itimestep, &
+ domain % blocklist % time_levs(1) % state % xtime % scalar, dt
+ write (56,'(100es24.16)') colmin_stat_all(1:nvars)
+ write (57,'(100es24.16)') colmax_stat_all(1:nvars)
+ close (51)
+ close (52)
+ close (53)
+ close (54)
+ close (55)
+ close (56)
+ close (57)
+
+ end subroutine write_stats
+! mrp 100120 end
+
subroutine compute_output_diagnostics(state, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields for a domain
Modified: branches/ocean/src/module_test_cases.F
===================================================================
--- branches/ocean/src/module_test_cases.F 2010-02-05 23:51:55 UTC (rev 114)
+++ branches/ocean/src/module_test_cases.F 2010-02-08 17:44:38 UTC (rev 115)
@@ -43,7 +43,8 @@
else if (config_test_case == 2) then
write(0,*) 'Setting up shallow water test case 2'
- write(0,*) ' -- Setup shallow water test case 2: Global Steady State Nonlinear Zonal Geostrophic Flow'
+ write(0,*) ' -- Setup shallow water test case 2: '// &
+ 'Global Steady State Nonlinear Zonal Geostrophic Flow'
block_ptr => domain % blocklist
do while (associated(block_ptr))
@@ -57,7 +58,8 @@
else if (config_test_case == 5) then
write(0,*) 'Setting up shallow water test case 5'
- write(0,*) ' -- Setup shallow water test case 5: Zonal Flow over an Isolated Mountain'
+ write(0,*) ' -- Setup shallow water test case 5:'// &
+ ' Zonal Flow over an Isolated Mountain'
block_ptr => domain % blocklist
do while (associated(block_ptr))
@@ -84,7 +86,8 @@
end do
else
- write(0,*) 'Only test case 1, 2, 5, and 6 are currently supported.'
+ write(0,*) &
+ 'Only test case 1, 2, 5, and 6 are currently supported.'
stop
end if
@@ -111,7 +114,7 @@
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
real (kind=RKIND), parameter :: alpha = pii/4.0
- integer :: iCell, iEdge, iVtx
+ integer :: iCell, iEdge, iVtx, iLevel !mrp on last one
real (kind=RKIND) :: r, u, v
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -165,10 +168,32 @@
state % h % array(1,iCell) = 0.0
end if
#ifdef EXPAND_LEVELS
- state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
+ ! mrp 100125:
+ !state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
+ ! make the total thickness equal to the single-layer thickness:
+ state % h % array(1:EXPAND_LEVELS, iCell) = state % h % array(1,iCell) &
+ / grid % nVertLevels
+ ! mrp 100125 end
+#else
+ ! mrp 100125:
+ ! make the total thickness equal to the single-layer thickness:
+ !state % h % array(:, iCell) = state % h % array(1,iCell) / 4.0
+ ! mrp 100125 end
#endif
end do
+ ! mrp 100121:
+ ! define the density
+ do iLevel = 1,grid % nVertLevels
+ grid % rho % array(iLevel) = (iLevel-1)*10
+ enddo
+ grid % rho % array = grid % rho % array &
+ - sum(grid % rho % array)/grid % nVertLevels + 1000.0
+print *, 'rho', grid % rho % array, sum(grid % rho % array)/grid % nVertLevels
+ grid % alpha % array = 1/grid % rho % array
+ grid % u_src % array = 0.0
+ ! mrp 100121 end
+
end subroutine sw_test_case_1
@@ -192,6 +217,7 @@
real (kind=RKIND), parameter :: alpha = 0.0
integer :: iCell, iEdge, iVtx
+ integer :: iLevel !mrp
real (kind=RKIND) :: u, v
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -263,10 +289,32 @@
) / &
gravity
#ifdef EXPAND_LEVELS
- state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
+ ! mrp 100125:
+ !state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
+ ! make the total thickness equal to the single-layer thickness:
+ state % h % array(1:EXPAND_LEVELS, iCell) = state % h % array(1,iCell) &
+ / grid % nVertLevels
+ ! mrp 100125 end
+#else
+ ! mrp 100125:
+ ! make the total thickness equal to the single-layer thickness:
+ !state % h % array(:, iCell) = state % h % array(1,iCell) / 4.0
+ ! mrp 100125 end
#endif
end do
+ ! mrp 100121:
+ ! define the density
+ do iLevel = 1,grid % nVertLevels
+ grid % rho % array(iLevel) = (iLevel-1)*10
+ enddo
+ grid % rho % array = grid % rho % array &
+ - sum(grid % rho % array)/grid % nVertLevels + 1000.0
+print *, 'rho', grid % rho % array, sum(grid % rho % array)/grid % nVertLevels
+ grid % alpha % array = 1/grid % rho % array
+ grid % u_src % array = 0.0
+ ! mrp 100121 end
+
end subroutine sw_test_case_2
@@ -286,13 +334,14 @@
real (kind=RKIND), parameter :: u0 = 20.
real (kind=RKIND), parameter :: gh0 = 5960.0*gravity
- real (kind=RKIND), parameter :: hs0 = 2000.
+! real (kind=RKIND), parameter :: hs0 = 2000. original
+ real (kind=RKIND), parameter :: hs0 = 250. !mrp 100204
real (kind=RKIND), parameter :: theta_c = pii/6.0
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
real (kind=RKIND), parameter :: rr = pii/9.0
real (kind=RKIND), parameter :: alpha = 0.0
- integer :: iCell, iEdge, iVtx
+ integer :: iCell, iEdge, iVtx, iLevel
real (kind=RKIND) :: r, u, v
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -360,6 +409,8 @@
r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
grid % h_s % array(iCell) = hs0 * (1.0 - r/rr)
end do
+! output about mountain
+print *, 'h_s',minval(grid % h_s % array),sum(grid % h_s % array)/grid % nCells, maxval(grid % h_s % array)
!
! Initialize tracer fields
@@ -388,10 +439,32 @@
gravity
state % h % array(1,iCell) = state % h % array(1,iCell) - grid % h_s % array(iCell)
#ifdef EXPAND_LEVELS
- state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
+ ! mrp 100125:
+ !state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
+ ! make the total thickness equal to the single-layer thickness:
+ state % h % array(1:EXPAND_LEVELS, iCell) = state % h % array(1,iCell) &
+ / grid % nVertLevels
+ ! mrp 100125 end
+#else
+ ! mrp 100125:
+ ! make the total thickness equal to the single-layer thickness:
+ !state % h % array(:, iCell) = state % h % array(1,iCell) / 4.0
+ ! mrp 100125 end
#endif
end do
+ ! mrp 100121:
+ ! define the density
+ do iLevel = 1,grid % nVertLevels
+ grid % rho % array(iLevel) = 0.0 + (iLevel-1)*10
+ enddo
+ grid % rho % array = grid % rho % array &
+ - sum(grid % rho % array)/grid % nVertLevels + 1000.0
+print *, 'rho', grid % rho % array, sum(grid % rho % array)/grid % nVertLevels
+ grid % alpha % array = 1/grid % rho % array
+ grid % u_src % array = 0.0
+ ! mrp 100121 end
+
end subroutine sw_test_case_5
@@ -414,7 +487,7 @@
real (kind=RKIND), parameter :: K = 7.848e-6
real (kind=RKIND), parameter :: R = 4.0
- integer :: iCell, iEdge, iVtx
+ integer :: iCell, iEdge, iVtx, iLevel
real (kind=RKIND) :: u, v
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -466,10 +539,32 @@
a*a*CC(grid%latCell%array(iCell)) * cos(2.0*R*grid%lonCell%array(iCell)) &
) / gravity
#ifdef EXPAND_LEVELS
- state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
+ ! mrp 100125:
+ !state % h % array(2:EXPAND_LEVELS, iCell) = state % h % array(1,iCell)
+ ! make the total thickness equal to the single-layer thickness:
+ state % h % array(1:EXPAND_LEVELS, iCell) = state % h % array(1,iCell) &
+ / grid % nVertLevels
+ ! mrp 100125 end
+#else
+ ! mrp 100125:
+ ! make the total thickness equal to the single-layer thickness:
+ !state % h % array(:, iCell) = state % h % array(1,iCell) / 4.0
+ ! mrp 100125 end
#endif
end do
+ ! mrp 100121:
+ ! define the density
+ do iLevel = 1,grid % nVertLevels
+ grid % rho % array(iLevel) = (iLevel-1)*10
+ enddo
+ grid % rho % array = grid % rho % array &
+ - sum(grid % rho % array)/grid % nVertLevels + 1000.0
+print *, 'rho', grid % rho % array, sum(grid % rho % array)/grid % nVertLevels
+ grid % alpha % array = 1/grid % rho % array
+ grid % u_src % array = 0.0
+ ! mrp 100121 end
+
end subroutine sw_test_case_6
Modified: branches/ocean/src/module_time_integration.F
===================================================================
--- branches/ocean/src/module_time_integration.F 2010-02-05 23:51:55 UTC (rev 114)
+++ branches/ocean/src/module_time_integration.F 2010-02-08 17:44:38 UTC (rev 115)
@@ -256,6 +256,9 @@
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: u_diffusion, visc
+ !mrp 100112:
+ real (kind=RKIND), dimension(:,:), pointer :: MontPot
+ !mrp 100112 end
!ocean
real (kind=RKIND), dimension(:,:), pointer :: u_src
@@ -274,6 +277,9 @@
ke => s % ke % array
pv_edge => s % pv_edge % array
vh => s % vh % array
+ !mrp 100112:
+ MontPot => s % MontPot % array
+ !mrp 100112 end
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -359,12 +365,21 @@
workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
end do
+ ! mrp 100112, this is the original shallow water formulation with grad H:
+ !tend_u(k,iEdge) = &
+ ! q &
+ ! + u_diffusion &
+ ! - ( ke(k,cell2) - ke(k,cell1) &
+ ! gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
+ ! ) / dcEdge(iEdge)
+ ! mrp 100112, changed to grad Montgomery potential:
tend_u(k,iEdge) = &
- q &
- + u_diffusion &
- - ( ke(k,cell2) - ke(k,cell1) + &
- gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ) / dcEdge(iEdge)
+ q &
+ + u_diffusion &
+ - ( ke(k,cell2) - ke(k,cell1) &
+ + MontPot(k,cell2) - MontPot(k,cell1) &
+ ) / dcEdge(iEdge)
+ ! mrp 100112 end
!ocean
tend_u(k,iEdge) = tend_u(k,iEdge) + u_src(k,iEdge)/rho_ref/h_edge(k,iEdge)
@@ -471,6 +486,14 @@
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence
+ !mrp 100112:
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ zmid, zbot, pmid, pbot, MontPot
+ real (kind=RKIND), dimension(:), pointer :: zSurface, rho, alpha
+ real (kind=RKIND) :: delta_p
+ character :: c1*6
+ !mrp 100112 end
+
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, uBC
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -492,6 +515,14 @@
pv_cell => s % pv_cell % array
gradPVn => s % gradPVn % array
gradPVt => s % gradPVt % array
+ !mrp 100112:
+ zmid => s % zmid % array
+ zbot => s % zbot % array
+ pmid => s % pmid % array
+ pbot => s % pbot % array
+ MontPot => s % MontPot % array
+ zSurface => s % zSurface % array
+ !mrp 100112 end
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -510,6 +541,10 @@
h_s => grid % h_s % array
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
+ !mrp 100112:
+ rho => grid % rho % array
+ alpha => grid % alpha % array
+ !mrp 100112 end
nCells = grid % nCells
nEdges = grid % nEdges
@@ -748,7 +783,42 @@
enddo
endif
+ !mrp 100112:
+ !
+ ! Compute mid- and bottom-depth of each layer, from bottom up.
+ ! Then compute mid- and bottom-pressure of each layer, and
+ ! Montgomery Potential, from top down
+ !
+ do iCell=1,nCells
+ ! h_s is ocean topography: elevation above lowest point,
+ ! and is positive. z coordinates are pos upward, and z=0
+ ! is at lowest ocean point.
+ zbot(nVertLevels,iCell) = h_s(iCell)
+ zmid(nVertLevels,iCell) = zbot(nVertLevels,iCell) + 0.5*h(nVertLevels,iCell)
+ do k=nVertLevels-1,1,-1
+ zbot(k,iCell) = zbot(k+1,iCell) + h(k+1,iCell)
+ zmid(k,iCell) = zbot(k,iCell) + 0.5*h(k,iCell)
+ end do
+ ! rather than using zbot(0,iCell), I am adding another 2D variable.
+ zSurface(iCell) = zbot(1,iCell) + h(1,iCell)
+
+ ! assume pressure at the surface is zero for now.
+ pmid(1,iCell) = 0.5*rho(1)*gravity* h(1,iCell) ! + psurf(iCell)
+ pbot(1,iCell) = rho(1)*gravity* h(1,iCell) ! + psurf(iCell)
+ MontPot(1,iCell) = gravity * zSurface(iCell)
+ do k=2,nVertLevels
+ delta_p = rho(k)*gravity* h(k,iCell)
+ pmid(k,iCell) = pbot(k-1,iCell) + 0.5*delta_p
+ pbot(k,iCell) = pbot(k-1,iCell) + delta_p
+
+ ! from delta M = p delta alpha
+ MontPot(k,iCell) = MontPot(k-1,iCell) &
+ + pbot(k-1,iCell)*(alpha(k) - alpha(k-1))
+ end do
+ end do
+ !mrp 100112 end
+
end subroutine compute_solve_diagnostics
</font>
</pre>