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