<p><b>lipscomb</b> 2012-01-10 14:05:05 -0700 (Tue, 10 Jan 2012)</p><p><br>
Commented out computations unlikely to be used by land ice model.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_global_diagnostics.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2012-01-10 20:27:25 UTC (rev 1338)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_global_diagnostics.F 2012-01-10 21:05:05 UTC (rev 1339)
@@ -48,7 +48,9 @@
! 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 :: h, u, v, h_edge, &
+!!! pv_edge, pv_vertex, pv_cell, &
+ h_vertex, weightsOnEdge
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -106,9 +108,9 @@
tracers => state % tracers % array
h_edge => state % h_edge % array
h_vertex => state % h_vertex % array
- pv_edge => state % pv_edge % array
- pv_vertex => state % pv_vertex % array
- pv_cell => state % pv_cell % array
+!!! pv_edge => state % pv_edge % array
+!!! pv_vertex => state % pv_vertex % array
+!!! pv_cell => state % pv_cell % array
! Step 2
! 2. Allocate the array with the correct dimensions.
@@ -159,41 +161,41 @@
cellVolume(iLevel,:) = h(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) &
- *h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
- volumeWeightedPotentialEnstrophy(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
- *pv_vertex(iLevel,1:nVerticesSolve)*h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+!!! volumeWeightedPotentialVorticity(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
+!!! *h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+!!! volumeWeightedPotentialEnstrophy(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
+!!! *pv_vertex(iLevel,1:nVerticesSolve)*h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
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)
+!!! 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))
do iEdge = 1,nEdgesSolve
q = 0.0
do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- workpv = 0.5 * (pv_edge(iLevel,iEdge) + pv_edge(iLevel,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * h_edge(iLevel,eoe)
+!!! eoe = edgesOnEdge(j,iEdge)
+!!! workpv = 0.5 * (pv_edge(iLevel,iEdge) + pv_edge(iLevel,eoe))
+!!! q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * h_edge(iLevel,eoe)
end do
- keTend_CoriolisForce(iLevel,iEdge) = h_edge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
+!!! keTend_CoriolisForce(iLevel,iEdge) = h_edge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
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)))
- 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)
- 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) &
- - h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
+!!! 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)
+!!! 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) &
+!!! - h_edge(iLevel,iEdge)*u(iLevel,iEdge)*dvEdge(iEdge)
end do
- peTend_DivThickness(iLevel,:) = peTend_DivThickness(iLevel,1:nCells)*gravity &
- *(h(iLevel,1:nCells)+h_s(1:nCells))
+!!! peTend_DivThickness(iLevel,:) = peTend_DivThickness(iLevel,1:nCells)*gravity &
+!!! *(h(iLevel,1:nCells)+h_s(1:nCells))
end do
do iEdge = 1,nEdgesSolve
@@ -206,15 +208,15 @@
! Step 4
! 4. Call Function to compute Global Stat that you want.
! Computing Kinetic and Potential Energy Tendency Terms
- call land_ice_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, keTend_PressureGradient, globalKineticEnergyTendency)
- call land_ice_compute_global_sum(dminfo, nVertLevels, nCells, peTend_DivThickness, globalPotentialEnergyTendency)
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, keTend_PressureGradient, globalKineticEnergyTendency)
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nCells, peTend_DivThickness, globalPotentialEnergyTendency)
! Computing top and bottom of global mass integral
call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, cellVolume, sumCellVolume)
call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, cellArea, sumCellArea)
- globalKineticEnergyTendency = globalKineticEnergyTendency / sumCellVolume
- globalPotentialEnergyTendency = globalPotentialEnergyTendency / sumCellVolume
+!!! globalKineticEnergyTendency = globalKineticEnergyTendency / sumCellVolume
+!!! globalPotentialEnergyTendency = globalPotentialEnergyTendency / sumCellVolume
! Step 5
! 5. Finish computing the global stat/integral
@@ -227,42 +229,42 @@
averageThickness(:) = (sumrefAreaWeightedSurfaceHeight/sumCellArea)-h_s(1:nCellsSolve)
! Compute Volume Weighted Averages of Potential Vorticity and Potential Enstrophy
- call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
- call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialEnstrophy, globalPotentialEnstrophy)
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialEnstrophy, globalPotentialEnstrophy)
call land_ice_compute_global_sum(dminfo, nVertLevels, nVerticesSolve, vertexVolume, sumVertexVolume)
- globalPotentialVorticity = globalPotentialVorticity/sumVertexVolume
- globalPotentialEnstrophy = globalPotentialEnstrophy/sumVertexVolume
+!!! globalPotentialVorticity = globalPotentialVorticity/sumVertexVolume
+!!! globalPotentialEnstrophy = globalPotentialEnstrophy/sumVertexVolume
- ! Compte Potential Enstrophy Reservior
- potentialEnstrophyReservior(:) = areaCell(:)*fCell(:)*fCell(:)/averageThickness
- call land_ice_compute_global_sum(dminfo, 1, nCellsSolve, potentialEnstrophyReservior, globalPotentialEnstrophyReservoir)
- globalPotentialEnstrophyReservoir = globalPotentialEnstrophyReservoir/sumCellVolume
+ ! Compute Potential Enstrophy Reservior
+!!! potentialEnstrophyReservior(:) = areaCell(:)*fCell(:)*fCell(:)/averageThickness
+!!! call land_ice_compute_global_sum(dminfo, 1, nCellsSolve, potentialEnstrophyReservior, globalPotentialEnstrophyReservoir)
+!!! globalPotentialEnstrophyReservoir = globalPotentialEnstrophyReservoir/sumCellVolume
- globalPotentialEnstrophy = globalPotentialEnstrophy - globalPotentialEnstrophyReservoir
+!!! globalPotentialEnstrophy = globalPotentialEnstrophy - globalPotentialEnstrophyReservoir
! Compute Kinetic and Potential Energy terms to be combined into total energy
- call land_ice_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, volumeWeightedKineticEnergy, globalKineticEnergy)
- call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergy, globalPotentialEnergy)
- call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyTopography, global_temp)
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, volumeWeightedKineticEnergy, globalKineticEnergy)
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergy, globalPotentialEnergy)
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyTopography, global_temp)
- globalKineticEnergy = globalKineticEnergy/sumCellVolume
- globalPotentialEnergy = (globalPotentialEnergy + global_temp)/sumCellVolume
+!!! globalKineticEnergy = globalKineticEnergy/sumCellVolume
+!!! globalPotentialEnergy = (globalPotentialEnergy + global_temp)/sumCellVolume
! Compute Potential energy reservoir to be subtracted from potential energy term
- volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*averageThickness*gravity*0.5
- call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, globalPotentialEnergyReservoir)
- volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*h_s(1:nCellsSolve)*gravity
- call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, global_temp)
+!!! volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*averageThickness*gravity*0.5
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, globalPotentialEnergyReservoir)
+!!! volumeWeightedPotentialEnergyReservoir(1:nCellsSolve) = areaCell(1:nCellsSolve)*averageThickness*h_s(1:nCellsSolve)*gravity
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnergyReservoir, global_temp)
- globalPotentialEnergyReservoir = (globalPotentialEnergyReservoir + global_temp)/sumCellVolume
+!!! globalPotentialEnergyReservoir = (globalPotentialEnergyReservoir + global_temp)/sumCellVolume
- globalPotentialEnergy = globalPotentialEnergy - globalPotentialEnergyReservoir
- globalEnergy = globalKineticEnergy + globalPotentialEnergy
+!!! globalPotentialEnergy = globalPotentialEnergy - globalPotentialEnergyReservoir
+!!! globalEnergy = globalKineticEnergy + globalPotentialEnergy
! Compute Coriolis energy tendency term
- call land_ice_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, keTend_CoriolisForce, globalCoriolisEnergyTendency)
- globalCoriolisEnergyTendency = globalCoriolisEnergyTendency/sumCellVolume
+!!! call land_ice_compute_global_sum(dminfo, nVertLevels, nEdgesSolve, keTend_CoriolisForce, globalCoriolisEnergyTendency)
+!!! globalCoriolisEnergyTendency = globalCoriolisEnergyTendency/sumCellVolume
! Step 6
! 6. Write out your global stat to the file
Modified: branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F
===================================================================
--- branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F 2012-01-10 20:27:25 UTC (rev 1338)
+++ branches/land_ice/mpas/src/core_land_ice/mpas_land_ice_time_integration.F 2012-01-10 21:05:05 UTC (rev 1339)
@@ -88,8 +88,9 @@
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
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(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)
end do
end do
@@ -118,18 +119,18 @@
block => domain % blocklist
do while (associated(block))
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % pv_edge % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+!!! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % pv_edge % array(:,:), &
+!!! block % mesh % nVertLevels, block % mesh % nEdges, &
+!!! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- if (config_h_mom_eddy_visc4 > 0.0) then
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % divergence % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % vorticity % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nVertices, &
- block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- end if
+!!! if (config_h_mom_eddy_visc4 > 0.0) then
+!!! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % divergence % array(:,:), &
+!!! block % mesh % nVertLevels, block % mesh % nCells, &
+!!! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+!!! call mpas_dmpar_exch_halo_field2d_real(domain % dminfo, provis % vorticity % array(:,:), &
+!!! block % mesh % nVertLevels, block % mesh % nVertices, &
+!!! block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+!!! end if
block => block % next
end do
@@ -197,8 +198,8 @@
do iCell=1,block % mesh % nCells
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) &
- + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ + rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
end do
end do
block => block % next
@@ -267,15 +268,18 @@
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 :: vh, weightsOnEdge, kiteAreasOnVertex, &
+ h_edge, h, u, v, tend_h, tend_u, &
+!!! circulation, vorticity, ke, pv_edge, divergence, &
+ h_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r, u_diffusion
+ real (kind=RKIND) :: r
+!!! real (kind=RKIND) :: 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
+!!! 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
@@ -286,11 +290,11 @@
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
- pv_edge => s % pv_edge % array
+!!! circulation => s % circulation % array
+!!! vorticity => s % vorticity % array
+!!! divergence => s % divergence % array
+!!! ke => s % ke % array
+!!! pv_edge => s % pv_edge % array
vh => s % vh % array
weightsOnEdge => grid % weightsOnEdge % array
@@ -348,28 +352,28 @@
!
! Compute u (normal) velocity tendency for each edge (cell face)
!
- tend_u(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+!!! tend_u(:,:) = 0.0
+!!! do iEdge=1,grid % nEdgesSolve
+!!! cell1 = cellsOnEdge(1,iEdge)
+!!! cell2 = cellsOnEdge(2,iEdge)
+!!! vertex1 = verticesOnEdge(1,iEdge)
+!!! 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
+!!! 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)) &
- ) / dcEdge(iEdge)
- end do
- 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)) &
+!!! ) / dcEdge(iEdge)
+!!! end do
+!!! end do
#endif
@@ -378,45 +382,45 @@
!
! 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)
+!!! 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))
+!!! 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))
+!!! 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
+!!! 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="red">abla vorticity
! only valid for visc == constant
- if (config_h_mom_eddy_visc2 > 0.0) then
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+!!! if (config_h_mom_eddy_visc2 > 0.0) then
+!!! do iEdge=1,grid % nEdgesSolve
+!!! cell1 = cellsOnEdge(1,iEdge)
+!!! cell2 = cellsOnEdge(2,iEdge)
+!!! vertex1 = verticesOnEdge(1,iEdge)
+!!! 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 = meshScalingDel2(iEdge) * config_h_mom_eddy_visc2 * u_diffusion
- tend_u(k,iEdge) = tend_u(k,iEdge) + u_diffusion
- end do
- end do
- end if
+!!! do k=1,nVertLevels
+!!! u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+!!! -(vorticity(k,vertex2) - vorticity(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
+!!! end do
+!!! end do
+!!! end if
!
! velocity tendency: del4 dissipation, -</font>
<font color="black">u_4 </font>
<font color="gray">abla^4 u
@@ -424,116 +428,116 @@
! applied recursively.
! 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))
+!!! 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))
- delsq_u(:,:) = 0.0
+!!! delsq_u(:,:) = 0.0
! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+!!! do iEdge=1,grid % nEdges
+!!! cell1 = cellsOnEdge(1,iEdge)
+!!! cell2 = cellsOnEdge(2,iEdge)
+!!! vertex1 = verticesOnEdge(1,iEdge)
+!!! vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+!!! 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(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
+!!! -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
- end do
- end do
+!!! end do
+!!! end do
! vorticity using </font>
<font color="red">abla^2 u
- delsq_circulation(:,:) = 0.0
- do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
- delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
- - dcEdge(iEdge) * delsq_u(k,iEdge)
- delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
- + dcEdge(iEdge) * delsq_u(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
- end do
- end do
+!!! delsq_circulation(:,:) = 0.0
+!!! do iEdge=1,nEdges
+!!! vertex1 = verticesOnEdge(1,iEdge)
+!!! vertex2 = verticesOnEdge(2,iEdge)
+!!! do k=1,nVertLevels
+!!! delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
+!!! - dcEdge(iEdge) * delsq_u(k,iEdge)
+!!! delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
+!!! + dcEdge(iEdge) * delsq_u(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
+!!! end do
+!!! end do
! Divergence using </font>
<font color="red">abla^2 u
- delsq_divergence(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
- delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
- + delsq_u(k,iEdge)*dvEdge(iEdge)
- delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
- - delsq_u(k,iEdge)*dvEdge(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
- end do
- end do
+!!! delsq_divergence(:,:) = 0.0
+!!! do iEdge=1,nEdges
+!!! cell1 = cellsOnEdge(1,iEdge)
+!!! cell2 = cellsOnEdge(2,iEdge)
+!!! do k=1,nVertLevels
+!!! delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
+!!! + delsq_u(k,iEdge)*dvEdge(iEdge)
+!!! delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
+!!! - delsq_u(k,iEdge)*dvEdge(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
+!!! 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
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
+!!! do iEdge=1,grid % nEdgesSolve
+!!! cell1 = cellsOnEdge(1,iEdge)
+!!! cell2 = cellsOnEdge(2,iEdge)
+!!! vertex1 = verticesOnEdge(1,iEdge)
+!!! vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
+!!! 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(k,cell2) &
+!!! - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
+!!! -( delsq_vorticity(k,vertex2) &
+!!! - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
- tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
+!!! u_diffusion = meshScalingDel4(iEdge) * config_h_mom_eddy_visc4 * u_diffusion
+!!! tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
- end do
- end do
+!!! end do
+!!! end do
- deallocate(delsq_divergence)
- deallocate(delsq_u)
- deallocate(delsq_circulation)
- deallocate(delsq_vorticity)
+!!! deallocate(delsq_divergence)
+!!! deallocate(delsq_u)
+!!! deallocate(delsq_circulation)
+!!! deallocate(delsq_vorticity)
- end if
+!!! 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_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
+!!! 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)))
+!!! 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
+!!! 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 land_ice_compute_tend
@@ -715,104 +719,104 @@
!
! 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
+!!! 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
+!!! 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 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)
+!!! do k=1,grid % nVertLevels
+!!! do iTracer=1, grid % nTracers
+!!! ! \kappa_2 </font>
<font color="blue">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
+!!! 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
+!!! end do
- deallocate(boundaryMask)
+!!! deallocate(boundaryMask)
- end if
+!!! 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
+!!! 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(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))
+!!! allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
- delsq_tracer(:,:,:) = 0.
+!!! 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 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
+!!! 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
+!!! 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
+!!! 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 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
+!!! 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
+!!! end do
- deallocate(delsq_tracer)
- deallocate(boundaryMask)
+!!! deallocate(delsq_tracer)
+!!! deallocate(boundaryMask)
- end if
+!!! end if
end subroutine land_ice_compute_scalar_tend
@@ -834,14 +838,17 @@
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, workpv
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
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
+!!! circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, &
+!!! gradPVn, gradPVt, divergence, vorticity_cell, &
+ h_vertex
+ integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, &
+ edgesOnCell, edgesOnEdge, edgesOnVertex, &
+ boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
@@ -856,16 +863,16 @@
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
- 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
+!!! circulation => s % circulation % array
+!!! vorticity => s % vorticity % array
+!!! divergence => s % divergence % array
+!!! ke => s % ke % 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
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -1035,73 +1042,73 @@
!
! Compute circulation and relative vorticity at each vertex
!
- circulation(:,:) = 0.0
- do iEdge=1,nEdges
- 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)
- end do
- end do
- do iVertex=1,nVertices
- do k=1,nVertLevels
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
- end do
- end do
+!!! circulation(:,:) = 0.0
+!!! do iEdge=1,nEdges
+!!! 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)
+!!! end do
+!!! end do
+!!! do iVertex=1,nVertices
+!!! do k=1,nVertLevels
+!!! vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+!!! end do
+!!! end do
!
! Compute the divergence at each cell center
!
- divergence(:,:) = 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)
- enddo
- endif
- if(cell2 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell2) = divergence(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
- enddo
- enddo
+!!! divergence(:,:) = 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)
+!!! enddo
+!!! endif
+!!! if(cell2 <= nCells) then
+!!! do k=1,nVertLevels
+!!! divergence(k,cell2) = divergence(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
+!!! enddo
+!!! enddo
!
! Compute kinetic energy in each cell
!
- ke(:,:) = 0.0
- do iCell=1,nCells
- do i=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i,iCell)
- 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
- end do
+!!! ke(:,:) = 0.0
+!!! do iCell=1,nCells
+!!! do i=1,nEdgesOnCell(iCell)
+!!! iEdge = edgesOnCell(i,iCell)
+!!! 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
+!!! end do
!
! Compute v (tangential) velocities
!
- v(:,:) = 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
- end do
- end do
+!!! v(:,:) = 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
+!!! end do
+!!! end do
#ifdef NCAR_FORMULATION
!
@@ -1131,7 +1138,7 @@
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(k,iVertex)) / h_vertex(k,iVertex)
end do
end do
@@ -1140,78 +1147,78 @@
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
- dvEdge(iEdge)
- enddo
- enddo
+!!! do iEdge = 1,nEdges
+!!! do k = 1,nVertLevels
+!!! gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
+!!! dvEdge(iEdge)
+!!! enddo
+!!! 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
+!!! 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
+!!! 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 )
!
- gradPVn(:,:) = 0.0
- do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
- do k = 1,nVertLevels
- gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
- dcEdge(iEdge)
- enddo
- endif
- enddo
+!!! gradPVn(:,:) = 0.0
+!!! do iEdge = 1,nEdges
+!!! if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
+!!! do k = 1,nVertLevels
+!!! gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
+!!! dcEdge(iEdge)
+!!! enddo
+!!! endif
+!!! enddo
! 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
</font>
</pre>