<p><b>mpetersen@lanl.gov</b> 2011-06-02 10:41:17 -0600 (Thu, 02 Jun 2011)</p><p>More changes to the Higdon code. Unsplit results are identical. Split not yet working.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/Makefile
===================================================================
--- branches/ocean_projects/timesplitting_mrp/Makefile 2011-06-02 13:06:27 UTC (rev 866)
+++ branches/ocean_projects/timesplitting_mrp/Makefile 2011-06-02 16:41:17 UTC (rev 867)
@@ -86,7 +86,7 @@
"CC = gcc" \
"SFC = ifort" \
"SCC = gcc" \
- "FFLAGS = -real-size 64 -convert big_endian -check all" \
+ "FFLAGS = -real-size 64 -O0 -convert big_endian " \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-06-02 13:06:27 UTC (rev 866)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-06-02 16:41:17 UTC (rev 867)
@@ -176,6 +176,7 @@
# Tendency variables: neither read nor written to any files
var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
var persistent real tend_h ( nVertLevels nCells Time ) 1 - h tend - -
+var persistent real tend_ssh ( nCells Time ) 1 - ssh tend - -
var persistent real tend_temperature ( nVertLevels nCells Time ) 1 - temperature tend tracers dynamics
var persistent real tend_salinity ( nVertLevels nCells Time ) 1 - salinity tend tracers dynamics
var persistent real tend_tracer1 ( nVertLevels nCells Time ) 1 - tracer1 tend tracers testing
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-02 13:06:27 UTC (rev 866)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-06-02 16:41:17 UTC (rev 867)
@@ -385,7 +385,7 @@
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
- real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, tend_ssh, tend_uBtr
+ real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
!print *, '1'
@@ -604,7 +604,7 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END baroclinic iterations on linear Coriolis term
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-print *, '5'
+!print *, '5'
@@ -645,7 +645,7 @@
block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(1) % state % ssh % array(iCell)
enddo
-print *, '6'
+!print *, '6'
do iEdge=1,block % mesh % nEdges
@@ -663,8 +663,8 @@
block => block % next
end do ! block
-print *, '7'
-
+!print *, '7'
+print *, 'Barotropic subcycle'
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN Barotropic subcycle loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -674,9 +674,9 @@
block => domain % blocklist
do while (associated(block))
- allocate(tend_ssh(block % mesh % nCells))
+!! mrp del allocate(tend_ssh(block % mesh % nCells))
-print *, '7.1'
+!print *, '7.1'
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -697,15 +697,15 @@
flux = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
* block % mesh % dvEdge % array(iEdge) &
* (sshEdge + hSum)
- tend_ssh(cell1) = tend_ssh(cell1) - flux
- tend_ssh(cell2) = tend_ssh(cell2) + flux
+ block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux
+ block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux
block % state % time_levs(1) % state % FBtr % array(iEdge) &
= block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ flux
end do
-print *, '7.2'
-
+!print *, '7.2'
+
do iCell=1,block % mesh % nCells
! \zeta_{n+(j+1)/J} = \zeta_{n+j/J} + \frac{\Delta t}{J} \left( -
@@ -714,7 +714,7 @@
block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
= block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ dt/config_n_btr_subcycles &
- * tend_ssh(iCell) / block % mesh % areaCell % array (iCell)
+ * block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
! sshNew(iCell) = sshNew(iCell) + sshSubcycleNew(iCell)
block % state % time_levs(2) % state % ssh % array(iCell) &
@@ -723,12 +723,12 @@
end do
-print *, '7.3'
- deallocate(tend_ssh)
-print *, '7.4'
+!print *, '7.3'
+! mrp del deallocate(tend_ssh)
+!print *, '7.4'
block => block % next
end do ! block
-print *, '8'
+!print *, '8'
! boundary update on \zeta_{n+(j+1)/J}
block => domain % blocklist
@@ -743,7 +743,7 @@
block => block % next
end do ! block
-print *, '9'
+!print *, '9'
! compute new barotropic velocity
@@ -803,9 +803,17 @@
enddo !do BtrCorIter=1,2
+! 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'
+!print *, '10'
! boundary update on \ubtr_{n+(j+1)/J}
block => domain % blocklist
@@ -822,7 +830,7 @@
! advance time pointers
oldBtrSubcycleTime = mod(oldBtrSubcycleTime,2)+1
newBtrSubcycleTime = mod(newBtrSubcycleTime,2)+1
-print *, '11 oldBtrSubcycleTime, newBtrSubcycleTime',oldBtrSubcycleTime, newBtrSubcycleTime
+!print *, '11 oldBtrSubcycleTime, newBtrSubcycleTime',oldBtrSubcycleTime, newBtrSubcycleTime
end do ! j=1,config_n_btr_subcycles
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -871,14 +879,16 @@
block => domain % blocklist
do while (associated(block))
- allocate(tend_ssh(block % mesh % nCells))
+ allocate(uTemp(block % mesh % nVertLevels))
+! allocate(tend_ssh(block % mesh % nCells))
+
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- tend_ssh(cell1) = tend_ssh(cell1) - block % state % time_levs(1) % state % FBtr % array(iEdge)
- tend_ssh(cell2) = tend_ssh(cell2) + block % state % time_levs(1) % state % FBtr % array(iEdge)
+ block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - block % state % time_levs(1) % state % FBtr % array(iEdge)
+ block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + block % state % time_levs(1) % state % FBtr % array(iEdge)
end do
@@ -887,7 +897,7 @@
block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
= block % state % time_levs(1) % state % ssh % array(iCell) &
+ dt &
- * tend_ssh(iCell) / block % mesh % areaCell % array (iCell)
+ * block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
end do
! Now can compare sshSubcycleNEW (big step using summed fluxes) with
! sshSubcycleOLD (individual steps to get there)
@@ -950,7 +960,9 @@
- deallocate(tend_ssh)
+! deallocate(tend_ssh)
+ deallocate(uTemp)
+
block => block % next
end do ! block
@@ -1065,6 +1077,20 @@
end do ! iCell
endif ! higdon_unsplit
+
+ ! compute u*, the velocity for tendency terms. Put in uNew.
+ ! uBclNew is at time n+1/2 here.
+ ! This overwrites u^{tr}, the tracer transport velocity, which was in uNew.
+ ! The following must occur after call compute_scalar_tend
+ do iEdge=1,block % mesh % nEdges
+ 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)
+ end do ! iEdge
+
+ ! mrp 110512 I really only need this to compute h_edge, density, pressure.
+ ! I can par this down later.
+ call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
elseif (higdon_step == config_n_ts_iter) then
@@ -1085,21 +1111,6 @@
endif ! higdon_step
-
- ! compute u*, the velocity for tendency terms. Put in uNew.
- ! uBclNew is at time n+1/2 here.
- ! This overwrites u^{tr}, the tracer transport velocity, which was in uNew.
- ! The following must occur after call compute_scalar_tend
- do iEdge=1,block % mesh % nEdges
- 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)
- end do ! iEdge
-
- ! mrp 110512 I really only need this to compute h_edge, density, pressure.
- ! I can par this down later.
- call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
-
block => block % next
end do
</font>
</pre>