<p><b>mpetersen@lanl.gov</b> 2012-02-02 15:42:04 -0700 (Thu, 02 Feb 2012)</p><p>BRANCH COMMIT Mostly spacing changes. Tested successfully with zstar and weighted zstar using split explicit and comparing to RK4.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/ale_split_exp/Makefile
===================================================================
--- branches/ocean_projects/ale_split_exp/Makefile 2012-02-02 20:52:11 UTC (rev 1455)
+++ branches/ocean_projects/ale_split_exp/Makefile 2012-02-02 22:42:04 UTC (rev 1456)
@@ -111,7 +111,7 @@
"CC = gcc" \
"SFC = ifort" \
"SCC = gcc" \
- "FFLAGS = -real-size 64 -O0 -check all -convert big_endian -FR" \
+ "FFLAGS = -real-size 64 -O3 -convert big_endian -FR" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
@@ -149,7 +149,7 @@
"CC = gcc" \
"SFC = ifort" \
"SCC = gcc" \
- "FFLAGS = -real-size 64 -O0 -check all -convert big_endian -FR" \
+ "FFLAGS = -real-size 64 -O3 -convert big_endian -FR" \
"CFLAGS = -O3 -m64" \
"LDFLAGS = -O3" \
"CORE = $(CORE)" \
Modified: branches/ocean_projects/ale_split_exp/namelist.input.ocean
===================================================================
--- branches/ocean_projects/ale_split_exp/namelist.input.ocean 2012-02-02 20:52:11 UTC (rev 1455)
+++ branches/ocean_projects/ale_split_exp/namelist.input.ocean 2012-02-02 22:42:04 UTC (rev 1456)
@@ -35,8 +35,6 @@
config_btr_mom_decay_time = 3600.0
config_btr_mom_eddy_visc2 = 0.0
config_btr_subcycle_loop_factor = 2
- config_SSH_from = 'avg_flux'
- config_new_btr_variables_from = 'btr_avg'
config_btr_gam1_uWt1 = 0.5
config_btr_gam2_SSHWt1 = 1.0
config_btr_gam3_uWt2 = 1.0
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-02-02 20:52:11 UTC (rev 1455)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_tendency.F 2012-02-02 22:42:04 UTC (rev 1456)
@@ -236,7 +236,9 @@
! velocity tendency: pressure gradient
!
call mpas_timer_start("pressure grad", .false., velPgradTimer)
- if (config_vert_grid_type.eq.'isopycnal') then
+! mrp changing isopycnal formulation for testing
+ if (1==2) then
+! if (config_vert_grid_type.eq.'isopycnal') then
call ocn_vel_pressure_grad_tend(grid, MontPot, zMid, rho, tend_u, err)
else
call ocn_vel_pressure_grad_tend(grid, pressure, zMid, rho, tend_u, err)
@@ -785,7 +787,9 @@
! This section must be after computing rho
!
! dwj: 10/25/2011 - Need to explore isopycnal vs zlevel flags
- if (config_vert_grid_type.eq.'isopycnal') then
+! mrp changing isopycnal formulation for testing
+ if (1==2) then
+! if (config_vert_grid_type.eq.'isopycnal') then
! For Isopycnal model.
! Compute pressure at top of each layer, and then
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-02-02 20:52:11 UTC (rev 1455)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-02 22:42:04 UTC (rev 1456)
@@ -89,7 +89,7 @@
n_bcl_iter(config_n_ts_iter), &
vertex1, vertex2, iVertex
type (block_type), pointer :: block
- real (kind=RKIND) :: uhSum, hSum, sshEdge, flux, &
+ real (kind=RKIND) :: uhSum, hSum, flux, sshEdge, &
uPerp, uCorr, temp, temp_h, coef, FBtr_coeff, sshCell1, sshCell2
integer :: num_tracers, ucorr_coef, err
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -319,26 +319,29 @@
if (config_filter_btr_mode) then
block % state % time_levs(1) % state % GBtrForcing % array(:) = 0.0
- endif
+ endif
- do iCell=1,block % mesh % nCells
- ! sshSubcycleOld = sshOld
- block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) = block % state % time_levs(1) % state % ssh % array(iCell)
- enddo
+ do iCell=1,block % mesh % nCells
+ ! sshSubcycleOld = sshOld
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ = block % state % time_levs(1) % state % ssh % array(iCell)
+ end do
- do iEdge=1,block % mesh % nEdges
+ do iEdge=1,block % mesh % nEdges
- ! uBtrSubcycleOld = uBtrOld
- block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) = block % state % time_levs(1) % state % uBtr % array(iEdge)
+ ! uBtrSubcycleOld = uBtrOld
+ block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ = block % state % time_levs(1) % state % uBtr % array(iEdge)
- ! uBtrNew = BtrOld This is the first for the summation
- block % state % time_levs(2) % state % uBtr % array(iEdge) = block % state % time_levs(1) % state % uBtr % array(iEdge)
+ ! uBtrNew = BtrOld This is the first for the summation
+ block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ = block % state % time_levs(1) % state % uBtr % array(iEdge)
- ! FBtr = 0
- block % state % time_levs(1) % state % FBtr % array(iEdge) = 0.0
- enddo
+ ! FBtr = 0
+ block % state % time_levs(1) % state % FBtr % array(iEdge) = 0.0
+ end do
- block => block % next
+ block => block % next
end do ! block
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -346,41 +349,41 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do j=1,config_n_btr_subcycles*config_btr_subcycle_loop_factor
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Barotropic subcycle: initial solve for velecity
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- if (config_btr_gam1_uWt1>1.0e-12) then ! only do this part if it is needed in next SSH solve
- uPerpTime = oldBtrSubcycleTime
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: initial solve for velecity
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ if (config_btr_gam1_uWt1>1.0e-12) then ! only do this part if it is needed in next SSH solve
+ uPerpTime = oldBtrSubcycleTime
- block => domain % blocklist
- do while (associated(block))
+ block => domain % blocklist
+ do while (associated(block))
- do iEdge=1,block % mesh % nEdges
+ do iEdge=1,block % mesh % nEdges
- cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
- cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
+ cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
+ cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
- ! 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(uPerpTime) % state % uBtrSubcycle % array(eoe) &
+ ! 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(uPerpTime) % state % uBtrSubcycle % array(eoe) &
* block % mesh % fEdge % array(eoe)
- end do
+ end do
- ! uBtrNew = uBtrOld + dt*(-f*uBtroldPerp - g*grad(SSH) + G)
- block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- = (block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- + dt / config_n_btr_subcycles * (uPerp - gravity &
- * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
- - block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
- / block % mesh % dcEdge % array(iEdge) &
- + block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1, iEdge)
- end do
+ ! uBtrNew = uBtrOld + dt*(-f*uBtroldPerp - g*grad(SSH) + G)
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ = (block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ + dt / config_n_btr_subcycles * (uPerp - gravity &
+ * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) &
+ - block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) ) &
+ / block % mesh % dcEdge % array(iEdge) &
+ + block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1, iEdge)
+ end do
- block => block % next
- end do ! block
+ block => block % next
+ end do ! block
! boundary update on uBtrNew
call mpas_timer_start("se halo ubtr", .false., timer_halo_ubtr)
@@ -396,7 +399,6 @@
call mpas_timer_stop("se halo ubtr", timer_halo_ubtr)
endif ! config_btr_gam1_uWt1>1.0e-12
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Barotropic subcycle: Compute thickness flux and new SSH: PREDICTOR
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -428,7 +430,8 @@
hSum = sum(block % mesh % hZLevel % array (1:block % mesh % maxLevelEdgeTop % array(iEdge)))
flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
- + config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) * (sshEdge + hSum)
+ + config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
+ * (hSum + sshEdge)
block % tend % ssh % array(cell1) = block % tend % ssh % array(cell1) - flux * block % mesh % dvEdge % array(iEdge)
block % tend % ssh % array(cell2) = block % tend % ssh % array(cell2) + flux * block % mesh % dvEdge % array(iEdge)
@@ -563,54 +566,54 @@
block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(:), &
block % mesh % nCells, block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do ! block
- call mpas_timer_stop("se halo ssh", timer_halo_ssh)
- endif ! config_btr_solve_SSH2
+ block => block % next
+ end do ! block
+ call mpas_timer_stop("se halo ssh", timer_halo_ssh)
+ endif ! config_btr_solve_SSH2
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Barotropic subcycle: Accumulate running sums, advance timestep pointers
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Barotropic subcycle: Accumulate running sums, advance timestep pointers
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- block => domain % blocklist
- do while (associated(block))
+ block => domain % blocklist
+ do while (associated(block))
- ! uBtrNew = uBtrNew + uBtrSubcycleNEW
- ! This accumulates the sum.
- ! If the Barotropic Coriolis iteration is limited to one, this could
- ! be merged with the above code.
- do iEdge=1,block % mesh % nEdges
+ ! uBtrNew = uBtrNew + uBtrSubcycleNEW
+ ! This accumulates the sum.
+ ! If the Barotropic Coriolis iteration is limited to one, this could
+ ! be merged with the above code.
+ do iEdge=1,block % mesh % nEdges
- 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)
+ 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 ! iEdge
- block => block % next
- end do ! block
+ end do ! iEdge
+ block => block % next
+ end do ! block
- ! advance time pointers
- oldBtrSubcycleTime = mod(oldBtrSubcycleTime,2)+1
- newBtrSubcycleTime = mod(newBtrSubcycleTime,2)+1
+ ! advance time pointers
+ oldBtrSubcycleTime = mod(oldBtrSubcycleTime,2)+1
+ newBtrSubcycleTime = mod(newBtrSubcycleTime,2)+1
end do ! j=1,config_n_btr_subcycles
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END Barotropic subcycle loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ! Normalize Barotropic subcycle sums: ssh, uBtr, and F
+ ! Normalize Barotropic subcycle sums: ssh, uBtr, and F
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*config_btr_subcycle_loop_factor)
+ 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*config_btr_subcycle_loop_factor)
- block % state % time_levs(2) % state % uBtr % array(iEdge) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
- / (config_n_btr_subcycles*config_btr_subcycle_loop_factor + 1)
- end do
+ block % state % time_levs(2) % state % uBtr % array(iEdge) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ / (config_n_btr_subcycles*config_btr_subcycle_loop_factor + 1)
+ end do
- block => block % next
+ block => block % next
end do ! block
@@ -634,57 +637,56 @@
block => domain % blocklist
do while (associated(block))
- allocate(uTemp(block % mesh % nVertLevels))
+ allocate(uTemp(block % mesh % nVertLevels))
- ! Correction velocity uCorr = (Flux - Sum(h u*))/H
- ! or, for the full latex version:
- !{\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.
+ ! Correction velocity uCorr = (Flux - Sum(h u*))/H
+ ! or, for the full latex version:
+ !{\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
- else
- ucorr_coef = 0
- endif
+ if (config_u_correction) then
+ ucorr_coef = 1
+ else
+ ucorr_coef = 0
+ endif
- do iEdge=1,block % mesh % nEdges
+ do iEdge=1,block % mesh % nEdges
- ! This is u^{avg}
- uTemp(:) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
- + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
+ ! This is u^{avg}
+ uTemp(:) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
+ + block % state % time_levs(2) % state % uBcl % array(:,iEdge)
- uhSum = block % state % time_levs(2) % state % h_edge % array(1,iEdge) * uTemp(1)
- hSum = block % state % time_levs(2) % state % h_edge % array(1,iEdge)
+ uhSum = block % state % time_levs(2) % state % h_edge % array(1,iEdge) * uTemp(1)
+ hSum = block % state % time_levs(2) % state % h_edge % array(1,iEdge)
- do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + block % state % time_levs(2) % state % h_edge % array(k,iEdge) * uTemp(k)
- hSum = hSum + block % state % time_levs(2) % state % h_edge % array(k,iEdge)
- enddo
+ do k=2,block % mesh % maxLevelEdgeTop % array(iEdge)
+ uhSum = uhSum + block % state % time_levs(2) % state % h_edge % array(k,iEdge) * uTemp(k)
+ hSum = hSum + block % state % time_levs(2) % state % h_edge % array(k,iEdge)
+ enddo
- uCorr = ucorr_coef*(( block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum)
+ uCorr = ucorr_coef*(( block % state % time_levs(1) % state % FBtr % array(iEdge) - uhSum)/hSum)
- ! put u^{tr}, the velocity for tracer transport, in uNew
- ! mrp 060611 not sure if boundary enforcement is needed here.
- if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
- block % state % time_levs(2) % state % u % array(:,iEdge) = 0.0
- else
- do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
- block % state % time_levs(2) % state % u % array(k,iEdge) = uTemp(k) + uCorr
- enddo
- do k=block % mesh % maxLevelEdgeTop % array(iEdge)+1,block % mesh % nVertLevels
- block % state % time_levs(2) % state % u % array(k,iEdge) = 0.0
- enddo
- endif
+ ! put u^{tr}, the velocity for tracer transport, in uNew
+ ! mrp 060611 not sure if boundary enforcement is needed here.
+ if (block % mesh % boundaryEdge % array(1,iEdge).eq.1) then
+ block % state % time_levs(2) % state % u % array(:,iEdge) = 0.0
+ else
+ do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
+ block % state % time_levs(2) % state % u % array(k,iEdge) = uTemp(k) + uCorr
+ enddo
+ do k=block % mesh % maxLevelEdgeTop % array(iEdge)+1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % u % array(k,iEdge) = 0.0
+ end do
+ endif
- end do ! iEdge
+ end do ! iEdge
- deallocate(uTemp)
+ deallocate(uTemp)
- block => block % next
+ block => block % next
end do ! block
-
endif ! split_explicit
call mpas_timer_stop("se btr vel", timer_btr_vel)
@@ -697,12 +699,12 @@
block => domain % blocklist
do while (associated(block))
- call ocn_wtop(block % state % time_levs(1) % state,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)
- 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)
+ 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
+ block => block % next
end do
! update halo for thickness and tracer tendencies
@@ -710,12 +712,12 @@
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 % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
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 % 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)
@@ -723,96 +725,99 @@
block => domain % blocklist
do while (associated(block))
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! If iterating, reset variables for next iteration
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- if (split_explicit_step < config_n_ts_iter) then
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! If iterating, 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.
+ do iCell=1,block % mesh % nCells
+ ! sshNew is a pointer, defined above.
+ do k=1,block % mesh % maxLevelCell % array(iCell)
- ! Only need T & S for earlier iterations,
- ! then all the tracers needed the last time through.
- do iCell=1,block % mesh % nCells
- ! sshNew is a pointer, defined above.
- do k=1,block % mesh % maxLevelCell % array(iCell)
+ ! this is h_{n+1}
+ temp_h &
+ = block % state % time_levs(1) % state % h % array(k,iCell) &
+ + dt* block % tend % h % array(k,iCell)
- ! this is h_{n+1}
- temp_h &
- = 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_h)
- ! 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_h)
+ do i=1,2
+ ! This is Phi at n+1
+ 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)) &
+ / temp_h
+
+ ! 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) &
+ + temp )
+ end do
+ end do
+ end do ! iCell
- do i=1,2
- ! This is Phi at n+1
- 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)) &
- / temp_h
+ ! uBclNew is u'_{n+1/2}
+ ! uBtrNew is {\bar u}_{avg}
+ ! uNew is u^{tr}
- ! 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) + temp )
- enddo
- end do
- end do ! iCell
+ ! mrp 110512 I really only need this to compute h_edge, density, pressure, and SSH
+ ! I can par this down later.
+ call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- ! uBclNew is u'_{n+1/2}
- ! uBtrNew is {\bar u}_{avg}
- ! uNew is u^{tr}
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! If large iteration complete, compute all variables at time n+1
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ elseif (split_explicit_step == config_n_ts_iter) then
- ! mrp 110512 I really only need this to compute h_edge, density, pressure, and SSH
- ! I can par this down later.
- call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % maxLevelCell % array(iCell)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! If large iteration complete, compute all variables at time n+1
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- elseif (split_explicit_step == config_n_ts_iter) then
+ ! 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 iCell=1,block % mesh % nCells
- do k=1,block % mesh % maxLevelCell % array(iCell)
+ ! This is Phi at n+1
+ do i=1,block % state % time_levs(1) % state % num_tracers
+ 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)) &
+ / block % state % time_levs(2) % state % h % array(k,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)
+ enddo
+ end do
+ end do
- 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)) &
- / block % state % time_levs(2) % state % h % array(k,iCell)
-
- 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.
+ ! 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
+ 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)
+ end do
+ end do ! iEdges
endif ! split_explicit_step
@@ -824,7 +829,6 @@
! END large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
block => domain % blocklist
do while (associated(block))
@@ -895,98 +899,42 @@
type (diagnostics_type), intent(in) :: d
type (mesh_type), intent(in) :: grid
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-! Some of these variables can be removed, but at a later time.
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j
+ integer :: iEdge, k
integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: vertSum, uhSum, hSum, sshEdge
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, dvEdge, dcEdge, areaCell, areaTriangle, &
- meshScalingDel2, meshScalingDel4
+ real (kind=RKIND) :: vertSum, uhSum, hSum
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
- tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- MontPot, wTop, divergence, vertViscTopOfEdge
+ h_edge, h, u,tend_u
type (dm_info) :: dminfo
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
- edgesOnEdge, edgesOnVertex
- real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+ integer, dimension(:), pointer :: maxLevelEdgeTop
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-
- real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
-
call mpas_timer_start("filter_btr_mode_tend_u")
h => s % h % array
u => s % u % array
- v => s % v % array
- wTop => s % wTop % array
h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- ke_edge => s % ke_edge % array
- pv_edge => s % pv_edge % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
- vertViscTopOfEdge => d % vertViscTopOfEdge % array
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
tend_u => tend % u % array
nCells = grid % nCells
nEdges = grid % nEdges
- nEdgesSolve = grid % nEdgesSolve
- nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- u_src => grid % u_src % array
+ do iEdge=1,nEdges
- do iEdge=1,grid % nEdges
- ! I am using hZLevel here. This assumes that SSH is zero everywhere already,
- ! which should be the case if the barotropic mode is filtered.
- ! The more general case is to use sshEdge or h_edge.
- uhSum = (grid % hZLevel % array(1)) * tend_u(1,iEdge)
- hSum = grid % hZLevel % array(1)
+ uhSum = h_edge(1,iEdge) * tend_u(1,iEdge)
+ hSum = h_edge(1,iEdge)
- do k=2,grid % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + grid % hZLevel % array(k) * tend_u(k,iEdge)
- hSum = hSum + grid % hZLevel % array(k)
+ do k=2,maxLevelEdgeTop(iEdge)
+ uhSum = uhSum + h_edge(k,iEdge) * tend_u(k,iEdge)
+ hSum = hSum + h_edge(k,iEdge)
enddo
vertSum = uhSum/hSum
- do k=1,grid % maxLevelEdgeTop % array(iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) - vertSum
enddo
enddo ! iEdge
@@ -1010,95 +958,40 @@
type (state_type), intent(inout) :: s
type (mesh_type), intent(in) :: grid
-! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
-! Some of these variables can be removed, but at a later time.
- integer :: iEdge, iCell, iVertex, k, cell1, cell2, &
- vertex1, vertex2, eoe, i, j
+ integer :: iEdge, k
integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
- real (kind=RKIND) :: vertSum, uhSum, hSum, sshEdge
- real (kind=RKIND), dimension(:), pointer :: &
- h_s, dvEdge, dcEdge, areaCell, areaTriangle, &
- meshScalingDel2, meshScalingDel4
+ real (kind=RKIND) :: vertSum, uhSum, hSum
real (kind=RKIND), dimension(:,:), pointer :: &
- weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, pressure, &
- tend_u, circulation, vorticity, ke, ke_edge, pv_edge, &
- MontPot, wTop, divergence, vertViscTopOfEdge
+ h_edge, h, u
type (dm_info) :: dminfo
- integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge, &
- maxLevelCell, maxLevelEdgeTop, maxLevelVertexBot
- integer, dimension(:,:), pointer :: &
- cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, &
- edgesOnEdge, edgesOnVertex
- real (kind=RKIND) :: u_diffusion
- real (kind=RKIND), dimension(:), allocatable:: fluxVertTop,w_dudzTopEdge
+ integer, dimension(:), pointer :: maxLevelEdgeTop
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
-
-
- real (kind=RKIND), dimension(:,:), pointer :: u_src
- real (kind=RKIND), parameter :: rho_ref = 1000.0
-
call mpas_timer_start("filter_btr_mode_u")
h => s % h % array
u => s % u % array
- v => s % v % array
- wTop => s % wTop % array
h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
- ke_edge => s % ke_edge % array
- pv_edge => s % pv_edge % array
- MontPot => s % MontPot % array
- pressure => s % pressure % array
- weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
- edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
- dcEdge => grid % dcEdge % array
- dvEdge => grid % dvEdge % array
- areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
- maxLevelCell => grid % maxLevelCell % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
- maxLevelVertexBot => grid % maxLevelVertexBot % array
nCells = grid % nCells
nEdges = grid % nEdges
- nEdgesSolve = grid % nEdgesSolve
- nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
- u_src => grid % u_src % array
+ do iEdge=1,nEdges
- do iEdge=1,grid % nEdges
- ! I am using hZLevel here. This assumes that SSH is zero everywhere already,
- ! which should be the case if the barotropic mode is filtered.
- ! The more general case is to use sshedge or h_edge.
- uhSum = (grid % hZLevel % array(1)) * u(1,iEdge)
- hSum = grid % hZLevel % array(1)
+ uhSum = h_edge(1,iEdge) * u(1,iEdge)
+ hSum = h_edge(1,iEdge)
- do k=2,grid % maxLevelEdgeTop % array(iEdge)
- uhSum = uhSum + grid % hZLevel % array(k) * u(k,iEdge)
- hSum = hSum + grid % hZLevel % array(k)
+ do k=2,maxLevelEdgeTop(iEdge)
+ uhSum = uhSum + h_edge(k,iEdge) * u(k,iEdge)
+ hSum = hSum + h_edge(k,iEdge)
enddo
vertSum = uhSum/hSum
- do k=1,grid % maxLevelEdgeTop % array(iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
u(k,iEdge) = u(k,iEdge) - vertSum
enddo
enddo ! iEdge
Modified: branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vel_pressure_grad.F
===================================================================
--- branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2012-02-02 20:52:11 UTC (rev 1455)
+++ branches/ocean_projects/ale_split_exp/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2012-02-02 22:42:04 UTC (rev 1456)
@@ -187,8 +187,12 @@
err = 0
if (config_vert_grid_type.eq.'isopycnal') then
- rho0Inv = 1.0
- grho0Inv = 0.0
+! mrp changing isopycnal formulation for testing
+ rho0Inv = 1.0/config_rho0
+ grho0Inv = gravity/config_rho0
+! rho0Inv = 1.0
+! grho0Inv = 0.0
+! mrp changing isopycnal formulation for testing end.
else
rho0Inv = 1.0/config_rho0
grho0Inv = gravity/config_rho0
</font>
</pre>