<p><b>mpetersen@lanl.gov</b> 2012-02-01 11:47:09 -0700 (Wed, 01 Feb 2012)</p><p>BRANCH COMMIT<br>
First cut at revising split explicit to use ale coordinate. This revision has:<br>
<br>
before: In stage 3, h is computed from barotropic SSH computed in<br>
stage 2. This can only be done with z-level coordinates, and can not<br>
be generalized to ALE.<br>
<br>
now: Whereever h is needed, use h^*, known from the previous step or<br>
iteration.<br>
<br>
This revision includes all deleted lines simply commented out. It<br>
appears to work for single processor, but has some boundary update<br>
bug.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_tendency.F 2012-01-31 21:12:35 UTC (rev 1443)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_tendency.F 2012-02-01 18:47:09 UTC (rev 1444)
@@ -137,8 +137,6 @@
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
! for explanation of divergence operator.
!
- ! for z-level, only compute height tendency for top layer.
-
call mpas_timer_start("hadv", .false., thickHadvTimer)
call ocn_thick_hadv_tend(grid, u, h_edge, tend_h, err)
call mpas_timer_stop("hadv", thickHadvTimer)
@@ -146,7 +144,6 @@
!
! height tendency: vertical advection term -d/dz(hw)
!
- ! Vertical advection computed for top layer of a z grid only.
call mpas_timer_start("vadv", .false., thickVadvTimer)
call ocn_thick_vadv_tend(grid, wtop, tend_h, err)
call mpas_timer_stop("vadv", thickVadvTimer)
@@ -440,7 +437,7 @@
real (kind=RKIND), dimension(:), pointer :: &
h_s, fVertex, dvEdge, dcEdge, areaCell, areaTriangle, &
- zTopZLevel
+ zTopZLevel, ssh
real (kind=RKIND), dimension(:,:), pointer :: &
weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure,&
circulation, vorticity, ke, ke_edge, MontPot, wTop, zMid, &
@@ -470,6 +467,7 @@
MontPot => s % MontPot % array
pressure => s % pressure % array
zMid => s % zMid % array
+ ssh => s % ssh % array
tracers => s % tracers % array
weightsOnEdge => grid % weightsOnEdge % array
@@ -843,9 +841,19 @@
endif
-! mrp 111206 remove
-! call ocn_wtop(s,grid)
+ !
+ ! Sea Surface Height
+ !
+ do iCell=1,nCells
+ ! Start at the bottom where we know the depth, and go up.
+ ! The bottom depth for this cell is
+ ! zTopZLevel(maxLevelCell(iCell)+1)
+ ssh(iCell) = zTopZLevel(maxLevelCell(iCell)+1) &
+ + sum(h(1:maxLevelCell(iCell),iCell))
+
+ end do
+
end subroutine ocn_diagnostic_solve!}}}
!***********************************************************************
@@ -861,7 +869,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_wtop(s, grid)!{{{
+ subroutine ocn_wtop(s1,s2, grid)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
!
@@ -872,7 +880,8 @@
implicit none
- type (state_type), intent(inout) :: s
+ type (state_type), intent(in) :: s1
+ type (state_type), intent(inout) :: s2
type (mesh_type), intent(in) :: grid
! mrp 110512 could clean this out, remove pointers?
@@ -895,10 +904,10 @@
maxLevelCell, maxLevelEdgeTop, maxLevelEdgeBot, &
maxLevelVertexBot, maxLevelVertexTop
- u => s % u % array
- h => s % h % array
- wTop => s % wTop % array
- h_edge => s % h_edge % array
+ h => s1 % h % array
+ h_edge => s1 % h_edge % array
+ u => s2 % u % array
+ wTop => s2 % wTop % array
areaCell => grid % areaCell % array
cellsOnEdge => grid % cellsOnEdge % array
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-01-31 21:12:35 UTC (rev 1443)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-01 18:47:09 UTC (rev 1444)
@@ -168,7 +168,7 @@
do while (associated(block))
! mrp 111206 put ocn_wtop call at top for ALE
- call ocn_wtop(provis, block % mesh)
+ call ocn_wtop(provis, provis, block % mesh)
if (.not.config_implicit_vertical_mix) then
call ocn_vmix_coefs(block % mesh, provis, block % diagnostics, err)
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-01-31 21:12:35 UTC (rev 1443)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-01 18:47:09 UTC (rev 1444)
@@ -87,12 +87,15 @@
integer :: iCell, i,k,j, iEdge, cell1, cell2, split_explicit_step, split, &
eoe, oldBtrSubcycleTime, newBtrSubcycleTime, uPerpTime, BtrCorIter, &
n_bcl_iter(config_n_ts_iter), &
- vertex1, vertex2, iVertex, sshSwapFlag
+ vertex1, vertex2, iVertex
+! mrp 120131 remove h computed from SSH in current step
+!, sshSwapFlag
type (block_type), pointer :: block
real (kind=RKIND) :: uhSum, hSum, sshEdge, flux, &
- uPerp, uCorr, tracerTemp, coef, FBtr_coeff, sshCell1, sshCell2
- real (kind=RKIND), dimension(:), pointer :: sshNew
+ uPerp, uCorr, temp, coef, FBtr_coeff, sshCell1, sshCell2
+! mrp 120131 remove h computed from SSH in current step
+! real (kind=RKIND), dimension(:), pointer :: sshNew
integer :: num_tracers, ucorr_coef, err
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -100,10 +103,11 @@
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
integer, dimension(:), pointer :: &
maxLevelCell, maxLevelEdgeTop
- real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp, hNew
+ real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
- sshSwapFlag = 1
+! mrp 120131 remove h computed from SSH in current step
+! sshSwapFlag = 1
call mpas_timer_start("se timestep", .false., timer_main)
@@ -133,10 +137,13 @@
enddo ! iEdge
- ! Initialize * variables that are used compute baroclinic tendencies below.
- block % state % time_levs(2) % state % ssh % array(:) &
+ ! Initialize * variables that are used to compute baroclinic tendencies below.
+ block % state % time_levs(2) % state % ssh % array(:) &
= block % state % time_levs(1) % state % ssh % array(:)
+ 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(:,:)
@@ -239,16 +246,24 @@
/block % mesh % dcEdge % array(iEdge) )
enddo
- ! Compute GBtrForcing, the vertically averaged forcing
- sshEdge = 0.5*( block % state % time_levs(1) % state % ssh % array(cell1) &
- + block % state % time_levs(1) % state % ssh % array(cell2) )
+! mrp 120131 remove h computed from SSH in current step
+! ! Compute GBtrForcing, the vertically averaged forcing
+! sshEdge = 0.5*( block % state % time_levs(1) % state % ssh % array(cell1) &
+! + block % state % time_levs(1) % state % ssh % array(cell2) )
+! uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * uTemp(1)
+! hSum = sshEdge + block % mesh % hZLevel % array(1)
- uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * uTemp(1)
- hSum = sshEdge + block % mesh % hZLevel % array(1)
+ uhSum = block % state % time_levs(2) % state % h_edge % array(1,iCell) * uTemp(1)
+ hSum = block % state % time_levs(2) % state % h_edge % array(1,iCell)
do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + block % mesh % hZLevel % array(k) * uTemp(k)
- hSum = hSum + block % mesh % hZLevel % array(k)
+! mrp 120131 remove h computed from SSH in current step
+! uhSum = uhSum + block % mesh % hZLevel % array(k) * uTemp(k)
+! hSum = hSum + block % mesh % hZLevel % array(k)
+
+ uhSum = uhSum + block % state % time_levs(2) % state % h_edge % array(k,iCell) * uTemp(k)
+ hSum = hSum + block % state % time_levs(2) % state % h_edge % array(k,iCell)
+
enddo
block % state % time_levs(1) % state % GBtrForcing % array(iEdge) = split*uhSum/hSum/dt
@@ -329,8 +344,9 @@
! sshSubcycleOld = sshOld
block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) = block % state % time_levs(1) % state % ssh % array(iCell)
- ! sshNew = sshOld This is the first for the summation
- block % state % time_levs(2) % state % ssh % array(iCell) = block % state % time_levs(1) % state % ssh % array(iCell)
+! mrp 120131 remove h computed from SSH in current step
+! ! sshNew = sshOld This is the first for the summation
+! block % state % time_levs(2) % state % ssh % array(iCell) = block % state % time_levs(1) % state % ssh % array(iCell)
enddo
do iEdge=1,block % mesh % nEdges
@@ -425,6 +441,7 @@
! config_btr_gam1_uWt1= 1 flux = uBtrNew*H
! config_btr_gam1_uWt1=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam1_uWt1= 0 flux = uBtrOld*H
+ ! mrp 120201 efficiency: could we combint the following edge and cell loops?
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -533,7 +550,7 @@
! config_btr_gam3_uWt2= 1 flux = uBtrNew*H
! config_btr_gam3_uWt2=0.5 flux = 1/2*(uBtrNew+uBtrOld)*H
! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
-
+ ! mrp 120201 efficiency: could we combint the following edge and cell loops?
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -582,11 +599,12 @@
block => domain % blocklist
do while (associated(block))
- ! Accumulate SSH in running sum over the subcycles.
- do iCell=1,block % mesh % nCells
- block % state % time_levs(2) % state % ssh % array(iCell) = block % state % time_levs(2) % state % ssh % array(iCell) &
- + block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell)
- end do
+! mrp 120131 remove h computed from SSH in current step
+! ! Accumulate SSH in running sum over the subcycles.
+! do iCell=1,block % mesh % nCells
+! block % state % time_levs(2) % state % ssh % array(iCell) = block % state % time_levs(2) % state % ssh % array(iCell) &
+! + block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell)
+! end do
! uBtrNew = uBtrNew + uBtrSubcycleNEW
! This accumulates the sum.
@@ -623,17 +641,18 @@
/ (config_n_btr_subcycles*config_btr_subcycle_loop_factor + 1)
end do
- if (config_SSH_from=='avg_of_SSH_subcycles') then
- 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*config_btr_subcycle_loop_factor + 1)
- end do
- elseif (config_SSH_from=='avg_flux') then
- ! see below
- else
- write(0,*) 'Abort: Unknown config_SSH_from option: ' //trim(config_SSH_from)
- call mpas_dmpar_abort(dminfo)
- endif
+! mrp 120131 remove h computed from SSH in current step
+! if (config_SSH_from=='avg_of_SSH_subcycles') then
+! 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*config_btr_subcycle_loop_factor + 1)
+! end do
+! elseif (config_SSH_from=='avg_flux') then
+! ! see below
+! else
+! write(0,*) 'Abort: Unknown config_SSH_from option: ' //trim(config_SSH_from)
+! call mpas_dmpar_abort(dminfo)
+! endif
block => block % next
end do ! block
@@ -661,37 +680,37 @@
allocate(uTemp(block % mesh % nVertLevels))
- if (config_SSH_from=='avg_flux') then
- ! Accumulate fluxes in the tend % ssh variable
- block % tend % ssh % array(:) = 0.0
- do iEdge=1,block % mesh % nEdges
- cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
- cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+!!!!!!!!!!!!!!!! do not delete this completely yet.
+!!!!!!!!! use this as a check to see if ssh from diag matches this from flux.
+! would need to be put at the end of the loop, after the ocn diag call.
+! if (config_SSH_from=='avg_flux') then
+! ! Accumulate fluxes in the tend % ssh variable
+! block % tend % ssh % array(:) = 0.0
+! do iEdge=1,block % mesh % nEdges
+! cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+! cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+! block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) &
+! - block % state % time_levs(1) % state % FBtr % array(iEdge) &
+! * block % mesh % dvEdge % array(iEdge)
+! block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) &
+! + block % state % time_levs(1) % state % FBtr % array(iEdge) &
+! * block % mesh % dvEdge % array(iEdge)
+! end do
+! do iCell=1,block % mesh % nCells
+! ! SSHnew = SSHold + dt*(-div(Flux))
+! block % state % time_levs(2) % state % ssh % array(iCell) = block % state % time_levs(1) % state % ssh % array(iCell) &
+! + dt * block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
+! end do
+!!!!!!!!!!!!!!!! END do not delete this completely yet.
- block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) &
- - block % state % time_levs(1) % state % FBtr % array(iEdge) &
- * block % mesh % dvEdge % array(iEdge)
+! mrp 120131 remove h computed from SSH in current step
+! endif
-
- block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) &
- + block % state % time_levs(1) % state % FBtr % array(iEdge) &
- * block % mesh % dvEdge % array(iEdge)
-
- end do
-
- do iCell=1,block % mesh % nCells
-
- ! SSHnew = SSHold + dt*(-div(Flux))
- block % state % time_levs(2) % state % ssh % array(iCell) = block % state % time_levs(1) % state % ssh % array(iCell) &
- + dt * block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
- end do
- endif
-
! Correction velocity uCorr = (Flux - Sum(h u*))/H
! or, for the full latex version:
- !u^{corr} = \left( {\overline {\bf F}}
- ! - \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right) u_k^* \right)
- !\left/ \sum_{k=1}^{N^{edge}} \left(\zeta_{k,n}^{*\;edge}+\Delta z_k\right) \right.
+ !{\bf u}^{corr} = \left( {\overline {\bf F}}
+ ! - \sum_{k=1}^{N^{edge}} h_{k,*}^{edge} {\bf u}_k^{avg} \right)
+ ! \left/ \sum_{k=1}^{N^{edge}} h_{k,*}^{edge} \right.
if (config_u_correction) then
ucorr_coef = 1
@@ -700,22 +719,29 @@
endif
do iEdge=1,block % mesh % nEdges
- cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
- cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+! mrp 120131 remove h computed from SSH in current step
+! cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+! cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- sshEdge = 0.5 * (block % state % time_levs(2) % state % ssh % array(cell1) &
- + block % state % time_levs(2) % state % ssh % array(cell2))
+! mrp 120131 remove h computed from SSH in current step
+! sshEdge = 0.5 * (block % state % time_levs(2) % state % ssh % array(cell1) &
+! + block % state % time_levs(2) % state % ssh % array(cell2))
- ! This is u*
+ ! This is u^{avg}
uTemp(:) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ block % state % time_levs(2) % state % uBcl % array(:,iEdge)
- uhSum = (sshEdge + block % mesh % hZLevel % array(1)) * uTemp(1)
- hSum = sshEdge + block % mesh % hZLevel % array(1)
+ uhSum = block % state % time_levs(2) % state % h_edge % array(1,iCell) * uTemp(1)
+ hSum = block % state % time_levs(2) % state % h_edge % array(1,iCell)
do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + block % mesh % hZLevel % array(k) * uTemp(k)
- hSum = hSum + block % mesh % hZLevel % array(k)
+! mrp 120131 remove h computed from SSH in current step
+! uhSum = uhSum + block % mesh % hZLevel % array(k) * uTemp(k)
+! hSum = hSum + block % mesh % hZLevel % array(k)
+
+ uhSum = uhSum + block % state % time_levs(2) % state % h_edge % array(k,iCell) * uTemp(k)
+ hSum = hSum + block % state % time_levs(2) % state % h_edge % array(k,iCell)
+
enddo
uCorr = ucorr_coef*(( block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum)
@@ -733,32 +759,34 @@
enddo
endif
- ! Put new sshEdge values in h_edge array, for the OcnTendScalar call below.
- block % state % time_levs(2) % state % h_edge % array(1,iEdge) = sshEdge + block % mesh % hZLevel % array(1)
+! mrp 120131 remove h computed from SSH in current step
+! ! Put new sshEdge values in h_edge array, for the OcnTendScalar call below.
+! block % state % time_levs(2) % state % h_edge % array(1,iEdge) = sshEdge + block % mesh % hZLevel % array(1)
+! do k=2,block % mesh % nVertLevels
+! block % state % time_levs(2) % state % h_edge % array(k,iEdge) = block % mesh % hZLevel % array(k)
+! enddo
- do k=2,block % mesh % nVertLevels
- block % state % time_levs(2) % state % h_edge % array(k,iEdge) = block % mesh % hZLevel % array(k)
- enddo
end do ! iEdge
- ! Put new SSH values in h array, for the OcnTendScalar call below.
- do iCell=1,block % mesh % nCells
- block % state % time_levs(2) % state % h % array(1,iCell) = block % state % time_levs(2) % state % ssh % array(iCell) &
- + block % mesh % hZLevel % array(1)
+! mrp 120131 remove h computed from SSH in current step
+! ! Put new SSH values in h array, for the OcnTendScalar call below.
+! do iCell=1,block % mesh % nCells
+! block % state % time_levs(2) % state % h % array(1,iCell) = block % state % time_levs(2) % state % ssh % array(iCell) &
+! + block % mesh % hZLevel % array(1)
+! ! mrp 110601 efficiency note: Since h just moves back and forth between pointers,
+! ! this is not necessary once initialized.
+! do k=2,block % mesh % nVertLevels * sshSwapFlag
+! block % state % time_levs(2) % state % h % array(k,iCell) = block % mesh % hZLevel % array(k)
+! enddo
+! enddo ! iCell
- ! mrp 110601 efficiency note: Since h just moves back and forth between pointers,
- ! this is not necessary once initialized.
- do k=2,block % mesh % nVertLevels * sshSwapFlag
- block % state % time_levs(2) % state % h % array(k,iCell) = block % mesh % hZLevel % array(k)
- enddo
- enddo ! iCell
-
deallocate(uTemp)
block => block % next
end do ! block
- sshSwapFlag = 0
- endif ! split_explicit
+! mrp 120131 remove h computed from SSH in current step
+! sshSwapFlag = 0
+ endif ! split_explicit
call mpas_timer_stop("se btr vel", timer_btr_vel)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -769,91 +797,139 @@
block => domain % blocklist
do while (associated(block))
- call ocn_wtop(block % state % time_levs(2) % state, block % mesh)
+ call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
- if (trim(config_time_integration) == 'unsplit_explicit') then
- call ocn_tend_h(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
- endif
-
+! mrp 120131 remove h computed from SSH in current step
+! if (trim(config_time_integration) == 'unsplit_explicit') then
+! call ocn_tend_h(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+! endif
+ call ocn_tend_h (block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
block => block % next
end do
- ! update halo for thicknes for unsplit only
- if (trim(config_time_integration) == 'unsplit_explicit') then
- call mpas_timer_start("se halo h", .false., timer_halo_h)
- block => domain % blocklist
- do while (associated(block))
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
- call mpas_timer_stop("se halo h", timer_halo_h)
- endif ! unsplit_explicit
+! mrp 120131 remove h computed from SSH in current step
+! ! update halo for thickness for unsplit only
+! if (trim(config_time_integration) == 'unsplit_explicit') then
+! call mpas_timer_start("se halo h", .false., timer_halo_h)
+! block => domain % blocklist
+! do while (associated(block))
+! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % h % array(:,:), &
+! block % mesh % nVertLevels, block % mesh % nCells, block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! block => block % next
+! end do
+! call mpas_timer_stop("se halo h", timer_halo_h)
+! endif ! unsplit_explicit
+ ! update halo for thickness and tracers
+ call mpas_timer_start("se halo h", .false., timer_halo_h)
+ block => domain % blocklist
+ do while (associated(block))
+ call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, block % tend % h % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+! mrp 120131 note: need to remove halo update on tracers themselves below, but later.
+ call mpas_dmpar_exch_halo_field3d_real(domain % dminfo, block % tend % tracers % array(:,:,:), &
+ block % tend % num_tracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
+ call mpas_timer_stop("se halo h", timer_halo_h)
+
block => domain % blocklist
do while (associated(block))
- allocate(hNew(block % mesh % nVertLevels))
+! mrp 120131 remove h computed from SSH in current step
+! allocate(hNew(block % mesh % nVertLevels))
- if (trim(config_new_btr_variables_from) == 'last_subcycle') then
- ! This points to the last barotropic SSH subcycle
- sshNew => block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array
- elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
- ! This points to the tendency variable SSH*
- sshNew => block % state % time_levs(2) % state % ssh % array
- endif
+! mrp 120131 remove h computed from SSH in current step
+! if (trim(config_new_btr_variables_from) == 'last_subcycle') then
+! ! This points to the last barotropic SSH subcycle
+! sshNew => block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array
+! elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
+! ! This points to the tendency variable SSH*
+! sshNew => block % state % time_levs(2) % state % ssh % array
+! endif
- if (trim(config_time_integration) == 'unsplit_explicit') then
+! mrp 120131 remove h computed from SSH in current step
+! if (trim(config_time_integration) == 'unsplit_explicit') then
+! do iCell=1,block % mesh % nCells
+! ! this is h_{n+1}
+! do k = 1,maxLevelCell(iCell)
+! block % state % time_levs(2) % state % h % array(k,iCell) &
+! = block % state % time_levs(1) % state % h % array(k,iCell) &
+! + dt* block % tend % h % array(k,iCell)
+! enddo
- do iCell=1,block % mesh % nCells
- ! this is h_{n+1}
- block % state % time_levs(2) % state % h % array(:,iCell) &
- = block % state % time_levs(1) % state % h % array(:,iCell) &
- + dt* block % tend % h % array(:,iCell)
+! mrp 120131 remove h computed from SSH in current step
+! ! this is only for the hNew computation below, so there is the correct
+! ! value in the ssh variable for unsplit_explicit case.
+! block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+! = block % state % time_levs(2) % state % h % array(1,iCell) &
+! - block % mesh % hZLevel % array(1)
+! end do ! iCell
+
+! mrp 120131 remove h computed from SSH in current step
+! endif ! unsplit_explicit
- ! this is only for the hNew computation below, so there is the correct
- ! value in the ssh variable for unsplit_explicit case.
- block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
- = block % state % time_levs(2) % state % h % array(1,iCell) &
- - block % mesh % hZLevel % array(1)
- end do ! iCell
- endif ! unsplit_explicit
+ !!!!!!!!!!!!!!! not on last iteration: reset variables for next iteration
+ if (split_explicit_step < config_n_ts_iter) then
+
+
! Only need T & S for earlier iterations,
! then all the tracers needed the last time through.
- if (split_explicit_step < config_n_ts_iter) then
- hNew(:) = block % mesh % hZLevel % array(:)
+! mrp 120131 remove h computed from SSH in current step
+! hNew(:) = block % mesh % hZLevel % array(:)
do iCell=1,block % mesh % nCells
! sshNew is a pointer, defined above.
- hNew(1) = sshNew(iCell) + block % mesh % hZLevel % array(1)
+! mrp 120131 remove h computed from SSH in current step
+! hNew(1) = sshNew(iCell) + block % mesh % hZLevel % array(1)
do k=1,block % mesh % maxLevelCell % array(iCell)
+
+ ! this is h_{n+1}
+ temp &
+ = block % state % time_levs(1) % state % h % array(k,iCell) &
+ + dt* block % tend % h % array(k,iCell)
+
+ ! this is h_{n+1/2}
+ block % state % time_levs(2) % state % h % array(k,iCell) &
+ = 0.5*( block % state % time_levs(1) % state % h % array(k,iCell) &
+ + temp)
+
do i=1,2
! This is Phi at n+1
- tracerTemp = ( block % state % time_levs(1) % state % tracers % array(i,k,iCell) &
+ temp = ( block % state % time_levs(1) % state % tracers % array(i,k,iCell) &
* block % state % time_levs(1) % state % h % array(k,iCell) &
- + dt * block % tend % tracers % array(i,k,iCell)) / hNew(k)
+ + dt * block % tend % tracers % array(i,k,iCell)) &
+ / block % state % time_levs(2) % state % h % array(k,iCell)
+! mrp 120131 remove h computed from SSH in current step
+! + dt * block % tend % tracers % array(i,k,iCell)) / hNew(k)
+
! This is Phi at n+1/2
block % state % time_levs(2) % state % tracers % array(i,k,iCell) = 0.5*( &
- block % state % time_levs(1) % state % tracers % array(i,k,iCell) + tracerTemp )
+ block % state % time_levs(1) % state % tracers % array(i,k,iCell) + temp )
enddo
end do
end do ! iCell
- if (trim(config_time_integration) == 'unsplit_explicit') then
+! mrp 120131 remove h computed from SSH in current step
+! if (trim(config_time_integration) == 'unsplit_explicit') then
+! ! compute h*, which is h at n+1/2 and put into array hNew
+! ! on last iteration, hNew remains at n+1
+! do iCell=1,block % mesh % nCells
+! block % state % time_levs(2) % state % h % array(1,iCell) = 0.5*( &
+! block % state % time_levs(2) % state % h % array(1,iCell) &
+! + block % state % time_levs(1) % state % h % array(1,iCell) )
+! end do ! iCell
+! mrp 120131 remove h computed from SSH in current step
+! endif ! unsplit_explicit
- ! compute h*, which is h at n+1/2 and put into array hNew
- ! on last iteration, hNew remains at n+1
- do iCell=1,block % mesh % nCells
- block % state % time_levs(2) % state % h % array(1,iCell) = 0.5*( &
- block % state % time_levs(2) % state % h % array(1,iCell) &
- + block % state % time_levs(1) % state % h % array(1,iCell) )
- end do ! iCell
- endif ! unsplit_explicit
-
+! mrp 120131 remove h computed from SSH in current step
+! this next section is probably not needed - just u^{tr} should be fine.
! 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.
@@ -869,29 +945,67 @@
end do
end do ! iEdge
+ ! uBclNew is u'_{n+1/2}
+ ! uBtrNew is {\bar u}_{avg}
+ ! uNew is u^{tr}
+
+ !!!!!!!!!!!!!!! On last iteration: Compute all variables at time n+1
elseif (split_explicit_step == config_n_ts_iter) then
- hNew(:) = block % mesh % hZLevel % array(:)
+! mrp 120131 remove h computed from SSH in current step
+! hNew(:) = block % mesh % hZLevel % array(:)
do iCell=1,block % mesh % nCells
- ! sshNew is a pointer, defined above.
- hNew(1) = sshNew(iCell) + block % mesh % hZLevel % array(1)
+! mrp 120131 remove h computed from SSH in current step
+! ! sshNew is a pointer, defined above.
+! hNew(1) = sshNew(iCell) + block % mesh % hZLevel % array(1)
do k=1,block % mesh % maxLevelCell % array(iCell)
+
+ ! this is h_{n+1}
+ block % state % time_levs(2) % state % h % array(k,iCell) &
+ = block % state % time_levs(1) % state % h % array(k,iCell) &
+ + dt* block % tend % h % array(k,iCell)
+
do i=1,block % state % time_levs(1) % state % num_tracers
! This is Phi at n+1
block % state % time_levs(2) % state % tracers % array(i,k,iCell) &
= (block % state % time_levs(1) % state % tracers % array(i,k,iCell) &
- * block % state % time_levs(1) % state % h % array(k,iCell) &
- + dt * block % tend % tracers % array(i,k,iCell)) / hNew(k)
+ * block % state % time_levs(1) % state % h % array(k,iCell) &
+ + dt * block % tend % tracers % array(i,k,iCell)) &
+ / block % state % time_levs(2) % state % h % array(k,iCell)
+! mrp 120131 remove h computed from SSH in current step
+! + dt * block % tend % tracers % array(i,k,iCell)) / hNew(k)
+
enddo
end do
end do
+
+ ! Recompute final u to go on to next step.
+ ! u_{n+1} = uBtr_{n+1} + uBcl_{n+1}
+ ! Right now uBclNew is at time n+1/2, so back compute to get uBcl at time n+1
+ ! using uBcl_{n+1/2} = 1/2*(uBcl_n + u_Bcl_{n+1})
+ ! so the following lines are
+ ! u_{n+1} = uBtr_{n+1} + 2*uBcl_{n+1/2} - uBcl_n
+ ! note that uBcl is recomputed at the beginning of the next timestep due to Imp Vert mixing,
+ ! so uBcl does not have to be recomputed here.
+
+ do iEdge=1,block % mesh % nEdges
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ block % state % time_levs(2) % state % u % array(k,iEdge) &
+ = block % state % time_levs(2) % state % uBtr % array( iEdge) &
+ +2*block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
+ - block % state % time_levs(1) % state % uBcl % array(k,iEdge)
+ enddo
+ enddo ! iEdges
+
endif ! split_explicit_step
- deallocate(hNew)
+! mrp 120131 remove h computed from SSH in current step
+ ! deallocate(hNew)
block => block % next
end do
+! mrp 120131 note: need to remove halo update on tracers themselves below, but later.
! Boundary update on tracers. This is placed here, rather than
! on tend % tracers as in RK4, because I needed to update
! afterwards for the del4 diffusion operator.
@@ -906,7 +1020,7 @@
if (split_explicit_step < config_n_ts_iter) then
- ! mrp 110512 I really only need this to compute h_edge, density, pressure.
+ ! mrp 110512 I really only need this to compute h_edge, density, pressure, and SSH
! I can par this down later.
block => domain % blocklist
do while (associated(block))
@@ -922,66 +1036,54 @@
! END large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
- !
+! mrp 120131 remove h computed from SSH in current step
+! !
+! ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
+! !
+
block => domain % blocklist
do while (associated(block))
- if (trim(config_new_btr_variables_from) == 'last_subcycle') then
- do iEdge=1,block % mesh % nEdges
- ! uBtrNew = uBtrSubcycleNew (old here is because counter already flipped)
- ! This line is not needed if u is resplit at the beginning of the timestep.
- block % state % time_levs(2) % state % uBtr % array(iEdge) &
- = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
- enddo ! iEdges
- elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
- ! uBtrNew from u*. this is done above, so u* is already in
- ! block % state % time_levs(2) % state % uBtr % array(iEdge)
- else
- write(0,*) 'Abort: Unknown config_new_btr_variables_from: '//trim(config_time_integration)
- call mpas_dmpar_abort(dminfo)
- endif
- ! Recompute final u to go on to next step.
- ! u_{n+1} = uBtr_{n+1} + uBcl_{n+1}
- ! Right now uBclNew is at time n+1/2, so back compute to get uBcl at time n+1
- ! using uBcl_{n+1/2} = 1/2*(uBcl_n + u_Bcl_{n+1})
- ! so the following lines are
- ! u_{n+1} = uBtr_{n+1} + 2*uBcl_{n+1/2} - uBcl_n
- ! note that uBcl is recomputed at the beginning of the next timestep due to Imp Vert mixing,
- ! so uBcl does not have to be recomputed here.
-
- do iEdge=1,block % mesh % nEdges
- do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
- block % state % time_levs(2) % state % u % array(k,iEdge) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
- + 2*block % state % time_levs(2) % state % uBcl % array(k,iEdge) - block % state % time_levs(1) % state % uBcl % array(k,iEdge)
- enddo
- enddo ! iEdges
+! mrp 120131 remove h computed from SSH in current step
+! if (trim(config_new_btr_variables_from) == 'last_subcycle') then
+! do iEdge=1,block % mesh % nEdges
+! ! uBtrNew = uBtrSubcycleNew (old here is because counter already flipped)
+! ! This line is not needed if u is resplit at the beginning of the timestep.
+! block % state % time_levs(2) % state % uBtr % array(iEdge) &
+! = block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
+! enddo ! iEdges
+! elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
+! ! uBtrNew from u*. this is done above, so u* is already in
+! ! block % state % time_levs(2) % state % uBtr % array(iEdge)
+! else
+! write(0,*) 'Abort: Unknown config_new_btr_variables_from: '//trim(config_time_integration)
+! call mpas_dmpar_abort(dminfo)
+! endif
- if (trim(config_time_integration) == 'split_explicit') then
- if (trim(config_new_btr_variables_from) == 'last_subcycle') then
- do iCell=1,block % mesh % nCells
- ! SSH for the next step is from the end of the barotropic subcycle.
- block % state % time_levs(2) % state % ssh % array(iCell) = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell)
- end do ! iCell
- elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
- ! sshNew from ssh*. This is done above, so ssh* is already in
- ! block % state % time_levs(2) % state % ssh % array(iCell)
- endif
- do iCell=1,block % mesh % nCells
- ! Put new SSH values in h array, for the OcnTendScalar call below.
- block % state % time_levs(2) % state % h % array(1,iCell) = block % state % time_levs(2) % state % ssh % array(iCell) + block % mesh % hZLevel % array(1)
+! mrp 120131 remove h computed from SSH in current step
+! if (trim(config_time_integration) == 'split_explicit') then
+! if (trim(config_new_btr_variables_from) == 'last_subcycle') then
+! do iCell=1,block % mesh % nCells
+! ! SSH for the next step is from the end of the barotropic subcycle.
+! block % state % time_levs(2) % state % ssh % array(iCell) = block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell)
+! end do ! iCell
+! elseif (trim(config_new_btr_variables_from) == 'btr_avg') then
+! ! sshNew from ssh*. This is done above, so ssh* is already in
+! ! block % state % time_levs(2) % state % ssh % array(iCell)
+! endif
+! do iCell=1,block % mesh % nCells
+! ! Put new SSH values in h array, for the OcnTendScalar call below.
+! block % state % time_levs(2) % state % h % array(1,iCell) = block % state % time_levs(2) % state % ssh % array(iCell) + block % mesh % hZLevel % array(1)
+! ! mrp 110601 efficiency note: Since h just moves back and forth between pointers,
+! ! this is not necessary once initialized.
+! do k=2,block % mesh % nVertLevels * sshSwapFlag
+! block % state % time_levs(2) % state % h % array(k,iCell) = block % mesh % hZLevel % array(k)
+! end do
+! end do ! iCell
+! sshSwapFlag = 0
+! end if ! split_explicit
- ! mrp 110601 efficiency note: Since h just moves back and forth between pointers,
- ! this is not necessary once initialized.
- do k=2,block % mesh % nVertLevels * sshSwapFlag
- block % state % time_levs(2) % state % h % array(k,iCell) = block % mesh % hZLevel % array(k)
- end do
- end do ! iCell
- sshSwapFlag = 0
- end if ! split_explicit
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Implicit vertical mixing, done after timestep is complete
@@ -1005,7 +1107,6 @@
!
! Implicit vertical solve for momentum
!
-
call ocn_vel_vmix_tend_implicit(block % mesh, dt, ke_edge, vertvisctopofedge, h, h_edge, u, err)
!
</font>
</pre>