var/da/da_tools/da_llxy_default_new.inc

   1 subroutine da_llxy_default_new (info)
   2
   3    !----------------------------------------------------------------------------
   4    !
   5    !                 routine llxy
   6    !                **************
   7    !
   8    !
   9    ! Purpose:  calculates the (x,y) location (dot) in the mesoscale grids
  10    ! -------   from latitudes and longitudes
  11    !
  12    !           for global domain co-ordinates
  13    !
  14    !  input:
  15    !  -----
  16    !   xlat:    latitudes
  17    !   xlon:    longitudes
  18    !
  19    ! output:
  20    ! -----
  21    !   x:        the coordinate in x (i)-direction.
  22    !   y:        the coordinate in y (j)-direction.
  23    !
  24    !----------------------------------------------------------------------------
  25
  26    implicit none
  27
  28    type(infa_type), intent(inout) :: info
  29
  30    real              :: dxlon
  31    real              :: xlat, xlon
  32    real              :: xx, yy, xc, yc
  33    real              :: cell, psi0, psx, r, flp
  34    real              :: centri, centrj
  35    real              :: ratio
  36    real              :: bb
  37    real, parameter   :: conv = 180.0 / pi
  38
  39    integer :: n
  40
  41    if (trace_use) call da_trace_entry("da_llxy_default_new")
  42
  43    ! Slow, but I can't be arsed to do all the temporary arrays
  44
  45    do n=1,ubound(info%lat,2)
  46       xlon = info%lon(1,n)
  47       xlat = info%lat(1,n)
  48
  49       xlat = max (xlat, -89.95)
  50       xlat = min (xlat, +89.95)
  51
  52       dxlon = xlon - xlonc
  53       if (dxlon >  180) dxlon = dxlon - 360.0
  54       if (dxlon < -180) dxlon = dxlon + 360.0
  55
  56       if (map_projection == 3) then
  57          xc = 0.0
  58          yc = YCNTR
  59
  60          cell = cos(xlat/conv)/(1.0+sin(xlat/conv))
  61          yy = -c2*alog(cell)
  62          xx = c2*dxlon/conv
  63       else
  64          psi0 = (pole - phic)/conv
  65          xc = 0.0
  66
  67          ! calculate x,y coords. relative to pole
  68
  69          flp = cone_factor*dxlon/conv
  70
  71          psx = (pole - xlat)/conv
  72
  73          if (map_projection == 2) then
  74             ! Polar stereographics:
  75             bb = 2.0*(cos(psi1/2.0)**2)
  76             yc = -earth_radius*bb*tan(psi0/2.0)
  77              r = -earth_radius*bb*tan(psx/2.0)
  78          else
  79             ! Lambert conformal:
  80             bb = -earth_radius/cone_factor*sin(psi1)
  81             yc = bb*(tan(psi0/2.0)/tan(psi1/2.0))**cone_factor
  82              r = bb*(tan(psx /2.0)/tan(psi1/2.0))**cone_factor
  83          end if
  84
  85          if (phic < 0.0) then
  86             xx = r*sin(flp)
  87             yy = r*cos(flp)
  88          else
  89             xx = -r*sin(flp)
  90             yy =  r*cos(flp)
  91          end if
  92
  93       end if
  94
  95       ! transform (1,1) to the origin
  96       ! the location of the center in the coarse domain
  97
  98       centri = real (coarse_ix + 1)/2.0
  99       centrj = real (coarse_jy + 1)/2.0
 100
 101       ! the (x,y) coordinates in the coarse domain
 102
 103       info%x(1,n) = (xx - xc)/coarse_ds + centri
 104       info%y(1,n) = (yy - yc)/coarse_ds + centrj
 105
 106       ratio = coarse_ds / dsm
 107
 108       ! only add 0.5 so that x/y is relative to first cross points:
 109
 110       info%x(:,n) = (info%x(1,n) - start_x)*ratio + 0.5
 111       info%y(:,n) = (info%y(1,n) - start_y)*ratio + 0.5
 112    end do
 113
 114    if (trace_use) call da_trace_exit("da_llxy_default_new")
 115
 116 end subroutine da_llxy_default_new
 117
 118