<p><b>mpetersen@lanl.gov</b> 2011-10-17 13:17:20 -0600 (Mon, 17 Oct 2011)</p><p>Merging bug fixes from split_explicit_timestepping branch onto trunk.<br>
</p><hr noshade><pre><font color="gray">
Property changes on: trunk/mpas
___________________________________________________________________
Modified: svn:mergeinfo
- /branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/time_manager:924-962
+ /branches/ocean_projects/imp_vert_mix_mrp:754-986
/branches/ocean_projects/split_explicit_timestepping:1044-1097
/branches/ocean_projects/vert_adv_mrp:704-745
/branches/time_manager:924-962
Modified: trunk/mpas/Makefile
===================================================================
--- trunk/mpas/Makefile 2011-10-17 18:57:45 UTC (rev 1097)
+++ trunk/mpas/Makefile 2011-10-17 19:17:20 UTC (rev 1098)
@@ -1,4 +1,3 @@
-CORE=ocean
#MODEL_FORMULATION = -DNCAR_FORMULATION
MODEL_FORMULATION = -DLANL_FORMULATION
@@ -83,7 +82,7 @@
ifort-serial:
( make all \
- "FC = ifort" \
+ "FC = mpif90" \
"CC = gcc" \
"SFC = ifort" \
"SCC = gcc" \
@@ -91,7 +90,7 @@
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
- "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
ifort-papi:
( make all \
Modified: trunk/mpas/namelist.input.ocean
===================================================================
--- trunk/mpas/namelist.input.ocean 2011-10-17 18:57:45 UTC (rev 1097)
+++ trunk/mpas/namelist.input.ocean 2011-10-17 19:17:20 UTC (rev 1098)
@@ -29,17 +29,18 @@
config_n_bcl_iter_end = 2
config_n_btr_subcycles = 10
config_n_btr_cor_iter = 2
- config_compute_tr_midstage = .true.
config_u_correction = .true.
config_filter_btr_mode = .false.
config_btr_mom_decay = .false.
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'
-/
+ config_btr_gam1_uWt1 = 0.5
+ config_btr_gam2_SSHWt1 = 1.0
+ config_btr_gam3_uWt2 = 1.0
+ config_btr_solve_SSH2 = .false.
&hmix
config_h_mom_eddy_visc2 = 1.0e5
config_h_mom_eddy_visc4 = 0.0
Modified: trunk/mpas/src/core_ocean/Registry
===================================================================
--- trunk/mpas/src/core_ocean/Registry 2011-10-17 18:57:45 UTC (rev 1097)
+++ trunk/mpas/src/core_ocean/Registry 2011-10-17 19:17:20 UTC (rev 1098)
@@ -26,16 +26,18 @@
namelist integer split_explicit_ts config_n_bcl_iter_end 4
namelist integer split_explicit_ts config_n_btr_subcycles 10
namelist integer split_explicit_ts config_n_btr_cor_iter 1
-namelist logical split_explicit_ts config_compute_tr_midstage true
namelist logical split_explicit_ts config_u_correction true
namelist logical split_explicit_ts config_filter_btr_mode false
namelist logical split_explicit_ts config_btr_mom_decay false
namelist real split_explicit_ts config_btr_mom_decay_time 3600.0
-namelist real split_explicit_ts config_btr_flux_coef 1.0
namelist real split_explicit_ts config_btr_mom_eddy_visc2 0.0
namelist integer split_explicit_ts config_btr_subcycle_loop_factor 2
namelist character split_explicit_ts config_SSH_from avg_flux
namelist character split_explicit_ts config_new_btr_variables_from btr_avg
+namelist real split_explicit_ts config_btr_gam1_uWt1 0.5
+namelist real split_explicit_ts config_btr_gam2_SSHWt1 1.0
+namelist real split_explicit_ts config_btr_gam3_uWt2 1.0
+namelist logical split_explicit_ts config_btr_solve_SSH2 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
@@ -186,9 +188,9 @@
var persistent real h ( nVertLevels nCells Time ) 2 ir h state - -
var persistent real rho ( nVertLevels nCells Time ) 2 ir rho state - -
var persistent real temperature ( nVertLevels nCells Time ) 2 iro temperature state tracers dynamics
-var persistent real salinity ( nVertLevels nCells Time ) 2 ir salinity state tracers dynamics
-var persistent real tracer1 ( nVertLevels nCells Time ) 2 ir tracer1 state tracers testing
-var persistent real tracer2 ( nVertLevels nCells Time ) 2 ir tracer2 state tracers testing
+var persistent real salinity ( nVertLevels nCells Time ) 2 iro salinity state tracers dynamics
+var persistent real tracer1 ( nVertLevels nCells Time ) 2 iro tracer1 state tracers testing
+var persistent real tracer2 ( nVertLevels nCells Time ) 2 iro tracer2 state tracers testing
# Tendency variables: neither read nor written to any files
var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
Modified: trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2011-10-17 18:57:45 UTC (rev 1097)
+++ trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2011-10-17 19:17:20 UTC (rev 1098)
@@ -249,7 +249,7 @@
! mrp changed to 0
! :block % mesh % nVertLevels,iCell) = -1e34
-! mrp 110516 temp, added just to test for conservation of tracers
+! mrp 110516, added just to test for conservation of tracers
block % state % time_levs(1) % state % tracers % array(3,:,iCell) = 1.0
end do
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2011-10-17 18:57:45 UTC (rev 1097)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2011-10-17 19:17:20 UTC (rev 1098)
@@ -148,10 +148,10 @@
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
if (config_h_mom_eddy_visc4 > 0.0) then
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % divergence % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % vorticity % array(:,:), &
block % mesh % nVertLevels, block % mesh % nVertices, &
block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
end if
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2011-10-17 18:57:45 UTC (rev 1097)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2011-10-17 19:17:20 UTC (rev 1098)
@@ -87,7 +87,7 @@
type (block_type), pointer :: block
real (kind=RKIND) :: uhSum, hSum, sshEdge, flux, &
- uPerp, uCorr, tracerTemp, coef
+ uPerp, uCorr, tracerTemp, coef, FBtr_coeff, sshCell1, sshCell2
real (kind=RKIND), dimension(:), pointer :: sshNew
integer :: num_tracers, ucorr_coef, err
@@ -159,7 +159,7 @@
block => domain % blocklist
do while (associated(block))
-! mrp 110512 not sure if I need the following three. Leave be, assume I need it.
+
call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % pv_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
@@ -350,6 +350,7 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Barotropic subcycle: initial solve for velecity
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if (config_btr_gam1_uWt1>1.0e-12) then ! only do this part if it is needed in next SSH solve
uPerpTime = oldBtrSubcycleTime
block => domain % blocklist
@@ -423,19 +424,30 @@
block => block % next
end do ! block
+ endif ! config_btr_gam1_uWt1>1.0e-12
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Barotropic subcycle: Compute thickness flux and new SSH
+ ! Barotropic subcycle: Compute thickness flux and new SSH: PREDICTOR
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
block => domain % blocklist
do while (associated(block))
block % tend % ssh % array(:) = 0.0
- ! config_btr_flux_coef sets the forward weighting of velocity in the SSH computation
- ! config_btr_flux_coef= 1 flux = uBtrNew*H
- ! config_btr_flux_coef=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
- ! config_btr_flux_coef= 0 flux = uBtrOld*H
+ if (config_btr_solve_SSH2) then
+ ! If config_btr_solve_SSH2=.true., then do NOT accumulate FBtr in this SSH predictor
+ ! section, because it will be accumulated in the SSH corrector section.
+ FBtr_coeff = 0.0
+ else
+ ! otherwise, DO accumulate FBtr in this SSH predictor section
+ FBtr_coeff = 1.0
+ endif
+
+ ! config_btr_gam1_uWt1 sets the forward weighting of velocity in the SSH computation
+ ! 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
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -445,19 +457,20 @@
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
- flux = ((1.0-config_btr_flux_coef) &
+ flux = ((1.0-config_btr_gam1_uWt1) &
* block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- + config_btr_flux_coef &
+ + config_btr_gam1_uWt1 &
* block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
- * block % mesh % dvEdge % array(iEdge) &
* (sshEdge + hSum)
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux
+ block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) &
+ - flux * block % mesh % dvEdge % array(iEdge)
+ block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) &
+ + flux * block % mesh % dvEdge % array(iEdge)
block % state % time_levs(1) % state % FBtr % array(iEdge) &
= block % state % time_levs(1) % state % FBtr % array(iEdge) &
- + flux
+ + FBtr_coeff*flux
end do
! SSHnew = SSHold + dt/J*(-div(Flux))
@@ -473,7 +486,7 @@
block => block % next
end do ! block
- ! boundary update on SSNnew
+ ! boundary update on SSHnew
block => domain % blocklist
do while (associated(block))
@@ -487,21 +500,7 @@
block => block % next
end do ! block
- block => domain % blocklist
- do while (associated(block))
- do iCell=1,block % mesh % nCells
-
- ! Accumulate SSH in running sum over the subcycles.
- block % state % time_levs(2) % state % ssh % array(iCell) &
- = block % state % time_levs(2) % state % ssh % array(iCell) &
- + block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell)
-
- end do
-
- block => block % next
- end do ! block
-
! mrp 110801 begin
! This whole section, bounded by 'mrp 110801', may be deleted later if it is found
! that barotropic del2 is not useful.
@@ -612,13 +611,22 @@
else
! uBtrNew = uBtrOld + dt*(-f*uBtroldPerp - g*grad(SSH) + G)
+
+ sshCell1 = &
+ (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1)
+
+ sshCell2 = &
+ (1-config_btr_gam2_SSHWt1)*block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ + config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
+
block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
= block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ dt/config_n_btr_subcycles *( &
uPerp &
- gravity &
- *( block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
- - block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
+ *( sshCell2 &
+ - sshCell1 )&
/block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge) &
+ block % state % time_levs(1) % state % u_diffusionBtr % array(iEdge))
@@ -662,13 +670,93 @@
end do !do BtrCorIter=1,config_n_btr_cor_iter
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: Compute thickness flux and new SSH: CORRECTOR
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if (config_btr_solve_SSH2) then
+ block => domain % blocklist
+ do while (associated(block))
+
+ block % tend % ssh % array(:) = 0.0
+
+ ! config_btr_gam3_uWt2 sets the forward weighting of velocity in the SSH computation
+ ! 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
+
+ do iEdge=1,block % mesh % nEdges
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
+ sshEdge = 0.5 &
+ *( block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ + block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
+ hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
+
+ flux = ((1.0-config_btr_gam3_uWt2) &
+ * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + config_btr_gam3_uWt2 &
+ * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * (sshEdge + hSum)
+
+ block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) &
+ - flux * block % mesh % dvEdge % array(iEdge)
+ block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) &
+ + flux * block % mesh % dvEdge % array(iEdge)
+
+ block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ = block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ + flux
+
+
+ end do
+
+ ! SSHnew = SSHold + dt/J*(-div(Flux))
+ do iCell=1,block % mesh % nCells
+
+ block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ + dt/config_n_btr_subcycles &
+ * block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
+
+ end do
+
+ block => block % next
+ end do ! block
+
+ ! boundary update on SSHnew
+ block => domain % blocklist
+ do while (associated(block))
+
+ call dmpar_exch_halo_field1dReal(domain % dminfo, &
+ block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
+ block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
+ block => block % next
+ end do ! block
+
+ endif ! config_btr_solve_SSH2
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: Accumulate running sums, advance timestep pointers
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ ! Accumulate SSH in running sum over the subcycles.
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(2) % state % ssh % array(iCell) &
+ = block % state % time_levs(2) % state % ssh % array(iCell) &
+ + block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell)
+ end do
+
! uBtrNew = uBtrNew + uBtrSubcycleNEW
! This accumulates the sum.
! If the Barotropic Coriolis iteration is limited to one, this could
! be merged with the above code.
- block => domain % blocklist
- do while (associated(block))
do iEdge=1,block % mesh % nEdges
block % state % time_levs(2) % state % uBtr % array(iEdge) &
@@ -752,11 +840,14 @@
block % tend % ssh % array(cell1) &
= block % tend % ssh % array(cell1) &
- - block % state % time_levs(1) % state % FBtr % array(iEdge)
+ - block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ * block % mesh % dvEdge % array(iEdge)
+
block % tend % ssh % array(cell2) &
= block % tend % ssh % array(cell2) &
- + block % state % time_levs(1) % state % FBtr % array(iEdge)
+ + block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ * block % mesh % dvEdge % array(iEdge)
end do
@@ -769,12 +860,6 @@
* 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)), &
-! 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 uCorr = (Flux - Sum(h u*))/H
! or, for the full latex version:
@@ -810,7 +895,6 @@
enddo
uCorr = ucorr_coef*(( block % state % time_levs(1) % state % FBtr % array(iEdge) &
- /block % mesh % dvEdge % array(iEdge) &
- uhSum)/hSum)
! put u^{tr}, the velocity for tracer transport, in uNew
@@ -818,7 +902,12 @@
if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
block % state % time_levs(2) % state % u % array(:,iEdge) = 0.0
else
- block % state % time_levs(2) % state % u % array(:,iEdge) = uTemp(:) + uCorr
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ block % state % time_levs(2) % state % u % array(k,iEdge) = uTemp(k) + uCorr
+ enddo
+ do k=block % mesh % maxLevelEdgeTop % array(iEdge)+1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % u % array(k,iEdge) = 0.0
+ enddo
endif
! Put new sshEdge values in h_edge array, for the OcnTendScalar call below.
Modified: trunk/mpas/src/core_ocean/mpas_ocn_tracer_vadv_spline.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_tracer_vadv_spline.F 2011-10-17 18:57:45 UTC (rev 1097)
+++ trunk/mpas/src/core_ocean/mpas_ocn_tracer_vadv_spline.F 2011-10-17 19:17:20 UTC (rev 1098)
@@ -168,7 +168,7 @@
if(config_vert_tracer_adv.eq.'spline') then
splineOn = .true.
- call ocn_tracer_vadv_spline2_init(err2)
+ call ocn_tracer_vadv_spline2_init(err1)
call ocn_tracer_vadv_spline3_init(err2)
err = err1 .or. err2
</font>
</pre>