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