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