da_llxy_ps_new.inc

1 subroutine da_llxy_ps_new(proj,info)
2 
3    !-----------------------------------------------------------------------
4    ! Purpose: Given latitude (-90 to 90), longitude (-180 to 180), and the
5    ! standard polar-stereographic projection information via the 
6    ! public proj structure, this routine returns the x/y indices which
7    ! if within the domain range from 1->nx and 1->ny, respectively.
8    !-----------------------------------------------------------------------
9 
10    implicit none
11 
12    type(proj_info), intent(in)    :: proj
13    type(infa_type), intent(inout) :: info
14    
15    real    :: reflon
16    integer :: n
17    real    :: scale_top, ala, rm, alo
18 
19    if (trace_use) call da_trace_entry("da_llxy_ps_new")
20 
21    reflon = proj%stdlon + 90.0
22 
23    ! FAST
24 
25 !   x(:,:) = proj%polei + proj%rebydx * COS(lat(:,:) * rad_per_deg) &
26 !      * (1.0 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) &
27 !      / (1.0 + proj%hemi *SIN(lat(:,:) * rad_per_deg)) &
28 !      * COS((lon(:,:) - reflon) * rad_per_deg)
29 
30 !   y(:,:) = proj%polej + proj%hemi * proj%rebydx * COS(lat(:,:) * rad_per_deg) &
31 !      * (1.0 + proj%hemi * SIN(proj%truelat1 * rad_per_deg) &
32 !      / (1.0 + proj%hemi *SIN(lat(:,:) * rad_per_deg)) &
33 !      * SIN((lon(:,:) - reflon) * rad_per_deg)
34 
35 ! SLOW
36 
37    do n=lbound(info%lat,2),ubound(info%lat,2)
38       ! compute numerator term of map scale factor
39 
40       scale_top = 1.0 + proj%hemi * sin(proj%truelat1 * rad_per_deg)
41 
42       ! find radius to desired point
43       ala = info%lat(1,n) * rad_per_deg
44       rm = proj%rebydx * cos(ala) * scale_top/(1.0 + proj%hemi *sin(ala))
45       alo = (info%lon(1,n) - reflon) * rad_per_deg
46       info%x(:,n) = proj%polei + rm * cos(alo)
47       info%y(:,n) = proj%polej + proj%hemi * rm * sin(alo)
48    end do
49 
50    if (trace_use) call da_trace_exit("da_llxy_ps_new")
51  
52 end subroutine da_llxy_ps_new
53 
54