NCEPLIBS-w3emc 2.12.0
Loading...
Searching...
No Matches
w3ft203.f
Go to the documentation of this file.
1C> @file
2C> @brief Convert (361,91) grid to (45,39) n. hemi. grid
3C> @author Ralph Jones @date 1994-05-18
4
5C> Convert a northern hemisphere 1.0 degree lat.,lon. 361 by
6C> 91 grid to a polar stereographic 45 by 39 grid. The polar
7C> stereographic map projection is true at 60 deg. n. , The mesh
8C> length is 190.5 km. and the oriention is 150 deg. w.
9C>
10C> ### Program History Log:
11C> Date | Programmer | Comment
12C> -----|------------|--------
13C> 1994-05-18 | Ralph Jones | Initial.
14C>
15C> @param[in] ALOLA 361*91 grid 1.0 lat,lon grid n. hemisphere
16C> 32851 point grid is o.n. 84 type ?? or ?? hex
17C> @param[in] INTERP 1 linear interpolation , ne.1 biquadratic
18C> @param[out] APOLA 45*39 grid of northern hemisphere. 1755 point grid is
19C> awips grid type 203
20C>
21C> @note
22C> - 1. W1 and w2 are used to store sets of constants which are
23C> reusable for repeated calls to the subroutine.
24C> - 2. Wind components are not rotated to the 45*39 grid orientation
25C> after interpolation. You may use w3fc08() to do this.
26C>
27C> @author Ralph Jones @date 1994-05-18
28 SUBROUTINE w3ft203(ALOLA,APOLA,INTERP)
29C
30 parameter(npts=1755,ii=45,jj=39)
31 parameter(orient=150.0,ipole=27,jpole=37)
32 parameter(xmesh=190.5)
33C
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)
41C
42 INTEGER IV(NPTS), JV(NPTS), JY(NPTS,4)
43 INTEGER IM1(NPTS), IP1(NPTS), IP2(NPTS)
44C
45 LOGICAL LIN
46C
47 SAVE
48C
49 equivalence(xi(1,1),xii(1)),(xj(1,1),xjj(1))
50C
51 DATA degprd/57.2957795/
52 DATA earthr/6371.2/
53 DATA intrpo/99/
54 DATA iswt /0/
55C
56 lin = .false.
57 IF (interp.EQ.1) lin = .true.
58C
59 IF (iswt.EQ.1) GO TO 900
60C
61 deg = 1.0
62 gi2 = (1.86603 * earthr) / xmesh
63 gi2 = gi2 * gi2
64C
65C NEXT 32 LINES OF CODE PUTS SUBROUTINE W3FB01 IN LINE
66C
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
73C
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
79C
80 xii(1647) = 1.0
81 DO 300 kk = 1,npts
82 angle(kk) = degprd * atan2(xjj(kk),xii(kk))
83 300 CONTINUE
84C
85 DO 400 kk = 1,npts
86 IF (angle(kk).LT.0.0) angle(kk) = angle(kk) + 360.0
87 400 CONTINUE
88C
89 DO 500 kk = 1,npts
90 wlon(kk) = 270.0 + orient - angle(kk)
91 500 CONTINUE
92C
93 DO 600 kk = 1,npts
94 IF (wlon(kk).LT.0.0) wlon(kk) = wlon(kk) + 360.0
95 600 CONTINUE
96C
97 DO 700 kk = 1,npts
98 IF (wlon(kk).GE.360.0) wlon(kk) = wlon(kk) - 360.0
99 700 CONTINUE
100C
101 xlat(1647) = 90.0
102 wlon(1647) = 0.0
103C
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
108C
109 iswt = 1
110 intrpo = interp
111 GO TO 1000
112C
113C AFTER THE 1ST CALL TO W3FT203 TEST INTERP, IF IT HAS
114C CHANGED RECOMPUTE SOME CONSTANTS
115C
116 900 CONTINUE
117 IF (interp.EQ.intrpo) GO TO 2100
118 intrpo = interp
119C
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
130C
131 IF (lin) GO TO 1400
132C
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
141C
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
151C
152 1400 CONTINUE
153C
154 IF (lin) GO TO 1700
155C
156 DO 1500 kk = 1,npts
157 IF (jv(kk).LT.2.OR.jv(kk).GT.89) xj2tm(kk) = 0.0
158 1500 CONTINUE
159C
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
166C
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
174C
175C LINEAR INTERPOLATION
176C
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
183C
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
192C
193 2100 CONTINUE
194 IF (lin) THEN
195C
196C LINEAR INTERPOLATION
197C
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
204C
205 DO 2300 kk = 1,npts
206 apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
207 & * xdelj(kk)
208 2300 CONTINUE
209C
210 ELSE
211C
212C QUADRATIC INTERPOLATION
213C
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
236C
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
242C
243C SET POLE POINT , WMO STANDARD FOR U OR V
244C
245 apola(1647) = alola(181,91)
246C
247 ENDIF
248C
249 RETURN
250 END
w3ft203
subroutine w3ft203(alola, apola, interp)
Convert a northern hemisphere 1.0 degree lat.,lon.
Definition w3ft203.f:29