<p><b>duda</b> 2010-04-23 14:43:27 -0600 (Fri, 23 Apr 2010)</p><p>BRANCH COMMIT<br>
<br>
In MPAS_INIT1(...), add computation of unit vectors in the east and north directions<br>
at each cell, and also add computation of the unit normal vector for each edge. These<br>
unit vectors will be needed when transforming wind vectors in the MPAS_TO_CAM(...) <br>
and CAM_TO_MPAS(...) routines.<br>
<br>
<br>
M core_hyd_atmos/Registry<br>
M driver_cam_interface/module_mpas_cam_interface.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/mpas_cam_coupling/src/core_hyd_atmos/Registry
===================================================================
--- branches/mpas_cam_coupling/src/core_hyd_atmos/Registry 2010-04-23 20:28:13 UTC (rev 207)
+++ branches/mpas_cam_coupling/src/core_hyd_atmos/Registry 2010-04-23 20:43:27 UTC (rev 208)
@@ -159,3 +159,9 @@
# 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 - -
Modified: branches/mpas_cam_coupling/src/driver_cam_interface/module_mpas_cam_interface.F
===================================================================
--- branches/mpas_cam_coupling/src/driver_cam_interface/module_mpas_cam_interface.F 2010-04-23 20:28:13 UTC (rev 207)
+++ branches/mpas_cam_coupling/src/driver_cam_interface/module_mpas_cam_interface.F 2010-04-23 20:43:27 UTC (rev 208)
@@ -33,7 +33,15 @@
implicit none
integer, intent(in) :: mpi_comm
-
+
+ integer :: iCell, iEdge, j
+ real (kind=8) :: m, angle
+ real (kind=8), dimension(3) :: X1, X2, n_rot_vector
+ integer, dimension(:), pointer :: nEdgesOnCell
+ real (kind=8), dimension(:), pointer :: xCell, yCell, zCell
+ integer, dimension(:,:), pointer :: edgesOnCell, cellsOnEdge
+ real (kind=8), dimension(:,:), pointer :: east, north, edge_normal
+
write(0,*) 'Called MPAS_INIT1'
allocate(dminfo)
@@ -44,9 +52,76 @@
call allocate_domain(domain, dminfo)
! NOTE: Really, we only want to read in the grid information, since any state
- ! variables will be overwritten by cam_inidat_to_mpas() or cam_restart_to_mpas();
- ! however, with a suitable registry, this shouldn't be a problem
+ ! variables will be overwritten by cam_inidat_to_mpas() or cam_restart_to_mpas();
+ ! however, with a suitable registry, this shouldn't be a problem
call input_state_for_domain(domain)
+
+ xCell => domain % blocklist % mesh % xCell % array
+ yCell => domain % blocklist % mesh % yCell % array
+ zCell => domain % blocklist % mesh % zCell % array
+ nEdgesOnCell => domain % blocklist % mesh % nEdgesOnCell % array
+ edgesOnCell => domain % blocklist % mesh % edgesOnCell % array
+ cellsOnEdge => domain % blocklist % mesh % cellsOnEdge % array
+ east => domain % blocklist % mesh % east % array
+ north => domain % blocklist % mesh % north % array
+ edge_normal => domain % blocklist % mesh % edge_normal % array
+
+ angle = 0.0003
+
+
+ !
+ ! Compute unit vectors in east and north directions for each cell
+ !
+ do iCell = 1,domain % blocklist % mesh % 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 mpas_init1
@@ -601,4 +676,77 @@
end subroutine mpas_final
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine rotate_about_vector(x, y, z, theta, a, b, c, u, v, w, xp, yp, zp)
+
+ 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
+ !
+ ! Computes the cross product of (ax, ay, az) and (bx, by, bz).
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine r3_cross(ax, ay, az, bx, by, bz, cx, cy, cz)
+
+ 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
+ !
+ ! Normalizes the vector (ax, ay, az)
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine r3_normalize(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 mpas_cam_interface
</font>
</pre>