NCEPLIBS-ip 4.0.0
ip_rot_equid_cylind_egrid_mod.f90
Go to the documentation of this file.
1
5
10module ip_rot_equid_cylind_egrid_mod
11 use iso_fortran_env, only: real64
12 use ip_grid_descriptor_mod
13 use ip_grid_mod
14 use constants_mod, only: dpr, pi
15 use earth_radius_mod
16 implicit none
17
18 private
19 public :: ip_rot_equid_cylind_egrid
20
21 integer, parameter :: kd = real64
22
23 type, extends(ip_grid) :: ip_rot_equid_cylind_egrid
24 real(kd) :: rlon0, rlon1, rlat1, clat0, slat0
25 real(kd) :: dlats, dlons, hi
26 integer :: irot
27 contains
28 procedure :: init_grib1
29 procedure :: init_grib2
30 procedure :: gdswzd => gdswzd_rot_equid_cylind_egrid
31 end type ip_rot_equid_cylind_egrid
32
33 INTEGER :: IROT
34
35 REAL(KIND=kd) :: clat, clat0, clatr
36 REAL(KIND=kd) :: clon, dlats, dlons, rerth
37 REAL(KIND=kd) :: rlon0, slat, slat0, slatr
38
39contains
40
48 subroutine init_grib1(self, g1_desc)
49 class(ip_rot_equid_cylind_egrid), intent(inout) :: self
50 type(grib1_descriptor), intent(in) :: g1_desc
51
52 integer :: iscan
53 real(kd) :: rlat0
54
55 real(kd) :: rlat1, rlon1, rlon0, slat1, clat1, slat0, clat0, clon1
56 real(kd) :: slatr, clatr, clonr, rlatr, rlonr, dlats, dlons, hs, hi
57 integer :: im, jm
58
59 integer :: is1, kscan, irot
60
61 associate(kgds => g1_desc%gds)
62 self%rerth = 6.3712e6_kd
63 self%eccen_squared = 0.0
64
65 im=kgds(2)
66 jm=kgds(3)
67
68 self%nscan_field_pos = 3
69 self%nscan = mod(kgds(11)/32,2)
70
71 rlat1=kgds(4)*1.e-3_kd
72 rlon1=kgds(5)*1.e-3_kd
73 rlat0=kgds(7)*1.e-3_kd
74 rlon0=kgds(8)*1.e-3_kd
75
76 irot=mod(kgds(6)/8,2)
77 kscan=mod(kgds(11)/256,2)
78 iscan=mod(kgds(11)/128,2)
79 hi=(-1.)**iscan
80 slat1=sin(rlat1/dpr)
81 clat1=cos(rlat1/dpr)
82 slat0=sin(rlat0/dpr)
83 clat0=cos(rlat0/dpr)
84 hs=sign(1._kd,mod(rlon1-rlon0+180+3600,360._kd)-180)
85 clon1=cos((rlon1-rlon0)/dpr)
86 slatr=clat0*slat1-slat0*clat1*clon1
87 clatr=sqrt(1-slatr**2)
88 clonr=(clat0*clat1*clon1+slat0*slat1)/clatr
89 rlatr=dpr*asin(slatr)
90 rlonr=hs*dpr*acos(clonr)
91 dlats=rlatr/(-(jm-1)/2)
92 dlons=rlonr/(-((im * 2 - 1) -1)/2)
93
94 IF(kscan.EQ.0) THEN
95 is1=(jm+1)/2
96 ELSE
97 is1=jm/2
98 ENDIF
99
100 self%im = im
101 self%jm = jm
102 self%rlon0 = rlon0
103 self%rlon1 = rlon1
104 self%rlat1 = rlat1
105 self%clat0 = clat0
106 self%slat0 = slat0
107 self%dlats = dlats
108 self%dlons = dlons
109 self%hi = hi
110 self%irot = irot
111 self%kscan = kscan
112
113
114 end associate
115
116 end subroutine init_grib1
117
118
125 subroutine init_grib2(self, g2_desc)
126 class(ip_rot_equid_cylind_egrid), intent(inout) :: self
127 type(grib2_descriptor), intent(in) :: g2_desc
128
129 integer :: iscale, iscan
130 real(kd) :: rlat0
131 integer :: i_offset_odd!, i_offset_even
132
133 associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
134 CALL earth_radius(igdtmpl,igdtlen,self%rerth,self%eccen_squared)
135
136 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137 ! ROUTINE ONLY WORKS FOR "E"-STAGGER GRIDS.
138 ! "V" GRID WHEN BIT 5 IS '1' AND BIT 6 IS '0'.
139 ! "H" GRID WHEN BIT 5 IS '0' AND BIT 6 IS '1'.
140 ! I_OFFSET_ODD=MOD(IGDTMPL(19)/8,2)
141 ! I_OFFSET_EVEN=MOD(IGDTMPL(19)/4,2)
142 ! IF(I_OFFSET_ODD==I_OFFSET_EVEN) THEN
143 ! CALL ROT_EQUID_CYLIND_EGRID_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
144 ! RETURN
145 ! ENDIF
146 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147
148 self%IM=igdtmpl(8)
149 self%JM=igdtmpl(9)
150
151 self%NSCAN=mod(igdtmpl(16)/32,2)
152 self%nscan_field_pos = 3
153
154 iscale=igdtmpl(10)*igdtmpl(11)
155 IF(iscale==0) iscale=10**6
156
157 self%RLON0=float(igdtmpl(21))/float(iscale)
158 self%DLATS=float(igdtmpl(18))/float(iscale)
159 ! THE GRIB2 CONVENTION FOR "I" RESOLUTION IS TWICE WHAT THIS ROUTINE ASSUMES.
160 self%DLONS=float(igdtmpl(17))/float(iscale) * 0.5_kd
161
162 self%IROT=mod(igdtmpl(14)/8,2)
163
164 i_offset_odd=mod(igdtmpl(19)/8,2)
165 self%KSCAN=i_offset_odd
166 iscan=mod(igdtmpl(19)/128,2)
167
168 self%HI=(-1.)**iscan
169
170 rlat0=float(igdtmpl(20))/float(iscale)
171 rlat0=rlat0+90.0_kd
172
173 self%SLAT0=sin(rlat0/dpr)
174 self%CLAT0=cos(rlat0/dpr)
175
176 self%RLAT1=float(igdtmpl(12))/float(iscale)
177 self%RLON1=float(igdtmpl(13))/float(iscale)
178 end associate
179 end subroutine init_grib2
180
181
229 SUBROUTINE gdswzd_rot_equid_cylind_egrid(self,IOPT,NPTS,&
230 FILL,XPTS,YPTS,RLON,RLAT,NRET, &
231 CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
232 IMPLICIT NONE
233 !
234 class(ip_rot_equid_cylind_egrid), intent(in) :: self
235
236 INTEGER, INTENT(IN ) :: IOPT, NPTS
237 INTEGER, INTENT( OUT) :: NRET
238 !
239 REAL, INTENT(IN ) :: FILL
240 REAL, INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
241 REAL, INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
242 REAL, OPTIONAL, INTENT( OUT) :: CROT(NPTS),SROT(NPTS)
243 REAL, OPTIONAL, INTENT( OUT) :: XLON(NPTS),XLAT(NPTS)
244 REAL, OPTIONAL, INTENT( OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
245 !
246 INTEGER :: IM, JM, IS1, N
247 INTEGER :: KSCAN
248 ! INTEGER :: I_OFFSET_ODD, I_OFFSET_EVEN
249 !
250 LOGICAL :: LROT, LMAP, LAREA
251 !
252 REAL(KIND=kd) :: rlat1, rlon1
253 REAL(KIND=kd) :: clonr
254 REAL(KIND=kd) :: rlatr, rlonr, sbd, wbd, hs
255 REAL :: HI
256 REAL :: XMAX, XMIN, YMAX, YMIN, XPTF, YPTF
257 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258 IF(PRESENT(crot)) crot=fill
259 IF(PRESENT(srot)) srot=fill
260 IF(PRESENT(xlon)) xlon=fill
261 IF(PRESENT(xlat)) xlat=fill
262 IF(PRESENT(ylon)) ylon=fill
263 IF(PRESENT(ylat)) ylat=fill
264 IF(PRESENT(area)) area=fill
265
266 rlon0=self%rlon0
267 irot=self%irot
268 im=self%im * 2 - 1
269 jm=self%jm
270 dlats=self%dlats
271 dlons=self%dlons
272 kscan=self%kscan
273 hi=self%hi
274 slat0=self%slat0
275 clat0=self%clat0
276 rlat1=self%rlat1
277 rlon1=self%rlon1
278
279 rerth = self%rerth
280
281 ! IS THE EARTH RADIUS DEFINED?
282 IF(rerth<0.)THEN
283 CALL rot_equid_cylind_egrid_error(iopt,fill,rlat,rlon,xpts,ypts,npts)
284 RETURN
285 ENDIF
286
287 sbd=rlat1
288 wbd=rlon1
289
290 IF (wbd > 180.0) wbd = wbd - 360.0
291 IF(kscan.EQ.0) THEN
292 is1=(jm+1)/2
293 ELSE
294 is1=jm/2
295 ENDIF
296
297 xmin=0
298 xmax=im+2
299 ymin=0
300 ymax=jm+1
301 nret=0
302
303 IF(PRESENT(crot).AND.PRESENT(srot))THEN
304 lrot=.true.
305 ELSE
306 lrot=.false.
307 ENDIF
308 IF(PRESENT(xlon).AND.PRESENT(xlat).AND.PRESENT(ylon).AND.PRESENT(ylat))THEN
309 lmap=.true.
310 ELSE
311 lmap=.false.
312 ENDIF
313 IF(PRESENT(area))THEN
314 larea=.true.
315 ELSE
316 larea=.false.
317 ENDIF
318
319 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
320 ! TRANSLATE GRID COORDINATES TO EARTH COORDINATES
321 IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
322 DO n=1,npts
323 xptf=ypts(n)+(xpts(n)-is1)
324 yptf=ypts(n)-(xpts(n)-is1)+kscan
325 IF(xptf.GE.xmin.AND.xptf.LE.xmax.AND. &
326 yptf.GE.ymin.AND.yptf.LE.ymax) THEN
327 hs=hi*sign(1.,xptf-(im+1)/2)
328 select type(desc => self%descriptor)
329 type is(grib1_descriptor)
330 rlonr=(xptf-(im+1)/2)*dlons
331 rlatr=(yptf-(jm+1)/2)*dlats
332 type is(grib2_descriptor)
333 rlonr=(xptf-1.0_kd)*dlons + wbd
334 rlatr=(yptf-1.0_kd)*dlats + sbd
335 end select
336 clonr=cos(rlonr/dpr)
337 slatr=sin(rlatr/dpr)
338 clatr=cos(rlatr/dpr)
339 slat=clat0*slatr+slat0*clatr*clonr
340 IF(slat.LE.-1) THEN
341 clat=0.
342 clon=cos(rlon0/dpr)
343 rlon(n)=0
344 rlat(n)=-90
345 ELSEIF(slat.GE.1) THEN
346 clat=0.
347 clon=cos(rlon0/dpr)
348 rlon(n)=0
349 rlat(n)=90
350 ELSE
351 clat=sqrt(1-slat**2)
352 clon=(clat0*clatr*clonr-slat0*slatr)/clat
353 clon=min(max(clon,-1._kd),1._kd)
354 rlon(n)=mod(rlon0+hs*dpr*acos(clon)+3600,360._kd)
355 rlat(n)=dpr*asin(slat)
356 ENDIF
357 nret=nret+1
358 IF(lrot) CALL rot_equid_cylind_egrid_vect_rot(rlon(n),crot(n),srot(n))
359 IF(lmap) CALL rot_equid_cylind_egrid_map_jacob(fill, rlon(n), &
360 xlon(n),xlat(n),ylon(n),ylat(n))
361 IF(larea) CALL rot_equid_cylind_egrid_grid_area(fill, area(n))
362 ELSE
363 rlon(n)=fill
364 rlat(n)=fill
365 ENDIF
366 ENDDO
367 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
368 ! TRANSLATE EARTH COORDINATES TO GRID COORDINATES
369 ELSEIF(iopt.EQ.-1) THEN
370 DO n=1,npts
371 IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90) THEN
372 hs=sign(1._kd,mod(rlon(n)-rlon0+180+3600,360._kd)-180)
373 clon=cos((rlon(n)-rlon0)/dpr)
374 slat=sin(rlat(n)/dpr)
375 clat=cos(rlat(n)/dpr)
376 slatr=clat0*slat-slat0*clat*clon
377 IF(slatr.LE.-1) THEN
378 clatr=0.
379 rlonr=0
380 rlatr=-90
381 ELSEIF(slatr.GE.1) THEN
382 clatr=0.
383 rlonr=0
384 rlatr=90
385 ELSE
386 clatr=sqrt(1-slatr**2)
387 clonr=(clat0*clat*clon+slat0*slat)/clatr
388 clonr=min(max(clonr,-1._kd),1._kd)
389 rlonr=hs*dpr*acos(clonr)
390 rlatr=dpr*asin(slatr)
391 ENDIF
392 select type(desc => self%descriptor)
393 type is(grib1_descriptor)
394 xptf=((rlonr-wbd)/dlons)+1.0_kd
395 yptf=((rlatr-sbd)/dlats)+1.0_kd
396 type is(grib2_descriptor)
397 xptf=(im+1)/2+rlonr/dlons
398 yptf=(jm+1)/2+rlatr/dlats
399 end select
400
401 IF(xptf.GE.xmin.AND.xptf.LE.xmax.AND. &
402 yptf.GE.ymin.AND.yptf.LE.ymax) THEN
403 xpts(n)=is1+(xptf-(yptf-kscan))/2
404 ypts(n)=(xptf+(yptf-kscan))/2
405 nret=nret+1
406 IF(lrot) CALL rot_equid_cylind_egrid_vect_rot(rlon(n),crot(n),srot(n))
407 IF(lmap) CALL rot_equid_cylind_egrid_map_jacob(fill, rlon(n), &
408 xlon(n),xlat(n),ylon(n),ylat(n))
409 IF(larea) CALL rot_equid_cylind_egrid_grid_area(fill, area(n))
410 ELSE
411 xpts(n)=fill
412 ypts(n)=fill
413 ENDIF
414 ELSE
415 xpts(n)=fill
416 ypts(n)=fill
417 ENDIF
418 ENDDO
419 ENDIF
420 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
421 END SUBROUTINE gdswzd_rot_equid_cylind_egrid
422 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
423 SUBROUTINE rot_equid_cylind_egrid_error(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
424 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
425 !
426 ! SUBPROGRAM: ROT_EQUID_CYLIND_EGRID_ERROR ERROR HANDLER
427 ! PRGMMR: GAYNO ORG: W/NMC23 DATE: 2015-07-13
428 !
429 ! ABSTRACT: UPON AN ERROR, THIS SUBPROGRAM ASSIGNS
430 ! A "FILL" VALUE TO THE OUTPUT FIELDS.
431
432 ! PROGRAM HISTORY LOG:
433 ! 2015-07-13 GAYNO INITIAL VERSION
434 ! 2015-09-17 GAYNO RENAME AS "ROT_EQUID_CYLIND_EGRID_ERROR"
435 !
436 ! USAGE: CALL ROT_EQUID_CYLIND_EGRID_ERROR(IOPT,FILL,RLAT,RLON,
437 ! XPTS,YPTS,NPTS)
438 !
439 ! INPUT ARGUMENT LIST:
440 ! IOPT - INTEGER OPTION FLAG
441 ! (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
442 ! (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
443 ! NPTS - INTEGER MAXIMUM NUMBER OF COORDINATES
444 ! FILL - REAL FILL VALUE TO SET INVALID OUTPUT DATA
445 ! (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
446 ! OUTPUT ARGUMENT LIST:
447 ! RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
448 ! RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
449 ! XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
450 ! YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
451 !
452 ! ATTRIBUTES:
453 ! LANGUAGE: FORTRAN 90
454 !
455 !$$$
456 IMPLICIT NONE
457 !
458 INTEGER, INTENT(IN ) :: IOPT, NPTS
459 !
460 REAL, INTENT(IN ) :: FILL
461 REAL, INTENT( OUT) :: RLAT(NPTS),RLON(NPTS)
462 REAL, INTENT( OUT) :: XPTS(NPTS),YPTS(NPTS)
463 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
464 IF(iopt>=0) THEN
465 rlon=fill
466 rlat=fill
467 ENDIF
468 IF(iopt<=0) THEN
469 xpts=fill
470 ypts=fill
471 ENDIF
472 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
473 END SUBROUTINE rot_equid_cylind_egrid_error
474
488 SUBROUTINE rot_equid_cylind_egrid_vect_rot(RLON, CROT, SROT)
489 IMPLICIT NONE
490
491 REAL , INTENT(IN ) :: RLON
492 REAL , INTENT( OUT) :: CROT, SROT
493
494 REAL(KIND=kd) :: slon
495
496 IF(irot.EQ.1) THEN
497 IF(clatr.LE.0) THEN
498 crot=-sign(1._kd,slatr*slat0)
499 srot=0.
500 ELSE
501 slon=sin((rlon-rlon0)/dpr)
502 crot=(clat0*clat+slat0*slat*clon)/clatr
503 srot=slat0*slon/clatr
504 ENDIF
505 ELSE
506 crot=1.
507 srot=0.
508 ENDIF
509
510 END SUBROUTINE rot_equid_cylind_egrid_vect_rot
511
523 SUBROUTINE rot_equid_cylind_egrid_map_jacob(FILL, RLON, &
524 XLON, XLAT, YLON, YLAT)
525 IMPLICIT NONE
526
527 REAL , INTENT(IN ) :: FILL, RLON
528 REAL , INTENT( OUT) :: XLON, XLAT, YLON, YLAT
529
530 REAL(KIND=kd) :: slon, term1, term2
531 REAL(KIND=kd) :: xlatf, xlonf, ylatf, ylonf
532
533 IF(clatr.LE.0._kd) THEN
534 xlon=fill
535 xlat=fill
536 ylon=fill
537 ylat=fill
538 ELSE
539 slon=sin((rlon-rlon0)/dpr)
540 term1=(clat0*clat+slat0*slat*clon)/clatr
541 term2=slat0*slon/clatr
542 xlonf=term1*clat/(dlons*clatr)
543 xlatf=-term2/(dlons*clatr)
544 ylonf=term2*clat/dlats
545 ylatf=term1/dlats
546 xlon=xlonf-ylonf
547 xlat=xlatf-ylatf
548 ylon=xlonf+ylonf
549 ylat=xlatf+ylatf
550 ENDIF
551
552 END SUBROUTINE rot_equid_cylind_egrid_map_jacob
553
561 SUBROUTINE rot_equid_cylind_egrid_grid_area(FILL, AREA)
562 IMPLICIT NONE
563
564 REAL, INTENT(IN ) :: FILL
565 REAL, INTENT( OUT) :: AREA
566
567 IF(clatr.LE.0._kd) THEN
568 area=fill
569 ELSE
570 area=rerth**2*clatr*dlats*dlons*2/dpr**2
571 ENDIF
572
573 END SUBROUTINE rot_equid_cylind_egrid_grid_area
574
575end module ip_rot_equid_cylind_egrid_mod
576
577
constants_mod
Module containing common constants.
Definition: constants_mod.f90:9
constants_mod::pi
real(real64), parameter pi
PI.
Definition: constants_mod.f90:15
constants_mod::dpr
real(real64), parameter dpr
Radians to degrees.
Definition: constants_mod.f90:16
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_rot_equid_cylind_egrid_mod
Rotated equidistant cylindrical grib decoder and grid coordinate transformations.
Definition: ip_rot_equid_cylind_egrid_mod.f90:10
ip_rot_equid_cylind_egrid_mod::rot_equid_cylind_egrid_map_jacob
subroutine rot_equid_cylind_egrid_map_jacob(FILL, RLON, XLON, XLAT, YLON, YLAT)
Computes the map jacobians for a rotated equidistant cylindrical grid.
Definition: ip_rot_equid_cylind_egrid_mod.f90:525
ip_rot_equid_cylind_egrid_mod::rot_equid_cylind_egrid_grid_area
subroutine rot_equid_cylind_egrid_grid_area(FILL, AREA)
Computes the grid box area for a rotated equidistant cylindrical grid.
Definition: ip_rot_equid_cylind_egrid_mod.f90:562
ip_rot_equid_cylind_egrid_mod::rot_equid_cylind_egrid_vect_rot
subroutine rot_equid_cylind_egrid_vect_rot(RLON, CROT, SROT)
Computes the vector rotation sines and cosines for a rotated equidistant cylindrical grid.
Definition: ip_rot_equid_cylind_egrid_mod.f90:489
ip_rot_equid_cylind_egrid_mod::gdswzd_rot_equid_cylind_egrid
subroutine gdswzd_rot_equid_cylind_egrid(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 rotated equidistant cylin...
Definition: ip_rot_equid_cylind_egrid_mod.f90:232
ip_rot_equid_cylind_egrid_mod::init_grib1
subroutine init_grib1(self, g1_desc)
Initializes a rotated equidistant cylindrical grid given a grib1_descriptor object.
Definition: ip_rot_equid_cylind_egrid_mod.f90:49
ip_rot_equid_cylind_egrid_mod::init_grib2
subroutine init_grib2(self, g2_desc)
Initializes a rotated equidistant cylindrical grid given a grib2_descriptor object.
Definition: ip_rot_equid_cylind_egrid_mod.f90:126
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
ip_rot_equid_cylind_egrid_mod::ip_rot_equid_cylind_egrid
Definition: ip_rot_equid_cylind_egrid_mod.f90:23