<p><b>qchen3@fsu.edu</b> 2013-03-21 10:36:37 -0600 (Thu, 21 Mar 2013)</p><p>BRANCH COMMIT<br>
<br>
In preparation for a merge of the trunk to the ocean_projects/gm_implement branch.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/gm_implement/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/Registry 2013-03-21 16:08:51 UTC (rev 2642)
+++ branches/ocean_projects/gm_implement/src/core_ocean/Registry 2013-03-21 16:36:37 UTC (rev 2643)
@@ -393,9 +393,14 @@
var scratch real grad_rho ( nVertLevels nEdges ) 0 - grad_rho scratch - -
var scratch real grad_rho_interface ( nVertLevelsP1 nEdges ) 0 - grad_rho_interface scratch - -
+var scratch real grad_rho_constZ ( nVertLevelsP1 nEdges ) 0 - grad_rho_constZ scratch - -
var scratch real rhoz ( nVertLevelsP1 nCells ) 0 - rhoz scratch - -
var scratch real rhoz_edge ( nVertLevelsP1 nEdges ) 0 - rhoz_edge scratch - -
var scratch real rhoz_center ( nVertLevels nCells ) 0 - rhoz_center scratch - -
var scratch real grad_zMid ( nVertLevels nEdges ) 0 - grad_zMid scratch - -
var scratch real grad_zMid_interface ( nVertLevelsP1 nEdges ) 0 - grad_zMid_interface scratch - -
+var scratch real gmStreamFunc ( nVertLevelsP1 nEdges ) 0 - gmStreamFunc scratch - -
+var scratch real tridiagA ( nVertLevels ) 0 - tridiagA scratch - -
+var scratch real tridiagB ( nVertLevels ) 0 - tridiagB scratch - -
+var scratch real tridiagC ( nVertLevels ) 0 - tridiagC scratch - -
Modified: branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F
===================================================================
--- branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F 2013-03-21 16:08:51 UTC (rev 2642)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F 2013-03-21 16:36:37 UTC (rev 2643)
@@ -27,6 +27,7 @@
! Private module variables
!
!--------------------------------------------------------------------
+ private :: tridiagonal_solve
contains
@@ -37,9 +38,9 @@
type(scratch_type), intent(inout) :: scratch
real(kind=RKIND), dimension(:,:), pointer :: rho, zMid, uBolusGM, hEddyFlux, h_edge, &
- grad_rho, grad_rho_interface, rhoz_edge, rhoz, grad_zMid, grad_zMid_interface, &
- slopeRelative, k33, rhoz_center
- real(kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge
+ grad_rho, grad_rho_interface, grad_rho_constZ, rhoz_edge, rhoz, grad_zMid, &
+ grad_zMid_interface, slopeRelative, k33, rhoz_center, gmStreamFunc
+ real(kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, tridiagA, tridiagB, tridiagC
integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
integer, dimension(:,:), pointer :: cellsOnEdge
integer :: k, iEdge, nEdges, cell1, cell2, iCell, nCells
@@ -47,12 +48,18 @@
call mpas_allocate_scratch_field(scratch % grad_rho, .True.)
call mpas_allocate_scratch_field(scratch % grad_rho_interface, .True.)
+ call mpas_allocate_scratch_field(scratch % grad_rho_constZ, .True.)
call mpas_allocate_scratch_field(scratch % rhoz, .True.)
call mpas_allocate_scratch_field(scratch % rhoz_edge, .True.)
call mpas_allocate_scratch_field(scratch % rhoz_center, .True.)
call mpas_allocate_scratch_field(scratch % grad_zMid, .True.)
call mpas_allocate_scratch_field(scratch % grad_zMid_interface, .True.)
+ call mpas_allocate_scratch_field(scratch % gmStreamFunc, .True.)
+ call mpas_allocate_scratch_field(scratch % tridiagA, .True.)
+ call mpas_allocate_scratch_field(scratch % tridiagB, .True.)
+ call mpas_allocate_scratch_field(scratch % tridiagC, .True.)
+
rho => s % rho % array
zMid => s % zMid % array
uBolusGM => s % uBolusGM % array
@@ -61,13 +68,19 @@
h_edge => s % h_edge % array
hEddyFlux => s % hEddyFlux % array
zMid => s % zMid % array
+
grad_rho => scratch % grad_rho % array
grad_rho_interface => scratch % grad_rho_interface % array
+ grad_rho_constZ => scratch % grad_rho_constZ % array
rhoz => scratch % rhoz % array
rhoz_edge => scratch % rhoz_edge % array
rhoz_center => scratch % rhoz_center % array
grad_zMid => scratch % grad_zMid % array
grad_zMid_interface => scratch % grad_zMid_interface % array
+ gmStreamFunc => scratch % gmStreamFunc % array
+ tridiagA => scratch % tridiagA % array
+ tridiagB => scratch % tridiagB % array
+ tridiagC => scratch % tridiagC % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelCell => grid % maxLevelCell % array
@@ -155,6 +168,9 @@
end do
end do
+ ! Compute the density gradient along constant z surfaces.
+ grad_rho_constZ(:,:) = grad_rho_interface(:,:) - rhoz_edge(:,:) * grad_zMid_interface(:,:)
+
! Compute slopeRelative at cell edge and layer interface
! QC Note: We may need to overwrite the invalid or unrealistically large values in the variable.
slopeRelative(:,:) = -(1.0 - config_diapycnal_diff) * grad_rho_interface(:,:) / rhoz_edge(:,:)
@@ -179,9 +195,16 @@
end do
end do
-
- call ocn_gm_compute_hEddyFlux(s, grid)
+ !
+ !Compute the stream function
+ !
+ ! Construct the tridiagonal matrices
+
+ ! Call the tridiagonal solver
+
+
+
if (config_vert_coord_movement .EQ. 'isopycnal') then
do iEdge = 1, nEdges
@@ -289,4 +312,49 @@
end subroutine ocn_get_h_kappa_q!}}}
+ subroutine tridiagonal_solve(a,b,c,r,x,n) !{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Solve the matrix equation Ax=r for x, where A is tridiagonal.
+ ! A is an nxn matrix, with:
+ ! a sub-diagonal, filled from 1:n-1 (a(1) appears on row 2)
+ ! b diagonal, filled from 1:n
+ ! c sup-diagonal, filled from 1:n-1 (c(1) apears on row 1)
+ !
+ ! Input: a,b,c,r,n
+ !
+ ! Output: x
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ integer,intent(in) :: n
+ real (KIND=RKIND), dimension(n), intent(in) :: a,b,c,r
+ real (KIND=RKIND), dimension(n), intent(out) :: x
+ real (KIND=RKIND), dimension(n) :: bTemp,rTemp
+ real (KIND=RKIND) :: m
+ integer i
+
+ call mpas_timer_start("tridiagonal_solve")
+
+ ! Use work variables for b and r
+ bTemp(1) = b(1)
+ rTemp(1) = r(1)
+
+ ! First pass: set the coefficients
+ do i = 2,n
+ m = a(i-1)/bTemp(i-1)
+ bTemp(i) = b(i) - m*c(i-1)
+ rTemp(i) = r(i) - m*rTemp(i-1)
+ end do
+
+ x(n) = rTemp(n)/bTemp(n)
+ ! Second pass: back-substition
+ do i = n-1, 1, -1
+ x(i) = (rTemp(i) - c(i)*x(i+1))/bTemp(i)
+ end do
+
+ call mpas_timer_stop("tridiagonal_solve")
+
+ end subroutine tridiagonal_solve !}}}
+
end module ocn_gm
</font>
</pre>