<p><b>qchen3@fsu.edu</b> 2013-03-25 10:55:55 -0600 (Mon, 25 Mar 2013)</p><p>BRANCH COMMIT<br>
<br>
Branch commit for gm_implement.<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-22 21:31:57 UTC (rev 2665)
+++ branches/ocean_projects/gm_implement/src/core_ocean/Registry 2013-03-25 16:55:55 UTC (rev 2666)
@@ -57,6 +57,7 @@
namelist real standard_GM config_h_kappa 0.0
namelist real standard_GM config_h_kappa_q 0.0
namelist real standard_GM config_diapycnal_diff 0.0000001
+namelist real standard_GM config_gravWaveSpeed_trunc 0.3
namelist logical Rayleigh_damping config_Rayleigh_friction .false.
namelist real Rayleigh_damping config_Rayleigh_damping_coeff 0.0
@@ -403,4 +404,5 @@
var scratch real tridiagA ( nVertLevels ) 0 - tridiagA scratch - -
var scratch real tridiagB ( nVertLevels ) 0 - tridiagB scratch - -
var scratch real tridiagC ( nVertLevels ) 0 - tridiagC scratch - -
+var scratch real rightHandSide ( nVertLevels ) 0 - rightHandSide 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-22 21:31:57 UTC (rev 2665)
+++ branches/ocean_projects/gm_implement/src/core_ocean/mpas_ocn_gm.F 2013-03-25 16:55:55 UTC (rev 2666)
@@ -3,6 +3,7 @@
use mpas_grid_types
use mpas_configure
use mpas_timer
+ use mpas_constants
implicit none
private
@@ -39,12 +40,12 @@
real(kind=RKIND), dimension(:,:), pointer :: rho, zMid, uBolusGM, hEddyFlux, h_edge, &
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
+ grad_zMid_interface, slopeRelative, k33, rhoz_center, gmStreamFunc, BruntVaisalaFreqTop
+ real(kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, tridiagA, tridiagB, tridiagC, rightHandSide
integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
integer, dimension(:,:), pointer :: cellsOnEdge
- integer :: k, iEdge, nEdges, cell1, cell2, iCell, nCells
- real(kind=RKIND) :: h1, h2, areaEdge
+ integer :: k, iEdge, nEdges, cell1, cell2, iCell, nCells, N
+ real(kind=RKIND) :: h1, h2, areaEdge, c, BruntVaisalaFreqTopEdge
call mpas_allocate_scratch_field(scratch % grad_rho, .True.)
call mpas_allocate_scratch_field(scratch % grad_rho_interface, .True.)
@@ -55,6 +56,7 @@
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 % rightHandSide, .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.)
@@ -68,6 +70,7 @@
h_edge => s % h_edge % array
hEddyFlux => s % hEddyFlux % array
zMid => s % zMid % array
+ BruntVaisalaFreqTop => s % BruntVaisalaFreqTop % array
grad_rho => scratch % grad_rho % array
grad_rho_interface => scratch % grad_rho_interface % array
@@ -78,6 +81,7 @@
grad_zMid => scratch % grad_zMid % array
grad_zMid_interface => scratch % grad_zMid_interface % array
gmStreamFunc => scratch % gmStreamFunc % array
+ rightHandSide => scratch % rightHandSide % array
tridiagA => scratch % tridiagA % array
tridiagB => scratch % tridiagB % array
tridiagC => scratch % tridiagC % array
@@ -198,37 +202,66 @@
!
!Compute the stream function
!
+ gmStreamFunc(:,:) = 0.0
+ c = config_gravWaveSpeed_trunc**2
+ do iEdge = 1, nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
- ! Construct the tridiagonal matrices
+ ! Construct the tridiagonal matrix
+ ! First row
+ k = 2
+ BruntVaisalaFreqTopEdge = 0.5 * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
+ tridiagB(k-1) = - 0.5*c/(h_edge(k-1,iEdge)*h_edge(k,iEdge)) - BruntVaisalaFreqTopEdge
+ tridiagC(k-1) = 2.*c/h_edge(k,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
+ rightHandSide(k-1) = config_h_kappa * gravity / config_rho0 * grad_rho_constZ(k,iEdge)
- ! Call the tridiagonal solver
-
+ ! Second to next to the last rows
+ do k = 3, maxLevelEdgeTop(iEdge)-1
+ BruntVaisalaFreqTopEdge = 0.5 * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2))
+ tridiagA(k-2) = 2.*c/h_edge(k-1,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
+ tridiagB(k-1) = - 0.5*c/(h_edge(k-1,iEdge)*h_edge(k,iEdge)) - BruntVaisalaFreqTopEdge
+ tridiagC(k-1) = 2.*c/h_edge(k,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
+ rightHandSide(k-1) = config_h_kappa * gravity / config_rho0 * grad_rho_constZ(k,iEdge)
+ end do
-
- if (config_vert_coord_movement .EQ. 'isopycnal') then
+ ! Last row
+ k = maxLevelEdgeTop(iEdge)
+ tridiagA(k-2) = 2.*c/h_edge(k-1,iEdge)/(h_edge(k-1,iEdge)+h_edge(k,iEdge))
+ tridiagB(k-1) = - 0.5*c/(h_edge(k-1,iEdge)*h_edge(k,iEdge)) - BruntVaisalaFreqTopEdge
+ rightHandSide(k-1) = config_h_kappa * gravity / config_rho0 * grad_rho_constZ(k,iEdge)
- do iEdge = 1, nEdges
- do k = 1, maxLevelEdgeTop(iEdge)
- uBolusGM(k,iEdge) = hEddyFlux(k,iEdge)/h_edge(k,iEdge)
- end do
- end do
+ ! Total number of rows
+ N = maxLevelEdgeTop(iEdge) - 1
+
+ ! Call the tridiagonal solver
+ call tridiagonal_solve(tridiagA, tridiagB, tridiagC, rightHandSide, gmStreamFunc(2:maxLevelEdgeTop(iEdge),iEdge), N)
- else
+ end do
- ! Nothing for now for all other grid types (zlevel, zstar, ztilde)
- uBolusGM(:,:) = 0.000000001 ! A tiny number just for testing for the moment
+ ! Compute uBolusGM from the stream function
+ do iEdge = 1, nEdges
+ do k = 1, maxLevelEdgeTop(iEdge)
+ uBolusGM(k,iEdge) = (gmStreamFunc(k,iEdge) - gmStreamFunc(k+1,iEdge)) / h_edge(k,iEdge)
+ end do
+ end do
- end if
+
call mpas_deallocate_scratch_field(scratch % grad_rho)
call mpas_deallocate_scratch_field(scratch % grad_rho_interface)
call mpas_deallocate_scratch_field(scratch % rhoz)
call mpas_deallocate_scratch_field(scratch % rhoz_edge)
call mpas_deallocate_scratch_field(scratch % rhoz_center)
call mpas_deallocate_scratch_field(scratch % grad_zMid)
- call mpas_deallocate_scratch_field(scratch % grad_zMid_interface)
+ call mpas_deallocate_scratch_field(scratch % gmStreamFunc)
+ call mpas_deallocate_scratch_field(scratch % rightHandSide)
+ call mpas_deallocate_scratch_field(scratch % tridiagA)
+ call mpas_deallocate_scratch_field(scratch % tridiagB)
+ call mpas_deallocate_scratch_field(scratch % tridiagC)
+
end subroutine ocn_gm_compute_uBolus!}}}
subroutine ocn_gm_compute_slopeRelative(s, grid, scratch)
</font>
</pre>