<p><b>mpetersen@lanl.gov</b> 2011-05-19 07:32:11 -0600 (Thu, 19 May 2011)</p><p>BRANCH COMMIT: Begin to add things for split version. Added initial barotropic/baroclinic splitting of u, splitting of Gbar in stage 1, and addition of g/rho0 grad(h_1^*) in stage 1.<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 00:34:30 UTC (rev 844)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-19 13:32:11 UTC (rev 845)
@@ -369,8 +369,9 @@
type (domain_type), intent(inout) :: domain
real (kind=RKIND), intent(in) :: dt
- integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step
+ integer :: iCell, i,k,j, iEdge, cell1, cell2, higdon_step, split
type (block_type), pointer :: block
+ real (kind=RKIND) :: uSum, hSum, GSum, grav_rho0Inv
integer :: nCells, nEdges, nVertLevels, num_tracers
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -390,23 +391,27 @@
block => domain % blocklist
do while (associated(block))
-! printing:
-! nCells = block % mesh % nCells
-! nEdges = block % mesh % nEdges
-! nVertLevels = block % mesh % nVertLevels
-!print *, 'uold',minval(block % state % time_levs(1) % state % u % array(:,1:nEdges)),&
-! maxval(block % state % time_levs(1) % state % u % array(:,1:nEdges))
-!print *, 'hold',minval(block % state % time_levs(1) % state % h % array(:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % h % array(:,1:nCells))
-!print *, 'Told',minval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells))
-!print *, 'Sold',minval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells)),&
-! maxval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells))
-! printing end
+ ! Compute barotropic velocity at the old timestep
+ if (trim(config_time_integration) == 'higdon_unsplit') then
+ block % state % time_levs(1) % state % uBtr % array(:) = 0.0
+ elseif (trim(config_time_integration) == 'higdon_split') then
+ do iEdge=1,block % mesh % nEdges
+ uSum = 0.0
+ hSum = 0.0
+ ! 110518 mrp efficiency note:
+ ! in a z-level model, could change hSum to use sum of hZLevel, or have another
+ ! 1-D z array which is total thickness between the surface and that level.
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ uSum = uSum + block % state % time_levs(1) % state % u % array(k,iEdge)
+ hSum = hSum + block % state % time_levs(1) % state % h_edge % array(k,iEdge)
+ enddo
+ block % state % time_levs(1) % state % uBtr % array(iEdge) = uSum/hSum
+ enddo
+ endif
+!mrp 110512 need to add w, p, and ssh here?
+
do iEdge=1,block % mesh % nEdges
- ! for unsplit only:
- block % state % time_levs(1) % state % uBtr % array(iEdge) = 0.0
block % state % time_levs(1) % state % uBcl % array(:,iEdge) &
= block % state % time_levs(1) % state % u % array(:,iEdge) &
@@ -435,7 +440,23 @@
end do
end do
-!mrp 110512 need to add w, p, and ssh here?
+! printing:
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertLevels = block % mesh % nVertLevels
+print *, 'uold ',minval(block % state % time_levs(1) % state % u % array(:,1:nEdges)),&
+ maxval(block % state % time_levs(1) % state % u % array(:,1:nEdges))
+print *, 'uBtrOld',minval(block % state % time_levs(1) % state % uBtr % array(1:nEdges)),&
+ maxval(block % state % time_levs(1) % state % uBtr % array(1:nEdges))
+print *, 'uBclOld',minval(block % state % time_levs(1) % state % uBcl % array(:,1:nEdges)),&
+ maxval(block % state % time_levs(1) % state % uBcl % array(:,1:nEdges))
+!print *, 'hold',minval(block % state % time_levs(1) % state % h % array(:,1:nCells)),&
+! maxval(block % state % time_levs(1) % state % h % array(:,1:nCells))
+!print *, 'Told',minval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells)),&
+! maxval(block % state % time_levs(1) % state % tracers % array(1,:,1:nCells))
+!print *, 'Sold',minval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells)),&
+! maxval(block % state % time_levs(1) % state % tracers % array(2,:,1:nCells))
+! printing end
block => block % next
end do
@@ -496,6 +517,14 @@
!print *, '4'
+ ! Use this G coefficient to avoid an if statement within the iEdge loop.
+ if (trim(config_time_integration) == 'higdon_unsplit') then
+ split = 0
+ elseif (trim(config_time_integration) == 'higdon_split') then
+ split = 1
+ endif
+ grav_rho0Inv = gravity/config_rho0
+
block => domain % blocklist
do while (associated(block))
allocate(G(block % mesh % nVertLevels))
@@ -504,16 +533,24 @@
call compute_f_uperp(block % state % time_levs(2) % state , block % mesh)
do iEdge=1,block % mesh % nEdges
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+
G = 0.0 ! could put this after with G(maxleveledgetop+1:nvertlevels)=0
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
!{\bf G}_k = -f{\bf u}_{k,n+1}^{'\;\perp} + {\bf T}({\bf u}_k^*,w_k^*,p_k^*) +\frac{g}{\rho_0}</font>
<font color="gray">abla \ssh^*
! This soon, with coriolis:
-! G(k) = -fEdge(iEdge)*uBclPerp(k,iEdge) + tend_u(k,iEdge)
+! G(k) = tend_u(k,iEdge) -fEdge(iEdge)*uBclPerp(k,iEdge) - g/rho0 grad (h_1^*)
! Here uNew is a work variable containing -fEdge(iEdge)*uBclPerp(k,iEdge)
! mrp temp new one:
- G(k) = block % state % time_levs(2) % state % u % array (k,iEdge) &
- + block % tend % u % array (k,iEdge)
+ G(k) = 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 % mesh % dcEdge % array(iEdge)
+
! mrp old one:
! G(k) = block % tend % u % array (k,iEdge)
@@ -521,13 +558,21 @@
enddo
-! set GBtrForcing to zero for unsplit.
- !{\overline {\bf G}} = \left. \sum_{k=1}^{N^{edge}} h_{k,n}^{edge} {\bf G}_k \right/\sum_{k=1}^{N^{edge}} h_{k,n}^{edge}
-! GBtrForcing(iEdge) = GBtrForcing(iEdge) /(hOld(1,iEdge) + h2toNZLevel)
-! mrp this is for unsplit only:
- block % state % time_levs(1) % state % GBtrForcing % array(iEdge) = 0.0
+ GSum = 0.0
+ hSum = 0.0
+ ! 110518 mrp efficiency note:
+ ! in a z-level model, could change hSum to use sum of hZLevel, or have another
+ ! 1-D z array which is total thickness between the surface and that level.
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ GSum = GSum + G(k)
+ hSum = hSum + block % state % time_levs(1) % state % h_edge % array(k,iEdge)
+ enddo
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge) = split*GSum/hSum
+
+
+
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
!{\bf u}'_{k,n+1} = {\bf u}'_{k,n} + \Delta t \left({\bf G}_k - {\overline {\bf G}}\right)
block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
= block % state % time_levs(1) % state % uBcl % array(k,iEdge) &
</font>
</pre>