<p><b>mpetersen@lanl.gov</b> 2011-05-18 09:02:58 -0600 (Wed, 18 May 2011)</p><p>Higdon unsplit now works. This is without baroclinic iterations. Tracer3=1.0 and is conserved to 14 digits. Works on hex 120km mesh and quad x3 mesh. Tested with config_n_ts_iter=1,10,100.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-05-17 19:25:13 UTC (rev 839)
+++ branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-05-18 15:02:58 UTC (rev 840)
@@ -1,16 +1,13 @@
&sw_model
config_test_case = 0
- config_time_integration = 'RK4'
- config_dt = 120.0
- config_ntimesteps = 100
+ config_time_integration = 'higdon_unsplit'
+ config_dt = 10.0
+ config_ntimesteps =100
config_output_interval = 100
config_stats_interval = 100
/
&io
- config_input_name = grid.nc
- config_output_name = output.nc
- config_restart_name = restart.nc
/
&restart
@@ -23,6 +20,14 @@
config_vert_grid_type = 'zlevel'
config_rho0 = 1000
/
+
+&timestep_higdon
+ config_n_ts_iter = 10
+ config_n_bcl_iter_beg = 1
+ config_n_bcl_iter_mid = 1
+ config_n_bcl_iter_end = 1
+ config_compute_tr_midstage = .true.
+/
&hmix
config_h_mom_eddy_visc2 = 1.0e5
config_h_mom_eddy_visc4 = 0.0
@@ -30,15 +35,15 @@
config_h_tracer_eddy_diff4 = 0.0
/
&vmix
- config_vert_visc_type = 'rich'
- config_vert_diff_type = 'rich'
- config_implicit_vertical_mix = .true.
+ config_vert_visc_type = 'const'
+ config_vert_diff_type = 'const'
+ config_implicit_vertical_mix = .false.
config_convective_visc = 1.0
config_convective_diff = 1.0
/
&vmix_const
- config_vert_visc = 2.5e-5
- config_vert_diff = 2.5e-5
+ config_vert_visc = 2.5e-5
+ config_vert_diff = 2.5e-5
/
&vmix_rich
config_bkrd_vert_visc = 1.0e-4
@@ -57,15 +62,15 @@
config_eos_type = 'jm'
/
&advection
- config_vert_tracer_adv = 'spline'
- config_vert_tracer_adv_order = 3
+ config_vert_tracer_adv = 'stencil'
+ config_vert_tracer_adv_order = 2
config_tracer_adv_order = 2
config_thickness_adv_order = 2
config_positive_definite = .false.
config_monotonic = .false.
/
&restore
-config_restoreTS = .false.
-config_restoreT_timescale = 90.0
-config_restoreS_timescale = 90.0
+ config_restoreTS = .false.
+ config_restoreT_timescale = 90.0
+ config_restoreS_timescale = 90.0
/
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-17 19:25:13 UTC (rev 839)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-18 15:02:58 UTC (rev 840)
@@ -430,13 +430,6 @@
! change to maxLevelCell % array(iCell) ?
do k=1,block % mesh % nVertLevels
- ! mrp 110516 efficiency note: Can we get rid of all this mult and divide
- ! of tracers by h here? Is it needed?
- ! here the tracer array now has tracer*h.
- block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % h % array(k,iCell)
-
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
= block % state % time_levs(1) % state % tracers % array(:,k,iCell)
end do
@@ -449,9 +442,9 @@
!print *, '2'
- !
- ! Begin iteration
- !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! BEGIN large iteration loop
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do higdon_step = 1, config_n_ts_iter
! --- update halos for diagnostic variables
@@ -528,9 +521,9 @@
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
!{\bf u}'_{k,n+1} = {\bf u}'_{k,n} + \Delta t \left({\bf G}_k - {\overline {\bf G}}\right)
- block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
- = block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
- + dt * (G(k) - block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
+ block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
+ = block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
+ + dt * (G(k) - block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
enddo
enddo
@@ -576,9 +569,11 @@
do while (associated(block))
do iEdge=1,block % mesh % nEdges
- block % state % time_levs(2) % state % u % array(:,iEdge) &
+ block % state % time_levs(2) % state % u % array(:,iEdge) &
= block % state % time_levs(2) % state % uBtr % array(iEdge) &
- + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
+ + 0.5*( &
+ block % state % time_levs(1) % state % uBcl % array(:,iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(:,iEdge) )
end do ! iEdge
@@ -591,7 +586,9 @@
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)),&
maxval(block % tend % tracers % array(3,1,1:nCells))
@@ -639,8 +636,9 @@
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- = (block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- + dt * block % tend % tracers % array(:,k,iCell) &
+ = ( block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ * block % state % time_levs(1) % state % h % array(k,iCell) &
+ + dt * block % tend % tracers % array(:,k,iCell) &
) / block % state % time_levs(2) % state % h % array(k,iCell)
end do
end do
@@ -652,7 +650,7 @@
block => block % next
end do
- end do
+ end do ! higdon_step = 1, config_n_ts_iter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -662,7 +660,7 @@
!
block => domain % blocklist
do while (associated(block))
-print *, '9'
+!print *, '9'
! printing:
nCells = block % mesh % nCells
nEdges = block % mesh % nEdges
@@ -679,18 +677,18 @@
! 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 *, 'hnew ',minval(block % state % time_levs(2) % state % h % array(:,1:nCells)),&
- maxval(block % state % time_levs(2) % state % h % array(:,1:nCells))
+!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 *, '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)),&
- maxval(block % state % time_levs(2) % state % tracers % array(3,1,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)),&
+! maxval(block % state % time_levs(2) % state % tracers % array(3,1,1:nCells))
! block % state % time_levs(2) % state % tracers % array(3,1,1:nCells) = domain % dminfo % my_proc_id
! printing end
@@ -931,13 +929,11 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1,min(1,maxLevelEdgeTop(iEdge))
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
-! mrp temp
- tend_h(k,cell1) = tend_h(k,cell1) - flux
+ tend_h(k,cell1) = tend_h(k,cell1) - flux
tend_h(k,cell2) = tend_h(k,cell2) + flux
end do
end do
do iCell=1,nCells
-! mrp temp
tend_h(1,iCell) = tend_h(1,iCell) / areaCell(iCell)
end do
@@ -949,11 +945,10 @@
! Vertical advection computed for top layer of a z grid only.
if (config_vert_grid_type.eq.'zlevel') then
do iCell=1,nCells
-! mrp temp
-! tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
+ tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
end do
endif ! coordinate type
-
+
end subroutine compute_tend_h
subroutine compute_tend_u(tend, s, d, grid)
@@ -1326,6 +1321,8 @@
!
! Input: s - current model state
! grid - grid metadata
+ ! note: the variable s % tracers really contains the tracers,
+ ! not tracers*h
!
! Output: tend - computed scalar tendencies
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1432,7 +1429,6 @@
do iTracer=1,num_tracers
tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
-! mrp temp
tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
end do
@@ -1709,10 +1705,9 @@
do iCell=1,nCellsSolve
do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
-! mrp temp
-! tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
-! - ( wTop(k ,iCell)*tracerTop(iTracer,k ,iCell) &
-! - wTop(k+1,iCell)*tracerTop(iTracer,k+1,iCell))
+ tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
+ - ( wTop(k ,iCell)*tracerTop(iTracer,k ,iCell) &
+ - wTop(k+1,iCell)*tracerTop(iTracer,k+1,iCell))
end do
end do
end do
@@ -1742,11 +1737,6 @@
do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
! \kappa_2 </font>
<font color="gray">abla \phi on edge
-! my new one:
-! tracer_turb_flux = config_h_tracer_eddy_diff2 &
-! *( tracers(iTracer,k,cell2)/h(k,cell2) &
-! - tracers(iTracer,k,cell1)/h(k,cell1))/dcEdge(iEdge)
-! orig:
tracer_turb_flux = config_h_tracer_eddy_diff2 &
*( tracers(iTracer,k,cell2) &
- tracers(iTracer,k,cell1))/dcEdge(iEdge)
@@ -1791,11 +1781,11 @@
do iTracer=1,num_tracers
delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
+ dvEdge(iEdge)*h_edge(k,iEdge) &
- *(tracers(iTracer,k,cell2)/h(k,cell2) - tracers(iTracer,k,cell1)/h(k,cell1)) &
+ *(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) &
/dcEdge(iEdge) * boundaryMask(k,iEdge)
delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
- dvEdge(iEdge)*h_edge(k,iEdge) &
- *(tracers(iTracer,k,cell2)/h(k,cell2) - tracers(iTracer,k,cell1)/h(k,cell1)) &
+ *(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) &
/dcEdge(iEdge) * boundaryMask(k,iEdge)
end do
end do
@@ -1857,9 +1847,10 @@
do iTracer=1,num_tracers
! compute \kappa_v d\phi/dz
fluxVertTop(iTracer,k) = vertDiffTopOfCell(k,iCell) &
- * ( tracers(iTracer,k-1,iCell)/h(k-1,iCell) &
- - tracers(iTracer,k ,iCell)/h(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
</font>
</pre>