NCEPLIBS-ip 4.0.0
spectral_interp_mod.f90
Go to the documentation of this file.
1module spectral_interp_mod
2 use gdswzd_mod
3 use ip_grid_mod
4 use ip_grid_descriptor_mod
5 use ip_grid_factory_mod
6 implicit none
7
8 private
9 public :: interpolate_spectral
10
11 interface interpolate_spectral
12 module procedure interpolate_spectral_scalar
13 module procedure interpolate_spectral_vector
14 end interface interpolate_spectral
15
16 interface polates4
17 module procedure polates4_grib1
18 module procedure polates4_grib2
19 end interface polates4
20
21 interface polatev4
22 module procedure polatev4_grib1
23 module procedure polatev4_grib2
24 end interface polatev4
25
26contains
27
28 subroutine interpolate_spectral_scalar(IPOPT,grid_in,grid_out, &
29 MI,MO,KM,IBI,GI, &
30 NO,RLAT,RLON,IBO,LO,GO,IRET)
31 INTEGER, INTENT(IN ) :: IPOPT(20)
32 class(ip_grid), intent(in) :: grid_in, grid_out
33 INTEGER, INTENT(IN ) :: MI, MO
34 INTEGER, INTENT(IN ) :: IBI(KM), KM
35 INTEGER, INTENT( OUT) :: IBO(KM), IRET, NO
36 !
37 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
38 !
39 REAL, INTENT(IN ) :: GI(MI,KM)
40 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
41 REAL, INTENT( OUT) :: GO(MO,KM)
42
43
44 select type(desc_in => grid_in%descriptor)
45 type is(grib1_descriptor)
46 select type(desc_out => grid_out%descriptor)
47 type is(grib1_descriptor)
48 CALL polates4(ipopt,desc_in%gds,desc_out%gds,mi,mo,km,ibi,gi,no,rlat,rlon,ibo,lo,go,iret)
49 end select
50
51 type is(grib2_descriptor)
52 select type(desc_out => grid_out%descriptor)
53 type is(grib2_descriptor)
54 CALL polates4(ipopt,desc_in%gdt_num,desc_in%gdt_tmpl,desc_in%gdt_len, &
55 desc_out%gdt_num,desc_out%gdt_tmpl,desc_out%gdt_len, &
56 mi,mo,km,ibi,gi,no,rlat,rlon,ibo,lo,go,iret)
57 end select
58 end select
59
60
61 end subroutine interpolate_spectral_scalar
62
63
64 subroutine interpolate_spectral_vector(IPOPT,grid_in,grid_out, &
65 MI,MO,KM,IBI,UI,VI, &
66 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
67 class(ip_grid), intent(in) :: grid_in, grid_out
68 INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
69 INTEGER, INTENT(IN ) :: KM, MI, MO
70 INTEGER, INTENT( OUT) :: IRET, IBO(KM), NO
71 !
72 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
73 !
74 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
75 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
76 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
77 REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
78
79
80 select type(desc_in => grid_in%descriptor)
81 type is(grib1_descriptor)
82 select type(desc_out => grid_out%descriptor)
83 type is(grib1_descriptor)
84 CALL polatev4_grib1(ipopt,desc_in%gds,desc_out%gds,mi,mo,km,ibi,ui,vi,&
85 no,rlat,rlon,crot,srot,ibo,lo,uo,vo,iret)
86 end select
87
88 type is(grib2_descriptor)
89 select type(desc_out => grid_out%descriptor)
90 type is(grib2_descriptor)
91 CALL polatev4(ipopt,desc_in%gdt_num,desc_in%gdt_tmpl,desc_in%gdt_len, &
92 desc_out%gdt_num,desc_out%gdt_tmpl,desc_out%gdt_len, &
93 mi,mo,km,ibi,ui,vi,&
94 no,rlat,rlon,crot,srot,ibo,lo,uo,vo,iret)
95 end select
96 end select
97
98
99 end subroutine interpolate_spectral_vector
100
101
102 SUBROUTINE polates4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
103 IGDTNUMO,IGDTMPLO,IGDTLENO, &
104 MI,MO,KM,IBI,GI, &
105 NO,RLAT,RLON,IBO,LO,GO,IRET)
106 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
107 !
108 ! SUBPROGRAM: POLATES4 INTERPOLATE SCALAR FIELDS (SPECTRAL)
109 ! PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
110 !
111 ! ABSTRACT: THIS SUBPROGRAM PERFORMS SPECTRAL INTERPOLATION
112 ! FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
113 ! IT REQUIRES THAT THE INPUT FIELDS BE UNIFORMLY GLOBAL.
114 ! OPTIONS ALLOW CHOICES BETWEEN TRIANGULAR SHAPE (IPOPT(1)=0)
115 ! AND RHOMBOIDAL SHAPE (IPOPT(1)=1) WHICH HAS NO DEFAULT;
116 ! A SECOND OPTION IS THE TRUNCATION (IPOPT(2)) WHICH DEFAULTS
117 ! TO A SENSIBLE TRUNCATION FOR THE INPUT GRID (IF OPT(2)=-1).
118 ! NOTE THAT IF THE OUTPUT GRID IS NOT FOUND IN A SPECIAL LIST,
119 ! THEN THE TRANSFORM BACK TO GRID IS NOT VERY FAST.
120 ! THIS SPECIAL LIST CONTAINS GLOBAL CYLINDRICAL GRIDS,
121 ! POLAR STEREOGRAPHIC GRIDS CENTERED AT THE POLE
122 ! AND MERCATOR GRIDS. ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
123 !
124 ! THE CODE RECOGNIZES THE FOLLOWING PROJECTIONS, WHERE
125 ! "IGDTNUMI/O" IS THE GRIB 2 GRID DEFINTION TEMPLATE NUMBER
126 ! FOR THE INPUT AND OnUTPUT GRIDS, RESPECTIVELY:
127 ! (IGDTNUMI/O=00) EQUIDISTANT CYLINDRICAL
128 ! (IGDTNUMO =01) ROTATED EQUIDISTANT CYLINDRICAL. "E" AND
129 ! NON-"E" STAGGERED
130 ! (IGDTNUMO =10) MERCATOR CYLINDRICAL
131 ! (IGDTNUMO =20) POLAR STEREOGRAPHIC AZIMUTHAL
132 ! (IGDTNUMO =30) LAMBERT CONFORMAL CONICAL
133 ! (IGDTNUMI/O=40) GAUSSIAN CYLINDRICAL
134 ! AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
135 ! AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
136 ! ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
137 ! IF IGDTNUMO<0, IN WHICH CASE THE NUMBER OF POINTS
138 ! AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT.
139 ! OUTPUT BITMAPS WILL NOT BE CREATED.
140 !
141 ! PROGRAM HISTORY LOG:
142 ! 96-04-10 IREDELL
143 ! 2001-06-18 IREDELL IMPROVE DETECTION OF SPECIAL FAST TRANSFORM
144 ! 2015-01-27 GAYNO REPLACE CALLS TO GDSWIZ WITH NEW MERGED
145 ! VERSION OF GDSWZD.
146 ! 2015-07-13 GAYNO CONVERT TO GRIB 2. REPLACE GRIB 1 KGDS ARRAYS
147 ! WITH GRIB 2 GRID DEFINITION TEMPLATE ARRAYS.
148 !
149 ! USAGE: CALL POLATES4(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
150 ! IGDTNUMO,IGDTMPLO,IGDTLENO, &
151 ! MI,MO,KM,IBI,LI,GI, &
152 ! NO,RLAT,RLON,IBO,LO,GO,IRET)
153 !
154 ! INPUT ARGUMENT LIST:
155 ! IPOPT - INTEGER (20) INTERPOLATION OPTIONS
156 ! IPOPT(1)=0 FOR TRIANGULAR, IPOPT(1)=1 FOR RHOMBOIDAL;
157 ! IPOPT(2) IS TRUNCATION NUMBER
158 ! (DEFAULTS TO SENSIBLE IF IPOPT(2)=-1).
159 ! IGDTNUMI - INTEGER GRID DEFINITION TEMPLATE NUMBER - INPUT GRID.
160 ! CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
161 ! NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
162 ! 00 - EQUIDISTANT CYLINDRICAL
163 ! 01 - ROTATED EQUIDISTANT CYLINDRICAL. "E"
164 ! AND NON-"E" STAGGERED
165 ! 10 - MERCATOR CYCLINDRICAL
166 ! 20 - POLAR STEREOGRAPHIC AZIMUTHAL
167 ! 30 - LAMBERT CONFORMAL CONICAL
168 ! 40 - GAUSSIAN EQUIDISTANT CYCLINDRICAL
169 ! IGDTMPLI - INTEGER (IGDTLENI) GRID DEFINITION TEMPLATE ARRAY -
170 ! INPUT GRID. CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
171 ! OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE:
172 ! (SECTION 3 INFO). SEE COMMENTS IN ROUTINE
173 ! IPOLATES FOR COMPLETE DEFINITION.
174 ! IGDTLENI - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
175 ! TEMPLATE ARRAY - INPUT GRID. CORRESPONDS TO THE GFLD%IGDTLEN
176 ! COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
177 ! IGDTNUMO - INTEGER GRID DEFINITION TEMPLATE NUMBER - OUTPUT GRID.
178 ! CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
179 ! NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE. IGDTNUMO<0
180 ! MEANS INTERPOLATE TO RANDOM STATION POINTS.
181 ! OTHERWISE, SAME DEFINITION AS "IGDTNUMI".
182 ! IGDTMPLO - INTEGER (IGDTLENO) GRID DEFINITION TEMPLATE ARRAY -
183 ! OUTPUT GRID. CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
184 ! OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
185 ! (SECTION 3 INFO). SEE COMMENTS IN ROUTINE
186 ! IPOLATES FOR COMPLETE DEFINITION.
187 ! IGDTLENO - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
188 ! TEMPLATE ARRAY - OUTPUT GRID. CORRESPONDS TO THE GFLD%IGDTLEN
189 ! COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
190 ! MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
191 ! OR DIMENSION OF INPUT GRID FIELDS IF KM=1
192 ! MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
193 ! OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
194 ! KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE
195 ! IBI - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
196 ! GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
197 ! RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF IGDTNUMO<0)
198 ! RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF IGDTNUMO<0)
199 !
200 ! OUTPUT ARGUMENT LIST:
201 ! NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF IGDTNUMO>=0)
202 ! RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF IGDTNUMO>=0)
203 ! RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF IGDTNUMO>=0)
204 ! IBO - INTEGER (KM) OUTPUT BITMAP FLAGS
205 ! LO - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
206 ! GO - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
207 ! IRET - INTEGER RETURN CODE
208 ! 0 SUCCESSFUL INTERPOLATION
209 ! 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
210 ! 3 UNRECOGNIZED OUTPUT GRID
211 ! 41 INVALID NONGLOBAL INPUT GRID
212 ! 42 INVALID SPECTRAL METHOD PARAMETERS
213 !
214 ! SUBPROGRAMS CALLED:
215 ! EARTH_RADIUS DETERMINE SIZE/SHAPE OF EARTH
216 ! GDSWZD GRID DESCRIPTION SECTION WIZARD
217 ! SPTRUN SPECTRALLY TRUNCATE GRIDDED SCALAR FIELDS
218 ! SPTRUNS SPECTRALLY INTERPOLATE SCALARS TO POLAR STEREO.
219 ! SPTRUNM SPECTRALLY INTERPOLATE SCALARS TO MERCATOR
220 ! SPTRUNG SPECTRALLY INTERPOLATE SCALARS TO STATIONS
221 !
222 ! ATTRIBUTES:
223 ! LANGUAGE: FORTRAN 90
224 !
225 !$$$
226 INTEGER, INTENT(IN ) :: IGDTNUMI, IGDTLENI
227 INTEGER, INTENT(IN ) :: IGDTMPLI(IGDTLENI)
228 INTEGER, INTENT(IN ) :: IGDTNUMO, IGDTLENO
229 INTEGER, INTENT(IN ) :: IGDTMPLO(IGDTLENO)
230 INTEGER, INTENT(IN ) :: IPOPT(20)
231 INTEGER, INTENT(IN ) :: MI, MO
232 INTEGER, INTENT(IN ) :: IBI(KM), KM
233 INTEGER, INTENT( OUT) :: IBO(KM), IRET, NO
234 !
235 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
236 !
237 REAL, INTENT(IN ) :: GI(MI,KM)
238 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
239 REAL, INTENT( OUT) :: GO(MO,KM)
240 !
241 REAL, PARAMETER :: FILL=-9999.
242 REAL, PARAMETER :: PI=3.14159265358979
243 REAL, PARAMETER :: DPR=180./pi
244 !
245 INTEGER :: IDRTI, IDRTO, IG, JG, IM, JM
246 INTEGER :: IGO, JGO, IMO, JMO
247 INTEGER :: ISCAN, JSCAN, NSCAN
248 INTEGER :: ISCANO, JSCANO, NSCANO
249 INTEGER :: ISKIPI, JSKIPI, ISCALE
250 INTEGER :: IMAXI, JMAXI, ISPEC
251 INTEGER :: IP, IPRIME, IPROJ, IROMB, K
252 INTEGER :: MAXWV, N, NI, NJ, NPS
253 !
254 REAL :: DE, DR, DY
255 REAL :: DLAT, DLON, DLATO, DLONO
256 REAL :: GO2(MO,KM), H, HI, HJ
257 REAL :: ORIENT, SLAT, RERTH, E2
258 REAL :: RLAT1, RLON1, RLAT2, RLON2, RLATI
259 REAL :: XMESH, XP, YP
260 REAL :: XPTS(MO), YPTS(MO)
261
262 type(grib2_descriptor) :: desc_in, desc_out
263 class(ip_grid), allocatable :: grid_in, grid_out
264
265 desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
266 desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
267
268 call init_grid(grid_in, desc_in)
269 call init_grid(grid_out, desc_out)
270 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
272 iret=0
273 IF(igdtnumo.GE.0) THEN
274 !CALL GDSWZD(IGDTNUMO,IGDTMPLO,IGDTLENO, 0,MO,FILL,XPTS,YPTS,RLON,RLAT,NO)
275 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
276 IF(no.EQ.0) iret=3
277 ENDIF
278 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
279 ! AFFIRM APPROPRIATE INPUT GRID
280 ! LAT/LON OR GAUSSIAN
281 ! NO BITMAPS
282 ! FULL ZONAL COVERAGE
283 ! FULL MERIDIONAL COVERAGE
284 idrti=igdtnumi
285 IF(idrti==40) idrti=4
286 IF(idrti==0.OR.idrti==4)THEN
287 im=igdtmpli(8)
288 jm=igdtmpli(9)
289 iscale=igdtmpli(10)*igdtmpli(11)
290 IF(iscale==0) iscale=10**6
291 rlon1=float(igdtmpli(13))/float(iscale)
292 rlon2=float(igdtmpli(16))/float(iscale)
293 iscan=mod(igdtmpli(19)/128,2)
294 jscan=mod(igdtmpli(19)/64,2)
295 nscan=mod(igdtmpli(19)/32,2)
296 ELSE
297 iret=41
298 ENDIF
299 DO k=1,km
300 IF(ibi(k).NE.0) iret=41
301 ENDDO
302 IF(iret.EQ.0) THEN
303 IF(iscan.EQ.0) THEN
304 dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
305 ELSE
306 dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
307 ENDIF
308 ig=nint(360/abs(dlon))
309 iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
310 imaxi=ig
311 jmaxi=jm
312 IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
313 ENDIF
314 IF(iret.EQ.0.AND.idrti.EQ.0) THEN
315 iscale=igdtmpli(10)*igdtmpli(11)
316 IF(iscale==0) iscale=10**6
317 rlat1=float(igdtmpli(12))/float(iscale)
318 rlat2=float(igdtmpli(15))/float(iscale)
319 dlat=(rlat2-rlat1)/(jm-1)
320 jg=nint(180/abs(dlat))
321 IF(jm.EQ.jg) idrti=256
322 IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
323 ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
324 jg=igdtmpli(18)*2
325 IF(jm.NE.jg) iret=41
326 ENDIF
327 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
328 ! SET PARAMETERS
329 IF(iret.EQ.0) THEN
330 iromb=ipopt(1)
331 maxwv=ipopt(2)
332 IF(maxwv.EQ.-1) THEN
333 IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
334 IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
335 IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
336 IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
337 IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
338 IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
339 ENDIF
340 IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
341 ENDIF
342 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
343 ! INTERPOLATE
344 IF(iret.EQ.0) THEN
345 IF(nscan.EQ.0) THEN
346 iskipi=1
347 jskipi=im
348 ELSE
349 iskipi=jm
350 jskipi=1
351 ENDIF
352 IF(iscan.EQ.1) iskipi=-iskipi
353 IF(jscan.EQ.0) jskipi=-jskipi
354 ispec=0
355 ! SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
356 IF((igdtnumo.EQ.0.OR.igdtnumo.EQ.40).AND. &
357 mod(igdtmplo(8),2).EQ.0.AND.igdtmplo(13).EQ.0.AND.igdtmplo(19).EQ.0) THEN
358 idrto=igdtnumo
359 IF(idrto==40)idrto=4
360 imo=igdtmplo(8)
361 jmo=igdtmplo(9)
362 iscale=igdtmplo(10)*igdtmplo(11)
363 IF(iscale==0) iscale=10**6
364 rlon2=float(igdtmplo(16))/float(iscale)
365 dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
366 igo=nint(360/abs(dlono))
367 IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
368 rlat1=float(igdtmplo(12))/float(iscale)
369 rlat2=float(igdtmplo(15))/float(iscale)
370 dlat=(rlat2-rlat1)/(jmo-1)
371 jgo=nint(180/abs(dlat))
372 IF(jmo.EQ.jgo) idrto=256
373 IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
374 ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
375 jgo=igdtmplo(18)*2
376 IF(jmo.EQ.jgo) ispec=1
377 ENDIF
378 IF(ispec.EQ.1) THEN
379 CALL sptrun(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
380 km,iprime,iskipi,jskipi,mi,0,0,mo,0,gi,go)
381 ENDIF
382 ! SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
383 ELSEIF(igdtnumo.EQ.20.AND. &
384 igdtmplo(8).EQ.igdtmplo(9).AND.mod(igdtmplo(8),2).EQ.1.AND. &
385 igdtmplo(15).EQ.igdtmplo(16).AND.igdtmplo(18).EQ.64) THEN
386 nps=igdtmplo(8)
387 rlat1=float(igdtmplo(10))*1.e-6
388 rlon1=float(igdtmplo(11))*1.e-6
389 orient=float(igdtmplo(14))*1.e-6
390 xmesh=float(igdtmplo(15))*1.e-3
391 iproj=mod(igdtmplo(17)/128,2)
392 ip=(nps+1)/2
393 h=(-1.)**iproj
394 slat=float(abs(igdtmplo(13)))*1.e-6
395 CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
396 de=(1.+sin(slat/dpr))*rerth
397 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
398 xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
399 yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
400 IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
401 IF(iproj.EQ.0) THEN
402 CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
403 iprime,iskipi,jskipi,mi,mo,0,0,0, &
404 slat,xmesh,orient,gi,go,go2)
405 ELSE
406 CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
407 iprime,iskipi,jskipi,mi,mo,0,0,0, &
408 slat,xmesh,orient,gi,go2,go)
409 ENDIF
410 ispec=1
411 ENDIF
412 ! SPECIAL CASE OF MERCATOR GRID
413 ELSEIF(igdtnumo.EQ.10) THEN
414 ni=igdtmplo(8)
415 nj=igdtmplo(9)
416 rlat1=float(igdtmplo(10))*1.0e-6
417 rlon1=float(igdtmplo(11))*1.0e-6
418 rlon2=float(igdtmplo(15))*1.0e-6
419 rlati=float(igdtmplo(13))*1.0e-6
420 iscano=mod(igdtmplo(16)/128,2)
421 jscano=mod(igdtmplo(16)/64,2)
422 nscano=mod(igdtmplo(16)/32,2)
423 dy=float(igdtmplo(19))*1.0e-3
424 hi=(-1.)**iscano
425 hj=(-1.)**(1-jscano)
426 CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
427 dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
428 dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
429 IF(nscano.EQ.0) THEN
430 CALL sptrunm(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
431 iprime,iskipi,jskipi,mi,mo,0,0,0, &
432 rlat1,rlon1,dlato,dlono,gi,go)
433 ispec=1
434 ENDIF
435 ENDIF
436 ! GENERAL SLOW CASE
437 IF(ispec.EQ.0) THEN
438 CALL sptrung(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
439 iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon,gi,go)
440 ENDIF
441 DO k=1,km
442 ibo(k)=0
443 DO n=1,no
444 lo(n,k)=.true.
445 ENDDO
446 ENDDO
447 ELSE
448 DO k=1,km
449 ibo(k)=1
450 DO n=1,no
451 lo(n,k)=.false.
452 go(n,k)=0.
453 ENDDO
454 ENDDO
455 ENDIF
456 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
457 END SUBROUTINE polates4_grib2
458
459
463 !
536 SUBROUTINE polates4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,GI, &
537 NO,RLAT,RLON,IBO,LO,GO,IRET)
538 INTEGER, INTENT(IN ) :: IPOPT(20), KGDSI(200)
539 INTEGER, INTENT(IN ) :: KGDSO(200), MI, MO
540 INTEGER, INTENT(IN ) :: IBI(KM), KM
541 INTEGER, INTENT( OUT) :: IBO(KM), IRET
542 !
543 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
544 !
545 REAL, INTENT(IN ) :: GI(MI,KM)
546 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
547 REAL, INTENT( OUT) :: GO(MO,KM)
548 !
549 REAL, PARAMETER :: FILL=-9999.
550 REAL, PARAMETER :: RERTH=6.3712e6
551 REAL, PARAMETER :: PI=3.14159265358979
552 REAL, PARAMETER :: DPR=180./pi
553 !
554 INTEGER :: IDRTI, IDRTO, IG, JG, IM, JM
555 INTEGER :: IGO, JGO, IMO, JMO
556 INTEGER :: ISCAN, JSCAN, NSCAN
557 INTEGER :: ISCANO, JSCANO, NSCANO
558 INTEGER :: ISKIPI, JSKIPI
559 INTEGER :: IMAXI, JMAXI, ISPEC
560 INTEGER :: IP, IPRIME, IPROJ, IROMB, K
561 INTEGER :: MAXWV, N, NI, NJ, NPS, NO
562 !
563 REAL :: DE, DR, DY
564 REAL :: DLAT, DLON, DLATO, DLONO
565 REAL :: GO2(MO,KM), H, HI, HJ
566 REAL :: ORIENT
567 REAL :: RLAT1, RLON1, RLAT2, RLON2, RLATI
568 REAL :: XMESH, XP, YP
569 REAL :: XPTS(MO), YPTS(MO)
570
571 type(grib1_descriptor) :: desc_in, desc_out
572 class(ip_grid), allocatable :: grid_in, grid_out
573
574 desc_in = init_descriptor(kgdsi)
575 desc_out = init_descriptor(kgdso)
576
577 call init_grid(grid_in, desc_in)
578 call init_grid(grid_out, desc_out)
579
580 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
581 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
582 iret=0
583 IF(kgdso(1).GE.0) THEN
584 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no)
585 IF(no.EQ.0) iret=3
586 ENDIF
587 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
588 ! AFFIRM APPROPRIATE INPUT GRID
589 ! LAT/LON OR GAUSSIAN
590 ! NO BITMAPS
591 ! FULL ZONAL COVERAGE
592 ! FULL MERIDIONAL COVERAGE
593 idrti=kgdsi(1)
594 im=kgdsi(2)
595 jm=kgdsi(3)
596 rlon1=kgdsi(5)*1.e-3
597 rlon2=kgdsi(8)*1.e-3
598 iscan=mod(kgdsi(11)/128,2)
599 jscan=mod(kgdsi(11)/64,2)
600 nscan=mod(kgdsi(11)/32,2)
601 IF(idrti.NE.0.AND.idrti.NE.4) iret=41
602 DO k=1,km
603 IF(ibi(k).NE.0) iret=41
604 ENDDO
605 IF(iret.EQ.0) THEN
606 IF(iscan.EQ.0) THEN
607 dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
608 ELSE
609 dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
610 ENDIF
611 ig=nint(360/abs(dlon))
612 iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
613 imaxi=ig
614 jmaxi=jm
615 IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
616 ENDIF
617 IF(iret.EQ.0.AND.idrti.EQ.0) THEN
618 rlat1=kgdsi(4)*1.e-3
619 rlat2=kgdsi(7)*1.e-3
620 dlat=(rlat2-rlat1)/(jm-1)
621 jg=nint(180/abs(dlat))
622 IF(jm.EQ.jg) idrti=256
623 IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
624 ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
625 jg=kgdsi(10)*2
626 IF(jm.NE.jg) iret=41
627 ENDIF
628 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
629 ! SET PARAMETERS
630 IF(iret.EQ.0) THEN
631 iromb=ipopt(1)
632 maxwv=ipopt(2)
633 IF(maxwv.EQ.-1) THEN
634 IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
635 IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
636 IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
637 IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
638 IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
639 IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
640 ENDIF
641 IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
642 ENDIF
643 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
644 ! INTERPOLATE
645 IF(iret.EQ.0) THEN
646 IF(nscan.EQ.0) THEN
647 iskipi=1
648 jskipi=im
649 ELSE
650 iskipi=jm
651 jskipi=1
652 ENDIF
653 IF(iscan.EQ.1) iskipi=-iskipi
654 IF(jscan.EQ.0) jskipi=-jskipi
655 ispec=0
656 ! SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
657 IF((kgdso(1).EQ.0.OR.kgdso(1).EQ.4).AND. &
658 mod(kgdso(2),2).EQ.0.AND.kgdso(5).EQ.0.AND.kgdso(11).EQ.0) THEN
659 idrto=kgdso(1)
660 imo=kgdso(2)
661 jmo=kgdso(3)
662 rlon2=kgdso(8)*1.e-3
663 dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
664 igo=nint(360/abs(dlono))
665 IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
666 rlat1=kgdso(4)*1.e-3
667 rlat2=kgdso(7)*1.e-3
668 dlat=(rlat2-rlat1)/(jmo-1)
669 jgo=nint(180/abs(dlat))
670 IF(jmo.EQ.jgo) idrto=256
671 IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
672 ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
673 jgo=kgdso(10)*2
674 IF(jmo.EQ.jgo) ispec=1
675 ENDIF
676 IF(ispec.EQ.1) THEN
677 CALL sptrun(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
678 km,iprime,iskipi,jskipi,mi,0,0,mo,0,gi,go)
679 ENDIF
680 ! SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
681 ELSEIF(kgdso(1).EQ.5.AND. &
682 kgdso(2).EQ.kgdso(3).AND.mod(kgdso(2),2).EQ.1.AND. &
683 kgdso(8).EQ.kgdso(9).AND.kgdso(11).EQ.64) THEN
684 nps=kgdso(2)
685 rlat1=kgdso(4)*1.e-3
686 rlon1=kgdso(5)*1.e-3
687 orient=kgdso(7)*1.e-3
688 xmesh=kgdso(8)
689 iproj=mod(kgdso(10)/128,2)
690 ip=(nps+1)/2
691 h=(-1.)**iproj
692 de=(1.+sin(60./dpr))*rerth
693 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
694 xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
695 yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
696 IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
697 IF(iproj.EQ.0) THEN
698 CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
699 iprime,iskipi,jskipi,mi,mo,0,0,0, &
700 60.,xmesh,orient,gi,go,go2)
701 ELSE
702 CALL sptruns(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
703 iprime,iskipi,jskipi,mi,mo,0,0,0, &
704 60.,xmesh,orient,gi,go2,go)
705 ENDIF
706 ispec=1
707 ENDIF
708 ! SPECIAL CASE OF MERCATOR GRID
709 ELSEIF(kgdso(1).EQ.1) THEN
710 ni=kgdso(2)
711 nj=kgdso(3)
712 rlat1=kgdso(4)*1.e-3
713 rlon1=kgdso(5)*1.e-3
714 rlon2=kgdso(8)*1.e-3
715 rlati=kgdso(9)*1.e-3
716 iscano=mod(kgdso(11)/128,2)
717 jscano=mod(kgdso(11)/64,2)
718 nscano=mod(kgdso(11)/32,2)
719 dy=kgdso(13)
720 hi=(-1.)**iscano
721 hj=(-1.)**(1-jscano)
722 dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
723 dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
724 IF(nscano.EQ.0) THEN
725 CALL sptrunm(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
726 iprime,iskipi,jskipi,mi,mo,0,0,0, &
727 rlat1,rlon1,dlato,dlono,gi,go)
728 ispec=1
729 ENDIF
730 ENDIF
731 ! GENERAL SLOW CASE
732 IF(ispec.EQ.0) THEN
733 CALL sptrung(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
734 iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon,gi,go)
735 ENDIF
736 DO k=1,km
737 ibo(k)=0
738 DO n=1,no
739 lo(n,k)=.true.
740 ENDDO
741 ENDDO
742 ELSE
743 DO k=1,km
744 ibo(k)=1
745 DO n=1,no
746 lo(n,k)=.false.
747 go(n,k)=0.
748 ENDDO
749 ENDDO
750 ENDIF
751 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
752 END SUBROUTINE polates4_grib1
753
754
755 SUBROUTINE polatev4_grib2(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
756 IGDTNUMO,IGDTMPLO,IGDTLENO, &
757 MI,MO,KM,IBI,UI,VI, &
758 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
759 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
760 !
761 ! SUBPROGRAM: POLATEV4 INTERPOLATE VECTOR FIELDS (SPECTRAL)
762 ! PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
763 !
764 ! ABSTRACT: THIS SUBPROGRAM PERFORMS SPECTRAL INTERPOLATION
765 ! FROM ANY GRID TO ANY GRID FOR VECTOR FIELDS.
766 ! IT REQUIRES THAT THE INPUT FIELDS BE UNIFORMLY GLOBAL.
767 ! OPTIONS ALLOW CHOICES BETWEEN TRIANGULAR SHAPE (IPOPT(1)=0)
768 ! AND RHOMBOIDAL SHAPE (IPOPT(1)=1) WHICH HAS NO DEFAULT;
769 ! A SECOND OPTION IS THE TRUNCATION (IPOPT(2)) WHICH DEFAULTS
770 ! TO A SENSIBLE TRUNCATION FOR THE INPUT GRID (IF OPT(2)=-1).
771 ! NOTE THAT IF THE OUTPUT GRID IS NOT FOUND IN A SPECIAL LIST,
772 ! THEN THE TRANSFORM BACK TO GRID IS NOT VERY FAST.
773 ! THIS SPECIAL LIST CONTAINS GLOBAL CYLINDRICAL GRIDS,
774 ! POLAR STEREOGRAPHIC GRIDS CENTERED AT THE POLE
775 ! AND MERCATOR GRIDS. ONLY HORIZONTAL INTERPOLATION
776 ! IS PERFORMED.
777 !
778 ! THE INPUT AND OUTPUT GRIDS ARE DEFINED BY THEIR GRIB 2 GRID
779 ! DEFINITION TEMPLATE AS DECODED BY THE NCEP G2 LIBRARY. THE
780 ! CODE RECOGNIZES THE FOLLOWING PROJECTIONS, WHERE
781 ! "IGDTNUMI/O" IS THE GRIB 2 GRID DEFINTION TEMPLATE NUMBER
782 ! FOR THE INPUT AND OUTPUT GRIDS, RESPECTIVELY:
783 ! (IGDTNUMI/O=00) EQUIDISTANT CYLINDRICAL
784 ! (IGDTNUMO =01) ROTATED EQUIDISTANT CYLINDRICAL. "E" AND
785 ! NON-"E" STAGGERED
786 ! (IGDTNUMO =10) MERCATOR CYLINDRICAL
787 ! (IGDTNUMO =20) POLAR STEREOGRAPHIC AZIMUTHAL
788 ! (IGDTNUMO =30) LAMBERT CONFORMAL CONICAL
789 ! (IGDTNUMI/O=40) GAUSSIAN CYLINDRICAL
790 !
791 ! THE INPUT AND OUTPUT VECTORS ARE ROTATED SO THAT THEY ARE
792 ! EITHER RESOLVED RELATIVE TO THE DEFINED GRID
793 ! IN THE DIRECTION OF INCREASING X AND Y COORDINATES
794 ! OR RESOLVED RELATIVE TO EASTERLY AND NORTHERLY DIRECTIONS,
795 ! AS DESIGNATED BY THEIR RESPECTIVE GRID DESCRIPTION SECTIONS.
796 !
797 ! AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
798 ! AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED
799 ! ALONG WITH THEIR VECTOR ROTATION PARAMETERS.
800 ! ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
801 ! IF IGDTNUMO<0, IN WHICH CASE THE NUMBER OF POINTS
802 ! AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT
803 ! ALONG WITH THEIR VECTOR ROTATION PARAMETERS.
804 !
805 ! OUTPUT BITMAPS WILL ONLY BE CREATED WHEN THE OUTPUT GRID
806 ! EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
807 ! THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
808 !
809 ! PROGRAM HISTORY LOG:
810 ! 96-04-10 IREDELL
811 ! 2001-06-18 IREDELL IMPROVE DETECTION OF SPECIAL FAST TRANSFORM
812 ! 2015-01-27 GAYNO REPLACE CALLS TO GDSWIZ WITH NEW MERGED
813 ! ROUTINE GDSWZD.
814 ! 2015-07-13 GAYNO CONVERT TO GRIB 2. REPLACE GRIB 1 KGDS ARRAYS
815 ! WITH GRIB 2 GRID DEFINITION TEMPLATE ARRAYS.
816 !
817 ! USAGE: CALL POLATEV4(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
818 ! IGDTNUMO,IGDTMPLO,IGDTLENO, &
819 ! MI,MO,KM,IBI,LI,UI,VI, &
820 ! NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
821 !
822 ! INPUT ARGUMENT LIST:
823 ! IPOPT - INTEGER (20) INTERPOLATION OPTIONS
824 ! IPOPT(1)=0 FOR TRIANGULAR, IPOPT(1)=1 FOR RHOMBOIDAL;
825 ! IPOPT(2) IS TRUNCATION NUMBER
826 ! (DEFAULTS TO SENSIBLE IF IPOPT(2)=-1).
827 ! IGDTNUMI - INTEGER GRID DEFINITION TEMPLATE NUMBER - INPUT GRID.
828 ! CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
829 ! NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE:
830 ! 00 - EQUIDISTANT CYLINDRICAL
831 ! 01 - ROTATED EQUIDISTANT CYLINDRICAL. "E"
832 ! AND NON-"E" STAGGERED
833 ! 10 - MERCATOR CYCLINDRICAL
834 ! 20 - POLAR STEREOGRAPHIC AZIMUTHAL
835 ! 30 - LAMBERT CONFORMAL CONICAL
836 ! 40 - GAUSSIAN EQUIDISTANT CYCLINDRICAL
837 ! IGDTMPLI - INTEGER (IGDTLENI) GRID DEFINITION TEMPLATE ARRAY -
838 ! INPUT GRID. CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
839 ! OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE
840 ! (SECTION 3 INFO). SEE COMMENTS IN ROUTINE
841 ! IPOLATEV FOR COMPLETE DEFINITION.
842 ! IGDTLENI - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
843 ! TEMPLATE ARRAY - INPUT GRID. CORRESPONDS TO THE GFLD%IGDTLEN
844 ! COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
845 ! IGDTNUMO - INTEGER GRID DEFINITION TEMPLATE NUMBER - OUTPUT GRID.
846 ! CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
847 ! NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE. IGDTNUMO<0
848 ! MEANS INTERPOLATE TO RANDOM STATION POINTS.
849 ! OTHERWISE, SAME DEFINITION AS "IGDTNUMI".
850 ! IGDTMPLO - INTEGER (IGDTLENO) GRID DEFINITION TEMPLATE ARRAY -
851 ! OUTPUT GRID. CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
852 ! OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE
853 ! (SECTION 3 INFO). SEE COMMENTS IN ROUTINE
854 ! IPOLATEV FOR COMPLETE DEFINITION.
855 ! IGDTLENO - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
856 ! TEMPLATE ARRAY - OUTPUT GRID. CORRESPONDS TO THE GFLD%IGDTLEN
857 ! COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
858 ! MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
859 ! OR DIMENSION OF INPUT GRID FIELDS IF KM=1
860 ! MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
861 ! OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
862 ! KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE
863 ! IBI - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
864 ! UI - REAL (MI,KM) INPUT U-COMPONENT FIELDS TO INTERPOLATE
865 ! VI - REAL (MI,KM) INPUT V-COMPONENT FIELDS TO INTERPOLATE
866 ! RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF IGDTNUMO<0)
867 ! RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF IGDTNUMO<0)
868 ! CROT - REAL (MO) VECTOR ROTATION COSINES (IF IGDTNUMO<0)
869 ! SROT - REAL (MO) VECTOR ROTATION SINES (IF IGDTNUMO<0)
870 ! (UGRID=CROT*UEARTH-SROT*VEARTH;
871 ! VGRID=SROT*UEARTH+CROT*VEARTH)
872 !
873 ! OUTPUT ARGUMENT LIST:
874 ! NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF IGDTNUMO>=0)
875 ! RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF IGDTNUMO>=0)
876 ! RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF IGDTNUMO>=0)
877 ! CROT - REAL (MO) VECTOR ROTATION COSINES (IF IGDTNUMO>=0)
878 ! SROT - REAL (MO) VECTOR ROTATION SINES (IF IGDTNUMO>=0)
879 ! (UGRID=CROT*UEARTH-SROT*VEARTH;
880 ! VGRID=SROT*UEARTH+CROT*VEARTH)
881 ! IBO - INTEGER (KM) OUTPUT BITMAP FLAGS
882 ! LO - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
883 ! UO - REAL (MO,KM) OUTPUT U-COMPONENT FIELDS INTERPOLATED
884 ! VO - REAL (MO,KM) OUTPUT V-COMPONENT FIELDS INTERPOLATED
885 ! IRET - INTEGER RETURN CODE
886 ! 0 SUCCESSFUL INTERPOLATION
887 ! 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
888 ! 3 UNRECOGNIZED OUTPUT GRID
889 ! 41 INVALID NONGLOBAL INPUT GRID
890 ! 42 INVALID SPECTRAL METHOD PARAMETERS
891 !
892 ! SUBPROGRAMS CALLED:
893 ! EARTH_RADIUS DETERMINE SIZE/SHAPE OF EARTH
894 ! GDSWZD GRID DESCRIPTION SECTION WIZARD
895 ! SPTRUNV SPECTRALLY TRUNCATE GRIDDED VECTOR FIELDS
896 ! SPTRUNSV SPECTRALLY INTERPOLATE VECTORS TO POLAR STEREO.
897 ! SPTRUNMV SPECTRALLY INTERPOLATE VECTORS TO MERCATOR
898 ! SPTRUNGV SPECTRALLY INTERPOLATE VECTORS TO STATIONS
899 !
900 ! ATTRIBUTES:
901 ! LANGUAGE: FORTRAN 90
902 !
903 !$$$
904 INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
905 INTEGER, INTENT(IN ) :: KM, MI, MO
906 INTEGER, INTENT( OUT) :: IRET, IBO(KM), NO
907 INTEGER, INTENT(IN ) :: IGDTNUMI, IGDTLENI
908 INTEGER, INTENT(IN ) :: IGDTMPLI(IGDTLENI)
909 INTEGER, INTENT(IN ) :: IGDTNUMO, IGDTLENO
910 INTEGER, INTENT(IN ) :: IGDTMPLO(IGDTLENO)
911 !
912 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
913 !
914 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
915 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
916 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
917 REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
918 !
919 REAL, PARAMETER :: FILL=-9999.
920 REAL, PARAMETER :: PI=3.14159265358979
921 REAL, PARAMETER :: DPR=180./pi
922 !
923 INTEGER :: IDRTO, IROMB, ISKIPI, ISPEC
924 INTEGER :: IDRTI, IMAXI, JMAXI, IM, JM
925 INTEGER :: IPRIME, IG, IMO, JMO, IGO, JGO
926 INTEGER :: ISCAN, JSCAN, NSCAN
927 INTEGER :: ISCANO, JSCANO, NSCANO
928 INTEGER :: ISCALE, IP, IPROJ, JSKIPI, JG
929 INTEGER :: K, MAXWV, N, NI, NJ, NPS
930 !
931 REAL :: DLAT, DLON, DLATO, DLONO, DE, DR, DY
932 REAL :: DUM, E2, H, HI, HJ
933 REAL :: ORIENT, RERTH, SLAT
934 REAL :: RLAT1, RLON1, RLAT2, RLON2, RLATI
935 REAL :: UROT, VROT, UO2(MO,KM),VO2(MO,KM)
936 REAL :: XMESH, X, XP, YP, XPTS(MO),YPTS(MO)
937
938 type(grib2_descriptor) :: desc_in, desc_out
939 class(ip_grid), allocatable :: grid_in, grid_out
940
941 desc_in = init_descriptor(igdtnumi, igdtleni, igdtmpli)
942 desc_out = init_descriptor(igdtnumo, igdtleno, igdtmplo)
943
944 call init_grid(grid_in, desc_in)
945 call init_grid(grid_out, desc_out)
946
947 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
948 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
949 iret=0
950 IF(igdtnumo.GE.0) THEN
951 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts, &
952 rlon,rlat,no,crot,srot)
953 IF(no.EQ.0) iret=3
954 ENDIF
955 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
956 ! AFFIRM APPROPRIATE INPUT GRID
957 ! LAT/LON OR GAUSSIAN
958 ! NO BITMAPS
959 ! FULL ZONAL COVERAGE
960 ! FULL MERIDIONAL COVERAGE
961 idrti=igdtnumi
962 IF(idrti==40) idrti=4
963 IF(idrti==0.OR.idrti==4)THEN
964 im=igdtmpli(8)
965 jm=igdtmpli(9)
966 iscale=igdtmpli(10)*igdtmpli(11)
967 IF(iscale==0) iscale=10**6
968 rlon1=float(igdtmpli(13))/float(iscale)
969 rlon2=float(igdtmpli(16))/float(iscale)
970 iscan=mod(igdtmpli(19)/128,2)
971 jscan=mod(igdtmpli(19)/64,2)
972 nscan=mod(igdtmpli(19)/32,2)
973 ELSE
974 iret=41
975 ENDIF
976 DO k=1,km
977 IF(ibi(k).NE.0) iret=41
978 ENDDO
979 IF(iret.EQ.0) THEN
980 IF(iscan.EQ.0) THEN
981 dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
982 ELSE
983 dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
984 ENDIF
985 ig=nint(360/abs(dlon))
986 iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
987 imaxi=ig
988 jmaxi=jm
989 IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
990 ENDIF
991 IF(iret.EQ.0.AND.idrti.EQ.0) THEN
992 iscale=igdtmpli(10)*igdtmpli(11)
993 IF(iscale==0) iscale=10**6
994 rlat1=float(igdtmpli(12))/float(iscale)
995 rlat2=float(igdtmpli(15))/float(iscale)
996 dlat=(rlat2-rlat1)/(jm-1)
997 jg=nint(180/abs(dlat))
998 IF(jm.EQ.jg) idrti=256
999 IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
1000 ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
1001 jg=igdtmpli(18)*2
1002 IF(jm.NE.jg) iret=41
1003 ENDIF
1004 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1005 ! SET PARAMETERS
1006 IF(iret.EQ.0) THEN
1007 iromb=ipopt(1)
1008 maxwv=ipopt(2)
1009 IF(maxwv.EQ.-1) THEN
1010 IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
1011 IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
1012 IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
1013 IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
1014 IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
1015 IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
1016 ENDIF
1017 IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
1018 ENDIF
1019 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1020 ! INTERPOLATE
1021 IF(iret.EQ.0) THEN
1022 IF(nscan.EQ.0) THEN
1023 iskipi=1
1024 jskipi=im
1025 ELSE
1026 iskipi=jm
1027 jskipi=1
1028 ENDIF
1029 IF(iscan.EQ.1) iskipi=-iskipi
1030 IF(jscan.EQ.0) jskipi=-jskipi
1031 ispec=0
1032 ! SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
1033 IF((igdtnumo.EQ.0.OR.igdtnumo.EQ.40).AND. &
1034 mod(igdtmplo(8),2).EQ.0.AND.igdtmplo(13).EQ.0.AND. &
1035 igdtmplo(19).EQ.0) THEN
1036 idrto=igdtnumo
1037 IF(idrto==40)idrto=4
1038 imo=igdtmplo(8)
1039 jmo=igdtmplo(9)
1040 iscale=igdtmplo(10)*igdtmplo(11)
1041 IF(iscale==0) iscale=10**6
1042 rlon2=float(igdtmplo(16))/float(iscale)
1043 dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
1044 igo=nint(360/abs(dlono))
1045 IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
1046 rlat1=float(igdtmplo(12))/float(iscale)
1047 rlat2=float(igdtmplo(15))/float(iscale)
1048 dlat=(rlat2-rlat1)/(jmo-1)
1049 jgo=nint(180/abs(dlat))
1050 IF(jmo.EQ.jgo) idrto=256
1051 IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
1052 ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
1053 jgo=igdtmplo(18)*2
1054 IF(jmo.EQ.jgo) ispec=1
1055 ENDIF
1056 IF(ispec.EQ.1) THEN
1057 CALL sptrunv(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
1058 km,iprime,iskipi,jskipi,mi,0,0,mo,0,ui,vi, &
1059 .true.,uo,vo,.false.,dum,dum,.false.,dum,dum)
1060 ENDIF
1061 ! SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
1062 ELSEIF(igdtnumo.EQ.20.AND. &
1063 igdtmplo(8).EQ.igdtmplo(9).AND.mod(igdtmplo(8),2).EQ.1.AND. &
1064 igdtmplo(15).EQ.igdtmplo(16).AND.igdtmplo(18).EQ.64.AND. &
1065 mod(igdtmplo(12)/8,2).EQ.1) THEN
1066 nps=igdtmplo(8)
1067 rlat1=float(igdtmplo(10))*1.e-6
1068 rlon1=float(igdtmplo(11))*1.e-6
1069 orient=float(igdtmplo(14))*1.e-6
1070 xmesh=float(igdtmplo(15))*1.e-3
1071 iproj=mod(igdtmplo(17)/128,2)
1072 ip=(nps+1)/2
1073 h=(-1.)**iproj
1074 slat=float(abs(igdtmplo(13)))*1.e-6
1075 CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
1076 de=(1.+sin(slat/dpr))*rerth
1077 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
1078 xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
1079 yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
1080 IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
1081 IF(iproj.EQ.0) THEN
1082 CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
1083 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1084 slat,xmesh,orient,ui,vi,.true.,uo,vo,uo2,vo2, &
1085 .false.,dum,dum,dum,dum, &
1086 .false.,dum,dum,dum,dum)
1087 ELSE
1088 CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
1089 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1090 slat,xmesh,orient,ui,vi,.true.,uo2,vo2,uo,vo, &
1091 .false.,dum,dum,dum,dum, &
1092 .false.,dum,dum,dum,dum)
1093 ENDIF
1094 ispec=1
1095 ENDIF
1096 ! SPECIAL CASE OF MERCATOR GRID
1097 ELSEIF(igdtnumo.EQ.10) THEN
1098 ni=igdtmplo(8)
1099 nj=igdtmplo(9)
1100 rlat1=float(igdtmplo(10))*1.0e-6
1101 rlon1=float(igdtmplo(11))*1.0e-6
1102 rlon2=float(igdtmplo(15))*1.0e-6
1103 rlati=float(igdtmplo(13))*1.0e-6
1104 iscano=mod(igdtmplo(16)/128,2)
1105 jscano=mod(igdtmplo(16)/64,2)
1106 nscano=mod(igdtmplo(16)/32,2)
1107 dy=float(igdtmplo(19))*1.0e-3
1108 hi=(-1.)**iscano
1109 hj=(-1.)**(1-jscano)
1110 CALL earth_radius(igdtmplo,igdtleno,rerth,e2)
1111 dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
1112 dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
1113 IF(nscano.EQ.0) THEN
1114 CALL sptrunmv(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
1115 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1116 rlat1,rlon1,dlato,dlono,ui,vi, &
1117 .true.,uo,vo,.false.,dum,dum,.false.,dum,dum)
1118 ispec=1
1119 ENDIF
1120 ENDIF
1121 ! GENERAL SLOW CASE
1122 IF(ispec.EQ.0) THEN
1123 CALL sptrungv(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
1124 iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon, &
1125 ui,vi,.true.,uo,vo,.false.,x,x,.false.,x,x)
1126 DO k=1,km
1127 ibo(k)=0
1128 DO n=1,no
1129 lo(n,k)=.true.
1130 urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
1131 vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
1132 uo(n,k)=urot
1133 vo(n,k)=vrot
1134 ENDDO
1135 ENDDO
1136 ENDIF
1137 ELSE
1138 DO k=1,km
1139 ibo(k)=1
1140 DO n=1,no
1141 lo(n,k)=.false.
1142 uo(n,k)=0.
1143 vo(n,k)=0.
1144 ENDDO
1145 ENDDO
1146 ENDIF
1147 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1148 END SUBROUTINE polatev4_grib2
1149
1153 !
1241 SUBROUTINE polatev4_grib1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,UI,VI, &
1242 NO,RLAT,RLON,CROT,SROT,IBO,LO,UO,VO,IRET)
1243 INTEGER, INTENT(IN ) :: IPOPT(20), IBI(KM)
1244 INTEGER, INTENT(IN ) :: KM, MI, MO
1245 INTEGER, INTENT( OUT) :: IRET, IBO(KM)
1246 INTEGER, INTENT(IN) :: KGDSI(200),KGDSO(200)
1247 !
1248 LOGICAL*1, INTENT( OUT) :: LO(MO,KM)
1249 !
1250 REAL, INTENT(IN ) :: UI(MI,KM),VI(MI,KM)
1251 REAL, INTENT( OUT) :: UO(MO,KM),VO(MO,KM)
1252 REAL, INTENT(INOUT) :: RLAT(MO),RLON(MO)
1253 REAL, INTENT( OUT) :: CROT(MO),SROT(MO)
1254 !
1255 REAL, PARAMETER :: FILL=-9999.
1256 REAL, PARAMETER :: RERTH=6.3712e6
1257 REAL, PARAMETER :: PI=3.14159265358979
1258 REAL, PARAMETER :: DPR=180./pi
1259 !
1260 INTEGER :: IDRTO, IROMB, ISKIPI, ISPEC
1261 INTEGER :: IDRTI, IMAXI, JMAXI, IM, JM
1262 INTEGER :: IPRIME, IG, IMO, JMO, IGO, JGO
1263 INTEGER :: ISCAN, JSCAN, NSCAN
1264 INTEGER :: ISCANO, JSCANO, NSCANO
1265 INTEGER :: IP, IPROJ, JSKIPI, JG
1266 INTEGER :: K, MAXWV, N, NI, NJ, NO, NPS
1267 !
1268 REAL :: DLAT, DLON, DLATO, DLONO, DE, DR, DY
1269 REAL :: DUM, H, HI, HJ
1270 REAL :: ORIENT
1271 REAL :: RLAT1, RLON1, RLAT2, RLON2, RLATI
1272 REAL :: UROT, VROT, UO2(MO,KM),VO2(MO,KM)
1273 REAL :: XMESH, X, XP, YP, XPTS(MO),YPTS(MO)
1274
1275 type(grib1_descriptor) :: desc_in, desc_out
1276 class(ip_grid), allocatable :: grid_in, grid_out
1277
1278 desc_in = init_descriptor(kgdsi)
1279 desc_out = init_descriptor(kgdso)
1280
1281 call init_grid(grid_in, desc_in)
1282 call init_grid(grid_out, desc_out)
1283 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1284 ! COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
1285 iret=0
1286 IF(kgdso(1).GE.0) THEN
1287 CALL gdswzd(grid_out, 0,mo,fill,xpts,ypts,rlon,rlat,no,crot,srot)
1288 IF(no.EQ.0) iret=3
1289 ENDIF
1290 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1291 ! AFFIRM APPROPRIATE INPUT GRID
1292 ! LAT/LON OR GAUSSIAN
1293 ! NO BITMAPS
1294 ! FULL ZONAL COVERAGE
1295 ! FULL MERIDIONAL COVERAGE
1296 idrti=kgdsi(1)
1297 im=kgdsi(2)
1298 jm=kgdsi(3)
1299 rlon1=kgdsi(5)*1.e-3
1300 rlon2=kgdsi(8)*1.e-3
1301 iscan=mod(kgdsi(11)/128,2)
1302 jscan=mod(kgdsi(11)/64,2)
1303 nscan=mod(kgdsi(11)/32,2)
1304 IF(idrti.NE.0.AND.idrti.NE.4) iret=41
1305 DO k=1,km
1306 IF(ibi(k).NE.0) iret=41
1307 ENDDO
1308 IF(iret.EQ.0) THEN
1309 IF(iscan.EQ.0) THEN
1310 dlon=(mod(rlon2-rlon1-1+3600,360.)+1)/(im-1)
1311 ELSE
1312 dlon=-(mod(rlon1-rlon2-1+3600,360.)+1)/(im-1)
1313 ENDIF
1314 ig=nint(360/abs(dlon))
1315 iprime=1+mod(-nint(rlon1/dlon)+ig,ig)
1316 imaxi=ig
1317 jmaxi=jm
1318 IF(mod(ig,2).NE.0.OR.im.LT.ig) iret=41
1319 ENDIF
1320 IF(iret.EQ.0.AND.idrti.EQ.0) THEN
1321 rlat1=kgdsi(4)*1.e-3
1322 rlat2=kgdsi(7)*1.e-3
1323 dlat=(rlat2-rlat1)/(jm-1)
1324 jg=nint(180/abs(dlat))
1325 IF(jm.EQ.jg) idrti=256
1326 IF(jm.NE.jg.AND.jm.NE.jg+1) iret=41
1327 ELSEIF(iret.EQ.0.AND.idrti.EQ.4) THEN
1328 jg=kgdsi(10)*2
1329 IF(jm.NE.jg) iret=41
1330 ENDIF
1331 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1332 ! SET PARAMETERS
1333 IF(iret.EQ.0) THEN
1334 iromb=ipopt(1)
1335 maxwv=ipopt(2)
1336 IF(maxwv.EQ.-1) THEN
1337 IF(iromb.EQ.0.AND.idrti.EQ.4) maxwv=(jmaxi-1)
1338 IF(iromb.EQ.1.AND.idrti.EQ.4) maxwv=(jmaxi-1)/2
1339 IF(iromb.EQ.0.AND.idrti.EQ.0) maxwv=(jmaxi-3)/2
1340 IF(iromb.EQ.1.AND.idrti.EQ.0) maxwv=(jmaxi-3)/4
1341 IF(iromb.EQ.0.AND.idrti.EQ.256) maxwv=(jmaxi-1)/2
1342 IF(iromb.EQ.1.AND.idrti.EQ.256) maxwv=(jmaxi-1)/4
1343 ENDIF
1344 IF((iromb.NE.0.AND.iromb.NE.1).OR.maxwv.LT.0) iret=42
1345 ENDIF
1346 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1347 ! INTERPOLATE
1348 IF(iret.EQ.0) THEN
1349 IF(nscan.EQ.0) THEN
1350 iskipi=1
1351 jskipi=im
1352 ELSE
1353 iskipi=jm
1354 jskipi=1
1355 ENDIF
1356 IF(iscan.EQ.1) iskipi=-iskipi
1357 IF(jscan.EQ.0) jskipi=-jskipi
1358 ispec=0
1359 ! SPECIAL CASE OF GLOBAL CYLINDRICAL GRID
1360 IF((kgdso(1).EQ.0.OR.kgdso(1).EQ.4).AND. &
1361 mod(kgdso(2),2).EQ.0.AND.kgdso(5).EQ.0.AND. &
1362 kgdso(11).EQ.0) THEN
1363 idrto=kgdso(1)
1364 imo=kgdso(2)
1365 jmo=kgdso(3)
1366 rlon2=kgdso(8)*1.e-3
1367 dlono=(mod(rlon2-1+3600,360.)+1)/(imo-1)
1368 igo=nint(360/abs(dlono))
1369 IF(imo.EQ.igo.AND.idrto.EQ.0) THEN
1370 rlat1=kgdso(4)*1.e-3
1371 rlat2=kgdso(7)*1.e-3
1372 dlat=(rlat2-rlat1)/(jmo-1)
1373 jgo=nint(180/abs(dlat))
1374 IF(jmo.EQ.jgo) idrto=256
1375 IF(jmo.EQ.jgo.OR.jmo.EQ.jgo+1) ispec=1
1376 ELSEIF(imo.EQ.igo.AND.idrto.EQ.4) THEN
1377 jgo=kgdso(10)*2
1378 IF(jmo.EQ.jgo) ispec=1
1379 ENDIF
1380 IF(ispec.EQ.1) THEN
1381 CALL sptrunv(iromb,maxwv,idrti,imaxi,jmaxi,idrto,imo,jmo, &
1382 km,iprime,iskipi,jskipi,mi,0,0,mo,0,ui,vi, &
1383 .true.,uo,vo,.false.,dum,dum,.false.,dum,dum)
1384 ENDIF
1385 ! SPECIAL CASE OF POLAR STEREOGRAPHIC GRID
1386 ELSEIF(kgdso(1).EQ.5.AND. &
1387 kgdso(2).EQ.kgdso(3).AND.mod(kgdso(2),2).EQ.1.AND. &
1388 kgdso(8).EQ.kgdso(9).AND.kgdso(11).EQ.64.AND. &
1389 mod(kgdso(6)/8,2).EQ.1) THEN
1390 nps=kgdso(2)
1391 rlat1=kgdso(4)*1.e-3
1392 rlon1=kgdso(5)*1.e-3
1393 orient=kgdso(7)*1.e-3
1394 xmesh=kgdso(8)
1395 iproj=mod(kgdso(10)/128,2)
1396 ip=(nps+1)/2
1397 h=(-1.)**iproj
1398 de=(1.+sin(60./dpr))*rerth
1399 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
1400 xp=1-h*sin((rlon1-orient)/dpr)*dr/xmesh
1401 yp=1+cos((rlon1-orient)/dpr)*dr/xmesh
1402 IF(nint(xp).EQ.ip.AND.nint(yp).EQ.ip) THEN
1403 IF(iproj.EQ.0) THEN
1404 CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
1405 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1406 60.,xmesh,orient,ui,vi,.true.,uo,vo,uo2,vo2, &
1407 .false.,dum,dum,dum,dum, &
1408 .false.,dum,dum,dum,dum)
1409 ELSE
1410 CALL sptrunsv(iromb,maxwv,idrti,imaxi,jmaxi,km,nps, &
1411 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1412 60.,xmesh,orient,ui,vi,.true.,uo2,vo2,uo,vo, &
1413 .false.,dum,dum,dum,dum, &
1414 .false.,dum,dum,dum,dum)
1415 ENDIF
1416 ispec=1
1417 ENDIF
1418 ! SPECIAL CASE OF MERCATOR GRID
1419 ELSEIF(kgdso(1).EQ.1) THEN
1420 ni=kgdso(2)
1421 nj=kgdso(3)
1422 rlat1=kgdso(4)*1.e-3
1423 rlon1=kgdso(5)*1.e-3
1424 rlon2=kgdso(8)*1.e-3
1425 rlati=kgdso(9)*1.e-3
1426 iscano=mod(kgdso(11)/128,2)
1427 jscano=mod(kgdso(11)/64,2)
1428 nscano=mod(kgdso(11)/32,2)
1429 dy=kgdso(13)
1430 hi=(-1.)**iscano
1431 hj=(-1.)**(1-jscano)
1432 dlono=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(ni-1)
1433 dlato=hj*dy/(rerth*cos(rlati/dpr))*dpr
1434 IF(nscano.EQ.0) THEN
1435 CALL sptrunmv(iromb,maxwv,idrti,imaxi,jmaxi,km,ni,nj, &
1436 iprime,iskipi,jskipi,mi,mo,0,0,0, &
1437 rlat1,rlon1,dlato,dlono,ui,vi, &
1438 .true.,uo,vo,.false.,dum,dum,.false.,dum,dum)
1439 ispec=1
1440 ENDIF
1441 ENDIF
1442 ! GENERAL SLOW CASE
1443 IF(ispec.EQ.0) THEN
1444 CALL sptrungv(iromb,maxwv,idrti,imaxi,jmaxi,km,no, &
1445 iprime,iskipi,jskipi,mi,mo,0,0,0,rlat,rlon, &
1446 ui,vi,.true.,uo,vo,.false.,x,x,.false.,x,x)
1447 DO k=1,km
1448 ibo(k)=0
1449 DO n=1,no
1450 lo(n,k)=.true.
1451 urot=crot(n)*uo(n,k)-srot(n)*vo(n,k)
1452 vrot=srot(n)*uo(n,k)+crot(n)*vo(n,k)
1453 uo(n,k)=urot
1454 vo(n,k)=vrot
1455 ENDDO
1456 ENDDO
1457 ENDIF
1458 ELSE
1459 DO k=1,km
1460 ibo(k)=1
1461 DO n=1,no
1462 lo(n,k)=.false.
1463 uo(n,k)=0.
1464 vo(n,k)=0.
1465 ENDDO
1466 ENDDO
1467 ENDIF
1468 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1469 END SUBROUTINE polatev4_grib1
1470
1471
1472
1473end module spectral_interp_mod
gdswzd_mod::gdswzd
Definition: gdswzd_mod.f90:37
ip_grid_descriptor_mod::init_descriptor
Definition: ip_grid_descriptor_mod.f90:57
ip_grid_factory_mod::init_grid
Definition: ip_grid_factory_mod.f90:20
spectral_interp_mod::interpolate_spectral
Definition: spectral_interp_mod.f90:11
spectral_interp_mod::polates4
Definition: spectral_interp_mod.f90:16
spectral_interp_mod::polatev4
Definition: spectral_interp_mod.f90:21
gdswzd_mod
Driver module for gdswzd routines.
Definition: gdswzd_mod.f90:25
ip_grid_descriptor_mod
Uses derived type grid descriptor objects to abstract away the raw Grib-1 and Grib-2 grid definitions...
Definition: ip_grid_descriptor_mod.f90:15
ip_grid_factory_mod
Routines for creating an ip_grid given a Grib descriptor.
Definition: ip_grid_factory_mod.f90:11
ip_grid_mod
Abstract ip_grid type.
Definition: ip_grid_mod.f90:8
ip_grid_descriptor_mod::grib1_descriptor
Descriptor representing a grib1 grib descriptor section (GDS) with an integer array.
Definition: ip_grid_descriptor_mod.f90:37
ip_grid_descriptor_mod::grib2_descriptor
Grib-2 descriptor containing a grib2 GDT represented by an integer array.
Definition: ip_grid_descriptor_mod.f90:45
ip_grid_mod::ip_grid
Abstract grid that holds fields and methods common to all grids.
Definition: ip_grid_mod.f90:45