da_xyll.inc

1 subroutine da_xyll(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    conv = 180.0 / pi
25 
26    ! seperate treatment for global fields
27    if (fg_format == fg_format_kma_global) then
28       xlat = yy * 180.0 /(coarse_jy-1)  -  90.0    
29       xlon = xx * 360.0 /(coarse_ix-1)  - 180.0    
30       return
31    end if
32 
33    CNTRI = real(coarse_ix+1)/2.
34    CNTRJ = real(coarse_jy+1)/2. 
35 
36    xcntr = 0.0
37 
38    ! CALCULATE X AND Y POSITIONS OF GRID
39 
40    x = (xx - 0.5)*dsm/coarse_ds + start_x
41    y = (yy - 0.5)*dsm/coarse_ds + start_y
42    x = xcntr + (x-cntri)*coarse_ds
43    y = ycntr + (y-cntrj)*coarse_ds
44 
45    ! NOW CALCULATE LAT AND LON OF THIS POinT
46 
47    if (map_projection.NE.3) then
48       if(Y.EQ.0.) then      
49          if (X.GE.0.0) FLP =  90.0/CONV 
50          if (X.LT.0.0) FLP = -90.0/CONV
51       else
52          if (PHIC.LT.0.0)then
53             FLP = ATAN2(X,Y)   
54          else
55             FLP = ATAN2(X,-Y) 
56          end if
57       end if 
58       FLPP = (FLP/cone_factor)*CONV+XLONC  
59       if (FLPP.LT.-180.) FLPP = FLPP + 360    
60       if (FLPP.GT.180.)  FLPP = FLPP - 360.  
61       XLON = FLPP 
62       ! NOW SOLVE FOR LATITUDE
63       R = sqrt(X*X+Y*Y)  
64       if (PHIC.LT.0.0) R = -R  
65       if (map_projection.EQ.1) then   
66          CELL = (R*cone_factor)/(earth_radius*Sin(PSI1))    
67          Rcone_factor  = 1.0/cone_factor   
68          CEL1 = TAN(PSI1/2.)*(CELL)**Rcone_factor    
69       end if 
70       if (map_projection.EQ.2) then
71          CELL = R/earth_radius        
72          CEL1 = CELL/(1.0+COS(PSI1))  
73       end if 
74       CEL2 = ATAN(CEL1)    
75       PSX  = 2.*CEL2*CONV
76       XLAT = POLE-PSX 
77    end if   
78    ! CALCULATIONS FOR MERCATOR LAT,LON    
79    if (map_projection.EQ.3) then   
80       XLON = XLONC + ((X-XCNTR)/C2)*CONV 
81       if (XLON.LT.-180.) XLON = XLON + 360
82       if (XLON.GT.180.)  XLON = XLON - 360.
83       CELL = EXP(Y/C2)  
84       XLAT = 2.*(CONV*ATAN(CELL))-90.0 
85    end if  
86 
87 end subroutine da_xyll
88 
89