<p><b>mpetersen@lanl.gov</b> 2012-01-31 14:08:15 -0700 (Tue, 31 Jan 2012)</p><p>Fixing some bugs on the ocean_core trunk:<br>
<br>
1. A variable was not initialized, <br>
areaCell(nCells+1) = -1.0e34<br>
<br>
2. The first barotropic splitting would get hSum=0 on land edges,<br>
so changed so hSum is nonzero.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-01-31 20:00:57 UTC (rev 1441)
+++ trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-01-31 21:08:15 UTC (rev 1442)
@@ -243,6 +243,8 @@
! The reconstructed velocity on land will have values not exactly
! -1e34 due to the interpolation of reconstruction.
+ block % mesh % areaCell % array(block % mesh % nCells+1) = -1.0e34
+
do iEdge=1,block % mesh % nEdges
! mrp 101115 note: in order to include flux boundary conditions, the following
! line will need to change. Right now, set boundary edges between land and
@@ -558,10 +560,19 @@
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
! uBtr = sum(u)/sum(h) on each column
- uhSum = 0.0
- hSum = 0.0
+ ! ocn_diagnostic_solve has not yet been called, so compute hEdge
+ ! just for this edge.
- do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ ! I took out k=1 so that hSum always has a nonzero value, to avoid
+ ! nans in uBtr.
+ 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
+
+ do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
! ocn_diagnostic_solve has not yet been called, so compute hEdge
! just for this edge.
hEdge1 = 0.5*( &
@@ -571,6 +582,7 @@
uhSum = uhSum &
+ hEdge1*block % state % time_levs(1) % state % u % array(k,iEdge)
hSum = hSum + hEdge1
+
enddo
block % state % time_levs(1) % state % uBtr % array(iEdge) = uhSum/hsum
Modified: trunk/mpas/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-01-31 20:00:57 UTC (rev 1441)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tendency.F 2012-01-31 21:08:15 UTC (rev 1442)
@@ -512,8 +512,6 @@
! initialize h_edge to avoid divide by zero and NaN problems.
h_edge = -1.0e34
- h_edge = -1.0e34
-
do iEdge=1,nEdges*hadv2nd
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -593,7 +591,6 @@
! set the velocity and height at dummy address
! used -1e34 so error clearly occurs if these values are used.
!
-!mrp 110516 change to zero, change back later:
u(:,nEdges+1) = -1e34
h(:,nCells+1) = -1e34
tracers(s % index_temperature,:,nCells+1) = -1e34
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-01-31 20:00:57 UTC (rev 1441)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-01-31 21:08:15 UTC (rev 1442)
@@ -859,9 +859,14 @@
! This overwrites u^{tr}, the tracer transport velocity, which was in uNew.
! The following must occur after call OcnTendScalar
do iEdge=1,block % mesh % nEdges
- block % state % time_levs(2) % state % u % array(:,iEdge) &
- = block % state % time_levs(2) % state % uBtr % array(iEdge) &
- + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ block % state % time_levs(2) % state % u % array(k,iEdge) &
+ = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(k,iEdge)
+ end do
+ do k=block % mesh % maxLevelEdgeTop % array(iEdge)+1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % u % array(k,iEdge) = 0.0
+ end do
end do ! iEdge
elseif (split_explicit_step == config_n_ts_iter) then
@@ -1012,7 +1017,6 @@
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
-
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
@@ -1024,6 +1028,7 @@
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
+
block => block % next
end do
call mpas_timer_stop("se timestep", timer_main)
Modified: trunk/mpas/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-01-31 20:00:57 UTC (rev 1441)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-01-31 21:08:15 UTC (rev 1442)
@@ -459,8 +459,8 @@
areaCell => grid % areaCell % array
allocate( &
- drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges+1), &
- du2TopOfCell(nVertLevels+1,nCells+1), du2TopOfEdge(nVertLevels+1,nEdges+1))
+ drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges), &
+ du2TopOfCell(nVertLevels+1,nCells+1), du2TopOfEdge(nVertLevels+1,nEdges))
! compute density of parcel displaced to next deeper z-level,
! in state % rhoDisplaced
</font>
</pre>