diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index fb8cfe959a..afd92b94ad 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -545,11 +545,6 @@
-
-
-
-
-
#ifdef MPAS_CAM_DYCORE
diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F
index 35c4034815..df9a800292 100644
--- a/src/core_atmosphere/mpas_atm_core.F
+++ b/src/core_atmosphere/mpas_atm_core.F
@@ -419,6 +419,10 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
logical, pointer :: config_do_restart, config_do_DAcycling
+ logical, pointer :: on_a_sphere
+ real (kind=RKIND), pointer :: sphere_radius
+
+
call atm_compute_signs(mesh)
call mpas_pool_get_subpool(block % structs, 'diag', diag)
@@ -470,6 +474,10 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
!!!!! End compute inverses
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ call mpas_pool_get_config(mesh, 'on_a_sphere', on_a_sphere)
+ call mpas_pool_get_config(mesh, 'sphere_radius', sphere_radius)
+ call atm_initialize_deformation_weights(mesh, nCells, on_a_sphere, sphere_radius)
+
call atm_adv_coef_compression(mesh)
call atm_couple_coef_3rd_order(mesh, block % configs)
@@ -1600,5 +1608,238 @@ subroutine mpas_atm_run_compatibility(dminfo, blockList, streamManager, ierr)
end subroutine mpas_atm_run_compatibility
+
+ subroutine atm_initialize_deformation_weights(mesh, nCells, on_a_sphere, sphere_radius)
+
+!
+! compute the cell coefficients for the deformation calculations
+! WCS, 13 July 2010
+!
+
+ use mpas_vector_operations, only : mpas_fix_periodicity
+ use mpas_timer, only : mpas_timer_start, mpas_timer_stop
+ use mpas_geometry_utils, only : mpas_sphere_angle, mpas_plane_angle, mpas_arc_length
+
+ implicit none
+
+ type (mpas_pool_type), intent(inout) :: mesh
+ integer, intent(in) :: nCells
+ logical, intent(in) :: on_a_sphere
+ real (kind=RKIND), intent(in) :: sphere_radius
+
+! local variables
+
+ real (kind=RKIND), dimension(:,:), pointer :: deformation_coef_c2, deformation_coef_s2, deformation_coef_cs
+ real (kind=RKIND), dimension(:,:), pointer :: deformation_coef_c, deformation_coef_s
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, cellsOnCell, verticesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+ real (kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
+ real (kind=RKIND), dimension(:), pointer :: xEdge, yEdge, zEdge, angleEdge
+
+ real (kind=RKIND), dimension(nCells) :: theta_abs
+
+ real (kind=RKIND), dimension(25) :: xc, yc, zc ! cell center coordinates
+ real (kind=RKIND), dimension(25) :: thetav, thetat, dl_sphere
+ real (kind=RKIND) :: dl
+ integer :: i, ip1, ip2, n
+ integer :: iCell
+ real (kind=RKIND) :: pii
+ real (kind=RKIND), dimension(25) :: xp, yp
+ real (kind=RKIND) :: xe, ye
+
+ integer, dimension(25) :: cell_list
+
+ integer :: iv, ie
+ logical :: do_the_cell
+ real (kind=RKIND) :: area_cell, sint2, cost2, sint_cost, dx, dy
+
+ real(kind=RKIND), pointer :: x_period, y_period
+
+
+ call mpas_pool_get_config(mesh, 'x_period', x_period)
+ call mpas_pool_get_config(mesh, 'y_period', y_period)
+
+ call mpas_pool_get_array(mesh, 'deformation_coef_c2', deformation_coef_c2)
+ call mpas_pool_get_array(mesh, 'deformation_coef_s2', deformation_coef_s2)
+ call mpas_pool_get_array(mesh, 'deformation_coef_cs', deformation_coef_cs)
+ call mpas_pool_get_array(mesh, 'deformation_coef_c', deformation_coef_c)
+ call mpas_pool_get_array(mesh, 'deformation_coef_s', deformation_coef_s)
+ call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
+ call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
+ call mpas_pool_get_array(mesh, 'edgesOnCell', edgesOnCell)
+ call mpas_pool_get_array(mesh, 'cellsOnCell', cellsOnCell)
+ call mpas_pool_get_array(mesh, 'verticesOnCell', verticesOnCell)
+ call mpas_pool_get_array(mesh, 'xCell', xCell)
+ call mpas_pool_get_array(mesh, 'yCell', yCell)
+ call mpas_pool_get_array(mesh, 'zCell', zCell)
+ call mpas_pool_get_array(mesh, 'xVertex', xVertex)
+ call mpas_pool_get_array(mesh, 'yVertex', yVertex)
+ call mpas_pool_get_array(mesh, 'zVertex', zVertex)
+ call mpas_pool_get_array(mesh, 'xEdge', xEdge)
+ call mpas_pool_get_array(mesh, 'yEdge', yEdge)
+ call mpas_pool_get_array(mesh, 'zEdge', zEdge)
+ call mpas_pool_get_array(mesh, 'angleEdge', angleEdge)
+
+ deformation_coef_c2(:,:) = 0.
+ deformation_coef_s2(:,:) = 0.
+ deformation_coef_cs(:,:) = 0.
+ deformation_coef_c(:,:) = 0.
+ deformation_coef_s(:,:) = 0.
+
+ pii = 2.*asin(1.0)
+
+ do iCell = 1, nCells
+
+ cell_list(1) = iCell
+ do i=2,nEdgesOnCell(iCell)+1
+ cell_list(i) = cellsOnCell(i-1,iCell)
+ end do
+ n = nEdgesOnCell(iCell) + 1
+
+! check to see if we are reaching outside the halo
+
+ do_the_cell = .true.
+ do i=1,n
+ if (cell_list(i) > nCells) do_the_cell = .false.
+ end do
+
+
+ if (.not. do_the_cell) cycle
+
+ ! compute poynomial fit for this cell if all needed neighbors exist
+
+ if (on_a_sphere) then
+
+ ! xc holds the center point and the vertex points of the cell,
+ ! normalized to a sphere or radius 1.
+
+ xc(1) = xCell(iCell)/sphere_radius
+ yc(1) = yCell(iCell)/sphere_radius
+ zc(1) = zCell(iCell)/sphere_radius
+
+ do i=2,n
+ iv = verticesOnCell(i-1,iCell)
+ xc(i) = xVertex(iv)/sphere_radius
+ yc(i) = yVertex(iv)/sphere_radius
+ zc(i) = zVertex(iv)/sphere_radius
+ end do
+
+ !
+ ! In case the current cell center lies at exactly z=1.0, the sphere_angle() routine
+ ! may generate an FPE since the triangle it is given will have a zero side length
+ ! adjacent to the vertex whose angle we are trying to find; in this case, simply
+ ! set the value of theta_abs directly
+ !
+ if (zc(1) == 1.0) then
+ theta_abs(iCell) = pii/2.
+ else
+ ! theta_abs is the angle to the first vertex from the center, normalized so that
+ ! an eastward pointing vector has a angle of 0.
+ theta_abs(iCell) = pii/2. - mpas_sphere_angle( xc(1), yc(1), zc(1), &
+ xc(2), yc(2), zc(2), &
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND )
+ end if
+
+ ! here we are constructing the tangent-plane cell.
+ ! thetat is the angle in the (x,y) tangent-plane coordinate from
+ ! the cell center to each vertex, normalized so that an
+ ! eastward pointing vector has a angle of 0.
+
+ ! dl_sphere is the spherical distance from the cell center
+ ! to the sphere vertex points for the cell.
+
+ thetat(1) = theta_abs(iCell)
+ do i=1,n-1
+
+ ip2 = i+2
+ if (ip2 > n) ip2 = 2
+
+ thetav(i) = mpas_sphere_angle( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1), &
+ xc(ip2), yc(ip2), zc(ip2) )
+ dl_sphere(i) = sphere_radius*mpas_arc_length( xc(1), yc(1), zc(1), &
+ xc(i+1), yc(i+1), zc(i+1) )
+ if(i.gt.1) thetat(i) = thetat(i-1)+thetav(i-1)
+ end do
+
+ ! xp and yp are the tangent-plane vertex points with the cell center at (0,0)
+
+ do i=1,n-1
+ xp(i) = cos(thetat(i)) * dl_sphere(i)
+ yp(i) = sin(thetat(i)) * dl_sphere(i)
+ end do
+
+ else ! On an x-y plane
+
+ do i=1,n-1
+ iv = verticesOnCell(i,iCell)
+ xp(i) = mpas_fix_periodicity(xVertex(iv),xCell(iCell),x_period) - xCell(iCell)
+ yp(i) = mpas_fix_periodicity(yVertex(iv),yCell(iCell),y_period) - yCell(iCell)
+ end do
+
+ do i=1,n-1
+ ie = edgesOnCell(i,iCell)
+ xe = mpas_fix_periodicity(xEdge(ie),xCell(iCell),x_period) - xCell(iCell)
+ ye = mpas_fix_periodicity(yEdge(ie),yCell(iCell),y_period) - yCell(iCell)
+ thetat(i) = atan2(ye,xe)
+ end do
+
+ theta_abs(iCell) = thetat(1)
+
+ end if
+
+ ! (1) compute cell area on the tangent plane used in the integrals
+ ! (2) compute angle of cell edge normal vector. here we are repurposing thetat
+ thetat(1) = theta_abs(iCell)
+
+ do i=2,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ thetat(i) = mpas_plane_angle( 0.0_RKIND, 0.0_RKIND, 0.0_RKIND, &
+ xp(i)-xp(i-1), yp(i)-yp(i-1), 0.0_RKIND, &
+ xp(ip1)-xp(i), yp(ip1)-yp(i), 0.0_RKIND, &
+ 0.0_RKIND, 0.0_RKIND, 1.0_RKIND)
+ thetat(i) = thetat(i) + thetat(i-1)
+ end do
+
+ area_cell = 0.
+ do i=1,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ dx = xp(ip1)-xp(i)
+ dy = yp(ip1)-yp(i)
+ area_cell = area_cell + 0.25*(xp(i)+xp(ip1))*(yp(ip1)-yp(i)) - 0.25*(yp(i)+yp(ip1))*(xp(ip1)-xp(i))
+ thetat(i) = atan2(dy,dx)-pii/2.
+ end do
+
+ ! coefficients - see documentation for the formulas.
+
+ do i=1,n-1
+ ip1 = i+1
+ if (ip1 == n) ip1 = 1
+ dl = sqrt((xp(ip1)-xp(i))**2 + (yp(ip1)-yp(i))**2)
+ sint2 = (sin(thetat(i)))**2
+ cost2 = (cos(thetat(i)))**2
+ sint_cost = sin(thetat(i))*cos(thetat(i))
+ deformation_coef_c2(i,iCell) = dl*cost2/area_cell
+ deformation_coef_s2(i,iCell) = dl*sint2/area_cell
+ deformation_coef_cs(i,iCell) = dl*sint_cost/area_cell
+ deformation_coef_c(i,iCell) = dl*cos(thetat(i))/area_cell
+ deformation_coef_s(i,iCell) = dl*sin(thetat(i))/area_cell
+ if (cellsOnEdge(1,EdgesOnCell(i,iCell)) /= iCell) then
+ deformation_coef_c2(i,iCell) = - deformation_coef_c2(i,iCell)
+ deformation_coef_s2(i,iCell) = - deformation_coef_s2(i,iCell)
+ deformation_coef_cs(i,iCell) = - deformation_coef_cs(i,iCell)
+! deformation_coef_c(i,iCell) = - deformation_coef_c(i,iCell)
+! deformation_coef_s(i,iCell) = - deformation_coef_s(i,iCell)
+ end if
+
+ end do
+
+ end do
+
+ end subroutine atm_initialize_deformation_weights
+
end module atm_core
diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index 938b8dde37..0a80532158 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -478,11 +478,6 @@
-
-
-
-
-
@@ -588,11 +583,6 @@
-
-
-
-
-
@@ -1134,21 +1124,6 @@
-
-
-
-
-
-
-
-
-
-
diff --git a/src/core_init_atmosphere/mpas_atm_advection.F b/src/core_init_atmosphere/mpas_atm_advection.F
index 4852b17117..f4d44c984e 100644
--- a/src/core_init_atmosphere/mpas_atm_advection.F
+++ b/src/core_init_atmosphere/mpas_atm_advection.F
@@ -11,7 +11,6 @@ module atm_advection
use mpas_derived_types
use mpas_pool_routines
use mpas_constants
- use mpas_vector_operations
use mpas_abort, only : mpas_dmpar_global_abort
use mpas_log, only : mpas_log_write
@@ -759,14 +758,10 @@ subroutine atm_initialize_deformation_weights( mesh, nCells, on_a_sphere, sphere
real (kind=RKIND), dimension(:,:), pointer :: defc_a, defc_b
real (kind=RKIND), dimension(:,:), pointer :: cell_gradient_coef_x, cell_gradient_coef_y
- real (kind=RKIND), dimension(:,:), pointer :: deformation_coef_c2, deformation_coef_s2, deformation_coef_cs
- real (kind=RKIND), dimension(:,:), pointer :: deformation_coef_c, deformation_coef_s
integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, cellsOnCell, verticesOnCell
integer, dimension(:), pointer :: nEdgesOnCell
real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
real (kind=RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
- real (kind=RKIND), dimension(:), pointer :: xEdge, yEdge, zEdge, angleEdge
- real (kind=RKIND), dimension(:), pointer :: areaCell
real (kind=RKIND), dimension(nCells) :: theta_abs
@@ -777,32 +772,19 @@ subroutine atm_initialize_deformation_weights( mesh, nCells, on_a_sphere, sphere
integer :: iCell
real (kind=RKIND) :: pii
real (kind=RKIND), dimension(25) :: xp, yp
- real (kind=RKIND) :: xe, ye
real (kind=RKIND) :: length_scale
integer, dimension(25) :: cell_list
- integer :: iv, ie
+ integer :: iv
logical :: do_the_cell
real (kind=RKIND) :: area_cell, sint2, cost2, sint_cost, dx, dy
- logical, pointer :: is_periodic
- real(kind=RKIND), pointer :: x_period, y_period
-
-
- call mpas_pool_get_config(mesh, 'is_periodic', is_periodic)
- call mpas_pool_get_config(mesh, 'x_period', x_period)
- call mpas_pool_get_config(mesh, 'y_period', y_period)
call mpas_pool_get_array(mesh, 'defc_a', defc_a)
call mpas_pool_get_array(mesh, 'defc_b', defc_b)
call mpas_pool_get_array(mesh, 'cell_gradient_coef_x', cell_gradient_coef_x)
call mpas_pool_get_array(mesh, 'cell_gradient_coef_y', cell_gradient_coef_y)
- call mpas_pool_get_array(mesh, 'deformation_coef_c2', deformation_coef_c2)
- call mpas_pool_get_array(mesh, 'deformation_coef_s2', deformation_coef_s2)
- call mpas_pool_get_array(mesh, 'deformation_coef_cs', deformation_coef_cs)
- call mpas_pool_get_array(mesh, 'deformation_coef_c', deformation_coef_c)
- call mpas_pool_get_array(mesh, 'deformation_coef_s', deformation_coef_s)
call mpas_pool_get_array(mesh, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(mesh, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(mesh, 'edgesOnCell', edgesOnCell)
@@ -814,19 +796,9 @@ subroutine atm_initialize_deformation_weights( mesh, nCells, on_a_sphere, sphere
call mpas_pool_get_array(mesh, 'xVertex', xVertex)
call mpas_pool_get_array(mesh, 'yVertex', yVertex)
call mpas_pool_get_array(mesh, 'zVertex', zVertex)
- call mpas_pool_get_array(mesh, 'xEdge', xEdge)
- call mpas_pool_get_array(mesh, 'yEdge', yEdge)
- call mpas_pool_get_array(mesh, 'zEdge', zEdge)
- call mpas_pool_get_array(mesh, 'angleEdge', angleEdge)
- call mpas_pool_get_array(mesh, 'areaCell', areaCell)
defc_a(:,:) = 0.
defc_b(:,:) = 0.
- deformation_coef_c2(:,:) = 0.
- deformation_coef_s2(:,:) = 0.
- deformation_coef_cs(:,:) = 0.
- deformation_coef_c(:,:) = 0.
- deformation_coef_s(:,:) = 0.
cell_gradient_coef_x(:,:) = 0.
cell_gradient_coef_y(:,:) = 0.
@@ -916,58 +888,21 @@ subroutine atm_initialize_deformation_weights( mesh, nCells, on_a_sphere, sphere
else ! On an x-y plane
- do i=1,n-1
- iv = verticesOnCell(i,iCell)
- xp(i) = mpas_fix_periodicity(xVertex(iv),xCell(iCell),x_period) - xCell(iCell)
- yp(i) = mpas_fix_periodicity(yVertex(iv),yCell(iCell),y_period) - yCell(iCell)
- end do
+ theta_abs(iCell) = 0.0
- ! if(iCell.lt.11) then
- ! call mpas_log_write(' setting defc coefs, cell $i', intArgs=(/iCell/))
- ! do i=1,n-1
- ! iv = verticesOnCell(i,iCell)
- ! call mpas_log_write(' xp,yp,xvc,yvc, $r $r $r $r', realArgs=(/xp(i),yp(i),xVertex(iv)-xCell(iCell),yVertex(iv)-yCell(iCell)/))
- ! end do
- ! end if
+ xp(1) = xCell(iCell)
+ yp(1) = yCell(iCell)
- do i=1,n-1
- ie = edgesOnCell(i,iCell)
- xe = mpas_fix_periodicity(xEdge(ie),xCell(iCell),x_period) - xCell(iCell)
- ye = mpas_fix_periodicity(yEdge(ie),yCell(iCell),y_period) - yCell(iCell)
- thetat(i) = atan2(ye,xe)
+ do i=2,n
+ iv = verticesOnCell(i-1,iCell)
+ xp(i) = xVertex(iv)
+ yp(i) = yVertex(iv)
end do
- ! if(iCell .lt. 11) then
- ! call mpas_log_write(' edge angles, plane calc, cell $i', intArgs=(/iCell/))
- ! do i=1,n-1
- ! call mpas_log_write(' edge angle $r', realArgs=(/thetat(i)*180./3.1415926/))
- ! end do
- ! end if
-
- theta_abs(iCell) = thetat(1)
-
end if
! (1) compute cell area on the tangent plane used in the integrals
! (2) compute angle of cell edge normal vector. here we are repurposing thetat
- thetat(1) = theta_abs(iCell)
-
- do i=2,n-1
- ip1 = i+1
- if (ip1 == n) ip1 = 1
- thetat(i) = plane_angle( 0.0_RKIND, 0.0_RKIND, 0.0_RKIND, &
- xp(i)-xp(i-1), yp(i)-yp(i-1), 0.0_RKIND, &
- xp(ip1)-xp(i), yp(ip1)-yp(i), 0.0_RKIND, &
- 0.0_RKIND, 0.0_RKIND, 1.0_RKIND)
- thetat(i) = thetat(i) + thetat(i-1)
- end do
-
- ! if(iCell .lt. 11) then
- ! call mpas_log_write(' edge angles, generic calc, cell $i', intArgs=(/iCell/))
- ! do i=1,n-1
- ! call mpas_log_write(' edge angle $r', realArgs=(/thetat(i)*180./3.1415926/))
- ! end do
- ! end if
area_cell = 0.
do i=1,n-1
@@ -992,241 +927,15 @@ subroutine atm_initialize_deformation_weights( mesh, nCells, on_a_sphere, sphere
defc_b(i,iCell) = dl*2.*sint_cost/area_cell
cell_gradient_coef_x(i,iCell) = dl*cos(thetat(i))/area_cell
cell_gradient_coef_y(i,iCell) = dl*sin(thetat(i))/area_cell
- deformation_coef_c2(i,iCell) = dl*cost2/area_cell
- deformation_coef_s2(i,iCell) = dl*sint2/area_cell
- deformation_coef_cs(i,iCell) = dl*sint_cost/area_cell
- deformation_coef_c(i,iCell) = dl*cos(thetat(i))/area_cell
- deformation_coef_s(i,iCell) = dl*sin(thetat(i))/area_cell
if (cellsOnEdge(1,EdgesOnCell(i,iCell)) /= iCell) then
defc_a(i,iCell) = - defc_a(i,iCell)
defc_b(i,iCell) = - defc_b(i,iCell)
- deformation_coef_c2(i,iCell) = - deformation_coef_c2(i,iCell)
- deformation_coef_s2(i,iCell) = - deformation_coef_s2(i,iCell)
- deformation_coef_cs(i,iCell) = - deformation_coef_cs(i,iCell)
-! deformation_coef_c(i,iCell) = - deformation_coef_c(i,iCell)
-! deformation_coef_s(i,iCell) = - deformation_coef_s(i,iCell)
end if
end do
end do
- call atm_init_test_coefs( deformation_coef_c2, deformation_coef_s2, &
- deformation_coef_cs, deformation_coef_c, &
- deformation_coef_s, &
- is_periodic, on_a_sphere, &
- x_period, y_period, &
- xEdge, yEdge, zEdge, &
- xCell, yCell, zCell, nCells, &
- angleEdge, nEdgesOnCell, edgesOnCell )
-
-
end subroutine atm_initialize_deformation_weights
-
- subroutine atm_init_test_coefs( deformation_coef_c2, deformation_coef_s2, &
- deformation_coef_cs, deformation_coef_c, &
- deformation_coef_s, &
- is_periodic, on_a_sphere, &
- x_period, y_period, &
- xEdge, yEdge, zEdge, &
- xCell, yCell, zCell, nCells, &
- angleEdge, nEdgesOnCell, edgesOnCell )
-
- implicit none
-
- logical :: is_periodic, on_a_sphere
- integer :: nCells
- integer, dimension(:) :: nEdgesOnCell
- real (kind=RKIND) :: x_period, y_period
- real (kind=RKIND), dimension(:,:) :: deformation_coef_c2, deformation_coef_s2
- real (kind=RKIND), dimension(:,:) :: deformation_coef_cs
- real (kind=RKIND), dimension(:,:) :: deformation_coef_c, deformation_coef_s
- integer, dimension(:,:) :: edgesOnCell
- real (kind=RKIND), dimension(:) :: angleEdge, xEdge, yEdge, zEdge
- real (kind=RKIND), dimension(:) :: xCell, yCell, zCell
-
- ! local variables
-
- integer :: iCell, iEdge, ie
- real (kind=RKIND) :: cos_edge, sin_edge, ux, uy, vx, vy, wx, wy
- real (kind=RKIND) :: xc, yc, xe, ye
- real (kind=RKIND) :: angle_e, ue, ve, we, e_int
- real (kind=RKIND) :: dudx, dudy, dvdx, dvdy, dwdx, dwdy
- real (kind=RKIND) :: dudx_c, dudy_c, dvdx_c, dvdy_c, dwdx_c, dwdy_c
-
- real (kind=RKIND) :: dudx_err_max, dudy_err_max, dvdx_err_max, dvdy_err_max, dwdx_err_max, dwdy_err_max
- real (kind=RKIND) :: dudx_err_tot, dudy_err_tot, dvdx_err_tot, dvdy_err_tot, dwdx_err_tot, dwdy_err_tot
- real (kind=RKIND) :: dudx_max, dudy_max, dvdx_max, dvdy_max, dwdx_max, dwdy_max
-
- real (kind=RKIND) :: ang
- real (kind=RKIND), parameter :: x_vel= 1.0, y_vel=1.0, w_vel=1.0
- real (kind=RKIND) :: u_edge, v_edge, w_edge, x, y, angle, xl, yl
- real (kind=RKIND) :: dudx_cell, dudy_cell, dvdx_cell, dvdy_cell, dwdx_cell, dwdy_cell
-
- ! Test tunction definitions
- !
- ! here are the velocity field functions and their derivatives.
- ! First a simple test: U = x_vel*(-x+y), V = y_vel * (-x+y), W = w_vel*(-x+y)
-
- u_edge(x,y,ang,xl,yl) = (x_vel*(x+y)) * cos(ang) + (y_vel * (x+y) * sin(ang))
- v_edge(x,y,ang,xl,yl) = -(x_vel*(x+y)) * sin(ang) + (y_vel * (x+y) * cos(ang))
- w_edge(x,y,xl,yl) = w_vel * (x+y)
-
- dudx_cell(x,y,xl,yl) = x_vel
- dudy_cell(x,y,xl,yl) = x_vel
- dvdx_cell(x,y,xl,yl) = y_vel
- dvdy_cell(x,y,xl,yl) = y_vel
- dwdx_cell(x,y,xl,yl) = w_vel
- dwdy_cell(x,y,xl,yl) = w_vel
-
- ! -----------------
-
- if ( (.not. on_a_sphere) .and. (is_periodic) ) then ! test is for doubly-periodic Cartesian plane only
-
- dudx_err_max = 0.
- dudy_err_max = 0.
- dvdx_err_max = 0.
- dvdy_err_max = 0.
- dwdx_err_max = 0.
- dwdy_err_max = 0.
-
- dudx_err_tot = 0.
- dudy_err_tot = 0.
- dvdx_err_tot = 0.
- dvdy_err_tot = 0.
- dwdx_err_tot = 0.
- dwdy_err_tot = 0.
-
- dudx_max = 0.
- dudy_max = 0.
- dvdx_max = 0.
- dvdy_max = 0.
- dwdx_max = 0.
- dwdy_max = 0.
-
- do iCell = 1, nCells
-
- dudx = 0.
- dudy = 0.
- dvdx = 0.
- dvdy = 0.
- dwdx = 0.
- dwdy = 0.
-
- xc = xCell(iCell)
- yc = yCell(iCell)
-
- dudx_c = dudx_cell(xc,yc,x_period,y_period)
- dudy_c = dudy_cell(xc,yc,x_period,y_period)
- dvdx_c = dvdx_cell(xc,yc,x_period,y_period)
- dvdy_c = dvdy_cell(xc,yc,x_period,y_period)
- dwdx_c = dwdx_cell(xc,yc,x_period,y_period)
- dwdy_c = dwdy_cell(xc,yc,x_period,y_period)
-
- do iEdge = 1, nEdgesOnCell(iCell)
-
- ie = edgesOnCell(iEdge,iCell)
- angle_e = angleEdge(ie)
- xe = xEdge(ie)
- ye = yEdge(ie)
-
- xe = mpas_fix_periodicity(xe,xc,x_period)
- ye = mpas_fix_periodicity(ye,yc,y_period)
-
- ue = u_edge(xe,ye,angle_e,x_period,y_period)
- ve = v_edge(xe,ye,angle_e,x_period,y_period)
- we = w_edge(xe,ye,x_period,y_period)
-
- dudx = dudx + deformation_coef_c2(iEdge,iCell)*ue &
- - deformation_coef_cs(iEdge,iCell)*ve
- dudy = dudy + deformation_coef_cs(iEdge,iCell)*ue &
- - deformation_coef_s2(iEdge,iCell)*ve
- dvdx = dvdx + deformation_coef_cs(iEdge,iCell)*ue &
- + deformation_coef_c2(iEdge,iCell)*ve
- dvdy = dvdy + deformation_coef_s2(iEdge,iCell)*ue &
- + deformation_coef_cs(iEdge,iCell)*ve
-
- dwdx = dwdx + deformation_coef_c(iEdge,iCell)*we
- dwdy = dwdy + deformation_coef_s(iEdge,iCell)*we
-
- end do
-
- ! call mpas_log_write(' u_x, u_y, $r, $r ', realArgs=(/dudx, dudy/))
- ! call mpas_log_write(' v_x, v_y, $r, $r ', realArgs=(/dvdx, dvdy/))
- ! call mpas_log_write(' w_x, w_y, $r, $r ', realArgs=(/dwdx, dwdy/))
-
- ! check result for cell
-
- e_int = abs(dudx_c - dudx)
- dudx_err_tot = dudx_err_tot + e_int
- dudx_err_max = max(dudx_err_max, e_int)
-
- e_int = abs(dudy_c - dudy)
- dudy_err_tot = dudy_err_tot + e_int
- dudy_err_max = max(dudy_err_max, e_int)
-
- e_int = abs(dvdx_c - dvdx)
- dvdx_err_tot = dvdx_err_tot + e_int
- dvdx_err_max = max(dvdx_err_max, e_int)
-
- e_int = abs(dvdy_c - dvdy)
- dvdy_err_tot = dvdy_err_tot + e_int
- dvdy_err_max = max(dvdy_err_max, e_int)
-
- e_int = abs(dwdx_c - dwdx)
- dwdx_err_tot = dwdx_err_tot + e_int
- dwdx_err_max = max(dwdx_err_max, e_int)
-
- e_int = abs(dwdy_c - dwdy)
- dwdy_err_tot = dwdy_err_tot + e_int
- dwdy_err_max = max(dwdy_err_max, e_int)
-
- dudx_max = max(dudx_max, abs(dudx_c))
- dudy_max = max(dudy_max, abs(dudy_c))
- dvdx_max = max(dvdx_max, abs(dvdx_c))
- dvdy_max = max(dvdy_max, abs(dvdy_c))
- dwdx_max = max(dwdx_max, abs(dwdx_c))
- dwdy_max = max(dwdy_max, abs(dwdy_c))
-
- end do
-
- ! scale errors
-
- dudx_err_max = dudx_err_max/dudx_max
- dudy_err_max = dudy_err_max/dudy_max
- dvdx_err_max = dvdx_err_max/dvdx_max
- dvdy_err_max = dvdy_err_max/dvdy_max
- dwdx_err_max = dwdx_err_max/dwdx_max
- dwdy_err_max = dwdy_err_max/dwdy_max
-
- dudx_err_tot = dudx_err_tot/dudx_max/real(nCells)
- dudy_err_tot = dudy_err_tot/dudy_max/real(nCells)
- dvdx_err_tot = dvdx_err_tot/dvdx_max/real(nCells)
- dvdy_err_tot = dvdy_err_tot/dvdy_max/real(nCells)
- dwdx_err_tot = dwdx_err_tot/dwdx_max/real(nCells)
- dwdy_err_tot = dwdy_err_tot/dwdy_max/real(nCells)
-
- ! output
-
- call mpas_log_write(' ')
- call mpas_log_write(' deformation coefficients check ')
- call mpas_log_write(' dudx check, max abs(dudx), max and avg error $r, $r, $r', &
- realArgs=(/dudx_max, dudx_err_max, dudx_err_tot/))
- call mpas_log_write(' dudy check, max abs(dudy), max and avg error $r, $r, $r', &
- realArgs=(/dudy_max, dudy_err_max, dudy_err_tot/))
- call mpas_log_write(' dvdx check, max abs(dvdx), max and avg error $r, $r, $r', &
- realArgs=(/dvdx_max, dvdx_err_max, dvdx_err_tot/))
- call mpas_log_write(' dvdy check, max abs(dvdy), max and avg error $r, $r, $r', &
- realArgs=(/dvdy_max, dvdy_err_max, dvdy_err_tot/))
- call mpas_log_write(' dwdx check, max abs(dwdx), max and avg error $r, $r, $r', &
- realArgs=(/dwdx_max, dwdx_err_max, dwdx_err_tot/))
- call mpas_log_write(' dwdy check, max abs(dwdy), max and avg error $r, $r, $r', &
- realArgs=(/dwdy_max, dwdy_err_max, dwdy_err_tot/))
- call mpas_log_write(' ')
-
- end if
-
- end subroutine atm_init_test_coefs
-
-end module atm_advection
+ end module atm_advection