<p><b>ringler@lanl.gov</b> 2011-04-18 15:39:04 -0600 (Mon, 18 Apr 2011)</p><p><br>
the pv-based shallow-water model is working<br>
for the barotropic vorticity equation.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/namelist.input.pvsw
===================================================================
--- branches/pv_based_swm/mpas/namelist.input.pvsw 2011-04-08 22:42:48 UTC (rev 793)
+++ branches/pv_based_swm/mpas/namelist.input.pvsw 2011-04-18 21:39:04 UTC (rev 794)
@@ -1,9 +1,10 @@
&sw_model
- config_test_case = 5
+ configure_bve = .true.
+ config_test_case = 0
config_time_integration = 'RK4'
- config_dt = 172.8
- config_ntimesteps = 7500
- config_output_interval = 500
+ config_dt = 2400
+ config_ntimesteps = 3600
+ config_output_interval = 36
config_stats_interval = 0
config_h_mom_eddy_visc2 = 0.0
config_h_mom_eddy_visc4 = 0.0
Modified: branches/pv_based_swm/mpas/src/core_pvsw/Registry
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-04-08 22:42:48 UTC (rev 793)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-04-18 21:39:04 UTC (rev 794)
@@ -13,6 +13,7 @@
namelist real sw_model config_h_tracer_eddy_diff4 0.0
namelist integer sw_model config_thickness_adv_order 2
namelist integer sw_model config_tracer_adv_order 2
+namelist logical sw_model config_bve false
namelist logical sw_model config_positive_definite false
namelist logical sw_model config_monotonic false
namelist logical sw_model config_wind_stress false
@@ -127,6 +128,8 @@
var persistent real div ( nVertLevels nCells Time ) 2 iro div state - -
var persistent real psi ( nVertLevels nCells Time ) 2 iro psi state - -
var persistent real chi ( nVertLevels nCells Time ) 2 iro chi state - -
+var persistent real psiV ( nVertLevels nVertices Time ) 2 iro psiV state - -
+var persistent real chiV ( nVertLevels nVertices Time ) 2 iro chiV state - -
# Tendency variables
var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
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-04-08 22:42:48 UTC (rev 793)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-04-18 21:39:04 UTC (rev 794)
@@ -27,12 +27,6 @@
write(6,*) 'calling sw_test_case'
if (.not. config_do_restart) call setup_sw_test_case(domain)
- write(6,*) 'calling init_poisson'
- call init_poisson(domain)
-
- write(6,*) 'calling poisson'
- call poisson(domain)
-
!
! Initialize core
!
@@ -43,6 +37,15 @@
block => block % next
end do
+ write(6,*) 'calling init_poisson'
+ call init_poisson(domain)
+
+ write(6,*) 'calling poisson'
+ call poisson(domain)
+
+ write(6,*) ' deriving velocity'
+ call derive_velocity(domain)
+
restart_frame = 1
end subroutine mpas_core_init
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-04-08 22:42:48 UTC (rev 793)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-04-18 21:39:04 UTC (rev 794)
@@ -7,7 +7,8 @@
implicit none
- real (KIND=RKIND), parameter :: threshold = 1.0e-7
+ real (KIND=RKIND), parameter :: threshold = 1.0e-8
+ integer, parameter :: nIterPerGS = 2
public :: init_poisson, gs_iteration, poisson
@@ -46,13 +47,6 @@
block % mesh % poisson_weights % array(j,iCell) = block % mesh % dvEdge % array(jEdge) / &
block % mesh % dcEdge % array(jEdge)
- ! dump weights
- write(10,*) iCell, j, jEdge, jCell, block % mesh % nEdgesOnCell % array(iCell)
- write(10,*) block % mesh % poisson_weights % array(j,iCell)
- write(10,*) block % mesh % dvEdge % array(jEdge)
- write(10,*) block % mesh % dcEdge % array(jEdge)
- write(10,*)
-
if(block % mesh % poisson_weights % array(j,iCell).le.0.0) write(6,*) ' too low '
if(block % mesh % poisson_weights % array(j,iCell).gt.1.0) write(6,*) ' too high'
@@ -119,64 +113,105 @@
integer :: iteration
logical :: converge
- real (kind=RKIND) :: l2_psi_error, l2_chi_error
+ real (kind=RKIND) :: l2_psi_error, l2_chi_error, globalPsiError, globalChiError, globalArea, localArea, rflag, gflag
type (block_type), pointer :: block
iteration = 0; l2_psi_error = 1.0; l2_chi_error = 1.0
converge = .false.
-! update vor and div to get started
-! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vor % array(:,:), &
-! block % mesh % nVertLevels, block % mesh % nCells, &
-! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ ! make sure that mean vor and div are zero
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_mean_vor_and_div(block % state % time_levs(2) % state, block % mesh, l2_psi_error, l2_chi_error, localArea)
+ call dmpar_sum_real(domain % dminfo, l2_psi_error, globalPsiError)
+ call dmpar_sum_real(domain % dminfo, l2_chi_error, globalChiError)
+ call dmpar_sum_real(domain % dminfo, localArea, globalArea)
+ block => block % next
+ end do
+ l2_psi_error = globalPsiError/globalArea
+ l2_chi_error = globalChiError/globalArea
+ write(6,*) ' subtracting mean ', l2_psi_error, l2_chi_error
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(2) % state % vor % array(:,:) = block % state % time_levs(2) % state % vor % array(:,:) - l2_psi_error
+ if(.not.config_bve) block % state % time_levs(2) % state % div % array(:,:) = block % state % time_levs(2) % state % div % array(:,:) - l2_chi_error
+ block => block % next
+ end do
-! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % div % array(:,:), &
-! block % mesh % nVertLevels, block % mesh % nCells, &
-! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
-
+ ! iterate until we converge
do while (.not.converge)
iteration = iteration + 1
+ ! do a few Guass-Seidel iterations
block => domain % blocklist
do while (associated(block))
-
call gs_iteration(block % state % time_levs(2) % state, block % mesh)
- if(mod(iteration,100).eq.0) call compute_poisson_error(block % state % time_levs(2) % state, block % mesh, l2_psi_error, l2_chi_error)
+ block => block % next
+ end do
-! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % psi % 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 % chi % array(:,:), &
-! block % mesh % nVertLevels, block % mesh % nCells, &
-! block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ ! do a boundary update on psi and chi
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % psi % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ if(.not.config_bve) then
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % chi % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ endif
+ block => block % next
+ end do
-! sum error
+ ! check the global error every so often
+ if(mod(iteration,100).eq.0) then
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_poisson_error(block % state % time_levs(2) % state, block % mesh, l2_psi_error, l2_chi_error, localArea)
+ call dmpar_sum_real(domain % dminfo, l2_psi_error, globalPsiError)
+ call dmpar_sum_real(domain % dminfo, l2_chi_error, globalChiError)
+ call dmpar_sum_real(domain % dminfo, localArea, globalArea)
block => block % next
-
end do
-
+ l2_psi_error = sqrt(globalPsiError/globalArea)
+ l2_chi_error = sqrt(globalChiError/globalArea)
if(l2_psi_error.lt.threshold.and.l2_chi_error.lt.threshold) converge=.true.
if(iteration.gt.10000) converge=.true.
+ endif
+
+
+ ! make sure that everyone agrees on convergence
+ if(mod(iteration,100).eq.0) then
block => domain % blocklist
do while (associated(block))
- if(mod(iteration,100).eq.0) write(6,10) iteration, l2_psi_error, l2_chi_error, maxval(block % state % time_levs(2) % state % chi % array(:,:))
- block => block % next
+ rflag = 0.0
+ if(converge) rflag=1.0
+ call dmpar_min_real(domain % dminfo, rflag, gflag)
+ block => block % next
end do
+ converge=.false.
+ if(gflag.gt.0.5) converge=.true.
+ write(6,10) iteration, l2_psi_error, l2_chi_error, rflag, gflag
+ 10 format(i8,6e20.10)
+ endif
- 10 format(i8,3e20.10)
-
+ enddo ! do while not converged
+
+ ! copy psi and chi into the other time level
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % psi % array(:,:) = block % state % time_levs(2) % state % psi % array(:,:)
+ block % state % time_levs(1) % state % chi % array(:,:) = block % state % time_levs(2) % state % chi % array(:,:)
+ block => block % next
end do
- write(6,*) ' exiting poisson'
-
end subroutine poisson
subroutine gs_iteration(state, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: do one iteration using Guass Seidel to solve
+ ! Purpose: do one or more iteration(s) using Guass Seidel to solve
! psi and chi where del2(psi)=vor and del2(chi)=div
!
! Input: vor, div, initial guess for psi and chi
@@ -189,7 +224,7 @@
type (state_type), intent(inout) :: state
type (mesh_type), intent(in) :: grid
- integer :: nCellsSolve
+ integer :: nCells
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: cellsOnCell
integer :: iCell, jCell, iter, k, j, nVertLevels
@@ -197,24 +232,28 @@
real (kind=RKIND), dimension(:,:), pointer :: vor, div, psi, chi
real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
real (kind=RKIND) :: sum1, sum2, sum3
+ real (kind=RKIND), parameter :: sor = 1.5
cellsOnCell => grid % cellsOnCell % array
nEdgesOnCell => grid % nEdgesOnCell % array
areaCell => grid % areaCell % array
- nCellsSolve = grid % nCellsSolve
+ poisson_weights => grid % poisson_weights % array
+ nCells = grid % nCells
nVertLevels = grid % nVertLevels
- poisson_weights => grid % poisson_weights % array
- vor => state % vor % array
+ vor => state % vor % array
div => state % div % array
- psi => state % psi % array
- chi => state % vor % array
+ psi => state % psi % array
+ chi => state % chi % array
+ ! do nIter GS iterations
+ do iter=1,nIterPerGS
+
! loop over the vertical
do k = 1, nVertLevels
! loop over cell centers
- do iCell=1,nCellsSolve
+ do iCell=1,nCells
sum1=0.0; sum2=0.0; sum3=0.0
! loop over neighboring cells
@@ -227,16 +266,18 @@
enddo ! loop over edges
- psi(k,iCell) = (sum2 - vor (k,iCell)*areaCell(iCell))/sum1
- chi(k,iCell) = (sum3 - div(k,iCell)*areaCell(iCell))/sum1
+ psi(k,iCell) = (1.0-sor)*psi(k,iCell) + sor*(sum2 - vor(k,iCell)*areaCell(iCell))/sum1
+ chi(k,iCell) = (1.0-sor)*chi(k,iCell) + sor*(sum3 - div(k,iCell)*areaCell(iCell))/sum1
enddo ! loop over cell centers
enddo ! loop over vertical levels
+ enddo ! loop over iterations
+
end subroutine gs_iteration
- subroutine compute_poisson_error(state, grid, l2_psi_error, l2_chi_error)
+ subroutine compute_poisson_error(state, grid, l2_psi_error, l2_chi_error, localArea)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: compute error in poisson solver
!
@@ -249,7 +290,7 @@
type (state_type), intent(in) :: state
type (mesh_type), intent(in) :: grid
- real (kind=RKIND), intent(out) :: l2_psi_error, l2_chi_error
+ real (kind=RKIND), intent(out) :: l2_psi_error, l2_chi_error, localArea
integer :: nCellsSolve
integer, dimension(:), pointer :: nEdgesOnCell
@@ -267,12 +308,11 @@
nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
- vor => state % vor % array
+ vor => state % vor % array
div => state % div % array
- psi => state % psi % array
- chi => state % vor % array
+ psi => state % psi % array
+ chi => state % chi % array
-
do k=1,nVertLevels
sum_error_psi = 0.0; sum_error_chi = 0.0; sum_area = 0.0
@@ -295,11 +335,194 @@
enddo
- l2_psi_error = sqrt(sum_error_psi/sum_area)
- l2_chi_error = sqrt(sum_error_chi/sum_area)
+ l2_psi_error = sum_error_psi
+ l2_chi_error = sum_error_chi
+ localArea = sum_area
enddo
end subroutine compute_poisson_error
+
+ subroutine compute_mean_vor_and_div(state, grid, l2_psi_error, l2_chi_error, localArea)
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: compute mean vor and div
+ !
+ ! not correct for multiple vertical levels
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(in) :: state
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), intent(out) :: l2_psi_error, l2_chi_error, localArea
+
+ integer :: nCellsSolve
+ integer :: iCell, nVertLevels, k
+ real (kind=RKIND), dimension(:), pointer :: areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: vor, div
+ real (kind=RKIND) :: sum_error_psi, sum_error_chi, sum_area
+
+ areaCell => grid % areaCell % array
+ nCellsSolve = grid % nCellsSolve
+ nVertLevels = grid % nVertLevels
+
+ vor => state % vor % array
+ div => state % div % array
+
+ do k=1,nVertLevels
+
+ sum_error_psi = 0.0; sum_error_chi = 0.0; sum_area = 0.0
+ do iCell=1,nCellsSolve
+
+ sum_error_psi = sum_error_psi + areaCell(iCell)* vor(k,iCell)
+ sum_error_chi = sum_error_chi + areaCell(iCell)* div(k,iCell)
+ sum_area = sum_area + areaCell(iCell)
+
+ enddo
+
+ enddo
+
+ l2_psi_error = sum_error_psi
+ l2_chi_error = sum_error_chi
+ localArea = sum_area
+
+ end subroutine compute_mean_vor_and_div
+
+
+ subroutine derive_velocity(domain)
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ type (block_type), pointer :: block
+
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ call compute_psiV_and_chiV(block % state % time_levs(2) % state, block % mesh)
+ block => block % next
+
+ enddo ! block loop
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ call compute_velocity(block % state % time_levs(2) % state, block % mesh)
+ block => block % next
+
+ enddo
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % v % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+
+ block => block % next
+ end do
+
+
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:)
+ block % state % time_levs(1) % state % v % array(:,:) = block % state % time_levs(2) % state % v % array(:,:)
+ block => block % next
+ end do
+
+ end subroutine derive_velocity
+
+
+ subroutine compute_velocity(state, grid)
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND) :: dcInverse, dvInverse
+ integer :: k, iEdge, cell1, cell2, vertex1, vertex2
+
+ integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge
+ real (kind=RKIND), dimension(:,:), pointer :: u, v, psi, chi, psiV, chiV
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ verticesOnEdge => grid % verticesOnEdge % array
+
+ psi => state % psi % array
+ chi => state % chi % array
+ psiV => state % psiV % array
+ chiV => state % chiV % array
+ u => state % u % array
+ v => state % v % array
+
+ do k=1,grid % nVertLevels
+ do iEdge=1,grid % nEdgesSolve
+
+ dcInverse = 1.0/(grid % dcEdge % array(iEdge))
+ dvInverse = 1.0/(grid % dvEdge % array(iEdge))
+
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ vertex1 = verticesOnEdge(1,iEdge)
+ vertex2 = verticesOnEdge(2,iEdge)
+
+ u(k,iEdge) = (chi(k,cell2)-chi(k,cell1))*dcInverse - (psiV(k,vertex2)-psiV(k,vertex1))*dvInverse
+ v(k,iEdge) = (psi(k,cell2)-psi(k,cell1))*dcInverse + (chiV(k,vertex2)-chiV(k,vertex1))*dvInverse
+
+ enddo
+ enddo
+
+ end subroutine compute_velocity
+
+
+ subroutine compute_psiV_and_chiV(state,grid)
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ integer :: k, iVertex, cell1, cell2, cell3
+
+ integer, dimension(:,:), pointer :: cellsOnVertex
+ real (kind=RKIND), dimension(:,:), pointer :: psi, chi, psiV, chiV
+
+ cellsOnVertex => grid % cellsOnVertex % array
+
+ psi => state % psi % array
+ chi => state % chi % array
+ psiV => state % psiV % array
+ chiV => state % chiV % array
+
+ psiV(:,:) = 0.0
+ chiV(:,:) = 0.0
+
+ do k=1,grid % nVertLevels
+
+ do iVertex=1,grid % nVertices
+
+ cell1 = cellsOnVertex(1,iVertex)
+ cell2 = cellsOnVertex(2,iVertex)
+ cell3 = cellsOnVertex(3,iVertex)
+
+ ! OK for perfect hexagons, can do better for distorted grids
+ psiV(k,iVertex) = (psi(k,cell1)+psi(k,cell2)+psi(k,cell3))/3.0
+ chiV(k,iVertex) = (chi(k,cell1)+chi(k,cell2)+chi(k,cell3))/3.0
+
+ enddo ! k loop
+ enddo ! iVertex loop
+
+
+ end subroutine compute_psiV_and_chiV
+
end module poisson_solver
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_test_cases.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_test_cases.F 2011-04-08 22:42:48 UTC (rev 793)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_test_cases.F 2011-04-18 21:39:04 UTC (rev 794)
@@ -26,7 +26,15 @@
if (config_test_case == 0) then
write(0,*) 'Using initial conditions supplied in input file'
+ block_ptr => domain % blocklist
+ do while (associated(block_ptr))
+ do i=2,nTimeLevs
+ call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
+ end do
+ block_ptr => block_ptr % next
+ end do
+
else if (config_test_case == 1) then
write(0,*) 'Setting up shallow water test case 1'
write(0,*) ' -- Advection of Cosine Bell over the Pole'
@@ -262,20 +270,20 @@
gravity
end do
- ! Initialize psi, chi, pv, vor and div (pv solver)
- do iCell=1,grid % nCells
- state % psi % array(1,iCell) = -a * u0 * ( &
- sin(grid%latCell%array(iCell)) * cos(alpha) - &
- cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) &
- )
- state % chi % array(1,iCell) = 0.0
- state % div % array(1,iCell) = 0.0
- state % pv % array(1,iCell) = ( (2.0*u0/a + 2.0*omega)* &
- (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
- sin(grid%latCell%array(iCell)) * cos(alpha) ) + grid%fCell%array(iCell) ) / &
- state % h % array(1,iCell)
- state % vor % array(1,iCell) = state % h % array(1,iCell) * state % pv % array(1,iCell) - grid%fCell%array(iCell)
- end do
+ ! ! Initialize psi, chi, pv, vor and div (pv solver)
+ ! do iCell=1,grid % nCells
+ ! state % psi % array(1,iCell) = -a * u0 * ( &
+ ! sin(grid%latCell%array(iCell)) * cos(alpha) - &
+ ! cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) &
+ ! )
+ ! state % chi % array(1,iCell) = 0.0
+ ! state % div % array(1,iCell) = 0.0
+ ! state % pv % array(1,iCell) = ( (2.0*u0/a + 2.0*omega)* &
+ ! (-cos(grid%lonCell%array(iCell)) * cos(grid%latCell%array(iCell)) * sin(alpha) + &
+ ! sin(grid%latCell%array(iCell)) * cos(alpha) ) + grid%fCell%array(iCell) ) / &
+ ! state % h % array(1,iCell)
+ ! state % vor % array(1,iCell) = state % h % array(1,iCell) * state % pv % array(1,iCell) - grid%fCell%array(iCell)
+ ! end do
end subroutine sw_test_case_2
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-04-08 22:42:48 UTC (rev 793)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-04-18 21:39:04 UTC (rev 794)
@@ -86,6 +86,13 @@
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(:,:)
+
+ block % state % time_levs(2) % state % pv % array(:,:) = block % state % time_levs(1) % state % pv % array(:,:)
+ block % state % time_levs(2) % state % vor % array(:,:) = block % state % time_levs(1) % state % vor % array(:,:)
+ block % state % time_levs(2) % state % div % array(:,:) = block % state % time_levs(1) % state % div % array(:,:)
+ block % state % time_levs(2) % state % psi % array(:,:) = block % state % time_levs(1) % state % psi % array(:,:)
+ block % state % time_levs(2) % state % chi % array(:,:) = block % state % time_levs(1) % state % chi % 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) &
@@ -116,23 +123,42 @@
! --- update halos for diagnostic variables
- block => domain % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_edge % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ ! block => domain % blocklist
+ ! do while (associated(block))
+ ! call dmpar_exch_halo_field2dReal(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 dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % 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(:,:), &
- block % mesh % nVertLevels, block % mesh % nVertices, &
- block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
- end if
+ ! if (config_h_mom_eddy_visc4 > 0.0) then
+ ! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % 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(:,:), &
+ ! block % mesh % nVertLevels, block % mesh % nVertices, &
+ ! block % parinfo % verticesToSend, block % parinfo % verticesToRecv)
+ ! end if
- block => block % next
+ ! block => block % next
+ ! end do
+
+! --- invert elliptic equations for psi and chi and then derive normal and tangent velocity
+
+ block => domain % blocklist
+ do while (associated(block))
+ do k=1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % vor % array(k,:) = provis % tracers % array(1,k,:) * provis % h % array(k,:) - block % mesh % fCell % array(:)
+ block % state % time_levs(2) % state % div % array(k,:) = 0.0
+ enddo
+ block => block % next
end do
+ call poisson(domain)
+ call derive_velocity(domain)
+ block => domain % blocklist
+ do while (associated(block))
+ provis % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:)
+ provis % v % array(:,:) = block % state % time_levs(2) % state % v % array(:,:)
+ block => block % next
+ end do
! --- compute tendencies
@@ -148,9 +174,9 @@
block => domain % blocklist
do while (associated(block))
- 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 % 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(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
@@ -160,23 +186,11 @@
block => block % next
end do
-! --- invert elliptic equations for psi and chi
- block => domain % blocklist
- do while (associated(block))
- block % state % time_levs(2) % state % vor % array(:,:) = provis % vorticity_cell % array(:,:)
- block % state % time_levs(2) % state % div % array(:,:) = provis % divergence % array(:,:)
- block => block % next
- end do
- write(11,*) 'calling poisson from module_time_integration'
- call poisson(domain)
-
! --- compute next substep state
if (rk_step < 4) then
block => domain % blocklist
do while (associated(block))
- provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
- + rk_substep_weights(rk_step) * block % tend % u % array(:,:)
provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
+ rk_substep_weights(rk_step) * block % tend % h % array(:,:)
do iCell=1,block % mesh % nCells
@@ -184,13 +198,11 @@
provis % tracers % array(:,k,iCell) = ( &
block % state % time_levs(1) % state % h % 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)
+ + rk_substep_weights(rk_step) * block % tend % tracers % array(:,k,iCell) &
+ ) / provis % h % array(k,iCell)
+ provis % pv % array(k,iCell) = provis % tracers % array(1,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(:,:)
- end if
call compute_solve_diagnostics(dt, provis, block % mesh)
block => block % next
end do
@@ -200,14 +212,12 @@
block => domain % blocklist
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(:,:)
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) &
+ block % state % time_levs(2) % state % tracers % array(:,k,iCell) &
+ rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
end do
end do
@@ -223,6 +233,7 @@
!
! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
!
+
block => domain % blocklist
do while (associated(block))
do iCell=1,block % mesh % nCells
@@ -230,13 +241,36 @@
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 % pv % array(k,iCell) = block % state % time_levs(2) % state % tracers % array(1,k,iCell)
end do
end do
+ block => block % next
+ 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(:,:)
- end if
+ ! --- invert elliptic equations for psi and chi and then derive normal and tangent velocity
+ block => domain % blocklist
+ do while (associated(block))
+ do k=1,block % mesh % nVertLevels
+ block % state % time_levs(2) % state % vor % array(k,:) = block % state % time_levs(2) % state % pv % array(k,:) * &
+ block % state % time_levs(2) % state % h % array(k,:) &
+ - block % mesh % fCell % array(:)
+ block % state % time_levs(2) % state % div % array(k,:) = 0.0
+ enddo
+ block => block % next
+ end do
+ call poisson(domain)
+ call derive_velocity(domain)
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:)
+ block % state % time_levs(1) % state % v % array(:,:) = block % state % time_levs(2) % state % v % array(:,:)
+ block => block % next
+ end do
+
+ block => domain % blocklist
+ do while (associated(block))
+
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
call reconstruct(block % state % time_levs(2) % state, block % mesh)
@@ -327,6 +361,7 @@
! Compute height tendency for each cell
!
tend_h(:,:) = 0.0
+ if(config_bve) return
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -342,197 +377,8 @@
end do
end do
-#ifdef LANL_FORMULATION
- !
- ! 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)
-
- 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
-
-
-#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="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)
-
- do k=1,nVertLevels
- u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -(vorticity(k,vertex2) - vorticity(k,vertex1) ) / dvEdge(iEdge)
- u_diffusion = 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="red">abla^4 u
- ! computed as </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="red">abla vorticity
- ! 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))
-
- 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 k=1,nVertLevels
-
- delsq_u(k,iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
-
- 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
-
- ! 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
-
- ! Compute - \kappa </font>
<font color="red">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="red">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 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 = config_h_mom_eddy_visc4 * u_diffusion
- tend_u(k,iEdge) = tend_u(k,iEdge) - u_diffusion
-
- end do
- end do
-
- deallocate(delsq_divergence)
- deallocate(delsq_u)
- deallocate(delsq_circulation)
- deallocate(delsq_vorticity)
-
- 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
</font>
</pre>