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