<p><b>duda</b> 2010-05-25 15:58:11 -0600 (Tue, 25 May 2010)</p><p>BRANCH COMMIT<br>
<br>
Add new module_earth_winds.F to compute zonal and meridional components<br>
of winds (Ux and Uy, respectively) at cell centers given reconstructed wind<br>
vector components in Cartesian space.<br>
<br>
<br>
A src/core_hyd_atmos/module_earth_winds.F<br>
M src/core_hyd_atmos/mpas_interface.F<br>
M src/core_hyd_atmos/Registry<br>
M src/core_hyd_atmos/module_time_integration.F<br>
M src/core_hyd_atmos/Makefile<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_hyd_atmos/Makefile
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/Makefile 2010-05-25 21:07:20 UTC (rev 312)
+++ branches/atmos_physics/src/core_hyd_atmos/Makefile 2010-05-25 21:58:11 UTC (rev 313)
@@ -3,6 +3,7 @@
OBJS = module_test_cases.o \
module_time_integration.o \
module_advection.o \
+ module_earth_winds.o \
mpas_interface.o
all: core_hyd
@@ -12,12 +13,14 @@
module_test_cases.o:
-module_time_integration.o:
+module_time_integration.o: module_earth_winds.o
module_advection.o:
-mpas_interface.o: module_advection.o module_test_cases.o module_time_integration.o
+module_earth_winds.o:
+mpas_interface.o: module_advection.o module_earth_winds.o module_test_cases.o module_time_integration.o
+
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/atmos_physics/src/core_hyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/Registry 2010-05-25 21:07:20 UTC (rev 312)
+++ branches/atmos_physics/src/core_hyd_atmos/Registry 2010-05-25 21:58:11 UTC (rev 313)
@@ -123,6 +123,8 @@
var real pressure ( nVertLevelsP1 nCells Time ) ro pressure - -
var real geopotential ( nVertLevelsP1 nCells Time ) ro geopotential - -
var real alpha ( nVertLevels nCells Time ) iro alpha - -
+var real Ux ( nVertLevels nCells Time ) o Ux - -
+var real Uy ( nVertLevels nCells Time ) o Uy - -
# Diagnostic fields: only written to output
var real v ( nVertLevels nEdges Time ) o v - -
@@ -173,6 +175,13 @@
# Arrays required for reconstruction of velocity field
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
+# Unit vectors in east and north directions for each cell
+var real east ( R3 nCells ) - east - -
+var real north ( R3 nCells ) - north - -
+
+# Unit vectors in the positive normal direction for each edge
+var real edge_normal ( R3 nEdges ) - edge_normal - -
+
#==================================================================================================
#... PARAMETERIZATION OF CONVECTION:
Added: branches/atmos_physics/src/core_hyd_atmos/module_earth_winds.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/module_earth_winds.F (rev 0)
+++ branches/atmos_physics/src/core_hyd_atmos/module_earth_winds.F 2010-05-25 21:58:11 UTC (rev 313)
@@ -0,0 +1,221 @@
+module earth_winds
+
+ use grid_types
+ use configure
+ use constants
+
+
+ contains
+
+
+ subroutine initialize_earth_winds( grid )
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute unit vectors in east and north directions, which are used by
+ ! the subroutine compute_earth_winds
+ !
+ ! Input: grid - grid meta-data fields
+ !
+ ! Output: grid - upon exit, the fields 'east' and 'north' are filled in with
+ ! components of the unit east and north vectors for each cell.
+ ! Also, edge_normal contains components of the unit vectors in
+ ! the positive normal direction for each edge.
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_meta), intent(inout) :: grid
+
+ integer :: iCell, iEdge, j
+ real (kind=RKIND) :: angle
+ real (kind=RKIND), dimension(3) :: X1, X2, n_rot_vector
+ integer, dimension(:), pointer :: nEdgesOnCell
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+ integer, dimension(:,:), pointer :: edgesOnCell, cellsOnEdge
+ real (kind=RKIND), dimension(:,:), pointer :: east, north, edge_normal
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ east => grid % east % array
+ north => grid % north % array
+ edge_normal => grid % edge_normal % array
+
+ angle = 0.0003
+
+
+ !
+ ! Compute unit vectors in east and north directions for each cell
+ !
+ do iCell = 1,grid % nCellsSolve
+ call rotate_about_vector(xCell(iCell), yCell(iCell), zCell(iCell), &
+ angle, &
+ 0., 0., 0., &
+ 0., 0., 1., &
+ east(1,iCell), east(2,iCell), east(3,iCell))
+ call rotate_about_vector(xCell(iCell), yCell(iCell), zCell(iCell), &
+ -angle, &
+ 0., 0., 0., &
+ 0., 0., 1., &
+ X1(1), X1(2), X1(3))
+ east(1,iCell) = east(1,iCell) - X1(1)
+ east(2,iCell) = east(2,iCell) - X1(2)
+ east(3,iCell) = east(3,iCell) - X1(3)
+ call r3_normalize(east(1,iCell), east(2,iCell), east(3,iCell))
+
+ call r3_cross(xCell(iCell), yCell(iCell), zCell(iCell), &
+ 0., 0., 1., &
+ n_rot_vector(1), n_rot_vector(2), n_rot_vector(3))
+ call rotate_about_vector(xCell(iCell), yCell(iCell), zCell(iCell), &
+ angle, &
+ 0., 0., 0., &
+ n_rot_vector(1), n_rot_vector(2), n_rot_vector(3), &
+ north(1,iCell), north(2,iCell), north(3,iCell))
+ call rotate_about_vector(xCell(iCell), yCell(iCell), zCell(iCell), &
+ -angle, &
+ 0., 0., 0., &
+ n_rot_vector(1), n_rot_vector(2), n_rot_vector(3), &
+ X1(1), X1(2), X1(3))
+ north(1,iCell) = north(1,iCell) - X1(1)
+ north(2,iCell) = north(2,iCell) - X1(2)
+ north(3,iCell) = north(3,iCell) - X1(3)
+ call r3_normalize(north(1,iCell), north(2,iCell), north(3,iCell))
+
+ X1(1) = xCell(iCell)
+ X1(2) = yCell(iCell)
+ X1(3) = zCell(iCell)
+
+ do j=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(j,iCell)
+ if (iCell == cellsOnEdge(1,iEdge)) then
+ X2(1) = xCell(cellsOnEdge(2,iEdge))
+ X2(2) = yCell(cellsOnEdge(2,iEdge))
+ X2(3) = zCell(cellsOnEdge(2,iEdge))
+ edge_normal(:,iEdge) = X2(:) - X1(:)
+ call r3_normalize(edge_normal(1,iEdge), edge_normal(2,iEdge), edge_normal(3,iEdge))
+ endif
+ enddo
+
+ end do
+
+ end subroutine initialize_earth_winds
+
+
+ subroutine compute_earth_winds( state, grid )
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Compute unit vectors in east and north directions, which are used by
+ ! the subroutine compute_earth_winds
+ !
+ ! Input: state - reconstructed velocity vector fields uReconstruct{x,y,z}
+ ! grid - 'east' and 'north' fields
+ !
+ ! Output: state - upon exit, the 'Ux' and 'Uy' fields will contain the zonal
+ ! and meridional components, respectively, of the interpolated
+ ! wind field given in uReconstruct{x,y,z}
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_state), intent(inout) :: state
+ type (grid_meta), intent(in) :: grid
+
+ integer :: iCell, k
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstX, uReconstY, uReconstZ, east, north, Ux, Uy
+
+ east => grid % east % array
+ north => grid % north % array
+ uReconstX => state % uReconstructX % array
+ uReconstY => state % uReconstructY % array
+ uReconstZ => state % uReconstructZ % array
+ Ux => state % Ux % array
+ Uy => state % Uy % array
+
+ do iCell=1,grid % nCellsSolve
+ do k=1,grid % nVertLevels
+ Ux(k,iCell) = uReconstX(k,iCell) * east(1,iCell) &
+ + uReconstY(k,iCell) * east(2,iCell) &
+ + uReconstZ(k,iCell) * east(3,iCell)
+ Uy(k,iCell) = uReconstX(k,iCell) * north(1,iCell) &
+ + uReconstY(k,iCell) * north(2,iCell) &
+ + uReconstZ(k,iCell) * north(3,iCell)
+ end do
+ end do
+
+ end subroutine compute_earth_winds
+
+
+ subroutine rotate_about_vector(x, y, z, theta, a, b, c, u, v, w, xp, yp, zp)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE ROTATE_ABOUT_VECTOR
+ !
+ ! Rotates the point (x,y,z) through an angle theta about the vector
+ ! originating at (a, b, c) and having direction (u, v, w).
+ !
+ ! Reference: http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/ArbitraryAxisRotation.html
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: x, y, z, theta, a, b, c, u, v, w
+ real (kind=RKIND), intent(out) :: xp, yp, zp
+
+ real (kind=RKIND) :: vw2, uw2, uv2
+ real (kind=RKIND) :: m
+
+ vw2 = v**2.0 + w**2.0
+ uw2 = u**2.0 + w**2.0
+ uv2 = u**2.0 + v**2.0
+ m = sqrt(u**2.0 + v**2.0 + w**2.0)
+
+ xp = (a*vw2 + u*(-b*v-c*w+u*x+v*y+w*z) + ((x-a)*vw2+u*(b*v+c*w-v*y-w*z))*cos(theta) + m*(-c*v+b*w-w*y+v*z)*sin(theta))/m**2.0
+ yp = (b*uw2 + v*(-a*u-c*w+u*x+v*y+w*z) + ((y-b)*uw2+v*(a*u+c*w-u*x-w*z))*cos(theta) + m*( c*u-a*w+w*x-u*z)*sin(theta))/m**2.0
+ zp = (c*uv2 + w*(-a*u-b*v+u*x+v*y+w*z) + ((z-c)*uv2+w*(a*u+b*v-u*x-v*y))*cos(theta) + m*(-b*u+a*v-v*x+u*y)*sin(theta))/m**2.0
+
+ end subroutine rotate_about_vector
+
+
+ subroutine r3_cross(ax, ay, az, bx, by, bz, cx, cy, cz)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE R3_CROSS
+ !
+ ! Computes the cross product of (ax, ay, az) and (bx, by, bz).
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: ax, ay, az
+ real (kind=RKIND), intent(in) :: bx, by, bz
+ real (kind=RKIND), intent(out) :: cx, cy, cz
+
+ cx = ay * bz - az * by
+ cy = az * bx - ax * bz
+ cz = ax * by - ay * bx
+
+ end subroutine r3_cross
+
+
+ subroutine r3_normalize(ax, ay, az)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! SUBROUTINE R3_NORMALIZE
+ !
+ ! Normalizes the vector (ax, ay, az)
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ real (kind=RKIND), intent(inout) :: ax, ay, az
+
+ real (kind=RKIND) :: mi
+
+ mi = 1.0 / sqrt(ax**2 + ay**2 + az**2)
+ ax = ax * mi
+ ay = ay * mi
+ az = az * mi
+
+ end subroutine r3_normalize
+
+end module earth_winds
Modified: branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F 2010-05-25 21:07:20 UTC (rev 312)
+++ branches/atmos_physics/src/core_hyd_atmos/module_time_integration.F 2010-05-25 21:58:11 UTC (rev 313)
@@ -5,6 +5,7 @@
use constants
use dmpar
use vector_reconstruction
+ use earth_winds
#ifdef DO_PHYSICS
use module_microphysics_driver
@@ -355,6 +356,7 @@
block => domain % blocklist
do while (associated(block))
call reconstruct(block % time_levs(2) % state, block % mesh)
+ call compute_earth_winds( block % time_levs(2) % state, block % mesh )
block => block % next
end do
Modified: branches/atmos_physics/src/core_hyd_atmos/mpas_interface.F
===================================================================
--- branches/atmos_physics/src/core_hyd_atmos/mpas_interface.F 2010-05-25 21:07:20 UTC (rev 312)
+++ branches/atmos_physics/src/core_hyd_atmos/mpas_interface.F 2010-05-25 21:58:11 UTC (rev 313)
@@ -17,6 +17,7 @@
use grid_types
use advection
+ use earth_winds
use time_integration
#ifdef DO_PHYSICS
use module_physics_init
@@ -35,6 +36,8 @@
call initialize_advection_rk(mesh)
call init_reconstruct(mesh)
call reconstruct(block % time_levs(1) % state, mesh)
+ call initialize_earth_winds(mesh)
+ call compute_earth_winds(block % time_levs(1) % state, mesh)
#ifdef DO_PHYSICS
call physics_wrf_allocate(mesh)
</font>
</pre>