[18] | 1 | !********************************************************************** |
---|
| 2 | ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * |
---|
| 3 | ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * |
---|
| 4 | ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * |
---|
| 5 | ! * |
---|
| 6 | ! This file is part of FLEXPART. * |
---|
| 7 | ! * |
---|
| 8 | ! FLEXPART is free software: you can redistribute it and/or modify * |
---|
| 9 | ! it under the terms of the GNU General Public License as published by* |
---|
| 10 | ! the Free Software Foundation, either version 3 of the License, or * |
---|
| 11 | ! (at your option) any later version. * |
---|
| 12 | ! * |
---|
| 13 | ! FLEXPART is distributed in the hope that it will be useful, * |
---|
| 14 | ! but WITHOUT ANY WARRANTY; without even the implied warranty of * |
---|
| 15 | ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * |
---|
| 16 | ! GNU General Public License for more details. * |
---|
| 17 | ! * |
---|
| 18 | ! You should have received a copy of the GNU General Public License * |
---|
| 19 | ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * |
---|
| 20 | !********************************************************************** |
---|
| 21 | |
---|
| 22 | ! Changes to the routines by A. Stohl |
---|
| 23 | ! xi,xi0,eta,eta0 are double precision variables to avoid problems |
---|
| 24 | ! at poles |
---|
| 25 | |
---|
| 26 | module cmapf_mod |
---|
| 27 | |
---|
| 28 | use par_mod, only: dp |
---|
| 29 | |
---|
| 30 | implicit none |
---|
| 31 | private |
---|
| 32 | |
---|
| 33 | public :: cc2gll, cll2xy, cgszll, cxy2ll, stlmbr, stcm2p |
---|
| 34 | |
---|
| 35 | real,parameter :: rearth=6371.2, almst1=.9999999 |
---|
| 36 | |
---|
| 37 | real,parameter :: pi=3.14159265358979 |
---|
| 38 | real,parameter :: radpdg=pi/180., dgprad=180./pi |
---|
| 39 | |
---|
| 40 | contains |
---|
| 41 | |
---|
| 42 | subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg) |
---|
| 43 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 44 | |
---|
| 45 | use par_mod, only: dp |
---|
| 46 | |
---|
| 47 | implicit none |
---|
| 48 | |
---|
| 49 | real :: strcmp(9), xlat, xlong, ue, vn, ug, vg |
---|
| 50 | |
---|
| 51 | real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot |
---|
| 52 | |
---|
| 53 | along = cspanf( xlong - strcmp(2), -180., 180.) |
---|
| 54 | if (xlat.gt.89.985) then |
---|
| 55 | !* North polar meteorological orientation: "north" along prime meridian |
---|
| 56 | rot = - strcmp(1) * along + xlong - 180. |
---|
| 57 | elseif (xlat.lt.-89.985) then |
---|
| 58 | !* South polar meteorological orientation: "north" along prime meridian |
---|
| 59 | rot = - strcmp(1) * along - xlong |
---|
| 60 | else |
---|
| 61 | rot = - strcmp(1) * along |
---|
| 62 | endif |
---|
| 63 | slong = sin( radpdg * rot ) |
---|
| 64 | clong = cos( radpdg * rot ) |
---|
| 65 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
| 66 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
| 67 | ug = ypolg * ue + xpolg * vn |
---|
| 68 | vg = ypolg * vn - xpolg * ue |
---|
| 69 | return |
---|
| 70 | end subroutine cc2gll |
---|
| 71 | |
---|
| 72 | subroutine ccrvll (strcmp, xlat,xlong, gx,gy) |
---|
| 73 | !* Written on 9/20/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 74 | |
---|
| 75 | use par_mod, only: dp |
---|
| 76 | |
---|
| 77 | implicit none |
---|
| 78 | |
---|
| 79 | real(kind=dp) :: xpolg,ypolg,temp,along,slong,clong,ctemp, curv |
---|
| 80 | real :: strcmp(9), xlat, xlong, gx, gy |
---|
| 81 | |
---|
| 82 | along = cspanf( xlong - strcmp(2), -180., 180.) |
---|
| 83 | slong = sin( radpdg * strcmp(1) * along) |
---|
| 84 | clong = cos( radpdg * strcmp(1) * along) |
---|
| 85 | xpolg = - slong * strcmp(5) + clong * strcmp(6) |
---|
| 86 | ypolg = clong * strcmp(5) + slong * strcmp(6) |
---|
| 87 | temp = sin(radpdg * xlat) |
---|
| 88 | ctemp = cos(radpdg * xlat) |
---|
| 89 | curv = (strcmp(1) - temp) / ctemp / rearth |
---|
| 90 | gx = curv * xpolg |
---|
| 91 | gy = curv * ypolg |
---|
| 92 | return |
---|
| 93 | end subroutine ccrvll |
---|
| 94 | |
---|
| 95 | subroutine ccrvxy (strcmp, x,y, gx,gy) |
---|
| 96 | !* Written on 9/20/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 97 | |
---|
| 98 | use par_mod, only: dp |
---|
| 99 | |
---|
| 100 | implicit none |
---|
| 101 | |
---|
| 102 | real :: strcmp(9), x, y, gx, gy |
---|
| 103 | real(kind=dp) :: xpolg,ypolg,temp,ymerc,efact,curv |
---|
| 104 | |
---|
| 105 | temp = strcmp(1) * strcmp(7) /rearth |
---|
| 106 | xpolg = strcmp(6) + temp * (strcmp(3) - x) |
---|
| 107 | ypolg = strcmp(5) + temp * (strcmp(4) - y) |
---|
| 108 | temp = sqrt ( xpolg ** 2 + ypolg ** 2 ) |
---|
| 109 | if (temp.gt.0.) then |
---|
| 110 | ymerc = - log( temp) /strcmp(1) |
---|
| 111 | efact = exp(ymerc) |
---|
| 112 | curv = ( (strcmp(1) - 1.d0) * efact + & |
---|
| 113 | (strcmp(1) + 1.d0) / efact ) & |
---|
| 114 | * .5d0 / rearth |
---|
| 115 | gx = xpolg * curv / temp |
---|
| 116 | gy = ypolg * curv / temp |
---|
| 117 | else |
---|
| 118 | if (abs(strcmp(1)) .eq. 1.) then |
---|
| 119 | gx = 0. |
---|
| 120 | gy = 0. |
---|
| 121 | else |
---|
| 122 | gx = 1./rearth |
---|
| 123 | gy = 1./rearth |
---|
| 124 | endif |
---|
| 125 | endif |
---|
| 126 | return |
---|
| 127 | end subroutine ccrvxy |
---|
| 128 | |
---|
| 129 | subroutine cg2cll (strcmp, xlat,xlong, ug,vg, ue,vn) |
---|
| 130 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 131 | |
---|
| 132 | use par_mod, only: dp |
---|
| 133 | |
---|
| 134 | implicit none |
---|
| 135 | |
---|
| 136 | real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot |
---|
| 137 | real :: strcmp(9), xlat, xlong, ug, vg, ue, vn |
---|
| 138 | |
---|
| 139 | along = cspanf( xlong - strcmp(2), -180., 180.) |
---|
| 140 | if (xlat.gt.89.985) then |
---|
| 141 | !* North polar meteorological orientation: "north" along prime meridian |
---|
| 142 | rot = - strcmp(1) * along + xlong - 180. |
---|
| 143 | elseif (xlat.lt.-89.985) then |
---|
| 144 | !* South polar meteorological orientation: "north" along prime meridian |
---|
| 145 | rot = - strcmp(1) * along - xlong |
---|
| 146 | else |
---|
| 147 | rot = - strcmp(1) * along |
---|
| 148 | endif |
---|
| 149 | slong = sin( radpdg * rot ) |
---|
| 150 | clong = cos( radpdg * rot ) |
---|
| 151 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
| 152 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
| 153 | ue = ypolg * ug - xpolg * vg |
---|
| 154 | vn = ypolg * vg + xpolg * ug |
---|
| 155 | return |
---|
| 156 | end subroutine cg2cll |
---|
| 157 | |
---|
| 158 | subroutine cg2cxy (strcmp, x,y, ug,vg, ue,vn) |
---|
| 159 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 160 | |
---|
| 161 | use par_mod, only: dp |
---|
| 162 | |
---|
| 163 | implicit none |
---|
| 164 | |
---|
| 165 | real :: strcmp(9) , x, y, ug, vg, ue, vn |
---|
| 166 | |
---|
| 167 | real :: clong, radial, rot, slong, xlat, xlong |
---|
| 168 | real(kind=dp) :: xpolg,ypolg,temp,xi0,eta0,xi,eta |
---|
| 169 | |
---|
| 170 | xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth |
---|
| 171 | eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth |
---|
| 172 | xi = xi0 * strcmp(5) - eta0 * strcmp(6) |
---|
| 173 | eta = eta0 * strcmp(5) + xi0 * strcmp(6) |
---|
| 174 | radial = 2. * eta - strcmp(1) * (xi*xi + eta*eta) |
---|
| 175 | if (radial.gt.strcmp(8)) then |
---|
| 176 | !* Case north of 89 degrees. Meteorological wind direction definition |
---|
| 177 | !* changes. |
---|
| 178 | call cnxyll(strcmp, xi,eta, xlat,xlong) |
---|
| 179 | !* North polar meteorological orientation: "north" along prime meridian |
---|
| 180 | rot = strcmp(1) * (xlong - strcmp(2)) - xlong - 180. |
---|
| 181 | slong = - sin( radpdg * rot ) |
---|
| 182 | clong = cos( radpdg * rot ) |
---|
| 183 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
| 184 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
| 185 | else if (radial.lt.strcmp(9)) then |
---|
| 186 | !* Case south of -89 degrees. Meteorological wind direction definition |
---|
| 187 | !* changes. |
---|
| 188 | call cnxyll(strcmp, xi,eta, xlat,xlong) |
---|
| 189 | !* South polar meteorological orientation: "north" along prime meridian |
---|
| 190 | rot = strcmp(1) * (xlong - strcmp(2)) + xlong |
---|
| 191 | slong = - sin( radpdg * rot ) |
---|
| 192 | clong = cos( radpdg * rot ) |
---|
| 193 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
| 194 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
| 195 | else |
---|
| 196 | !* Normal case. Meteorological direction of wind related to true north. |
---|
| 197 | xpolg = strcmp(6) - strcmp(1) * xi0 |
---|
| 198 | ypolg = strcmp(5) - strcmp(1) * eta0 |
---|
| 199 | temp = sqrt ( xpolg ** 2 + ypolg ** 2 ) |
---|
| 200 | xpolg = xpolg / temp |
---|
| 201 | ypolg = ypolg / temp |
---|
| 202 | end if |
---|
| 203 | ue = ( ypolg * ug - xpolg * vg ) |
---|
| 204 | vn = ( ypolg * vg + xpolg * ug ) |
---|
| 205 | return |
---|
| 206 | end subroutine cg2cxy |
---|
| 207 | |
---|
| 208 | real function cgszll (strcmp, xlat,xlong) |
---|
| 209 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 210 | |
---|
| 211 | use par_mod, only: dp |
---|
| 212 | |
---|
| 213 | implicit none |
---|
| 214 | |
---|
| 215 | real :: strcmp(9), xlat, xlong |
---|
| 216 | |
---|
| 217 | real(kind=dp) :: slat,ymerc,efact |
---|
| 218 | |
---|
| 219 | if (xlat .gt. 89.985) then |
---|
| 220 | !* Close to north pole |
---|
| 221 | if (strcmp(1) .gt. 0.9999) then |
---|
| 222 | !* and to gamma == 1. |
---|
| 223 | cgszll = 2. * strcmp(7) |
---|
| 224 | return |
---|
| 225 | endif |
---|
| 226 | efact = cos(radpdg * xlat) |
---|
| 227 | if (efact .le. 0.) then |
---|
| 228 | cgszll = 0. |
---|
| 229 | return |
---|
| 230 | else |
---|
| 231 | ymerc = - log( efact /(1. + sin(radpdg * xlat))) |
---|
| 232 | endif |
---|
| 233 | else if (xlat .lt. -89.985) then |
---|
| 234 | !* Close to south pole |
---|
| 235 | if (strcmp(1) .lt. -0.9999) then |
---|
| 236 | !* and to gamma == -1.0 |
---|
| 237 | cgszll = 2. * strcmp(7) |
---|
| 238 | return |
---|
| 239 | endif |
---|
| 240 | efact = cos(radpdg * xlat) |
---|
| 241 | if (efact .le. 0.) then |
---|
| 242 | cgszll = 0. |
---|
| 243 | return |
---|
| 244 | else |
---|
| 245 | ymerc = log( efact /(1. - sin(radpdg * xlat))) |
---|
| 246 | endif |
---|
| 247 | else |
---|
| 248 | slat = sin(radpdg * xlat) |
---|
| 249 | ymerc = log((1. + slat) / (1. - slat))/2. |
---|
| 250 | !efact = exp(ymerc) |
---|
| 251 | !cgszll = 2. * strcmp(7) * exp (strcmp(1) * ymerc) |
---|
| 252 | !c / (efact + 1./efact) |
---|
| 253 | endif |
---|
| 254 | cgszll = strcmp(7) * cos(radpdg * xlat) * exp(strcmp(1) *ymerc) |
---|
| 255 | return |
---|
| 256 | end function cgszll |
---|
| 257 | |
---|
| 258 | real function cgszxy (strcmp, x,y) |
---|
| 259 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 260 | |
---|
| 261 | use par_mod, only: dp |
---|
| 262 | |
---|
| 263 | implicit none |
---|
| 264 | |
---|
| 265 | real :: strcmp(9) , x, y |
---|
| 266 | real(kind=dp) :: ymerc,efact, radial, temp |
---|
| 267 | real(kind=dp) :: xi0,eta0,xi,eta |
---|
| 268 | |
---|
| 269 | |
---|
| 270 | xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth |
---|
| 271 | eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth |
---|
| 272 | xi = xi0 * strcmp(5) - eta0 * strcmp(6) |
---|
| 273 | eta = eta0 * strcmp(5) + xi0 * strcmp(6) |
---|
| 274 | radial = 2. * eta - strcmp(1) * (xi*xi + eta*eta) |
---|
| 275 | efact = strcmp(1) * radial |
---|
| 276 | if (efact .gt. almst1) then |
---|
| 277 | if (strcmp(1).gt.almst1) then |
---|
| 278 | cgszxy = 2. * strcmp(7) |
---|
| 279 | else |
---|
| 280 | cgszxy = 0. |
---|
| 281 | endif |
---|
| 282 | return |
---|
| 283 | endif |
---|
| 284 | if (abs(efact) .lt. 1.e-2) then |
---|
| 285 | temp = (efact / (2. - efact) )**2 |
---|
| 286 | ymerc = radial / (2. - efact) * (1. + temp * & |
---|
| 287 | (1./3. + temp * & |
---|
| 288 | (1./5. + temp * & |
---|
| 289 | (1./7. )))) |
---|
| 290 | else |
---|
| 291 | ymerc = - log( 1. - efact ) /2. /strcmp(1) |
---|
| 292 | endif |
---|
| 293 | if (ymerc .gt. 6.) then |
---|
| 294 | if (strcmp(1) .gt. almst1) then |
---|
| 295 | cgszxy = 2. * strcmp(7) |
---|
| 296 | else |
---|
| 297 | cgszxy = 0. |
---|
| 298 | endif |
---|
| 299 | else if (ymerc .lt. -6.) then |
---|
| 300 | if (strcmp(1) .lt. -almst1) then |
---|
| 301 | cgszxy = 2. * strcmp(7) |
---|
| 302 | else |
---|
| 303 | cgszxy = 0. |
---|
| 304 | endif |
---|
| 305 | else |
---|
| 306 | efact = exp(ymerc) |
---|
| 307 | cgszxy = 2. * strcmp(7) * exp (strcmp(1) * ymerc) & |
---|
| 308 | / (efact + 1./efact) |
---|
| 309 | endif |
---|
| 310 | return |
---|
| 311 | end function cgszxy |
---|
| 312 | |
---|
| 313 | subroutine cll2xy (strcmp, xlat,xlong, x,y) |
---|
| 314 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 315 | |
---|
| 316 | implicit none |
---|
| 317 | |
---|
| 318 | real :: strcmp(9) , xlat, xlong, x, y, xi, eta |
---|
| 319 | |
---|
| 320 | call cnllxy(strcmp, xlat,xlong, xi,eta) |
---|
| 321 | x = strcmp(3) + rearth/strcmp(7) * & |
---|
| 322 | (xi * strcmp(5) + eta * strcmp(6) ) |
---|
| 323 | y = strcmp(4) + rearth/strcmp(7) * & |
---|
| 324 | (eta * strcmp(5) - xi * strcmp(6) ) |
---|
| 325 | return |
---|
| 326 | end subroutine cll2xy |
---|
| 327 | |
---|
| 328 | subroutine cnllxy (strcmp, xlat,xlong, xi,eta) |
---|
| 329 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 330 | ! main transformation routine from latitude-longitude to |
---|
| 331 | ! canonical (equator-centered, radian unit) coordinates |
---|
| 332 | |
---|
| 333 | use par_mod, only: dp |
---|
| 334 | |
---|
| 335 | implicit none |
---|
| 336 | |
---|
| 337 | real :: strcmp(9), xlat, xlong, xi, eta, & |
---|
| 338 | gdlong, sndgam, csdgam, rhog1 |
---|
| 339 | real(kind=dp) :: gamma |
---|
| 340 | real(kind=dp) :: dlong,dlat,slat,mercy,gmercy |
---|
| 341 | |
---|
| 342 | gamma = strcmp(1) |
---|
| 343 | dlat = xlat |
---|
| 344 | dlong = cspanf(xlong - strcmp(2), -180., 180.) |
---|
| 345 | dlong = dlong * radpdg |
---|
| 346 | gdlong = gamma * dlong |
---|
| 347 | if (abs(gdlong) .lt. .01) then |
---|
| 348 | ! Code for gamma small or zero. This avoids round-off error or divide- |
---|
| 349 | ! by zero in the case of mercator or near-mercator projections. |
---|
| 350 | gdlong = gdlong * gdlong |
---|
| 351 | sndgam = dlong * (1. - 1./6. * gdlong * & |
---|
| 352 | (1. - 1./20. * gdlong * & |
---|
| 353 | (1. - 1./42. * gdlong ))) |
---|
| 354 | csdgam = dlong * dlong * .5 * & |
---|
| 355 | (1. - 1./12. * gdlong * & |
---|
| 356 | (1. - 1./30. * gdlong * & |
---|
| 357 | (1. - 1./56. * gdlong ))) |
---|
| 358 | else |
---|
| 359 | ! Code for moderate values of gamma |
---|
| 360 | sndgam = sin (gdlong) /gamma |
---|
| 361 | csdgam = (1. - cos(gdlong) )/gamma /gamma |
---|
| 362 | endif |
---|
| 363 | slat = sin(radpdg * dlat) |
---|
| 364 | if ((slat .ge. almst1) .or. (slat .le. -almst1)) then |
---|
| 365 | eta = 1./strcmp(1) |
---|
| 366 | xi = 0. |
---|
| 367 | return |
---|
| 368 | endif |
---|
| 369 | mercy = .5 * log( (1. + slat) / (1. - slat) ) |
---|
| 370 | gmercy = gamma * mercy |
---|
| 371 | if (abs(gmercy) .lt. .001) then |
---|
| 372 | ! Code for gamma small or zero. This avoids round-off error or divide- |
---|
| 373 | ! by zero in the case of mercator or near-mercator projections. |
---|
| 374 | rhog1 = mercy * (1. - .5 * gmercy * & |
---|
| 375 | (1. - 1./3. * gmercy * & |
---|
| 376 | (1. - 1./4. * gmercy ) ) ) |
---|
| 377 | else |
---|
| 378 | ! Code for moderate values of gamma |
---|
| 379 | rhog1 = (1. - exp(-gmercy)) / gamma |
---|
| 380 | endif |
---|
| 381 | eta = rhog1 + (1. - gamma * rhog1) * gamma * csdgam |
---|
| 382 | xi = (1. - gamma * rhog1 ) * sndgam |
---|
| 383 | end subroutine cnllxy |
---|
| 384 | |
---|
| 385 | subroutine cnxyll (strcmp, xi,eta, xlat,xlong) |
---|
| 386 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 387 | ! main transformation routine from canonical (equator-centered, |
---|
| 388 | ! radian unit) coordinates |
---|
| 389 | |
---|
| 390 | use par_mod, only: dp |
---|
| 391 | |
---|
| 392 | implicit none |
---|
| 393 | |
---|
| 394 | real :: strcmp(9), xlat, xlong, odist |
---|
| 395 | real(kind=dp) :: gamma,temp,arg1,arg2,ymerc,along,gxi,cgeta |
---|
| 396 | real(kind=dp) :: xi,eta |
---|
| 397 | |
---|
| 398 | gamma = strcmp(1) |
---|
| 399 | ! Calculate equivalent mercator coordinate |
---|
| 400 | odist = xi*xi + eta*eta |
---|
| 401 | arg2 = 2. * eta - gamma * (xi*xi + eta*eta) |
---|
| 402 | arg1 = gamma * arg2 |
---|
| 403 | ! Change by A. Stohl to avoid problems close to the poles |
---|
| 404 | ! if (arg1 .ge. almst1) then |
---|
| 405 | ! distance to north (or south) pole is zero (or imaginary ;) ) |
---|
| 406 | ! xlat = sign(90.,strcmp(1)) |
---|
| 407 | ! xlong = strcmp(2) |
---|
| 408 | ! return |
---|
| 409 | ! endif |
---|
| 410 | if (abs(arg1) .lt. .01) then |
---|
| 411 | ! Code for gamma small or zero. This avoids round-off error or divide- |
---|
| 412 | ! by zero in the case of mercator or near-mercator projections. |
---|
| 413 | temp = (arg1 / (2. - arg1) )**2 |
---|
| 414 | ymerc = arg2 / (2. - arg1) * (1. + temp * & |
---|
| 415 | (1./3. + temp * & |
---|
| 416 | (1./5. + temp * & |
---|
| 417 | (1./7. )))) |
---|
| 418 | else |
---|
| 419 | ! Code for moderate values of gamma |
---|
| 420 | ymerc = - log ( 1. - arg1 ) /2. / gamma |
---|
| 421 | endif |
---|
| 422 | ! Convert ymerc to latitude |
---|
| 423 | temp = exp( - abs(ymerc) ) |
---|
| 424 | xlat = sign(atan2((1. - temp) * (1. + temp), 2. * temp), ymerc) |
---|
| 425 | ! Find longitudes |
---|
| 426 | gxi = gamma*xi |
---|
| 427 | cgeta = 1. - gamma * eta |
---|
| 428 | if ( abs(gxi) .lt. .01*cgeta ) then |
---|
| 429 | ! Code for gamma small or zero. This avoids round-off error or divide- |
---|
| 430 | ! by zero in the case of mercator or near-mercator projections. |
---|
| 431 | temp = ( gxi /cgeta )**2 |
---|
| 432 | along = xi / cgeta * (1. - temp * & |
---|
| 433 | (1./3. - temp * & |
---|
| 434 | (1./5. - temp * & |
---|
| 435 | (1./7. )))) |
---|
| 436 | else |
---|
| 437 | ! Code for moderate values of gamma |
---|
| 438 | along = atan2( gxi, cgeta) / gamma |
---|
| 439 | endif |
---|
| 440 | xlong = sngl(strcmp(2) + dgprad * along) |
---|
| 441 | xlat = xlat * dgprad |
---|
| 442 | return |
---|
| 443 | end subroutine cnxyll |
---|
| 444 | |
---|
| 445 | subroutine cpolll (strcmp, xlat,xlong, enx,eny,enz) |
---|
| 446 | !* Written on 11/23/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 447 | |
---|
| 448 | use par_mod, only: dp |
---|
| 449 | |
---|
| 450 | implicit none |
---|
| 451 | |
---|
| 452 | real(kind=dp) :: xpolg,ypolg,along,slong,clong,rot |
---|
| 453 | real :: strcmp(9), xlat, xlong, enx, eny, enz, clat |
---|
| 454 | |
---|
| 455 | along = cspanf( xlong - strcmp(2), -180., 180.) |
---|
| 456 | rot = - strcmp(1) * along |
---|
| 457 | slong = sin( radpdg * rot ) |
---|
| 458 | clong = cos( radpdg * rot ) |
---|
| 459 | xpolg = slong * strcmp(5) + clong * strcmp(6) |
---|
| 460 | ypolg = clong * strcmp(5) - slong * strcmp(6) |
---|
| 461 | clat = cos(radpdg * xlat) |
---|
| 462 | enx = clat * xpolg |
---|
| 463 | eny = clat * ypolg |
---|
| 464 | enz = sin(radpdg * xlat) |
---|
| 465 | return |
---|
| 466 | end subroutine cpolll |
---|
| 467 | |
---|
| 468 | subroutine cpolxy (strcmp, x,y, enx,eny,enz) |
---|
| 469 | !* Written on 11/26/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 470 | |
---|
| 471 | use par_mod, only: dp |
---|
| 472 | |
---|
| 473 | implicit none |
---|
| 474 | |
---|
| 475 | real :: strcmp(9) , x, y, enx, eny, enz |
---|
| 476 | real(kind=dp) :: xpol,ypol,temp,xi0,eta0,xi,eta,radial |
---|
| 477 | real(kind=dp) :: temp2,ymerc,arg,oarg,clat |
---|
| 478 | |
---|
| 479 | xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth |
---|
| 480 | eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth |
---|
| 481 | xi = xi0 * strcmp(5) - eta0 * strcmp(6) |
---|
| 482 | eta = eta0 * strcmp(5) + xi0 * strcmp(6) |
---|
| 483 | radial = 2. * eta - strcmp(1) * (xi*xi + eta*eta) |
---|
| 484 | temp = strcmp(1) * radial |
---|
| 485 | if (temp .ge. 1.) then |
---|
| 486 | enx = 0. |
---|
| 487 | eny = 0. |
---|
| 488 | enz = sign(1.,strcmp(1)) |
---|
| 489 | return |
---|
| 490 | endif |
---|
| 491 | if (abs(temp).lt.1.e-2) then |
---|
| 492 | temp2 = (temp / (2. - temp))**2 |
---|
| 493 | ymerc = radial / (2. - temp) * (1. + temp2 * & |
---|
| 494 | (1./3. + temp2 * & |
---|
| 495 | (1./5. + temp2 * & |
---|
| 496 | (1./7.)))) |
---|
| 497 | else |
---|
| 498 | ymerc = -.5 * log(1. - temp) / strcmp(1) |
---|
| 499 | endif |
---|
| 500 | arg = exp( ymerc ) |
---|
| 501 | oarg = 1./arg |
---|
| 502 | clat = 2./(arg + oarg) |
---|
| 503 | enz = (arg - oarg) * clat /2. |
---|
| 504 | temp = clat / sqrt(1. - temp) |
---|
| 505 | xpol = - xi * strcmp(1) * temp |
---|
| 506 | ypol = (1. - eta * strcmp(1) ) * temp |
---|
| 507 | enx = xpol * strcmp(5) + ypol * strcmp(6) |
---|
| 508 | eny = ypol * strcmp(5) - xpol * strcmp(6) |
---|
| 509 | return |
---|
| 510 | end subroutine cpolxy |
---|
| 511 | |
---|
| 512 | real function cspanf (value, begin, end) |
---|
| 513 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 514 | !* real function cspanf returns a value in the interval (begin,end] |
---|
| 515 | !* which is equivalent to value, mod (end - begin). It is used to |
---|
| 516 | !* reduce periodic variables to a standard range. It adjusts for the |
---|
| 517 | !* behavior of the mod function which provides positive results for |
---|
| 518 | !* positive input, and negative results for negative input |
---|
| 519 | !* input: |
---|
| 520 | !* value - real number to be reduced to the span |
---|
| 521 | !* begin - first value of the span |
---|
| 522 | !* end - last value of the span |
---|
| 523 | !* returns: |
---|
| 524 | !* the reduced value |
---|
| 525 | !* examples: |
---|
| 526 | !* along = cspanf(xlong, -180., +180.) |
---|
| 527 | !* dir = cspanf(angle, 0., 360.) |
---|
| 528 | |
---|
| 529 | implicit none |
---|
| 530 | |
---|
| 531 | real :: first,last, value, begin, end, val |
---|
| 532 | |
---|
| 533 | first = min(begin,end) |
---|
| 534 | last = max(begin,end) |
---|
| 535 | val = mod( value - first , last - first) |
---|
| 536 | if ( val .le. 0.) then |
---|
| 537 | cspanf = val + last |
---|
| 538 | else |
---|
| 539 | cspanf = val + first |
---|
| 540 | endif |
---|
| 541 | return |
---|
| 542 | end function cspanf |
---|
| 543 | |
---|
| 544 | subroutine cxy2ll (strcmp, x,y, xlat,xlong) |
---|
| 545 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 546 | |
---|
| 547 | use par_mod, only: dp |
---|
| 548 | |
---|
| 549 | implicit none |
---|
| 550 | |
---|
| 551 | real(kind=dp) :: xi0,eta0,xi,eta |
---|
| 552 | real :: strcmp(9), x, y, xlat, xlong |
---|
| 553 | |
---|
| 554 | xi0 = ( x - strcmp(3) ) * strcmp(7) / rearth |
---|
| 555 | eta0 = ( y - strcmp(4) ) * strcmp(7) /rearth |
---|
| 556 | xi = xi0 * strcmp(5) - eta0 * strcmp(6) |
---|
| 557 | eta = eta0 * strcmp(5) + xi0 * strcmp(6) |
---|
| 558 | call cnxyll(strcmp, xi,eta, xlat,xlong) |
---|
| 559 | xlong = cspanf(xlong, -180., 180.) |
---|
| 560 | return |
---|
| 561 | end subroutine cxy2ll |
---|
| 562 | |
---|
| 563 | real function eqvlat (xlat1,xlat2) |
---|
| 564 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 565 | |
---|
| 566 | implicit none |
---|
| 567 | |
---|
| 568 | real :: xlat1, xlat2, x, ssind, sinl1, sinl2, al1, al2, tau |
---|
| 569 | |
---|
| 570 | ssind(x) = sin (radpdg*x) |
---|
| 571 | sinl1 = ssind (xlat1) |
---|
| 572 | sinl2 = ssind (xlat2) |
---|
| 573 | if (abs(sinl1 - sinl2) .gt. .001) then |
---|
| 574 | al1 = log((1. - sinl1)/(1. - sinl2)) |
---|
| 575 | al2 = log((1. + sinl1)/(1. + sinl2)) |
---|
| 576 | else |
---|
| 577 | ! Case lat1 near or equal to lat2 |
---|
| 578 | tau = - (sinl1 - sinl2)/(2. - sinl1 - sinl2) |
---|
| 579 | tau = tau*tau |
---|
| 580 | al1 = 2. / (2. - sinl1 - sinl2) * (1. + tau * & |
---|
| 581 | (1./3. + tau * & |
---|
| 582 | (1./5. + tau * & |
---|
| 583 | (1./7.)))) |
---|
| 584 | tau = (sinl1 - sinl2)/(2. + sinl1 + sinl2) |
---|
| 585 | tau = tau*tau |
---|
| 586 | al2 = -2. / (2. + sinl1 + sinl2) * (1. + tau * & |
---|
| 587 | (1./3. + tau * & |
---|
| 588 | (1./5. + tau * & |
---|
| 589 | (1./7.)))) |
---|
| 590 | endif |
---|
| 591 | eqvlat = asin((al1 + al2) / (al1 - al2))/radpdg |
---|
| 592 | return |
---|
| 593 | end function eqvlat |
---|
| 594 | |
---|
| 595 | subroutine stcm1p(strcmp, x1,y1, xlat1,xlong1, & |
---|
| 596 | xlatg,xlongg, gridsz, orient) |
---|
| 597 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 598 | |
---|
| 599 | implicit none |
---|
| 600 | |
---|
| 601 | integer :: k |
---|
| 602 | real :: strcmp(9), x1, y1, xlat1, xlong1, turn, orient, & |
---|
| 603 | xlatg, xlongg, gridsz, x1a, y1a |
---|
| 604 | |
---|
| 605 | do k=3,4 |
---|
| 606 | strcmp (k) = 0. |
---|
| 607 | enddo |
---|
| 608 | turn = radpdg * (orient - strcmp(1) * & |
---|
| 609 | cspanf(xlongg - strcmp(2), -180., 180.) ) |
---|
| 610 | strcmp (5) = cos (turn) |
---|
| 611 | strcmp (6) = - sin (turn) |
---|
| 612 | strcmp (7) = 1. |
---|
| 613 | strcmp (7) = gridsz * strcmp(7) & |
---|
| 614 | / cgszll(strcmp, xlatg, strcmp(2)) |
---|
| 615 | call cll2xy (strcmp, xlat1,xlong1, x1a,y1a) |
---|
| 616 | strcmp(3) = strcmp(3) + x1 - x1a |
---|
| 617 | strcmp(4) = strcmp(4) + y1 - y1a |
---|
| 618 | return |
---|
| 619 | end subroutine stcm1p |
---|
| 620 | |
---|
| 621 | subroutine stcm2p(strcmp, x1,y1, xlat1,xlong1, & |
---|
| 622 | x2,y2, xlat2,xlong2) |
---|
| 623 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 624 | |
---|
| 625 | implicit none |
---|
| 626 | |
---|
| 627 | real :: strcmp(9), x1, y1, xlat1, xlong1, & |
---|
| 628 | x2, y2, xlat2, xlong2 |
---|
| 629 | |
---|
| 630 | integer :: k |
---|
| 631 | real :: x1a, y1a, x2a, y2a, den, dena |
---|
| 632 | |
---|
| 633 | do k=3,6 |
---|
| 634 | strcmp (k) = 0. |
---|
| 635 | enddo |
---|
| 636 | strcmp (5) = 1. |
---|
| 637 | strcmp (7) = 1. |
---|
| 638 | call cll2xy (strcmp, xlat1,xlong1, x1a,y1a) |
---|
| 639 | call cll2xy (strcmp, xlat2,xlong2, x2a,y2a) |
---|
| 640 | den = sqrt( (x1 - x2)**2 + (y1 - y2)**2 ) |
---|
| 641 | dena = sqrt( (x1a - x2a)**2 + (y1a - y2a)**2 ) |
---|
| 642 | strcmp(5) = ((x1a - x2a)*(x1 - x2) + (y1a - y2a) * (y1 - y2)) & |
---|
| 643 | /den /dena |
---|
| 644 | strcmp(6) = ((y1a - y2a)*(x1 - x2) - (x1a - x2a) * (y1 - y2)) & |
---|
| 645 | /den /dena |
---|
| 646 | strcmp (7) = strcmp(7) * dena / den |
---|
| 647 | call cll2xy (strcmp, xlat1,xlong1, x1a,y1a) |
---|
| 648 | strcmp(3) = strcmp(3) + x1 - x1a |
---|
| 649 | strcmp(4) = strcmp(4) + y1 - y1a |
---|
| 650 | return |
---|
| 651 | end subroutine stcm2p |
---|
| 652 | |
---|
| 653 | !* General conformal map routines for meteorological modelers |
---|
| 654 | !* written on 3/31/94 by |
---|
| 655 | |
---|
| 656 | !* Dr. Albion Taylor |
---|
| 657 | !* NOAA / OAR / ARL phone: (301) 713-0295 x 132 |
---|
| 658 | !* rm. 3151, 1315 east-west highway fax: (301) 713-0119 |
---|
| 659 | !* silver spring, md 20910 e-mail: adtaylor@arlrisc.ssmc.noaa.gov |
---|
| 660 | |
---|
| 661 | !* subroutine stlmbr (strcmp, tnglat, clong) |
---|
| 662 | !* This routine initializes the map structure array strcmp to |
---|
| 663 | !* the form of a specific map projection |
---|
| 664 | !* inputs: |
---|
| 665 | !* tnglat - the latitude at which the projection will be tangent |
---|
| 666 | !* to the earth. +90. For north polar stereographic, |
---|
| 667 | !* -90. for south polar stereographic, 0. For mercator, |
---|
| 668 | !* and other values for lambert conformal. |
---|
| 669 | !* -90 <= tnglat <= 90. |
---|
| 670 | !* clong - a longitude in the region under consideration. Longitudes |
---|
| 671 | !* between clong-180. and clong+180. Will be mapped in one |
---|
| 672 | !* connected region |
---|
| 673 | !* outputs: |
---|
| 674 | !* strcmp - a 9-value map structure array for use with subsequent |
---|
| 675 | !* calls to the coordinate transform routines. |
---|
| 676 | !* |
---|
| 677 | !* real function eqvlat (xlat1,xlat2) |
---|
| 678 | !* This function is provided to assist in finding the tangent latitude |
---|
| 679 | !* equivalent to the 2-reference latitude specification in the legend |
---|
| 680 | !* of most lambert conformal maps. If the map specifies "scale |
---|
| 681 | !* 1:xxxxx true at 40n and 60n", then eqvlat(40.,60.) will return the |
---|
| 682 | !* equivalent tangent latitude. |
---|
| 683 | !* inputs: |
---|
| 684 | !* xlat1,xlat2: the two latitudes specified in the map legend |
---|
| 685 | !* returns: |
---|
| 686 | !* the equivalent tangent latitude |
---|
| 687 | !* example: call stlmbr(strcmp, eqvlat(40.,60.), 90.) |
---|
| 688 | |
---|
| 689 | !* subroutine stcm2p (strcmp, x1,y1, xlat1,xlong1, |
---|
| 690 | !* x2,y2, xlat2,xlong2) |
---|
| 691 | !* subroutine stcm1p (strcmp, x1,y1, xlat1,xlong1, |
---|
| 692 | !* xlatg,xlongg, gridsz, orient) |
---|
| 693 | !* These routines complete the specification of the map structure |
---|
| 694 | !* array by conforming the map coordinates to the specifications |
---|
| 695 | !* of a particular grid. Either stcm1p or stcm2p must be called, |
---|
| 696 | !* but not both |
---|
| 697 | !* inputs: |
---|
| 698 | !* strcmp - a 9-value map structure array, set to a particular map |
---|
| 699 | !* form by a previous call to stlmbr |
---|
| 700 | !* for stcm2p: |
---|
| 701 | !* x1,y1, x2,y2 - the map coordinates of two points on the grid |
---|
| 702 | !* xlat1,xlong1, xlat2,xlong2 - the geographic coordinates of the |
---|
| 703 | !* same two points |
---|
| 704 | !* for stcm1p: |
---|
| 705 | !* x1,y1 - the map coordinates of one point on the grid |
---|
| 706 | !* xlat1,xlong1 - the geographic coordinates of the same point |
---|
| 707 | !* xlatg,xlongg - latitude and longitude of reference point for |
---|
| 708 | !* gridsz and orientation specification. |
---|
| 709 | !* gridsz - the desired grid size in kilometers, at xlatg,xlongg |
---|
| 710 | !* orient - the angle, with respect to north, of a y-grid line, at |
---|
| 711 | !* the point xlatg,xlongg |
---|
| 712 | !* outputs: |
---|
| 713 | !* strcmp - a 9-value map structure array, fully set for use by |
---|
| 714 | !* other subroutines in this system |
---|
| 715 | |
---|
| 716 | !* subroutine cll2xy (strcmp, xlat,xlong, x,y) |
---|
| 717 | !* subroutine cxy2ll (strcmp, x,y, xlat,xlong) |
---|
| 718 | !* these routines convert between map coordinates x,y |
---|
| 719 | !* and geographic coordinates xlat,xlong |
---|
| 720 | !* inputs: |
---|
| 721 | !* strcmp(9) - 9-value map structure array |
---|
| 722 | !* for cll2xy: xlat,xlong - geographic coordinates |
---|
| 723 | !* for cxy2ll: x,y - map coordinates |
---|
| 724 | !* outputs: |
---|
| 725 | !* for cll2xy: x,y - map coordinates |
---|
| 726 | !* for cxy2ll: xlat,xlong - geographic coordinates |
---|
| 727 | |
---|
| 728 | !* subroutine cc2gxy (strcmp, x,y, ue,vn, ug,vg) |
---|
| 729 | !* subroutine cg2cxy (strcmp, x,y, ug,vg, ue,vn) |
---|
| 730 | !* subroutine cc2gll (strcmp, xlat,xlong, ue,vn, ug,vg) |
---|
| 731 | !* subroutine cg2cll (strcmp, xlat,xlong, ug,vg, ue,vn) |
---|
| 732 | !* These subroutines convert vector wind components from |
---|
| 733 | !* geographic, or compass, coordinates, to map or |
---|
| 734 | !* grid coordinates. The site of the wind to be |
---|
| 735 | !* converted may be given either in geographic or |
---|
| 736 | !* map coordinates. Wind components are all in kilometers |
---|
| 737 | !* per hour, whether geographic or map coordinates. |
---|
| 738 | !* inputs: |
---|
| 739 | !* strcmp(9) - 9-value map structure array |
---|
| 740 | !* for cc2gxy and cg2cxy: x,y - map coordinates of site |
---|
| 741 | !* for cc2gll and cg2cll: xlat,xlong - geographic coordinates of site |
---|
| 742 | !* for cc2gxy and cc2gll: ue,vn - east and north wind components |
---|
| 743 | !* for cg2cxy and cg2cll: ug,vg - x- and y- direction wind components |
---|
| 744 | !* outputs: |
---|
| 745 | !* for cc2gxy and cc2gll: ug,vg - x- and y- direction wind components |
---|
| 746 | !* for cg2cxy and cg2cll: ue,vn - east and north wind components |
---|
| 747 | |
---|
| 748 | !* subroutine ccrvxy (strcmp, x, y, gx,gy) |
---|
| 749 | !* subroutine ccrvll (strcmp, xlat,xlong, gx,gy) |
---|
| 750 | !* These subroutines return the curvature vector (gx,gy), as referenced |
---|
| 751 | !* to map coordinates, induced by the map transformation. When |
---|
| 752 | !* non-linear terms in wind speed are important, a "geodesic" force |
---|
| 753 | !* should be included in the vector form [ (u,u) g - (u,g) u ] where the |
---|
| 754 | !* inner product (u,g) is defined as ux*gx + uy*gy. |
---|
| 755 | !* inputs: |
---|
| 756 | !* strcmp(9) - 9-value map structure array |
---|
| 757 | !* for ccrvxy: x,y - map coordinates of site |
---|
| 758 | !* for ccrvll: xlat,xlong - geographic coordinates of site |
---|
| 759 | !* outputs: |
---|
| 760 | !* gx,gy - vector coefficients of curvature, in units radians |
---|
| 761 | !* per kilometer |
---|
| 762 | |
---|
| 763 | !* real function cgszll (strcmp, xlat,xlong) |
---|
| 764 | !* real function cgszxy (strcmp, x,y) |
---|
| 765 | !* These functions return the size, in kilometers, of each unit of |
---|
| 766 | !* motion in map coordinates (grid size). The grid size at any |
---|
| 767 | !* location depends on that location; the position may be given in |
---|
| 768 | !* either map or geographic coordinates. |
---|
| 769 | !* inputs: |
---|
| 770 | !* strcmp(9) - 9-value map structure array |
---|
| 771 | !* for cgszxy: x,y - map coordinates of site |
---|
| 772 | !* for cgszll: xlat,xlong - geographic coordinates of site |
---|
| 773 | !* returns: |
---|
| 774 | !* gridsize in kilometers at given site. |
---|
| 775 | |
---|
| 776 | !* subroutine cpolxy (strcmp, x,y, enx,eny,enz) |
---|
| 777 | !* subroutine cpolll (strcmp, xlat,xlong, enx,eny,enz) |
---|
| 778 | !* These subroutines provide 3-d vector components of a unit vector |
---|
| 779 | !* in the direction of the north polar axis. When multiplied |
---|
| 780 | !* by twice the rotation rate of the earth (2 * pi/24 hr), the |
---|
| 781 | !* vertical component yields the coriolis factor. |
---|
| 782 | !* inputs: |
---|
| 783 | !* strcmp(9) - 9-value map structure array |
---|
| 784 | !* for cpolxy: x,y - map coordinates of site |
---|
| 785 | !* for cpolll: xlat,xlong - geographic coordinates of site |
---|
| 786 | !* returns: |
---|
| 787 | !* enx,eny,enz the direction cosines of a unit vector in the |
---|
| 788 | !* direction of the rotation axis of the earth |
---|
| 789 | |
---|
| 790 | !* subroutine cnllxy (strcmp, xlat,xlong, xi,eta) |
---|
| 791 | !* subroutine cnxyll (strcmp, xi,eta, xlat,xlong) |
---|
| 792 | !* These subroutines perform the underlying transformations from |
---|
| 793 | !* geographic coordinates to and from canonical (equator centered) |
---|
| 794 | !* coordinates. They are called by cxy2ll and cll2xy, but are not |
---|
| 795 | !* intended to be called directly |
---|
| 796 | |
---|
| 797 | !* real function cspanf (value, begin, end) |
---|
| 798 | !* This function assists other routines in providing a longitude in |
---|
| 799 | !* the proper range. It adds to value whatever multiple of |
---|
| 800 | !* (end - begin) is needed to return a number begin < cspanf <= end |
---|
| 801 | |
---|
| 802 | subroutine stlmbr(strcmp, tnglat, xlong) |
---|
| 803 | !* Written on 3/31/94 by Dr. Albion Taylor NOAA / OAR / ARL |
---|
| 804 | |
---|
| 805 | implicit none |
---|
| 806 | |
---|
| 807 | real :: strcmp(9), tnglat, xlong |
---|
| 808 | |
---|
| 809 | real :: eta, xi |
---|
| 810 | |
---|
| 811 | strcmp(1) = sin(radpdg * tnglat) |
---|
| 812 | !* gamma = sine of the tangent latitude |
---|
| 813 | strcmp(2) = cspanf( xlong, -180., +180.) |
---|
| 814 | !* lambda_0 = reference longitude |
---|
| 815 | strcmp(3) = 0. |
---|
| 816 | !* x_0 = x- grid coordinate of origin (xi,eta) = (0.,0.) |
---|
| 817 | strcmp(4) = 0. |
---|
| 818 | !* y_0 = y-grid coordinate of origin (xi,eta) = (0.,0.) |
---|
| 819 | strcmp(5) = 1. |
---|
| 820 | !* Cosine of rotation angle from xi,eta to x,y |
---|
| 821 | strcmp(6) = 0. |
---|
| 822 | !* Sine of rotation angle from xi,eta to x,y |
---|
| 823 | strcmp(7) = rearth |
---|
| 824 | !* Gridsize in kilometers at the equator |
---|
| 825 | call cnllxy(strcmp, 89.,xlong, xi,eta) |
---|
| 826 | strcmp(8) = 2. * eta - strcmp(1) * eta * eta |
---|
| 827 | !* Radial coordinate for 1 degree from north pole |
---|
| 828 | call cnllxy(strcmp, -89.,xlong, xi,eta) |
---|
| 829 | strcmp(9) = 2. * eta - strcmp(1) * eta * eta |
---|
| 830 | !* Radial coordinate for 1 degree from south pole |
---|
| 831 | return |
---|
| 832 | end subroutine stlmbr |
---|
| 833 | |
---|
| 834 | end module cmapf_mod |
---|