<p><b>qchen3@fsu.edu</b> 2011-11-06 20:19:21 -0700 (Sun, 06 Nov 2011)</p><p>Branch commit: an augmented FV scheme for SWE; serial version working.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/fvswe/Makefile
===================================================================
--- branches/fvswe/Makefile 2011-11-04 21:22:27 UTC (rev 1173)
+++ branches/fvswe/Makefile 2011-11-07 03:19:21 UTC (rev 1174)
@@ -92,6 +92,18 @@
"CORE = $(CORE)" \
"CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ifort-serial:
+ ( make all \
+ "FC = ifort" \
+ "CC = gcc" \
+ "SFC = ifort" \
+ "SCC = gcc" \
+ "FFLAGS = -real-size 64 -O3 -convert big_endian" \
+ "CFLAGS = -O3 -m64" \
+ "LDFLAGS = -O3" \
+ "CORE = $(CORE)" \
+ "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+
gfortran:
( make all \
"FC = mpif90" \
Copied: branches/fvswe/namelist.input.swe (from rev 1110, branches/fvswe/namelist.input.sw)
===================================================================
--- branches/fvswe/namelist.input.swe (rev 0)
+++ branches/fvswe/namelist.input.swe 2011-11-07 03:19:21 UTC (rev 1174)
@@ -0,0 +1,31 @@
+&sw_model
+ config_test_case = 5
+ config_time_integration = 'RK4'
+ config_dt = 172.8
+ config_ntimesteps = 7500
+ config_output_interval = 500
+ config_stats_interval = 0
+ config_h_ScaleWithMesh = .false.
+ config_h_mom_eddy_visc2 = 0.0
+ config_h_mom_eddy_visc4 = 0.0
+ config_h_tracer_eddy_diff2 = 0.0
+ config_h_tracer_eddy_diff4 = 0.0
+ config_thickness_adv_order = 2
+ config_tracer_adv_order = 2
+ config_positive_definite = .false.
+ config_monotonic = .false.
+ config_wind_stress = .false.
+ config_bottom_drag = .false.
+/
+
+&io
+ config_input_name = 'grid.nc'
+ config_output_name = 'output.nc'
+ config_restart_name = 'restart.nc'
+/
+
+&restart
+ config_restart_interval = 3000
+ config_do_restart = .false.
+ config_restart_time = 1036800.0
+/
Modified: branches/fvswe/src/core_swe/Registry
===================================================================
--- branches/fvswe/src/core_swe/Registry 2011-11-04 21:22:27 UTC (rev 1173)
+++ branches/fvswe/src/core_swe/Registry 2011-11-07 03:19:21 UTC (rev 1174)
@@ -99,7 +99,8 @@
var persistent real fEdge ( nEdges ) 0 iro fEdge mesh - -
var persistent real fVertex ( nVertices ) 0 iro fVertex mesh - -
var persistent real fCell ( nCells ) 0 iro fCell mesh - -
-var persistent real h_s ( nCells ) 0 iro h_s mesh - -
+var persistent real h_s_cell ( nCells ) 0 iro h_s_cell mesh - -
+var persistent real h_s_vertex ( nVertices ) 0 iro h_s_vertex mesh - -
# Space needed for advection
var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 o deriv_two mesh - -
Modified: branches/fvswe/src/core_swe/module_global_diagnostics.F
===================================================================
--- branches/fvswe/src/core_swe/module_global_diagnostics.F 2011-11-04 21:22:27 UTC (rev 1173)
+++ branches/fvswe/src/core_swe/module_global_diagnostics.F 2011-11-07 03:19:21 UTC (rev 1174)
@@ -47,8 +47,8 @@
! Step 1
! 1. Define the array to integrate, and the variable for the value to be stored in after the integration
- real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, h_s, fCell, fEdge
- real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, pv_edge, pv_vertex, pv_cell, h_vertex, weightsOnEdge
+ real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, h_s_cell, fCell, fEdge
+ real (kind=RKIND), dimension(:,:), pointer :: h_cell, u, v, h_edge, pv_edge, pv_vertex, pv_cell, h_vertex, weightsOnEdge
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -86,7 +86,7 @@
nVerticesSolve = grid % nVerticesSolve
nCells = grid % nCells
- h_s => grid % h_s % array
+ h_s_cell => grid % h_s_cell % array
areaCell => grid % areaCell % array
dcEdge => grid % dcEdge % array
dvEdge => grid % dvEdge % array
@@ -100,7 +100,7 @@
areaEdge = dcEdge(1:nEdgesSolve)*dvEdge(1:nEdgesSolve)
weightsOnEdge => grid % weightsOnEdge % array
- h => state % h % array
+ h_cell => state % h_cell % array
u => state % u % array
v => state % v % array
tracers => state % tracers % array
@@ -156,7 +156,7 @@
! eg. GlobalFluidThickness = Sum(h dA)/Sum(dA), See below for array filling
do iLevel = 1,nVertLevels
! eg. GlobalFluidThickness top (Sum( h dA)) = Sum(cellVolume)
- cellVolume(iLevel,:) = h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
+ cellVolume(iLevel,:) = h_cell(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
! eg. GlobalFluidThickness bot (Sum(dA)) = Sum(cellArea)
cellArea(iLevel,:) = areaCell(1:nCellsSolve)
volumeWeightedPotentialVorticity(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
@@ -166,9 +166,9 @@
vertexVolume(iLevel,:) = h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
volumeWeightedKineticEnergy(iLevel,:) = u(iLevel,1:nEdgesSolve)*u(iLevel,1:nEdgesSolve) &
*h_edge(iLevel,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
- volumeWeightedPotentialEnergy(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
- volumeWeightedPotentialEnergyTopography(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h_s(1:nCellsSolve)*areaCell(1:nCellsSolve)
- refAreaWeightedSurfaceHeight(iLevel,:) = areaCell(1:nCellsSolve)*(h(iLevel,1:nCellsSolve)+h_s(1:nCellsSolve))
+ volumeWeightedPotentialEnergy(iLevel,:) = gravity*h_cell(iLevel,1:nCellsSolve)*h_cell(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
+ volumeWeightedPotentialEnergyTopography(iLevel,:) = gravity*h_cell(iLevel,1:nCellsSolve)*h_s_cell(1:nCellsSolve)*areaCell(1:nCellsSolve)
+ refAreaWeightedSurfaceHeight(iLevel,:) = areaCell(1:nCellsSolve)*(h_cell(iLevel,1:nCellsSolve)+h_s_cell(1:nCellsSolve))
do iEdge = 1,nEdgesSolve
q = 0.0
@@ -182,10 +182,10 @@
iCell1 = cellsOnEdge(1,iEdge)
iCell2 = cellsOnEdge(2,iEdge)
- refAreaWeightedSurfaceHeight_edge(iLevel,iEdge) = areaEdge(iEdge)*(h_edge(iLevel,iEdge) + 0.5*(h_s(iCell1) + h_s(iCell2)))
+ refAreaWeightedSurfaceHeight_edge(iLevel,iEdge) = areaEdge(iEdge)*(h_edge(iLevel,iEdge) + 0.5*(h_s_cell(iCell1) + h_s_cell(iCell2)))
keTend_PressureGradient(iLevel,iEdge) = areaEdge(iEdge)*h_edge(iLevel,iEdge)*u(iLevel,iEdge) &
- *gravity*(h(iLevel,iCell2)+h_s(iCell2) - h(iLevel,iCell1)-h_s(iCell1))/dcEdge(iEdge)
+ *gravity*(h_cell(iLevel,iCell2)+h_s_cell(iCell2) - h_cell(iLevel,iCell1)-h_s_cell(iCell1))/dcEdge(iEdge)
peTend_DivThickness(iLevel,iCell1) = peTend_DivThickness(iLevel,iCell1) &
+ h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
peTend_DivThickness(iLevel,iCell2) = peTend_DivThickness(iLevel,iCell2) &
@@ -193,14 +193,14 @@
end do
peTend_DivThickness(iLevel,:) = peTend_DivThickness(iLevel,1:nCells)*gravity &
- *(h(iLevel,1:nCells)+h_s(1:nCells))
+ *(h_cell(iLevel,1:nCells)+h_s_cell(1:nCells))
end do
do iEdge = 1,nEdgesSolve
iCell1 = cellsOnEdge(1,iEdge)
iCell2 = cellsOnEdge(2,iEdge)
- h_s_edge(iEdge) = 0.5*(h_s(iCell1) + h_s(iCell2))
+ h_s_edge(iEdge) = 0.5*(h_s_cell(iCell1) + h_s_cell(iCell2))
end do
! Step 4
@@ -224,7 +224,7 @@
! Reservoir computations
call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, refAreaWeightedSurfaceHeight, sumrefAreaWeightedSurfaceHeight)
- averageThickness(:) = (sumrefAreaWeightedSurfaceHeight/sumCellArea)-h_s(1:nCellsSolve)
+ averageThickness(:) = (sumrefAreaWeightedSurfaceHeight/sumCellArea)-h_s_cell(1:nCellsSolve)
! Compute Volume Weighted Averages of Potential Vorticity and Potential Enstrophy
call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
@@ -252,7 +252,7 @@
! Compute Potential energy reservoir to be subtracted from potential energy term
volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*averageThickness*gravity*0.5
call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, globalPotentialEnergyReservoir)
- volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*h_s(1:nCellsSolve)*gravity
+ volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*h_s_cell(1:nCellsSolve)*gravity
call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, global_temp)
globalPotentialEnergyReservoir = (globalPotentialEnergyReservoir + global_temp)/sumCellVolume
Modified: branches/fvswe/src/core_swe/module_test_cases.F
===================================================================
--- branches/fvswe/src/core_swe/module_test_cases.F 2011-11-04 21:22:27 UTC (rev 1173)
+++ branches/fvswe/src/core_swe/module_test_cases.F 2011-11-07 03:19:21 UTC (rev 1174)
@@ -157,9 +157,9 @@
do iCell=1,grid % nCells
r = sphere_distance(theta_c, lambda_c, grid % latCell % array(iCell), grid % lonCell % array(iCell), a)
if (r < a/3.0) then
- state % h % array(1,iCell) = (h0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
+ state % h_cell % array(1,iCell) = (h0 / 2.0) * (1.0 + cos(pii*r*3.0/a))
else
- state % h % array(1,iCell) = 0.0
+ state % h_cell % array(1,iCell) = 0.0
end if
end do
@@ -187,7 +187,7 @@
integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex, psiCell
!
@@ -213,19 +213,32 @@
! Initialize wind field
!
allocate(psiVertex(grid % nVertices))
+ allocate(psiCell(grid % nCells))
do iVtx=1,grid % nVertices
psiVertex(iVtx) = -a * u0 * ( &
sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
)
end do
+ do iCell=1,grid % nCells
+ psiCell(iCell) = -a * u0 * ( &
+ sin(grid%latCell%array(iCell)) * cos(alpha) - &
+ cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) &
+ )
+ end do
+
do iEdge=1,grid % nEdges
state % u % array(1,iEdge) = -1.0 * ( &
psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
) / grid%dvEdge%array(iEdge)
+ state % v % array(1,iEdge) = 1.0 * ( &
+ psiCell(grid%cellsOnEdge%array(2,iEdge)) - &
+ psiCell(grid%cellsOnEdge%array(1,iEdge)) &
+ ) / grid%dcEdge%array(iEdge)
end do
deallocate(psiVertex)
+ deallocate(psiCell)
!
! Generate rotated Coriolis field
@@ -242,12 +255,18 @@
sin(grid%latVertex%array(iVtx)) * cos(alpha) &
)
end do
+ do iCell=1,grid % nCells
+ grid % fCell % array(iCell) = 2.0 * omega * &
+ (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
+ sin(grid%latCell%array(iCell)) * cos(alpha) &
+ )
+ end do
!
! Initialize height field (actually, fluid thickness field)
!
do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
+ state % h_cell % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
(-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
sin(grid%latCell%array(iCell)) * cos(alpha) &
)**2.0 &
@@ -255,6 +274,15 @@
gravity
end do
+ do iVtx = 1,grid % nVertices
+ state % h_vertex % array(1,iVtx) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
+ (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
+ sin(grid%latVertex%array(iVtx)) * cos(alpha) &
+ )**2.0 &
+ ) / &
+ gravity
+ end do
+
end subroutine sw_test_case_2
@@ -282,7 +310,7 @@
integer :: iCell, iEdge, iVtx
real (kind=RKIND) :: r, u, v
- real (kind=RKIND), allocatable, dimension(:) :: psiVertex
+ real (kind=RKIND), allocatable, dimension(:) :: psiVertex, psiCell
!
@@ -307,19 +335,33 @@
! Initialize wind field
!
allocate(psiVertex(grid % nVertices))
+ allocate(psiCell(grid % nCells))
do iVtx=1,grid % nVertices
psiVertex(iVtx) = -a * u0 * ( &
sin(grid%latVertex%array(iVtx)) * cos(alpha) - &
cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) &
)
end do
+ do iCell=1,grid % nCells
+ psiCell(iCell) = -a * u0 * ( &
+ sin(grid%latCell%array(iCell)) * cos(alpha) - &
+ cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) &
+ )
+ end do
+
+
do iEdge=1,grid % nEdges
state % u % array(1,iEdge) = -1.0 * ( &
psiVertex(grid%verticesOnEdge%array(2,iEdge)) - &
psiVertex(grid%verticesOnEdge%array(1,iEdge)) &
) / grid%dvEdge%array(iEdge)
+ state % v % array(1,iEdge) = 1.0 * ( &
+ psiCell(grid%cellsOnEdge%array(2,iEdge)) - &
+ psiCell(grid%cellsOnEdge%array(1,iEdge)) &
+ ) / grid%dcEdge%array(iEdge)
end do
deallocate(psiVertex)
+ deallocate(psiCell)
!
! Generate rotated Coriolis field
@@ -336,15 +378,26 @@
sin(grid%latVertex%array(iVtx)) * cos(alpha) &
)
end do
+ do iCell=1,grid % nCells
+ grid % fCell % array(iCell) = 2.0 * omega * &
+ (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
+ sin(grid%latCell%array(iCell)) * cos(alpha) &
+ )
+ end do
!
! Initialize mountain
!
- do iCell=1,grid % nCells
+ do iCell = 1,grid % nCells
if (grid % lonCell % array(iCell) < 0.0) grid % lonCell % array(iCell) = grid % lonCell % array(iCell) + 2.0 * pii
r = sqrt(min(rr**2.0, (grid % lonCell % array(iCell) - lambda_c)**2.0 + (grid % latCell % array(iCell) - theta_c)**2.0))
- grid % h_s % array(iCell) = hs0 * (1.0 - r/rr)
+ grid % h_s_cell % array(iCell) = hs0 * (1.0 - r/rr)
end do
+ do iVtx = 1,grid % nVertices
+ if (grid % lonVertex % array(iVtx) < 0.0) grid % lonVertex % array(iVtx) = grid % lonVertex % array(iVtx) + 2.0 * pii
+ r = sqrt(min(rr**2.0, (grid % lonVertex % array(iVtx) - lambda_c)**2.0 + (grid % latVertex % array(iVtx) - theta_c)**2.0))
+ grid % h_s_vertex % array(iVtx) = hs0 * (1.0 - r/rr)
+ end do
!
! Initialize tracer fields
@@ -367,15 +420,25 @@
! Initialize height field (actually, fluid thickness field)
!
do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
+ state % h_cell % array(1,iCell) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
(-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
sin(grid%latCell%array(iCell)) * cos(alpha) &
)**2.0 &
) / &
gravity
- state % h % array(1,iCell) = state % h % array(1,iCell) - grid % h_s % array(iCell)
+ state % h_cell % array(1,iCell) = state % h_cell % array(1,iCell) - grid % h_s_cell % array(iCell)
end do
+ do iVtx = 1,grid % nVertices
+ state % h_vertex % array(1,iVtx) = (gh0 - (a * omega * u0 + 0.5 * u0**2.0) * &
+ (-cos(grid%lonVertex%array(iVtx)) * cos(grid%latVertex%array(iVtx)) * sin(alpha) + &
+ sin(grid%latVertex%array(iVtx)) * cos(alpha) &
+ )**2.0 &
+ ) / &
+ gravity
+ state % h_vertex % array(1,iVtx) = state % h_vertex % array(1,iVtx) - grid % h_s_vertex % array(iVtx)
+ end do
+
end subroutine sw_test_case_5
@@ -442,7 +505,7 @@
! Initialize height field (actually, fluid thickness field)
!
do iCell=1,grid % nCells
- state % h % array(1,iCell) = (gravity * h0 + a*a*AA(grid%latCell%array(iCell)) + &
+ state % h_cell % array(1,iCell) = (gravity * h0 + a*a*AA(grid%latCell%array(iCell)) + &
a*a*BB(grid%latCell%array(iCell)) * cos(R*grid%lonCell%array(iCell)) + &
a*a*CC(grid%latCell%array(iCell)) * cos(2.0*R*grid%lonCell%array(iCell)) &
) / gravity
Modified: branches/fvswe/src/core_swe/module_time_integration.F
===================================================================
--- branches/fvswe/src/core_swe/module_time_integration.F 2011-11-04 21:22:27 UTC (rev 1173)
+++ branches/fvswe/src/core_swe/module_time_integration.F 2011-11-07 03:19:21 UTC (rev 1174)
@@ -181,14 +181,14 @@
if (rk_step < 4) then
block => domain % blocklist
do while (associated(block))
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
- provis % v % array(:,:) = block % state % time_levs(1) % state % v % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % v % array(:,:)
- provis % h_cell % array(:,:) = block % state % time_levs(1) % state % h_cell % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % h_cell % array(:,:)
- provis % h_vertex % array(:,:) = block % state % time_levs(1) % state % h_vertex % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % h_vertex % array(:,:)
+ provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
+ provis % v % array(:,:) = block % state % time_levs(1) % state % v % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % v % array(:,:)
+ provis % h_cell % array(:,:) = block % state % time_levs(1) % state % h_cell % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % h_cell % array(:,:)
+ provis % h_vertex % array(:,:) = block % state % time_levs(1) % state % h_vertex % array(:,:) &
+ + rk_substep_weights(rk_step) * block % tend % h_vertex % array(:,:)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
provis % tracers % array(:,k,iCell) = ( &
@@ -282,18 +282,18 @@
type (mesh_type), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
+ real (kind=RKIND) :: flux, flux_cell, flux_vert, vorticity_abs, workpv, q, upstream_bias
- integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve
+ real (kind=RKIND), dimension(:), pointer :: h_s, h_s_cell, h_s_vertex, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
meshScalingDel2, meshScalingDel4
- real (kind=RKIND), dimension(:,:), pointer :: kiteAreasOnVertex, h_edge, h_cell, h_vertex, u, v, tend_h_cell, tend_h_vertex &
+ real (kind=RKIND), dimension(:,:), pointer :: kiteAreasOnVertex, h_edge, h_cell, h_vertex, u, v, tend_h_cell, tend_h_vertex, &
tend_u, tend_v, circulation_cell, circulation_vertex, &
vorticity_vertex, vorticity_cell, ke_cell, ke_vertex, &
pv_vertex, pv_cell, pv_edge, divergence_vertex, divergence_cell
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r, u_diffusion
+ real (kind=RKIND) :: r, u_diffusion, v_diffusion
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence_cell, delsq_divergence_vertex
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u, delsq_v
@@ -302,7 +302,7 @@
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
- real (kind=RKIND) :: ke_edge
+ !real (kind=RKIND) :: ke_edge
h_cell => s % h_cell % array
@@ -316,12 +316,12 @@
vorticity_vertex => s % vorticity_vertex % array
divergence_cell => s % divergence_cell % array
divergence_vertex => s % divergence_vertex % array
- ke_edge => s % ke_edge % array
+ !ke_edge => s % ke_edge % array
ke_cell => s % ke_cell % array
ke_vertex => s % ke_vertex % array
pv_edge => s % pv_edge % array
- weightsOnEdge => grid % weightsOnEdge % array
+ !weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
cellsOnEdge => grid % cellsOnEdge % array
cellsOnVertex => grid % cellsOnVertex % array
@@ -335,7 +335,8 @@
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
+ h_s_cell => grid % h_s_cell % array
+ h_s_vertex => grid % h_s_vertex % array
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
@@ -348,6 +349,7 @@
nEdges = grid % nEdges
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ nEdgesSolve = grid % nEdgesSolve
u_src => grid % u_src % array
@@ -363,8 +365,8 @@
do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- vert1 = verticesOnEdge(1,iEdge)
- vert2 = verticesOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
do k=1,nVertLevels
flux_cell = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
@@ -372,8 +374,8 @@
tend_h_cell(k,cell1) = tend_h_cell(k,cell1) - flux_cell
tend_h_cell(k,cell2) = tend_h_cell(k,cell2) + flux_cell
- tend_h_vertex(k,vert1) = tend_h_vertex(k,vert1) - flux_vert
- tend_h_vertex(k,vert2) = tend_h_vertex(k,vert2) + flux_vert
+ tend_h_vertex(k,vertex1) = tend_h_vertex(k,vertex1) - flux_vert
+ tend_h_vertex(k,vertex2) = tend_h_vertex(k,vertex2) + flux_vert
end do
end do
@@ -403,11 +405,11 @@
vertex2 = verticesOnEdge(2,iEdge)
do k=1,nVertLevels
- tend_u(k,iEdge) = pv_edge(k,iEdge)*h_edge(k,iEdge)*v(k,iEdge)
+ tend_u(k,iEdge) = pv_edge(k,iEdge)*h_edge(k,iEdge)*v(k,iEdge) &
- ( ke_cell(k,cell2) - ke_cell(k,cell1) + &
gravity * (h_cell(k,cell2) + h_s_cell(cell2) - h_cell(k,cell1) - h_s_cell(cell1)) &
) / dcEdge(iEdge)
- tend_v(k,iEdge) = -pv_edge(k,iEdge)*h_edge(k,iEdge)*u(k,iEdge)
+ tend_v(k,iEdge) = -pv_edge(k,iEdge)*h_edge(k,iEdge)*u(k,iEdge) &
- ( ke_vertex(k,vertex2) - ke_vertex(k,vertex1) + &
gravity * (h_vertex(k,vertex2) + h_s_vertex(vertex2) - h_vertex(k,vertex1) - h_s_vertex(vertex1)) &
) / dvEdge(iEdge)
@@ -599,6 +601,9 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
+ type (tend_type), intent(inout) :: tend
+ type (state_type), intent(in) :: s
+ type (mesh_type), intent(in) :: grid
end subroutine compute_scalar_tend
@@ -623,9 +628,9 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, ke_edge, ke_edge_weighted
- integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fCell, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, kiteAreasOnVertex, h_edge, h_cell, h_vertex, u, v, &
+ integer :: nCells, nEdges, nVertices, nVertLevels, nEdgesSolve, nVerticesSolve, nCellsSolve
+ real (kind=RKIND), dimension(:), pointer :: h_s_cell, fVertex, fCell, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
+ real (kind=RKIND), dimension(:,:), pointer :: kiteAreasOnVertex, h_edge, h_cell, h_vertex, u, v, &
circulation_vertex, circulation_cell, vorticity_vertex, vorticity_cell, &
ke_cell, ke_vertex, pv_edge, pv_vertex, pv_cell, &
gradPVn, gradPVt, divergence_cell, divergence_vertex
@@ -655,7 +660,7 @@
gradPVn => s % gradPVn % array
gradPVt => s % gradPVt % array
- weightsOnEdge => grid % weightsOnEdge % array
+ !weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
cellsOnEdge => grid % cellsOnEdge % array
cellsOnVertex => grid % cellsOnVertex % array
@@ -669,7 +674,7 @@
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
areaTriangle => grid % areaTriangle % array
- h_s => grid % h_s % array
+ h_s_cell => grid % h_s_cell % array
fVertex => grid % fVertex % array
fCell => grid % fCell % array
fEdge => grid % fEdge % array
@@ -679,6 +684,9 @@
nEdges = grid % nEdges
nVertices = grid % nVertices
nVertLevels = grid % nVertLevels
+ nEdgesSolve = grid % nEdgesSolve
+ nVerticesSolve = grid % nVerticesSolve
+ nCellsSolve = grid % nCellsSolve
boundaryEdge => grid % boundaryEdge % array
boundaryCell => grid % boundaryCell % array
@@ -828,13 +836,13 @@
do iCell = 1,nCells
do k = 1,nVertLevels
- ke_cell(k,iCell) = 0.25 * ke_cell(k,iCell)/area_cell(iCell)
+ ke_cell(k,iCell) = 0.25 * ke_cell(k,iCell)/areaCell(iCell)
end do
end do
do iVertex = 1,nVertices
do k = 1,nVertLevels
- ke_vertex(k,iVertex) = 0.25 * ke_vertex(k,iVertex)/area_triangle(iVertex)
+ ke_vertex(k,iVertex) = 0.25 * ke_vertex(k,iVertex)/areaTriangle(iVertex)
end do
end do
</font>
</pre>