<p><b>mpetersen@lanl.gov</b> 2011-05-18 15:01:25 -0600 (Wed, 18 May 2011)</p><p>Linear Coriolis term, f*u^{perp} was separated out of the eta*u term. This revision works with config_n_bcl_iter_beg>1. Tested using config_n_bcl_iter_beg = 1, 2, 10.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-18 20:48:10 UTC (rev 842)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-18 21:01:25 UTC (rev 843)
@@ -500,6 +500,9 @@
do while (associated(block))
allocate(G(block % mesh % nVertLevels))
+ ! Put f*uBcl^{perp} in uNew as a work variable
+ call compute_f_uperp(block % state % time_levs(2) % state , block % mesh)
+
do iEdge=1,block % mesh % nEdges
G = 0.0 ! could put this after with G(maxleveledgetop+1:nvertlevels)=0
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
@@ -507,7 +510,12 @@
! This soon, with coriolis:
! G(k) = -fEdge(iEdge)*uBclPerp(k,iEdge) + tend_u(k,iEdge)
- G(k) = block % tend % u % array (k,iEdge)
+ ! Here uNew is a work variable containing -fEdge(iEdge)*uBclPerp(k,iEdge)
+! mrp temp new one:
+ G(k) = block % state % time_levs(2) % state % u % array (k,iEdge) &
+ + block % tend % u % array (k,iEdge)
+! mrp old one:
+! G(k) = block % tend % u % array (k,iEdge)
! GBtrForcing(iEdge) = GBtrForcing(iEdge) + hOld(k,iEdge)* G(k))
@@ -528,6 +536,7 @@
enddo
deallocate(G)
+
block => block % next
end do
@@ -536,10 +545,17 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % uBcl % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+! printing:
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+print *, 'ubcl, j= ',j,minval(block % state % time_levs(2) % state % uBcl % array(:,1:nEdges)),&
+ maxval(block % state % time_levs(2) % state % uBcl % array(:,1:nEdges))
+
block => block % next
end do
- enddo ! config_n_bcl_iter
+ enddo ! do j=1,config_n_bcl_iter
!print *, '5'
@@ -586,7 +602,7 @@
call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
! printing:
-! nCells = block % mesh % nCells
+ nCells = block % mesh % nCells
print *, 'h_tend ',minval(block % tend % h % array(1,1:nCells)),&
maxval(block % tend % h % array(1,1:nCells))
print *, 'tr_tend',minval(block % tend % tracers % array(3,1,1:nCells)),&
@@ -1315,7 +1331,119 @@
end subroutine compute_tend_u
+ ! Put f*uBcl^{perp} in u as a work variable
+ subroutine compute_f_uperp(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute height and normal wind tendencies, as well as diagnostic variables
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tendencies for prognostic variables
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
+! Some of these variables can be removed, but at a later time.
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
+ vertex1, vertex2, eoe, i, j
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, q, &
+ upstream_bias, wTopEdge, rho0Inv, r
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ zMidZLevel, zTopZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, uBcl, v, pressure, &
+ tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
+ MontPot, wTop, divergence, vertViscTopOfEdge
+ type (dm_info) :: dminfo
+
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
+ edgesOnEdge, edgesOnVertex
+ real (kind=RKIND) :: u_diffusion
+ real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+
+
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+ h => s % h % array
+ u => s % u % array
+ uBcl => s % uBcl % array
+ v => s % v % array
+ wTop => s % wTop % array
+ h_edge => s % h_edge % array
+ circulation => s % circulation % array
+ vorticity => s % vorticity % array
+ divergence => s % divergence % array
+ ke => s % ke % array
+ ke_edge => s % ke_edge % array
+ pv_edge => s % pv_edge % array
+ MontPot => s % MontPot % array
+ pressure => s % pressure % array
+
+ weightsOnEdge => grid % weightsOnEdge % array
+ kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ cellsOnVertex => grid % cellsOnVertex % array
+ verticesOnEdge => grid % verticesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnVertex => grid % edgesOnVertex % array
+ dcEdge => grid % dcEdge % array
+ dvEdge => grid % dvEdge % array
+ areaCell => grid % areaCell % array
+ areaTriangle => grid % areaTriangle % array
+ h_s => grid % h_s % array
+ fVertex => grid % fVertex % array
+ fEdge => grid % fEdge % array
+ zMidZLevel => grid % zMidZLevel % array
+ zTopZLevel => grid % zTopZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelVertexBot => grid % maxLevelVertexBot % array
+
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nEdgesSolve = grid % nEdgesSolve
+ nVertices = grid % nVertices
+ nVertLevels = grid % nVertLevels
+
+ !
+ ! Put f*uBcl^{perp} in u as a work variable
+ !
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ u(k,iEdge) = 0.0
+ do j = 1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ u(k,iEdge) = u(k,iEdge) + weightsOnEdge(j,iEdge) * uBcl(k,eoe) * fEdge(eoe)
+ end do
+ end do
+ end do
+
+ end subroutine compute_f_uperp
+
subroutine compute_scalar_tend(tend, s, d, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
@@ -1934,7 +2062,7 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
- integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree, fCoef
real (kind=RKIND), dimension(:), pointer :: &
@@ -2226,6 +2354,18 @@
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
+ if (trim(config_time_integration) == 'RK4') then
+ ! for RK4, PV is really PV = (eta+f)/h
+ fCoef = 1
+ elseif (trim(config_time_integration) == 'higdon_split' &
+ .or.trim(config_time_integration) == 'higdon_unsplit') then
+ ! for higdon, PV is eta/h because f is added separately to the momentum forcing.
+! mrp temp, new should be:
+ fCoef = 0
+! old, for testing:
+! fCoef = 1
+ end if
+
do iVertex = 1,nVertices
do k=1,maxLevelVertexBot(iVertex)
h_vertex = 0.0
@@ -2234,7 +2374,7 @@
end do
h_vertex = h_vertex / areaTriangle(iVertex)
- pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
+ pv_vertex(k,iVertex) = (fCoef*fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
end do
</font>
</pre>