<p><b>mpetersen@lanl.gov</b> 2011-06-10 07:09:06 -0600 (Fri, 10 Jun 2011)</p><p>Correct the initialization for normal and config_filter_btr_mode = .true.. This conducts a splitting before the first step, in module_mpas_core.F. For config_filter_btr_mode = .true., the barotropic mode is simply filtered out.<br>
</p><hr noshade><pre><font color="gray">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 20:46:23 UTC (rev 892)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-06-10 13:09:06 UTC (rev 893)
@@ -26,6 +26,8 @@
if (.not. config_do_restart) call setup_sw_test_case(domain)
+ call compute_maxLevel(domain)
+
if (config_vert_grid_type.eq.'isopycnal') then
print *, ' Using isopycnal coordinates'
elseif (config_vert_grid_type.eq.'zlevel') then
@@ -37,8 +39,6 @@
call dmpar_abort(dminfo)
endif
- call compute_maxLevel(domain)
-
!
! Initialize core
!
@@ -327,20 +327,11 @@
hRatioZLevelKm1(k) = 0.5*hZLevel(k-1)/hMeanTopZLevel(k)
end do
-! mrp delete later - probably dont need HTotZLevelEdge variable
-! do iEdge=1,block % mesh % nEdges
-! block % mesh % HTotZLevelEdge % array(iEdge) = 0.0
-! do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
-! block % mesh % HTotZLevelEdge % array(iEdge) &
-! = block % mesh % HTotZLevelEdge % array(iEdge) &
-! + block % mesh % hZLevel % array(k)
-! enddo
-! enddo
-
! mrp 110601 For now, h is the variable saved in the restart file
! I am computing SSH here. In the future, could make smaller
! restart files for z-Level runs by saving SSH only.
do iCell=1,block % mesh % nCells
+
block % state % time_levs(1) % state % ssh % array(iCell) &
= block % state % time_levs(1) % state % h % array(1,iCell) &
- block % mesh % hZLevel % array(1)
@@ -350,7 +341,21 @@
! This is only done upon start-up.
if (trim(config_time_integration) == 'higdon_unsplit') then
block % state % time_levs(1) % state % uBtr % array(:) = 0.0
+
+ block % state % time_levs(1) % state % uBcl % array(:,:) &
+ = block % state % time_levs(1) % state % u % array(:,:)
+
elseif (trim(config_time_integration) == 'higdon_split') then
+
+ if (config_filter_btr_mode) then
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(1) % state % h % array(1,iCell) &
+ = block % mesh % hZLevel % array(1)
+
+ block % state % time_levs(1) % state % ssh % array(iCell) = 0.0
+ enddo
+ endif
+
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -371,9 +376,56 @@
+ block % mesh % hZLevel % array(k)
enddo
block % state % time_levs(1) % state % uBtr % array(iEdge) = uhSum/hsum
+
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
+ = block % state % time_levs(1) % state % u % array(k,iEdge) &
+ - block % state % time_levs(1) % state % uBtr % array(iEdge)
+ enddo
+
+ do k=block % mesh % maxLevelEdgeTop % array(iEdge)+1, block % mesh % nVertLevels
+ block % state % time_levs(1) % state % uBcl % array(k,iEdge) = 0.0
+ block % state % time_levs(1) % state % u % array(k,iEdge) = 0.0
+ enddo
enddo
+
+ if (config_filter_btr_mode) then
+ block % state % time_levs(1) % state % uBtr % array(:) = 0.0
+
+ block % state % time_levs(1) % state % u % array(:,:) &
+ = block % state % time_levs(1) % state % uBcl % array(:,:)
+ endif
+
endif
+
+!print *, '11 u ',minval(domain % blocklist % state % time_levs(1) % state % u % array(:,1:domain % blocklist % mesh % nEdgesSolve)), &
+! maxval(domain % blocklist % state % time_levs(1) % state % u % array(:,1:domain % blocklist % mesh % nEdgesSolve))
+!print *, '11 uBtr ',minval(domain % blocklist % state % time_levs(1) % state % uBtr % array(1:domain % blocklist % mesh % nEdgesSolve)), &
+! maxval(domain % blocklist % state % time_levs(1) % state % uBtr % array(1:domain % blocklist % mesh % nEdgesSolve))
+!print *, '11 uBcl ',minval(domain % blocklist % state % time_levs(1) % state % uBcl % array(:,1:domain % blocklist % mesh % nEdgesSolve)), &
+! maxval(domain % blocklist % state % time_levs(1) % state % uBcl % array(:,1:domain % blocklist % mesh % nEdgesSolve))
+
+
+! mrp temp testing - is uBcl vert sum zero?
+ do iEdge=1,block % mesh % nEdges
+ uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * block % state % time_levs(1) % state % uBcl % array(1,iEdge)
+ hSum = sshEdge + block % mesh % hZLevel % array(1)
+
+ do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
+ uhSum = uhSum + block % mesh % hZLevel % array(k) * block % state % time_levs(1) % state % uBcl % array(k,iEdge)
+ hSum = hSum + block % mesh % hZLevel % array(k)
+ enddo
+ block % state % time_levs(1) % state % FBtr % array(iEdge) = uhSum/hSum
+
+ enddo ! iEdge
+
+!print *, 'uBcl vert sum IC',minval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve)), &
+! maxval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve))
+
+! mrp temp testing - is uBcl vert sum zero? end
+
+
block => block % next
end do
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 20:46:23 UTC (rev 892)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-06-10 13:09:06 UTC (rev 893)
@@ -402,9 +402,10 @@
do iEdge=1,block % mesh % nEdges
- block % state % time_levs(1) % state % uBcl % array(:,iEdge) &
- = block % state % time_levs(1) % state % u % array(:,iEdge) &
- - block % state % time_levs(1) % state % uBtr % array(iEdge)
+ ! mrp delete soon
+ ! block % state % time_levs(1) % state % uBcl % array(:,iEdge) &
+ ! = block % state % time_levs(1) % state % u % array(:,iEdge) &
+ ! - block % state % time_levs(1) % state % uBtr % array(iEdge)
block % state % time_levs(2) % state % u % array(:,iEdge) &
= block % state % time_levs(1) % state % u % array(:,iEdge)
@@ -544,15 +545,33 @@
!{\bf u}'_{k,n+1} = {\bf u}'_{k,n} - \Delta t {\overline {\bf G}}
!{\bf u}'_{k,n+1/2} = \frac{1}{2}\left({\bf u}^{'}_{k,n} +{\bf u}'_{k,n+1}\right)
! so that uBclNew is at time n+1/2
- block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
- = 0.5*( &
- block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
- + uTemp(k) - dt * block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
-
+ block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
+ = 0.5*( &
+ block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
+ + uTemp(k) - dt * block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
enddo
+! mrp temp testing - is uBcl vert sum zero?
+! uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * block % state % time_levs(1) % state % uBcl % array(1,iEdge)
+! hSum = sshEdge + block % mesh % hZLevel % array(1)
+! do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
+! uhSum = uhSum + block % mesh % hZLevel % array(k) * block % state % time_levs(1) % state % uBcl % array(k,iEdge)
+! hSum = hSum + block % mesh % hZLevel % array(k)
+! enddo
+! block % state % time_levs(1) % state % FBtr % array(iEdge) = uhSum/hSum
+! mrp temp testing end
+
enddo ! iEdge
+! mrp temp testing - is uBcl vert sum zero?
+!print *, '12 uBcl ',minval(domain % blocklist % state % time_levs(1) % state % uBcl % array(:,1:domain % blocklist % mesh % nEdgesSolve)), &
+! maxval(domain % blocklist % state % time_levs(1) % state % uBcl % array(:,1:domain % blocklist % mesh % nEdgesSolve))
+
+!print *, 'uBcl vert sum ',minval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve)), &
+! maxval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdgesSolve))
+! mrp temp testing - is uBcl vert sum zero? end
+
+
deallocate(uTemp)
block => block % next
@@ -603,8 +622,6 @@
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
@@ -744,10 +761,6 @@
block => domain % blocklist
do while (associated(block))
-
- ! Iterate on u two times.
-
-
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
@@ -899,6 +912,10 @@
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}}
@@ -911,6 +928,14 @@
ucorr_coef = 0
endif
+
+!print *, 'uBtr ',minval(block % state % time_levs(2) % state % uBtr % array(1:block % mesh % nEdges)), &
+! maxval(block % state % time_levs(2) % state % uBtr % array(1:block % mesh % nEdges))
+!print *, 'ssh ',minval(block % state % time_levs(2) % state % ssh % array(1:block % mesh % nEdges)), &
+! maxval(block % state % time_levs(2) % state % ssh % array(1:block % mesh % nEdges))
+!print *, 'FBtr ',minval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdges)), &
+! maxval(block % state % time_levs(1) % state % FBtr % array(1:block % mesh % nEdges))
+
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -935,6 +960,8 @@
uCorr = ucorr_coef*(( block % state % time_levs(1) % state % FBtr % array(iEdge) &
/block % mesh % dvEdge % array(iEdge) &
- uhSum)/hSum)
+! mrp temp for printing
+!block % state % time_levs(1) % state % GBtrForcing % array(iEdge) = uCorr
! put u^{tr}, the velocity for tracer transport, in uNew
! mrp 060611 not sure if boundary enforcement is needed here.
@@ -953,9 +980,11 @@
= block % mesh % hZLevel % array(k)
enddo
-
end do ! iEdge
+!print *, 'uCorr ',minval(block % state % time_levs(1) % state % GBtrForcing % array(1:block % mesh % nEdgesSolve)), &
+! maxval(block % state % time_levs(1) % state % GBtrForcing % array(1:block % mesh % nEdgesSolve))
+
! Put new SSH values in h array, for the compute_scalar_tend call below.
do iCell=1,block % mesh % nCells
block % state % time_levs(2) % state % h % array(1,iCell) &
</font>
</pre>