!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! rot_mod - rotates a spherical coordinate system, from M. Eby via MOM ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! module rot_mod use mpp_mod, only: stdout use axis_utils_mod, only: nearest_index implicit none private public rot_intrp_sclr, rot_intrp_vctr, rotate, psir, thetar, phir ! rotation angles set in ice_grid real :: psir = 0.0, thetar = 0.0, phir = 0.0 integer :: warns = 0 ! > 0 allows #warns warnings from interpolator contains subroutine rot_intrp_vctr (g, xg, yg, ig, jg, r, xr, yr, ir, jr, & psir, thetar, phir) ! !======================================================================= ! interpolate vector data from an geographic data grid to a ! rotated model grid ! ! input ! psir, thetar, phir = Euler angles defining rotation ! g = vector on geographic data grid ! xg = longitude of data points on geographic data grid ! yg = latitude of data points on geographic data grid ! ig = number of longitudes in geographic data grid ! jg = number of latitudes in geographic data grid ! xr = longitude of points on rotated model grid ! yr = latitude of points on rotated model grid ! ir = number of longitudes in rotated model grid ! jr = number of latitudes in rotated model grid ! ! output ! r = vector on rotated model grid ! ! internal ! (rln,rlt) = (longitude,latitude) in rotated coordinates ! (gln,glt) = (longitude,latitude) in geographic coordinates ! xg(iw) = point on the geographic grid to the west of (gln,glt) ! xg(ie) = point on the geographic grid to the east of (gln,glt) ! yg(js) = point on the geographic grid to the south of (gln,glt) ! yg(jn) = point on the geographic grid to the north of (gln,glt) ! !======================================================================= ! integer :: ig, jg, ir, jr real, dimension(ig,jg,2) :: g real, dimension(ig ) :: xg real, dimension( jg ) :: yg real, dimension(ir,jr,2) :: r real, dimension(ir ) :: xr real, dimension( jr ) :: yr real :: psir, thetar, phir integer :: i, j real :: rad, a, angle, vmag, glt, gln, rlt, rln ! rad = acos(-1.)/180. ! ! interpolate vector components as scalers on rotated model grid ! call rot_intrp_sclr (g(1,1,1), xg, yg, ig, jg, r(1,1,1), xr, yr, ir, jr, & psir, thetar, phir) call rot_intrp_sclr (g(1,1,2), xg, yg, ig, jg, r(1,1,2), xr, yr, ir, jr, & psir, thetar, phir) ! ! correct vector direction ! do j=1,jr do i=1,ir vmag = sqrt(r(i,j,1)**2 + r(i,j,2)**2) if (vmag .gt. 0.) then a = r(i,j,1)/vmag a = min(a, 1.) a = max(a, -1.) a = acos(a) if (r(i,j,2) .lt. 0.) a = -a call rotate (yr(j), xr(i), -psir, -thetar, -phir, glt, gln) call rotvec(glt, gln, phir, thetar, psir, angle) a = a + angle*rad r(i,j,1) = vmag*cos(a) r(i,j,2) = vmag*sin(a) else r(i,j,1) = 0. r(i,j,2) = 0. endif enddo enddo ! return end subroutine rot_intrp_vctr subroutine rot_intrp_sclr (g, xg, yg, ig, jg, r, xr, yr, & ir, jr, psir, thetar, phir ) ! !======================================================================= ! interpolate scaler data from an geographic data grid to a ! rotated model grid ! ! input ! psir, thetar, phir = Euler angles defining rotation ! g = scaler on geographic data grid ! xg = longitude of data points on geographic data grid ! yg = latitude of data points on geographic data grid ! ig = number of longitudes in on geographic data grid ! jg = number of latitudes in on geographic data grid ! xr = longitude of points on rotated model grid ! yr = latitude of points on rotated model grid ! ir = number of longitudes in rotated model grid ! jr = number of latitudes in rotated model grid ! ! output ! r = scaler on rotated model grid ! ! internal ! (rln,rlt) = (longitude,latitude) in rotated coordinates ! (gln,glt) = (longitude,latitude) in geographic coordinates ! xg(iw) = point on the geographic grid to the west of (gln,glt) ! xg(ie) = point on the geographic grid to the east of (gln,glt) ! yg(js) = point on the geographic grid to the south of (gln,glt) ! yg(jn) = point on the geographic grid to the north of (gln,glt) ! !======================================================================= ! integer :: ig, jg, ir, jr, outunit real :: psir, thetar, phir real, dimension(ig,jg) :: g real, dimension(ig ) :: xg real, dimension( jg) :: yg real, dimension(ir,jr) :: r real, dimension(ir ) :: xr real, dimension( jr) :: yr real :: rln, rlt, gln, glt integer :: i, j, iw, ie, js, jn, istrt, jstrt, iend, jend, ln_err, lt_err real :: epsln, glt_min, glt_max, gln_min, gln_max, del, wtw, wte, wts, wtn ! epsln = 1.e-10 glt_min = 90. glt_max = -90. gln_min = 360. gln_max = -360. ln_err = 0 lt_err = 0 ! ! find longitude points of data within interval [0., 360.] istrt = 1 do i=2,ig if (xg(i-1) .lt. 0. .and. xg(i) .ge. 0.) istrt = i enddo iend = ig do i=2,ig if (xg(i-1) .lt. 360. .and. xg(i) .ge. 360.) iend = i enddo ! ! find latitude points of data within interval [-90., 90.] jstrt = 1 do j=2,jg if (yg(j-1) .lt. -90. .and. yg(j) .ge. -90.) jstrt = j enddo jend = jg do j=2,jg if (yg(j-1) .lt. 90. .and. yg(j) .ge. 90.) jend = j enddo ! ! interpolate data to model grid ! do j=1,jr do i=1,ir call rotate (yr(j), xr(i), -psir, -thetar, -phir, glt, gln) if (gln .lt. 0.) gln = gln + 360. if (gln .ge. 360.) gln = gln - 360. glt_min = min(glt,glt_min) glt_max = max(glt,glt_max) gln_min = min(gln,gln_min) gln_max = max(gln,gln_max) ! iw = nearest_index (gln, xg(istrt:iend) ) + istrt - 1 if (xg(iw) .gt. gln) iw = iw - 1 ie = iw + 1 if (iw .ge. istrt .and. ie .le. iend) then del = xg(ie) - xg(iw) wtw = (xg(ie) - gln)/del else ! east or west of the last data value. this could be because a ! cyclic condition is needed or the dataset is too small. in either ! case apply a cyclic condition ln_err = 1 iw = iend ie = istrt del = xg(ie) + 360. + epsln - xg(iw) if (xg(ie) .ge. gln) then wtw = (xg(ie) - gln)/del else wtw = (xg(ie) + 360. + epsln - gln)/del endif endif wte = 1. - wtw ! js = nearest_index (glt, yg(jstrt:jend) ) + jstrt - 1 if (yg(js) .gt. glt) js = max(js - 1,jstrt) jn = min(js + 1,jend) if (yg(jn) .ne. yg(js) .and. yg(js) .le. glt) then wts = (yg(jn) - glt)/(yg(jn) - yg(js)) else ! north or south of the last data value. this could be because a ! pole is not included in the data set or the dataset is too small. ! in either case extrapolate north or south lt_err = 1 wts = 1. endif wtn = 1. - wts ! r(i,j) = g(ie,jn)*wte*wtn + g(ie,js)*wte*wts & + g(iw,jn)*wtw*wtn + g(iw,js)*wtw*wts ! enddo enddo ! outunit = stdout() if (ln_err .eq. 1 .and. warns > 0) then write (outunit,'(/,(1x,a))') & '==> Warning: the geographic data set does not extend far ', & ' enough east or west - a cyclic boundary ', & ' condition was applied. check if appropriate ' write (outunit,'(/,(1x,a,2f8.2))') & ' data required between longitudes:', gln_min, gln_max, & ' data set is between longitudes:', xg(istrt), xg(iend) warns = warns - 1 endif ! if (lt_err .eq. 1 .and. warns > 0) then write (outunit,'(/,(1x,a))') & '==> Warning: the geographic data set does not extend far ',& ' enough north or south - extrapolation from ',& ' the nearest data was applied. this may create ',& ' artificial gradients near a geographic pole ' write (outunit,'(/,(1x,a,2f8.2))') & ' data required between latitudes:', glt_min, glt_max, & ' data set is between latitudes:', yg(jstrt), yg(jend) warns = warns - 1 endif ! return end subroutine rot_intrp_sclr subroutine rotate (glt, gln, phir, thetar, psir, rlt, rln) real :: glt, gln, phir, thetar, psir, rlt, rln ! !======================================================================= ! subroutine rotate takes a geographic latitude and longitude and ! finds the the equivalent latitude and longitude on a rotated grid. ! when going from a geographic grid to a rotated grid, all of the ! defined rotation angles given to rotate by the calling program ! are positive, but when going from a rotated grid back to the ! geographic grid, the calling program must reverse the angle order ! (phir and psir are switched) and all of the angles made negative. ! ! the first rotation angle phir is defined as a rotation about the ! original z axis. the second rotation angle thetar is defined as a ! rotation about the new x axis. the final rotation angle psir is ! defined as a rotation about the new z axis. these rotation angles ! are just the Euler angles as defined in "classical mechanics" ! Goldstein (1951). ! ! author: M. Eby e-mail eby@uvic.ca !======================================================================= ! ! g... = geographic value ! r... = rotated value ! ...lt = latitude (or equivalent spherical coordinate) ! ...ln = longitude (or equivalent spherical coordinate) ! ...x = x coordinate ! ...y = y coordinate ! ...z = z coordinate ! psir, thetar, phir = Euler angles defining rotation ! real rad, gx, gy, gz, phis, thetas, rx, ry, rz ! define rad for conversion to radians. rad = acos(-1.)/180. ! ! convert latitude and longitude to spherical coordinates thetas = gln if (thetas .gt. 180.) thetas = thetas - 360. if (thetas .lt. -180.) thetas = thetas + 360. phis = (90. - glt)*rad thetas = thetas*rad ! ! translate point into Cartesian coordinates for rotation. gx = sin(phis)*cos(thetas) gy = sin(phis)*sin(thetas) gz = cos(phis) ! ! rotate the point (gx, gy, gz) about the z axis by phir then the x ! axis by thetar and finally about the z axis by psir. ! rx = gx*(cos(psir)*cos(phir) - cos(thetar)*sin(phir)*sin(psir)) + & gy*(cos(psir)*sin(phir) + cos(thetar)*cos(phir)*sin(psir)) + & gz*sin(psir)*sin(thetar) ! ry = gx*(-sin(psir)*cos(phir) - cos(thetar)*sin(phir)*cos(psir)) + & gy*(-sin(psir)*sin(phir) + cos(thetar)*cos(phir)*cos(psir)) + & gz*(cos(psir)*sin(thetar)) ! rz = gx*(sin(thetar)*sin(phir)) + gy*(-sin(thetar)*cos(phir)) + & gz*(cos(thetar)) ! ! convert rotated point back to spherical coordinates ! ! check for rounding error (arccos(x): abs(x) must be .le. 1) rz = min(rz, 1.) rz = max(rz, -1.) rlt = acos(rz) ! if point is at a pole set rotated longitude equal to initial. if (rlt .le. 0. .or. rlt .ge. 180.*rad) then rln = thetas else ! if rln lies between -135 and -45 or between 45 and 135 degrees ! it is more accurate to use an arccos calculation. if (abs(rx/sin(rlt)) .lt. cos(45.*rad)) then rln = rx/sin(rlt) ! check for rounding error (arccos(x): abs(x) must be .le. 1) rln = min(rln, 1.) rln = max(rln, -1.) rln = acos(rln) ! arccos will give rln between 0 and 180 degrees. if the point ! is negative in y, rln must be equal to negative rln. if (ry .lt. 0.) rln = -rln else ! if rln lies between -45 and 45 or between 135 and -135 degrees ! it is more accurate to use an arcsin calculation. rln = ry/sin(rlt) ! check for rounding error (arcsin(x): abs(x) must be .le. 1) rln = min(rln, 1.) rln = max(rln, -1.) rln = asin(rln) ! arcsin will give rln between -90 and 90 degrees. if the point ! is negative in x, rln must be equal to 180 degrees minus rln. if (rx .lt. 0.) rln = 180.*rad - rln endif endif ! ! convert back to degrees of latitude and longitude. rlt = 90. - rlt/rad rln = rln/rad if (rln .gt. 180.) rln = rln - 360. if (rln .le. -180.) rln = rln + 360. ! return end subroutine rotate subroutine rotvec (glt, gln, phir, thetar, psir, angle) real :: glt, gln, phir, thetar, psir, angle real rad, delta, rlt, rln, rlth, rlnh, glth, glnh, dst, t ! !======================================================================= ! subroutine rotvec takes a geographic latitude and longitude and ! finds the the vector rotation angle angle (in degrees) for a ! vector on the rotated grid. when going from the geographic to a ! rotated grid, all of the defined rotation angles given to rotvec ! by the calling program are positive, but when going from a ! rotated grid back to the geographic grid, the calling program ! must reverse the angle order (phir and psir are switched) and all ! of the angles made negative. if a pole is detected then an angle ! of zero is returned. ! ! rotvec rotates the point defining the head of a very short north ! or south pointing geographic vector. the angle between this ! vector and a similar direction vector defined in the new grid is ! calculated using the law of cosines. the accuracy of this ! calculation depends on the size of the direction vector (delta) ! and the precision of the computation. the smaller the vector the ! more accurate the calculation but then the more precision ! required. double precision is strongly recommended. ! ! author: M. Eby e-mail eby@uvic.ca !======================================================================= ! ! include "stdunits.h" ! if gx and gy are vector components at glt and gln, the corrected ! vector components rx and ry at rlt and rln can be calculated as ! follows: ! ! rad = acos(-1.)/180. ! call rotate (rlt, rln, -psir, -thetar, -phir, glt, gln) ! call get_vector (glt, gln, gx, gy) ! some routine to get vector ! r = sqrt(gx**2 + gy**2) ! if (r .gt. 0.) then ! a = gx/r ! a = min(a, 1.) ! a = max(a, -1.) ! a = acos(a) ! if (gy .lt. 0.) a = -a ! call rotvec(glt, gln, phir, thetar, psir, angle) ! a = a + angle*rad ! rx = r*cos(a) ! ry = r*sin(a) ! endif ! ! g... = geographic value ! r... = rotated value ! ...lt = latitude ! ...ln = longitude ! ...lth = latitude of head of vector ! ...lnh = longitude of head of vector ! dst = distance between heads of vectors ! delta = length of vector ! angle = angle to rotate vectors back to original orientation ! psir, thetar, phir = Euler angles defining rotation ! ! define multiplier rad for conversion to radians. rad = acos(-1.)/180. ! ! define length of direction vector (single precision may require a ! longer and thus less accurate vector length). delta = 0.001 ! ! if the base is in the north of the geographic grid use a south ! pointing vector to avoid any possibility of going over the pole. if (glt .ge. 0.) delta = -delta ! ! find the base of the direction vectors in the rotated grid. call rotate (glt, gln, phir, thetar, psir, rlt, rln) ! ! if base in the rotated grid is near a pole return an angle of zero if (abs(rlt) .ge. 90.-abs(delta)) then angle = 0. return endif ! ! find the head of the geographic grid direction vector in the ! rotated grid. call rotate (glt+delta, gln, phir, thetar, psir, glth, glnh) ! ! if the base is in opposite hemispheres switch the vector ! direction for better accuracy. if (glt*rlt .lt. 0) delta = -delta ! ! find the head of the rotated grid direction vector. rlth = rlt + delta rlnh = rln ! ! find the distance between the heads of the direction vectors. call dist (glth, glnh, rlth, rlnh, dst) ! ! find the angle between direction vectors with the law of cosines. delta = abs(delta) angle = (cos(dst)-cos(delta)**2)/(sin(delta)**2) angle = min(angle, 1.) angle = max(angle, -1.) angle = acos(angle)/rad t = abs(delta) ! ! adjust the angle if the direction vectors are opposite. if (glt*rlt .lt. 0) angle = 180. - angle ! determine the sign of the angle by checking the offset longitudes. if (glnh+360. .gt. rlnh+360.) angle = -angle ! change sign if the original direction vector was pointing south. if (glt .ge. 0.) angle = -angle ! return end subroutine rotvec subroutine dist (lat1, lng1, lat2, lng2, dst) ! !======================================================================= ! subroutine dist calculates the arc distance between two ! points given their latitudes and longitudes !======================================================================= ! real lat1, lng1, lat2, lng2, dst, rad ! ! define multiplier rad for conversion to radians. rad = acos(-1.)/180. ! ! check input. lat1 = min(lat1, 90.) lat1 = max(lat1, -90.) if (lng1 .lt. 0.) lng1 = lng1 + 360. if (lng1 .gt. 360.) lng1 = lng1 - 360. lat2 = min(lat2, 90.) lat2 = max(lat2, -90.) if (lng2 .lt. 0.) lng2 = lng2 + 360. if (lng2 .gt. 360.) lng2 = lng2 - 360. ! dst = sin(lat1*rad)*sin(lat2*rad)+cos(lat1*rad) & *cos(lat2*rad)*cos((lng1-lng2)*rad) dst = min(dst, 1.) dst = max(dst, -1.) dst = (acos(dst)/rad) return end subroutine dist end module rot_mod