ip_polar_stereo_grid_mod.f90

module ip_polar_stereo_grid_mod
  use ip_grid_descriptor_mod
  use ip_grid_mod
  use constants_mod, only: dpr, pi, pi2, pi4, rerth_wgs84, e2_wgs84
  use earth_radius_mod
  implicit none
 
  private
  public :: ip_polar_stereo_grid
 
  type, extends(ip_grid) :: ip_polar_stereo_grid
     logical :: elliptical
     real :: rlat1, rlon1, orient, h, dxs, dys, slatr
     integer :: irot
   contains
     procedure :: init_grib1
     procedure :: init_grib2
     procedure :: gdswzd => gdswzd_polar_stereo
  end type ip_polar_stereo_grid
 
  INTEGER                         :: IROT
 
  REAL                            :: DE2, DXS, DYS
  REAL                            :: E2, RERTH , H, ORIENT
 
CONTAINS
 
  subroutine init_grib1(self, g1_desc)
    class(ip_polar_stereo_grid), intent(inout) :: self
    type(grib1_descriptor), intent(in) :: g1_desc
 
    REAL, PARAMETER :: SLAT=60.0  ! standard latitude according grib1 standard
 
    real :: dx, dy, hi, hj
    integer :: iproj, iscan, jscan
 
    associate(kgds => g1_desc%gds)
      self%ELLIPTICAL=mod(kgds(6)/64,2).EQ.1
      
      if (.not. self%elliptical) then
         self%rerth = 6.3712e6
         self%eccen_squared = 0d0
      else
         self%rerth = rerth_wgs84
         self%eccen_squared = e2_wgs84 !wgs84 datum
      end if
 
      self%IM=kgds(2)
      self%JM=kgds(3)
      
      self%RLAT1=kgds(4)*1.e-3
      self%RLON1=kgds(5)*1.e-3
 
      self%IROT=mod(kgds(6)/8,2)
 
      self%SLATR=slat/dpr
 
      self%ORIENT=kgds(7)*1.e-3
 
      dx=kgds(8)
      dy=kgds(9)
 
      iproj=mod(kgds(10)/128,2)
      iscan=mod(kgds(11)/128,2)
      jscan=mod(kgds(11)/64,2)
 
      self%H=(-1.)**iproj
      hi=(-1.)**iscan
      hj=(-1.)**(1-jscan)
 
      IF(self%H.EQ.-1)self%ORIENT=self%ORIENT+180.
 
      self%DXS=dx*hi
      self%DYS=dy*hj
 
      self%iwrap= 0
      self%jwrap1 = 0
      self%jwrap2 = 0
      self%nscan = mod(kgds(11) / 32, 2)
      self%nscan_field_pos = self%nscan
      self%kscan = 0
    end associate
 
  end subroutine init_grib1
 
  subroutine init_grib2(self, g2_desc)
    class(ip_polar_stereo_grid), intent(inout) :: self
    type(grib2_descriptor), intent(in) :: g2_desc
 
    real :: slat, dx, dy, hi, hj
    integer :: iproj, iscan, jscan
 
    associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
      call earth_radius(igdtmpl, igdtlen, self%rerth, self%eccen_squared)
 
      self%ELLIPTICAL = self%eccen_squared > 0.0
 
      self%IM=igdtmpl(8)
      self%JM=igdtmpl(9)
 
      self%RLAT1=float(igdtmpl(10))*1.e-6
      self%RLON1=float(igdtmpl(11))*1.e-6
 
      self%IROT=mod(igdtmpl(12)/8,2)
 
      slat=float(abs(igdtmpl(13)))*1.e-6
      self%SLATR=slat/dpr
 
      self%ORIENT=float(igdtmpl(14))*1.e-6
 
      dx=float(igdtmpl(15))*1.e-3
      dy=float(igdtmpl(16))*1.e-3
 
      iproj=mod(igdtmpl(17)/128,2)
      iscan=mod(igdtmpl(18)/128,2)
      jscan=mod(igdtmpl(18)/64,2)
 
      self%H=(-1.)**iproj
      hi=(-1.)**iscan
      hj=(-1.)**(1-jscan)
 
      self%DXS=dx*hi
      self%DYS=dy*hj
 
      self%nscan = mod(igdtmpl(18) / 32, 2)
      self%nscan_field_pos = self%nscan
      self%iwrap = 0
      self%jwrap1 = 0
      self%jwrap2 = 0
      self%kscan = 0
    end associate
  end subroutine init_grib2
 
 
  SUBROUTINE gdswzd_polar_stereo(self,IOPT,NPTS, &
       FILL,XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  GDSWZD_POLAR_STEREO   GDS WIZARD FOR POLAR STEREOGRAPHIC 
    !                                    AZIMUTHAL
    !   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
    !
    ! ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB 2 GRID DEFINITION
    !           TEMPLATE (PASSED IN INTEGER FORM AS DECODED BY THE
    !           NCEP G2 LIBRARY) AND RETURNS ONE OF THE FOLLOWING:
    !             (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES
    !             (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES
    !           WORKS FOR POLAR STEREOGRAPHIC AZIMUTHAL PROJECTIONS.
    !           IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT
    !           BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT
    !           OUTPUT ELEMENTS ARE SET TO FILL VALUES.
    !           THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO.
    !           OPTIONALLY, THE VECTOR ROTATIONS, MAP JACOBIANS, AND
    !           GRID BOX AREAS MAY BE RETURNED AS WELL.  ROUTINE WORKS
    !           FOR BOTH SPHERICAL AND ELLIPTICAL EARTHS WITH THE
    !           EXCEPTION OF THE MAP JACOBIANS AND GRID BOX AREAS, WHICH
    !           ARE ONLY COMPUTED FOR SPHERICAL EARTHS.  TO COMPUTE
    !           THE VECTOR ROTATIONS, THE OPTIONAL ARGUMENTS 'SROT' AND 'CROT'
    !           MUST BE PRESENT.  TO COMPUTE THE MAP JACOBIANS, THE
    !           OPTIONAL ARGUMENTS 'XLON', 'XLAT', 'YLON', 'YLAT' MUST BE PRESENT.
    !           TO COMPUTE THE GRID BOX AREAS, THE OPTIONAL ARGUMENT
    !           'AREA' MUST BE PRESENT.
    !
    ! PROGRAM HISTORY LOG:
    !   96-04-10  IREDELL
    !   97-10-20  IREDELL  INCLUDE MAP OPTIONS
    !   09-05-13  GAYNO    ENSURE AREA ALWAYS POSITIVE
    ! 2015-01-21  GAYNO    MERGER OF GDSWIZ05 AND GDSWZD05.  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.
    !                      INCLUDE OPTION FOR ELLIPTICAL EARTHS.
    ! 2015-07-13  GAYNO    CONVERT TO GRIB 2. REPLACE GRIB 1 KGDS ARRAY
    !                      WITH GRIB 2 GRID DEFINITION TEMPLATE ARRAY.
    !                      RENAME ROUTINE AS "GDSWZD_POLAR_STEREO".
    ! 2018-07-20  WESLEY   ADD THREADING.
    !
    ! USAGE:   CALL GDSWZD_POLAR_STEREO(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,
    !    &                              FILL,XPTS,YPTS,RLON,RLAT,NRET,
    !    &                              CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    !
    !   INPUT ARGUMENT LIST:
    !     IGDTNUM  - INTEGER GRID DEFINITION TEMPLATE NUMBER.
    !                CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
    !                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
    !                MUST BE "20" FOR POLAR STEREOGRAPHIC GRIDS.
    !     IGDTMPL  - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY.
    !                CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT OF THE
    !                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE FOR SECTION
    !                THREE:
    !                 (1):  SHAPE OF EARTH, OCTET 15
    !                 (2):  SCALE FACTOR OF SPHERICAL EARTH RADIUS,
    !                       OCTET 16
    !                 (3):  SCALED VALUE OF RADIUS OF SPHERICAL EARTH,
    !                       OCTETS 17-20
    !                 (4):  SCALE FACTOR OF MAJOR AXIS OF ELLIPTICAL EARTH,
    !                       OCTET 21
    !                 (5):  SCALED VALUE OF MAJOR AXIS OF ELLIPTICAL EARTH,
    !                       OCTETS 22-25
    !                 (6):  SCALE FACTOR OF MINOR AXIS OF ELLIPTICAL EARTH,
    !                       OCTET 26
    !                 (7):  SCALED VALUE OF MINOR AXIS OF ELLIPTICAL EARTH,
    !                       OCTETS 27-30
    !                 (8):  NUMBER OF POINTS ALONG X-AXIS, OCTETS 31-34
    !                 (9):  NUMBER OF POINTS ALONG Y-AXIS, OCTETS 35-38
    !                 (10): LATITUDE OF FIRST GRID POINT, OCTETS 39-42
    !                 (11): LONGITUDE OF FIRST GRID POINT, OCTETS 43-46
    !                 (12): RESOLUTION AND COMPONENT FLAGS, OCTET 47
    !                 (13): TRUE LATITUDE, OCTETS 48-51
    !                 (14): ORIENTATION LONGITUDE, OCTETS 52-55
    !                 (15): X-DIRECTION GRID LENGTH, OCTETS 56-59
    !                 (16): Y-DIRECTION GRID LENGTH, OCTETS 60-63
    !                 (17): PROJECTION CENTER FLAG, OCTET 64
    !                 (18): SCANNING MODE FLAGS, OCTET 65
    !     IGDTLEN  - INTEGER NUMBER OF ELEMENTS (18) OF THE GRID DEFINITION
    !                TEMPLATE ARRAY.  CORRESPONDS TO THE GFLD%IGDTLEN
    !                COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
    !     IOPT     - INTEGER OPTION FLAG
    !                (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
    !                (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
    !     NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
    !     FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
    !                (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
    !     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
    !     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
    !     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
    !                (ACCEPTABLE RANGE: -360. TO 360.)
    !     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
    !                (ACCEPTABLE RANGE: -90. TO 90.)
    !
    !   OUTPUT ARGUMENT LIST:
    !     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
    !     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
    !     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
    !     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
    !     NRET     - INTEGER NUMBER OF VALID POINTS COMPUTED
    !     CROT     - REAL, OPTIONAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES
    !     SROT     - REAL, OPTIONAL (NPTS) CLOCKWISE VECTOR ROTATION SINES
    !                (UGRID=CROT*UEARTH-SROT*VEARTH;
    !                 VGRID=SROT*UEARTH+CROT*VEARTH)
    !     XLON     - REAL, OPTIONAL (NPTS) DX/DLON IN 1/DEGREES
    !     XLAT     - REAL, OPTIONAL (NPTS) DX/DLAT IN 1/DEGREES
    !     YLON     - REAL, OPTIONAL (NPTS) DY/DLON IN 1/DEGREES
    !     YLAT     - REAL, OPTIONAL (NPTS) DY/DLAT IN 1/DEGREES
    !     AREA     - REAL, OPTIONAL (NPTS) AREA WEIGHTS IN M**2
    !                (PROPORTIONAL TO THE SQUARE OF THE MAP FACTOR)
    !
    ! ATTRIBUTES:
    !   LANGUAGE: FORTRAN 90
    !
    !$$$
    IMPLICIT NONE
    !
 
    class(ip_polar_stereo_grid), intent(in) :: self
    INTEGER,          INTENT(IN   ) :: IOPT, NPTS
    INTEGER,          INTENT(  OUT) :: NRET
    !
    REAL,             INTENT(IN   ) :: FILL
    REAL,             INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
    REAL,             INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
    REAL, OPTIONAL,   INTENT(  OUT) :: CROT(NPTS),SROT(NPTS)
    REAL, OPTIONAL,   INTENT(  OUT) :: XLON(NPTS),XLAT(NPTS)
    REAL, OPTIONAL,   INTENT(  OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
    !
    INTEGER                         :: IM, JM
    INTEGER                         :: ITER, N
    !
    LOGICAL                         :: ELLIPTICAL, LROT, LMAP, LAREA
    !
    REAL                            :: ALAT, ALAT1, ALONG, DIFF
    REAL                            :: DI, DJ, DE
    REAL                            :: DR, E, E_OVER_2
    REAL                            :: MC, SLAT, SLATR
    REAL                            :: RLAT1, RLON1, RHO, T, TC
    REAL                            :: XMAX, XMIN, YMAX, YMIN
    REAL                            :: XP, YP, DR2
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    IF(PRESENT(crot)) crot=fill
    IF(PRESENT(srot)) srot=fill
    IF(PRESENT(xlon)) xlon=fill
    IF(PRESENT(xlat)) xlat=fill
    IF(PRESENT(ylon)) ylon=fill
    IF(PRESENT(ylat)) ylat=fill
    IF(PRESENT(area)) area=fill
 
    elliptical = self%elliptical
    im=self%im
    jm=self%jm
 
    rlat1=self%rlat1
    rlon1=self%rlon1
 
    irot=self%irot
    slatr=self%slatr
    orient=self%orient
 
    h=self%h
    dxs=self%dxs
    dys=self%dys
 
    rerth = self%rerth
    e2 = self%eccen_squared
    !
    ! FIND X/Y OF POLE
    IF (.NOT.elliptical) THEN
       de=(1.+sin(slatr))*rerth
       dr=de*cos(rlat1/dpr)/(1+h*sin(rlat1/dpr))
       xp=1-h*sin((rlon1-orient)/dpr)*dr/dxs
       yp=1+cos((rlon1-orient)/dpr)*dr/dys
       de2=de**2
    ELSE
       e=sqrt(e2)
       e_over_2=e*0.5
       alat=h*rlat1/dpr
       along = (rlon1-orient)/dpr
       t=tan(pi4-alat/2.)/((1.-e*sin(alat))/  &
            (1.+e*sin(alat)))**(e_over_2)
       tc=tan(pi4-slatr/2.)/((1.-e*sin(slatr))/  &
            (1.+e*sin(slatr)))**(e_over_2)
       mc=cos(slatr)/sqrt(1.0-e2*(sin(slatr)**2))
       rho=rerth*mc*t/tc
       yp = 1.0 + rho*cos(h*along)/dys
       xp = 1.0 - rho*sin(h*along)/dxs
    ENDIF ! ELLIPTICAL
    xmin=0
    xmax=im+1
    ymin=0
    ymax=jm+1
    nret=0
    IF(PRESENT(crot).AND.PRESENT(srot))THEN
       lrot=.true.
    ELSE
       lrot=.false.
    ENDIF
    IF(PRESENT(xlon).AND.PRESENT(xlat).AND.PRESENT(ylon).AND.PRESENT(ylat))THEN
       lmap=.true.
    ELSE
       lmap=.false.
    ENDIF
    IF(PRESENT(area))THEN
       larea=.true.
    ELSE
       larea=.false.
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  TRANSLATE GRID COORDINATES TO EARTH COORDINATES
    IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
       IF(.NOT.elliptical)THEN
          !$OMP PARALLEL DO PRIVATE(N,DI,DJ,DR2) REDUCTION(+:NRET) SCHEDULE(STATIC)
          DO n=1,npts
             IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
                  ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
                di=(xpts(n)-xp)*dxs
                dj=(ypts(n)-yp)*dys
                dr2=di**2+dj**2
                IF(dr2.LT.de2*1.e-6) THEN
                   rlon(n)=0.
                   rlat(n)=h*90.
                ELSE
                   rlon(n)=mod(orient+h*dpr*atan2(di,-dj)+3600,360.)
                   rlat(n)=h*dpr*asin((de2-dr2)/(de2+dr2))
                ENDIF
                nret=nret+1
                IF(lrot) CALL polar_stereo_vect_rot(rlon(n),crot(n),srot(n))
                IF(lmap) CALL polar_stereo_map_jacob(rlon(n),rlat(n),dr2, &
                     xlon(n),xlat(n),ylon(n),ylat(n))
                IF(larea) CALL polar_stereo_grid_area(rlat(n),dr2,area(n))
             ELSE
                rlon(n)=fill
                rlat(n)=fill
             ENDIF
          ENDDO
          !$OMP END PARALLEL DO
       ELSE ! ELLIPTICAL
          !$OMP PARALLEL DO PRIVATE(N,DI,DJ,RHO,T,ALONG,ALAT1,ALAT,DIFF) &
          !$OMP& REDUCTION(+:NRET) SCHEDULE(STATIC)
          DO n=1,npts
             IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND.  &
                  ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
                di=(xpts(n)-xp)*dxs
                dj=(ypts(n)-yp)*dys
                rho=sqrt(di*di+dj*dj)
                t=(rho*tc)/(rerth*mc)
                IF(abs(ypts(n)-yp)<0.01)THEN
                   IF(di>0.0) along=orient+h*90.0
                   IF(di<=0.0) along=orient-h*90.0
                ELSE
                   along=orient+h*atan(di/(-dj))*dpr
                   IF(dj>0) along=along+180.
                END IF
                alat1=pi2-2.0*atan(t)
                DO iter=1,10
                   alat = pi2 - 2.0*atan(t*(((1.0-e*sin(alat1))/  &
                        (1.0+e*sin(alat1)))**(e_over_2)))
                   diff = abs(alat-alat1)*dpr
                   IF (diff < 0.000001) EXIT
                   alat1=alat
                ENDDO
                rlat(n)=h*alat*dpr
                rlon(n)=along
                IF(rlon(n)<0.0) rlon(n)=rlon(n)+360.
                IF(rlon(n)>360.0) rlon(n)=rlon(n)-360.0
                nret=nret+1
                IF(lrot) CALL polar_stereo_vect_rot(rlon(n),crot(n),srot(n))
             ELSE
                rlon(n)=fill
                rlat(n)=fill
             ENDIF
          ENDDO
          !$OMP END PARALLEL DO
       ENDIF ! ELLIPTICAL
       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !  TRANSLATE EARTH COORDINATES TO GRID COORDINATES
    ELSEIF(iopt.EQ.-1) THEN
       IF(.NOT.elliptical)THEN
          !$OMP PARALLEL DO PRIVATE(N,DR,DR2) REDUCTION(+:NRET) SCHEDULE(STATIC)
          DO n=1,npts
             IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90.AND. &
                  h*rlat(n).NE.-90) THEN
                dr=de*tan((90-h*rlat(n))/2/dpr)
                dr2=dr**2
                xpts(n)=xp+h*sin((rlon(n)-orient)/dpr)*dr/dxs
                ypts(n)=yp-cos((rlon(n)-orient)/dpr)*dr/dys
                IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
                     ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
                   nret=nret+1
                   IF(lrot) CALL polar_stereo_vect_rot(rlon(n),crot(n),srot(n))
                   IF(lmap) CALL polar_stereo_map_jacob(rlon(n),rlat(n),dr2, &
                        xlon(n),xlat(n),ylon(n),ylat(n))
                   IF(larea) CALL polar_stereo_grid_area(rlat(n),dr2,area(n))
                ELSE
                   xpts(n)=fill
                   ypts(n)=fill
                ENDIF
             ELSE
                xpts(n)=fill
                ypts(n)=fill
             ENDIF
          ENDDO
          !$OMP END PARALLEL DO
       ELSE  ! ELLIPTICAL CASE
          !$OMP PARALLEL DO PRIVATE(N,ALAT,ALONG,T,RHO) REDUCTION(+:NRET) SCHEDULE(STATIC)
          DO n=1,npts
             IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90.AND.  &
                  h*rlat(n).NE.-90) THEN
                alat = h*rlat(n)/dpr
                along = (rlon(n)-orient)/dpr
                t=tan(pi4-alat*0.5)/((1.-e*sin(alat))/  &
                     (1.+e*sin(alat)))**(e_over_2)
                rho=rerth*mc*t/tc
                xpts(n)= xp + rho*sin(h*along) / dxs
                ypts(n)= yp - rho*cos(h*along) / dys
                IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND.  &
                     ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
                   nret=nret+1
                   IF(lrot) CALL polar_stereo_vect_rot(rlon(n),crot(n),srot(n))
                ELSE
                   xpts(n)=fill
                   ypts(n)=fill
                ENDIF
             ELSE
                xpts(n)=fill
                ypts(n)=fill
             ENDIF
          ENDDO
          !$OMP END PARALLEL DO
       ENDIF
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  END SUBROUTINE gdswzd_polar_stereo
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  SUBROUTINE polar_stereo_error(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  POLAR_STEREO_ERROR    ERROR HANDLER
    !   PRGMMR: GAYNO       ORG: W/NMC23       DATE: 2015-07-13
    !
    ! ABSTRACT: UPON AN ERROR, THIS SUBPROGRAM ASSIGNS
    !           A "FILL" VALUE TO THE OUTPUT FIELDS.
 
    ! PROGRAM HISTORY LOG:
    ! 2015-07-13  GAYNO     INITIAL VERSION
    !
    ! USAGE:    CALL POLAR_STEREO_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
    !
    !   INPUT ARGUMENT LIST:
    !     IOPT     - INTEGER OPTION FLAG
    !                (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
    !                (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
    !     NPTS     - INTEGER MAXIMUM NUMBER OF COORDINATES
    !     FILL     - REAL FILL VALUE TO SET INVALID OUTPUT DATA
    !                (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
    !   OUTPUT ARGUMENT LIST:
    !     RLON     - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
    !     RLAT     - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
    !     XPTS     - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
    !     YPTS     - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
    !
    ! ATTRIBUTES:
    !   LANGUAGE: FORTRAN 90
    !
    !$$$
    IMPLICIT NONE
    !
    INTEGER, INTENT(IN   ) :: IOPT, NPTS
    !
    REAL,    INTENT(IN   ) :: FILL
    REAL,    INTENT(  OUT) :: RLAT(NPTS),RLON(NPTS)
    REAL,    INTENT(  OUT) :: XPTS(NPTS),YPTS(NPTS)
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    IF(iopt>=0) THEN
       rlon=fill
       rlat=fill
    ENDIF
    IF(iopt<=0) THEN
       xpts=fill
       ypts=fill
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  END SUBROUTINE polar_stereo_error
  !
  SUBROUTINE polar_stereo_vect_rot(RLON, CROT, SROT)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  POLAR_STEREO_VECT_ROT   VECTOR ROTATION FIELDS FOR
    !                                      POLAR STEREOGRAPHIC GRIDS.
    !
    !   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
    !
    ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE VECTOR ROTATION SINES AND
    !           COSINES FOR A POLAR STEREOGRAPHIC AZIMUTHAL GRID.
    !
    ! PROGRAM HISTORY LOG:
    ! 2015-01-21  GAYNO    INITIAL VERSION
    ! 2015-09-17  GAYNO    RENAME AS "POLAR_STEREO_VECT_ROT"
    !
    ! USAGE:    CALL POLAR_STEREO_VECT_ROT(RLON,CROT,SROT)
    !
    !   INPUT ARGUMENT LIST:
    !     RLON     - GRID POINT LONGITUDE IN DEGREES (REAL)
    !
    !   OUTPUT ARGUMENT LIST:
    !     CROT     - CLOCKWISE VECTOR ROTATION COSINES (REAL)
    !     SROT     - CLOCKWISE VECTOR ROTATION SINES (REAL)
    !                (UGRID=CROT*UEARTH-SROT*VEARTH;
    !                 VGRID=SROT*UEARTH+CROT*VEARTH)
    !
    ! ATTRIBUTES:
    !   LANGUAGE: FORTRAN 90
    !
    !$$$
    !
    IMPLICIT NONE
 
    REAL,             INTENT(IN   ) :: RLON
    REAL,             INTENT(  OUT) :: CROT, SROT
 
    IF(irot.EQ.1) THEN
       crot=h*cos((rlon-orient)/dpr)
       srot=sin((rlon-orient)/dpr)
    ELSE
       crot=1.
       srot=0.
    ENDIF
 
  END SUBROUTINE polar_stereo_vect_rot
  !
  SUBROUTINE polar_stereo_map_jacob(RLON,RLAT,DR2,XLON,XLAT,YLON,YLAT)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  POLAR_STEREO_MAP_JACOB  MAP JACOBIANS FOR
    !                                      POLAR STEREOGRAPHIC GRIDS.
    !
    !   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
    !
    ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE MAP JACOBIANS FOR
    !           A POLAR STEREOGRAPHIC AZIMUTHAL GRID (SPHERICAL
    !           EARTH).
    !
    ! PROGRAM HISTORY LOG:
    ! 2015-01-21  GAYNO    INITIAL VERSION
    ! 2015-09-17  GAYNO    RENAME AS "POLAR_STEREO_MAP_JACOB"
    ! 2018-07-20  WESLEY   PASS IN DR2 FOR THREADING.
    !
    ! USAGE:  CALL POLAR_STEREO_MAP_JACOB(RLON,RLAT,DR2,XLON,XLAT,YLON,YLAT)
    !
    !   INPUT ARGUMENT LIST:
    !     RLON     - LONGITUDE IN DEGREES (REAL)
    !     RLAT     - LATITUDE IN DEGREES (REAL)
    !     DR2      - SQUARED DISTANCE FROM POLE (REAL)
    !
    !   OUTPUT ARGUMENT LIST:
    !     XLON     - DX/DLON IN 1/DEGREES (REAL)
    !     XLAT     - DX/DLAT IN 1/DEGREES (REAL)
    !     YLON     - DY/DLON IN 1/DEGREES (REAL)
    !     YLAT     - DY/DLAT IN 1/DEGREES (REAL)
    !
    ! ATTRIBUTES:
    !   LANGUAGE: FORTRAN 90
    !
    !$$$
    !
    IMPLICIT NONE
 
    REAL,             INTENT(IN   ) :: RLON, RLAT, DR2
    REAL,             INTENT(  OUT) :: XLON, XLAT, YLON, YLAT
 
    REAL                            :: CLAT, DE, DR
 
    IF(dr2.LT.de2*1.e-6) THEN
       de=sqrt(de2)
       xlon=0.
       xlat=-sin((rlon-orient)/dpr)/dpr*de/dxs/2
       ylon=0.
       ylat=h*cos((rlon-orient)/dpr)/dpr*de/dys/2
    ELSE
       dr=sqrt(dr2)
       clat=cos(rlat/dpr)
       xlon=h*cos((rlon-orient)/dpr)/dpr*dr/dxs
       xlat=-sin((rlon-orient)/dpr)/dpr*dr/dxs/clat
       ylon=sin((rlon-orient)/dpr)/dpr*dr/dys
       ylat=h*cos((rlon-orient)/dpr)/dpr*dr/dys/clat
    ENDIF
 
  END SUBROUTINE polar_stereo_map_jacob
  !
  SUBROUTINE polar_stereo_grid_area(RLAT, DR2, AREA)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  POLAR_STEREO_GRID_AREA  GRID BOX AREA FOR
    !                                      POLAR STEREOGRAPHIC GRIDS
    !
    !   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
    !
    ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE GRID BOX AREA FOR
    !           A POLAR STEREOGRAPHIC AZIMUTHAL GRID (SPHERICAL
    !           EARTH).
    !
    ! PROGRAM HISTORY LOG:
    ! 2015-01-21  GAYNO    INITIAL VERSION
    ! 2015-09-17  GAYNO    RENAME AS "POLAR_STEREO_GRID_AREA".
    ! 2018-07-20  WESLEY   PASS IN DR2 FOR THREADING.
    !
    ! USAGE:  CALL POLAR_STEREO_GRID_AREA(RLAT,DR2,AREA)
    !
    !   INPUT ARGUMENT LIST:
    !     RLAT     - LATITUDE OF GRID POINT IN DEGREES (REAL)
    !     DR2      - SQUARED DISTANCE FROM POLE (REAL)
    !
    !   OUTPUT ARGUMENT LIST:
    !     AREA     - AREA WEIGHTS IN M**2 (REAL)
    !
    ! ATTRIBUTES:
    !   LANGUAGE: FORTRAN 90
    !
    !$$$
    !
    IMPLICIT NONE
 
    REAL,             INTENT(IN   ) :: RLAT, DR2
    REAL,             INTENT(  OUT) :: AREA
 
    REAL                            :: CLAT
 
    IF(dr2.LT.de2*1.e-6) THEN
       area=rerth**2*abs(dxs)*abs(dys)*4/de2
    ELSE
       clat=cos(rlat/dpr)
       area=rerth**2*clat**2*abs(dxs)*abs(dys)/dr2
    ENDIF
 
  END SUBROUTINE polar_stereo_grid_area
 
 
end module ip_polar_stereo_grid_mod