NCEPLIBS-ip 5.3.0
All Data Structures Namespaces Files Functions Variables Pages
ip_rot_equid_cylind_grid_mod.F90
Go to the documentation of this file.
1!> @file
2!> @brief Rotated equidistant cylindrical GRIB decoder and grid
3!> coordinate transformations for Arakawa grids A through D.
4!>
5!> @author Mark Iredell, George Gayno, Kyle Gerheiser
6!> @date July 2021
7
8!> Rotated equidistant cylindrical GRIB decoder and grid coordinate
9!> transformations for Arakawa grids A through D. (To handle the E
10!> grid, see ip_rot_equid_cylind_egrid_mod).
11!>
12!> See more info about [Awakawa
13!> grids](https://en.wikipedia.org/wiki/Arakawa_grids).
14!>
15!> Octet numbers refer to [GRIB2 - GRID DEFINITION TEMPLATE 3.1 Rotate
16!> Latitude/Longitude](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-1.shtml).
17!>
18!> @author Gayno @date 2007-NOV-15
20 use iso_fortran_env, only: real64
22 use ip_grid_mod
23 use ip_constants_mod, only: dpr, pi
25 implicit none
26
27 private
29
30 integer, parameter :: kd = real64 !< Fortran kind for reals.
31
33 real(kd) :: clat0 !< Cosine of the latitude of the southern pole of projection.
34 real(kd) :: dlats !< 'J'-direction grid increment.
35 real(kd) :: dlons !< 'I'-direction grid increment.
36 real(kd) :: rlon0 !< Longitude of southern pole of projection.
37 real(kd) :: slat0 !< Sine of the latitude of the southern pole of projection.
38 real(kd) :: wbd !< Longitude of the western boundary of the grid before rotation.
39 real(kd) :: sbd !< Latitude of the southern boundary of the grid before rotation.
40 !> Rotation flag. When '0' the u/v vector components are relative
41 !> to north/east. When '1' the u/v vector components are grid
42 !> relative.
43 integer :: irot
44 contains
45 !> Initializes a Rotated equidistant cylindrical grid given a
46 !> grib1_descriptor object.
47 !> @return N/A @copydoc ip_rot_equid_cylind_grid_mod::init_grib1
48 procedure :: init_grib1
49 !> Initializes a Rotated equidistant cylindrical given a
50 !> grib2_descriptor object
51 !> @return N/A @copydoc ip_rot_equid_cylind_grid_mod::init_grib2
52 procedure :: init_grib2
53 !> Calculates Earth coordinates (iopt = 1) or grid coorindates (iopt = -1)
54 !> for Gaussian grids.
55 !> @return N/A @copydoc ip_rot_equid_cylind_grid_mod::gdswzd_rot_equid_cylind
58
59 INTEGER :: irot !< Local copy of irot.
60 REAL(kind=kd) :: rerth !< Radius of the Earth.
61 REAL(kind=kd) :: clat0 !< Local copy of clat0.
62 REAL(kind=kd) :: dlats !< Local copy of dlats.
63 REAL(kind=kd) :: dlons !< Local copy of dlons.
64 REAL(kind=kd) :: rlon0 !< Local copy of rlon0.
65 REAL(kind=kd) :: slat0 !< Local copy of slat0.
66
67CONTAINS
68
69 !> Initializes a Rotated equidistant cylindrical grid given a
70 !> grib1_descriptor object.
71 !>
72 !> @param[inout] self The grid to initialize
73 !> @param[in] g1_desc A grib1_descriptor
74 !>
75 !> @author Gayno @date 2007-NOV-15
76 subroutine init_grib1(self, g1_desc)
77 class(ip_rot_equid_cylind_grid), intent(inout) :: self
78 type(grib1_descriptor), intent(in) :: g1_desc
79
80 real(kd) :: rlat1, rlon1, rlat0, rlat2, rlon2, nbd, ebd
81 real(kd) :: hs, hs2, slat1, slat2, slatr, clon1, clon2, clat1, clat2, clatr, clonr, rlonr, rlatr
82
83 associate(kgds => g1_desc%gds)
84 self%rerth = 6.3712e6_kd
85 self%eccen_squared = 0d0
86
87 rlat1=kgds(4)*1.e-3_kd
88 rlon1=kgds(5)*1.e-3_kd
89 rlat0=kgds(7)*1.e-3_kd
90 self%RLON0=kgds(8)*1.e-3_kd
91 rlat2=kgds(12)*1.e-3_kd
92 rlon2=kgds(13)*1.e-3_kd
93
94 self%IROT=mod(kgds(6)/8,2)
95 self%IM=kgds(2)
96 self%JM=kgds(3)
97
98 slat1=sin(rlat1/dpr)
99 clat1=cos(rlat1/dpr)
100 self%SLAT0=sin(rlat0/dpr)
101 self%CLAT0=cos(rlat0/dpr)
102
103 hs=sign(1._kd,mod(rlon1-self%RLON0+180+3600,360._kd)-180)
104 clon1=cos((rlon1-self%RLON0)/dpr)
105 slatr=self%CLAT0*slat1-self%SLAT0*clat1*clon1
106 clatr=sqrt(1-slatr**2)
107 clonr=(self%CLAT0*clat1*clon1+self%SLAT0*slat1)/clatr
108 rlatr=dpr*asin(slatr)
109 rlonr=hs*dpr*acos(clonr)
110
111 self%WBD=rlonr
112 self%SBD=rlatr
113 slat2=sin(rlat2/dpr)
114 clat2=cos(rlat2/dpr)
115 hs2=sign(1._kd,mod(rlon2-self%RLON0+180+3600,360._kd)-180)
116 clon2=cos((rlon2-self%RLON0)/dpr)
117 slatr=self%CLAT0*slat2-self%SLAT0*clat2*clon2
118 clatr=sqrt(1-slatr**2)
119 clonr=(self%CLAT0*clat2*clon2+self%SLAT0*slat2)/clatr
120 nbd=dpr*asin(slatr)
121 ebd=hs2*dpr*acos(clonr)
122 self%DLATS=(nbd-self%SBD)/float(self%JM-1)
123 self%DLONS=(ebd-self%WBD)/float(self%IM-1)
124
125 self%iwrap = 0
126 self%jwrap1 = 0
127 self%jwrap2 = 0
128 self%nscan = mod(kgds(11) / 32, 2)
129 self%nscan_field_pos = self%nscan
130 self%kscan = 0
131 end associate
132
133 end subroutine init_grib1
134
135 !> Initializes a Rotated equidistant cylindrical grid given a
136 !> grib2_descriptor object. Call 'use_ncep_post_arakawa()' before
137 !> using this subroutine to use ncep_post-compatible grid definition.
138 !>
139 !> @param[inout] self The grid to initialize
140 !> @param[in] g2_desc A grib2_descriptor
141 !>
142 !> @author Alex Richert @date 2024-MAY-20
143 subroutine init_grib2(self, g2_desc)
144 class(ip_rot_equid_cylind_grid), intent(inout) :: self
145 type(grib2_descriptor), intent(in) :: g2_desc
146
147 if (ncep_post_arakawa.and.(g2_desc%gdt_num.eq.32769)) then
148 call init_grib2_ncep_post(self, g2_desc)
149 else
150 call init_grib2_default(self, g2_desc)
151 endif
152
153 end subroutine init_grib2
154
155 !> Initializes a Rotated equidistant cylindrical grid given a
156 !> grib2_descriptor object. Uses template definitions consistent
157 !> with WMO standards.
158 !>
159 !> @param[inout] self The grid to initialize
160 !> @param[in] g2_desc A grib2_descriptor
161 !>
162 !> @author Gayno @date 2007-NOV-15
163 subroutine init_grib2_default(self, g2_desc)
164 class(ip_rot_equid_cylind_grid), intent(inout) :: self
165 type(grib2_descriptor), intent(in) :: g2_desc
166
167 real(kd) :: rlat1, rlon1, rlat0, rlat2, rlon2, nbd, ebd
168 integer :: iscale
169 integer :: i_offset_odd, i_offset_even, j_offset
170
171 associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
172
173 CALL earth_radius(igdtmpl,igdtlen,self%rerth,self%eccen_squared)
174
175 i_offset_odd=mod(igdtmpl(19)/8,2)
176 i_offset_even=mod(igdtmpl(19)/4,2)
177 j_offset=mod(igdtmpl(19)/2,2)
178
179 iscale=igdtmpl(10)*igdtmpl(11)
180 IF(iscale==0) iscale=10**6
181
182 rlat1=float(igdtmpl(12))/float(iscale)
183 rlon1=float(igdtmpl(13))/float(iscale)
184 rlat0=float(igdtmpl(20))/float(iscale)
185 rlat0=rlat0+90.0_kd
186
187 self%RLON0=float(igdtmpl(21))/float(iscale)
188
189 rlat2=float(igdtmpl(15))/float(iscale)
190 rlon2=float(igdtmpl(16))/float(iscale)
191
192 self%IROT=mod(igdtmpl(14)/8,2)
193 self%IM=igdtmpl(8)
194 self%JM=igdtmpl(9)
195
196 self%SLAT0=sin(rlat0/dpr)
197 self%CLAT0=cos(rlat0/dpr)
198
199 self%WBD=rlon1
200 IF (self%WBD > 180.0) self%WBD = self%WBD - 360.0
201 self%SBD=rlat1
202
203 nbd=rlat2
204 ebd=rlon2
205
206 self%DLATS=(nbd-self%SBD)/float(self%JM-1)
207 self%DLONS=(ebd-self%WBD)/float(self%IM-1)
208
209 IF(i_offset_odd==1) self%WBD=self%WBD+(0.5_kd*self%DLONS)
210 IF(j_offset==1) self%SBD=self%SBD+(0.5_kd*self%DLATS)
211
212 self%iwrap = 0
213 self%jwrap1 = 0
214 self%jwrap2 = 0
215 self%kscan = 0
216 self%nscan = mod(igdtmpl(19) / 32, 2)
217 self%nscan_field_pos = self%nscan
218 end associate
219 end subroutine init_grib2_default
220
221 !> Initializes a Rotated equidistant cylindrical grid given a
222 !> grib2_descriptor object. Uses template number definitions in
223 !> a manner compatible with wgrib2 and ncep_post, which is based
224 !> on grib1.
225 !>
226 !> @param[inout] self The grid to initialize
227 !> @param[in] g2_desc A grib2_descriptor
228 !>
229 !> @author Alex Richert, George Gayno @date 2024-MAY-20
230 subroutine init_grib2_ncep_post(self, g2_desc)
231 class(ip_rot_equid_cylind_grid), intent(inout) :: self
232 type(grib2_descriptor), intent(in) :: g2_desc
233
234 real(kd) :: rlat1, rlon1, rlat0, rlat2, rlon2, nbd, ebd
235 integer :: iscale
236 integer :: i_offset_odd, i_offset_even, j_offset
237 real(kd) :: hs, hs2, slat1, slat2, slatr, clon1, clon2, clat1, clat2, clatr, clonr, rlonr, rlatr
238
239 associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
240
241 CALL earth_radius(igdtmpl,igdtlen,self%rerth,self%eccen_squared)
242
243 i_offset_odd=mod(igdtmpl(19)/8,2)
244 i_offset_even=mod(igdtmpl(19)/4,2)
245 j_offset=mod(igdtmpl(19)/2,2)
246
247 iscale=igdtmpl(10)*igdtmpl(11)
248 IF(iscale==0) iscale=10**6
249
250 rlat1=float(igdtmpl(12))/float(iscale)
251 rlon1=float(igdtmpl(13))/float(iscale)
252 rlat0=float(igdtmpl(15))/float(iscale)
253
254 self%RLON0=float(igdtmpl(16))/float(iscale)
255
256 rlat2=float(igdtmpl(20))/float(iscale)
257 rlon2=float(igdtmpl(21))/float(iscale)
258
259 self%IROT=mod(igdtmpl(14)/8,2)
260 self%IM=igdtmpl(8)
261 self%JM=igdtmpl(9)
262
263 slat1=sin(rlat1/dpr)
264 clat1=cos(rlat1/dpr)
265
266 self%SLAT0=sin(rlat0/dpr)
267 self%CLAT0=cos(rlat0/dpr)
268
269 hs=sign(1._kd,mod(rlon1-self%RLON0+180+3600,360._kd)-180)
270
271 clon1=cos((rlon1-self%RLON0)/dpr)
272 slatr=self%CLAT0*slat1-self%SLAT0*clat1*clon1
273 clatr=sqrt(1-slatr**2)
274 clonr=(self%CLAT0*clat1*clon1+self%SLAT0*slat1)/clatr
275 rlatr=dpr*asin(slatr)
276 rlonr=hs*dpr*acos(clonr)
277
278 self%WBD=rlonr
279 IF (self%WBD > 180.0) self%WBD = self%WBD - 360.0
280 self%SBD=rlatr
281
282 slat2=sin(rlat2/dpr)
283 clat2=cos(rlat2/dpr)
284 hs2=sign(1._kd,mod(rlon2-self%RLON0+180+3600,360._kd)-180)
285 clon2=cos((rlon2-self%RLON0)/dpr)
286 slatr=self%CLAT0*slat2-self%SLAT0*clat2*clon2
287 clatr=sqrt(1-slatr**2)
288 clonr=(self%CLAT0*clat2*clon2+self%SLAT0*slat2)/clatr
289 nbd=dpr*asin(slatr)
290 ebd=hs2*dpr*acos(clonr)
291 self%DLATS=(nbd-self%SBD)/float(self%JM-1)
292 self%DLONS=(ebd-self%WBD)/float(self%IM-1)
293
294 IF(i_offset_odd==1) self%WBD=self%WBD+(0.5_kd*self%DLONS)
295 IF(j_offset==1) self%SBD=self%SBD+(0.5_kd*self%DLATS)
296
297 self%iwrap = 0
298 self%jwrap1 = 0
299 self%jwrap2 = 0
300 self%kscan = 0
301 self%nscan = mod(igdtmpl(19) / 32, 2)
302 self%nscan_field_pos = self%nscan
303 end associate
304 end subroutine init_grib2_ncep_post
305
306 !> GDS wizard for rotated equidistant cylindrical.
307 !>
308 !> This subprogram decodes the grib 2 grid definition template
309 !> (passed in integer form as decoded by the ncep g2 library) and
310 !> returns one of the following:
311 !> - (iopt=+1) earth coordinates of selected grid coordinates
312 !> - (iopt=-1) grid coordinates of selected earth coordinates
313 !>
314 !> Works for non-"e" staggered rotated equidistant cylindrical
315 !> projections. the scan mode (section 3, octet 72, bits 5-6)
316 !> determine whether this is an "h" or "v" grid.
317 !>
318 !> If the selected coordinates are more than one gridpoint beyond
319 !> the the edges of the grid domain, then the relevant output
320 !> elements are set to fill values. The actual number of valid
321 !> points computed is returned too.
322 !>
323 !> Optionally, the vector rotations, the map jacobians and the grid
324 !> box areas may be returned as well.
325 !>
326 !> To compute the vector rotations, the optional arguments 'srot'
327 !> and 'crot' must be present. To compute the map jacobians, the
328 !> optional arguments 'xlon', 'xlat', 'ylon', 'ylat' must be
329 !> present. To compute the grid box areas, the optional argument
330 !> 'area' must be present.
331 !>
332 !> ### Program History Log
333 !> Date | Programmer | Comments
334 !> -----|------------|---------
335 !> 2010-jan-15 | gayno | based on routines gdswzdcb and gdswzdca
336 !> 2015-jan-21 | gayno | merger of gdswizcd and gdswzdcd. make crot,sort,xlon,xlat,ylon,ylat and area optional arguments. make part of a module. move vector rotation, map jacobian and grid box area computations to separate subroutines.
337 !> 2015-jul-13 | gayno | convert to grib 2. replace grib 1 kgds array with grib 2 grid definition template array. rename as "gdswzd_rot_equid_cylind."
338 !> 2018-07-20 | wesley | add threads.
339 !>
340 !> @param[in] self Module reference.
341 !> @param[in] iopt integer option flag
342 !> - 1 to compute earth coords of selected grid coords
343 !> - -1 to compute grid coords of selected earth coords
344 !> @param[in] npts integer maximum number of coordinates
345 !> @param[in] fill real fill value to set invalid output data
346 !> (must be impossible value; suggested value: -9999.)
347 !> @param[inout] xpts real (npts) grid x point coordinates if iopt>0
348 !> @param[inout] ypts real (npts) grid y point coordinates if iopt>0
349 !> @param[inout] rlon real (npts) earth longitudes in degrees e if iopt<0
350 !> (acceptable range: -360. to 360.)
351 !> @param[inout] rlat real (npts) earth latitudes in degrees n if iopt<0
352 !> (acceptable range: -90. to 90.)
353 !> @param[out] nret integer number of valid points computed
354 !> @param[out] crot real, optional (npts) clockwise vector rotation cosines
355 !> @param[out] srot real, optional (npts) clockwise vector rotation sines
356 !> (ugrid=crot*uearth-srot*vearth;
357 !> vgrid=srot*uearth+crot*vearth)
358 !> @param[out] xlon real, optional (npts) dx/dlon in 1/degrees
359 !> @param[out] xlat real, optional (npts) dx/dlat in 1/degrees
360 !> @param[out] ylon real, optional (npts) dy/dlon in 1/degrees
361 !> @param[out] ylat real, optional (npts) dy/dlat in 1/degrees
362 !> @param[out] area real, optional (npts) area weights in m**2
363 !>
364 !> @author Gayno @date 2007-NOV-15
365 SUBROUTINE gdswzd_rot_equid_cylind(self,IOPT,NPTS, &
366 FILL,XPTS,YPTS,RLON,RLAT,NRET, &
367 CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
368 IMPLICIT NONE
369
370 class(ip_rot_equid_cylind_grid), intent(in) :: self
371 INTEGER, INTENT(IN ) :: IOPT, NPTS
372 INTEGER, INTENT( OUT) :: NRET
373 !
374 REAL, INTENT(IN ) :: FILL
375 REAL, INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
376 REAL, INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
377 REAL, OPTIONAL, INTENT( OUT) :: CROT(NPTS),SROT(NPTS)
378 REAL, OPTIONAL, INTENT( OUT) :: XLON(NPTS),XLAT(NPTS)
379 REAL, OPTIONAL, INTENT( OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
380 !
381 INTEGER :: IM,JM,N
382 !
383 LOGICAL :: LROT, LMAP, LAREA
384 !
385 REAL(KIND=kd) :: hs
386 REAL(KIND=kd) :: clonr,clatr,slatr
387 REAL(KIND=kd) :: clat,slat,clon
388 REAL(KIND=kd) :: rlatr,rlonr
389 REAL(KIND=kd) :: wbd,sbd
390 REAL :: XMIN,XMAX,YMIN,YMAX
391 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
392 IF(PRESENT(crot)) crot=fill
393 IF(PRESENT(srot)) srot=fill
394 IF(PRESENT(xlon)) xlon=fill
395 IF(PRESENT(xlat)) xlat=fill
396 IF(PRESENT(ylon)) ylon=fill
397 IF(PRESENT(ylat)) ylat=fill
398 IF(PRESENT(area)) area=fill
399 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
400 ! IS THE EARTH RADIUS DEFINED?
401 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
402 ! IS THIS AN "E"-STAGGER GRID? ROUTINE CAN'T PROCESS THOSE.
403 ! I_OFFSET_ODD=MOD(IGDTMPL(19)/8,2)
404 ! I_OFFSET_EVEN=MOD(IGDTMPL(19)/4,2)
405 ! J_OFFSET=MOD(IGDTMPL(19)/2,2)
406 ! IF(I_OFFSET_ODD/=I_OFFSET_EVEN) THEN
407 ! CALL ROT_EQUID_CYLIND_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
408 ! RETURN
409 ! ENDIF
410 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
411
412
413 rlon0=self%rlon0
414 irot=self%irot
415
416 im=self%im
417 jm=self%jm
418
419 slat0=self%slat0
420 clat0=self%clat0
421
422 wbd=self%wbd
423 sbd=self%sbd
424
425 dlats=self%dlats
426 dlons=self%dlons
427
428 xmin=0
429 xmax=im+1
430 ymin=0
431 ymax=jm+1
432 nret=0
433
434 rerth = self%rerth
435 IF(rerth<0.)THEN
436 CALL rot_equid_cylind_error(iopt,fill,rlat,rlon,xpts,ypts,npts)
437 RETURN
438 ENDIF
439
440 IF(PRESENT(crot).AND.PRESENT(srot))THEN
441 lrot=.true.
442 ELSE
443 lrot=.false.
444 ENDIF
445 IF(PRESENT(xlon).AND.PRESENT(xlat).AND.PRESENT(ylon).AND.PRESENT(ylat))THEN
446 lmap=.true.
447 ELSE
448 lmap=.false.
449 ENDIF
450 IF(PRESENT(area))THEN
451 larea=.true.
452 ELSE
453 larea=.false.
454 ENDIF
455 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
456 ! TRANSLATE GRID COORDINATES TO EARTH COORDINATES
457 IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
458 !$OMP PARALLEL DO PRIVATE(N,RLONR,RLATR,HS,CLONR,SLATR,CLATR,SLAT,CLAT,CLON) &
459 !$OMP& REDUCTION(+:NRET) SCHEDULE(STATIC)
460 DO n=1,npts
461 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
462 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
463 rlonr=wbd+(xpts(n)-1._kd)*dlons
464 rlatr=sbd+(ypts(n)-1._kd)*dlats
465 IF(rlonr <= 0._kd) THEN
466 hs=-1.0_kd
467 ELSE
468 hs=1.0_kd
469 ENDIF
470 clonr=cos(rlonr/dpr)
471 slatr=sin(rlatr/dpr)
472 clatr=cos(rlatr/dpr)
473 slat=clat0*slatr+slat0*clatr*clonr
474 IF(slat.LE.-1) THEN
475 clat=0.
476 clon=cos(rlon0/dpr)
477 rlon(n)=0.
478 rlat(n)=-90.
479 ELSEIF(slat.GE.1) THEN
480 clat=0.
481 clon=cos(rlon0/dpr)
482 rlon(n)=0.
483 rlat(n)=90.
484 ELSE
485 clat=sqrt(1-slat**2)
486 clon=(clat0*clatr*clonr-slat0*slatr)/clat
487 clon=min(max(clon,-1._kd),1._kd)
488 rlon(n)=real(mod(rlon0+hs*dpr*acos(clon)+3600,360._kd))
489 rlat(n)=real(dpr*asin(slat))
490 ENDIF
491 nret=nret+1
492 IF(lrot) CALL rot_equid_cylind_vect_rot(rlon(n), clatr, slatr, &
493 clat, slat, clon, crot(n), srot(n))
494 IF(lmap) CALL rot_equid_cylind_map_jacob(fill, rlon(n), clatr, &
495 clat, slat, clon, xlon(n), xlat(n), ylon(n), ylat(n))
496 IF(larea) CALL rot_equid_cylind_grid_area(clatr, fill, area(n))
497 ELSE
498 rlon(n)=fill
499 rlat(n)=fill
500 ENDIF
501 ENDDO
502 !$OMP END PARALLEL DO
503 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
504 ! TRANSLATE EARTH COORDINATES TO GRID COORDINATES
505 ELSEIF(iopt.EQ.-1) THEN
506 !$OMP PARALLEL DO PRIVATE(N,HS,CLON,SLAT,CLAT,SLATR,CLATR,CLONR,RLONR,RLATR) &
507 !$OMP& REDUCTION(+:NRET) SCHEDULE(STATIC)
508 DO n=1,npts
509 IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90) THEN
510 hs=sign(1._kd,mod(rlon(n)-rlon0+180+3600,360._kd)-180)
511 clon=cos((rlon(n)-rlon0)/dpr)
512 slat=sin(rlat(n)/dpr)
513 clat=cos(rlat(n)/dpr)
514 slatr=clat0*slat-slat0*clat*clon
515 IF(slatr.LE.-1) THEN
516 clatr=0._kd
517 rlonr=0.
518 rlatr=-90.
519 ELSEIF(slatr.GE.1) THEN
520 clatr=0._kd
521 rlonr=0.
522 rlatr=90.
523 ELSE
524 clatr=sqrt(1-slatr**2)
525 clonr=(clat0*clat*clon+slat0*slat)/clatr
526 clonr=min(max(clonr,-1._kd),1._kd)
527 rlonr=hs*dpr*acos(clonr)
528 rlatr=dpr*asin(slatr)
529 ENDIF
530 xpts(n)=real((rlonr-wbd)/dlons+1._kd)
531 ypts(n)=real((rlatr-sbd)/dlats+1._kd)
532 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
533 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
534 nret=nret+1
535 IF(lrot) CALL rot_equid_cylind_vect_rot(rlon(n), clatr, slatr, &
536 clat, slat, clon, crot(n), srot(n))
537 IF(lmap) CALL rot_equid_cylind_map_jacob(fill, rlon(n), clatr, &
538 clat, slat, clon, xlon(n), xlat(n), ylon(n), ylat(n))
539 IF(larea) CALL rot_equid_cylind_grid_area(clatr, fill, area(n))
540 ELSE
541 xpts(n)=fill
542 ypts(n)=fill
543 ENDIF
544 ELSE
545 xpts(n)=fill
546 ypts(n)=fill
547 ENDIF
548 ENDDO
549 !$OMP END PARALLEL DO
550 ENDIF
551 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
552 END SUBROUTINE gdswzd_rot_equid_cylind
553
554 !> Error handler.
555 !>
556 !> Upon an error, this subprogram assigns a "fill" value to the
557 !> output fields.
558 !>
559 !> @param[in] iopt integer option flag
560 !> - +1 to compute earth coords of selected grid coords
561 !> - -1 to compute grid coords of selected earth coords
562 !> @param[in] fill real fill value to set invalid output data
563 !> (must be impossible value; suggested value: -9999.)
564 !> @param[out] rlat real (npts) earth latitudes in degrees n if iopt<0
565 !> @param[out] rlon real (npts) earth longitudes in degrees e if iopt<0
566 !> @param[out] xpts real (npts) grid x point coordinates if iopt>0
567 !> @param[out] ypts real (npts) grid y point coordinates if iopt>0
568 !> @param[in] npts integer maximum number of coordinates
569 !>
570 !> @author Gayno @date 2015-07-13
571 SUBROUTINE rot_equid_cylind_error(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
572 IMPLICIT NONE
573 !
574 INTEGER, INTENT(IN ) :: IOPT, NPTS
575 !
576 REAL, INTENT(IN ) :: FILL
577 REAL, INTENT( OUT) :: RLAT(NPTS),RLON(NPTS)
578 REAL, INTENT( OUT) :: XPTS(NPTS),YPTS(NPTS)
579 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
580 IF(iopt>=0) THEN
581 rlon=fill
582 rlat=fill
583 ENDIF
584 IF(iopt<=0) THEN
585 xpts=fill
586 ypts=fill
587 ENDIF
588 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
589 END SUBROUTINE rot_equid_cylind_error
590
591 !> Vector rotation fields for rotated equidistant cylindrical grids -
592 !> non "e" stagger.
593 !>
594 !> This subprogram computes the vector rotation sines and cosines
595 !> for a rotated equidistant cylindrical grid - non "e" stagger.
596 !>
597 !> ### Program History Log
598 !> Date | Programmer | Comments
599 !> -----|------------|---------
600 !> 2015-01-21 | gayno | initial version
601 !> 2015-07-19 | gayno | rename as "rot_equid_cylind_vect_rot."
602 !> 2018-07-20 | wesley | pass in clatr, slatr, clat, slat, clon for threading.
603 !>
604 !> @param[in] rlon longitude in degrees (real)
605 !> @param[in] clatr cosine of rotated latitude (real)
606 !> @param[in] slatr sine of rotated latitude (real)
607 !> @param[in] clat cosine of latitude (real)
608 !> @param[in] slat sine of latitude (real)
609 !> @param[in] clon cosine of longitude (real)
610 !> @param[out] crot clockwise vector rotation cosines (real)
611 !> @param[out] srot clockwise vector rotation sines (real)
612 !> (ugrid=crot*uearth-srot*vearth;
613 !> vgrid=srot*uearth+crot*vearth)
614 !>
615 !> @author Gayno @date 2015-01-21
616 SUBROUTINE rot_equid_cylind_vect_rot(RLON, CLATR, SLATR, CLAT, SLAT, &
617 CLON, CROT, SROT)
618 IMPLICIT NONE
619
620 REAL(KIND=kd), INTENT(IN ) :: clat, clatr, clon, slat, slatr
621 REAL , INTENT(IN ) :: RLON
622 REAL , INTENT( OUT) :: CROT, SROT
623
624 REAL(KIND=kd) :: slon
625
626 IF(irot.EQ.1) THEN
627 IF(clatr.LE.0._kd) THEN
628 crot=real(-sign(1._kd,slatr*slat0))
629 srot=0.
630 ELSE
631 slon=sin((rlon-rlon0)/dpr)
632 crot=real((clat0*clat+slat0*slat*clon)/clatr)
633 srot=real(slat0*slon/clatr)
634 ENDIF
635 ELSE
636 crot=1.
637 srot=0.
638 ENDIF
639
640 END SUBROUTINE rot_equid_cylind_vect_rot
641
642 !> Map jacobians for rotated equidistant cylindrical
643 !> grids - non "e" stagger.
644 !>
645 !> This subprogram computes the map jacobians for a rotated
646 !> equidistant cylindrical grid - non "e" stagger.
647 !>
648 !> ### Program History Log
649 !> Date | Programmer | Comments
650 !> -----|------------|---------
651 !> 2015-01-21 | gayno | initial version
652 !> 2015-09-17 | gayno | rename as "rot_equid_cylind_map_jacob".
653 !> 2018-07-20 | wesley | pass in clatr, clat, slat, clon to allow threading.
654 !>
655 !> @param[in] fill fill value for undefined points (real)
656 !> @param[in] rlon longitude in degrees (real)
657 !> @param[in] clatr cosine of unrotated latitude (real)
658 !> @param[in] clat cosine of latitude (real)
659 !> @param[in] slat sine of latitude (real)
660 !> @param[in] clon cosine of latitude (real)
661 !> @param[out] xlon dx/dlon in 1/degrees (real)
662 !> @param[out] xlat dx/dlat in 1/degrees (real)
663 !> @param[out] ylon dy/dlon in 1/degrees (real)
664 !> @param[out] ylat dy/dlat in 1/degrees (real)
665 !>
666 !> @author Gayno @date 2015-01-21
667 SUBROUTINE rot_equid_cylind_map_jacob(FILL, RLON, CLATR, CLAT, &
668 SLAT, CLON, XLON, XLAT, YLON, YLAT)
669 IMPLICIT NONE
670
671 REAL(KIND=kd), INTENT(IN ) :: clatr, clat, slat, clon
672 REAL , INTENT(IN ) :: FILL, RLON
673 REAL , INTENT( OUT) :: XLON, XLAT, YLON, YLAT
674
675 REAL(KIND=kd) :: slon, term1, term2
676
677 IF(clatr.LE.0._kd) THEN
678 xlon=fill
679 xlat=fill
680 ylon=fill
681 ylat=fill
682 ELSE
683 slon=sin((rlon-rlon0)/dpr)
684 term1=(clat0*clat+slat0*slat*clon)/clatr
685 term2=slat0*slon/clatr
686 xlon=real(term1*clat/(dlons*clatr))
687 xlat=real(-term2/(dlons*clatr))
688 ylon=real(term2*clat/dlats)
689 ylat=real(term1/dlats)
690 ENDIF
691
692 END SUBROUTINE rot_equid_cylind_map_jacob
693
694 !> Grid box area for rotated equidistant cylindrical grids - non "e"
695 !> stagger.
696 !>
697 !> This subprogram computes the grid box area for a rotated
698 !> equidistant cylindrical grid - non "e" stagger.
699 !>
700 !> ### Program History Log
701 !> Date | Programmer | Comments
702 !> -----|------------|---------
703 !> 2015-01-21 | Gayno | initial version
704 !> 2015-07-19 | gayno | rename as "rot_equid_cylind_grid_area."
705 !> 2018-07-20 | wesley | pass in clatr for threading
706 !>
707 !> @param[in] clatr cosine of unrotated latitude (real)
708 !> @param[in] fill fill value for undefined points (real)
709 !> @param[out] area area weights in m**2 (real)
710 !>
711 !> @author Gayno @date 2015-01-21
712 SUBROUTINE rot_equid_cylind_grid_area(CLATR, FILL, AREA)
713 IMPLICIT NONE
714
715 REAL(KIND=kd), INTENT(IN ) :: clatr
716 REAL, INTENT(IN ) :: FILL
717 REAL, INTENT( OUT) :: AREA
718
719 IF(clatr.LE.0._kd) THEN
720 area=fill
721 ELSE
722 area=real(2._kd*(rerth**2)*clatr*(dlons/dpr)*sin(0.5_kd*dlats/dpr))
723 ENDIF
724
725 END SUBROUTINE rot_equid_cylind_grid_area
726
728
void gdswzd(int igdtnum, int *igdtmpl, int igdtlen, int iopt, int npts, float fill, float *xpts, float *ypts, float *rlon, float *rlat, int *nret, float *crot, float *srot, float *xlon, float *xlat, float *ylon, float *ylat, float *area)
gdswzd() interface for C for _4 build of library.
Determine earth radius and shape.
Module containing common constants.
real, parameter pi
PI.
real, parameter dpr
Radians to degrees.
Users derived type grid descriptor objects to abstract away the raw GRIB1 and GRIB2 grid definitions.
Abstract ip_grid type.
Rotated equidistant cylindrical GRIB decoder and grid coordinate transformations for Arakawa grids A ...
subroutine gdswzd_rot_equid_cylind(self, iopt, npts, fill, xpts, ypts, rlon, rlat, nret, crot, srot, xlon, xlat, ylon, ylat, area)
GDS wizard for rotated equidistant cylindrical.
real(kind=kd) rlon0
Local copy of rlon0.
real(kind=kd) slat0
Local copy of slat0.
subroutine rot_equid_cylind_vect_rot(rlon, clatr, slatr, clat, slat, clon, crot, srot)
Vector rotation fields for rotated equidistant cylindrical grids - non "e" stagger.
subroutine init_grib2_ncep_post(self, g2_desc)
Initializes a Rotated equidistant cylindrical grid given a grib2_descriptor object.
subroutine rot_equid_cylind_error(iopt, fill, rlat, rlon, xpts, ypts, npts)
Error handler.
subroutine rot_equid_cylind_grid_area(clatr, fill, area)
Grid box area for rotated equidistant cylindrical grids - non "e" stagger.
subroutine init_grib2_default(self, g2_desc)
Initializes a Rotated equidistant cylindrical grid given a grib2_descriptor object.
real(kind=kd) dlons
Local copy of dlons.
real(kind=kd) dlats
Local copy of dlats.
subroutine init_grib1(self, g1_desc)
Initializes a Rotated equidistant cylindrical grid given a grib1_descriptor object.
subroutine init_grib2(self, g2_desc)
Initializes a Rotated equidistant cylindrical grid given a grib2_descriptor object.
subroutine rot_equid_cylind_map_jacob(fill, rlon, clatr, clat, slat, clon, xlon, xlat, ylon, ylat)
Map jacobians for rotated equidistant cylindrical grids - non "e" stagger.
real(kind=kd) rerth
Radius of the Earth.
integer, parameter kd
Fortran kind for reals.
real(kind=kd) clat0
Local copy of clat0.
Descriptor representing a grib1 grib descriptor section (GDS) with an integer array.
Abstract grid that holds fields and methods common to all grids.