<p><b>mpetersen@lanl.gov</b> 2012-03-26 11:07:56 -0600 (Mon, 26 Mar 2012)</p><p>Change split explicit to include uBolusGM velocity. Changed computation of wTop to use uTransport rather than u.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/gmvar/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/gmvar/src/core_ocean/mpas_ocn_tendency.F 2012-03-26 16:06:44 UTC (rev 1721)
+++ branches/ocean_projects/gmvar/src/core_ocean/mpas_ocn_tendency.F 2012-03-26 17:07:56 UTC (rev 1722)
@@ -107,13 +107,12 @@
type (state_type), intent(in) :: s !< Input: State information
type (mesh_type), intent(in) :: grid !< Input: Grid information
- real (kind=RKIND), dimension(:,:), pointer :: h_edge, u, wTop, tend_h, uTransport
+ real (kind=RKIND), dimension(:,:), pointer :: h_edge, wTop, tend_h, uTransport
integer :: err
call mpas_timer_start("ocn_tend_h")
- u => s % u % array
uTransport => s % uTransport % array
wTop => s % wTop % array
h_edge => s % h_edge % array
@@ -283,7 +282,7 @@
real (kind=RKIND), intent(in) :: dt !< Input: Time step
real (kind=RKIND), dimension(:,:), pointer :: &
- u, uTransport, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh
+ uTransport, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
@@ -291,7 +290,6 @@
call mpas_timer_start("ocn_tend_scalar")
- u => s % u % array
uTransport => s % uTransport % array
h => s % h % array
wTop => s % wTop % array
@@ -888,7 +886,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: u,h,wTop, h_edge
+ real (kind=RKIND), dimension(:,:), pointer :: uTransport,h,wTop, h_edge
real (kind=RKIND), dimension(:,:), allocatable:: div_hu
real (kind=RKIND), dimension(:), allocatable:: div_hu_btr, h_tend_col, h_weights
@@ -901,7 +899,7 @@
h => s1 % h % array
h_edge => s1 % h_edge % array
- u => s2 % u % array
+ uTransport => s2 % uTransport % array
wTop => s2 % wTop % array
areaCell => grid % areaCell % array
@@ -926,7 +924,7 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeBot(iEdge)
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ flux = uTransport(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
div_hu(k,cell1) = div_hu(k,cell1) + flux
div_hu(k,cell2) = div_hu(k,cell2) - flux
end do
Modified: branches/ocean_projects/gmvar/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/gmvar/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-26 16:06:44 UTC (rev 1721)
+++ branches/ocean_projects/gmvar/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-26 17:07:56 UTC (rev 1722)
@@ -663,9 +663,11 @@
do iEdge=1,block % mesh % nEdges
- ! This is u^{avg}
- uTemp(:) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
- + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
+ ! velocity for uCorrection is uBtr + uBcl + uBolus
+ uTemp(:) &
+ = block % state % time_levs(2) % state % uBtr % array( iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(:,iEdge) &
+ + block % state % time_levs(1) % state % uBolusGM % array(:,iEdge)
! hSum is initialized outside the loop because on land boundaries
! maxLevelEdgeTop=0, but I want to initialize hSum with a
@@ -680,19 +682,23 @@
uCorr = ucorr_coef*(( block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum)
- ! put u^{tr}, the velocity for tracer transport, in uNew
- ! mrp 060611 not sure if boundary enforcement is needed here.
- if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
- block % state % time_levs(2) % state % u % array(:,iEdge) = 0.0
- else
- do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
- block % state % time_levs(2) % state % u % array(k,iEdge) = uTemp(k) + uCorr
- enddo
- 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
- endif
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ ! uTranport = uBtr + uBcl + uBolus + uCorrection
+ ! This is u used in advective terms for h and tracers
+ ! in tendancy calls in stage 3.
+ block % state % time_levs(2) % state % uTransport % array(k,iEdge) &
+ = block % mesh % edgeMask % array(k,iEdge) &
+ *( block % state % time_levs(2) % state % uBtr % array( iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
+ + block % state % time_levs(1) % state % uBolusGM % array(k,iEdge) &
+ + uCorr )
+
+ enddo
+ do k=block % mesh % maxLevelEdgeTop % array(iEdge)+1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % uTransport % array(k,iEdge) = 0.0
+ end do
+
end do ! iEdge
deallocate(uTemp)
@@ -801,10 +807,25 @@
end do
end do ! iCell
- ! uBclNew is u'_{n+1/2}
- ! uBtrNew is {\bar u}_{avg}
- ! uNew is u^{tr}
+ do iEdge=1,block % mesh % nEdges
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+
+ ! u = uBtr + uBcl
+ ! here uBcl is at time n+1/2
+ ! This is u used in next iteration or step
+ block % state % time_levs(2) % state % u % array(k,iEdge) &
+ = block % mesh % edgeMask % array(k,iEdge) &
+ *( block % state % time_levs(2) % state % uBtr % array( iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(k,iEdge) )
+
+ enddo
+ 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
+
! mrp 110512 I really only need this to compute h_edge, density, pressure, and SSH
! I can par this down later.
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
</font>
</pre>