NCEPLIBS-w3emc 2.12.0
Loading...
Searching...
No Matches
w3ft43v.f
Go to the documentation of this file.
1C> @file
2C> @brief Convert (361,181) grid to (65,65) n. hemi. grid.
3C> @author Ralph Jones @date 1993-03-29
4
5C> Convert a global 1.0 degree lat.,lon. 361 by
6C> 181 grid to a polar stereographic 65 by 65 grid. the polar
7C> stereographic map projection is true at 60 deg. n. , the mesh
8C> length is 381 km. and the oriention is 80 deg. w.
9C>
10C> ### Program History Log:
11C> Date | Programmer | Comment
12C> -----|------------|--------
13C> 1993-03-29 | Ralph Jones | Add save statement.
14C>
15C> @param[in] ALOLA 361*181 grid 1.0 deg. lat,lon grid n. hemi.
16C> 65341 point grid. 360 * 181 one degree grib grid 3 was flipped, greenwish
17C> added to right side to make 361 * 181.
18C> @param[in] INTERP 1 linear interpolation , ne.1 biquadratic
19C> @param[out] APOLA 65*65 grid of northern hemisphere. 4225 point grid is
20C> o.n.84 type 27 or 1b hex
21C>
22C> @note
23C> - 1. W1 and w2 are used to store sets of constants which are
24C> reusable for repeated calls to the subroutine. 20 other arrays
25C> are saved and reused on the next calls to the subroutine.
26C> - 2. Wind components are not rotated to the 65*65 grid orientation
27C> after interpolation. You may use w3fc08 to do this.
28C> - 3. The about 1100 points below the equator will be in this map.
29C>
30C> @author Ralph Jones @date 1993-03-29
31 SUBROUTINE w3ft43v(ALOLA,APOLA,INTERP)
32C
33 parameter(npts=4225,ii=65,jj=65)
34 parameter(orient=80.0,ipole=33,jpole=33)
35 parameter(xmesh=381.0)
36C
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)
44C
45 INTEGER IV(NPTS), JV(NPTS), JY(NPTS,4)
46 INTEGER IM1(NPTS), IP1(NPTS), IP2(NPTS)
47C
48 LOGICAL LIN
49C
50 SAVE
51C
52 equivalence(xi(1,1),xii(1)),(xj(1,1),xjj(1))
53C
54 DATA degprd/57.2957795/
55 DATA earthr/6371.2/
56 DATA intrpo/99/
57 DATA iswt /0/
58C
59 lin = .false.
60 IF (interp.EQ.1) lin = .true.
61C
62 IF (iswt.EQ.1) GO TO 900
63C
64 deg = 1.0
65 gi2 = (1.86603 * earthr) / xmesh
66 gi2 = gi2 * gi2
67C
68C NEXT 32 LINES OF CODE PUTS SUBROUTINE W3FB05 IN LINE
69C
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
76C
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
82C
83 xii(2113) = 1.0
84 DO 300 kk = 1,npts
85 angle(kk) = degprd * atan2(xjj(kk),xii(kk))
86 300 CONTINUE
87C
88 DO 400 kk = 1,npts
89 IF (angle(kk).LT.0.0) angle(kk) = angle(kk) + 360.0
90 400 CONTINUE
91C
92 DO 500 kk = 1,npts
93 wlon(kk) = 270.0 + orient - angle(kk)
94 500 CONTINUE
95C
96 DO 600 kk = 1,npts
97 IF (wlon(kk).LT.0.0) wlon(kk) = wlon(kk) + 360.0
98 600 CONTINUE
99C
100 DO 700 kk = 1,npts
101 IF (wlon(kk).GE.360.0) wlon(kk) = wlon(kk) - 360.0
102 700 CONTINUE
103C
104 xlat(2113) = 90.0
105 wlon(2113) = 0.0
106C
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
111C
112 iswt = 1
113 intrpo = interp
114 GO TO 1000
115C
116C AFTER THE 1ST CALL TO W3FT43V TEST INTERP, IF IT HAS
117C CHANGED RECOMPUTE SOME CONSTANTS
118C
119 900 CONTINUE
120 IF (interp.EQ.intrpo) GO TO 2100
121 intrpo = interp
122C
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
133C
134 IF (lin) GO TO 1400
135C
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
144C
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
154C
155 1400 CONTINUE
156C
157 IF (lin) GO TO 1700
158C
159 DO 1500 kk = 1,npts
160 IF (jv(kk).GE.180) xj2tm(kk) = 0.0
161 1500 CONTINUE
162C
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
169C
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
177C
178C LINEAR INTERPOLATION
179C
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
184C
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
191C
192 2100 CONTINUE
193 IF (lin) THEN
194C
195C LINEAR INTERPOLATION
196C
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
203C
204 DO 2300 kk = 1,npts
205 apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
206 & * xdelj(kk)
207 2300 CONTINUE
208C
209 ELSE
210C
211C QUADRATIC INTERPOLATION
212C
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
235C
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
241C
242 ENDIF
243C
244C SET POLE POINT , WMO STANDARD FOR U OR V
245C
246 apola(2113) = alola(181,181)
247C
248 RETURN
249 END
w3ft43v
subroutine w3ft43v(alola, apola, interp)
Convert a global 1.0 degree lat.,lon.
Definition w3ft43v.f:32