<p><b>mpetersen@lanl.gov</b> 2011-05-31 16:39:08 -0600 (Tue, 31 May 2011)</p><p>BRANCH COMMIT: More progress on the Stage 2 Barotropic solver for split version. Higdon split not yet complete. Unsplit is working, as before.<br>
</p><hr noshade><pre><font color="gray">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-31 19:41:50 UTC (rev 864)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-31 22:39:08 UTC (rev 865)
@@ -371,9 +371,11 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
- integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step, split
+ integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step, split, &
+ eoe, oldBtrSubcycleTime, newBtrSubcycleTime
+
type (block_type), pointer :: block
- real (kind=RKIND) :: uSum, hSum, GSum, grav_rho0Inv, sshEdge, flux
+ real (kind=RKIND) :: uSum, hSum, GSum, grav_rho0Inv, sshEdge, flux, uPerp
integer :: nCells, nEdges, nVertLevels, num_tracers
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -381,11 +383,13 @@
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
- real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, aBtrSumCoef,bBtrSumCoef, tend_ssh, tend_uBtr
+ real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, tend_ssh, tend_uBtr
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
!print *, '1'
+ grav_rho0Inv = gravity/config_rho0
+
!
! Prep variables before first iteration
!
@@ -426,9 +430,12 @@
enddo ! iEdge
- block % state % time_levs(2) % state % h % array(:,:) &
- = block % state % time_levs(1) % state % h % array(:,:)
+ block % state % time_levs(2) % state % ssh % array(:) &
+ = block % state % time_levs(1) % state % ssh % array(:)
+! Are the following two needed?
+! block % state % time_levs(2) % state % ssh_edge % array(:,:) &
+! = block % state % time_levs(1) % state % ssh_edge % array(:,:)
block % state % time_levs(2) % state % h_edge % array(:,:) &
= block % state % time_levs(1) % state % h_edge % array(:,:)
@@ -524,7 +531,6 @@
elseif (trim(config_time_integration) == 'higdon_split') then
split = 1
endif
- grav_rho0Inv = gravity/config_rho0
block => domain % blocklist
do while (associated(block))
@@ -550,15 +556,10 @@
+ dt * (block % tend % u % array (k,iEdge) &
+ block % state % time_levs(2) % state % u % array (k,iEdge) &
+ split*grav_rho0Inv &
- *( block % state % time_levs(2) % state % h % array(1,cell2) &
- - block % state % time_levs(2) % state % h % array(1,cell1) ) &
+ *( block % state % time_levs(2) % state % ssh % array(cell2) &
+ - block % state % time_levs(2) % state % ssh % array(cell1) ) &
/block % mesh % dcEdge % array(iEdge) )
-! mrp old one:
-! G(k) = block % tend % u % array (k,iEdge)
-
-! GBtrForcing(iEdge) = GBtrForcing(iEdge) + hOld(k,iEdge)* G(k))
-
enddo
GSum = 0.0
@@ -615,66 +616,59 @@
! Stage 2: Barotropic velocity (2D) prediction, explicitly subcycled
!
+ oldBtrSubcycleTime = 1
+ newBtrSubcycleTime = 2
+
if (trim(config_time_integration) == 'higdon_unsplit') then
block => domain % blocklist
do while (associated(block))
+
block % state % time_levs(2) % state % uBtr % array(:) = 0.0
! mrp 110531 change following lines once I have pointers set up for subcycling
- block % state % time_levs(2) % state % uBtrSubcycle % array(:) = 0.0
- block % state % time_levs(1) % state % uBtrSubcycle % array(:) = 0.0
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:) = 0.0
+
block => block % next
end do ! block
elseif (trim(config_time_integration) == 'higdon_split') then
- ! note: eventually, put this in subroutine mpas_core_init
- ! using an input variable to choose the filter type.
- allocate(aBtrSumCoef(0:2*config_n_btr_subcycles), &
- bBtrSumCoef(0:2*config_n_btr_subcycles) )
- aBtrSumCoef = 1.0/(2.0*config_n_btr_subcycles+1)
- bBtrSumCoef(2*config_n_btr_subcycles) = aBtrSumCoef(2*config_n_btr_subcycles)
- do j=2*config_n_btr_subcycles-1,0,-1
- ! mrp not sure - do I need to divide by 2*config_n_btr_subcycles here?
- ! or do I divide by sum(bBtrSumCoef) at the end?
- bBtrSumCoef(j) = bBtrSumCoef(j+1) + aBtrSumCoef(j)
- end do
-
+ ! Initialize variables for barotropic subcycling
block => domain % blocklist
do while (associated(block))
do iCell=1,block % mesh % nCells
-! sshSubcycleOld = sshOld
- block % state % time_levs(1) % state % sshSubcycle % array(iCell) &
+ ! sshSubcycleOld = sshOld
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
= block % state % time_levs(1) % state % ssh % array(iCell)
-! sshNew = aBtrSumCoef(0)*sshOld
+ ! sshNew = sshOld This is the first for the summation
block % state % time_levs(2) % state % ssh % array(iCell) &
- = aBtrSumCoef(0) &
- * block % state % time_levs(1) % state % ssh % array(iCell)
+ = block % state % time_levs(1) % state % ssh % array(iCell)
enddo
do iEdge=1,block % mesh % nEdges
-! uBtrSubcycleOld = uBtrOld
- block % state % time_levs(1) % state % uBtrSubcycle % array(iEdge) &
+
+ ! uBtrSubcycleOld = uBtrOld
+ block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
= block % state % time_levs(1) % state % uBtr % array(iEdge)
-! uBtrNew = aBtrSumCoef(0)*uBtrOld
- block % state % time_levs(2) % state % uBtrSubcycle % array(iEdge) &
- = aBtrSumCoef(0) &
- * block % state % time_levs(1) % state % uBtr % array(iEdge)
+ ! uBtrNew = BtrOld This is the first for the summation
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ = block % state % time_levs(1) % state % uBtr % array(iEdge)
-! FBtr = 0
- block % state % time_levs(1) % state % FBtr % array(iEdge) = 0.0
+ ! FBtr = 0
+ block % state % time_levs(1) % state % FBtr % array(iEdge) = 0.0
enddo
block => block % next
end do ! block
- ! mrp 110519 change 2*config_n_btr_subcycles to be filter specific in the future.
- do j=1,2*config_n_btr_subcycles
+ ! Barotropic subcycle loop
+ do j=1,config_n_btr_subcycles
+ ! Compute thickness flux and new SSH
block => domain % blocklist
do while (associated(block))
@@ -685,24 +679,28 @@
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
! flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
-! FBtr(iEdge) = FBtr(iEdge) + bBtrSumCoef(j)*flux
+! FBtr(iEdge) = FBtr(iEdge) + flux
! tend_ssh(cell1) = tend_ssh(cell1) - flux
! tend_ssh(cell2) = tend_ssh(cell2) + flux
sshEdge = 0.5 &
- *( block % state % time_levs(1) % state % sshSubcycle % array(cell1) &
- + block % state % time_levs(1) % state % sshSubcycle % array(cell2) )
+ *( 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)))
! \zeta_{n+j/J}^{edge} = Interp(\zeta_{n+j/J})
! {\bf F}_j = \ubtr_{n+j/J}
! \left(\zeta_{n+j/J}^{edge} + \sum_{k=1}^{N^{edge}} \Delta z_k \right)
- flux = block % state % time_levs(1) % state % uBtrSubcycle % array(iEdge) &
+ flux = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
* block % mesh % dvEdge % array(iEdge) &
* h_edge(k,iEdge) &
* (hSum + sshEdge)
tend_ssh(cell1) = tend_ssh(cell1) - flux
tend_ssh(cell2) = tend_ssh(cell2) + flux
+
+ block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ = block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ + flux
end do
do iCell=1,block % mesh % nCells
@@ -710,16 +708,15 @@
! \zeta_{n+(j+1)/J} = \zeta_{n+j/J} + \frac{\Delta t}{J} \left( -
! </font>
<font color="gray">abla \cdot {\bf F}_j \right)
! sshSubcycleNew(iCell) = sshSubcycleOld(iCell) + dt/J*tend_ssh(iCell)/areaCell(iCell)
- block % state % time_levs(2) % state % sshSubcycle % array(iCell) &
- = block % state % time_levs(1) % state % sshSubcycle % array(iCell) &
+ 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)
-! sshNew(iCell) = sshNew(iCell) + aBtrSumCoef(j)*sshSubcycleNew(iCell)
+! sshNew(iCell) = sshNew(iCell) + sshSubcycleNew(iCell)
block % state % time_levs(2) % state % ssh % array(iCell) &
= block % state % time_levs(2) % state % ssh % array(iCell) &
- + aBtrSumCoef(j) &
- * block % state % time_levs(2) % state % sshSubcycle % array(iCell)
+ + block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell)
end do
@@ -733,7 +730,7 @@
do while (associated(block))
call dmpar_exch_halo_field1dReal(domain % dminfo, &
- block % state % time_levs(2) % state % sshSubcycle % array(:), &
+ block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
@@ -741,41 +738,58 @@
end do ! block
-! |compute \ubtr_{n+j/J}^{\perp} from \ubtr_{n+j/J}
-! \ubtr_{n+(j+1)/J} = \ubtr_{n+j/J} + \frac{\Delta t}{J} \left(
-! - f\ubtr_{n+j/J}^{\perp}
-! - \frac{g}{\rho_0}</font>
<font color="red">abla \zeta_{n+(j+1)/J} + {\overline {\bf G}}
-! \right),
+ ! compute new barotropic velocity
+
! do iEdge=1,nEdges
! cell1 = cellsOnEdge(1,iEdge)
! cell2 = cellsOnEdge(2,iEdge)
! grad_ssh = - gravity*rho0Inv*( sshSubcycleNew(cell2) &
! - sshSubcycleNew(cell1) )/dcEdge(iEdge)
-! uBtrSubCycleNew(iEdge) = uBtrSubcycleOld(iEdge) + dt*(fEdge(iEdge)*uBtrPerp(iEdge)
-! - grad_ssh + GBtrForcing(iEdge))
-! uBtrNew(iEdge) = uBtrNew(iEdge) + aBtrSumCoef(j)*uBtrSubcycleNew(iEdge)
+! uBtrNew(iEdge) = uBtrNew(iEdge) + uBtrSubcycleNew(iEdge)
! end do
block => domain % blocklist
do while (associated(block))
- allocate(tend_uBtr(block % mesh % nEdges))
- do iCell=1,block % mesh % nCells
+
+ do iEdge=1,block % mesh % nEdges
- block % state % time_levs(2) % state % sshSubcycle % array(iCell) &
- = block % state % time_levs(1) % state % sshSubcycle % array(iCell) &
- + dt/config_n_btr_subcycles &
- * tend_ssh(iCell) / block % mesh % areaCell % array (iCell)
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- block % state % time_levs(2) % state % ssh % array(iCell) &
- = block % state % time_levs(2) % state % ssh % array(iCell) &
- + aBtrSumCoef(j) &
- * block % state % time_levs(2) % state % sshSubcycle % array(iCell)
+ ! Compute -f*uPerp
+ uPerp = 0.0
+ do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
+ eoe = block % mesh % edgesOnEdge % array(i,iEdge)
+ uPerp = uPerp + block % mesh % weightsOnEdge % array(i,iEdge) &
+ * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(eoe) &
+ * block % mesh % fEdge % array(eoe)
+ end do
+ ! timestep in uBtrSubcycle:
+! \ubtr_{n+(j+1)/J} = \ubtr_{n+j/J} + \frac{\Delta t}{J} \left(
+! - f\ubtr_{n+j/J}^{\perp}
+! - \frac{g}{\rho_0}</font>
<font color="gray">abla \zeta_{n+(j+1)/J} + {\overline {\bf G}}
+! \right),
+! uBtrSubCycleNew(iEdge) = uBtrSubcycleOld(iEdge) + dt*(fEdge(iEdge)*uBtrPerp(iEdge)
+! - grad_ssh + GBtrForcing(iEdge))
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + dt/config_n_btr_subcycles *( &
+ uPerp &
+ + block % state % time_levs(2) % state % u % array (k,iEdge) &
+ - grav_rho0Inv &
+ *( block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ - block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
+ /block % mesh % dcEdge % array(iEdge) &
+ + block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
+
+ block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ + block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
end do
- deallocate(tend_uBtr)
block => block % next
end do ! block
@@ -785,16 +799,18 @@
do while (associated(block))
call dmpar_exch_halo_field1dReal(domain % dminfo, &
- block % state % time_levs(2) % state % uBtrSubcycle % array(:), &
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:), &
block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
block => block % next
end do ! block
- deallocate(aBtrSumCoef,bBtrSumCoef)
+ ! advance time pointers
+ oldBtrSubcycleTime = mod(oldBtrSubcycleTime,2)+1
+ newBtrSubcycleTime = mod(newBtrSubcycleTime,2)+1
- end do !j=1,2*config_n_btr_subcycles
+ end do ! barotropic subcycling on j
! Note: Normalize so that \sum_{j=0}^{2J} a_j=1 and \sum_{j=1}^{2J} b_j=1. Then the following
! three lines are already done, so that\\
@@ -808,6 +824,32 @@
block => domain % blocklist
do while (associated(block))
+
+
+
+
+
+
+ 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
+
+ block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ / (config_n_btr_subcycles + 1)
+ end do
+
+ 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)
+ end do
+
+
+
+
+
block => block % next
end do ! block
@@ -844,9 +886,9 @@
call compute_wTop(block % state % time_levs(2) % state, block % mesh)
-!loop on cells
-! usplit only:
+ if (trim(config_time_integration) == 'higdon_unsplit') then
call compute_tend_h(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ endif
call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
@@ -880,10 +922,14 @@
block => domain % blocklist
do while (associated(block))
- ! for unsplit only
+ if (trim(config_time_integration) == 'higdon_unsplit') then
+
block % state % time_levs(2) % state % h % array(:,:) &
= block % state % time_levs(1) % state % h % array(:,:) &
+ dt* block % tend % h % array(:,:)
+
+ endif
+
! printing:
nCells = block % mesh % nCells
nEdges = block % mesh % nEdges
@@ -964,14 +1010,18 @@
! so the following lines are
! u_{n+1} = uBtr_{n+1} + 2*uBcl_{n+1/2} - uBcl_n
- ! mrp 110531 change once I have pointers set up for subcycling
do iEdge=1,block % mesh % nEdges
block % state % time_levs(2) % state % u % array(:,iEdge) &
- = block % state % time_levs(2) % state % uBtrSubcycle % 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)
enddo
+ do iCell=1,block % mesh % nCells
+ block % state % time_levs(2) % state % ssh % array(iCell) &
+ = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell)
+ enddo
+
u => block % state % time_levs(2) % state % u % array
tracers => block % state % time_levs(2) % state % tracers % array
h => block % state % time_levs(2) % state % h % array
</font>
</pre>