NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
ip_polar_stereo_grid_mod.F90
Go to the documentation of this file.
1
5
13module ip_polar_stereo_grid_mod
14 use ip_grid_descriptor_mod
15 use ip_grid_mod
16 use ip_constants_mod, only: dpr, pi, pi2, pi4, rerth_wgs84, e2_wgs84
17 use earth_radius_mod
18 implicit none
19
20 private
21 public :: ip_polar_stereo_grid
22
23 type, extends(ip_grid) :: ip_polar_stereo_grid
24 logical :: elliptical
25 real :: rlat1
26 real :: rlon1
27 real :: orient
28 real :: h
29 real :: dxs
30 real :: dys
31 real :: slatr
35 integer :: irot
36 contains
39 procedure :: init_grib1
42 procedure :: init_grib2
46 procedure :: gdswzd => gdswzd_polar_stereo
47 end type ip_polar_stereo_grid
48
49 INTEGER :: irot
50 REAL :: de2
51 REAL :: dxs
52 REAL :: dys
53 REAL :: e2
54 REAL :: rerth
55 REAL :: h
56 REAL :: orient
57 REAL :: tinyreal=tiny(1.0)
58
59CONTAINS
60
68 subroutine init_grib1(self, g1_desc)
69 class(ip_polar_stereo_grid), intent(inout) :: self
70 type(grib1_descriptor), intent(in) :: g1_desc
71
72 REAL, PARAMETER :: SLAT=60.0 ! standard latitude according grib1 standard
73
74 real :: dx, dy, hi, hj
75 integer :: iproj, iscan, jscan
76
77 associate(kgds => g1_desc%gds)
78 self%ELLIPTICAL=mod(kgds(6)/64,2).EQ.1
79
80 if (.not. self%elliptical) then
81 self%rerth = 6.3712e6
82 self%eccen_squared = 0d0
83 else
84 self%rerth = rerth_wgs84
85 self%eccen_squared = e2_wgs84 !wgs84 datum
86 end if
87
88 self%IM=kgds(2)
89 self%JM=kgds(3)
90
91 self%RLAT1=kgds(4)*1.e-3
92 self%RLON1=kgds(5)*1.e-3
93
94 self%IROT=mod(kgds(6)/8,2)
95
96 self%SLATR=slat/dpr
97
98 self%ORIENT=kgds(7)*1.e-3
99
100 dx=kgds(8)
101 dy=kgds(9)
102
103 iproj=mod(kgds(10)/128,2)
104 iscan=mod(kgds(11)/128,2)
105 jscan=mod(kgds(11)/64,2)
106
107 self%H=(-1.)**iproj
108 hi=(-1.)**iscan
109 hj=(-1.)**(1-jscan)
110
111 IF(abs(self%H+1.).LT.tinyreal) self%ORIENT=self%ORIENT+180.
112
113 self%DXS=dx*hi
114 self%DYS=dy*hj
115
116 self%iwrap= 0
117 self%jwrap1 = 0
118 self%jwrap2 = 0
119 self%nscan = mod(kgds(11) / 32, 2)
120 self%nscan_field_pos = self%nscan
121 self%kscan = 0
122 end associate
123
124 end subroutine init_grib1
125
133 subroutine init_grib2(self, g2_desc)
134 class(ip_polar_stereo_grid), intent(inout) :: self
135 type(grib2_descriptor), intent(in) :: g2_desc
136
137 real :: slat, dx, dy, hi, hj
138 integer :: iproj, iscan, jscan
139
140 associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
141 call earth_radius(igdtmpl, igdtlen, self%rerth, self%eccen_squared)
142
143 self%ELLIPTICAL = self%eccen_squared > 0.0
144
145 self%IM=igdtmpl(8)
146 self%JM=igdtmpl(9)
147
148 self%RLAT1=float(igdtmpl(10))*1.e-6
149 self%RLON1=float(igdtmpl(11))*1.e-6
150
151 self%IROT=mod(igdtmpl(12)/8,2)
152
153 slat=float(abs(igdtmpl(13)))*1.e-6
154 self%SLATR=slat/dpr
155
156 self%ORIENT=float(igdtmpl(14))*1.e-6
157
158 dx=float(igdtmpl(15))*1.e-3
159 dy=float(igdtmpl(16))*1.e-3
160
161 iproj=mod(igdtmpl(17)/128,2)
162 iscan=mod(igdtmpl(18)/128,2)
163 jscan=mod(igdtmpl(18)/64,2)
164
165 self%H=(-1.)**iproj
166 hi=(-1.)**iscan
167 hj=(-1.)**(1-jscan)
168
169 self%DXS=dx*hi
170 self%DYS=dy*hj
171
172 self%nscan = mod(igdtmpl(18) / 32, 2)
173 self%nscan_field_pos = self%nscan
174 self%iwrap = 0
175 self%jwrap1 = 0
176 self%jwrap2 = 0
177 self%kscan = 0
178 end associate
179 end subroutine init_grib2
180
244 SUBROUTINE gdswzd_polar_stereo(self,IOPT,NPTS, &
245 FILL,XPTS,YPTS,RLON,RLAT,NRET, &
246 CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
247 IMPLICIT NONE
248 !
249
250 class(ip_polar_stereo_grid), intent(in) :: self
251 INTEGER, INTENT(IN ) :: IOPT, NPTS
252 INTEGER, INTENT( OUT) :: NRET
253 !
254 REAL, INTENT(IN ) :: FILL
255 REAL, INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
256 REAL, INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
257 REAL, OPTIONAL, INTENT( OUT) :: CROT(NPTS),SROT(NPTS)
258 REAL, OPTIONAL, INTENT( OUT) :: XLON(NPTS),XLAT(NPTS)
259 REAL, OPTIONAL, INTENT( OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
260 !
261 INTEGER :: IM, JM
262 INTEGER :: ITER, N
263 !
264 LOGICAL :: ELLIPTICAL, LROT, LMAP, LAREA
265 !
266 REAL :: ALAT, ALAT1, ALONG, DIFF
267 REAL :: DI, DJ, DE
268 REAL :: DR, E, E_OVER_2
269 REAL :: MC, SLATR
270 REAL :: RLAT1, RLON1, RHO, T, TC
271 REAL :: XMAX, XMIN, YMAX, YMIN
272 REAL :: XP, YP, DR2
273 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274 IF(PRESENT(crot)) crot=fill
275 IF(PRESENT(srot)) srot=fill
276 IF(PRESENT(xlon)) xlon=fill
277 IF(PRESENT(xlat)) xlat=fill
278 IF(PRESENT(ylon)) ylon=fill
279 IF(PRESENT(ylat)) ylat=fill
280 IF(PRESENT(area)) area=fill
281
282 elliptical = self%elliptical
283 im=self%im
284 jm=self%jm
285
286 rlat1=self%rlat1
287 rlon1=self%rlon1
288
289 irot=self%irot
290 slatr=self%slatr
291 orient=self%orient
292
293 h=self%h
294 dxs=self%dxs
295 dys=self%dys
296
297 rerth = self%rerth
298 e2 = self%eccen_squared
299 !
300 ! FIND X/Y OF POLE
301 IF (.NOT.elliptical) THEN
302 de=(1.+sin(slatr))*rerth
303 dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
304 xp=1-h*sin((rlon1-orient)/dpr)*dr/dxs
305 yp=1+cos((rlon1-orient)/dpr)*dr/dys
306 de2=de**2
307 ELSE
308 e=sqrt(e2)
309 e_over_2=e*0.5
310 alat=h*rlat1/dpr
311 along = (rlon1-orient)/dpr
312 t=tan(pi4-alat/2.)/((1.-e*sin(alat))/ &
313 (1.+e*sin(alat)))**(e_over_2)
314 tc=tan(pi4-slatr/2.)/((1.-e*sin(slatr))/ &
315 (1.+e*sin(slatr)))**(e_over_2)
316 mc=cos(slatr)/sqrt(1.0-e2*(sin(slatr)**2))
317 rho=rerth*mc*t/tc
318 yp = 1.0 + rho*cos(h*along)/dys
319 xp = 1.0 - rho*sin(h*along)/dxs
320 ENDIF ! ELLIPTICAL
321 xmin=0
322 xmax=im+1
323 ymin=0
324 ymax=jm+1
325 nret=0
326 IF(PRESENT(crot).AND.PRESENT(srot))THEN
327 lrot=.true.
328 ELSE
329 lrot=.false.
330 ENDIF
331 IF(PRESENT(xlon).AND.PRESENT(xlat).AND.PRESENT(ylon).AND.PRESENT(ylat))THEN
332 lmap=.true.
333 ELSE
334 lmap=.false.
335 ENDIF
336 IF(PRESENT(area))THEN
337 larea=.true.
338 ELSE
339 larea=.false.
340 ENDIF
341 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
342 ! TRANSLATE GRID COORDINATES TO EARTH COORDINATES
343 IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
344 IF(.NOT.elliptical)THEN
345 !$OMP PARALLEL DO PRIVATE(N,DI,DJ,DR2) REDUCTION(+:NRET) SCHEDULE(STATIC)
346 DO n=1,npts
347 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
348 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
349 di=(xpts(n)-xp)*dxs
350 dj=(ypts(n)-yp)*dys
351 dr2=di**2+dj**2
352 IF(dr2.LT.de2*1.e-6) THEN
353 rlon(n)=0.
354 rlat(n)=h*90.
355 ELSE
356 rlon(n)=mod(orient+h*dpr*atan2(di,-dj)+3600,360.)
357 rlat(n)=h*dpr*asin((de2-dr2)/(de2+dr2))
358 ENDIF
359 nret=nret+1
360 IF(lrot) CALL polar_stereo_vect_rot(rlon(n),crot(n),srot(n))
361 IF(lmap) CALL polar_stereo_map_jacob(rlon(n),rlat(n),dr2, &
362 xlon(n),xlat(n),ylon(n),ylat(n))
363 IF(larea) CALL polar_stereo_grid_area(rlat(n),dr2,area(n))
364 ELSE
365 rlon(n)=fill
366 rlat(n)=fill
367 ENDIF
368 ENDDO
369 !$OMP END PARALLEL DO
370 ELSE ! ELLIPTICAL
371 !$OMP PARALLEL DO PRIVATE(N,DI,DJ,RHO,T,ALONG,ALAT1,ALAT,DIFF) &
372 !$OMP& REDUCTION(+:NRET) SCHEDULE(STATIC)
373 DO n=1,npts
374 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
375 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
376 di=(xpts(n)-xp)*dxs
377 dj=(ypts(n)-yp)*dys
378 rho=sqrt(di*di+dj*dj)
379 t=(rho*tc)/(rerth*mc)
380 IF(abs(ypts(n)-yp)<0.01)THEN
381 IF(di>0.0) along=orient+h*90.0
382 IF(di<=0.0) along=orient-h*90.0
383 ELSE
384 along=orient+h*atan(di/(-dj))*dpr
385 IF(dj>0) along=along+180.
386 END IF
387 alat1=pi2-2.0*atan(t)
388 DO iter=1,10
389 alat = pi2 - 2.0*atan(t*(((1.0-e*sin(alat1))/ &
390 (1.0+e*sin(alat1)))**(e_over_2)))
391 diff = abs(alat-alat1)*dpr
392 IF (diff < 0.000001) EXIT
393 alat1=alat
394 ENDDO
395 rlat(n)=h*alat*dpr
396 rlon(n)=along
397 IF(rlon(n)<0.0) rlon(n)=rlon(n)+360.
398 IF(rlon(n)>360.0) rlon(n)=rlon(n)-360.0
399 nret=nret+1
400 IF(lrot) CALL polar_stereo_vect_rot(rlon(n),crot(n),srot(n))
401 ELSE
402 rlon(n)=fill
403 rlat(n)=fill
404 ENDIF
405 ENDDO
406 !$OMP END PARALLEL DO
407 ENDIF ! ELLIPTICAL
408 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
409 ! TRANSLATE EARTH COORDINATES TO GRID COORDINATES
410 ELSEIF(iopt.EQ.-1) THEN
411 IF(.NOT.elliptical)THEN
412 !$OMP PARALLEL DO PRIVATE(N,DR,DR2) REDUCTION(+:NRET) SCHEDULE(STATIC)
413 DO n=1,npts
414 IF(abs(rlon(n)).LT.(360.+tinyreal).AND.abs(rlat(n)).LT.(90.+tinyreal).AND. &
415 abs(h*rlat(n)+90).GT.tinyreal) THEN
416 dr=de*tan((90-h*rlat(n))/2/dpr)
417 dr2=dr**2
418 xpts(n)=xp+h*sin((rlon(n)-orient)/dpr)*dr/dxs
419 ypts(n)=yp-cos((rlon(n)-orient)/dpr)*dr/dys
420 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
421 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
422 nret=nret+1
423 IF(lrot) CALL polar_stereo_vect_rot(rlon(n),crot(n),srot(n))
424 IF(lmap) CALL polar_stereo_map_jacob(rlon(n),rlat(n),dr2, &
425 xlon(n),xlat(n),ylon(n),ylat(n))
426 IF(larea) CALL polar_stereo_grid_area(rlat(n),dr2,area(n))
427 ELSE
428 xpts(n)=fill
429 ypts(n)=fill
430 ENDIF
431 ELSE
432 xpts(n)=fill
433 ypts(n)=fill
434 ENDIF
435 ENDDO
436 !$OMP END PARALLEL DO
437 ELSE ! ELLIPTICAL CASE
438 !$OMP PARALLEL DO PRIVATE(N,ALAT,ALONG,T,RHO) REDUCTION(+:NRET) SCHEDULE(STATIC)
439 DO n=1,npts
440 IF(abs(rlon(n)).LT.(360+tinyreal).AND.abs(rlat(n)).LT.(90+tinyreal).AND. &
441 abs(h*rlat(n)+90).GT.tinyreal) THEN
442 alat = h*rlat(n)/dpr
443 along = (rlon(n)-orient)/dpr
444 t=tan(pi4-alat*0.5)/((1.-e*sin(alat))/ &
445 (1.+e*sin(alat)))**(e_over_2)
446 rho=rerth*mc*t/tc
447 xpts(n)= xp + rho*sin(h*along) / dxs
448 ypts(n)= yp - rho*cos(h*along) / dys
449 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
450 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
451 nret=nret+1
452 IF(lrot) CALL polar_stereo_vect_rot(rlon(n),crot(n),srot(n))
453 ELSE
454 xpts(n)=fill
455 ypts(n)=fill
456 ENDIF
457 ELSE
458 xpts(n)=fill
459 ypts(n)=fill
460 ENDIF
461 ENDDO
462 !$OMP END PARALLEL DO
463 ENDIF
464 ENDIF
465 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
466 END SUBROUTINE gdswzd_polar_stereo
467
487 SUBROUTINE polar_stereo_vect_rot(RLON, CROT, SROT)
488 IMPLICIT NONE
489
490 REAL, INTENT(IN ) :: RLON
491 REAL, INTENT( OUT) :: CROT, SROT
492
493 IF(irot.EQ.1) THEN
494 crot=h*cos((rlon-orient)/dpr)
495 srot=sin((rlon-orient)/dpr)
496 ELSE
497 crot=1.
498 srot=0.
499 ENDIF
500
501 END SUBROUTINE polar_stereo_vect_rot
502
525 SUBROUTINE polar_stereo_map_jacob(RLON,RLAT,DR2,XLON,XLAT,YLON,YLAT)
526 IMPLICIT NONE
527
528 REAL, INTENT(IN ) :: RLON, RLAT, DR2
529 REAL, INTENT( OUT) :: XLON, XLAT, YLON, YLAT
530
531 REAL :: CLAT, DE, DR
532
533 IF(dr2.LT.de2*1.e-6) THEN
534 de=sqrt(de2)
535 xlon=0.
536 xlat=-sin((rlon-orient)/dpr)/dpr*de/dxs/2
537 ylon=0.
538 ylat=h*cos((rlon-orient)/dpr)/dpr*de/dys/2
539 ELSE
540 dr=sqrt(dr2)
541 clat=cos(rlat/dpr)
542 xlon=h*cos((rlon-orient)/dpr)/dpr*dr/dxs
543 xlat=-sin((rlon-orient)/dpr)/dpr*dr/dxs/clat
544 ylon=sin((rlon-orient)/dpr)/dpr*dr/dys
545 ylat=h*cos((rlon-orient)/dpr)/dpr*dr/dys/clat
546 ENDIF
547
548 END SUBROUTINE polar_stereo_map_jacob
549
568 SUBROUTINE polar_stereo_grid_area(RLAT, DR2, AREA)
569 IMPLICIT NONE
570
571 REAL, INTENT(IN ) :: RLAT, DR2
572 REAL, INTENT( OUT) :: AREA
573
574 REAL :: CLAT
575
576 IF(dr2.LT.de2*1.e-6) THEN
577 area=rerth**2*abs(dxs)*abs(dys)*4/de2
578 ELSE
579 clat=cos(rlat/dpr)
580 area=rerth**2*clat**2*abs(dxs)*abs(dys)/dr2
581 ENDIF
582
583 END SUBROUTINE polar_stereo_grid_area
584
585end module ip_polar_stereo_grid_mod
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::pi2
real, parameter pi2
PI / 2.0.
Definition ip_constants_mod.F90:16
ip_constants_mod::pi
real, parameter pi
PI.
Definition ip_constants_mod.F90:14
ip_constants_mod::pi4
real, parameter pi4
PI / 4.0.
Definition ip_constants_mod.F90:17
ip_constants_mod::rerth_wgs84
real, parameter rerth_wgs84
Radius of the Earth defined by WGS-84.
Definition ip_constants_mod.F90:18
ip_constants_mod::dpr
real, parameter dpr
Radians to degrees.
Definition ip_constants_mod.F90:15
ip_constants_mod::e2_wgs84
real, parameter e2_wgs84
Eccentricity squared of Earth defined by WGS-84.
Definition ip_constants_mod.F90:19
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_polar_stereo_grid_mod
GDS wizard for polar stereographic azimuthal.
Definition ip_polar_stereo_grid_mod.F90:13
ip_polar_stereo_grid_mod::irot
integer irot
Local copy of irot.
Definition ip_polar_stereo_grid_mod.F90:49
ip_polar_stereo_grid_mod::rerth
real rerth
Radius of the Earth.
Definition ip_polar_stereo_grid_mod.F90:54
ip_polar_stereo_grid_mod::init_grib2
subroutine init_grib2(self, g2_desc)
Initializes a polar stereographic grid given a grib2_descriptor object.
Definition ip_polar_stereo_grid_mod.F90:134
ip_polar_stereo_grid_mod::h
real h
Local copy of h.
Definition ip_polar_stereo_grid_mod.F90:55
ip_polar_stereo_grid_mod::de2
real de2
Square of DE.
Definition ip_polar_stereo_grid_mod.F90:50
ip_polar_stereo_grid_mod::polar_stereo_vect_rot
subroutine polar_stereo_vect_rot(rlon, crot, srot)
Vector rotation fields for polar stereographic grids.
Definition ip_polar_stereo_grid_mod.F90:488
ip_polar_stereo_grid_mod::init_grib1
subroutine init_grib1(self, g1_desc)
Initializes a polar stereographic grid given a grib1_descriptor object.
Definition ip_polar_stereo_grid_mod.F90:69
ip_polar_stereo_grid_mod::polar_stereo_map_jacob
subroutine polar_stereo_map_jacob(rlon, rlat, dr2, xlon, xlat, ylon, ylat)
Map jacobians for polar stereographic grids.
Definition ip_polar_stereo_grid_mod.F90:526
ip_polar_stereo_grid_mod::orient
real orient
Local copy of orient.
Definition ip_polar_stereo_grid_mod.F90:56
ip_polar_stereo_grid_mod::e2
real e2
Eccentricity squared.
Definition ip_polar_stereo_grid_mod.F90:53
ip_polar_stereo_grid_mod::dxs
real dxs
Local copy of dxs.
Definition ip_polar_stereo_grid_mod.F90:51
ip_polar_stereo_grid_mod::polar_stereo_grid_area
subroutine polar_stereo_grid_area(rlat, dr2, area)
Grid box area for polar stereographic grids.
Definition ip_polar_stereo_grid_mod.F90:569
ip_polar_stereo_grid_mod::gdswzd_polar_stereo
subroutine gdswzd_polar_stereo(self, iopt, npts, fill, xpts, ypts, rlon, rlat, nret, crot, srot, xlon, xlat, ylon, ylat, area)
GDS wizard for polar stereographic azimuthal.
Definition ip_polar_stereo_grid_mod.F90:247
ip_polar_stereo_grid_mod::dys
real dys
Local copy of dys.
Definition ip_polar_stereo_grid_mod.F90:52
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:58
ip_polar_stereo_grid_mod::ip_polar_stereo_grid
Definition ip_polar_stereo_grid_mod.F90:23