<p><b>mpetersen@lanl.gov</b> 2011-03-24 16:25:27 -0600 (Thu, 24 Mar 2011)</p><p>BRANCH COMMIT Update to implicit vertical mix branch. Branch is now ready for review.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_mpas_core.F 2011-03-24 21:42:29 UTC (rev 768)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_mpas_core.F 2011-03-24 22:25:27 UTC (rev 769)
@@ -81,38 +81,10 @@
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
-! mrp causes error Subscript #1 of the array INDX has value 7 which is greater than the upper bound of 6
call rbfInterp_initialize(mesh)
call init_reconstruct(mesh)
call reconstruct(block % state % time_levs(1) % state, mesh)
- ! mrp 110131 testing only
- block % mesh % u_src % array(:,:) = &
- block % mesh % u_src % array(:,:) *100
- print *, 'min,maxval(u_src) ',minval( block % mesh % u_src % array(:,:)),maxval( block % mesh % u_src % array(:,:))
-
- do iCell=1,block % mesh % nCells
- do k=1,mesh % nVertLevels
- !block % state % time_levs(1) % state % tracers % array( &
- ! 1, k ,iCell) = 31-k/1000.0 !31-k
- !block % state % time_levs(1) % state % tracers % array( &
- ! 2, k ,iCell) = 32 !+ k*0.32
- !block % state % time_levs(1) % state % tracers % array( &
- ! 3, k ,iCell) = k
- enddo
- block % state % time_levs(1) % state % tracers % array( &
- 3, :,iCell) = 1.
- block % state % time_levs(1) % state % tracers % array( &
- 4, :,iCell) = 1.
-
- do k=1,mesh % nVertLevels,2
- block % state % time_levs(1) % state % tracers % array( &
- 4, k ,iCell) = 0.
- enddo
-
- end do
- ! mrp 110131 testing only end
-
! initialize velocities and tracers on land to be -1e34
! The reconstructed velocity on land will have values not exactly
! -1e34 due to the interpolation of reconstruction.
Modified: branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F 2011-03-24 21:42:29 UTC (rev 768)
+++ branches/ocean_projects/imp_vert_mix_mrp/src/core_ocean/module_time_integration.F 2011-03-24 22:25:27 UTC (rev 769)
@@ -266,8 +266,9 @@
call compute_vertical_mix_coefficients(provis, block % diagnostics, block % mesh)
-!print '(a,10es12.4)', 'u1', minval(u(:,1:nEdges)),maxval(u(:,1:nEdges))
-! should I use nEdges or nEdgesSolve here?
+ !
+ ! Implicit vertical solve for momentum
+ !
do iEdge=1,nEdges
if (maxLevelEdgeTop(iEdge).gt.0) then
@@ -296,25 +297,15 @@
call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
-!print '(a,2i4)', 'iEdge, maxLevelEdgeTop(iEdge))', iEdge,maxLevelEdgeTop(iEdge)
-!print '(a,100es12.4)', 'vertViscTopOfEdge(:,iEdge)', vertViscTopOfEdge(:,iEdge)
-!print '(a,100es12.4)', 'h_edge', h_edge(:,iEdge)
-!print '(a,100es12.4)', 'dt', dt
-!print '(a,100es12.4)', 'A', A
-!print '(a,100es12.4)', 'C', C
-!print '(a,100es12.4)', 'u(:,iEdge)', u(:,iEdge)
-!print '(a,100es12.4)', 'uTemp ', uTemp
-!stop
-
-! mrp temp, turn off implicit on u:
u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
end if
end do
-!print '(a,10es12.4)', 'u2', minval(u(:,1:nEdges)),maxval(u(:,1:nEdges))
-
+ !
+ ! Implicit vertical solve for tracers
+ !
do iCell=1,nCells
! Compute A(k), C(k) for tracers
! mrp 110315 efficiency note: for z-level, could precompute
@@ -336,19 +327,12 @@
tracersTemp, maxLevelCell(iCell), &
nVertLevels,num_tracers)
-!mrp temp
tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
-! mrp temp remove soon:
-! do k=1,maxLevelCell(iCell)
-! tracers(:,k,iCell) = tracers(:,k,iCell) / h(k,iCell)
-! end do
-! mrp temp end
! mrp set as follows
! tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = 0
end do
-!print '(a,10es12.4)', 'tracer', minval(tracers(:,:,1:nCells)),maxval(tracers(:,:,1:nCells))
deallocate(A,C,uTemp,tracersTemp)
end if
@@ -358,7 +342,6 @@
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
-! mrp causes error Subscript #1 of the array INDX has value 7 which is greater than the upper bound of 6
call reconstruct(block % state % time_levs(2) % state, block % mesh)
block => block % next
@@ -762,28 +745,28 @@
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
-!mrp indent later
if (.not.config_implicit_vertical_mix) then
- allocate(fluxVertTop(nVertLevels+1))
- fluxVertTop(1) = 0.0
- do iEdge=1,grid % nEdgesSolve
+ allocate(fluxVertTop(nVertLevels+1))
+ fluxVertTop(1) = 0.0
+ do iEdge=1,grid % nEdgesSolve
- do k=2,maxLevelEdgeTop(iEdge)
- fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &
- * ( u(k-1,iEdge) - u(k,iEdge) ) &
- * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
- enddo
- fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
+ do k=2,maxLevelEdgeTop(iEdge)
+ fluxVertTop(k) = vertViscTopOfEdge(k,iEdge) &
+ * ( u(k-1,iEdge) - u(k,iEdge) ) &
+ * 2 / (h_edge(k-1,iEdge) + h_edge(k,iEdge))
+ enddo
+ fluxVertTop(maxLevelEdgeTop(iEdge)+1) = 0.0
- do k=1,maxLevelEdgeTop(iEdge)
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- + (fluxVertTop(k) - fluxVertTop(k+1)) &
- / h_edge(k,iEdge)
- enddo
+ do k=1,maxLevelEdgeTop(iEdge)
+ tend_u(k,iEdge) = tend_u(k,iEdge) &
+ + (fluxVertTop(k) - fluxVertTop(k+1)) &
+ / h_edge(k,iEdge)
+ enddo
- end do
- deallocate(fluxVertTop)
+ end do
+ deallocate(fluxVertTop)
endif
+
end subroutine compute_tend
@@ -1290,33 +1273,32 @@
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
-!mrp indent later
if (.not.config_implicit_vertical_mix) then
- allocate(fluxVertTop(num_tracers,nVertLevels+1))
- fluxVertTop(:,1) = 0.0
- do iCell=1,nCellsSolve
+ allocate(fluxVertTop(num_tracers,nVertLevels+1))
+ fluxVertTop(:,1) = 0.0
+ do iCell=1,nCellsSolve
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! compute \kappa_v d\phi/dz
- fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
- * (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell) )&
- * 2 / (h(k-1,iCell) + h(k,iCell))
- enddo
- enddo
- fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ ! compute \kappa_v d\phi/dz
+ fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
+ * (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell) )&
+ * 2 / (h(k-1,iCell) + h(k,iCell))
+ enddo
+ enddo
+ fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! This is h d/dz( fluxVertTop) but h and dz cancel, so
- ! reduces to delta( fluxVertTop)
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
- enddo
- enddo
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,num_tracers
+ ! This is h d/dz( fluxVertTop) but h and dz cancel, so
+ ! reduces to delta( fluxVertTop)
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
+ enddo
+ enddo
- enddo ! iCell loop
- deallocate(fluxVertTop)
+ enddo ! iCell loop
+ deallocate(fluxVertTop)
endif
! print some diagnostics - for debugging
@@ -1725,9 +1707,8 @@
if (config_vert_grid_type.eq.'zlevel') then
call equation_of_state(s,grid,-1)
-! mrp temp, just to visualize rhoDisplaced
- call equation_of_state(s, grid, 1)
-! mrp temp, just to visualize rhoDisplaced end
+ ! mrp 110324 In order to vissualize rhoDisplaced, include the following
+ !call equation_of_state(s, grid, 1)
endif
@@ -1846,11 +1827,9 @@
! mrp clean out these after subroutine complete
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, i, j
+ integer :: nCells, nEdges, nVertices, nVertLevels
- integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
-
real (kind=RKIND), dimension(:,:), allocatable:: &
drhoTopOfCell, drhoTopOfEdge, &
du2TopOfCell, du2TopOfEdge
@@ -1863,22 +1842,13 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
hZLevel, zTopZLevel
- real (kind=RKIND), dimension(:,:), pointer :: &
- v, w, pressure,&
- circulation, vorticity, ke, ke_edge, MontPot, wTop, &
- pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- temperature, salinity
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- real (kind=RKIND), dimension(:), allocatable:: pTop
- real (kind=RKIND), dimension(:,:), allocatable:: div_u
character :: c1*6
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
- verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
- boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot, maxLevelVertexTop
+ integer, dimension(:,:), pointer :: &
+ cellsOnEdge
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
@@ -1894,67 +1864,20 @@
RiTopOfEdge => d % RiTopOfEdge % array
RiTopOfCell => d % RiTopOfCell % array
-
-!
-!
-! mmmmmmmmmmmmmmmmmmmmmmm mrp delete variables later
-!
-!
-
-
- v => s % v % array
- wTop => s % wTop % 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
- pv_vertex => s % pv_vertex % array
- pv_cell => s % pv_cell % array
- gradPVn => s % gradPVn % array
- gradPVt => s % gradPVt % array
- tracers => s % tracers % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
-
-
zTopZLevel => grid % zTopZLevel % array
-
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % 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
- hZLevel => grid % hZLevel % array
- deriv_two => grid % deriv_two % array
- maxLevelCell => grid % maxLevelCell % array
- maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelEdgeBot => grid % maxLevelEdgeBot % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
- maxLevelVertexTop => grid % maxLevelVertexTop % array
-
nCells = grid % nCells
nEdges = grid % nEdges
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- vertexDegree = grid % vertexDegree
- boundaryEdge => grid % boundaryEdge % array
- boundaryCell => grid % boundaryCell % array
-
!
! Compute Richardson Number
!
@@ -1967,11 +1890,6 @@
! compute density of parcel displaced to surface, $\rho^*$,
! in state % rhoDisplaced
call equation_of_state(s, grid, 1)
-! mrp delete:
-!print '(a,10es12.4)', 'temperature', minval(tracers(1,2:10,1:nCells)),maxval(tracers(1,2:10,1:nCells))
-!print '(a,10es12.4)', 'salinity', minval(tracers(2,2:10,1:nCells)),maxval(tracers(2,2:10,1:nCells))
-!print '(a,10es12.4)', 'rhoDisplaced', minval(rhoDisplaced(2:10,1:nCells)),maxval(rhoDisplaced(2:10,1:nCells))
-!print '(a,10es12.4)', 'h_edge', minval(h_edge(2:10,1:nEdges)),maxval(h_edge(2:10,1:nEdges))
! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
drhoTopOfCell = 0.0
@@ -2018,10 +1936,6 @@
end do
end do
-!test this interp
-
-! are there issues with updating proc boundaries here?
-
! compute RiTopOfEdge using drhoTopOfEdge and du2TopOfEdge
! coef = -g/rho_0/2
RiTopOfEdge = 0.0
@@ -2046,17 +1960,6 @@
end do
end do
-!print '(a,10es12.4)', 'RiTopOfCell', minval(RiTopOfCell(2:10,1:nCells)),maxval(RiTopOfCell(2:10,1:nCells))
-!print '(a,10es12.4)', 'RiTopOfEdge', minval(RiTopOfEdge(2:10,1:nEdges)),maxval(RiTopOfEdge(2:10,1:nEdges))
-!print *, 'gravity',gravity
-!print '(a,10es12.4)', 'drhoTopOfCell', minval(drhoTopOfCell(2:10,1:nCells)),maxval(drhoTopOfCell(2:10,1:nCells))
-!print '(a,10es12.4)', 'du2TopOfCell', minval(du2TopOfCell(2:10,1:nCells)),maxval(du2TopOfCell(2:10,1:nCells))
-!print '(a,10es12.4)', 'h', minval(h(2:10,1:nCells)),maxval(h(2:10,1:nCells))
-
-
-! what about boundary update?
-
-
deallocate(drhoTopOfCell, drhoTopOfEdge, &
du2TopOfCell, du2TopOfEdge)
endif
@@ -2162,9 +2065,6 @@
endif
-!print '(a,10es12.4)', 'vertDiffTopOfCell', minval(vertDiffTopOfCell(2:10,1:nCells)),maxval(vertDiffTopOfCell(2:10,1:nCells))
-!print '(a,10es12.4)', 'vertViscTopOfEdge', minval(vertViscTopOfEdge(2:10,1:nEdges)),maxval(vertViscTopOfEdge(2:10,1:nEdges))
-
end subroutine compute_vertical_mix_coefficients
@@ -2491,26 +2391,12 @@
DENOMK = 1.0/(BULK_MOD - p(k))
rhoPointer(k,iCell) = (unt0 + RHO_S)*BULK_MOD*DENOMK
-!mrp
-! if (k_displaced.gt.0) then
-!print '(a,2i5,10es22.14)', 'eos icell,k,T,S,rho', icell,k,TQ,SQ,rhoPointer(k,iCell)
-! endif
-!mrp
end do
end do
-!mrp
-! if (k_displaced.gt.0) then
-!print '(a,10es22.14)', 'eos temperature', minval(tracers(1,2:10,1:nCells)),maxval(tracers(1,2:10,1:nCells))
-!print '(a,10es22.14)', 'eos salinity', minval(tracers(2,2:10,1:nCells)),maxval(tracers(2,2:10,1:nCells))
-!print '(a,10es22.14)', 'eos rhoDisplaced', minval(rhoPointer(2:10,1:nCells)),maxval(rhoPointer(2:10,1:nCells))
-! stop
-! endif
-!mrp
+ deallocate(pRefEOS,p,p2)
- deallocate(pRefEOS,p,p2)
-
end subroutine equation_of_state_jm
</font>
</pre>