<p><b>qchen3@fsu.edu</b> 2011-10-19 22:12:06 -0600 (Wed, 19 Oct 2011)</p><p>Implementing an augmented FV scheme for the shallow water equations.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/fvswe/src/core_swe/Registry
===================================================================
--- trunk/mpas/src/core_sw/Registry 2011-08-15 23:30:26 UTC (rev 940)
+++ branches/fvswe/src/core_swe/Registry 2011-10-20 04:12:06 UTC (rev 1109)
@@ -123,24 +123,30 @@
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 iro u state - -
-var persistent real h ( nVertLevels nCells Time ) 2 iro h state - -
+var persistent real v ( nVertLevels nEdges Time ) 2 iro v state - -
+var persistent real h_cell ( nVertLevels nCells Time ) 2 iro h_cell state - -
+var persistent real h_vertex ( nVertLevels nVertices Time ) 2 iro h_vertex state - -
var persistent real tracers ( nTracers nVertLevels nCells Time ) 2 iro tracers state - -
# Tendency variables
var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
-var persistent real tend_h ( nVertLevels nCells Time ) 1 - h tend - -
+var persistent real tend_v ( nVertLevels nEdges Time ) 1 - v tend - -
+var persistent real tend_h_cell ( nVertLevels nCells Time ) 1 - h_cell tend - -
+var persistent real tend_h_vertex ( nVertLevels nVertices Time ) 1 - h_vertex tend - -
var persistent real tend_tracers ( nTracers nVertLevels nCells Time ) 1 - tracers tend - -
# Diagnostic fields: only written to output
-var persistent real v ( nVertLevels nEdges Time ) 2 o v state - -
-var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
-var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
+var persistent real h_edge ( nVertLevels nEdges Time ) 2 o h_edge state - -
+var persistent real divergence_cell ( nVertLevels nCells Time ) 2 o divergence_cell state - -
+var persistent real divergence_vertex ( nVertLevels nVertices Time ) 2 o divergence_vertex state - -
+var persistent real vorticity_vertex ( nVertLevels nVertices Time ) 2 o vorticity_vertex state - -
var persistent real vorticity_cell ( nVertLevels nCells Time ) 2 o vorticity_cell state - -
+var persistent real vorticity_edge ( nVertLevels nEdges Time ) 2 o vorticity_edge state - -
+var persistent real pv_cell ( nVertLevels nCells Time ) 2 o pv_cell state - -
+var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 o pv_vertex state - -
var persistent real pv_edge ( nVertLevels nEdges Time ) 2 o pv_edge state - -
-var persistent real h_edge ( nVertLevels nEdges Time ) 2 o h_edge state - -
-var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
-var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 o pv_vertex state - -
-var persistent real pv_cell ( nVertLevels nCells Time ) 2 o pv_cell state - -
+var persistent real ke_cell ( nVertLevels nCells Time ) 2 o ke_cell state - -
+var persistent real ke_vertex ( nVertLevels nVertices Time ) 2 o ke_vertex state - -
var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
var persistent real uReconstructY ( nVertLevels nCells Time ) 2 o uReconstructY state - -
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 o uReconstructZ state - -
@@ -148,9 +154,8 @@
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
# Other diagnostic variables: neither read nor written to any files
-var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
-var persistent real circulation ( nVertLevels nVertices Time ) 2 - circulation state - -
+var persistent real circulation_vertex ( nVertLevels nVertices Time ) 2 - circulation_vertex state - -
+var persistent real circulation_cell ( nVertLevels nCells Time ) 2 - circulation_cell state - -
var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
-var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
Modified: branches/fvswe/src/core_swe/module_time_integration.F
===================================================================
--- trunk/mpas/src/core_sw/module_time_integration.F 2011-08-15 23:30:26 UTC (rev 940)
+++ branches/fvswe/src/core_swe/module_time_integration.F 2011-10-20 04:12:06 UTC (rev 1109)
@@ -84,11 +84,13 @@
do while (associated(block))
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
- do iCell=1,block % mesh % nCells ! couple tracers to h
+ block % state % time_levs(2) % state % v % array(:,:) = block % state % time_levs(1) % state % v % array(:,:)
+ block % state % time_levs(2) % state % h_cell % array(:,:) = block % state % time_levs(1) % state % h_cell % array(:,:)
+ block % state % time_levs(2) % state % h_vertex % array(:,:) = block % state % time_levs(1) % state % h_vertex % array(:,:)
+ do iCell=1,block % mesh % nCells ! couple tracers to h_cell
do k=1,block % mesh % nVertLevels
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % state % time_levs(1) % state % h % array(k,iCell)
+ * block % state % time_levs(1) % state % h_cell % array(k,iCell)
end do
end do
@@ -120,14 +122,23 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_edge % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, provis % h_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
if (config_h_mom_eddy_visc4 > 0.0) then
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence_cell % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence_vertex % array(:,:), &
block % mesh % nVertLevels, block % mesh % nVertices, &
block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity_vertex % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorticity_cell % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
end if
block => block % next
@@ -150,9 +161,15 @@
call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % u % array(:,:), &
block % mesh % nVertLevels, block % mesh % nEdges, &
block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % h % array(:,:), &
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % v % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % h_cell % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % h_vertex % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nVertices, &
+ block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
call dmpar_exch_halo_field3dReal(domain % dminfo, block % tend % tracers % array(:,:,:), &
block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
@@ -166,19 +183,24 @@
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 % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % h % 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) = ( &
- block % state % time_levs(1) % state % h % array(k,iCell) * &
+ block % state % time_levs(1) % state % h_cell % array(k,iCell) * &
block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
+ rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
- ) / provis % h % array(k,iCell)
+ ) / provis % h_cell % array(k,iCell)
end do
end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ provis % v % array(:,:) = block % state % time_levs(1) % state % v % array(:,:)
end if
call compute_solve_diagnostics(dt, provis, block % mesh)
block => block % next
@@ -191,8 +213,12 @@
do while (associated(block))
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
+ rk_weights(rk_step) * block % tend % u % array(:,:)
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
- + rk_weights(rk_step) * block % tend % h % array(:,:)
+ block % state % time_levs(2) % state % v % array(:,:) = block % state % time_levs(2) % state % v % array(:,:) &
+ + rk_weights(rk_step) * block % tend % v % array(:,:)
+ block % state % time_levs(2) % state % h_cell % array(:,:) = block % state % time_levs(2) % state % h_cell % array(:,:) &
+ + rk_weights(rk_step) * block % tend % h_cell % array(:,:)
+ block % state % time_levs(2) % state % h_vertex % array(:,:) = block % state % time_levs(2) % state % h_vertex % array(:,:) &
+ + rk_weights(rk_step) * block % tend % h_vertex % array(:,:)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
@@ -218,12 +244,13 @@
do k=1,block % mesh % nVertLevels
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
- / block % state % time_levs(2) % state % h % array(k,iCell)
+ / block % state % time_levs(2) % state % h_cell % array(k,iCell)
end do
end do
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ block % state % time_levs(2) % state % v % array(:,:) = block % state % time_levs(1) % state % v % array(:,:)
end if
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
@@ -260,31 +287,39 @@
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle, &
meshScalingDel2, meshScalingDel4
- real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, divergence, 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), allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence_cell, delsq_divergence_vertex
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u, delsq_v
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation_cell, delsq_circulation_vertex
+ real (kind=RKIND), allocatable, dimension(:,:) :: delsq_vorticity_cell, delsq_vorticity_vertex
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
real (kind=RKIND) :: ke_edge
- h => s % h % array
+ h_cell => s % h_cell % array
+ h_vertex => s % h_vertex % array
+ h_edge => s % h_edge % array
u => s % u % array
v => s % v % array
- h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
+ circulation_cell => s % circulation_cell % array
+ circulation_vertex => s % circulation_vertex % array
+ vorticity_cell => s % vorticity_cell % array
+ 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_cell => s % ke_cell % array
+ ke_vertex => s % ke_vertex % array
pv_edge => s % pv_edge % array
- vh => s % vh % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -304,8 +339,10 @@
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
- tend_h => tend % h % array
+ tend_h_cell => tend % h_cell % array
+ tend_h_vertex => tend % h_vertex % array
tend_u => tend % u % array
+ tend_v => tend % v % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -319,29 +356,46 @@
!
- ! Compute height tendency for each cell
+ ! Compute height tendency for each primary and dual cell
!
- tend_h(:,:) = 0.0
- do iEdge=1,nEdges
+ tend_h_cell(:,:) = 0.0
+ tend_h_vertex(:,:) = 0.0
+ do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ vert1 = verticesOnEdge(1,iEdge)
+ vert2 = verticesOnEdge(2,iEdge)
+
do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell1) = tend_h(k,cell1) - flux
- tend_h(k,cell2) = tend_h(k,cell2) + flux
+ flux_cell = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ flux_vert = v(k,iEdge) * dcEdge(iEdge) * h_edge(k,iEdge)
+
+ 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
+
end do
- end do
+ end do
+
do iCell=1,grid % nCellsSolve
do k=1,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
+ tend_h_cell(k,iCell) = tend_h_cell(k,iCell) / areaCell(iCell)
end do
end do
-#ifdef LANL_FORMULATION
+ do iVertex=1,grid % nVerticesSolve
+ do k=1,nVertLevels
+ tend_h_vertex(k,iVertex) = tend_h_vertex(k,iVertex) / areaTriangle(iVertex)
+ end do
+ end do
+
!
- ! Compute u (normal) velocity tendency for each edge (cell face)
+ ! Compute u and v velocity tendency for each edge (cell face)
!
tend_u(:,:) = 0.0
+ tend_v(:,:) = 0.0
+
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -349,51 +403,20 @@
vertex2 = verticesOnEdge(2,iEdge)
do k=1,nVertLevels
- q = 0.0
- do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
- end do
-
- tend_u(k,iEdge) = &
- q &
- - ( ke(k,cell2) - ke(k,cell1) + &
- gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
+ 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)
- end do
- end do
+ 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)
-
-#endif
-
-#ifdef NCAR_FORMULATION
- !
- ! Compute u (normal) velocity tendency for each edge (cell face)
- !
- tend_u(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,nVertLevels
- vorticity_abs = fEdge(iEdge) + (circulation(k,vertex1) + circulation(k,vertex2)) / &
- (areaTriangle(vertex1) + areaTriangle(vertex2))
-
- workpv = 2.0 * vorticity_abs / (h(k,cell1) + h(k,cell2))
-
- tend_u(k,iEdge) = workpv * vh(k,iEdge) - &
- (ke(k,cell2) - ke(k,cell1) + &
- gravity * (h(k,cell2) + h_s(cell2) - h(k,cell1) - h_s(cell1)) &
- ) / &
- dcEdge(iEdge)
end do
end do
-#endif
- ! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="blue">abla vorticity
+
+ ! Compute diffusion, computed as (</font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity) \cdot n
! only valid for visc == constant
if (config_h_mom_eddy_visc2 > 0.0) then
do iEdge=1,grid % nEdgesSolve
@@ -403,10 +426,15 @@
vertex2 = verticesOnEdge(2,iEdge)
do k=1,nVertLevels
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = ( divergence_cell(k,cell2) - divergence_cell(k,cell1) ) / dcEdge(iEdge) &
+ -(vorticity_vertex(k,vertex2) - vorticity_vertex(k,vertex1) ) / dvEdge(iEdge)
u_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
+
+ v_diffusion = ( divergence_vertex(k,vertex2) - divergence_vertex(k,vertex1) ) / dvEdge(iEdge) &
+ + (vorticity_cell(k,cell2) - vorticity_cell(k,cell1) ) / dcEdge(iEdge)
+ v_diffusion = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * v_diffusion
+ tend_v(k,iEdge) = tend_v(k,iEdge) + v_diffusion
end do
end do
end if
@@ -418,12 +446,17 @@
! strictly only valid for h_mom_eddy_visc4 == constant
!
if (config_h_mom_eddy_visc4 > 0.0) then
- allocate(delsq_divergence(nVertLevels, nCells+1))
allocate(delsq_u(nVertLevels, nEdges+1))
- allocate(delsq_circulation(nVertLevels, nVertices+1))
- allocate(delsq_vorticity(nVertLevels, nVertices+1))
+ allocate(delsq_v(nVertLevels, nEdges+1))
+ allocate(delsq_divergence_cell(nVertLevels, nCells+1))
+ allocate(delsq_divergence_vertex(nVertLevels, nVertices+1))
+ allocate(delsq_vorticity_cell(nVertLevels, nCells+1))
+ allocate(delsq_vorticity_vertex(nVertLevels, nVertices+1))
+ allocate(delsq_circulation_cell(nVertLevels, nCells+1))
+ allocate(delsq_circulation_vertex(nVertLevels, nVertices+1))
delsq_u(:,:) = 0.0
+ delsq_v(:,:) = 0.0
! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity
do iEdge=1,grid % nEdges
@@ -434,50 +467,86 @@
do k=1,nVertLevels
- delsq_u(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+ delsq_u(k,iEdge) = ( divergence_cell(k,cell2) - divergence_cell(k,cell1) ) / dcEdge(iEdge) &
+ - ( vorticity_vertex(k,vertex2) - vorticity_vertex(k,vertex1)) / dvEdge(iEdge)
+ delsq_v(k,iEdge) = ( divergence_vertex(k,vertex2) - divergence_vertex(k,vertex1) ) / dvEdge(iEdge) &
+ + ( vorticity_cell(k,cell2) - vorticity_cell(k,cell1)) / dcEdge(iEdge)
end do
end do
- ! vorticity using </font>
<font color="red">abla^2 u
- delsq_circulation(:,:) = 0.0
+ ! vorticity_cell and vorticity_vertex using \Delta u and \Delta v
+ delsq_circulation_vertex(:,:) = 0.0
+ delsq_circulation_cell(:,:) = 0.0
do iEdge=1,nEdges
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
do k=1,nVertLevels
- delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
+ delsq_circulation_vertex(k,vertex1) = delsq_circulation_vertex(k,vertex1) &
- dcEdge(iEdge) * delsq_u(k,iEdge)
- delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
+ delsq_circulation_vertex(k,vertex2) = delsq_circulation_vertex(k,vertex2) &
+ dcEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_circulation_cell(k,cell1) = delsq_circulation_cell(k,cell1) &
+ + dvEdge(iEdge) * delsq_v(k,iEdge)
+ delsq_circulation_cell(k,cell2) = delsq_circulation_cell(k,cell2) &
+ - dvEdge(iEdge) * delsq_v(k,iEdge)
end do
end do
+
do iVertex=1,nVertices
r = 1.0 / areaTriangle(iVertex)
+
do k=1,nVertLevels
- delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
+ delsq_vorticity_vertex(k,iVertex) = delsq_circulation_vertex(k,iVertex) * r
end do
end do
+ do iCell = 1,nCells
+ r = 1.0 / areaCell(iCell)
+ do k = 1,nVertLevels
+ delsq_vorticity_cell(k,iCell) = delsq_circulation_cell(k,iCell) * r
+ end do
+ end do
+
! Divergence using </font>
<font color="red">abla^2 u
- delsq_divergence(:,:) = 0.0
+ delsq_divergence_cell(:,:) = 0.0
+ delsq_divergence_vertex(:,:) = 0.0
+
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
do k=1,nVertLevels
- delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+ delsq_divergence_cell(k,cell1) = delsq_divergence_cell(k,cell1) &
+ delsq_u(k,iEdge)*dvEdge(iEdge)
- delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
+ delsq_divergence_cell(k,cell2) = delsq_divergence_cell(k,cell2) &
- delsq_u(k,iEdge)*dvEdge(iEdge)
+ delsq_divergence_vertex(k,vertex1) = delsq_divergence_vertex(k,vertex1) &
+ + delsq_v(k,iEdge)*dcEdge(iEdge)
+ delsq_divergence_vertex(k,vertex2) = delsq_divergence_vertex(k,vertex2) &
+ - delsq_v(k,iEdge)*dcEdge(iEdge)
end do
end do
+
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
do k = 1,nVertLevels
- delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
+ delsq_divergence_cell(k,iCell) = delsq_divergence_cell(k,iCell) * r
end do
end do
+ do iVertex = 1,nVertices
+ r = 1.0 / areaTriangle(iVertex)
+ do k = 1,nVertLevels
+ delsq_divergence_vertex(k,iVertex) = delsq_divergence_vertex(k,iVertex) * r
+ end do
+ end do
+
! Compute - \kappa </font>
<font color="black">abla^4 u
! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="gray">abla^2 u) )
do iEdge=1,grid % nEdgesSolve
@@ -488,45 +557,34 @@
do k=1,nVertLevels
- u_diffusion = ( delsq_divergence(k,cell2) &
- - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
- -( delsq_vorticity(k,vertex2) &
- - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = ( delsq_divergence_cell(k,cell2) &
+ - delsq_divergence_cell(k,cell1) ) / dcEdge(iEdge) &
+ -( delsq_vorticity_vertex(k,vertex2) &
+ - delsq_vorticity_vertex(k,vertex1) ) / dvEdge(iEdge)
+ v_diffusion = ( delsq_divergence_vertex(k,vertex2) &
+ - delsq_divergence_vertex(k,vertex1) ) / dvEdge(iEdge) &
+ + ( delsq_vorticity_cell(k,cell2) &
+ - delsq_vorticity_cell(k,cell1) ) / dcEdge(iEdge)
u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
+ v_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * v_diffusion
+
tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
+ tend_v(k,iEdge) = tend_v(k,iEdge) - v_diffusion
end do
end do
- deallocate(delsq_divergence)
+ deallocate(delsq_divergence_cell)
+ deallocate(delsq_divergence_vertex)
deallocate(delsq_u)
- deallocate(delsq_circulation)
- deallocate(delsq_vorticity)
+ deallocate(delsq_v)
+ deallocate(delsq_circulation_cell)
+ deallocate(delsq_circulation_vertex)
+ deallocate(delsq_vorticity_cell)
+ deallocate(delsq_vorticity_vertex)
end if
-
- ! Compute u (velocity) tendency from wind stress (u_src)
- if(config_wind_stress) then
- do iEdge=1,grid % nEdges
- tend_u(1,iEdge) = tend_u(1,iEdge) &
- + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
- end do
- endif
-
- if (config_bottom_drag) then
- do iEdge=1,grid % nEdges
- ! bottom drag is the same as POP:
- ! -c |u| u where c is unitless and 1.0e-3.
- ! see POP Reference guide, section 3.4.4.
- ke_edge = 0.5 * ( ke(1,cellsOnEdge(1,iEdge)) &
- + ke(1,cellsOnEdge(2,iEdge)))
-
- tend_u(1,iEdge) = tend_u(1,iEdge) &
- - 1.0e-3*u(1,iEdge) &
- *sqrt(2.0*ke_edge)/h_edge(1,iEdge)
- end do
- endif
end subroutine compute_tend
@@ -542,271 +600,7 @@
implicit none
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (mesh_type), intent(in) :: grid
- integer :: iCell, iEdge, k, iTracer, cell1, cell2, i
- real (kind=RKIND) :: flux, tracer_edge, r
- real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
- integer, dimension(:,:), pointer :: boundaryEdge
- real (kind=RKIND), dimension(:,:), allocatable :: boundaryMask
- real (kind=RKIND), dimension(:,:,:), allocatable:: delsq_tracer
-
- real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
- integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND) :: coef_3rd_order
- real (kind=RKIND), dimension(:,:), pointer :: u, h_edge
-
- u => s % u % array
- h_edge => s % h_edge % array
- dcEdge => grid % dcEdge % array
- deriv_two => grid % deriv_two % array
- dvEdge => grid % dvEdge % array
- tracers => s % tracers % array
- cellsOnEdge => grid % cellsOnEdge % array
- boundaryCell=> grid % boundaryCell % array
- boundaryEdge=> grid % boundaryEdge % array
- areaCell => grid % areaCell % array
- tracer_tend => tend % tracers % array
-
- coef_3rd_order = 0.
- if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
- if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
-
-
- tracer_tend(:,:,:) = 0.0
-
- if (config_tracer_adv_order == 2) then
-
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- end do
- end do
- end if
- end do
-
- else if (config_tracer_adv_order == 3) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,grid % nTracers
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- !-- else u <= 0:
- else
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- end if
-
- !-- update tendency
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
- end do
-
- else if (config_tracer_adv_order == 4) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if an edge is not on the outer-most ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,grid % nTracers
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- !-- update tendency
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- enddo
- end do
- end if
- end do
-
- endif ! if (config_tracer_adv_order == 2 )
-
- !
- ! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="red">abla \phi)
- !
- if ( config_h_tracer_eddy_diff2 > 0.0 ) then
-
- !
- ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
- !
- allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
- boundaryMask = 1.0
- where(boundaryEdge.eq.1) boundaryMask=0.0
-
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0/areaCell(cell1)
- invAreaCell2 = 1.0/areaCell(cell2)
-
- do k=1,grid % nVertLevels
- do iTracer=1, grid % nTracers
- ! \kappa_2 </font>
<font color="red">abla \phi on edge
- tracer_turb_flux = config_h_tracer_eddy_diff2 &
- *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
-
- ! div(h \kappa_2 </font>
<font color="red">abla \phi) at cell center
- flux = dvEdge(iEdge) * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) - flux * invAreaCell2
- end do
- end do
-
- end do
-
- deallocate(boundaryMask)
-
- end if
-
- !
- ! tracer tendency: del4 horizontal tracer diffusion, &
- ! div(h \kappa_4 </font>
<font color="black">abla [div(h </font>
<font color="red">abla \phi)])
- !
- if ( config_h_tracer_eddy_diff4 > 0.0 ) then
-
- !
- ! compute a boundary mask to enforce insulating boundary conditions in the horizontal
- !
- allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
- boundaryMask = 1.0
- where(boundaryEdge.eq.1) boundaryMask=0.0
-
- allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
-
- delsq_tracer(:,:,:) = 0.
-
- ! first del2: div(h </font>
<font color="red">abla \phi) at cell center
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- do k=1,grid % nVertLevels
- do iTracer=1, grid % nTracers
- delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
- + dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
- delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
- - dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
- end do
- end do
-
- end do
-
- do iCell = 1, grid % nCells
- r = 1.0 / grid % areaCell % array(iCell)
- do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
- end do
- end do
- end do
-
- ! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- invAreaCell1 = 1.0 / grid % areaCell % array(cell1)
- invAreaCell2 = 1.0 / grid % areaCell % array(cell2)
-
- do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) / dcEdge(iEdge)
- flux = dvEdge(iEdge) * tracer_turb_flux
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1 * boundaryMask(k,iEdge)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2 * boundaryMask(k,iEdge)
- end do
- enddo
-
- end do
-
- deallocate(delsq_tracer)
- deallocate(boundaryMask)
-
- end if
-
end subroutine compute_scalar_tend
@@ -827,13 +621,14 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, workpv
+ real (kind=RKIND) :: flux, vorticity_abs, ke_edge, ke_edge_weighted
integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
- h_vertex, vorticity_cell
+ 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, &
+ circulation_vertex, circulation_cell, vorticity_vertex, vorticity_cell, &
+ ke_cell, ke_vertex, pv_edge, pv_vertex, pv_cell, &
+ gradPVn, gradPVt, divergence_cell, divergence_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
@@ -841,22 +636,22 @@
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
- h => s % h % array
+ h_cell => s % h_cell % array
+ h_vertex => s % h_vertex % array
+ h_edge => s % h_edge % array
u => s % u % array
v => s % v % array
- vh => s % vh % array
- h_edge => s % h_edge % array
- h_vertex => s % h_vertex % array
- tend_h => s % h % array
- tend_u => s % u % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
- divergence => s % divergence % array
- ke => s % ke % array
+ circulation_vertex => s % circulation_vertex % array
+ circulation_cell => s % circulation_cell % array
+ vorticity_vertex => s % vorticity_vertex % array
+ vorticity_cell => s % vorticity_cell % array
+ divergence_cell => s % divergence_cell % array
+ divergence_vertex => s % divergence_vertex % array
+ ke_cell => s % ke_cell % array
+ ke_vertex => s % ke_vertex % array
pv_edge => s % pv_edge % array
pv_vertex => s % pv_vertex % array
pv_cell => s % pv_cell % array
- vorticity_cell => s % vorticity_cell % array
gradPVn => s % gradPVn % array
gradPVt => s % gradPVt % array
@@ -876,6 +671,7 @@
areaTriangle => grid % areaTriangle % array
h_s => grid % h_s % array
fVertex => grid % fVertex % array
+ fCell => grid % fCell % array
fEdge => grid % fEdge % array
deriv_two => grid % deriv_two % array
@@ -907,118 +703,19 @@
! Namelist options control the order of accuracy of the reconstructed h_edge value
!
- coef_3rd_order = 0.
- if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
- if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
- if (config_thickness_adv_order == 2) then
+ if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
+ do k=1,grid % nVertLevels
+ h_edge(k,iEdge) = 0.25 * ( h_cell(k,cell1) + h_cell(k,cell2) + h_vertex(k,vertex1) + h_vertex(k,vertex2) )
+ end do
+ end if
+ end do
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end if
- end do
-
- else if (config_thickness_adv_order == 3) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- !-- else u <= 0:
- else
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end if
-
- end do ! do k
- end if ! if (cell1 <=
- end do ! do iEdge
-
- else if (config_thickness_adv_order == 4) then
-
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- end do ! do k
- end if ! if (cell1 <=
- end do ! do iEdge
-
- endif ! if(config_thickness_adv_order == 2)
-
!
! set the velocity in the nEdges+1 slot to zero, this is a dummy address
! used to when reading for edges that do not exist
@@ -1028,107 +725,147 @@
!
! Compute circulation and relative vorticity at each vertex
!
- circulation(:,:) = 0.0
- do iEdge=1,nEdges
+ circulation_vertex(:,:) = 0.0
+ circulation_cell(:,:) = 0.0
+
+ do iEdge=1,nEdgesSolve
do k=1,nVertLevels
- circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
- circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ circulation_vertex(k,vertex1) = circulation_vertex(k,vertex1) - dcEdge(iEdge) * u(k,iEdge)
+ circulation_vertex(k,vertex2) = circulation_vertex(k,vertex2) + dcEdge(iEdge) * u(k,iEdge)
+ circulation_cell(k,cell1) = circulation_cell(k,cell1) + v(k,iEdge) * dvEdge(iEdge)
+ circulation_cell(k,cell2) = circulation_cell(k,cell2) - v(k,iEdge) * dvEdge(iEdge)
end do
end do
- do iVertex=1,nVertices
+
+ do iVertex=1,nVerticesSolve
do k=1,nVertLevels
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+ vorticity_vertex(k,iVertex) = circulation_vertex(k,iVertex) / areaTriangle(iVertex)
end do
end do
+ do iCell = 1, nCellsSolve
+ do k = 1, nVertLevels
+ vorticity_cell(k,iCell) = circulation_cell(k,iCell) / areaCell(iCell)
+ end do
+ end do
+
+
!
- ! Compute the divergence at each cell center
+ ! Compute the divergence at each cell center and at each vertex
!
- divergence(:,:) = 0.0
+ divergence_cell(:,:) = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
if (cell1 <= nCells) then
do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ divergence_cell(k,cell1) = divergence_cell(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
enddo
endif
if(cell2 <= nCells) then
do k=1,nVertLevels
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ divergence_cell(k,cell2) = divergence_cell(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
enddo
end if
end do
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
do k = 1,nVertLevels
- divergence(k,iCell) = divergence(k,iCell) * r
+ divergence_cell(k,iCell) = divergence_cell(k,iCell) * r
enddo
enddo
- !
- ! Compute kinetic energy in each cell
- !
- ke(:,:) = 0.0
- do iCell=1,nCells
- do i=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i,iCell)
+ divergence_vertex(:,:) = 0.0
+ do iEdge=1,nEdges
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+ if (vertex1 <= nVertices) then
do k=1,nVertLevels
- ke(k,iCell) = ke(k,iCell) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
- end do
- end do
- do k=1,nVertLevels
- ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- end do
+ divergence_vertex(k,vertex1) = divergence_vertex(k,vertex1) + v(k,iEdge)*dcEdge(iEdge)
+ enddo
+ endif
+ if(vertex2 <= nVertices) then
+ do k=1,nVertLevels
+ divergence_vertex(k,vertex2) = divergence_vertex(k,vertex2) - v(k,iEdge)*dcEdge(iEdge)
+ enddo
+ end if
end do
+ do iVertex = 1,nVertices
+ r = 1.0 / areaTriangle(iVertex)
+ do k = 1,nVertLevels
+ divergence_vertex(k,iVertex) = divergence_vertex(k,iVertex) * r
+ enddo
+ enddo
+
!
- ! Compute v (tangential) velocities
+ ! Compute kinetic energy in each cell and in each dual cell
!
- v(:,:) = 0.0
+ ke_cell(:,:) = 0.0
+ ke_vertex(:,:) = 0.0
do iEdge = 1,nEdges
- do i=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(i,iEdge)
- do k = 1,nVertLevels
- v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
- end do
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ do k = 1,nVertLevels
+ ke_edge = 0.5*(u(k,iEdge)**2 + v(k,iEdge)**2)
+ ke_edge_weighted = ke_edge * dcEdge(iEdge)*dvEdge(iEdge)
+
+ ke_cell(k,cell1) = ke_cell(k,cell1) + ke_edge_weighted
+ ke_cell(k,cell2) = ke_cell(k,cell2) + ke_edge_weighted
+ ke_vertex(k,vertex1) = ke_vertex(k,vertex1) + ke_edge_weighted
+ ke_vertex(k,vertex2) = ke_vertex(k,vertex2) + ke_edge_weighted
end do
end do
-#ifdef NCAR_FORMULATION
- !
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
- vh(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- do j=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- do k=1,nVertLevels
- vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * u(k,eoe) * h_edge(k,eoe)
- end do
+ do iCell = 1,nCells
+ do k = 1,nVertLevels
+ ke_cell(k,iCell) = 0.25 * ke_cell(k,iCell)/area_cell(iCell)
end do
end do
-#endif
+ do iVertex = 1,nVertices
+ do k = 1,nVertLevels
+ ke_vertex(k,iVertex) = 0.25 * ke_vertex(k,iVertex)/area_triangle(iVertex)
+ end do
+ end do
+
!
- ! Compute height at vertices, pv at vertices, and average pv to edge locations
+ ! Compute pv at vertices and cells, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
- do iVertex = 1,nVertices
+ do iVertex = 1,nVerticesSolve
do k=1,nVertLevels
- h_vertex(k,iVertex) = 0.0
- do i=1,grid % vertexDegree
- h_vertex(k,iVertex) = h_vertex(k,iVertex) + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
- end do
- h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)
-
- pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
+ pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity_vertex(k,iVertex)) / h_vertex(k,iVertex)
end do
end do
+ do iCell = 1, nCellsSolve
+ do k = 1, nVertLevels
+ pv_cell(k,iCell) = (fCell(iCell) + vorticity_cell(k,iCell)) / h_cell(k,iCell)
+ enddo
+ end do
+ do iEdge = 1, nEdgesSolve
+ do k = 1, nVertLevels
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ pv_edge(k,iEdge) = 0.25*(pv_cell(k,cell1) + pv_cell(k,cell2) + pv_vertex(k,vertex1) + pv_vertex(k,vertex2) )
+ end do
+ end do
+
!
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
@@ -1141,50 +878,16 @@
enddo
!
- ! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells )
- !
- pv_edge(:,:) = 0.0
- do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
- do k=1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
- end do
- end do
- end do
-
-
- !
! Modify PV edge with upstream bias.
!
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
- enddo
- enddo
+! do iEdge = 1,nEdges
+! do k = 1,nVertLevels
+! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
+! enddo
+! enddo
!
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
- !
- pv_cell(:,:) = 0.0
- vorticity_cell(:,:) = 0.0
- do iVertex = 1, nVertices
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- if (iCell <= nCells) then
- do k = 1,nVertLevels
- pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
- vorticity_cell(k,iCell) = vorticity_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) / areaCell(iCell)
- enddo
- endif
- enddo
- enddo
-
-
- !
! Compute gradient of PV in normal direction
! ( this computes gradPVn for all edges bounding real cells )
!
@@ -1200,37 +903,13 @@
! Modify PV edge with upstream bias.
!
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
- enddo
- enddo
+! do iEdge = 1,nEdges
+! do k = 1,nVertLevels
+! pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
+! enddo
+! enddo
- !
- ! set pv_edge = fEdge / h_edge at boundary points
- !
- ! if (maxval(boundaryEdge).ge.0) then
- ! do iEdge = 1,nEdges
- ! cell1 = cellsOnEdge(1,iEdge)
- ! cell2 = cellsOnEdge(2,iEdge)
- ! do k = 1,nVertLevels
- ! if(boundaryEdge(k,iEdge).eq.1) then
- ! v(k,iEdge) = 0.0
- ! if(cell1.gt.0) then
- ! h1 = h(k,cell1)
- ! pv_edge(k,iEdge) = fEdge(iEdge) / h1
- ! h_edge(k,iEdge) = h1
- ! else
- ! h2 = h(k,cell2)
- ! pv_edge(k,iEdge) = fEdge(iEdge) / h2
- ! h_edge(k,iEdge) = h2
- ! endif
- ! endif
- ! enddo
- ! enddo
- ! endif
-
end subroutine compute_solve_diagnostics
</font>
</pre>