NCEPLIBS-ip 4.0.0
ip_gaussian_grid_mod.f90
Go to the documentation of this file.
1
5
10module ip_gaussian_grid_mod
11 use ip_grid_descriptor_mod
12 use ip_grid_mod
13 use earth_radius_mod
14 use constants_mod
15 implicit none
16
17 private
18 public :: ip_gaussian_grid
19
20 type, extends(ip_grid) :: ip_gaussian_grid
21 integer :: jh
22 real :: dlon, rlat1, rlon1, rlon2, hi
23 integer :: jg, jscan
24 contains
25 procedure :: init_grib1
26 procedure :: init_grib2
27 procedure :: gdswzd => gdswzd_gaussian
28 end type ip_gaussian_grid
29
30 INTEGER :: J1, JH
31 REAL, ALLOCATABLE :: BLAT(:)
32 REAL :: DLON, RERTH
33 REAL, ALLOCATABLE :: YLAT_ROW(:)
34
35contains
36
44 subroutine init_grib1(self, g1_desc)
45 class(ip_gaussian_grid), intent(inout) :: self
46 type(grib1_descriptor), intent(in) :: g1_desc
47
48 integer :: iscan, jg
49
50 associate(kgds => g1_desc%gds)
51 self%rerth = 6.3712e6
52 self%eccen_squared = 0.0
53
54 self%IM=kgds(2)
55 self%JM=kgds(3)
56 self%RLAT1=kgds(4)*1.e-3
57 self%RLON1=kgds(5)*1.e-3
58 self%RLON2=kgds(8)*1.e-3
59 self%JG=kgds(10)*2
60 iscan=mod(kgds(11)/128,2)
61 self%JSCAN=mod(kgds(11)/64,2)
62 self%HI=(-1.)**iscan
63 self%JH=(-1)**self%JSCAN
64 self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
65
66 self%iwrap = 0
67 self%jwrap1 = 0
68 self%jwrap2 = 0
69 self%nscan = mod(kgds(11) / 32, 2)
70 self%nscan_field_pos = self%nscan
71 self%kscan = 0
72
73 self%iwrap=nint(360 / abs(self%dlon))
74 if(self%im < self%iwrap) self%iwrap = 0
75
76 if(self%iwrap > 0 .and. mod(self%iwrap, 2) == 0) then
77 jg=kgds(10)*2
78 if(self%jm == self%jg) then
79 self%jwrap1 = 1
80 self%jwrap2 = 2 * self%jm + 1
81 endif
82 endif
83
84 end associate
85 end subroutine init_grib1
86
93 subroutine init_grib2(self, g2_desc)
94 class(ip_gaussian_grid), intent(inout) :: self
95 type(grib2_descriptor), intent(in) :: g2_desc
96
97 integer :: iscale, iscan, jg
98
99 associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
100 call earth_radius(igdtmpl, igdtlen, self%rerth, self%eccen_squared)
101
102 self%IM=igdtmpl(8)
103 self%JM=igdtmpl(9)
104 iscale=igdtmpl(10)*igdtmpl(11)
105 IF(iscale==0) iscale=10**6
106 self%RLAT1=float(igdtmpl(12))/float(iscale)
107 self%RLON1=float(igdtmpl(13))/float(iscale)
108 self%RLON2=float(igdtmpl(16))/float(iscale)
109 self%JG=igdtmpl(18)*2
110 iscan=mod(igdtmpl(19)/128,2)
111 self%JSCAN=mod(igdtmpl(19)/64,2)
112 self%HI=(-1.)**iscan
113 self%JH=(-1)**self%JSCAN
114 self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
115
116
117 self%iwrap = nint(360 / abs(self%dlon))
118 if(self%im < self%iwrap) self%iwrap = 0
119 self%jwrap1 = 0
120 self%jwrap2 = 0
121 if(self%iwrap > 0 .and. mod(self%iwrap, 2) == 0) then
122 jg = igdtmpl(18) * 2
123 if(self%jm == jg) then
124 self%jwrap1=1
125 self%jwrap2 = 2 * self%jm + 1
126 endif
127 endif
128 self%nscan = mod(igdtmpl(19) / 32, 2)
129 self%nscan_field_pos = self%nscan
130 self%kscan = 0
131 end associate
132
133 end subroutine init_grib2
134
179 SUBROUTINE gdswzd_gaussian(self,IOPT,NPTS,FILL, &
180 XPTS,YPTS,RLON,RLAT,NRET, &
181 CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
182 IMPLICIT NONE
183 !
184 class(ip_gaussian_grid), intent(in) :: self
185 INTEGER, INTENT(IN ) :: IOPT, NPTS
186 INTEGER, INTENT( OUT) :: NRET
187 !
188 REAL, INTENT(IN ) :: FILL
189 REAL, INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
190 REAL, INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
191 REAL, OPTIONAL, INTENT( OUT) :: CROT(NPTS),SROT(NPTS)
192 REAL, OPTIONAL, INTENT( OUT) :: XLON(NPTS),XLAT(NPTS)
193 REAL, OPTIONAL, INTENT( OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
194 !
195 INTEGER :: JSCAN, IM, JM
196 INTEGER :: J, JA, JG
197 INTEGER :: ISCALE, N
198 !
199 LOGICAL :: LROT, LMAP, LAREA
200 !
201 REAL, ALLOCATABLE :: ALAT(:), ALAT_JSCAN(:)
202 REAL, ALLOCATABLE :: ALAT_TEMP(:),BLAT_TEMP(:)
203 REAL :: HI, RLATA, RLATB, RLAT1, RLON1, RLON2
204 REAL :: XMAX, XMIN, YMAX, YMIN, YPTSA, YPTSB
205 REAL :: WB
206 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207 IF(PRESENT(crot)) crot=fill
208 IF(PRESENT(srot)) srot=fill
209 IF(PRESENT(xlon)) xlon=fill
210 IF(PRESENT(xlat)) xlat=fill
211 IF(PRESENT(ylon)) ylon=fill
212 IF(PRESENT(ylat)) ylat=fill
213 IF(PRESENT(area)) area=fill
214 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215
216 IF(PRESENT(crot).AND.PRESENT(srot))THEN
217 lrot=.true.
218 ELSE
219 lrot=.false.
220 ENDIF
221 IF(PRESENT(xlon).AND.PRESENT(xlat).AND.PRESENT(ylon).AND.PRESENT(ylat))THEN
222 lmap=.true.
223 ELSE
224 lmap=.false.
225 ENDIF
226 IF(PRESENT(area))THEN
227 larea=.true.
228 ELSE
229 larea=.false.
230 ENDIF
231
232 im=self%im
233 jm=self%jm
234
235 rlat1=self%rlat1
236 rlon1=self%rlon1
237 rlon2=self%rlon2
238
239 jg=self%jg
240 jscan=self%jscan
241 hi=self%hi
242
243 jh=self%jh
244 dlon=self%dlon
245 rerth = self%rerth
246
247 ALLOCATE(alat_temp(jg))
248 ALLOCATE(blat_temp(jg))
249 CALL splat(4,jg,alat_temp,blat_temp)
250 ALLOCATE(alat(0:jg+1))
251 ALLOCATE(blat(0:jg+1))
252 !$OMP PARALLEL DO PRIVATE(JA) SCHEDULE(STATIC)
253 DO ja=1,jg
254 alat(ja)=dpr*asin(alat_temp(ja))
255 blat(ja)=blat_temp(ja)
256 ENDDO
257 !$OMP END PARALLEL DO
258 DEALLOCATE(alat_temp,blat_temp)
259 alat(0)=180.-alat(1)
260 alat(jg+1)=-alat(0)
261 blat(0)=-blat(1)
262 blat(jg+1)=blat(0)
263 j1=1
264 DO WHILE(j1.LT.jg.AND.rlat1.LT.(alat(j1)+alat(j1+1))/2)
265 j1=j1+1
266 ENDDO
267 IF(lmap)THEN
268 ALLOCATE(alat_jscan(jg))
269 DO ja=1,jg
270 alat_jscan(j1+jh*(ja-1))=alat(ja)
271 ENDDO
272 ALLOCATE(ylat_row(0:jg+1))
273 DO ja=2,(jg-1)
274 ylat_row(ja)=2.0/(alat_jscan(ja+1)-alat_jscan(ja-1))
275 ENDDO
276 ylat_row(1)=1.0/(alat_jscan(2)-alat_jscan(1))
277 ylat_row(0)=ylat_row(1)
278 ylat_row(jg)=1.0/(alat_jscan(jg)-alat_jscan(jg-1))
279 ylat_row(jg+1)=ylat_row(jg)
280 DEALLOCATE(alat_jscan)
281 ENDIF
282 xmin=0
283 xmax=im+1
284 IF(im.EQ.nint(360/abs(dlon))) xmax=im+2
285 ymin=0.5
286 ymax=jm+0.5
287 nret=0
288 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289 ! TRANSLATE GRID COORDINATES TO EARTH COORDINATES
290 IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
291 !$OMP PARALLEL DO PRIVATE(N,J,WB,RLATA,RLATB) REDUCTION(+:NRET) SCHEDULE(STATIC)
292 DO n=1,npts
293 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
294 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
295 rlon(n)=mod(rlon1+dlon*(xpts(n)-1)+3600,360.)
296 j=ypts(n)
297 wb=ypts(n)-j
298 rlata=alat(j1+jh*(j-1))
299 rlatb=alat(j1+jh*j)
300 rlat(n)=rlata+wb*(rlatb-rlata)
301 nret=nret+1
302 IF(lrot) CALL gaussian_vect_rot(crot(n),srot(n))
303 IF(lmap) CALL gaussian_map_jacob(ypts(n),&
304 xlon(n),xlat(n),ylon(n),ylat(n))
305 IF(larea) CALL gaussian_grid_area(ypts(n),area(n))
306 ELSE
307 rlon(n)=fill
308 rlat(n)=fill
309 ENDIF
310 ENDDO
311 !$OMP END PARALLEL DO
312 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
313 ! TRANSLATE EARTH COORDINATES TO GRID COORDINATES
314 ELSEIF(iopt.EQ.-1) THEN
315 !$OMP PARALLEL DO PRIVATE(N,JA,YPTSA, YPTSB, WB) REDUCTION(+:NRET) SCHEDULE(STATIC)
316 DO n=1,npts
317 xpts(n)=fill
318 ypts(n)=fill
319 IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90) THEN
320 xpts(n)=1+hi*mod(hi*(rlon(n)-rlon1)+3600,360.)/dlon
321 ja=min(int((jg+1)/180.*(90-rlat(n))),jg)
322 IF(rlat(n).GT.alat(ja)) ja=max(ja-2,0)
323 IF(rlat(n).LT.alat(ja+1)) ja=min(ja+2,jg)
324 IF(rlat(n).GT.alat(ja)) ja=ja-1
325 IF(rlat(n).LT.alat(ja+1)) ja=ja+1
326 yptsa=1+jh*(ja-j1)
327 yptsb=1+jh*(ja+1-j1)
328 wb=(alat(ja)-rlat(n))/(alat(ja)-alat(ja+1))
329 ypts(n)=yptsa+wb*(yptsb-yptsa)
330 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
331 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
332 nret=nret+1
333 IF(lrot) CALL gaussian_vect_rot(crot(n),srot(n))
334 IF(lmap) CALL gaussian_map_jacob(ypts(n), &
335 xlon(n),xlat(n),ylon(n),ylat(n))
336 IF(larea) CALL gaussian_grid_area(ypts(n),area(n))
337 ELSE
338 xpts(n)=fill
339 ypts(n)=fill
340 ENDIF
341 ENDIF
342 ENDDO
343 !$OMP END PARALLEL DO
344 ENDIF
345 DEALLOCATE(alat, blat)
346 IF (ALLOCATED(ylat_row)) DEALLOCATE(ylat_row)
347 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348 END SUBROUTINE gdswzd_gaussian
349
350 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351 SUBROUTINE gaussian_error(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
352 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
353 !
354 ! SUBPROGRAM: GAUSSIAN_ERROR ERROR HANDLER
355 ! PRGMMR: GAYNO ORG: W/NMC23 DATE: 2015-07-13
356 !
357 ! ABSTRACT: UPON AN ERROR, THIS SUBPROGRAM ASSIGNS
358 ! A "FILL" VALUE TO THE OUTPUT FIELDS.
359 !
360 ! PROGRAM HISTORY LOG:
361 ! 2015-07-13 GAYNO INITIAL VERSION
362 !
363 ! USAGE: CALL GAUSSIAN_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
364 !
365 ! INPUT ARGUMENT LIST:
366 ! IOPT - INTEGER OPTION FLAG
367 ! (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
368 ! (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
369 ! NPTS - INTEGER MAXIMUM NUMBER OF COORDINATES
370 ! FILL - REAL FILL VALUE TO SET INVALID OUTPUT DATA
371 ! (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
372 ! OUTPUT ARGUMENT LIST:
373 ! RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
374 ! RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
375 ! XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
376 ! YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
377 !
378 ! ATTRIBUTES:
379 ! LANGUAGE: FORTRAN 90
380 !
381 !$$$
382 IMPLICIT NONE
383 !
384 INTEGER, INTENT(IN ) :: IOPT, NPTS
385 !
386 REAL, INTENT(IN ) :: FILL
387 REAL, INTENT( OUT) :: RLAT(NPTS),RLON(NPTS)
388 REAL, INTENT( OUT) :: XPTS(NPTS),YPTS(NPTS)
389 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
390 IF(iopt>=0) THEN
391 rlon=fill
392 rlat=fill
393 ENDIF
394 IF(iopt<=0) THEN
395 xpts=fill
396 ypts=fill
397 ENDIF
398 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
399 END SUBROUTINE gaussian_error
400
413 SUBROUTINE gaussian_vect_rot(CROT,SROT)
414 IMPLICIT NONE
415
416 REAL, INTENT( OUT) :: CROT, SROT
417
418 crot=1.0
419 srot=0.0
420
421 END SUBROUTINE gaussian_vect_rot
422
423
434 SUBROUTINE gaussian_map_jacob(YPTS, XLON, XLAT, YLON, YLAT)
435 IMPLICIT NONE
436
437 REAL, INTENT(IN ) :: YPTS
438 REAL, INTENT( OUT) :: XLON, XLAT, YLON, YLAT
439
440 xlon=1/dlon
441 xlat=0.
442 ylon=0.
443 ylat=ylat_row(nint(ypts))
444
445 END SUBROUTINE gaussian_map_jacob
446
454 SUBROUTINE gaussian_grid_area(YPTS,AREA)
455 IMPLICIT NONE
456
457 REAL, INTENT(IN ) :: YPTS
458 REAL, INTENT( OUT) :: AREA
459
460 INTEGER :: J
461
462 REAL :: WB, WLAT, WLATA, WLATB
463
464 j = ypts
465 wb=ypts-j
466 wlata=blat(j1+jh*(j-1))
467 wlatb=blat(j1+jh*j)
468 wlat=wlata+wb*(wlatb-wlata)
469 area=rerth**2*wlat*dlon/dpr
470
471 END SUBROUTINE gaussian_grid_area
472
473
474end module ip_gaussian_grid_mod
475
constants_mod
Module containing common constants.
Definition: constants_mod.f90:9
ip_gaussian_grid_mod
Gaussian grid coordinate transformations.
Definition: ip_gaussian_grid_mod.f90:10
ip_gaussian_grid_mod::gdswzd_gaussian
subroutine gdswzd_gaussian(self, IOPT, NPTS, FILL, XPTS, YPTS, RLON, RLAT, NRET, CROT, SROT, XLON, XLAT, YLON, YLAT, AREA)
Calculates Earth coordinates (iopt = 1) or grid coorindates (iopt = -1) for Gaussian grids.
Definition: ip_gaussian_grid_mod.f90:182
ip_gaussian_grid_mod::init_grib1
subroutine init_grib1(self, g1_desc)
Initializes a gaussian grid given a grib1_descriptor object.
Definition: ip_gaussian_grid_mod.f90:45
ip_gaussian_grid_mod::gaussian_grid_area
subroutine gaussian_grid_area(YPTS, AREA)
Computes the grid box area for a gaussian cylindrical grid.
Definition: ip_gaussian_grid_mod.f90:455
ip_gaussian_grid_mod::gaussian_vect_rot
subroutine gaussian_vect_rot(CROT, SROT)
Computes the vector rotation sines and cosines for a gaussian cylindrical grid.
Definition: ip_gaussian_grid_mod.f90:414
ip_gaussian_grid_mod::gaussian_map_jacob
subroutine gaussian_map_jacob(YPTS, XLON, XLAT, YLON, YLAT)
Computes the map jacobians for a gaussian cylindrical grid.
Definition: ip_gaussian_grid_mod.f90:435
ip_gaussian_grid_mod::init_grib2
subroutine init_grib2(self, g2_desc)
Initializes a gaussian grid given a grib2_descriptor object.
Definition: ip_gaussian_grid_mod.f90:94
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_mod
Abstract ip_grid type.
Definition: ip_grid_mod.f90:8
ip_gaussian_grid_mod::ip_gaussian_grid
Definition: ip_gaussian_grid_mod.f90:20
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_mod::ip_grid
Abstract grid that holds fields and methods common to all grids.
Definition: ip_grid_mod.f90:45