<p><b>ringler@lanl.gov</b> 2011-06-10 12:44:09 -0600 (Fri, 10 Jun 2011)</p><p><br>
add drag term to vorticity 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-06-10 17:15:11 UTC (rev 895)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-06-10 18:44:09 UTC (rev 896)
@@ -21,6 +21,7 @@
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 real sw_model config_dragCoefficient 0.001
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -148,6 +149,7 @@
# Diagnostic fields: only written to output
var persistent real v ( nVertLevels nEdges Time ) 2 - v state - -
+var persistent real gradKE ( nVertLevels nEdges Time ) 2 o gradKE state - -
var persistent real divergence ( nVertLevels nCells Time ) 2 - divergence state - -
var persistent real h_edge ( nVertLevels nEdges Time ) 2 - h_edge state - -
var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
@@ -156,9 +158,13 @@
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 - uReconstructZ state - -
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
+var persistent real gradKeReconstructX ( nVertLevels nCells Time ) 2 o gradKeReconstructX state - -
+var persistent real gradKeReconstructY ( nVertLevels nCells Time ) 2 o gradKeReconstructY state - -
+var persistent real gradKeReconstructZ ( nVertLevels nCells Time ) 2 o gradKeReconstructZ state - -
+var persistent real gradKeReconstructZonal ( nVertLevels nCells Time ) 2 o gradKeReconstructZonal state - -
+var persistent real gradKeReconstructMerid ( nVertLevels nCells Time ) 2 o gradKeReconstructMerid state - -
var persistent real dvordt ( nVertLevels nCells Time ) 2 o dvordt state - -
var persistent real pv_advection ( nVertLevels nCells Time ) 2 o pv_advection state - -
var persistent real pv_dissipation ( nVertLevels nCells Time ) 2 o pv_dissipation state - -
var persistent real pv_forcing ( nVertLevels nCells Time ) 2 o pv_forcing state - -
var persistent real pv_drag ( nVertLevels nCells Time ) 2 o pv_drag state - -
-
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-06-10 17:15:11 UTC (rev 895)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-06-10 18:44:09 UTC (rev 896)
@@ -225,6 +225,7 @@
end do
end do
call compute_solve_diagnostics(dt, provis, block % mesh)
+ call compute_drag(provis, block % mesh)
block => block % next
end do
end if
@@ -310,8 +311,10 @@
do while (associated(block))
call compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
+ call compute_drag(provis, block % mesh)
call reconstruct(block % state % time_levs(2) % state, block % mesh)
+ call reconstructGradKE(block % state % time_levs(2) % state, block % mesh)
block => block % next
end do
@@ -731,6 +734,75 @@
end subroutine compute_scalar_tend
+ subroutine compute_drag(s, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! drag on RHS of PV equation
+ !
+ ! Input: grid - grid metadata
+ !
+ ! Output: s - computed diagnostics
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: s
+ type (mesh_type), intent(in) :: grid
+
+ real (kind=RKIND), pointer, dimension(:) :: xCell, yCell, zCell
+ real (kind=RKIND), pointer, dimension(:,:) :: pv_drag, vor, ke, h
+ real (kind=RKIND), pointer, dimension(:,:) :: ux, uy, uz, gx, gy, gz
+ real (kind=RKIND) :: p(3), q(3), t(3), r
+ integer :: iCell, k
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+
+ pv_drag => s % pv_drag % array
+ vor => s % vor % array
+ ke => s % ke % array
+ h => s % h % array
+
+ ux => s % uReconstructX % array
+ uy => s % uReconstructY % array
+ uz => s % uReconstructZ % array
+
+ gx => s % gradKeReconstructX % array
+ gy => s % gradKeReconstructY % array
+ gz => s % gradKeReconstructZ % array
+
+
+ do iCell = 1,grid%nCells
+ do k=1,grid % nVertLevels
+ pv_drag(k,iCell) = -config_dragCoefficient * sqrt(2.0*ke(k,iCell)) * vor(k,iCell) / h(k,iCell)
+ enddo
+ enddo
+
+ do iCell = 1,grid%nCells
+ do k=1,grid % nVertLevels
+ p(1) = gx(k,iCell); p(2) = gy(k,iCell); p(3) = gz(k,iCell)
+ q(1) = ux(k,iCell); q(2) = uy(k,iCell); q(3) = uz(k,iCell)
+ call cross_product_in_R3(p,q,t)
+ r = xCell(iCell)*t(1) + yCell(iCell)*t(2) + zCell(iCell)*t(3)
+ pv_drag(k,iCell) = pv_drag(k,iCell) - config_dragCoefficient * r / h(k,iCell)
+ enddo
+ enddo
+
+ contains
+
+ subroutine cross_product_in_R3(p_1,p_2,p_out)
+ real (kind=RKIND), intent(in) :: p_1 (3), p_2 (3)
+ real (kind=RKIND), intent(out) :: p_out (3)
+
+ p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
+ p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
+ p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
+ end subroutine cross_product_in_R3
+
+
+ end subroutine compute_drag
+
+
subroutine compute_solve_diagnostics(dt, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
@@ -753,7 +825,7 @@
integer :: nCells, nEdges, nVertices, nVertLevels
real (kind=RKIND), dimension(:), pointer :: h_s, dvEdge, dcEdge, areaCell
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, h_edge, h, u, v, tend_h, tend_u, &
- ke, divergence
+ ke, divergence, gradKE
integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgesOnEdge, boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r
@@ -769,6 +841,7 @@
tend_u => s % u % array
divergence => s % divergence % array
ke => s % ke % array
+ gradKE => s % gradKE % array
weightsOnEdge => grid % weightsOnEdge % array
cellsOnEdge => grid % cellsOnEdge % array
@@ -849,6 +922,18 @@
end do
!
+ ! Compute kinetic energy at every edge
+ !
+ gradKE(:,:) = 0.0
+ do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,nVertLevels
+ gradKE(k,iEdge) = (ke(k,cell2) - ke(k,cell1)) / dcEdge(iEdge)
+ enddo
+ enddo
+
+ !
! Compute v (tangential) velocities
!
v(:,:) = 0.0
Modified: branches/pv_based_swm/mpas/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/pv_based_swm/mpas/src/operators/module_vector_reconstruction.F 2011-06-10 17:15:11 UTC (rev 895)
+++ branches/pv_based_swm/mpas/src/operators/module_vector_reconstruction.F 2011-06-10 18:44:09 UTC (rev 896)
@@ -7,7 +7,7 @@
implicit none
- public :: init_reconstruct, reconstruct
+ public :: init_reconstruct, reconstruct, reconstructGradKE
contains
@@ -200,4 +200,93 @@
end subroutine reconstruct
+ subroutine reconstructGradKE(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: reconstruct gradient of KE vector field at cell centers based on radial basis functions
+ !
+ ! Input: grid meta data and vector component data residing at cell edges
+ !
+ ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (state_type), intent(inout) :: state
+ type (mesh_type), intent(in) :: grid
+
+ ! temporary arrays needed in the compute procedure
+ integer :: nCellsSolve
+ integer, dimension(:,:), pointer :: edgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer :: iCell,iEdge, i
+ real (kind=RKIND), dimension(:,:), pointer :: u
+ real(kind=RKIND), dimension(:), pointer :: latCell, lonCell
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructZonal, uReconstructMeridional
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+
+ logical :: on_a_sphere
+
+ real (kind=RKIND) :: clat, slat, clon, slon
+
+
+ ! stored arrays used during compute procedure
+ coeffs_reconstruct => grid % coeffs_reconstruct % array
+
+ ! temporary variables
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ nCellsSolve = grid % nCellsSolve
+ u => state % u % array
+ uReconstructX => state % gradKeReconstructX % array
+ uReconstructY => state % gradKeReconstructY % array
+ uReconstructZ => state % gradKeReconstructZ % array
+
+ latCell => grid % latCell % array
+ lonCell => grid % lonCell % array
+ uReconstructZonal => state % gradKeReconstructZonal % array
+ uReconstructMeridional => state % gradKeReconstructMerid % array
+ on_a_sphere = grid % on_a_sphere
+
+ ! init the intent(out)
+ uReconstructX = 0.0
+ uReconstructY = 0.0
+ uReconstructZ = 0.0
+
+ ! loop over cell centers
+ do iCell=1,nCellsSolve
+ ! a more efficient reconstruction where rbf_values*matrix_reconstruct has been precomputed
+ ! in coeffs_reconstruct
+ do i=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i,iCell)
+ uReconstructX(:,iCell) = uReconstructX(:,iCell) &
+ + coeffs_reconstruct(1,i,iCell) * u(:,iEdge)
+ uReconstructY(:,iCell) = uReconstructY(:,iCell) &
+ + coeffs_reconstruct(2,i,iCell) * u(:,iEdge)
+ uReconstructZ(:,iCell) = uReconstructZ(:,iCell) &
+ + coeffs_reconstruct(3,i,iCell) * u(:,iEdge)
+
+ enddo
+ enddo ! iCell
+
+ if(on_a_sphere) then
+ do iCell=1,nCellsSolve
+ clat = cos(latCell(iCell))
+ slat = sin(latCell(iCell))
+ clon = cos(lonCell(iCell))
+ slon = sin(lonCell(iCell))
+ uReconstructZonal(:,iCell) = -uReconstructX(:,iCell)*slon + uReconstructY(:,iCell)*clon
+ uReconstructMeridional(:,iCell) = -(uReconstructX(:,iCell)*clon &
+ + uReconstructY(:,iCell)*slon)*slat &
+ + uReconstructZ(:,iCell)*clat
+ end do
+ else
+ uReconstructZonal = uReconstructX
+ uReconstructMeridional = uReconstructY
+ end if
+
+ end subroutine reconstructGradKE
+
+
end module vector_reconstruction
</font>
</pre>