NCEPLIBS-w3emc  2.11.0
w3ft213.f
Go to the documentation of this file.
1 C> @file
2 C> @brief Convert (361,91) grid to (129,85) n. hemi. grid
3 C> @author Ralph Jones @date 1993-10-23
4 
5 C> Convert a northern hemisphere 1.0 degree lat.,lon. 361 by
6 C> 91 grid to a polar stereographic 129 by 85 grid. The polar
7 C> stereographic map projection is true at 60 deg. n. , The mesh
8 C> length is 95.25 km. and the oriention is 105 deg. w.
9 C> awips grid 213 national - conus - double resolution
10 C>
11 C> ### Program History Log:
12 C> Date | Programmer | Comment
13 C> -----|------------|--------
14 C> 1993-10-23 | Ralph Jones | Initial.
15 C>
16 C> @param[in] ALOLA 361*91 grid 1.0 deg. lat,lon grid n. hemi.
17 C> 32851 point grid. 360 * 181 one degree grib grid 3 was flipped, greenwish added
18 C> to right side and cut to 361 * 91.
19 C> @param[in] INTERP 1 linear interpolation , ne.1 biquadratic
20 C> @param[out] APOLA 129*85 grid of northern hemisphere. 10965 point grid is
21 C> awips grid type 213
22 C> @note
23 C> - 1. W1 and w2 are used to store sets of constants which are
24 C> reusable for repeated calls to the subroutine.
25 C> - 2. Wind components are not rotated to the 129*85 grid orientation
26 C> after interpolation. You may use w3fc08() to do this.
27 C>
28 C> @author Ralph Jones @date 1993-10-23
29  SUBROUTINE w3ft213(ALOLA,APOLA,INTERP)
30 C
31  parameter(npts=10965,ii=129,jj=85)
32  parameter(orient=105.0,ipole=65,jpole=89)
33  parameter(xmesh=95.250)
34 C
35  REAL R2(NPTS), WLON(NPTS)
36  REAL XLAT(NPTS), XI(II,JJ), XJ(II,JJ)
37  REAL XII(NPTS), XJJ(NPTS), ANGLE(NPTS)
38  REAL ALOLA(361,91), APOLA(NPTS), ERAS(NPTS,4)
39  REAL W1(NPTS), W2(NPTS)
40  REAL XDELI(NPTS), XDELJ(NPTS)
41  REAL XI2TM(NPTS), XJ2TM(NPTS)
42 C
43  INTEGER IV(NPTS), JV(NPTS), JY(NPTS,4)
44  INTEGER IM1(NPTS), IP1(NPTS), IP2(NPTS)
45 C
46  LOGICAL LIN
47 C
48  SAVE
49 C
50  equivalence(xi(1,1),xii(1)),(xj(1,1),xjj(1))
51 C
52  DATA degprd/57.2957795/
53  DATA earthr/6371.2/
54  DATA intrpo/99/
55  DATA iswt /0/
56 C
57  lin = .false.
58  IF (interp.EQ.1) lin = .true.
59 C
60  IF (iswt.EQ.1) GO TO 900
61 C
62  deg = 1.0
63  gi2 = (1.86603 * earthr) / xmesh
64  gi2 = gi2 * gi2
65 C
66 C NEXT 32 LINES OF CODE PUTS SUBROUTINE W3FB05 IN LINE
67 C
68  DO 100 j = 1,jj
69  xj1 = j - jpole
70  DO 100 i = 1,ii
71  xi(i,j) = i - ipole
72  xj(i,j) = xj1
73  100 CONTINUE
74 C
75  DO 200 kk = 1,npts
76  r2(kk) = xjj(kk) * xjj(kk) + xii(kk) * xii(kk)
77  xlat(kk) = degprd *
78  & asin((gi2 - r2(kk)) / (gi2 + r2(kk)))
79  200 CONTINUE
80 C
81  DO 300 kk = 1,npts
82  angle(kk) = degprd * atan2(xjj(kk),xii(kk))
83  300 CONTINUE
84 C
85  DO 400 kk = 1,npts
86  IF (angle(kk).LT.0.0) angle(kk) = angle(kk) + 360.0
87  400 CONTINUE
88 C
89  DO 500 kk = 1,npts
90  wlon(kk) = 270.0 + orient - angle(kk)
91  500 CONTINUE
92 C
93  DO 600 kk = 1,npts
94  IF (wlon(kk).LT.0.0) wlon(kk) = wlon(kk) + 360.0
95  600 CONTINUE
96 C
97  DO 700 kk = 1,npts
98  IF (wlon(kk).GE.360.0) wlon(kk) = wlon(kk) - 360.0
99  700 CONTINUE
100 C
101  DO 800 kk = 1,npts
102  w1(kk) = (360.0 - wlon(kk)) / deg + 1.0
103  w2(kk) = xlat(kk) / deg + 1.0
104  800 CONTINUE
105 C
106  iswt = 1
107  intrpo = interp
108  GO TO 1000
109 C
110 C AFTER THE 1ST CALL TO W3FT213 TEST INTERP, IF IT HAS
111 C CHANGED RECOMPUTE SOME CONSTANTS
112 C
113  900 CONTINUE
114  IF (interp.EQ.intrpo) GO TO 2100
115  intrpo = interp
116 C
117  1000 CONTINUE
118  DO 1100 k = 1,npts
119  iv(k) = w1(k)
120  jv(k) = w2(k)
121  xdeli(k) = w1(k) - iv(k)
122  xdelj(k) = w2(k) - jv(k)
123  ip1(k) = iv(k) + 1
124  jy(k,3) = jv(k) + 1
125  jy(k,2) = jv(k)
126  1100 CONTINUE
127 C
128  IF (lin) GO TO 1400
129 C
130  DO 1200 k = 1,npts
131  ip2(k) = iv(k) + 2
132  im1(k) = iv(k) - 1
133  jy(k,1) = jv(k) - 1
134  jy(k,4) = jv(k) + 2
135  xi2tm(k) = xdeli(k) * (xdeli(k) - 1.0) * .25
136  xj2tm(k) = xdelj(k) * (xdelj(k) - 1.0) * .25
137  1200 CONTINUE
138 C
139  DO 1300 kk = 1,npts
140  IF (iv(kk).EQ.1) THEN
141  ip2(kk) = 3
142  im1(kk) = 360
143  ELSE IF (iv(kk).EQ.360) THEN
144  ip2(kk) = 2
145  im1(kk) = 359
146  ENDIF
147  1300 CONTINUE
148 C
149  1400 CONTINUE
150 C
151  IF (lin) GO TO 1700
152 C
153  DO 1500 kk = 1,npts
154  IF (jv(kk).LT.2.OR.jv(kk).GT.89) xj2tm(kk) = 0.0
155  1500 CONTINUE
156 C
157  DO 1600 kk = 1,npts
158  IF (ip2(kk).LT.1) ip2(kk) = 1
159  IF (im1(kk).LT.1) im1(kk) = 1
160  IF (ip2(kk).GT.361) ip2(kk) = 361
161  IF (im1(kk).GT.361) im1(kk) = 361
162  1600 CONTINUE
163 C
164  1700 CONTINUE
165  DO 1800 kk = 1,npts
166  IF (iv(kk).LT.1) iv(kk) = 1
167  IF (ip1(kk).LT.1) ip1(kk) = 1
168  IF (iv(kk).GT.361) iv(kk) = 361
169  IF (ip1(kk).GT.361) ip1(kk) = 361
170  1800 CONTINUE
171 C
172 C LINEAR INTERPOLATION
173 C
174  DO 1900 kk = 1,npts
175  IF (jy(kk,2).LT.1) jy(kk,2) = 1
176  IF (jy(kk,2).GT.91) jy(kk,2) = 91
177  IF (jy(kk,3).LT.1) jy(kk,3) = 1
178  IF (jy(kk,3).GT.91) jy(kk,3) = 91
179  1900 CONTINUE
180 C
181  IF (.NOT.lin) THEN
182  DO 2000 kk = 1,npts
183  IF (jy(kk,1).LT.1) jy(kk,1) = 1
184  IF (jy(kk,1).GT.91) jy(kk,1) = 91
185  IF (jy(kk,4).LT.1) jy(kk,4) = 1
186  IF (jy(kk,4).GT.91) jy(kk,4) = 91
187  2000 CONTINUE
188  ENDIF
189 C
190  2100 CONTINUE
191  IF (lin) THEN
192 C
193 C LINEAR INTERPOLATION
194 C
195  DO 2200 kk = 1,npts
196  eras(kk,2) = (alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
197  & * xdeli(kk) + alola(iv(kk),jy(kk,2))
198  eras(kk,3) = (alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
199  & * xdeli(kk) + alola(iv(kk),jy(kk,3))
200  2200 CONTINUE
201 C
202  DO 2300 kk = 1,npts
203  apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
204  & * xdelj(kk)
205  2300 CONTINUE
206 C
207  ELSE
208 C
209 C QUADRATIC INTERPOLATION
210 C
211  DO 2400 kk = 1,npts
212  eras(kk,1)=(alola(ip1(kk),jy(kk,1))-alola(iv(kk),jy(kk,1)))
213  & * xdeli(kk) + alola(iv(kk),jy(kk,1)) +
214  & ( alola(im1(kk),jy(kk,1)) - alola(iv(kk),jy(kk,1))
215  & - alola(ip1(kk),jy(kk,1))+alola(ip2(kk),jy(kk,1)))
216  & * xi2tm(kk)
217  eras(kk,2)=(alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
218  & * xdeli(kk) + alola(iv(kk),jy(kk,2)) +
219  & ( alola(im1(kk),jy(kk,2)) - alola(iv(kk),jy(kk,2))
220  & - alola(ip1(kk),jy(kk,2))+alola(ip2(kk),jy(kk,2)))
221  & * xi2tm(kk)
222  eras(kk,3)=(alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
223  & * xdeli(kk) + alola(iv(kk),jy(kk,3)) +
224  & ( alola(im1(kk),jy(kk,3)) - alola(iv(kk),jy(kk,3))
225  & - alola(ip1(kk),jy(kk,3))+alola(ip2(kk),jy(kk,3)))
226  & * xi2tm(kk)
227  eras(kk,4)=(alola(ip1(kk),jy(kk,4))-alola(iv(kk),jy(kk,4)))
228  & * xdeli(kk) + alola(iv(kk),jy(kk,4)) +
229  & ( alola(im1(kk),jy(kk,4)) - alola(iv(kk),jy(kk,4))
230  & - alola(ip1(kk),jy(kk,4))+alola(ip2(kk),jy(kk,4)))
231  & * xi2tm(kk)
232  2400 CONTINUE
233 C
234  DO 2500 kk = 1,npts
235  apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
236  & * xdelj(kk) + (eras(kk,1) - eras(kk,2)
237  & - eras(kk,3) + eras(kk,4)) * xj2tm(kk)
238  2500 CONTINUE
239 C
240  ENDIF
241 C
242  RETURN
243  END
w3ft213
subroutine w3ft213(ALOLA, APOLA, INTERP)
Convert a northern hemisphere 1.0 degree lat.,lon.
Definition: w3ft213.f:30