<p><b>mpetersen@lanl.gov</b> 2011-05-19 16:39:34 -0600 (Thu, 19 May 2011)</p><p>BRANCH COMMIT: Added half of the new code for the Higdon split method. Code compiles and works fine for config_time_integration = 'higdon_unsplit'. Does not yet work for 'higdon_split'.<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-19 21:41:14 UTC (rev 847)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-19 22:39:34 UTC (rev 848)
@@ -371,7 +371,7 @@
integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step, split
type (block_type), pointer :: block
- real (kind=RKIND) :: uSum, hSum, GSum, grav_rho0Inv
+ real (kind=RKIND) :: uSum, hSum, GSum, grav_rho0Inv, sshEdge, flux
integer :: nCells, nEdges, nVertLevels, num_tracers
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -379,7 +379,7 @@
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
- real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, G
+ real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, G, aBtrSumCoef,bBtrSumCoef, tend_ssh, tend_uBtr
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
!print *, '1'
@@ -609,17 +609,200 @@
! Stage 2: Barotropic velocity (2D) prediction, explicitly subcycled
!
- ! for unsplit version only:
- block => domain % blocklist
- do while (associated(block))
+ 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
+ block => block % next
+ end do ! block
- ! for unsplit only:
- block % state % time_levs(2) % state % uBtr % array(:) = 0.0
+ elseif (trim(config_time_integration) == 'higdon_split') then
- block => block % next
- end do ! block
+ block => domain % blocklist
+ do while (associated(block))
+
+ ! 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
+
+ do iCell=1,block % mesh % nCells
+! sshSubcycleOld = sshOld
+ block % state % time_levs(1) % state % sshSubcycle % array(iCell) &
+ = block % state % time_levs(1) % state % ssh % array(iCell)
+
+! sshNew = aBtrSumCoef(0)*sshOld
+ block % state % time_levs(2) % state % ssh % array(iCell) &
+ = aBtrSumCoef(0) &
+ * 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) &
+ = 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)
+
+! 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
+
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ 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)
+
+! flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+! FBtr(iEdge) = FBtr(iEdge) + bBtrSumCoef(j)*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) )
+ 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) &
+ * 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
+ end do
+ do iCell=1,nCells
+ tend_ssh(iCell) = tend_ssh(iCell) / block % mesh % areaCell % array (iCell)
+ end do
+ do iCell=1,block % mesh % nCells ! couple tracers to h
+
+! \zeta_{n+(j+1)/J} = \zeta_{n+j/J} + \frac{\Delta t}{J} \left( -
+! </font>
<font color="blue">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) &
+ + dt/config_n_btr_subcycles &
+ * tend_ssh(iCell) / block % mesh % areaCell % array (iCell)
+
+! sshNew(iCell) = sshNew(iCell) + aBtrSumCoef(j)*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)
+
+ end do
+
+ deallocate(tend_ssh)
+ block => block % next
+ end do ! block
+
+
+! |boundary update on \zeta_{n+(j+1)/J}
+ block => domain % blocklist
+ do while (associated(block))
+
+ call dmpar_exch_halo_field1dReal(domain % dminfo, &
+ block % state % time_levs(2) % state % sshSubcycle % array(:), &
+ block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+
+ block => block % next
+ 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="blue">abla \zeta_{n+(j+1)/J} + {\overline {\bf G}}
+! \right),
+! 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)
+! end do
+ block => domain % blocklist
+ do while (associated(block))
+ allocate(tend_uBtr(block % mesh % nEdges))
+
+ deallocate(tend_uBtr)
+ block => block % next
+ end do ! block
+
+
+! |boundary update on \ubtr_{n+(j+1)/J}
+ block => domain % blocklist
+ do while (associated(block))
+
+ call dmpar_exch_halo_field1dReal(domain % dminfo, &
+ block % state % time_levs(2) % state % uBtrSubcycle % array(:), &
+ block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ block => block % next
+ end do ! block
+
+ deallocate(aBtrSumCoef,bBtrSumCoef)
+
+ end do !j=1,2*config_n_btr_subcycles
+
+! 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\\
+! \zeta^* = \left. \sum_{j=0}^{2J} a_j \zeta_{n+j/J} \right/ \sum_{j=0}^{2J} a_j is ! sshNew
+! \ubtr^* = \left. \sum_{j=0}^{2J} a_j \ubtr_{n+j/J} \right/ \sum_{j=0}^{2J} a_j is ! uBtrNew
+! {\overline {\bf F}} = \left. \sum_{j=1}^{2J} b_j {\bf F}_j \right/ \sum_{j=1}^{2J} b_j is ! FBtr
+! {\bf u}^*_k = \ubtr^* + {\bf u}^{'*}_k
+! uNew(iEdge) = uBtrNew(iEdge) + 0.5*(uBclOld(iEdge) + uBclNew(iEdge))
+! h_1^* = \Delta z_1 + \zeta^*
+! hNew(iCell) = hZlevel(1) + sshNew(iCell)
+ block => domain % blocklist
+ do while (associated(block))
+
+ block => block % next
+ end do ! block
+
+
+! boundary update on {\overline {\bf F}}
+ block => domain % blocklist
+ do while (associated(block))
+
+ block => block % next
+ end do ! block
+
+
+ endif
+
+
+
!print *, '6'
!
</font>
</pre>