<p><b>mpetersen@lanl.gov</b> 2011-05-16 10:09:37 -0600 (Mon, 16 May 2011)</p><p>Added new variables and flags for time splitting, and higden_timestep subroutine. This version works for higdon unsplit, no iteration on baroclinic stage, on x3 quad to 1000 steps (p29d).<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-13 21:22:31 UTC (rev 834)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/Registry 2011-05-16 16:09:37 UTC (rev 835)
@@ -175,6 +175,15 @@
var persistent real tend_tracer1 ( nVertLevels nCells Time ) 1 - tracer1 tend tracers testing
var persistent real tend_tracer2 ( nVertLevels nCells Time ) 1 - tracer2 tend tracers testing
+# state variables for Higdon timesplitting
+var persistent real uBtr ( nEdges Time ) 2 o uBtr state - -
+var persistent real eta ( nCells Time ) 2 o eta state - -
+var persistent real uBtrSubcycle ( nEdges Time ) 2 o uBtrSubcycle state - -
+var persistent real etaSubcycle ( nCells Time ) 2 o etaSubcycle state - -
+var persistent real FBtr ( nEdges Time ) 1 o FBtr state - -
+var persistent real GBtrForcing ( nEdges Time ) 1 o GBtrForcing state - -
+var persistent real uBcl ( nVertLevels nEdges Time ) 2 o uBcl state - -
+
# Diagnostic fields: only written to output
var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
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-05-13 21:22:31 UTC (rev 834)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_mpas_core.F 2011-05-16 16:09:37 UTC (rev 835)
@@ -99,12 +99,16 @@
block % state % time_levs(1) % state % u % array( &
block % mesh % maxLevelEdgeBot % array(iEdge)+1: &
- block % mesh % nVertLevels,iEdge) = -1e34
+ block % mesh % nVertLevels,iEdge) = 0.0
+! mrp changed to 0
+! block % mesh % nVertLevels,iEdge) = -1e34
end do
do iCell=1,block % mesh % nCells
block % state % time_levs(1) % state % tracers % array( &
:, block % mesh % maxLevelCell % array(iCell)+1 &
- :block % mesh % nVertLevels,iCell) = -1e34
+ :block % mesh % nVertLevels,iCell) = 0.0
+! mrp changed to 0
+! :block % mesh % nVertLevels,iCell) = -1e34
end do
if (.not. config_do_restart) then
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-13 21:22:31 UTC (rev 834)
+++ branches/ocean_projects/timesplitting_mrp/src/core_ocean/module_time_integration.F 2011-05-16 16:09:37 UTC (rev 835)
@@ -356,7 +356,7 @@
subroutine higdon_timestep(domain, dt)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Advance model state forward in time by the specified time step using
- ! 4th order Runge-Kutta
+ ! Higdon timestepping scheme
!
! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
! plus grid meta-data
@@ -378,25 +378,67 @@
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
- real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp
+ real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, G
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
+print *, '1'
+
!
- ! Initialize time_levs(2) with state at current time
- ! Initialize first RK state
- ! Couple tracers time_levs(2) with h in time-levels
- ! Initialize RK weights
+ ! Prep variables before first iteration
!
block => domain % blocklist
do while (associated(block))
- block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+! 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
+
+ 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) &
+ - 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)
+
+ block % state % time_levs(2) % state % uBcl % array(:,iEdge) &
+ = block % state % time_levs(1) % state % uBcl % array(:,iEdge)
+
+ enddo ! iEdge
+
+ block % state % time_levs(2) % state % h % array(:,:) &
+ = block % state % time_levs(1) % state % h % array(:,:)
+
+ block % state % time_levs(2) % state % h_edge % array(:,:) &
+ = block % state % time_levs(1) % state % h_edge % array(:,:)
+
do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % maxLevelCell % array(iCell)
- block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % h % array(k,iCell)
+ ! change to maxLevelCell % array(iCell) ?
+ do k=1,block % mesh % nVertLevels
+
+ ! mrp 110516 efficiency note: Can we get rid of all this mult and divide
+ ! of tracers by h here? Is it needed?
+ ! here the tracer array now has tracer*h.
+ block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ * block % state % time_levs(1) % state % h % array(k,iCell)
+
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ = block % state % time_levs(1) % state % tracers % array(:,k,iCell)
end do
end do
@@ -404,11 +446,12 @@
block => block % next
end do
+print *, '2'
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! BEGIN large iteration (predictor/corrector when config_n_ts_iter=2)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Begin iteration
+ !
do higdon_step = 1, config_n_ts_iter
! --- update halos for diagnostic variables
@@ -430,48 +473,132 @@
block => block % next
end do
+print *, '3'
!
! Stage 1: Baroclinic velocity (3D) prediction, explicit with long timestep
!
-! --- compute tendencies
+ ! compute velocity tendencies, T(u*,w*,p*)
+ block => domain % blocklist
+ do while (associated(block))
+ if (.not.config_implicit_vertical_mix) then
+ call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
+ end if
+ call compute_tend_u(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+! mrp soon: for split version: add g/rho grad eta to tend_u
+ call enforce_boundaryEdge(block % tend, block % mesh)
+ block => block % next
+ end do
+
+ ! baroclinic iterations on linear Coriolis term
+! change this from beg later
+ do j=1,config_n_bcl_iter_beg
+
+! mrp soon: add linear Coriolis to tend_u maybe at end of this iteration?
+ !compute {\bf u}_{k,n+1}^{'\;\perp} from {\bf u}_{k,n+1}'
+! call compute_uPerp(uBclNew,uBclPerp)
+! need subroutine and variable just for perp, say uBclPerp - may already be there.
+
+print *, '4'
+
+ block => domain % blocklist
+ do while (associated(block))
+ allocate(G(block % mesh % nVertLevels))
+
+ do iEdge=1,block % mesh % nEdges
+ 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 \eta^*
+
+! This soon, with coriolis:
+! G(k) = -fEdge(iEdge)*uBclPerp(k,iEdge) + tend_u(k,iEdge)
+ G(k) = block % tend % u % array (k,iEdge)
+
+! GBtrForcing(iEdge) = GBtrForcing(iEdge) + hOld(k,iEdge)* G(k))
+
+ 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
+
+ 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) &
+ + dt * (G(k) - block % state % time_levs(1) % state % GBtrForcing % array(iEdge))
+ enddo
+ enddo
+
+ deallocate(G)
+ block => block % next
+ end do
+
block => domain % blocklist
do while (associated(block))
- if (.not.config_implicit_vertical_mix) then
- call compute_vertical_mix_coefficients(block % state % time_levs(2) % state, block % diagnostics, block % mesh)
- end if
- call compute_tend_u(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
- call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
- call enforce_boundaryEdge(block % tend, block % mesh)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % uBcl % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
block => block % next
end do
+ enddo ! config_n_bcl_iter
+print *, '5'
+
+
+
!
! Stage 2: Barotropic velocity (2D) prediction, explicitly subcycled
!
+ ! for unsplit version only:
+ block => domain % blocklist
+ do while (associated(block))
+ ! for unsplit only:
+ block % state % time_levs(2) % state % uBtr % array(:) = 0.0
+
+ block => block % next
+ end do ! block
+
+
+print *, '6'
+
!
! Stage 3: Tracer, density, pressure, vertical velocity prediction
!
- block => domain % blocklist
- do while (associated(block))
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ 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
+
+ call compute_wTop(block % state % time_levs(2) % state, block % mesh)
+
+!loop on cells
+! usplit only:
+ call compute_tend_h(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+
call compute_scalar_tend(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
- call compute_tend_h(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
- call enforce_boundaryEdge(block % tend, block % mesh)
+
block => block % next
- end do
+ end do
+print *, '7'
! --- update halos for prognostic variables
block => domain % blocklist
do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
@@ -480,27 +607,46 @@
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
+print *, '8'
-!--- accumulate update (for HIGDON_TIMESTEP)
block => domain % blocklist
do while (associated(block))
-! block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
-! + rk_weights(rk_step) * block % tend % u % array(:,:)
-! block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
-! + rk_weights(rk_step) * block % tend % h % array(:,:)
+ ! for unsplit only
+ block % state % time_levs(2) % state % h % array(:,:) &
+ = block % state % time_levs(1) % state % h % array(:,:) &
+ + dt* block % tend % h % array(:,:)
+! printing:
+ nCells = block % mesh % nCells
+ nEdges = block % mesh % nEdges
+ nVertLevels = block % mesh % nVertLevels
+!print *, 'hnew ',minval(block % tend % h % array(:,1:nCells)),&
+! maxval(block % tend % h % array(:,1:nCells))
+!print *, 'hnew1 ',minval(block % tend % h % array(1,1:nCells)),&
+! maxval(block % tend % h % array(1,1:nCells))
+! do iCell = 1,10 !nCells
+! print '(a,i5,50e12.2)', 'htend iCell=',iCell,block % tend % h % array(:,iCell)
+! enddo
+! printing end
+
+ ! mrp 110512 at a later time, only need T & S for iterations
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
- ! block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
- ! block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- ! + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ = (block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ + dt * block % tend % tracers % array(:,k,iCell) &
+ ) / block % state % time_levs(2) % state % h % array(k,iCell)
end do
end do
- block => block % next
- end do
+ ! 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
+
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END large iteration loop
@@ -511,6 +657,33 @@
!
block => domain % blocklist
do while (associated(block))
+print *, '9'
+! 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 *, 'hold1 ',minval(block % state % time_levs(1) % state % h % array(1,1:nCells)),&
+ maxval(block % state % time_levs(1) % state % h % array(1,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))
+print *, 'unew ',minval(block % state % time_levs(2) % state % u % array(:,1:nEdges)),&
+ maxval(block % state % time_levs(2) % state % u % array(:,1:nEdges))
+print *, 'hnew ',minval(block % state % time_levs(2) % state % h % array(:,1:nCells)),&
+ maxval(block % state % time_levs(2) % state % h % array(:,1:nCells))
+print *, 'hnew1 ',minval(block % state % time_levs(2) % state % h % array(1,1:nCells)),&
+ maxval(block % state % time_levs(2) % state % h % array(1,1:nCells))
+print *, 'Tnew ',minval(block % state % time_levs(2) % state % tracers % array(1,:,1:nCells)),&
+ maxval(block % state % time_levs(2) % state % tracers % array(1,:,1:nCells))
+print *, 'Snew ',minval(block % state % time_levs(2) % state % tracers % array(2,:,1:nCells)),&
+ maxval(block % state % time_levs(2) % state % tracers % array(2,:,1:nCells))
+ block % state % time_levs(2) % state % tracers % array(3,1,1:nCells) = domain % dminfo % my_proc_id
+! printing end
u => block % state % time_levs(2) % state % u % array
tracers => block % state % time_levs(2) % state % tracers % array
@@ -526,11 +699,12 @@
nEdges = block % mesh % nEdges
nVertLevels = block % mesh % nVertLevels
- do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- tracers(:,k,iCell) = tracers(:,k,iCell) / h(k,iCell)
- end do
- end do
+! dividing by h above for higdon.
+! do iCell=1,nCells
+! do k=1,maxLevelCell(iCell)
+! tracers(:,k,iCell) = tracers(:,k,iCell) / h(k,iCell)
+! end do
+! end do
if (config_implicit_vertical_mix) then
allocate(A(nVertLevels),C(nVertLevels),uTemp(nVertLevels), &
@@ -1821,6 +1995,9 @@
! mrp 101115 note: in order to include flux boundary conditions, we will need to
! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
+ ! mrp 110516 efficiency note: For z-level, only do this on level 1. h_edge for all
+ ! lower levels is defined by hZlevel.
+
coef_3rd_order = 0.
if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
@@ -1927,6 +2104,7 @@
! set the velocity and height at dummy address
! used -1e34 so error clearly occurs if these values are used.
!
+!mrp 110516 change to zero, change back later:
u(:,nEdges+1) = -1e34
h(:,nCells+1) = -1e34
tracers(s % index_temperature,:,nCells+1) = -1e34
@@ -2165,6 +2343,61 @@
endif
+ call compute_wtop(s,grid)
+
+
+ end subroutine compute_solve_diagnostics
+
+
+ subroutine compute_wtop(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute diagnostic fields used in the tendency computations
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ ! mrp 110512 could clean this out, remove pointers?
+ integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
+ real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv
+
+ integer :: nCells, nEdges, nVertices, nVertLevels, vertexDegree
+
+
+ real (kind=RKIND), dimension(:), pointer :: &
+ h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ hZLevel
+ real (kind=RKIND), dimension(:,:), pointer :: u,wTop
+ real (kind=RKIND), dimension(:,:), allocatable:: div_u
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, &
+ verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, &
+ boundaryEdge, boundaryCell
+ integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
+ maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
+ maxLevelVertexBot, maxLevelVertexTop
+
+
+ u => s % u % array
+ wTop => s % wTop % array
+
+ areaCell => grid % areaCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ hZLevel => grid % hZLevel % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeBot => grid % maxLevelEdgeBot % array
+ dvEdge => grid % dvEdge % array
+
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ nVertLevels = grid % nVertLevels
+
!
! vertical velocity through layer interface
!
@@ -2197,14 +2430,14 @@
wTop(maxLevelCell(iCell)+1,iCell) = 0.0
do k=maxLevelCell(iCell),2,-1
wTop(k,iCell) = wTop(k+1,iCell) &
- - div_u(k,iCell)/areaCell(iCell)*h(k,iCell)
+ - div_u(k,iCell)/areaCell(iCell)*hZLevel(k)
end do
end do
deallocate(div_u)
endif
- end subroutine compute_solve_diagnostics
+ end subroutine compute_wtop
subroutine compute_vertical_mix_coefficients(s, d, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
</font>
</pre>