<p><b>ringler@lanl.gov</b> 2011-07-28 16:46:56 -0600 (Thu, 28 Jul 2011)</p><p><br>
added global diagnostics to assess potential enstrophy budget<br>
<br>
cleaned up forcing in PV equation<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/src/core_pvsw/Registry
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-07-27 22:14:06 UTC (rev 931)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-07-28 22:46:56 UTC (rev 932)
@@ -25,7 +25,6 @@
namelist logical sw_model config_alpha_closure false
namelist real sw_model config_dragCoefficient 0.001
namelist real sw_model config_enstrophy_injection_rate 0.0
-namelist real sw_model config_enstrophy_forcing_amp 0.0
namelist real sw_model config_Lx 1.0
namelist real sw_model config_Ly 1.0
namelist logical sw_model config_forcing false
@@ -57,6 +56,7 @@
# var persistence type name_in_file ( dims ) time_levs iro- name_in_code struct super-array array_class
#
var persistent real xtime ( Time ) 2 ro xtime state - -
+var persistent real pvForcingAmp ( Time ) 2 ro pvForcingAmp state - -
var persistent real latCell ( nCells ) 0 iro latCell mesh - -
var persistent real lonCell ( nCells ) 0 iro lonCell mesh - -
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_global_diagnostics.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_global_diagnostics.F 2011-07-27 22:14:06 UTC (rev 931)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_global_diagnostics.F 2011-07-28 22:46:56 UTC (rev 932)
@@ -48,7 +48,7 @@
! 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, weightsOnEdge
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, pv, weightsOnEdge, pv_advection, pv_drag, pv_forcing, pv_dissipation
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -59,9 +59,11 @@
real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnstrophy, volumeWeightedKineticEnergy
real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnergy, volumeWeightedPotentialEnergyTopography
real (kind=RKIND), dimension(:,:), allocatable :: keTend_CoriolisForce, keTend_PressureGradient
- real (kind=RKIND), dimension(:,:), allocatable ::peTend_DivThickness, refAreaWeightedSurfaceHeight, refAreaWeightedSurfaceHeight_edge
+ real (kind=RKIND), dimension(:,:), allocatable :: peTend_DivThickness, refAreaWeightedSurfaceHeight, refAreaWeightedSurfaceHeight_edge
+ real (kind=RKIND), dimension(:,:), allocatable :: pe_forcing, pe_drag, pe_dissipation, pe_advection, pe
real (kind=RKIND) :: sumCellVolume, sumCellArea, sumrefAreaWeightedSurfaceHeight
+ real (kind=RKIND) :: sum_pe_forcing, sum_pe_drag, sum_pe_dissipation, sum_pe_advection, sum_pe
real (kind=RKIND) :: globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, globalEnergy
real (kind=RKIND) :: globalCoriolisEnergyTendency, globalKEPETendency, globalPotentialEnstrophyReservoir
@@ -107,6 +109,11 @@
h_edge => state % h_edge % array
pv => state % pv % array
+ pv_forcing => state % pv_forcing % array
+ pv_drag => state % pv_drag % array
+ pv_dissipation => state % pv_dissipation % array
+ pv_advection => state % pv_advection % array
+
! Step 2
! 2. Allocate the array with the correct dimensions.
allocate(cellVolume(nVertLevels,nCellsSolve))
@@ -125,10 +132,15 @@
allocate(peTend_DivThickness(nVertLevels,nCells))
allocate(averageThickness(nCellsSolve))
-
allocate(h_s_edge(nEdgesSOlve))
+ allocate(pe_forcing(nVertLevels,nCellsSolve))
+ allocate(pe_drag(nVertLevels,nCellsSolve))
+ allocate(pe_dissipation(nVertLevels,nCellsSolve))
+ allocate(pe_advection(nVertLevels,nCellsSolve))
+ allocate(pe(nVertLevels,nCellsSolve))
+
cellVolume = 0
refAreaWeightedSurfaceHeight = 0
refAreaWeightedSurfaceHeight_edge = 0
@@ -150,9 +162,7 @@
! 3. Fill the array with the data to be integrated.
! 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)
- ! eg. GlobalFluidThickness bot (Sum(dA)) = Sum(cellArea)
cellArea(iLevel,:) = areaCell(1:nCellsSolve)
volumeWeightedPotentialVorticity(iLevel,:) = pv(iLevel,1:nCellsSolve) &
*h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
@@ -163,6 +173,12 @@
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))
+ pe(iLevel,:) = pv(iLevel,1:nCellsSolve) * pv(iLevel,1:nCellsSolve) * cellVolume(iLevel,:)
+ pe_forcing(iLevel,:) = pv(iLevel,1:nCellsSolve) * pv_forcing(iLevel,1:nCellsSolve) * cellVolume(iLevel,:)
+ pe_drag(iLevel,:) = pv(iLevel,1:nCellsSolve) * pv_drag(iLevel,1:nCellsSolve) * cellVolume(iLevel,:)
+ pe_dissipation(iLevel,:) = pv(iLevel,1:nCellsSolve) * pv_dissipation(iLevel,1:nCellsSolve) * cellVolume(iLevel,:)
+ pe_advection(iLevel,:) = pv(iLevel,1:nCellsSolve) * pv_advection(iLevel,1:nCellsSolve) * cellVolume(iLevel,:)
+
do iEdge = 1,nEdgesSolve
q = 0.0
! do j = 1,nEdgesOnEdge(iEdge)
@@ -206,6 +222,19 @@
call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, cellVolume, sumCellVolume)
call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, cellArea, sumCellArea)
+ ! Computing global potential enstrophy diagnostics
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, pe, sum_pe)
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, pe_forcing, sum_pe_forcing)
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, pe_drag, sum_pe_drag)
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, pe_dissipation, sum_pe_dissipation)
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, pe_advection, sum_pe_advection)
+
+ sum_pe = sum_pe / sumCellVolume
+ sum_pe_forcing = sum_pe_forcing / sumCellVolume
+ sum_pe_drag = sum_pe_drag / sumCellVolume
+ sum_pe_dissipation = sum_pe_dissipation / sumCellVolume
+ sum_pe_advection = sum_pe_advection / sumCellVolume
+
globalKineticEnergyTendency = globalKineticEnergyTendency / sumCellVolume
globalPotentialEnergyTendency = globalPotentialEnergyTendency / sumCellVolume
@@ -232,7 +261,7 @@
call computeGlobalSum(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 computeGlobalSum(dminfo, nVertLevels, nEdgesSolve, volumeWeightedKineticEnergy, globalKineticEnergy)
@@ -268,12 +297,34 @@
open(fileID, file='GlobalIntegrals.txt',POSITION='append')
endif
write(fileID,'(1i0, 100es24.16)') timeIndex, timeIndex*dt, globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, &
- globalEnergy, globalCoriolisEnergyTendency, globalKineticEnergyTendency+globalPotentialEnergyTendency, &
- globalKineticEnergy, globalPotentialEnergy
+ sum_pe, sum_pe_forcing, sum_pe_drag, sum_pe_dissipation, sum_pe_advection, globalKineticEnergy, globalPotentialEnergy
close(fileID)
end if
+
deallocate(areaEdge)
+ deallocate(cellVolume)
+ deallocate(cellArea)
+ deallocate(refAreaWeightedSurfaceHeight)
+ deallocate(refAreaWeightedSurfaceHeight_edge)
+ deallocate(volumeWeightedPotentialVorticity)
+ deallocate(volumeWeightedPotentialEnstrophy)
+ deallocate(potentialEnstrophyReservior)
+ deallocate(volumeWeightedKineticEnergy)
+ deallocate(volumeWeightedPotentialEnergy)
+ deallocate(volumeWeightedPotentialEnergyTopography)
+ deallocate(volumeWeightedPotentialEnergyReservoir)
+ deallocate(keTend_CoriolisForce)
+ deallocate(keTend_PressureGradient)
+ deallocate(peTend_DivThickness)
+ deallocate(averageThickness)
+ deallocate(h_s_edge)
+ deallocate(pe_forcing)
+ deallocate(pe_drag)
+ deallocate(pe_dissipation)
+ deallocate(pe_advection)
+ deallocate(pe)
+
end subroutine computeGlobalDiagnostics
integer function getFreeUnit()
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-07-27 22:14:06 UTC (rev 931)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-07-28 22:46:56 UTC (rev 932)
@@ -34,7 +34,7 @@
dt = config_dt
block => domain % blocklist
do while (associated(block))
- call mpas_init_block(block, block % mesh, dt)
+ call mpas_init_block(domain, block, block % mesh, dt)
block => block % next
end do
@@ -57,20 +57,21 @@
end subroutine mpas_core_init
- subroutine mpas_init_block(block, mesh, dt)
+ subroutine mpas_init_block(domain,block, mesh, dt)
use grid_types
use time_integration
use RBF_interpolation
use vector_reconstruction
+ use global_diagnostics
implicit none
+ type (domain_type), intent(inout) :: domain
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
-
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
call rbfInterp_initialize(mesh)
@@ -81,6 +82,11 @@
call compute_drag(block % state % time_levs(1) % state, mesh)
call compute_forcing(block % state % time_levs(1) % state, mesh)
+ call timer_start("global_diagnostics")
+ call computeGlobalDiagnostics(domain % dminfo, &
+ block % state % time_levs(2) % state, block % mesh, &
+ 0, dt)
+ call timer_stop("global_diagnostics")
if (.not. config_do_restart) block % state % time_levs(1) % state % xtime % scalar = 0.0
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-07-27 22:14:06 UTC (rev 931)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-07-28 22:46:56 UTC (rev 932)
@@ -41,6 +41,7 @@
block => domain % blocklist
do while (associated(block))
block % state % time_levs(2) % state % xtime % scalar = block % state % time_levs(1) % state % xtime % scalar + dt
+ block % state % time_levs(2) % state % pvForcingAmp % scalar = block % state % time_levs(1) % state % pvForcingAmp % scalar
block => block % next
end do
@@ -895,33 +896,22 @@
pv_forcing => s % pv_forcing % array
pv_forcing(:,:) = 0.0
- if(config_forcing) then
+ rtime = s % xtime % scalar
+ amp = s % pvForcingAmp % scalar
- write(6,*) ' inside compute forcing'
+ amp = 1.0e-8
-! following lines need to be removed
+ if(config_forcing) then
- rtime = 0.0
- amp = 1.0e-12
+ wx = 2.0 * pii / (10.0*86400.0)
+ wy = 2.0 * 3.0 / (10.0*86400.0)
-! above lines need to be removed
-
- wx = 2.0*pii / (10.0*86400.0)
- wy = wx * 3.0 / pii
-
phasex = pii * sin(wx*rtime)
phasey = pii * sin(wy*rtime)
kwave = 2.0*pii*4.0 / config_Lx
lwave = 2.0*pii*4.0 / config_Ly
- write(6,*) kwave, lwave
- write(6,*) phasex, phasey
- write(6,*) wx, wy
- write(6,*) config_enstrophy_forcing_amp, amp
- write(6,*) pii
- write(6,*)
-
do iCell = 1,grid%nCells
do k=1,grid % nVertLevels
pv_forcing(k,iCell) = - amp * cos( kwave*xCell(iCell) + phasex)
</font>
</pre>