<p><b>mpetersen@lanl.gov</b> 2011-09-07 08:45:56 -0600 (Wed, 07 Sep 2011)</p><p>Added the ability to:<br>
1. Baroclinic subcycle to 2J (that is, to timestep n+2)<br>
2. Compute SSH from barotropic flux, rather than avg of barotropic SSH<br>
3. Step forward from the average barotropic velocity, rather than last <br>
subcycle step.<br>
<br>
The new flags for these these settings are:<br>
config_btr_subcycle_loop_factor = 2<br>
config_SSH_from = 'avg_flux'<br>
config_new_btr_variables_from = 'btr_avg'<br>
Using these settings, tracers are still conserved, and I can take a long<br>
timestep (5000s) in the 100km hex gridcell basin domain. Still testing<br>
the global domain.<br>
<br>
For the previous setting, use <br>
config_btr_subcycle_loop_factor = 1<br>
config_SSH_from = 'avg_of_SSH_subcycles'<br>
config_new_btr_variables_from = 'last_subcycle'<br>
This was limited to 300s timesteps for Higdon split in teh 100km basin<br>
domain.<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-09-07 00:23:57 UTC (rev 979)
+++ branches/ocean_projects/timesplitting_mrp/namelist.input.ocean 2011-09-07 14:45:56 UTC (rev 980)
@@ -32,6 +32,9 @@
config_btr_mom_decay_time = 3600.0
config_btr_flux_coef = 1.0
config_btr_mom_eddy_visc2 = 0.0
+ config_btr_subcycle_loop_factor = 2
+ config_SSH_from = 'avg_flux'
+ config_new_btr_variables_from = 'btr_avg'
/
&hmix
config_h_mom_eddy_visc2 = 1.0e5
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-09-07 00:23:57 UTC (rev 979)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-09-07 14:45:56 UTC (rev 980)
@@ -30,6 +30,9 @@
namelist real timestep_higdon config_btr_mom_decay_time 3600.0
namelist real timestep_higdon config_btr_flux_coef 1.0
namelist real timestep_higdon config_btr_mom_eddy_visc2 0.0
+namelist integer timestep_higdon config_btr_subcycle_loop_factor 2
+namelist character timestep_higdon config_SSH_from avg_flux
+namelist character timestep_higdon config_new_btr_variables_from btr_avg
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
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-09-07 00:23:57 UTC (rev 979)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-09-07 14:45:56 UTC (rev 980)
@@ -409,6 +409,7 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
+ type (dm_info) :: dminfo
integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step, split, &
eoe, oldBtrSubcycleTime, newBtrSubcycleTime, uPerpTime, BtrCorIter, &
n_bcl_iter(config_n_ts_iter), &
@@ -417,6 +418,7 @@
type (block_type), pointer :: block
real (kind=RKIND) :: uhSum, hSum, sshEdge, flux, &
uPerp, uCorr, tracerTemp, coef
+ real (kind=RKIND), dimension(:), pointer :: sshNew
integer :: num_tracers, ucorr_coef
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -427,7 +429,6 @@
real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, hNew
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Prep variables before first iteration
@@ -672,7 +673,7 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN Barotropic subcycle loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do j=1,config_n_btr_subcycles
+ do j=1,config_n_btr_subcycles*config_btr_subcycle_loop_factor
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Barotropic subcycle: initial solve for velecity
@@ -1023,18 +1024,26 @@
do iEdge=1,block % mesh % nEdges
block % state % time_levs(1) % state % FBtr % array(iEdge) &
= block % state % time_levs(1) % state % FBtr % array(iEdge) &
- / config_n_btr_subcycles
+ / (config_n_btr_subcycles*config_btr_subcycle_loop_factor)
block % state % time_levs(2) % state % uBtr % array(iEdge) &
= block % state % time_levs(2) % state % uBtr % array(iEdge) &
- / (config_n_btr_subcycles + 1)
+ / (config_n_btr_subcycles*config_btr_subcycle_loop_factor + 1)
end do
+ if (config_SSH_from=='avg_of_SSH_subcycles') then
do iCell=1,block % mesh % nCells
block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(2) % state % ssh % array(iCell) &
- / (config_n_btr_subcycles + 1)
+ / (config_n_btr_subcycles*config_btr_subcycle_loop_factor + 1)
end do
+ elseif (config_SSH_from=='avg_flux') then
+ ! see below
+ else
+ write(0,*) 'Abort: Unknown config_SSH_from option: '&
+ //trim(config_SSH_from)
+ call dmpar_abort(dminfo)
+ endif
block => block % next
end do ! block
@@ -1062,6 +1071,7 @@
allocate(uTemp(block % mesh % nVertLevels))
+ if (config_SSH_from=='computed_from_avg_flux') then
! Accumulate fluxes in the tend % ssh variable
block % tend % ssh % array(:) = 0.0
do iEdge=1,block % mesh % nEdges
@@ -1081,12 +1091,12 @@
do iCell=1,block % mesh % nCells
! SSHnew = SSHold + dt*(-div(Flux))
- block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(1) % state % ssh % array(iCell) &
+ dt &
* block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
end do
-
+ endif
! 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)), &
@@ -1210,6 +1220,14 @@
do while (associated(block))
allocate(hNew(block % mesh % nVertLevels))
+ if (trim(config_new_btr_variables_from) == 'last_subcycle') then
+ ! This points to the last barotropic SSH subcycle
+ sshNew => block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array
+ elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
+ ! This points to the tendency variable SSH*
+ sshNew => block % state % time_levs(2) % state % ssh % array
+ endif
+
if (trim(config_time_integration) == 'higdon_unsplit') then
do iCell=1,block % mesh % nCells
@@ -1233,8 +1251,8 @@
hNew(:) = block % mesh % hZLevel % array(:)
do iCell=1,block % mesh % nCells
- hNew(1) = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
- + block % mesh % hZLevel % array(1)
+ ! sshNew is a pointer, defined above.
+ hNew(1) = sshNew(iCell) + block % mesh % hZLevel % array(1)
do k=1,block % mesh % maxLevelCell % array(iCell)
do i=1,2
! This is Phi at n+1
@@ -1286,8 +1304,8 @@
hNew(:) = block % mesh % hZLevel % array(:)
do iCell=1,block % mesh % nCells
- hNew(1) = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
- + block % mesh % hZLevel % array(1)
+ ! sshNew is a pointer, defined above.
+ hNew(1) = sshNew(iCell) + block % mesh % hZLevel % array(1)
do k=1,block % mesh % maxLevelCell % array(iCell)
do i=1,block % state % time_levs(1) % state % num_tracers
! This is Phi at n+1
@@ -1318,6 +1336,22 @@
block => domain % blocklist
do while (associated(block))
+ if (trim(config_new_btr_variables_from) == 'last_subcycle') then
+ do iEdge=1,block % mesh % nEdges
+ ! uBtrNew = uBtrSubcycleNew (old here is because counter already flipped)
+ ! This line is not needed if u is resplit at the beginning of the timestep.
+ block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
+ enddo ! iEdges
+ elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
+ ! uBtrNew from u*. this is done above, so u* is already in
+ ! block % state % time_levs(2) % state % uBtr % array(iEdge)
+ else
+ write(0,*) 'Abort: Unknown config_new_btr_variables_from: '&
+ //trim(config_time_integration)
+ call dmpar_abort(dminfo)
+ endif
+
! Recompute final u to go on to next step.
! u_{n+1} = uBtr_{n+1} + uBcl_{n+1}
! Right now uBclNew is at time n+1/2, so back compute to get uBcl at time n+1
@@ -1330,17 +1364,11 @@
do iEdge=1,block % mesh % nEdges
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
- ! uBtrNew = uBtrSubcycleNew (old here is because counter already flipped)
- ! This line is not needed if u is resplit at the beginning of the timestep.
- block % state % time_levs(2) % state % uBtr % array(iEdge) &
- = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
-
block % state % time_levs(2) % state % u % array(k,iEdge) &
- = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+2*block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
- block % state % time_levs(1) % state % uBcl % array(k,iEdge)
enddo
-
! mrp 110607 zero out velocity below land edges. efficiency: this may not be required.
do k=block % mesh % maxLevelEdgeTop % array(iEdge) + 1, block % mesh % nVertLevels
block % state % time_levs(2) % state % u % array(k,iEdge) = 0.0
@@ -1348,19 +1376,20 @@
enddo ! iEdges
-! do iEdge=1,block % mesh % nEdges
-! block % state % time_levs(2) % state % u % array(:,iEdge) &
-! = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
-! +2*block % state % time_levs(2) % state % uBcl % array(:,iEdge) &
-! - block % state % time_levs(1) % state % uBcl % array(:,iEdge)
-! end do
+ if (trim(config_time_integration) == 'higdon_split') then
- if (trim(config_time_integration) == 'higdon_split') then
+ if (trim(config_new_btr_variables_from) == 'last_subcycle') then
do iCell=1,block % mesh % nCells
! SSH for the next step is from the end of the barotropic subcycle.
block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell)
+ end do ! iCell
+ elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
+ ! sshNew from ssh*. This is done above, so ssh* is already in
+ ! block % state % time_levs(2) % state % ssh % array(iCell)
+ endif
+ do iCell=1,block % mesh % nCells
! Put new SSH values in h array, for the compute_scalar_tend call below.
block % state % time_levs(2) % state % h % array(1,iCell) &
= block % state % time_levs(2) % state % ssh % array(iCell) &
</font>
</pre>