<p><b>ringler@lanl.gov</b> 2011-04-21 21:33:28 -0600 (Thu, 21 Apr 2011)</p><p><br>
a working version of the alpha closure<br>
(currently only working for the barotropic vorticity equation)<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/src/core_pvsw/Makefile
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Makefile 2011-04-21 19:55:42 UTC (rev 800)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Makefile 2011-04-22 03:33:28 UTC (rev 801)
@@ -4,6 +4,7 @@
module_test_cases.o \
module_advection.o \
module_poisson_solver.o \
+ module_alpha_solver.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -18,8 +19,10 @@
module_poisson_solver.o:
-module_time_integration.o: module_poisson_solver.o
+module_alpha_solver.o: module_poisson_solver.o
+module_time_integration.o: module_poisson_solver.o module_alpha_solver.o
+
module_global_diagnostics.o:
module_mpas_core.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o module_advection.o module_poisson_solver.o
Modified: branches/pv_based_swm/mpas/src/core_pvsw/Registry
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-04-21 19:55:42 UTC (rev 800)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-04-22 03:33:28 UTC (rev 801)
@@ -19,6 +19,8 @@
namelist logical sw_model config_wind_stress false
namelist logical sw_model config_bottom_drag false
namelist real sw_model config_apvm_upwinding 0.5
+namelist real sw_model config_alpha 0.0
+namelist logical sw_model config_alpha_closure false
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -126,6 +128,8 @@
var persistent real pv ( nVertLevels nCells Time ) 2 iro pv state - -
var persistent real vor ( nVertLevels nCells Time ) 2 iro vor state - -
var persistent real div ( nVertLevels nCells Time ) 2 iro div state - -
+var persistent real vorSmooth ( nVertLevels nCells Time ) 2 iro vorSmooth state - -
+var persistent real divSmooth ( nVertLevels nCells Time ) 2 iro divSmooth 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 - -
Deleted: branches/pv_based_swm/mpas/src/core_pvsw/module_alpha.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_alpha.F 2011-04-21 19:55:42 UTC (rev 800)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_alpha.F 2011-04-22 03:33:28 UTC (rev 801)
@@ -1,273 +0,0 @@
-module alpha_solver
-
- use grid_types
- use configure
- use constants
- use poisson_solver
- use dmpar
-
- implicit none
-
- real (KIND=RKIND), parameter :: threshold = 1.0e-8
- integer, parameter :: nIterPerGS = 2
-
- public :: init_alpha, gs_iteration_alpha, alpha
-
- contains
-
- subroutine init_alpha(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: put any init stuff in here
- !
- ! Input: grid meta data
- !
- ! Output:
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- type (block_type), pointer :: block
-
- block => domain % blocklist
- do while (associated(block))
-
- block => block % next
- enddo ! block
-
- end subroutine init_alpha
-
-
- subroutine alpha(domain)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: given the rough vorticity and divergence, find the smooth vorticity and divergence
- !
- ! Input: grid meta data
- !
- ! Output: vorSmooth, divSmooth
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
-
- integer :: iteration
- logical :: converge
- real (kind=RKIND) :: l2_vorSmooth_error, l2_divSmooth_error, globalVorSmoothError, globalDivSmoothError, globalArea, localArea, rflag, gflag
- type (block_type), pointer :: block
-
- iteration = 0; l2_vorSmooth_error = 1.0; l2_divSmooth_error = 1.0
- converge = .false.
-
- ! 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_alpha(block % state % time_levs(2) % state, block % mesh)
- block => block % next
- end do
-
- ! do a boundary update on vorSmooth and divSmooth
- block => domain % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorSmooth % 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 % divSmooth % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- endif
- block => block % next
- end do
-
-
- ! check the global error every so often
- if(mod(iteration,100).eq.0) then
- block => domain % blocklist
- do while (associated(block))
- call compute_alpha_error(block % state % time_levs(2) % state, block % mesh, l2_vorSmooth_error, l2_divSmooth_error, localArea)
- call dmpar_sum_real(domain % dminfo, l2_vorSmooth_error, globalVorSmoothError)
- call dmpar_sum_real(domain % dminfo, l2_divSmooth_error, globalDivSmoothError)
- call dmpar_sum_real(domain % dminfo, localArea, globalArea)
- block => block % next
- end do
- l2_vorSmooth_error = sqrt(globalVorSmoothError/globalArea)
- l2_divSmooth_error = sqrt(globalDivSmoothError/globalArea)
- if(l2_vorSmooth_error.lt.threshold.and.l2_divSmooth_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))
- 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_vorSmooth_error, l2_divSmooth_error, rflag, gflag
- 10 format(i8,6e20.10)
- endif
-
- enddo ! do while not converged
-
- ! copy vorSmooth and divSmooth into the other time level
- block => domain % blocklist
- do while (associated(block))
- block % state % time_levs(1) % state % vorSmooth % array(:,:) = block % state % time_levs(2) % state % vorSmooth % array(:,:)
- block % state % time_levs(1) % state % divSmooth % array(:,:) = block % state % time_levs(2) % state % divSmooth % array(:,:)
- block => block % next
- end do
-
- end subroutine alpha
-
-
- subroutine gs_iteration_alpha(state, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: do one or more iteration(s) using Guass Seidel to solve
- ! vorSmooth and divSmooth where del2(vorSmooth)=vor and del2(divSmooth)=div
- !
- ! Input: vor, div, initial guess for vorSmooth and divSmooth
- !
- ! Output: update vorSmooth and divSmooth
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(inout) :: state
- type (mesh_type), intent(in) :: grid
-
- integer :: nCells
- integer, dimension(:), pointer :: nEdgesOnCell
- integer, dimension(:,:), pointer :: cellsOnCell
- integer :: iCell, jCell, iter, k, j, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: areaCell
- real (kind=RKIND), dimension(:,:), pointer :: vor, div, vorSmooth, divSmooth
- 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
- poisson_weights => grid % poisson_weights % array
- nCells = grid % nCells
- nVertLevels = grid % nVertLevels
-
- vor => state % vor % array
- div => state % div % array
- vorSmooth => state % vorSmooth % array
- divSmooth => state % divSmooth % array
-
- ! do nIter GS iterations
- do iter=1,nIterPerGS
-
- ! loop over the vertical
- do k = 1, nVertLevels
-
- ! loop over cell centers
- do iCell=1,nCells
- sum1=0.0; sum2=0.0; sum3=0.0
-
- ! loop over neighboring cells
- do j=1,nEdgesOnCell(iCell)
- jCell = cellsOnCell(j,iCell)
-
- sum1 = sum1 + poisson_weights(j,iCell)
- sum2 = sum2 + poisson_weights(j,iCell) * vorSmooth(k,jCell)
- sum3 = sum3 + poisson_weights(j,iCell) * divSmooth(k,jCell)
-
- enddo ! loop over edges
-
- sum1 = sum1 * alphasq / areaCell(iCell)
- sum2 = sum2 * alphasq / areaCell(iCell)
- sum3 = sum3 * alphasq / areaCell(iCell)
-
- vorSmooth(k,iCell) = (1.0-sor)*vorSmooth(k,iCell) + sor*(sum2 + vor(k,iCell))/(1.0+sum1)
- divSmooth(k,iCell) = (1.0-sor)*divSmooth(k,iCell) + sor*(sum3 + div(k,iCell))/(1.0+sum1)
-
- enddo ! loop over cell centers
- enddo ! loop over vertical levels
-
- enddo ! loop over iterations
-
- end subroutine gs_iteration_alpha
-
-
- subroutine compute_alpha_error(state, grid, l2_vorSmooth_error, l2_divSmooth_error, localArea)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: compute error in poisson solver
- !
- ! Input: vorSmooth, divSmooth, vor and div
- !
- ! Output: converge (true or false)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (state_type), intent(in) :: state
- type (mesh_type), intent(in) :: grid
- real (kind=RKIND), intent(out) :: l2_vorSmooth_error, l2_divSmooth_error, localArea
-
- integer :: nCellsSolve
- integer, dimension(:), pointer :: nEdgesOnCell
- integer, dimension(:,:), pointer :: cellsOnCell
- integer :: iCell, jCell, iter, k, j, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: areaCell
- real (kind=RKIND), dimension(:,:), pointer :: vor, div, vorSmooth, divSmooth
- real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
- real (kind=RKIND) :: sum_error_vorSmooth, sum_error_divSmooth, sum_area, r1, r2
-
- poisson_weights => grid % poisson_weights % array
- cellsOnCell => grid % cellsOnCell % array
- nEdgesOnCell => grid % nEdgesOnCell % array
- areaCell => grid % areaCell % array
- nCellsSolve = grid % nCellsSolve
- nVertLevels = grid % nVertLevels
-
- vor => state % vor % array
- div => state % div % array
- vorSmooth => state % vorSmooth % array
- divSmooth => state % divSmooth % array
-
- do k=1,nVertLevels
-
- sum_error_vorSmooth = 0.0; sum_error_divSmooth = 0.0; sum_area = 0.0
- do iCell=1,nCellsSolve
-
- r1=0.0; r2=0.0
- do j=1,nEdgesOnCell(iCell)
-
- jCell = cellsOnCell(j,iCell)
- r1 = r1 + poisson_weights(j,iCell) * (vorSmooth(k,jCell)-vorSmooth(k,iCell))
- r2 = r2 + poisson_weights(j,iCell) * (divSmooth(k,jCell)-divSmooth(k,iCell))
-
- enddo
- r1 = vorSmooth(k,iCell) - r1 * alphasq / areaCell(iCell)
- r2 = divSmooth(k,iCell) - r2 * alphasq / areaCell(iCell)
-
- sum_error_vorSmooth = sum_error_vorSmooth + areaCell(iCell)*(r1 - vor (k,iCell))**2
- sum_error_divSmooth = sum_error_divSmooth + areaCell(iCell)*(r2 - div(k,iCell))**2
- sum_area = sum_area + areaCell(iCell)
-
- enddo
-
- l2_vorSmooth_error = sum_error_vorSmooth
- l2_divSmooth_error = sum_error_divSmooth
- localArea = sum_area
-
- enddo
-
- end subroutine compute_alpha_error
-
-end module alpha_solver
Added: branches/pv_based_swm/mpas/src/core_pvsw/module_alpha_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_alpha_solver.F (rev 0)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_alpha_solver.F 2011-04-22 03:33:28 UTC (rev 801)
@@ -0,0 +1,340 @@
+module alpha_solver
+
+ use grid_types
+ use configure
+ use constants
+ use poisson_solver
+ use dmpar
+
+ implicit none
+
+ real (KIND=RKIND), parameter :: threshold_alpha = 1.0e-8
+ integer, parameter :: nIterPerGS = 2
+
+ private
+ public :: init_alpha, alpha
+
+ contains
+
+ subroutine init_alpha(domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: put any init stuff in here
+ !
+ ! Input: grid meta data
+ !
+ ! Output:
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ type (block_type), pointer :: block
+
+ block => domain % blocklist
+ do while (associated(block))
+
+ block => block % next
+ enddo ! block
+
+ end subroutine init_alpha
+
+
+ subroutine alpha(domain)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: given the rough vorticity and divergence, find the smooth vorticity and divergence
+ !
+ ! Input: grid meta data
+ !
+ ! Output: vorSmooth, divSmooth
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+
+ integer :: iteration
+ logical :: converge
+ real (kind=RKIND) :: l2_vorSmooth_error, l2_divSmooth_error, globalVorSmoothError, globalDivSmoothError, globalArea, localArea, rflag, gflag
+ type (block_type), pointer :: block
+
+ iteration = 0; l2_vorSmooth_error = 1.0; l2_divSmooth_error = 1.0
+ converge = .false.
+
+ ! 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_vorSmooth_error, l2_divSmooth_error, localArea)
+ call dmpar_sum_real(domain % dminfo, l2_vorSmooth_error, globalVorSmoothError)
+ call dmpar_sum_real(domain % dminfo, l2_divSmooth_error, globalDivSmoothError)
+ call dmpar_sum_real(domain % dminfo, localArea, globalArea)
+ block => block % next
+ end do
+ l2_vorSmooth_error = globalVorSmoothError/globalArea
+ l2_divSmooth_error = globalDivSmoothError/globalArea
+ write(6,*) ' subtracting mean ', l2_vorSmooth_error, l2_divSmooth_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_vorSmooth_error
+ if(.not.config_bve) block % state % time_levs(2) % state % div % array(:,:) = block % state % time_levs(2) % state % div % array(:,:) - l2_divSmooth_error
+ block => block % next
+ end do
+
+ ! 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)
+ block => block % next
+ end do
+
+ ! do a boundary update on vorSmooth and divSmooth
+ block => domain % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % vorSmooth % 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 % divSmooth % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ endif
+ block => block % next
+ end do
+
+
+ ! check the global error every so often
+ if(mod(iteration,100).eq.0) then
+ block => domain % blocklist
+ do while (associated(block))
+ call compute_alpha_error(block % state % time_levs(2) % state, block % mesh, l2_vorSmooth_error, l2_divSmooth_error, localArea)
+ call dmpar_sum_real(domain % dminfo, l2_vorSmooth_error, globalVorSmoothError)
+ call dmpar_sum_real(domain % dminfo, l2_divSmooth_error, globalDivSmoothError)
+ call dmpar_sum_real(domain % dminfo, localArea, globalArea)
+ block => block % next
+ end do
+ l2_vorSmooth_error = sqrt(globalVorSmoothError/globalArea)
+ l2_divSmooth_error = sqrt(globalDivSmoothError/globalArea)
+ if(l2_vorSmooth_error.lt.threshold_alpha.and.l2_divSmooth_error.lt.threshold_alpha) 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))
+ 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) 'alpha ', iteration, l2_vorSmooth_error, l2_divSmooth_error, rflag, gflag
+ 10 format(a10,i8,6e20.10)
+ endif
+
+ enddo ! do while not converged
+
+ ! copy vorSmooth and divSmooth into the other time level
+ block => domain % blocklist
+ do while (associated(block))
+ block % state % time_levs(1) % state % vorSmooth % array(:,:) = block % state % time_levs(2) % state % vorSmooth % array(:,:)
+ block % state % time_levs(1) % state % divSmooth % array(:,:) = block % state % time_levs(2) % state % divSmooth % array(:,:)
+ block => block % next
+ end do
+
+ end subroutine alpha
+
+
+ subroutine gs_iteration(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: do one or more iteration(s) using Guass Seidel to solve
+ ! vorSmooth and divSmooth where del2(vorSmooth)=vor and del2(divSmooth)=div
+ !
+ ! Input: vor, div, initial guess for vorSmooth and divSmooth
+ !
+ ! Output: update vorSmooth and divSmooth
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ integer :: nCells
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnCell
+ integer :: iCell, jCell, iter, k, j, nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: vor, div, vorSmooth, divSmooth
+ 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
+ poisson_weights => grid % poisson_weights % array
+ nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
+
+ vor => state % vor % array
+ div => state % div % array
+ vorSmooth => state % vorSmooth % array
+ divSmooth => state % divSmooth % array
+
+ ! do nIter GS iterations
+ do iter=1,nIterPerGS
+
+ ! loop over the vertical
+ do k = 1, nVertLevels
+
+ ! loop over cell centers
+ do iCell=1,nCells
+ sum1=0.0; sum2=0.0; sum3=0.0
+
+ ! loop over neighboring cells
+ do j=1,nEdgesOnCell(iCell)
+ jCell = cellsOnCell(j,iCell)
+
+ sum1 = sum1 + poisson_weights(j,iCell)
+ sum2 = sum2 + poisson_weights(j,iCell) * vorSmooth(k,jCell)
+ sum3 = sum3 + poisson_weights(j,iCell) * divSmooth(k,jCell)
+
+ enddo ! loop over edges
+
+ sum1 = sum1 * (config_alpha)**2 / areaCell(iCell)
+ sum2 = sum2 * (config_alpha)**2 / areaCell(iCell)
+ sum3 = sum3 * (config_alpha)**2 / areaCell(iCell)
+
+ vorSmooth(k,iCell) = (1.0-sor)*vorSmooth(k,iCell) + sor*(sum2 + vor(k,iCell))/(1.0+sum1)
+ divSmooth(k,iCell) = (1.0-sor)*divSmooth(k,iCell) + sor*(sum3 + div(k,iCell))/(1.0+sum1)
+
+ enddo ! loop over cell centers
+ enddo ! loop over vertical levels
+
+ enddo ! loop over iterations
+
+ end subroutine gs_iteration
+
+
+ subroutine compute_alpha_error(state, grid, l2_vorSmooth_error, l2_divSmooth_error, localArea)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: compute error in poisson solver
+ !
+ ! Input: vorSmooth, divSmooth, vor and div
+ !
+ ! Output: converge (true or false)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(in) :: state
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), intent(out) :: l2_vorSmooth_error, l2_divSmooth_error, localArea
+
+ integer :: nCellsSolve
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer, dimension(:,:), pointer :: cellsOnCell
+ integer :: iCell, jCell, iter, k, j, nVertLevels
+ real (kind=RKIND), dimension(:), pointer :: areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: vor, div, vorSmooth, divSmooth
+ real (kind=RKIND), dimension(:,:), pointer :: poisson_weights
+ real (kind=RKIND) :: sum_error_vorSmooth, sum_error_divSmooth, sum_area, r1, r2
+
+ poisson_weights => grid % poisson_weights % array
+ cellsOnCell => grid % cellsOnCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ areaCell => grid % areaCell % array
+ nCellsSolve = grid % nCellsSolve
+ nVertLevels = grid % nVertLevels
+
+ vor => state % vor % array
+ div => state % div % array
+ vorSmooth => state % vorSmooth % array
+ divSmooth => state % divSmooth % array
+
+ do k=1,nVertLevels
+
+ sum_error_vorSmooth = 0.0; sum_error_divSmooth = 0.0; sum_area = 0.0
+ do iCell=1,nCellsSolve
+
+ r1=0.0; r2=0.0
+ do j=1,nEdgesOnCell(iCell)
+
+ jCell = cellsOnCell(j,iCell)
+ r1 = r1 + poisson_weights(j,iCell) * (vorSmooth(k,jCell)-vorSmooth(k,iCell))
+ r2 = r2 + poisson_weights(j,iCell) * (divSmooth(k,jCell)-divSmooth(k,iCell))
+
+ enddo
+ r1 = vorSmooth(k,iCell) - r1 * (config_alpha)**2 / areaCell(iCell)
+ r2 = divSmooth(k,iCell) - r2 * (config_alpha)**2 / areaCell(iCell)
+
+ sum_error_vorSmooth = sum_error_vorSmooth + areaCell(iCell)*(r1 - vor (k,iCell))**2
+ sum_error_divSmooth = sum_error_divSmooth + areaCell(iCell)*(r2 - div(k,iCell))**2
+ sum_area = sum_area + areaCell(iCell)
+
+ enddo
+
+ l2_vorSmooth_error = sum_error_vorSmooth
+ l2_divSmooth_error = sum_error_divSmooth
+ localArea = sum_area
+
+ enddo
+
+ end subroutine compute_alpha_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
+
+
+end module alpha_solver
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-21 19:55:42 UTC (rev 800)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_mpas_core.F 2011-04-22 03:33:28 UTC (rev 801)
@@ -15,6 +15,7 @@
use grid_types
use test_cases
use poisson_solver
+ use alpha_solver
implicit none
@@ -40,6 +41,11 @@
write(6,*) 'calling init_poisson'
call init_poisson(domain)
+ if(config_alpha_closure) then
+ write(6,*) 'calling alpha'
+ call alpha(domain)
+ endif
+
write(6,*) 'calling poisson'
call poisson(domain)
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-21 19:55:42 UTC (rev 800)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-04-22 03:33:28 UTC (rev 801)
@@ -10,7 +10,8 @@
real (KIND=RKIND), parameter :: threshold = 1.0e-8
integer, parameter :: nIterPerGS = 2
- public :: init_poisson, gs_iteration, poisson
+ private
+ public :: init_poisson, poisson, derive_velocity
contains
@@ -133,8 +134,13 @@
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
+ if(config_alpha_closure) then
+ block % state % time_levs(2) % state % vorSmooth % array(:,:) = block % state % time_levs(2) % state % vorSmooth % array(:,:) - l2_psi_error
+ if(.not.config_bve) block % state % time_levs(2) % state % divSmooth % array(:,:) = block % state % time_levs(2) % state % divSmooth % array(:,:) - l2_chi_error
+ else
+ 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
+ endif
block => block % next
end do
@@ -241,8 +247,14 @@
nCells = grid % nCells
nVertLevels = grid % nVertLevels
- vor => state % vor % array
- div => state % div % array
+ if(config_alpha_closure) then
+ vor => state % vorSmooth % array
+ div => state % divSmooth % array
+ else
+ vor => state % vor % array
+ div => state % div % array
+ endif
+
psi => state % psi % array
chi => state % chi % array
@@ -308,8 +320,14 @@
nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
- vor => state % vor % array
- div => state % div % array
+ if(config_alpha_closure) then
+ vor => state % vorSmooth % array
+ div => state % divSmooth % array
+ else
+ vor => state % vor % array
+ div => state % div % array
+ endif
+
psi => state % psi % array
chi => state % chi % array
@@ -368,8 +386,13 @@
nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
- vor => state % vor % array
- div => state % div % array
+ if(config_alpha_closure) then
+ vor => state % vorSmooth % array
+ div => state % divSmooth % array
+ else
+ vor => state % vor % array
+ div => state % div % array
+ endif
do k=1,nVertLevels
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-21 19:55:42 UTC (rev 800)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-04-22 03:33:28 UTC (rev 801)
@@ -5,6 +5,7 @@
use configure
use constants
use poisson_solver
+ use alpha_solver
use dmpar
@@ -141,7 +142,7 @@
! block => block % next
! end do
-! --- invert elliptic equations for psi and chi and then derive normal and tangent velocity
+! --- get div and vor
block => domain % blocklist
do while (associated(block))
@@ -151,6 +152,14 @@
enddo
block => block % next
end do
+
+! --- compute smoothed vorticity for alpha closure
+ if(config_alpha_closure) then
+ call alpha(domain)
+ endif
+
+! --- invert elliptic equations for psi and chi and then derive normal and tangent velocity
+
call poisson(domain)
call derive_velocity(domain)
block => domain % blocklist
</font>
</pre>