var/da/da_verif_obs/da_verif_tools.f90

   1
   2
   3 module da_verif_tools
   4
   5
   6    implicit none
   7
   8
   9    ! Define some private constants
  10
  11    real, private, parameter    :: pi = 3.1415926
  12    real, private, parameter    :: deg_per_rad = 180.0/pi
  13    real, private, parameter    :: rad_per_deg = pi / 180.0
  14    logical, private, parameter :: print_detail_map = .false.
  15
  16    integer, private, parameter :: stdout = 6
  17
  18    ! Mean Earth Radius in m.  The value below is consistent
  19    ! with NCEP's routines and grids.
  20    ! real, public , parameter    :: earth_radius_m = 6371200.0 ! from Brent
  21    real, public , parameter    :: earth_radius_m = 6370000.0
  22    real, public , parameter    :: radians_per_degree = pi / 180.0
  23
  24    ! Define public parameters
  25
  26    ! Projection codes for proj_info structure:
  27    integer, public, parameter  :: PROJ_LATLON = 0
  28    integer, public, parameter  :: PROJ_MERC = 1
  29    integer, public, parameter  :: PROJ_LC = 3
  30    integer, public, parameter  :: PROJ_PS = 5
  31
  32
  33   ! Define data structures to define various projections
  34
  35    type proj_info
  36       integer          :: code     ! integer code for projection type
  37       real             :: lat1     ! SW latitude (1,1) in degrees (-90->90N)
  38       real             :: lon1     ! SW longitude (1,1) in degrees (-180->180E)
  39       real             :: dx       ! Grid spacing in meters at truelats, used
  40                                    ! only for ps, lc, and merc projections
  41       real             :: dlat     ! Lat increment for lat/lon grids
  42       real             :: dlon     ! Lon increment for lat/lon grids
  43       real             :: stdlon   ! Longitude parallel to y-axis (-180->180E)
  44       real             :: truelat1 ! First true latitude (all projections)
  45       real             :: truelat2 ! Second true lat (LC only)
  46       real             :: hemi     ! 1 for NH, -1 for SH
  47       real             :: cone     ! Cone factor for LC projections
  48       real             :: polei    ! Computed i-location of pole point
  49       real             :: polej    ! Computed j-location of pole point
  50       real             :: rsw      ! Computed radius to SW corner
  51       real             :: rebydx   ! Earth radius divided by dx
  52       real             :: knowni   ! X-location of known lat/lon
  53       real             :: knownj   ! Y-location of known lat/lon
  54       real             :: latinc   ! Latitude increments in degrees
  55       real             :: loninc   ! Longitude increments in degrees
  56       logical          :: init     ! Flag to indicate if this struct is
  57                                  ! ready for use
  58    end type proj_info
  59
  60    type(proj_info) :: map_info
  61
  62    character(len=10000),private :: message(50)
  63
  64
  65 contains
  66
  67 subroutine da_llxy_latlon(lat, lon, proj, x, y)
  68
  69    !-----------------------------------------------------------------------
  70    ! Purpose: Compute the x/y location of a lat/lon on a LATLON grid.
  71    !-----------------------------------------------------------------------
  72
  73    implicit none
  74
  75    real, intent(in)             :: lat
  76    real, intent(in)             :: lon
  77    type(proj_info), intent(in)  :: proj
  78    real, intent(out)            :: x
  79    real, intent(out)            :: y
  80
  81    real                         :: deltalat
  82    real                         :: deltalon
  83    real                         :: lon360
  84    real                         :: latinc
  85    real                         :: loninc
  86
  87    ! Extract the latitude and longitude increments for this grid
  88    ! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
  89    ! loninc is saved in the stdlon tag and latinc is saved in truelat1
  90
  91    latinc = proj%truelat1
  92    loninc = proj%stdlon
  93
  94    ! Compute deltalat and deltalon as the difference between the input
  95    ! lat/lon and the origin lat/lon
  96
  97    deltalat = lat - proj%lat1
  98
  99    ! To account for issues around the dateline, convert the incoming
 100    ! longitudes to be 0->360.0
 101    if (lon < 0) then
 102       lon360 = lon + 360.0
 103    else
 104       lon360 = lon
 105    end if
 106    deltalon = lon360 - proj%lon1
 107
 108    ! Compute x/y
 109    x = deltalon/loninc + 1.0
 110    y = deltalat/latinc + 1.0
 111
 112 end subroutine da_llxy_latlon
 113
 114
 115 subroutine da_llxy_lc(lat, lon, proj, x, y)
 116
 117    !-----------------------------------------------------------------------
 118    ! Purpose: compute the geographical latitude and longitude values
 119    ! to the cartesian x/y on a Lambert Conformal projection.
 120    !-----------------------------------------------------------------------
 121
 122    implicit none
 123
 124    real, intent(in)              :: lat      ! Latitude (-90->90 deg N)
 125    real, intent(in)              :: lon      ! Longitude (-180->180 E)
 126    type(proj_info),intent(in)    :: proj     ! Projection info structure
 127
 128    real, intent(out)             :: x        ! Cartesian X coordinate
 129    real, intent(out)             :: y        ! Cartesian Y coordinate
 130
 131    real                          :: arg
 132    real                          :: deltalon
 133    real                          :: tl1r
 134    real                          :: rm
 135    real                          :: ctl1r
 136
 137    ! Compute deltalon between known longitude and standard lon and ensure
 138    ! it is not in the cut zone
 139    deltalon = lon - proj%stdlon
 140    if (deltalon > +180.0) deltalon = deltalon - 360.0
 141    if (deltalon < -180.0) deltalon = deltalon + 360.0
 142
 143    ! Convert truelat1 to radian and compute COS for later use
 144    tl1r = proj%truelat1 * rad_per_deg
 145    ctl1r = COS(tl1r)
 146
 147    ! Radius to desired point
 148    rm = proj%rebydx * ctl1r/proj%cone * &
 149        (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / &
 150         TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0))**proj%cone
 151
 152    arg = proj%cone*(deltalon*rad_per_deg)
 153    x = proj%polei + proj%hemi * rm * Sin(arg)
 154    y = proj%polej - rm * COS(arg)
 155
 156    ! Finally, if we are in the southern hemisphere, flip the i/j
 157    ! values to a coordinate system where (1,1) is the SW corner
 158    ! (what we assume) which is different than the original NCEP
 159    ! algorithms which used the NE corner as the origin in the
 160    ! southern hemisphere (left-hand vs. right-hand coordinate?)
 161    if (proj%hemi == -1.0) then
 162       x = 2.0 - x
 163       y = 2.0 - y
 164    end if
 165
 166 end subroutine da_llxy_lc
 167
 168
 169 subroutine da_llxy_merc(lat, lon, proj, x, y)
 170
 171    !-----------------------------------------------------------------------
 172    ! Purpose: Compute x,y coordinate from lat lon for mercator projection
 173    !-----------------------------------------------------------------------
 174
 175    implicit none
 176
 177    real, intent(in)              :: lat
 178    real, intent(in)              :: lon
 179    type(proj_info),intent(in)    :: proj
 180    real,intent(out)              :: x
 181    real,intent(out)              :: y
 182    real                          :: deltalon
 183
 184    deltalon = lon - proj%lon1
 185    if (deltalon < -180.0) deltalon = deltalon + 360.0
 186    if (deltalon > 180.0) deltalon = deltalon - 360.0
 187    x = 1.0 + (deltalon/(proj%dlon*deg_per_rad))
 188    y = 1.0 + (ALOG(TAN(0.5*((lat + 90.0) * rad_per_deg)))) / &
 189            proj%dlon - proj%rsw
 190
 191 end subroutine da_llxy_merc
 192
 193
 194 subroutine da_llxy_ps(lat,lon,proj,x,y)
 195
 196    !-----------------------------------------------------------------------
 197    ! Purpose: Given latitude (-90 to 90), longitude (-180 to 180), and the
 198    ! standard polar-stereographic projection information via the
 199    ! public proj structure, this routine returns the x/y indices which
 200    ! if within the domain range from 1->nx and 1->ny, respectively.
 201    !-----------------------------------------------------------------------
 202
 203    implicit none
 204
 205    real, intent(in)               :: lat
 206    real, intent(in)               :: lon
 207    type(proj_info),intent(in)     :: proj
 208
 209    real, intent(out)              :: x !(x-index)
 210    real, intent(out)              :: y !(y-index)
 211
 212    real                           :: reflon
 213    real                           :: scale_top
 214    real                           :: ala
 215    real                           :: alo
 216    real                           :: rm
 217
 218    reflon = proj%stdlon + 90.0
 219
 220    ! Compute numerator term of map scale factor
 221
 222    scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
 223
 224    ! Find radius to desired point
 225    ala = lat * rad_per_deg
 226    rm = proj%rebydx * COS(ala) * scale_top/(1.0 + proj%hemi *Sin(ala))
 227    alo = (lon - reflon) * rad_per_deg
 228    x = proj%polei + rm * COS(alo)
 229    y = proj%polej + proj%hemi * rm * Sin(alo)
 230
 231 end subroutine da_llxy_ps
 232
 233
 234 subroutine da_llxy_wrf(proj, lat, lon, x, y)
 235
 236    !-----------------------------------------------------------------------
 237    ! Purpose: Converts input lat/lon values to the cartesian (x, y) value
 238    ! for the given projection.
 239    !-----------------------------------------------------------------------
 240
 241    implicit none
 242
 243    type(proj_info), intent(in)  :: proj
 244    real,            intent(in)  :: lat
 245    real,            intent(in)  :: lon
 246    real,            intent(out) :: x
 247    real,            intent(out) :: y
 248
 249    if (.NOT.proj%init) then
 250       write(stdout,*) "da_llxy_wrf.inc"// &
 251            "You have not called map_set for this projection!"
 252       stop
 253    end if
 254
 255    select case(proj%code)
 256
 257       case(PROJ_LATLON)
 258          call da_llxy_latlon(lat,lon,proj,x,y)
 259
 260       case(PROJ_MERC)
 261          call da_llxy_merc(lat,lon,proj,x,y)
 262          x = x + proj%knowni - 1.0
 263          y = y + proj%knownj - 1.0
 264
 265       case(PROJ_PS)
 266          call da_llxy_ps(lat,lon,proj,x,y)
 267
 268       case(PROJ_LC)
 269          call da_llxy_lc(lat,lon,proj,x,y)
 270          x = x + proj%knowni - 1.0
 271          y = y + proj%knownj - 1.0
 272
 273       case default
 274          write(unit=message(1),fmt='(A,I2)') &
 275             'Unrecognized map projection code: ', proj%code
 276          write(stdout,*) "da_llxy_wrf.inc"//message(1:1)
 277          stop
 278    end select
 279
 280 end subroutine da_llxy_wrf
 281
 282
 283 subroutine da_xyll(proj, xx, yy, lat, lon)
 284
 285    !-----------------------------------------------------------------------
 286    ! Purpose: Computes geographical latitude and longitude for a given (i,j)
 287    ! point in a grid with a projection of proj
 288    !-----------------------------------------------------------------------
 289
 290    implicit none
 291
 292    type(proj_info), intent(in)  :: proj
 293    real,            intent(in)  :: xx
 294    real,            intent(in)  :: yy
 295    real,            intent(out) :: lat
 296    real,            intent(out) :: lon
 297
 298    real :: x, y
 299
 300    if (.NOT.proj%init) then
 301       write(stdout,*) "da_xyll.inc"// &
 302            "You have not called map_set for this projection!"
 303       stop
 304    end if
 305
 306    x = xx
 307    y = yy
 308
 309    select case (proj%code)
 310       case (PROJ_LATLON)
 311          call da_xyll_latlon(x, y, proj, lat, lon)
 312
 313       case (PROJ_MERC)
 314          x = xx - proj%knowni + 1.0
 315          y = yy - proj%knownj + 1.0
 316          call da_xyll_merc(x, y, proj, lat, lon)
 317
 318       case (PROJ_PS)
 319          call da_xyll_ps(x, y, proj, lat, lon)
 320
 321       case (PROJ_LC)
 322
 323          x = xx - proj%knowni + 1.0
 324          y = yy - proj%knownj + 1.0
 325          call da_xyll_lc(x, y, proj, lat, lon)
 326
 327       case default
 328          write(unit=message(1),fmt='(A,I2)') &
 329             "Unrecognized map projection code: ", proj%code
 330          write(stdout,*) "da_xyll.inc"//message(1:1)
 331          stop
 332
 333    end select
 334
 335 end subroutine da_xyll
 336
 337
 338 subroutine da_xyll_latlon(x, y, proj, lat, lon)
 339
 340    !-----------------------------------------------------------------------
 341    ! Purpose: Compute the lat/lon location of an x/y on a LATLON grid.
 342    !-----------------------------------------------------------------------
 343
 344    implicit none
 345
 346    real, intent(in)             :: x
 347    real, intent(in)             :: y
 348    type(proj_info), intent(in)  :: proj
 349    real, intent(out)            :: lat
 350    real, intent(out)            :: lon
 351
 352    real                         :: deltalat
 353    real                         :: deltalon
 354    real                         :: latinc
 355    real                         :: loninc
 356
 357    ! Extract the latitude and longitude increments for this grid
 358    ! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
 359    ! loninc is saved in the stdlon tag and latinc is saved in truelat1
 360
 361    latinc = proj%truelat1
 362    loninc = proj%stdlon
 363
 364    ! Compute deltalat and deltalon
 365
 366    deltalat = (x-1.0)*latinc
 367    deltalon = (y-1.0)*loninc
 368    lat = proj%lat1 + deltalat
 369    lon = proj%lon1 + deltalon
 370
 371    if ((ABS(lat) > 90.0).OR.(ABS(deltalon) > 360.0)) then
 372       ! Off the earth for this grid
 373       lat = -999.0
 374       lon = -999.0
 375    else
 376       lon = lon + 360.0
 377       lon = MOD(lon,360.0)
 378       if (lon > 180.0) lon = lon -360.0
 379    end if
 380
 381 end subroutine da_xyll_latlon
 382
 383
 384 subroutine da_xyll_lc(x, y, proj, lat, lon)
 385
 386    ! subroutine da_to convert from the (x,y) cartesian coordinate to the
 387    ! geographical latitude and longitude for a Lambert Conformal projection.
 388
 389    implicit none
 390
 391    real, intent(in)              :: x        ! Cartesian X coordinate
 392    real, intent(in)              :: y        ! Cartesian Y coordinate
 393    type(proj_info),intent(in)    :: proj     ! Projection info structure
 394
 395
 396    real, intent(out)             :: lat      ! Latitude (-90->90 deg N)
 397    real, intent(out)             :: lon      ! Longitude (-180->180 E)
 398
 399    real                          :: inew
 400    real                          :: jnew
 401    real                          :: r
 402    real                          :: chi,chi1,chi2
 403    real                          :: r2
 404    real                          :: xx
 405    real                          :: yy
 406
 407    chi1 = (90.0 - proj%hemi*proj%truelat1)*rad_per_deg
 408    chi2 = (90.0 - proj%hemi*proj%truelat2)*rad_per_deg
 409
 410    ! See if we are in the southern hemispere and flip the indices
 411    ! if we are.
 412    if (proj%hemi == -1.0) then
 413       inew = -x + 2.0
 414       jnew = -y + 2.0
 415    else
 416       inew = x
 417       jnew = y
 418    end if
 419
 420    ! Compute radius**2 to i/j location
 421    xx = inew - proj%polei
 422    yy = proj%polej - jnew
 423    r2 = (xx*xx + yy*yy)
 424    r = sqrt(r2)/proj%rebydx
 425
 426    ! Convert to lat/lon
 427    if (r2 == 0.0) then
 428       lat = proj%hemi * 90.0
 429       lon = proj%stdlon
 430    else
 431
 432       ! Longitude
 433       lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
 434       lon = MOD(lon+360.0, 360.0)
 435
 436       ! Latitude.  Latitude determined by solving an equation adapted
 437       ! from:
 438       !  Maling, D.H., 1973: Coordinate Systems and Map Projections
 439       ! Equations #20 in Appendix I.
 440
 441       if (chi1 == chi2) then
 442          chi = 2.0*ATAN((r/TAN(chi1))**(1.0/proj%cone) * TAN(chi1*0.5))
 443       else
 444          chi = 2.0*ATAN((r*proj%cone/Sin(chi1))**(1.0/proj%cone) * TAN(chi1*0.5))
 445       end if
 446       lat = (90.0-chi*deg_per_rad)*proj%hemi
 447    end if
 448
 449    if (lon > +180.0) lon = lon - 360.0
 450    if (lon < -180.0) lon = lon + 360.0
 451
 452 end subroutine da_xyll_lc
 453
 454
 455 subroutine da_xyll_merc(x, y, proj, lat, lon)
 456
 457    !-----------------------------------------------------------------------
 458    ! Compute the lat/lon from i/j for mercator projection
 459    !-----------------------------------------------------------------------
 460
 461    implicit none
 462
 463    real,intent(in)               :: x
 464    real,intent(in)               :: y
 465    type(proj_info),intent(in)    :: proj
 466    real, intent(out)             :: lat
 467    real, intent(out)             :: lon
 468
 469    lat = 2.0*ATAN(EXP(proj%dlon*(proj%rsw + y-1.0)))*deg_per_rad - 90.0
 470    lon = (x-1.0)*proj%dlon*deg_per_rad + proj%lon1
 471    if (lon > 180.0) lon = lon - 360.0
 472    if (lon < -180.0) lon = lon + 360.0
 473
 474 end subroutine da_xyll_merc
 475
 476
 477 subroutine da_xyll_ps(x, y, proj, lat, lon)
 478
 479    ! This is the inverse subroutine da_of llij_ps.  It returns the
 480    ! latitude and longitude of an x/y point given the projection info
 481    ! structure.
 482
 483    implicit none
 484
 485    real, intent(in)                    :: x    ! Column
 486    real, intent(in)                    :: y    ! Row
 487    type (proj_info), intent(in)        :: proj
 488
 489    real, intent(out)                   :: lat     ! -90 -> 90 North
 490    real, intent(out)                   :: lon     ! -180 -> 180 East
 491
 492    real                                :: reflon
 493    real                                :: scale_top
 494    real                                :: xx,yy
 495    real                                :: gi2, r2
 496    real                                :: arccos
 497
 498    ! Compute the reference longitude by rotating 90 degrees to the east
 499    ! to find the longitude line parallel to the positive x-axis.
 500    reflon = proj%stdlon + 90.0
 501
 502    ! Compute numerator term of map scale factor
 503    scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
 504
 505    ! Compute radius to point of interest
 506    xx = x - proj%polei
 507    yy = (y - proj%polej) * proj%hemi
 508    r2 = xx**2 + yy**2
 509
 510    ! Now the magic code
 511    if (r2 == 0.0) then
 512       lat = proj%hemi * 90.0
 513       lon = reflon
 514    else
 515       gi2 = (proj%rebydx * scale_top)**2.0
 516       lat = deg_per_rad * proj%hemi * ASin((gi2-r2)/(gi2+r2))
 517       arccos = ACOS(xx/sqrt(r2))
 518       if (yy > 0) then
 519          lon = reflon + deg_per_rad * arccos
 520       else
 521          lon = reflon - deg_per_rad * arccos
 522       end if
 523    end if
 524
 525    ! Convert to a -180 -> 180 East convention
 526    if (lon > 180.0) lon = lon - 360.0
 527    if (lon < -180.0) lon = lon + 360.0
 528
 529 end subroutine da_xyll_ps
 530
 531
 532 subroutine da_set_lc(proj)
 533
 534    !-----------------------------------------------------------------------
 535    ! Purpose: Initialize the remaining items in the proj structure for a
 536    ! lambert conformal grid.
 537    !-----------------------------------------------------------------------
 538
 539    implicit none
 540
 541    type(proj_info), intent(inout) :: proj
 542
 543    real :: arg
 544    real :: deltalon1
 545    real :: tl1r
 546    real :: ctl1r
 547
 548    ! Compute cone factor
 549    call da_lc_cone(proj%truelat1, proj%truelat2, proj%cone)
 550    if (print_detail_map) then
 551       write(unit=stdout, fmt='(A,F8.6)') 'Computed cone factor: ', proj%cone
 552    end if
 553    ! Compute longitude differences and ensure we stay out of the
 554    ! forbidden "cut zone"
 555    deltalon1 = proj%lon1 - proj%stdlon
 556    if (deltalon1 .gt. +180.0) deltalon1 = deltalon1 - 360.0
 557    if (deltalon1 .lt. -180.0) deltalon1 = deltalon1 + 360.0
 558
 559    ! Convert truelat1 to radian and compute COS for later use
 560    tl1r = proj%truelat1 * rad_per_deg
 561    ctl1r = COS(tl1r)
 562
 563    ! Compute the radius to our known lower-left (SW) corner
 564    proj%rsw = proj%rebydx * ctl1r/proj%cone * &
 565            (TAN((90.0*proj%hemi-proj%lat1)*rad_per_deg/2.0) / &
 566             TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0))**proj%cone
 567
 568    ! Find pole point
 569    arg = proj%cone*(deltalon1*rad_per_deg)
 570    proj%polei = 1.0 - proj%hemi * proj%rsw * Sin(arg)
 571    proj%polej = 1.0 + proj%rsw * COS(arg)
 572    if (print_detail_map) then
 573       write(unit=stdout,fmt='(A,2F10.3)') 'Computed pole i/j = ', proj%polei, proj%polej
 574    end if
 575
 576 end subroutine da_set_lc
 577
 578
 579 subroutine da_set_ps(proj)
 580
 581    ! Initializes a polar-stereographic map projection from the partially
 582    ! filled proj structure. This routine computes the radius to the
 583    ! southwest corner and computes the i/j location of the pole for use
 584    ! in llij_ps and ijll_ps.
 585
 586    implicit none
 587
 588    type(proj_info), intent(inout)    :: proj
 589
 590    real :: ala1
 591    real :: alo1
 592    real :: reflon
 593    real :: scale_top
 594
 595    ! To define the cone factor for polar stereographic projection
 596    proj%cone = 1.0
 597
 598    reflon = proj%stdlon + 90.0
 599
 600    ! Compute numerator term of map scale factor
 601    scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
 602
 603    ! Compute radius to lower-left (SW) corner
 604    ala1 = proj%lat1 * rad_per_deg
 605    proj%rsw = proj%rebydx*COS(ala1)*scale_top/(1.0+proj%hemi*Sin(ala1))
 606
 607    ! Find the pole point
 608    alo1 = (proj%lon1 - reflon) * rad_per_deg
 609    proj%polei = proj%knowni - proj%rsw * COS(alo1)
 610    proj%polej = proj%knownj - proj%hemi * proj%rsw * Sin(alo1)
 611    if (print_detail_map) then
 612       write(unit=stdout,fmt='(A,2F10.1)') 'Computed (I,J) of pole point: ',proj%polei,proj%polej
 613    end if
 614
 615 end subroutine da_set_ps
 616
 617
 618 subroutine da_map_init(proj)
 619
 620    !-----------------------------------------------------------------------
 621    ! Purpose: Initializes the map projection structure to missing values
 622    !-----------------------------------------------------------------------
 623
 624    implicit none
 625
 626    type(proj_info), intent(inout)  :: proj
 627
 628    proj%lat1     = -999.9
 629    proj%lon1     = -999.9
 630    proj%dx       = -999.9
 631    proj%stdlon   = -999.9
 632    proj%truelat1 = -999.9
 633    proj%truelat2 = -999.9
 634    proj%hemi     = 0.0
 635    proj%cone     = -999.9
 636    proj%polei    = -999.9
 637    proj%polej    = -999.9
 638    proj%rsw      = -999.9
 639    proj%knowni   = -999.9
 640    proj%knownj   = -999.9
 641    proj%latinc   = -999.9
 642    proj%loninc   = -999.9
 643    proj%init     = .false.
 644
 645 end subroutine da_map_init
 646
 647
 648 subroutine da_map_set(proj_code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,latinc,loninc,proj)
 649    ! Given a partially filled proj_info structure, this routine computes
 650    ! polei, polej, rsw, and cone (if LC projection) to complete the
 651    ! structure.  This allows us to eliminate redundant calculations when
 652    ! calling the coordinate conversion routines multiple times for the
 653    ! same map.
 654    ! This will generally be the first routine called when a user wants
 655    ! to be able to use the coordinate conversion routines, and it
 656    ! will call the appropriate subroutines based on the
 657    ! proj%code which indicates which projection type  this is.
 658
 659    implicit none
 660
 661    integer,         intent(in)  :: proj_code
 662    real,            intent(in)  :: lat1
 663    real,            intent(in)  :: lon1
 664    real,            intent(in)  :: dx
 665    real,            intent(in)  :: stdlon
 666    real,            intent(in)  :: truelat1
 667    real,            intent(in)  :: truelat2
 668    real,            intent(in)  :: knowni , knownj
 669    real,            intent(in)  :: latinc, loninc
 670    type(proj_info), intent(out) :: proj
 671
 672    ! First, check for validity of mandatory variables in proj
 673    if (ABS(lat1) > 90.0) then
 674       write(stdout,*) "da_map_set.inc"// &
 675          "Latitude of origin corner required as follows: -90N <= lat1 < = 90N"
 676       stop
 677    end if
 678    if (ABS(lon1) > 180.0) then
 679       write(stdout,*) "da_map_set.inc"// &
 680            "Longitude of origin required as follows: -180E <= lon1 <= 180W"
 681       stop
 682    end if
 683    if ((dx .LE. 0.0).AND.(proj_code .NE. PROJ_LATLON)) then
 684       write(stdout,*)  "da_map_set.inc"// &
 685            "Require grid spacing (dx) in meters be positive!"
 686       stop
 687    end if
 688    if ((ABS(stdlon) > 180.0).AND.(proj_code .NE. PROJ_MERC)) then
 689       write(stdout,*) "da_map_set.inc"// &
 690            "Need orientation longitude (stdlon) as: -180E <= lon1 <= 180W"
 691       stop
 692    end if
 693    if (proj%code .NE. PROJ_LATLON .and. ABS(truelat1)>90.0) then
 694       write(stdout,*) "da_map_set.inc"// &
 695            "Set true latitude 1 for all projections!"
 696       stop
 697    end if
 698
 699    call da_map_init(proj)
 700    proj%code  = proj_code
 701    proj%lat1 = lat1
 702    proj%lon1 = lon1
 703    proj%knowni = knowni
 704    proj%knownj = knownj
 705    proj%dx    = dx
 706    proj%stdlon = stdlon
 707    proj%truelat1 = truelat1
 708    proj%truelat2 = truelat2
 709    proj%latinc   = latinc
 710    proj%loninc   = loninc
 711    if (proj%code .NE. PROJ_LATLON) then
 712       proj%dx = dx
 713       if (truelat1 < 0.0) then
 714          proj%hemi = -1.0
 715       else
 716          proj%hemi = 1.0
 717       end if
 718       proj%rebydx = earth_radius_m / dx
 719    end if
 720    pick_proj: select case(proj%code)
 721
 722       case(PROJ_PS)
 723          if (print_detail_map) then
 724             write(unit=stdout,fmt='(A)') 'Setting up POLAR STEREOGRAPHIC map...'
 725          end if
 726          call da_set_ps(proj)
 727
 728       case(PROJ_LC)
 729          if (print_detail_map) then
 730             write(unit=stdout,fmt='(A)') 'Setting up LAMBERT CONFORMAL map...'
 731          end if
 732          if (ABS(proj%truelat2) > 90.0) then
 733             if (print_detail_map) then
 734                write(unit=stdout,fmt='(A)') 'Second true latitude not set, assuming a tangent'
 735                write(unit=stdout,fmt='(A,F10.3)') 'projection at truelat1: ', proj%truelat1
 736             end if
 737             proj%truelat2=proj%truelat1
 738          end if
 739          call da_set_lc(proj)
 740
 741       case (PROJ_MERC)
 742          if (print_detail_map) then
 743             write(unit=stdout,fmt='(A)') 'Setting up MERCATOR map...'
 744          end if
 745          call da_set_merc(proj)
 746
 747       case (PROJ_LATLON)
 748          if (print_detail_map) then
 749             write(unit=stdout,fmt='(A)') 'Setting up CYLinDRICAL EQUIDISTANT LATLON map...'
 750          end if
 751          ! Convert lon1 to 0->360 notation
 752          if (proj%lon1 < 0.0) proj%lon1 = proj%lon1 + 360.0
 753
 754          proj%cone = 1.0
 755       case default
 756          write(unit=message(1),fmt='(A,I2)') 'Unknown projection code: ', proj%code
 757          write(stdout,*) "da_map_set.inc"//message(1:1)
 758          stop
 759
 760    end select pick_proj
 761    proj%init = .true.
 762
 763 end subroutine da_map_set
 764
 765
 766 subroutine da_set_merc(proj)
 767
 768    !--------------------------------------------------------------------------
 769    ! Purpose: Sets up the remaining basic elements for the mercator projection
 770    !--------------------------------------------------------------------------
 771
 772    implicit none
 773
 774    type(proj_info), intent(inout)       :: proj
 775
 776    real :: clain
 777
 778    !  Preliminary variables
 779
 780    clain = COS(rad_per_deg*proj%truelat1)
 781    proj%dlon = proj%dx / (earth_radius_m * clain)
 782
 783    ! Compute distance from equator to origin, and store in the
 784    ! proj%rsw tag.
 785
 786    proj%rsw = 0.0
 787    if (proj%lat1 .NE. 0.0) then
 788       proj%rsw = (alog(tan(0.5*((proj%lat1+90.)*rad_per_deg))))/proj%dlon
 789    end if
 790
 791 end subroutine da_set_merc
 792
 793
 794 subroutine da_lc_cone(truelat1, truelat2, cone)
 795
 796    !------------------------------------------------------------------------
 797    ! Purpose: compute the cone factor of a Lambert Conformal projection
 798    !------------------------------------------------------------------------
 799
 800    implicit none
 801
 802    real, intent(in)             :: truelat1  ! (-90 -> 90 degrees N)
 803    real, intent(in)             :: truelat2  !   "   "  "   "     "
 804    real, intent(out)            :: cone
 805
 806    ! First, see if this is a secant or tangent projection.  For tangent
 807    ! projections, truelat1 = truelat2 and the cone is tangent to the
 808    ! Earth's surface at this latitude.  For secant projections, the cone
 809    ! intersects the Earth's surface at each of the distinctly different
 810    ! latitudes
 811    if (abs(truelat1-truelat2) > 0.1) then
 812       cone = alog10(cos(truelat1*rad_per_deg)) - &
 813              alog10(cos(truelat2*rad_per_deg))
 814       cone = cone /(alog10(tan((45.0 - abs(truelat1)/2.0) * rad_per_deg)) - &
 815              alog10(tan((45.0 - abs(truelat2)/2.0) * rad_per_deg)))
 816    else
 817       cone = sin(abs(truelat1)*rad_per_deg)
 818    end if
 819
 820 end subroutine da_lc_cone
 821
 822 end module da_verif_tools
 823