<p><b>mpetersen@lanl.gov</b> 2013-03-20 07:12:06 -0600 (Wed, 20 Mar 2013)</p><p>branch commit: remove coment lines, revert advection comments to previous version.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_advection.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_advection.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_advection.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -172,6 +172,9 @@
if ( iCell .ne. grid % CellsOnEdge % array(1,iEdge)) &
angle_2d(i) = angle_2d(i) - pii
+! xp(i) = grid % xCell % array(cell_list(i)) - grid % xCell % array(iCell)
+! yp(i) = grid % yCell % array(cell_list(i)) - grid % yCell % array(iCell)
+
xp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * cos(angle_2d(i))
yp(i) = grid % dcEdge % array(grid % EdgesOnCell % array(i,iCell)) * sin(angle_2d(i))
@@ -302,6 +305,10 @@
else
thetae(2,grid % EdgesOnCell % array (i,iCell)) = thetae_tmp
end if
+! else
+!
+! xe(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (xv1 + xv2)
+! ye(grid % EdgesOnCell % array (i,iCell)) = 0.5 * (yv1 + yv2)
end if
@@ -350,6 +357,13 @@
cos2t = cos2t**2
sin2t = sin2t**2
+! do j=1,n
+!
+! deriv_two(j,1,iEdge) = 2.*xe(iEdge)*xe(iEdge)*bmatrix(4,j) &
+! + 2.*xe(iEdge)*ye(iEdge)*bmatrix(5,j) &
+! + 2.*ye(iEdge)*ye(iEdge)*bmatrix(6,j)
+! end do
+
if (iCell == grid % cellsOnEdge % array(1,iEdge)) then
do j=1,n
deriv_two(j,1,iEdge) = 2.*cos2t*bmatrix(4,j) &
@@ -371,6 +385,25 @@
if (debug) stop
+
+! write(0,*) ' check for deriv2 coefficients, iEdge 4 '
+!
+! iEdge = 4
+! j = 1
+! iCell = grid % cellsOnEdge % array(1,iEdge)
+! write(0,*) ' j, icell, coef ',j,iCell,deriv_two(j,1,iEdge)
+! do j=2,7
+! write(0,*) ' j, icell, coef ',j,grid % CellsOnCell % array(j-1,iCell),deriv_two(j,1,iEdge)
+! end do
+!
+! j = 1
+! iCell = grid % cellsOnEdge % array(2,iEdge)
+! write(0,*) ' j, icell, coef ',j,iCell,deriv_two(j,2,iEdge)
+! do j=2,7
+! write(0,*) ' j, icell, coef ',j,grid % CellsOnCell % array(j-1,iCell),deriv_two(j,2,iEdge)
+! end do
+! stop
+
end subroutine ocn_initialize_advection_rk!}}}
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_diagnostics.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_diagnostics.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_diagnostics.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -176,8 +176,6 @@
! Compute height on cell edges at velocity locations
! Namelist options control the order of accuracy of the reconstructed h_edge value
!
- ! mrp 101115 note: in order to include flux boundary conditions, we will need to
- ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
! initialize h_edge to avoid divide by zero and NaN problems.
h_edge = -1.0e34
@@ -244,13 +242,10 @@
! Compute v (tangential) velocities
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
- ! mrp 101115 note: in order to include flux boundary conditions,
- ! the following loop may need to change to maxLevelEdgeBot
do k = 1,maxLevelEdgeTop(iEdge)
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
end do
-
end do
!
@@ -290,8 +285,6 @@
!
! Compute ke on cell edges at velocity locations for quadratic bottom drag.
!
- ! mrp 101025 efficiency note: we could get rid of ke_edge completely by
- ! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -393,7 +386,6 @@
! Pressure
! This section must be after computing rho
!
- ! dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
if (config_pressure_gradient_type.eq.'MontgomeryPotential') then
! For Isopycnal model.
@@ -484,7 +476,6 @@
if (config_h_kappa .GE. epsilon(0D0)) then
call ocn_gm_compute_uBolus(s,grid)
else
- ! mrp efficiency note: if uBolusGM is guaranteed to be zero, this can be removed.
uBolusGM = 0.0
end if
@@ -644,10 +635,7 @@
type (state_type), intent(inout) :: s !< Input/Output: State information
type (mesh_type), intent(in) :: grid !< Input: Grid information
-! 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, cell1, cell2, eoe, i, j, k
-
integer :: nEdgesSolve
real (kind=RKIND), dimension(:), pointer :: fEdge
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, u, uBcl
@@ -827,15 +815,13 @@
endif
if (trim(config_time_integrator) == 'RK4') then
- ! for RK4, PV is really PV = (eta+f)/h
+ ! For RK4, PV includes f: PV = (eta+f)/h.
fCoef = 1
elseif (trim(config_time_integrator) == 'split_explicit' &
.or.trim(config_time_integrator) == 'unsplit_explicit') then
- ! for split explicit, PV is eta/h because f is added separately to the momentum forcing.
- ! mrp temp, new should be:
+ ! For split explicit, PV is eta/h because the Coriolis term
+ ! is added separately to the momentum tendencies.
fCoef = 0
- ! old, for testing:
- ! fCoef = 1
end if
end subroutine ocn_diagnostics_init!}}}
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_mpas_core.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_mpas_core.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -169,16 +169,10 @@
block => block % next
end do
- ! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
- ! input arguement into mpas_init. Ask about that later. For now, there will be
- ! no initial statistics write.
-
if (config_write_stats_on_startup) then
call mpas_timer_start("global diagnostics", .false., globalDiagTimer)
call ocn_compute_global_diagnostics(domain, 1 , 0, dt)
call mpas_timer_stop("global diagnostics", globalDiagTimer)
-! call mpas_output_state_init(output_obj, domain, "OUTPUT")
-! call ocn_write_output_frame(output_obj, output_frame, domain)
endif
current_outfile_frames = 0
@@ -296,7 +290,6 @@
block % state % time_levs(1) % state % uReconstructMeridional % array &
)
-!TDR
call mpas_reconstruct(mesh, mesh % u_src % array, &
block % state % time_levs(1) % state % uSrcReconstructX % array, &
block % state % time_levs(1) % state % uSrcReconstructY % array, &
@@ -304,18 +297,12 @@
block % state % time_levs(1) % state % uSrcReconstructZonal % array, &
block % state % time_levs(1) % state % uSrcReconstructMeridional % array &
)
-!TDR
- ! 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.
+ ! initialize velocities and tracers on land to be zero.
block % mesh % areaCell % array(block % mesh % nCells+1) = -1.0e34
do iEdge=1,block % mesh % nEdges
- ! mrp 101115 note: in order to include flux boundary conditions, the following
- ! line will need to change. Right now, set boundary edges between land and
- ! water to have zero velocity.
block % state % time_levs(1) % state % u % array( &
block % mesh % maxLevelEdgeTop % array(iEdge)+1 &
:block % mesh % maxLevelEdgeBot % array(iEdge), iEdge) = 0.0
@@ -323,15 +310,11 @@
block % state % time_levs(1) % state % u % array( &
block % mesh % maxLevelEdgeBot % array(iEdge)+1: &
block % mesh % nVertLevels,iEdge) = 0.0
-! mrp changed to 0
-! block % mesh % nVertLevels,iEdge) = -1e34
end do
do iCell=1,block % mesh % nCells
block % state % time_levs(1) % state % tracers % array( &
:, block % mesh % maxLevelCell % array(iCell)+1 &
:block % mesh % nVertLevels,iCell) = 0.0
-! mrp changed to 0
-! :block % mesh % nVertLevels,iCell) = -1e34
end do
do i=2,nTimeLevs
@@ -598,11 +581,6 @@
nVertLevels = block % mesh % nVertLevels
num_tracers = size(tracers, dim=1)
- ! mrp 120208 right now hZLevel is in the grid.nc file.
- ! We would like to transition to using refBottomDepth
- ! as the defining variable instead, and will transition soon.
- ! When the transition is done, hZLevel can be removed from
- ! registry and the following four lines deleted.
refBottomDepth(1) = hZLevel(1)
do k = 2,nVertLevels
refBottomDepth(k) = refBottomDepth(k-1) + hZLevel(k)
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tendency.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tendency.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -187,7 +187,6 @@
!
! velocity tendency: start accumulating tendency terms
!
- ! mrp 110516 efficiency: could remove next line and have first tend_u operation not be additive
tend_u(:,:) = 0.0
if(config_disable_u_all_tend) return
@@ -230,8 +229,6 @@
!
! velocity tendency: forcing and bottom drag
!
- ! mrp 101115 note: in order to include flux boundary conditions, we will need to
- ! know the bottom edge with nonzero velocity and place the drag there.
call mpas_timer_start("forcings", .false., velForceTimer)
call ocn_vel_forcing_tend(grid, u, u_src, ke_edge, h_edge, tend_u, err)
@@ -304,10 +301,6 @@
!
! tracer tendency: horizontal advection term -div( h \phi u)
!
- ! mrp 101115 note: in order to include flux boundary conditions, we will need to
- ! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the compute_solve_diagnostics
- ! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
- ! tracer_edge at the boundary will also need to be defined for flux boundaries.
! Monotonoic Advection, or standard advection
call mpas_timer_start("adv", .false., tracerHadvTimer)
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -271,7 +271,7 @@
! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
! be computed. For kpp, more variables may be needed. Either way, this
! could be made more efficient by only computing what is needed for the
- ! implicit vmix routine that follows. mrp 121023.
+ ! implicit vmix routine that follows.
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
@@ -314,7 +314,6 @@
block % state % time_levs(2) % state % uReconstructMeridional % array &
)
-!TDR
call mpas_reconstruct(block % mesh, block % mesh % u_src % array, &
block % state % time_levs(2) % state % uSrcReconstructX % array, &
block % state % time_levs(2) % state % uSrcReconstructY % array, &
@@ -322,7 +321,6 @@
block % state % time_levs(2) % state % uSrcReconstructZonal % array, &
block % state % time_levs(2) % state % uSrcReconstructMeridional % array &
)
-!TDR
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_time_integration_split.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_time_integration_split.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -420,7 +420,6 @@
! config_btr_gam1_uWt1= 1 flux = uBtrNew*H
! config_btr_gam1_uWt1=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam1_uWt1= 0 flux = uBtrOld*H
- ! mrp 120201 efficiency: could we combine the following edge and cell loops?
do iCell = 1, block % mesh % nCells
do i = 1, block % mesh % nEdgesOnCell % array(iCell)
@@ -560,7 +559,6 @@
! config_btr_gam3_uWt2= 1 flux = uBtrNew*H
! config_btr_gam3_uWt2=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
- ! mrp 120201 efficiency: could we combine the following edge and cell loops?
do iCell = 1, block % mesh % nCells
do i = 1, block % mesh % nEdgesOnCell % array(iCell)
@@ -771,13 +769,8 @@
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !TDR: it seems almost trivial to hold off on doing T, S and rho updates until the
- !TDR: dycore time step is complete. we might want to take this opportunity to clean-up
- !TDR: Stage3 in order to faciliate the testing of not doing tracer updates after this code is committed to trunk.
- !TDR: at this point, I am suggesting just pushing some of this code into subroutines.
- !TDR: see comments farther down
-
- ! dwj: 02/22/12 splitting thickness and tracer tendency computations and halo updates to allow monotonic advection.
+ ! Thickness tendency computations and thickness halo updates are completed before tracer
+ ! tendency computations to allow monotonic advection.
block => domain % blocklist
do while (associated(block))
@@ -819,9 +812,6 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (split_explicit_step < config_n_ts_iter) then
- !TDR: should we move this code into a subroutine called "compute_intermediate_value_at_midtime"
- !TDR: this could be within a contains statement in this routine
-
! Only need T & S for earlier iterations,
! then all the tracers needed the last time through.
do iCell=1,block % mesh % nCells
@@ -872,8 +862,8 @@
end do ! iEdge
- ! mrp 110512 I really only need this to compute h_edge, density, pressure, and SSH
- ! I can par this down later.
+ ! Efficiency note: We really only need this to compute h_edge, density, pressure, and SSH
+ ! in this diagnostics solve.
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -883,9 +873,6 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
elseif (split_explicit_step == config_n_ts_iter) then
- !TDR: should we move this code into a subroutine called "compute_final_values_at_nplus1"?
- !TDR: this could be within a contains statement in this routine
-
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
@@ -945,7 +932,7 @@
! For Richardson vertical mixing, only rho, h_edge, and ke_edge need to
! be computed. For kpp, more variables may be needed. Either way, this
! could be made more efficient by only computing what is needed for the
- ! implicit vmix routine that follows. mrp 121023.
+ ! implicit vmix routine that follows.
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
@@ -990,7 +977,6 @@
block % state % time_levs(2) % state % uReconstructMeridional % array &
)
-!TDR
call mpas_reconstruct(block % mesh, block % mesh % u_src % array, &
block % state % time_levs(2) % state % uSrcReconstructX % array, &
block % state % time_levs(2) % state % uSrcReconstructY % array, &
@@ -998,7 +984,6 @@
block % state % time_levs(2) % state % uSrcReconstructZonal % array, &
block % state % time_levs(2) % state % uSrcReconstructMeridional % array &
)
-!TDR
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tracer_vadv.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tracer_vadv.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tracer_vadv.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -123,11 +123,6 @@
err = 0
- ! mrp 120202 efficiency note:
- ! The following if statement is not needed, since wTop is set to
- ! zero for isopycnal coordinates. This if statment saves flops
- ! for isopycnal coordinates. However, if the loops are pushed
- ! out, we could get rid of this if statement.
if(.not.vadvOn) return
call ocn_tracer_vadv_stencil_tend(grid, h, wTop, tracers, tend, err1)
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -150,8 +150,6 @@
depthTop(k+1) = depthTop(k) + h(k,iCell)
enddo
- ! mrp 110201 efficiency note: push tracer loop down
- ! into spline subroutines to improve efficiency
do iTracer=1,num_tracers
! Place data in arrays to avoid creating new temporary arrays for every
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_vel_vadv.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_vel_vadv.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_vel_vadv.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -112,11 +112,6 @@
real (kind=RKIND) :: wTopEdge
real (kind=RKIND), dimension(:), allocatable :: w_dudzTopEdge
- ! mrp 120202 efficiency note:
- ! The following if statement is not needed, since wTop is set to
- ! zero for isopycnal coordinates. This if statment saves flops
- ! for isopycnal coordinates. However, if the loops are pushed
- ! out, we could get rid of this if statement.
if(.not.velVadvOn) return
err = 0
Modified: branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-03-19 22:45:53 UTC (rev 2629)
+++ branches/ocean_projects/comment_cleanup/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-03-20 13:12:06 UTC (rev 2630)
@@ -223,13 +223,11 @@
do iEdge = 1,nEdges
do k = 2,maxLevelEdgeTop(iEdge)
- ! mrp 110324 efficiency note: this if is inside iEdge and k loops.
+ ! efficiency note: these if statements are inside iEdge and k loops.
! Perhaps there is a more efficient way to do this.
if (RiTopOfEdge(k,iEdge)>0.0) then
vertViscTopOfEdge(k,iEdge) = vertViscTopOfEdge(k, iEdge) + config_bkrd_vert_visc &
+ config_rich_mix / (1.0 + 5.0*RiTopOfEdge(k,iEdge))**2
- ! maltrud do limiting of coefficient--should not be necessary
- ! also probably better logic could be found
if (vertViscTopOfEdge(k,iEdge) > config_convective_visc) then
vertViscTopOfEdge(k,iEdge) = config_convective_visc
end if
@@ -314,15 +312,13 @@
coef = -gravity/config_rho0/2.0
do iCell = 1,nCells
do k = 2,maxLevelCell(iCell)
- ! mrp 110324 efficiency note: this if is inside iCell and k loops.
+ ! efficiency note: these if statements are inside iEdge and k loops.
! Perhaps there is a more efficient way to do this.
if (RiTopOfCell(k,iCell)>0.0) then
vertDiffTopOfCell(k,iCell) = vertDiffTopOfCell(k, iCell) + config_bkrd_vert_diff &
+ (config_bkrd_vert_visc &
+ config_rich_mix / (1.0 + 5.0*RiTopOfCell(k,iCell))**2) &
/ (1.0 + 5.0*RiTopOfCell(k,iCell))
- ! maltrud do limiting of coefficient--should not be necessary
- ! also probably better logic could be found
if (vertDiffTopOfCell(k,iCell) > config_convective_diff) then
vertDiffTopOfCell(k,iCell) = config_convective_diff
end if
@@ -432,16 +428,6 @@
drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges), &
du2TopOfCell(nVertLevels+1,nCells+1), du2TopOfEdge(nVertLevels+1,nEdges))
- ! compute density of parcel displaced to next deeper z-level,
- ! in state % rhoDisplaced
-!maltrud make sure rho is current--check this for redundancy
-! call OcnEquationOfStateRho(grid, 'relative', 0, indexT, indexS, &
-! tracers, rho, err)
- ! mrp 110324 In order to visualize rhoDisplaced, include the following
-! call OcnEquationOfStateRho(grid, 'relative', 1, indexT, indexS, &
-! tracers, rhoDisplaced, err)
-
-
! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
drhoTopOfCell = 0.0
do iCell=1,nCells
</font>
</pre>