<p><b>mpetersen@lanl.gov</b> 2012-02-08 08:59:05 -0700 (Wed, 08 Feb 2012)</p><p>BRANCH COMMIT: Corrections from Todd's review comments, some improved commenting. Todd recommended subdividing ocn_time_integrator_split with more subroutines, but I am postponing that work for now.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/Registry 2012-02-08 00:31:49 UTC (rev 1480)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/Registry 2012-02-08 15:59:05 UTC (rev 1481)
@@ -207,10 +207,6 @@
var persistent real FBtr ( nEdges Time ) 2 - FBtr state - -
var persistent real GBtrForcing ( nEdges Time ) 2 - GBtrForcing state - -
var persistent real uBcl ( nVertLevels nEdges Time ) 2 - uBcl state - -
-var persistent real circulationBtr ( nVertices Time ) 2 - circulationBtr state - -
-var persistent real divergenceBtr ( nCells Time ) 2 - divergenceBtr state - -
-var persistent real vorticityBtr ( nVertices Time ) 2 - vorticityBtr state - -
-var persistent real u_diffusionBtr ( nEdges Time ) 2 - u_diffusionBtr state - -
% Diagnostic fields: only written to output
var persistent real zMid ( nVertLevels nCells Time ) 2 io zMid state - -
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-08 00:31:49 UTC (rev 1480)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_mpas_core.F 2012-02-08 15:59:05 UTC (rev 1481)
@@ -225,6 +225,8 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(1) % state, mesh)
call mpas_timer_stop("diagnostic solve", initDiagSolveTimer)
+ call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(1) % state, mesh)
+
call ocn_compute_mesh_scaling(mesh)
call mpas_rbf_interp_initialize(mesh)
@@ -561,12 +563,12 @@
! ocn_diagnostic_solve has not yet been called, so compute hEdge
! just for this edge.
- ! I took out k=1 so that hSum always has a nonzero value, to avoid
- ! nans in uBtr.
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
hEdge1 = 0.5*( &
block % state % time_levs(1) % state % h % array(1,cell1) &
+ block % state % time_levs(1) % state % h % array(1,cell2) )
-
uhSum = hEdge1*block % state % time_levs(1) % state % u % array(1,iEdge)
hSum = hEdge1
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-08 00:31:49 UTC (rev 1480)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-08 15:59:05 UTC (rev 1481)
@@ -86,21 +86,17 @@
type (dm_info) :: dminfo
integer :: iCell, i,k,j, iEdge, cell1, cell2, split_explicit_step, split, &
eoe, oldBtrSubcycleTime, newBtrSubcycleTime, uPerpTime, BtrCorIter, &
- n_bcl_iter(config_n_ts_iter), &
-!TDR: vertex1, vertex2, iVertex unused
- vertex1, vertex2, iVertex
+ n_bcl_iter(config_n_ts_iter)
type (block_type), pointer :: block
real (kind=RKIND) :: uhSum, hSum, flux, sshEdge, &
- uPerp, uCorr, temp, temp_h, coef, FBtr_coeff, sshCell1, sshCell2
+ CoriolisTerm, uCorr, temp, temp_h, coef, FBtr_coeff, sshCell1, sshCell2
integer :: num_tracers, ucorr_coef, err
real (kind=RKIND), dimension(:,:), pointer :: &
u, h, h_edge, ke_edge, vertViscTopOfEdge, vertDiffTopOfCell
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
-!TDR: i think that A and C are unused. if they are used, the names should be
-!TDR: expanded in order to find via grep
- real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp
+ real (kind=RKIND), dimension(:), allocatable:: uTemp
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
call mpas_timer_start("se timestep", .false., timer_main)
@@ -114,40 +110,41 @@
block => domain % blocklist
do while (associated(block))
+ ! Initialize * variables that are used to compute baroclinic tendencies below.
do iEdge=1,block % mesh % nEdges
+ do k=1,block % mesh % nVertLevels !maxLevelEdgeTop % array(iEdge)
- ! The baroclinic velocity needs be recomputed at the beginning of a
- ! timestep because the implicit vertical mixing is conducted on the
- ! total u. We keep uBtr from the previous timestep.
- block % state % time_levs(1) % state % uBcl % array(:,iEdge) &
- = block % state % time_levs(1) % state % u % array(:,iEdge) &
- - block % state % time_levs(1) % state % uBtr % array( iEdge)
+ ! The baroclinic velocity needs be recomputed at the beginning of a
+ ! timestep because the implicit vertical mixing is conducted on the
+ ! total u. We keep uBtr from the previous timestep.
+ block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
+ = block % state % time_levs(1) % state % u % array(k,iEdge) &
+ - block % state % time_levs(1) % state % uBtr % array( iEdge)
- block % state % time_levs(2) % state % u % array(:,iEdge) &
- = block % state % time_levs(1) % state % u % array(:,iEdge)
+ block % state % time_levs(2) % state % u % array(k,iEdge) &
+ = block % state % time_levs(1) % state % u % array(k,iEdge)
- block % state % time_levs(2) % state % uBcl % array(:,iEdge) &
- = block % state % time_levs(1) % state % uBcl % array(:,iEdge)
+ block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
+ = block % state % time_levs(1) % state % uBcl % array(k,iEdge)
- enddo ! iEdge
+ block % state % time_levs(2) % state % h_edge % array(k,iEdge) &
+ = block % state % time_levs(1) % state % h_edge % array(k,iEdge)
- ! Initialize * variables that are used to compute baroclinic tendencies below.
+ end do
+ end do
+
block % state % time_levs(2) % state % ssh % array(:) &
= block % state % time_levs(1) % state % ssh % array(:)
- block % state % time_levs(2) % state % h % array(:,:) &
- = block % state % time_levs(1) % state % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % maxLevelCell % array(iCell)
- block % state % time_levs(2) % state % h_edge % array(:,:) &
- = block % state % time_levs(1) % state % h_edge % array(:,:)
+ block % state % time_levs(2) % state % h % array(k,iCell) &
+ = block % state % time_levs(1) % state % h % array(k,iCell)
-!TDR: either the following code or the comment "couple tracers to h" is incorrect below
- do iCell=1,block % mesh % nCells ! couple tracers to h
- ! change to maxLevelCell % array(iCell) ?
- do k=1,block % mesh % nVertLevels
-
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
= block % state % time_levs(1) % state % tracers % array(:,k,iCell)
+
end do
end do
@@ -243,17 +240,9 @@
/block % mesh % dcEdge % array(iEdge) )
enddo
-!TDR: it would be clean and less prone to error if the following action was computed via
-! uhSum = c0
-! hSum = c0
-! do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
-! uhSum = uhSum + block % state % time_levs(2) % state % h_edge % array(k,iEdge) * uTemp(k)
-! hSum = hSum + block % state % time_levs(2) % state % h_edge % array(k,iEdge)
-! enddo
-! block % state % time_levs(1) % state % GBtrForcing % array(iEdge) = split*uhSum/hSum/dt
-!TDR
-
-
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
uhSum = block % state % time_levs(2) % state % h_edge % array(1,iEdge) * uTemp(1)
hSum = block % state % time_levs(2) % state % h_edge % array(1,iEdge)
@@ -365,7 +354,7 @@
do j=1,config_n_btr_subcycles*config_btr_subcycle_loop_factor
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Barotropic subcycle: initial solve for velecity
+ ! Barotropic subcycle: VELOCITY PREDICTOR STEP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (config_btr_gam1_uWt1>1.0e-12) then ! only do this part if it is needed in next SSH solve
uPerpTime = oldBtrSubcycleTime
@@ -378,26 +367,22 @@
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
-!TDR: we say that we are computing -f*uPerp, but then we call it uPerp. this is confusing on several fronts
-!TDR: first, we are computing an acceleration (i.e. tendency to u) but give it a name that implies velocity
-!TDR: second, we all the total of -f*uPerp to be uPerp ...
-!TDR: I would replace all "uPerp" with "CoriolisAcceleration" or something a bit shorter but along those lines
-
- ! Compute -f*uPerp
- uPerp = 0.0
+ ! Compute the barotropic Coriolis term, -f*uPerp
+ CoriolisTerm = 0.0
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
eoe = block % mesh % edgesOnEdge % array(i,iEdge)
- uPerp = uPerp + block % mesh % weightsOnEdge % array(i,iEdge) &
+ CoriolisTerm = CoriolisTerm &
+ + block % mesh % weightsOnEdge % array(i,iEdge) &
* block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
- * block % mesh % fEdge % array(eoe)
+ * block % mesh % fEdge % array(eoe)
end do
- ! uBtrNew = uBtrOld + dt*(-f*uBtroldPerp - g*grad(SSH) + G)
+ ! uBtrNew = uBtrOld + dt/J*(-f*uBtroldPerp - g*grad(SSH) + G)
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
= (block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- + dt / config_n_btr_subcycles * (uPerp - gravity &
+ + dt / config_n_btr_subcycles * (CoriolisTerm - gravity &
* (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
- - block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
+ - block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
/ block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1, iEdge)
end do
@@ -420,7 +405,7 @@
endif ! config_btr_gam1_uWt1>1.0e-12
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Barotropic subcycle: Compute thickness flux and new SSH: PREDICTOR
+ ! Barotropic subcycle: SSH PREDICTOR STEP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
block => domain % blocklist
do while (associated(block))
@@ -482,7 +467,6 @@
call mpas_timer_start("se halo ssh", .false., timer_halo_ssh)
block => domain % blocklist
do while (associated(block))
- ! block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
call mpas_dmpar_exch_halo_field1d_real(domain % dminfo, &
block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
@@ -493,7 +477,7 @@
call mpas_timer_stop("se halo ssh", timer_halo_ssh)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Barotropic subcycle: Final solve for velocity. Iterate for Coriolis term.
+ ! Barotropic subcycle: VELOCITY CORRECTOR STEP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do BtrCorIter=1,config_n_btr_cor_iter
uPerpTime = newBtrSubcycleTime
@@ -504,27 +488,27 @@
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- ! Compute -f*uPerp
- uPerp = 0.0
+ ! Compute the barotropic Coriolis term, -f*uPerp
+ CoriolisTerm = 0.0
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
eoe = block % mesh % edgesOnEdge % array(i,iEdge)
- uPerp = uPerp + block % mesh % weightsOnEdge % array(i,iEdge) &
+ CoriolisTerm = CoriolisTerm + block % mesh % weightsOnEdge % array(i,iEdge) &
* block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
* block % mesh % fEdge % array(eoe)
end do
+ ! In this final solve for velocity, SSH is a linear
+ ! combination of SSHold and SSHnew.
sshCell1 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
- + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
-
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
sshCell2 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
- + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
-!TDR: I can not find where u_diffusionBtr is computed
+ ! uBtrNew = uBtrOld + dt/J*(-f*uBtroldPerp - g*grad(SSH) + G)
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
= (block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- + dt/config_n_btr_subcycles *(uPerp - gravity *(sshCell2 - sshCell1) /block % mesh % dcEdge % array(iEdge) &
- + block % state % time_levs(1) % state % GBtrForcing % array(iEdge) &
- + block % state % time_levs(1) % state % u_diffusionBtr % array(iEdge))) * block % mesh % edgeMask % array(1,iEdge)
+ + dt/config_n_btr_subcycles *(CoriolisTerm - gravity *(sshCell2 - sshCell1) /block % mesh % dcEdge % array(iEdge) &
+ + block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1,iEdge)
end do
block => block % next
@@ -544,7 +528,7 @@
end do !do BtrCorIter=1,config_n_btr_cor_iter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Barotropic subcycle: Compute thickness flux and new SSH: CORRECTOR
+ ! Barotropic subcycle: SSH CORRECTOR STEP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (config_btr_solve_SSH2) then
@@ -552,27 +536,32 @@
do while (associated(block))
block % tend % ssh % array(:) = 0.0
- ! config_btr_gam3_uWt2 sets the forward weighting of velocity in the SSH computation
+ ! config_btr_gam3_uWt2 sets the forward weighting of velocity in the SSH computation
! config_btr_gam3_uWt2= 1 flux = uBtrNew*H
! config_btr_gam3_uWt2=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
do iEdge=1,block % mesh % nEdges
- cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
- cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
- + block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2))
- hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
+ ! SSH is a linear combination of SSHold and SSHnew.
+ sshCell1 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
+ sshCell2 = (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+
+ sshEdge = 0.5 * (sshCell1 + sshCell2)
+ hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
- flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- + config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
- * (sshEdge + hSum)
+ flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * (sshEdge + hSum)
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
+ block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
+ block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
- block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) + flux
+ block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) + flux
end do
! SSHnew = SSHold + dt/J*(-div(Flux))
@@ -684,6 +673,9 @@
uTemp(:) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ block % state % time_levs(2) % state % uBcl % array(:,iEdge)
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
uhSum = block % state % time_levs(2) % state % h_edge % array(1,iEdge) * uTemp(1)
hSum = block % state % time_levs(2) % state % h_edge % array(1,iEdge)
@@ -724,16 +716,12 @@
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!TDR: it seems almost trivial to hold off on doing T, S and rho updates until the
-!TDR: dycore time step is complete. we might want to take this opportunity to clean-up
-!TDR: Stage3 in order to faciliate the testing of not doing tracer updates after this code is committed to trunk.
-!TDR: at this point, I am suggesting just pushing some of this code into subroutines.
-!TDR: see comments farther down
+ !TDR: it seems almost trivial to hold off on doing T, S and rho updates until the
+ !TDR: dycore time step is complete. we might want to take this opportunity to clean-up
+ !TDR: Stage3 in order to faciliate the testing of not doing tracer updates after this code is committed to trunk.
+ !TDR: at this point, I am suggesting just pushing some of this code into subroutines.
+ !TDR: see comments farther down
-!TDR: this is the first occurrence of computing wTop that I can find.
-!TDR: the call to ocn_wtop is not found at initialization and wTop is not written to restart
-!TDR: so upon restart, won't uBaroclinic be computed with wTop=0 on the first time step?
-
block => domain % blocklist
do while (associated(block))
call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
@@ -769,8 +757,8 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (split_explicit_step < config_n_ts_iter) then
-!TDR: should we move this code into a subroutine called "compute_intermediate_value_at_midtime"
-!TDR: this could be within a contains statement in this routine
+ !TDR: should we move this code into a subroutine called "compute_intermediate_value_at_midtime"
+ !TDR: this could be within a contains statement in this routine
! Only need T & S for earlier iterations,
! then all the tracers needed the last time through.
@@ -821,8 +809,8 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
elseif (split_explicit_step == config_n_ts_iter) then
-!TDR: should we move this code into a subroutine called "compute_final_values_at_nplus1"?
-!TDR: this could be within a contains statement in this routine
+ !TDR: should we move this code into a subroutine called "compute_final_values_at_nplus1"?
+ !TDR: this could be within a contains statement in this routine
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
@@ -968,6 +956,9 @@
do iEdge=1,nEdges
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
uhSum = h_edge(1,iEdge) * tend_u(1,iEdge)
hSum = h_edge(1,iEdge)
@@ -1025,6 +1016,9 @@
do iEdge=1,nEdges
+ ! hSum is initialized outside the loop because on land boundaries
+ ! maxLevelEdgeTop=0, but I want to initialize hSum with a
+ ! nonzero value to avoid a NaN.
uhSum = h_edge(1,iEdge) * u(1,iEdge)
hSum = h_edge(1,iEdge)
</font>
</pre>