NCEPLIBS-ip 5.3.0
•All Data Structures Namespaces Files Functions Variables Pages
ip_mercator_grid_mod.F90
Go to the documentation of this file.
1!> @file
2!> @brief GDS wizard for mercator cylindrical.
3!>
4!> @author Iredell @date 96-04-10
5
6!> @brief GDS wizard for mercator cylindrical.
7!>
8!> Octet numbers refer to [GRIB2 - GRID DEFINITION TEMPLATE 3.10 -
9!> Mercator](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-10.shtml).
10!>
11!> @author Iredell @date 96-04-10
12module ip_mercator_grid_mod
13 use ip_grid_descriptor_mod
14 use ip_grid_mod
15 use ip_constants_mod, only: dpr, pi
16 use earth_radius_mod
17 implicit none
18
19 private
20 public :: ip_mercator_grid
21
22 type, extends(ip_grid) :: ip_mercator_grid
23 real :: rlat1 !< Latitude of first grid point. Section 3, octets 39-42.
24 real :: rlon1 !< Longitude of first grid point. Section 3, octets 43-46.
25 real :: rlon2 !< Longitude of last grid point. Section 3, octets 56-59.
26 real :: rlati !< Latitude at which the Mercator projection intersects the Earth. Section 3, octets 48-51.
27 real :: hi !< Scan mode in the 'i' direction. Section 3, octet 60.
28 real :: dlon !< Longitudinal direction grid length. Section 3, octets 65-68.
29 real :: dphi !< Latitudinal direction grid length. Section 3, octets 69-72.
30 contains
31 !> Initializes a gaussian grid given a grib1_descriptor object.
32 !> @return N/A @copydoc ip_mercator_grid_mod::init_grib1
33 procedure :: init_grib1
34 !> Initializes a gaussian grid given a grib2_descriptor object.
35 !> @return N/A @copydoc ip_mercator_grid_mod::init_grib2
36 procedure :: init_grib2
37 !> Calculates Earth coordinates (iopt = 1) or grid coorindates (iopt = -1)
38 !> for Gaussian grids.
39 !> @return N/A @copydoc ip_mercator_grid_mod::gdswzd_mercator
40 procedure :: gdswzd => gdswzd_mercator
41 end type ip_mercator_grid
42
43 REAL :: dlon !< Longitudinal direction grid length.
44 REAL :: dphi !< Latitudinal direction grid length.
45 REAL :: rerth !< Radius of the Earth.
46
47CONTAINS
48
49 !> Initializes a mercator grid given a grib1_descriptor object.
50 !>
51 !> @param[inout] self ip_mercator_grid object.
52 !> @param[in] g1_desc GRIB1 descriptor.
53 !>
54 !> @author Iredell @date 96-04-10
55 subroutine init_grib1(self, g1_desc)
56 class(ip_mercator_grid), intent(inout) :: self
57 type(grib1_descriptor), intent(in) :: g1_desc
58
59 integer :: iscan, jscan
60 real :: dy, hj
61
62 associate(kgds => g1_desc%gds)
63 self%rerth = 6.3712e6
64 self%eccen_squared = 0.0
65
66 self%IM=kgds(2)
67 self%JM=kgds(3)
68
69 self%RLAT1=kgds(4)*1.e-3
70 self%RLON1=kgds(5)*1.e-3
71 self%RLON2=kgds(8)*1.e-3
72 self%RLATI=kgds(9)*1.e-3
73
74 iscan=mod(kgds(11)/128,2)
75 jscan=mod(kgds(11)/64,2)
76
77 dy=kgds(13)
78 self%HI=(-1.)**iscan
79 hj=(-1.)**(1-jscan)
80 self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
81 self%DPHI=hj*dy/(self%RERTH*cos(self%RLATI/dpr))
82
83 ! defaults
84 self%iwrap = 0
85 self%jwrap1 = 0
86 self%jwrap2 = 0
87 self%nscan = mod(kgds(11) / 32, 2)
88 self%nscan_field_pos = self%nscan
89 self%kscan = 0
90
91 self%iwrap = nint(360 / abs(self%dlon))
92 if (self%im < self%iwrap) self%iwrap = 0
93 end associate
94
95 end subroutine init_grib1
96
97 !> Init GRIB2.
98 !>
99 !> @param[inout] self ip_mercator_grid object.
100 !> @param[in] g2_desc GRIB2 descriptor.
101 !>
102 !> @author Iredell @date 96-04-10
103 subroutine init_grib2(self, g2_desc)
104 class(ip_mercator_grid), intent(inout) :: self
105 type(grib2_descriptor), intent(in) :: g2_desc
106
107 integer :: iscan, jscan
108 real :: hj, dy
109
110 associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
111
112 call earth_radius(igdtmpl, igdtlen, self%rerth, self%eccen_squared)
113
114 self%IM=igdtmpl(8)
115 self%JM=igdtmpl(9)
116
117 self%RLAT1=float(igdtmpl(10))*1.0e-6
118 self%RLON1=float(igdtmpl(11))*1.0e-6
119 self%RLON2=float(igdtmpl(15))*1.0e-6
120 self%RLATI=float(igdtmpl(13))*1.0e-6
121
122 iscan=mod(igdtmpl(16)/128,2)
123 jscan=mod(igdtmpl(16)/64,2)
124
125 dy=float(igdtmpl(19))*1.0e-3
126 self%HI=(-1.)**iscan
127 hj=(-1.)**(1-jscan)
128 self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
129 self%DPHI=hj*dy/(self%RERTH*cos(self%RLATI/dpr))
130
131 self%jwrap1 = 0
132 self%jwrap2 = 0
133 self%kscan = 0
134 self%nscan=mod(igdtmpl(16) / 32,2)
135 self%nscan_field_pos = self%nscan
136
137 self%iwrap = nint(360 / abs(self%dlon))
138 if(self%im < self%iwrap) self%iwrap = 0
139
140 end associate
141 end subroutine init_grib2
142
143 !> GDS wizard for mercator cylindrical.
144 !>
145 !> This routine decodes the grib 2 grid definition template (passed
146 !> in integer form as decoded by the ncep g2 library) and returns
147 !> one of the following:
148 !> - (iopt=+1) earth coordinates of selected grid coordinates
149 !> - (iopt=-1) grid coordinates of selected earth coordinates
150 !>
151 !> Works for mercator cylindrical projections.
152 !>
153 !> If the selected coordinates are more than one gridpoint beyond
154 !> the the edges of the grid domain, then the relevant output
155 !> elements are set to fill values.
156 !>
157 !> The actual number of valid points computed is returned too.
158 !>
159 !> Optionally, the vector rotations, map jacobians and the grid box
160 !> areas may be returned. To compute the vector rotations, the
161 !> optional arguments 'srot' and 'crot' must be present. To compute
162 !> the map jacobians, the optional arguments 'xlon', 'xlat', 'ylon',
163 !> 'ylat' must be present. to compute the grid box areas, the
164 !> optional argument 'area' must be present.
165 !>
166 !> ### Program History Log
167 !> Date | Programmer | Comments
168 !> -----|------------|---------
169 !> 96-04-10 | iredell | Initial
170 !> 96-10-01 | iredell | protected against unresolvable points
171 !> 97-10-20 | iredell | include map options
172 !> 2015-01-21 | gayno | merger of gdswiz01() and gdswzd01(). 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.
173 !> 2015-07-13 | gayno | convert to grib 2. replace grib 1 kgds array with grib 2 grid definition template array. Rename.
174 !> 2018-07-20 | wesley | add threads.
175 !>
176 !> @param[in] self grid descriptor.
177 !> @param[in] iopt option flag
178 !> - 1 to compute earth coords of selected grid coords
179 !> - -1 to compute grid coords of selected earth coords
180 !> @param[in] npts maximum number of coordinates
181 !> @param[in] fill fill value to set invalid output data (must be
182 !> impossible value; suggested value: -9999.)
183 !> @param[inout] xpts (npts) grid x point coordinates if iopt>0
184 !> @param[inout] ypts (npts) grid y point coordinates if iopt>0
185 !> @param[inout] rlon (npts) earth longitudes in degrees e if iopt<0
186 !> (acceptable range: -360. to 360.)
187 !> @param[inout] rlat (npts) earth latitudes in degrees n if iopt<0
188 !> (acceptable range: -90. to 90.)
189 !> @param[out] nret number of valid points computed
190 !> @param[out] crot optional (npts) clockwise vector rotation cosines
191 !> @param[out] srot optional (npts) clockwise vector rotation sines
192 !> (ugrid=crot*uearth-srot*vearth; vgrid=srot*uearth+crot*vearth)
193 !> @param[out] xlon optional (npts) dx/dlon in 1/degrees
194 !> @param[out] xlat optional (npts) dx/dlat in 1/degrees
195 !> @param[out] ylon optional (npts) dy/dlon in 1/degrees
196 !> @param[out] ylat optional (npts) dy/dlat in 1/degrees
197 !> @param[out] area optional (npts) area weights in m**2
198 !> (proportional to the square of the map factor)
199 !>
200 !> @author Iredell @date 96-04-10
201 SUBROUTINE gdswzd_mercator(self,IOPT,NPTS,FILL, &
202 XPTS,YPTS,RLON,RLAT,NRET, &
203 CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
204 IMPLICIT NONE
205 !
206 class(ip_mercator_grid), intent(in) :: self
207 INTEGER, INTENT(IN ) :: IOPT, NPTS
208 INTEGER, INTENT( OUT) :: NRET
209 !
210 REAL, INTENT(IN ) :: FILL
211 REAL, INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
212 REAL, INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
213 REAL, OPTIONAL, INTENT( OUT) :: CROT(NPTS),SROT(NPTS)
214 REAL, OPTIONAL, INTENT( OUT) :: XLON(NPTS),XLAT(NPTS)
215 REAL, OPTIONAL, INTENT( OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
216 !
217 INTEGER :: IM, JM, N
218 !
219 LOGICAL :: LROT, LMAP, LAREA
220 !
221 REAL :: HI
222 REAL :: RLAT1, RLON1, RLON2, RLATI
223 REAL :: XMAX, XMIN, YMAX, YMIN
224 REAL :: YE
225 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
226 IF(PRESENT(crot)) crot=fill
227 IF(PRESENT(srot)) srot=fill
228 IF(PRESENT(xlon)) xlon=fill
229 IF(PRESENT(xlat)) xlat=fill
230 IF(PRESENT(ylon)) ylon=fill
231 IF(PRESENT(ylat)) ylat=fill
232 IF(PRESENT(area)) area=fill
233 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234
235 im=self%im
236 jm=self%jm
237
238 rlat1=self%rlat1
239 rlon1=self%rlon1
240 rlon2=self%rlon2
241 rlati=self%rlati
242
243 hi=self%hi
244
245 dlon=self%dlon
246 dphi=self%dphi
247 rerth = self%rerth
248
249 ye=1-log(tan((rlat1+90)/2/dpr))/dphi
250 xmin=0
251 xmax=im+1
252 IF(im.EQ.nint(360/abs(dlon))) xmax=im+2
253 ymin=0
254 ymax=jm+1
255 nret=0
256 IF(PRESENT(crot).AND.PRESENT(srot))THEN
257 lrot=.true.
258 ELSE
259 lrot=.false.
260 ENDIF
261 IF(PRESENT(xlon).AND.PRESENT(xlat).AND.PRESENT(ylon).AND.PRESENT(ylat))THEN
262 lmap=.true.
263 ELSE
264 lmap=.false.
265 ENDIF
266 IF(PRESENT(area))THEN
267 larea=.true.
268 ELSE
269 larea=.false.
270 ENDIF
271 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272 ! TRANSLATE GRID COORDINATES TO EARTH COORDINATES
273 IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
274 !$OMP PARALLEL DO PRIVATE(N) REDUCTION(+:NRET) SCHEDULE(STATIC)
275 DO n=1,npts
276 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
277 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
278 rlon(n)=mod(rlon1+dlon*(xpts(n)-1)+3600,360.)
279 rlat(n)=2*atan(exp(dphi*(ypts(n)-ye)))*dpr-90
280 nret=nret+1
281 IF(lrot) CALL mercator_vect_rot(crot(n),srot(n))
282 IF(lmap) CALL mercator_map_jacob(rlat(n),xlon(n),xlat(n),ylon(n),ylat(n))
283 IF(larea) CALL mercator_grid_area(rlat(n),area(n))
284 ELSE
285 rlon(n)=fill
286 rlat(n)=fill
287 ENDIF
288 ENDDO
289 !$OMP END PARALLEL DO
290 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291 ! TRANSLATE EARTH COORDINATES TO GRID COORDINATES
292 ELSEIF(iopt.EQ.-1) THEN
293 !$OMP PARALLEL DO PRIVATE(N) REDUCTION(+:NRET) SCHEDULE(STATIC)
294 DO n=1,npts
295 IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LT.90) THEN
296 xpts(n)=1+hi*mod(hi*(rlon(n)-rlon1)+3600,360.)/dlon
297 ypts(n)=ye+log(tan((rlat(n)+90)/2/dpr))/dphi
298 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
299 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
300 nret=nret+1
301 IF(lrot) CALL mercator_vect_rot(crot(n),srot(n))
302 IF(lmap) CALL mercator_map_jacob(rlat(n),xlon(n),xlat(n),ylon(n),ylat(n))
303 IF(larea) CALL mercator_grid_area(rlat(n),area(n))
304 ELSE
305 xpts(n)=fill
306 ypts(n)=fill
307 ENDIF
308 ELSE
309 xpts(n)=fill
310 ypts(n)=fill
311 ENDIF
312 ENDDO
313 !$OMP END PARALLEL DO
314 ENDIF
315 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
316 END SUBROUTINE gdswzd_mercator
317
318 !> Vector rotation fields for mercator cylindrical grids.
319 !>
320 !> This subprogram computes the vector rotation sines and cosines
321 !> for a mercator cylindrical grid.
322 !>
323 !> ### Program History Log
324 !> Date | Programmer | Comments
325 !> -----|------------|---------
326 !> 2015-01-21 | gayno | initial version
327 !> 2015-09-17 | gayno | rename as "mercator_vect_rot".
328 !>
329 !> @param[in] crot clockwise vector rotation cosines (real)
330 !> @param[in] srot clockwise vector rotation sines (real)
331 !> (ugrid=crot*uearth-srot*vearth; vgrid=srot*uearth+crot*vearth)
332 !>
333 !> @author Gayno @date 2015-01-21
334 SUBROUTINE mercator_vect_rot(CROT,SROT)
335 IMPLICIT NONE
336
337 REAL, INTENT( OUT) :: CROT, SROT
338
339 crot=1.0
340 srot=0.0
341
342 END SUBROUTINE mercator_vect_rot
343
344 !> Map jacobians for mercator cylindrical grids.
345 !>
346 !> This subprogram computes the map jacobians for a mercator
347 !> cylindrical grid.
348 !>
349 !> ### Program History Log
350 !> Date | Programmer | Comments
351 !> -----|------------|---------
352 !> 2015-01-21 | gayno | initial version
353 !> 2015-09-17 | gayno | rename as "mercator_map_jacob"
354 !>
355 !> @param[in] rlat latitude in degrees (real)
356 !> @param[out] xlon dx/dlon in 1/degrees (real)
357 !> @param[out] xlat dx/dlat in 1/degrees (real)
358 !> @param[out] ylon dy/dlon in 1/degrees (real)
359 !> @param[out] ylat dy/dlat in 1/degrees (real)
360 !>
361 !> @author Gayno @date 2015-01-21
362 SUBROUTINE mercator_map_jacob(RLAT,XLON,XLAT,YLON,YLAT)
363 IMPLICIT NONE
364
365 REAL, INTENT(IN ) :: RLAT
366 REAL, INTENT( OUT) :: XLON, XLAT, YLON, YLAT
367
368 xlon=1./dlon
369 xlat=0.
370 ylon=0.
371 ylat=1./dphi/cos(rlat/dpr)/dpr
372
373 END SUBROUTINE mercator_map_jacob
374
375 !> Grid box area for mercator cylindrical grids.
376 !>
377 !> This subprogram computes the grid box area for a mercator
378 !> cylindrical grid.
379 !>
380 !> ### Program History Log
381 !> Date | Programmer | Comments
382 !> -----|------------|---------
383 !> 2015-01-21 | gayno | initial version
384 !> 2015-09-17 | gayno | rename as "mercator_grid_area"
385 !>
386 !> @param[in] rlat latitude of grid point in degrees (real)
387 !> @param[out] area area weights in m**2 (real)
388 !>
389 !> @author Gayno @date 2015-01-21
390 SUBROUTINE mercator_grid_area(RLAT,AREA)
391 IMPLICIT NONE
392
393 REAL, INTENT(IN ) :: RLAT
394 REAL, INTENT( OUT) :: AREA
395
396 area=rerth**2*cos(rlat/dpr)**2*dphi*dlon/dpr
397
398 END SUBROUTINE mercator_grid_area
399
400end module ip_mercator_grid_mod
401
gdswzd
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.
earth_radius_mod
Determine earth radius and shape.
Definition earth_radius_mod.F90:7
ip_constants_mod
Module containing common constants.
Definition ip_constants_mod.F90:9
ip_constants_mod::pi
real, parameter pi
PI.
Definition ip_constants_mod.F90:14
ip_constants_mod::dpr
real, parameter dpr
Radians to degrees.
Definition ip_constants_mod.F90:15
ip_grid_descriptor_mod
Users derived type grid descriptor objects to abstract away the raw GRIB1 and GRIB2 grid definitions.
Definition ip_grid_descriptor_mod.F90:15
ip_grid_mod
Abstract ip_grid type.
Definition ip_grid_mod.F90:10
ip_mercator_grid_mod
GDS wizard for mercator cylindrical.
Definition ip_mercator_grid_mod.F90:12
ip_mercator_grid_mod::dlon
real dlon
Longitudinal direction grid length.
Definition ip_mercator_grid_mod.F90:43
ip_mercator_grid_mod::dphi
real dphi
Latitudinal direction grid length.
Definition ip_mercator_grid_mod.F90:44
ip_mercator_grid_mod::init_grib1
subroutine init_grib1(self, g1_desc)
Initializes a mercator grid given a grib1_descriptor object.
Definition ip_mercator_grid_mod.F90:56
ip_mercator_grid_mod::mercator_grid_area
subroutine mercator_grid_area(rlat, area)
Grid box area for mercator cylindrical grids.
Definition ip_mercator_grid_mod.F90:391
ip_mercator_grid_mod::gdswzd_mercator
subroutine gdswzd_mercator(self, iopt, npts, fill, xpts, ypts, rlon, rlat, nret, crot, srot, xlon, xlat, ylon, ylat, area)
GDS wizard for mercator cylindrical.
Definition ip_mercator_grid_mod.F90:204
ip_mercator_grid_mod::mercator_vect_rot
subroutine mercator_vect_rot(crot, srot)
Vector rotation fields for mercator cylindrical grids.
Definition ip_mercator_grid_mod.F90:335
ip_mercator_grid_mod::init_grib2
subroutine init_grib2(self, g2_desc)
Init GRIB2.
Definition ip_mercator_grid_mod.F90:104
ip_mercator_grid_mod::mercator_map_jacob
subroutine mercator_map_jacob(rlat, xlon, xlat, ylon, ylat)
Map jacobians for mercator cylindrical grids.
Definition ip_mercator_grid_mod.F90:363
ip_mercator_grid_mod::rerth
real rerth
Radius of the Earth.
Definition ip_mercator_grid_mod.F90:45
ip_grid_descriptor_mod::grib1_descriptor
Descriptor representing a grib1 grib descriptor section (GDS) with an integer array.
Definition ip_grid_descriptor_mod.F90:39
ip_grid_mod::ip_grid
Abstract grid that holds fields and methods common to all grids.
Definition ip_grid_mod.F90:59
ip_mercator_grid_mod::ip_mercator_grid
Definition ip_mercator_grid_mod.F90:22