<p><b>mpetersen@lanl.gov</b> 2011-06-09 12:30:13 -0600 (Thu, 09 Jun 2011)</p><p>Added ability to filter Barotropic mode. Removed some print statements.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-06-09 16:46:08 UTC (rev 890)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-06-09 18:30:13 UTC (rev 891)
@@ -24,6 +24,7 @@
namelist integer timestep_higdon config_n_btr_cor_iter 1
namelist logical timestep_higdon config_compute_tr_midstage true
namelist logical timestep_higdon config_u_correction true
+namelist logical timestep_higdon config_filter_btr_mode false
namelist logical sw_model config_h_ScaleWithMesh false
namelist real hmix config_h_mom_eddy_visc2 0.0
namelist real hmix config_h_mom_eddy_visc4 0.0
@@ -203,16 +204,16 @@
var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
var persistent real ke_edge ( nVertLevels nEdges Time ) 2 - ke_edge state - -
var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 - pv_vertex state - -
-var persistent real pv_cell ( nVertLevels nCells Time ) 2 o pv_cell state - -
-var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
-var persistent real uReconstructY ( nVertLevels nCells Time ) 2 o uReconstructY state - -
-var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 o uReconstructZ state - -
+var persistent real pv_cell ( nVertLevels nCells Time ) 2 - pv_cell state - -
+var persistent real uReconstructX ( nVertLevels nCells Time ) 2 - uReconstructX state - -
+var persistent real uReconstructY ( nVertLevels nCells Time ) 2 - uReconstructY state - -
+var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 - uReconstructZ state - -
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
-var persistent real MontPot ( nVertLevels nCells Time ) 2 o MontPot state - -
-var persistent real pressure ( nVertLevels nCells Time ) 2 o pressure state - -
-var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 o wTop state - -
-var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 o rhoDisplaced state - -
+var persistent real MontPot ( nVertLevels nCells Time ) 2 - MontPot state - -
+var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
+var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
+var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
# Other diagnostic variables: neither read nor written to any files
var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
@@ -230,9 +231,9 @@
var persistent real CFLNumberGlobal ( Time ) 2 o CFLNumberGlobal state - -
# Diagnostics fields, only one time level required
-var persistent real RiTopOfCell ( nVertLevelsP1 nCells Time ) 1 o RiTopOfCell diagnostics - -
-var persistent real RiTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 o RiTopOfEdge diagnostics - -
-var persistent real vertViscTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 o vertViscTopOfEdge diagnostics - -
-var persistent real vertDiffTopOfCell ( nVertLevelsP1 nCells Time ) 1 o vertDiffTopOfCell diagnostics - -
+var persistent real RiTopOfCell ( nVertLevelsP1 nCells Time ) 1 - RiTopOfCell diagnostics - -
+var persistent real RiTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 - RiTopOfEdge diagnostics - -
+var persistent real vertViscTopOfEdge ( nVertLevelsP1 nEdges Time ) 1 - vertViscTopOfEdge diagnostics - -
+var persistent real vertDiffTopOfCell ( nVertLevelsP1 nCells Time ) 1 - vertDiffTopOfCell diagnostics - -
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-06-09 16:46:08 UTC (rev 890)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-06-09 18:30:13 UTC (rev 891)
@@ -149,8 +149,7 @@
dt = config_dt
ntimesteps = config_ntimesteps
-! mrp temp remove first output frame
-! call write_output_frame(output_obj, output_frame, domain)
+ call write_output_frame(output_obj, output_frame, domain)
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
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-06-09 16:46:08 UTC (rev 890)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-06-09 18:30:13 UTC (rev 891)
@@ -347,7 +347,7 @@
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
! mrp temp
-! call reconstruct(block % state % time_levs(2) % state, block % mesh)
+ call reconstruct(block % state % time_levs(2) % state, block % mesh)
block => block % next
end do
@@ -380,7 +380,7 @@
real (kind=RKIND) :: uhSum, hSum, grav_rho0Inv, sshEdge, flux, &
uPerp, uCorr, tracerTemp
- integer :: nCells, nEdges, nVertLevels, num_tracers, ucorr_coef
+ integer :: num_tracers, ucorr_coef
real (kind=RKIND), dimension(:,:), pointer :: &
u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -389,7 +389,6 @@
real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, hNew
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
-!print *, '1'
grav_rho0Inv = gravity/config_rho0
@@ -433,27 +432,8 @@
end do
end do
-! printing:
- nCells = block % mesh % nCells
- nEdges = block % mesh % nEdges
- nVertLevels = block % mesh % nVertLevels
-!print *, 'uold ',minval(block % state % time_levs(1) % state % u % array(:,1:nEdges)),&
-! maxval(block % state % time_levs(1) % state % u % array(:,1:nEdges))
-!print *, 'uBtrOld',minval(block % state % time_levs(1) % state % uBtr % array(1:nEdges)),&
-! maxval(block % state % time_levs(1) % state % uBtr % array(1:nEdges))
-!print *, 'uBclOld',minval(block % state % time_levs(1) % state % uBcl % array(:,1:nEdges)),&
-! maxval(block % state % time_levs(1) % state % uBcl % array(:,1:nEdges))
-!print *, 'hold',minval(block % state % time_levs(1) % state % h % array(:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % h % array(:,1:nCells))
-!print *, 'Told',minval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells))
-!print *, 'Sold',minval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells))
-! printing end
-
block => block % next
end do
-!print *, '2'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -484,7 +464,6 @@
block => block % next
end do
-!print *, '3'
!
! Stage 1: Baroclinic velocity (3D) prediction, explicit with long timestep
@@ -508,13 +487,6 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do j=1,n_bcl_iter(higdon_step)
-! mrp soon: add linear Coriolis to tend_u maybe at end of this iteration?
- !compute {\bf u}_{k,n+1}^{'\;\perp} from {\bf u}_{k,n+1}'
-! call compute_uPerp(uBclNew,uBclPerp)
-! need subroutine and variable just for perp, say uBclPerp - may already be there.
-
-!print *, '4'
-
! Use this G coefficient to avoid an if statement within the iEdge loop.
if (trim(config_time_integration) == 'higdon_unsplit') then
split = 0
@@ -592,12 +564,6 @@
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
@@ -605,10 +571,8 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END baroclinic iterations on linear Coriolis term
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!print *, '5'
-
!
! Stage 2: Barotropic velocity (2D) prediction, explicitly subcycled
!
@@ -637,6 +601,12 @@
block => domain % blocklist
do while (associated(block))
+ if (config_filter_btr_mode) then
+ block % state % time_levs(1) % state % GBtrForcing % array(:) = 0.0
+ block % state % time_levs(1) % state % ssh % array(:) = 0.0
+ block % state % time_levs(1) % state % uBtr % array(:) = 0.0
+ endif
+
do iCell=1,block % mesh % nCells
! sshSubcycleOld = sshOld
block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
@@ -646,7 +616,6 @@
block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(1) % state % ssh % array(iCell)
enddo
-!print *, '6'
do iEdge=1,block % mesh % nEdges
@@ -664,8 +633,7 @@
block => block % next
end do ! block
-!print *, '7'
-!print *, 'Barotropic subcycle'
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN Barotropic subcycle loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -677,7 +645,6 @@
!! mrp del allocate(tend_ssh(block % mesh % nCells))
-!print *, '7.1'
block % tend % ssh % array(:) = 0.0
do iEdge=1,block % mesh % nEdges
@@ -707,7 +674,6 @@
= block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ flux
end do
-!print *, '7.2'
do iCell=1,block % mesh % nCells
@@ -721,12 +687,9 @@
end do
-!print *, '7.3'
! mrp del deallocate(tend_ssh)
-!print *, '7.4'
block => block % next
end do ! block
-!print *, '8'
! boundary update on \zeta_{n+(j+1)/J}
block => domain % blocklist
@@ -741,7 +704,6 @@
block => block % next
end do ! block
-!print *, '9'
block => domain % blocklist
do while (associated(block))
@@ -756,12 +718,10 @@
end do
-!print *, '7.3'
! mrp del deallocate(tend_ssh)
-!print *, '7.4'
+
block => block % next
end do ! block
-!print *, '8'
! compute new barotropic velocity
@@ -788,13 +748,11 @@
! Iterate on u two times.
-!print *, '9.1'
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
-!print *, '9.2'
! Compute -f*uPerp
uPerp = 0.0
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
@@ -803,7 +761,6 @@
* block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
* block % mesh % fEdge % array(eoe)
end do
-!print *, '9.3'
! timestep in uBtrSubcycle:
! \ubtr_{n+(j+1)/J} = \ubtr_{n+j/J} + \frac{\Delta t}{J} \left(
@@ -830,24 +787,13 @@
endif
-!print *, '9.4'
end do
-!print *, '9.5'
-! printing:
- nCells = block % mesh % nCells
- nEdges = block % mesh % nEdges
-!print *, 'uBtr, j= ',j,minval(block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(1:nEdges)),&
-! maxval(block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(1:nEdges))
-!print *, 'ssh, j= ',j,minval(block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(1:nCells)),&
-! maxval(block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(1:nCells))
block => block % next
end do ! block
-!print *, '10'
-
! boundary update on \ubtr_{n+(j+1)/J}
block => domain % blocklist
do while (associated(block))
@@ -882,7 +828,6 @@
! advance time pointers
oldBtrSubcycleTime = mod(oldBtrSubcycleTime,2)+1
newBtrSubcycleTime = mod(newBtrSubcycleTime,2)+1
-!print *, '11 oldBtrSubcycleTime, newBtrSubcycleTime',oldBtrSubcycleTime, newBtrSubcycleTime
end do ! j=1,config_n_btr_subcycles
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -954,10 +899,6 @@
end do
! Now can compare sshSubcycleNEW (big step using summed fluxes) with
! sshSubcycleOLD (individual steps to get there)
-!print *, 'ssh, by substeps',minval(block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve)), &
-! maxval(block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve))
-!print *, 'ssh, by 1 step ',minval(block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve)), &
-! maxval(block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(1:block % mesh % nCellsSolve))
! Correction velocity, u^{corr}:
!u^{corr} = \left( {\overline {\bf F}}
@@ -1042,8 +983,6 @@
-!print *, '16'
-
!
! Stage 3: Tracer, density, pressure, vertical velocity prediction
!
@@ -1059,18 +998,9 @@
call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
-! printing:
- nCells = block % mesh % nCells
-!print *, '1h_tend ',minval(block % tend % h % array(1,1:nCells)),&
-! maxval(block % tend % h % array(1,1:nCells))
-!print *, '1tr_tend',minval(block % tend % tracers % array(3,1,1:nCells)),&
-! maxval(block % tend % tracers % array(3,1,1:nCells))
-
block => block % next
end do
-!print *, '17'
-
! --- update halos for prognostic variables
block => domain % blocklist
@@ -1083,12 +1013,11 @@
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
-!print *, '18'
block => domain % blocklist
do while (associated(block))
- allocate(hNew(nVertLevels))
+ allocate(hNew(block % mesh % nVertLevels))
if (trim(config_time_integration) == 'higdon_unsplit') then
@@ -1107,19 +1036,6 @@
endif ! higdon_unsplit
-! printing:
- nCells = block % mesh % nCells
- nEdges = block % mesh % nEdges
- nVertLevels = block % mesh % nVertLevels
-!print *, 'hnew ',minval(block % tend % h % array(:,1:nCells)),&
-! maxval(block % tend % h % array(:,1:nCells))
-!print *, 'hnew1 ',minval(block % tend % h % array(1,1:nCells)),&
-! maxval(block % tend % h % array(1,1:nCells))
-! do iCell = 1,10 !nCells
-! print '(a,i5,50e12.2)', 'htend iCell=',iCell,block % tend % h % array(:,iCell)
-! enddo
-! printing end
-
! Only need T & S for earlier iterations,
! then all the tracers needed the last time through.
if (higdon_step < config_n_ts_iter) then
@@ -1215,11 +1131,12 @@
!
block => domain % blocklist
do while (associated(block))
-!print *, '9'
+
+! nCells = block % mesh % nCells
+! nEdges = block % mesh % nEdges
+! nVertLevels = block % mesh % nVertLevels
+
! printing:
- nCells = block % mesh % nCells
- nEdges = block % mesh % nEdges
- nVertLevels = block % mesh % nVertLevels
!print *, 'uold ',minval(block % state % time_levs(1) % state % u % array(:,1:nEdges)),&
! maxval(block % state % time_levs(1) % state % u % array(:,1:nEdges))
!print *, 'hold ',minval(block % state % time_levs(1) % state % h % array(:,1:nCells)),&
@@ -1230,16 +1147,16 @@
! maxval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells))
!print *, 'Sold ',minval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells)),&
! maxval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells))
-print *, 'unew ',minval(block % state % time_levs(2) % state % u % array(:,1:nEdges)),&
- maxval(block % state % time_levs(2) % state % u % array(:,1:nEdges))
+!print *, 'unew ',minval(block % state % time_levs(2) % state % u % array(:,1:nEdges)),&
+! maxval(block % state % time_levs(2) % state % u % array(:,1:nEdges))
!print *, 'hnew ',minval(block % state % time_levs(2) % state % h % array(:,1:nCells)),&
! maxval(block % state % time_levs(2) % state % h % array(:,1:nCells))
-print *, 'hnew1 ',minval(block % state % time_levs(2) % state % h % array(1,1:nCells)),&
- maxval(block % state % time_levs(2) % state % h % array(1,1:nCells))
-print *, 'Tnew ',minval(block % state % time_levs(2) % state % tracers % array(1,:,1:nCells)),&
- maxval(block % state % time_levs(2) % state % tracers % array(1,:,1:nCells))
-print *, 'Snew ',minval(block % state % time_levs(2) % state % tracers % array(2,:,1:nCells)),&
- maxval(block % state % time_levs(2) % state % tracers % array(2,:,1:nCells))
+!print *, 'hnew1 ',minval(block % state % time_levs(2) % state % h % array(1,1:nCells)),&
+! maxval(block % state % time_levs(2) % state % h % array(1,1:nCells))
+!print *, 'Tnew ',minval(block % state % time_levs(2) % state % tracers % array(1,:,1:nCells)),&
+! maxval(block % state % time_levs(2) % state % tracers % array(1,:,1:nCells))
+!print *, 'Snew ',minval(block % state % time_levs(2) % state % tracers % array(2,:,1:nCells)),&
+! maxval(block % state % time_levs(2) % state % tracers % array(2,:,1:nCells))
!print *, 'Tr1Old1',minval(block % state % time_levs(1) % state % tracers % array(3,1,1:nCells)),&
! maxval(block % state % time_levs(1) % state % tracers % array(3,1,1:nCells))
!print *, 'Tr1New1',minval(block % state % time_levs(2) % state % tracers % array(3,1,1:nCells)),&
@@ -1305,6 +1222,7 @@
vertDiffTopOfCell => block % diagnostics % vertDiffTopOfCell % array
maxLevelCell => block % mesh % maxLevelCell % array
maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+
! dividing by h above for higdon.
! do iCell=1,nCells
@@ -1314,15 +1232,15 @@
! end do
if (config_implicit_vertical_mix) then
- allocate(A(nVertLevels),C(nVertLevels),uTemp(nVertLevels), &
- tracersTemp(num_tracers,nVertLevels))
+ allocate(A(block % mesh % nVertLevels),C(block % mesh % nVertLevels),uTemp(block % mesh % nVertLevels), &
+ tracersTemp(num_tracers,block % mesh % nVertLevels))
call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
!
! Implicit vertical solve for momentum
!
- do iEdge=1,nEdges
+ do iEdge=1,block % mesh % nEdges
if (maxLevelEdgeTop(iEdge).gt.0) then
! Compute A(k), C(k) for momentum
@@ -1351,7 +1269,7 @@
call tridiagonal_solve(A,C,A,u(:,iEdge),uTemp,maxLevelEdgeTop(iEdge))
u(1:maxLevelEdgeTop(iEdge),iEdge) = uTemp(1:maxLevelEdgeTop(iEdge))
- u(maxLevelEdgeTop(iEdge)+1:nVertLevels,iEdge) = 0.0
+ u(maxLevelEdgeTop(iEdge)+1:block % mesh % nVertLevels,iEdge) = 0.0
end if
end do
@@ -1359,7 +1277,7 @@
!
! Implicit vertical solve for tracers
!
- do iCell=1,nCells
+ do iCell=1,block % mesh % nCells
! Compute A(k), C(k) for tracers
! mrp 110315 efficiency note: for z-level, could precompute
! -2.0*dt/(h(k)_h(k+1))/h(k) in setup
@@ -1378,10 +1296,10 @@
call tridiagonal_solve_mult(A,C,A,tracers(:,:,iCell), &
tracersTemp, maxLevelCell(iCell), &
- nVertLevels,num_tracers)
+ block % mesh % nVertLevels,num_tracers)
tracers(:,1:maxLevelCell(iCell),iCell) = tracersTemp(:,1:maxLevelCell(iCell))
- tracers(:,maxLevelCell(iCell)+1:nVertLevels,iCell) = -1e34
+ tracers(:,maxLevelCell(iCell)+1:block % mesh % nVertLevels,iCell) = -1e34
end do
deallocate(A,C,uTemp,tracersTemp)
end if
@@ -1392,8 +1310,7 @@
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
-!
-! call reconstruct(block % state % time_levs(2) % state, block % mesh)
+ call reconstruct(block % state % time_levs(2) % state, block % mesh)
block => block % next
end do
</font>
</pre>