var/da/da_tools/da_xyll_default.inc

   1 subroutine da_xyll_default(XX,YY,XLAT,XLON)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: calculates the latitudes and longitudes from the
   5    !                  (x,y) location (dot) in the mesoscale grids.
   6    ! on entry     :
   7    ! x            : the coordinate in x (j)-direction.
   8    ! y            : the coordinate in y (i)-direction.
   9    !
  10    ! on exit      :
  11    ! xlat         : latitudes
  12    ! xlon         : longitudes
  13    !-----------------------------------------------------------------------
  14
  15    implicit none
  16
  17    real, intent(in)  :: XX, YY
  18    real, intent(out) :: XLAT,XLON
  19
  20    real              :: flp, flpp, r, cell, cel1, cel2
  21    real              :: rcone_factor, psx, conv
  22    real              :: cntri, cntrj, x, y, xcntr
  23
  24    if (trace_use) call da_trace_entry("da_xyll_default")
  25
  26    conv = 180.0 / pi
  27
  28    ! seperate treatment for global fields
  29    if (fg_format == fg_format_kma_global) then
  30       xlat = yy * 180.0 /(coarse_jy-1)  -  90.0
  31       xlon = xx * 360.0 /(coarse_ix-1)  - 180.0
  32       return
  33    end if
  34
  35    cntri = real(coarse_ix+1)/2.0
  36    cntrj = real(coarse_jy+1)/2.0
  37
  38    xcntr = 0.0
  39
  40    ! calculate x and y positions of grid
  41
  42    x = (xx - 0.5)*dsm/coarse_ds + start_x
  43    y = (yy - 0.5)*dsm/coarse_ds + start_y
  44    x = xcntr + (x-cntri)*coarse_ds
  45    y = ycntr + (y-cntrj)*coarse_ds
  46
  47    ! now calculate lat and lon of this point
  48
  49    if (map_projection.ne.3) then
  50       if(y.eq.0.0) then
  51          if (x.ge.0.0) flp =  90.0/conv
  52          if (x.lt.0.0) flp = -90.0/conv
  53       else
  54          if (phic.lt.0.0)then
  55             flp = atan2(x,y)
  56          else
  57             flp = atan2(x,-y)
  58          end if
  59       end if
  60       flpp = (flp/cone_factor)*conv+xlonc
  61       if (flpp.lt.-180.0) flpp = flpp + 360.0
  62       if (flpp.gt.180.0)  flpp = flpp - 360.0
  63       xlon = flpp
  64       ! now solve for latitude
  65       r = sqrt(x*x+y*y)
  66       if (phic.lt.0.0) r = -r
  67       if (map_projection.eq.1) then
  68          cell = (r*cone_factor)/(earth_radius*sin(psi1))
  69          rcone_factor  = 1.0/cone_factor
  70          cel1 = tan(psi1/2.0)*(cell)**rcone_factor
  71       end if
  72       if (map_projection.eq.2) then
  73          cell = r/earth_radius
  74          cel1 = cell/(1.0+cos(psi1))
  75       end if
  76       cel2 = atan(cel1)
  77       psx  = 2.0*cel2*conv
  78       xlat = pole-psx
  79    end if
  80    ! calculations for mercator lat,lon
  81    if (map_projection.eq.3) then
  82       xlon = xlonc + ((x-xcntr)/c2)*conv
  83       if (xlon.lt.-180.0) xlon = xlon + 360.0
  84       if (xlon.gt.180.0)  xlon = xlon - 360.0
  85       cell = exp(y/c2)
  86       xlat = 2.0*(conv*atan(cell))-90.0
  87    end if
  88
  89    if (trace_use) call da_trace_exit("da_xyll_default")
  90
  91 end subroutine da_xyll_default
  92
  93