ip_rot_equid_cylind_grid_mod.f90

module ip_rot_equid_cylind_grid_mod
  use iso_fortran_env, only: real64
  use ip_grid_descriptor_mod
  use ip_grid_mod
  use constants_mod, only: dpr, pi
  use earth_radius_mod
  implicit none
 
  private
  public :: ip_rot_equid_cylind_grid
 
  integer, parameter :: kd = real64
 
  type, extends(ip_grid) :: ip_rot_equid_cylind_grid
     real(kd) :: clat0, dlats, dlons, rlon0, slat0, wbd, sbd
     integer :: irot
   contains
     procedure :: init_grib1
     procedure :: init_grib2
     procedure :: gdswzd => gdswzd_rot_equid_cylind
  end type ip_rot_equid_cylind_grid
 
  INTEGER                                :: IROT
 
  REAL(KIND=kd)                          :: rerth
  REAL(KIND=kd)                          :: clat0, dlats, dlons
  REAL(KIND=kd)                          :: rlon0, slat0
 
 
CONTAINS
 
  subroutine init_grib1(self, g1_desc)
    class(ip_rot_equid_cylind_grid), intent(inout) :: self
    type(grib1_descriptor), intent(in) :: g1_desc
 
    real(kd) :: rlat1, rlon1, rlat0, rlat2, rlon2, nbd, ebd
    real(kd) :: hs, hs2, slat1, slat2, slatr, clon1, clon2, clat1, clat2, clatr, clonr, rlonr, rlatr
    integer :: iscale
 
    associate(kgds => g1_desc%gds)
      self%rerth = 6.3712e6_kd
      self%eccen_squared = 0d0
 
      rlat1=kgds(4)*1.e-3_kd
      rlon1=kgds(5)*1.e-3_kd
      rlat0=kgds(7)*1.e-3_kd
      self%RLON0=kgds(8)*1.e-3_kd
      rlat2=kgds(12)*1.e-3_kd
      rlon2=kgds(13)*1.e-3_kd
 
      self%IROT=mod(kgds(6)/8,2)
      self%IM=kgds(2)
      self%JM=kgds(3)
 
      slat1=sin(rlat1/dpr)
      clat1=cos(rlat1/dpr)
      self%SLAT0=sin(rlat0/dpr)
      self%CLAT0=cos(rlat0/dpr)
 
      hs=sign(1._kd,mod(rlon1-self%RLON0+180+3600,360._kd)-180)
      clon1=cos((rlon1-self%RLON0)/dpr)
      slatr=self%CLAT0*slat1-self%SLAT0*clat1*clon1
      clatr=sqrt(1-slatr**2)
      clonr=(self%CLAT0*clat1*clon1+self%SLAT0*slat1)/clatr
      rlatr=dpr*asin(slatr)
      rlonr=hs*dpr*acos(clonr)
 
      self%WBD=rlonr
      self%SBD=rlatr
      slat2=sin(rlat2/dpr)
      clat2=cos(rlat2/dpr)
      hs2=sign(1._kd,mod(rlon2-self%RLON0+180+3600,360._kd)-180)
      clon2=cos((rlon2-self%RLON0)/dpr)
      slatr=self%CLAT0*slat2-self%SLAT0*clat2*clon2
      clatr=sqrt(1-slatr**2)
      clonr=(self%CLAT0*clat2*clon2+self%SLAT0*slat2)/clatr
      nbd=dpr*asin(slatr)
      ebd=hs2*dpr*acos(clonr)
      self%DLATS=(nbd-self%SBD)/float(self%JM-1)
      self%DLONS=(ebd-self%WBD)/float(self%IM-1)
 
      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_rot_equid_cylind_grid), intent(inout) :: self
    type(grib2_descriptor), intent(in) :: g2_desc
 
    real(kd) :: rlat1, rlon1, rlat0, rlat2, rlon2, nbd, ebd
    integer :: iscale
    integer :: i_offset_odd, i_offset_even, j_offset
 
    associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
 
      CALL earth_radius(igdtmpl,igdtlen,self%rerth,self%eccen_squared)
 
      i_offset_odd=mod(igdtmpl(19)/8,2)
      i_offset_even=mod(igdtmpl(19)/4,2)
      j_offset=mod(igdtmpl(19)/2,2)
 
      iscale=igdtmpl(10)*igdtmpl(11)
      IF(iscale==0) iscale=10**6
 
      rlat1=float(igdtmpl(12))/float(iscale)
      rlon1=float(igdtmpl(13))/float(iscale)
      rlat0=float(igdtmpl(20))/float(iscale)
      rlat0=rlat0+90.0_kd
 
      self%RLON0=float(igdtmpl(21))/float(iscale)
 
      rlat2=float(igdtmpl(15))/float(iscale)
      rlon2=float(igdtmpl(16))/float(iscale)
 
      self%IROT=mod(igdtmpl(14)/8,2)
      self%IM=igdtmpl(8)
      self%JM=igdtmpl(9)
 
      self%SLAT0=sin(rlat0/dpr)
      self%CLAT0=cos(rlat0/dpr)
 
      self%WBD=rlon1
      IF (self%WBD > 180.0) self%WBD = self%WBD - 360.0
      self%SBD=rlat1
 
      nbd=rlat2
      ebd=rlon2
 
      self%DLATS=(nbd-self%SBD)/float(self%JM-1)
      self%DLONS=(ebd-self%WBD)/float(self%IM-1)
 
      IF(i_offset_odd==1) self%WBD=self%WBD+(0.5_kd*self%DLONS)
      IF(j_offset==1) self%SBD=self%SBD+(0.5_kd*self%DLATS)
 
      self%iwrap = 0
      self%jwrap1 = 0
      self%jwrap2 = 0
      self%kscan = 0
      self%nscan = mod(igdtmpl(19) / 32, 2)
      self%nscan_field_pos = self%nscan
    end associate
  end subroutine init_grib2
 
 
 
  SUBROUTINE gdswzd_rot_equid_cylind(self,IOPT,NPTS, &
       FILL,XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  GDSWZD_ROT_EQUID_CYLIND  GDS WIZARD FOR ROTATED 
    !                                       EQUIDISTANT CYLINDRICAL
    !   PRGMMR: GAYNO       ORG: W/NMC23       DATE: 2007-NOV-15
    !
    ! 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 NON-"E" STAGGERED ROTATED EQUIDISTANT CYLINDRICAL 
    !           PROJECTIONS. THE SCAN MODE (SECTION 3, OCTET 72, BITS 5-6)
    !           DETERMINE WHETHER THIS IS AN "H" OR "V" GRID.  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, THE MAP JACOBIANS AND
    !           THE GRID BOX AREAS MAY BE RETURNED AS WELL.  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:
    ! 2010-JAN-15  GAYNO     BASED ON ROUTINES GDSWZDCB AND GDSWZDCA
    ! 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.
    ! 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."
    ! 2018-07-20   WESLEY    ADD THREADS.
    !
    ! USAGE:   CALL GDSWZD_ROT_EQUID_CYLIND(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.
    !     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 A PARALLEL, OCTS 31-34
    !                 (9):  NUMBER OF POINTS ALONG A MERIDIAN, OCTS 35-38
    !                 (10): BASIC ANGLE OF INITIAL PRODUCTION DOMAIN,
    !                       OCTETS 39-42
    !                 (11): SUBDIVISIONS OF BASIC ANGLE, OCTETS 43-46
    !                 (12): LATITUDE OF FIRST GRID POINT IN X/Y SPACE
    !                       (BEFORE ROTATION), OCTETS 47-50
    !                 (13): LONGITUDE OF FIRST GRID POINT IN X/Y
    !                       SPACE (BEFORE ROTATION), OCTETS 51-54
    !                 (14): RESOLUTION AND COMPONENT FLAGS, OCTET 55
    !                 (15): LATITUDE OF LAST GRID POINT IN X/Y SPACE
    !                       (BEFORE ROTATION), OCTETS 56-59
    !                 (16): LONGITUDE OF LAST GRID POINT IN X/Y SPACE
    !                       (BEFORE ROTATION), OCTETS 60-63
    !                 (17): I-DIRECTION INCREMENT, OCTETS 64-67
    !                 (18): J-DIRECTION INCREMENT, OCTETS 68-71
    !                 (19): SCANNING MODE, OCTET 72
    !                 (20): LATITUDE OF SOUTHERN POLE OF PROJECTION,
    !                       OCTETS 73-76
    !                 (21): LONGITUDE OF SOUTHERN POLE OF PROJECTION,
    !                       OCTETS 77-80
    !                 (22): ANGLE OF ROTATION OF PROJECTION, OCTS 81-84
    !     IGDTLEN  - INTEGER NUMBER OF ELEMENTS (22) 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
    !
    ! ATTRIBUTES:
    !   LANGUAGE: FORTRAN 90
    !
    !$$$
    IMPLICIT NONE
 
    class(ip_rot_equid_cylind_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,N
    !
    LOGICAL                                :: LROT, LMAP, LAREA
    !
    REAL                                   :: DUM1,DUM2
    REAL(KIND=kd)                          :: hs
    REAL(KIND=kd)                          :: clonr,clatr,slatr
    REAL(KIND=kd)                          :: clat,slat,clon
    REAL(KIND=kd)                          :: rlatr,rlonr
    REAL(KIND=kd)                          :: wbd,sbd
    REAL                                   :: XMIN,XMAX,YMIN,YMAX
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    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
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! IS THE EARTH RADIUS DEFINED?
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  IS THIS AN "E"-STAGGER GRID?  ROUTINE CAN'T PROCESS THOSE.
    ! I_OFFSET_ODD=MOD(IGDTMPL(19)/8,2)
    ! I_OFFSET_EVEN=MOD(IGDTMPL(19)/4,2)
    ! J_OFFSET=MOD(IGDTMPL(19)/2,2)
    ! IF(I_OFFSET_ODD/=I_OFFSET_EVEN) THEN
    !    CALL ROT_EQUID_CYLIND_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
    !    RETURN
    ! ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
 
    rlon0=self%rlon0
    irot=self%irot
 
    im=self%im
    jm=self%jm
 
    slat0=self%slat0
    clat0=self%clat0
 
    wbd=self%wbd
    sbd=self%sbd
 
    dlats=self%dlats
    dlons=self%dlons
 
    xmin=0
    xmax=im+1
    ymin=0
    ymax=jm+1
    nret=0
 
    rerth = self%rerth
    IF(rerth<0.)THEN
       CALL rot_equid_cylind_error(iopt,fill,rlat,rlon,xpts,ypts,npts)
       RETURN
    ENDIF
 
    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
       !$OMP PARALLEL DO PRIVATE(N,RLONR,RLATR,HS,CLONR,SLATR,CLATR,SLAT,CLAT,CLON) &
       !$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
             rlonr=wbd+(xpts(n)-1._kd)*dlons
             rlatr=sbd+(ypts(n)-1._kd)*dlats
             IF(rlonr <= 0._kd) THEN
                hs=-1.0_kd
             ELSE
                hs=1.0_kd
             ENDIF
             clonr=cos(rlonr/dpr)
             slatr=sin(rlatr/dpr)
             clatr=cos(rlatr/dpr)
             slat=clat0*slatr+slat0*clatr*clonr
             IF(slat.LE.-1) THEN
                clat=0.
                clon=cos(rlon0/dpr)
                rlon(n)=0.
                rlat(n)=-90.
             ELSEIF(slat.GE.1) THEN
                clat=0.
                clon=cos(rlon0/dpr)
                rlon(n)=0.
                rlat(n)=90.
             ELSE
                clat=sqrt(1-slat**2)
                clon=(clat0*clatr*clonr-slat0*slatr)/clat
                clon=min(max(clon,-1._kd),1._kd)
                rlon(n)=mod(rlon0+hs*dpr*acos(clon)+3600,360._kd)
                rlat(n)=dpr*asin(slat)
             ENDIF
             nret=nret+1
             IF(lrot) CALL rot_equid_cylind_vect_rot(rlon(n), clatr, slatr, &
                  clat, slat, clon, crot(n), srot(n))
             IF(lmap) CALL rot_equid_cylind_map_jacob(fill, rlon(n), clatr, &
                  clat, slat, clon, xlon(n), xlat(n), ylon(n), ylat(n))
             IF(larea) CALL rot_equid_cylind_grid_area(clatr, fill, area(n))
          ELSE
             rlon(n)=fill
             rlat(n)=fill
          ENDIF
       ENDDO
       !$OMP END PARALLEL DO
       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !  TRANSLATE EARTH COORDINATES TO GRID COORDINATES
    ELSEIF(iopt.EQ.-1) THEN
       !$OMP PARALLEL DO PRIVATE(N,HS,CLON,SLAT,CLAT,SLATR,CLATR,CLONR,RLONR,RLATR) &
       !$OMP& REDUCTION(+:NRET) SCHEDULE(STATIC)
       DO n=1,npts
          IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90) THEN
             hs=sign(1._kd,mod(rlon(n)-rlon0+180+3600,360._kd)-180)
             clon=cos((rlon(n)-rlon0)/dpr)
             slat=sin(rlat(n)/dpr)
             clat=cos(rlat(n)/dpr)
             slatr=clat0*slat-slat0*clat*clon
             IF(slatr.LE.-1) THEN
                clatr=0._kd
                rlonr=0.
                rlatr=-90.
             ELSEIF(slatr.GE.1) THEN
                clatr=0._kd
                rlonr=0.
                rlatr=90.
             ELSE
                clatr=sqrt(1-slatr**2)
                clonr=(clat0*clat*clon+slat0*slat)/clatr
                clonr=min(max(clonr,-1._kd),1._kd)
                rlonr=hs*dpr*acos(clonr)
                rlatr=dpr*asin(slatr)
             ENDIF
             xpts(n)=(rlonr-wbd)/dlons+1._kd
             ypts(n)=(rlatr-sbd)/dlats+1._kd
             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 rot_equid_cylind_vect_rot(rlon(n), clatr, slatr, &
                     clat, slat, clon, crot(n), srot(n))
                IF(lmap) CALL rot_equid_cylind_map_jacob(fill, rlon(n), clatr, &
                     clat, slat, clon, xlon(n), xlat(n), ylon(n), ylat(n))
                IF(larea) CALL rot_equid_cylind_grid_area(clatr, fill, area(n))
             ELSE
                xpts(n)=fill
                ypts(n)=fill
             ENDIF
          ELSE
             xpts(n)=fill
             ypts(n)=fill
          ENDIF
       ENDDO
       !$OMP END PARALLEL DO
    ENDIF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  END SUBROUTINE gdswzd_rot_equid_cylind
  ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  SUBROUTINE rot_equid_cylind_error(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  ROT_EQUID_CYLIND_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 ROT_EQUID_CYLIND_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 rot_equid_cylind_error
  !
  SUBROUTINE rot_equid_cylind_vect_rot(RLON, CLATR, SLATR, CLAT, SLAT, &
       CLON, CROT, SROT)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  ROT_EQUID_CYLIND_VECT_ROT  VECTOR ROTATION FIELDS FOR
    !                                         ROTATED EQUIDISTANT CYLINDRICAL
    !                                         GRIDS - NON "E" STAGGER.
    !
    !   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
    !
    ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE VECTOR ROTATION SINES AND
    !           COSINES FOR A ROTATED EQUIDISTANT CYLINDRICAL GRID -
    !           NON "E" STAGGER.
    !
    ! PROGRAM HISTORY LOG:
    ! 2015-01-21  GAYNO    INITIAL VERSION
    ! 2015-07-19  GAYNO    RENAME AS "ROT_EQUID_CYLIND_VECT_ROT."
    ! 2018-07-20  WESLEY   PASS IN CLATR, SLATR, CLAT, SLAT, CLON
    !                      FOR THREADING.
    !
    ! USAGE: CALL ROT_EQUID_CYLIND_VECT_ROT(RLON, CLATR, SLATR, CLAT, SLAT, &
    !                                       CLON, CROT, SROT)
    !
    !   INPUT ARGUMENT LIST:
    !     RLON     - LONGITUDE IN DEGREES (REAL)
    !     CLATR    - COSINE OF ROTATED LATITUDE (REAL)
    !     SLATR    - SINE OF ROTATED LATITUDE (REAL)
    !     CLAT     - COSINE OF LATITUDE (REAL)
    !     SLAT     - SINE OF LATITUDE (REAL)
    !     CLON     - COSINE OF LONGITUDE (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(KIND=kd),    INTENT(IN   ) :: clat, clatr, clon, slat, slatr
    REAL         ,    INTENT(IN   ) :: RLON
    REAL         ,    INTENT(  OUT) :: CROT, SROT
 
    REAL(KIND=kd)                   :: slon
 
    IF(irot.EQ.1) THEN
       IF(clatr.LE.0._kd) THEN
          crot=-sign(1._kd,slatr*slat0)
          srot=0.
       ELSE
          slon=sin((rlon-rlon0)/dpr)
          crot=(clat0*clat+slat0*slat*clon)/clatr
          srot=slat0*slon/clatr
       ENDIF
    ELSE
       crot=1.
       srot=0.
    ENDIF
 
  END SUBROUTINE rot_equid_cylind_vect_rot
  !
  SUBROUTINE rot_equid_cylind_map_jacob(FILL, RLON, CLATR, CLAT, &
       SLAT, CLON, XLON, XLAT, YLON, YLAT)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  ROT_EQUID_CYLIND_MAP_JACOB 
    !     MAP JACOBIANS FOR ROTATED EQUIDISTANT CYLINDRICAL
    !     GRIDS - NON "E" STAGGER.
    !
    !   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
    !
    ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE MAP JACOBIANS FOR
    !           A ROTATED EQUIDISTANT CYLINDRICAL GRID -
    !           NON "E" STAGGER.
    !
    ! PROGRAM HISTORY LOG:
    ! 2015-01-21  GAYNO    INITIAL VERSION
    ! 2015-09-17  GAYNO    RENAME AS "ROT_EQUID_CYLIND_MAP_JACOB".
    ! 2018-07-20  WESLEY   PASS IN CLATR, CLAT, SLAT, CLON TO ALLOW
    !                      THREADING.
    !
    ! USAGE:  CALL ROT_EQUID_CYLIND_MAP_JACOB(FILL, RLON, CLATR, CLAT, &
    !                               SLAT, CLON, XLON, XLAT, YLON, YLAT)
    !
    !   INPUT ARGUMENT LIST:
    !     CLATR    - COSINE OF UNROTATED LATITUDE (REAL)
    !     CLAT     - COSINE OF LATITUDE (REAL)
    !     SLAT     - SINE OF LATITUDE (REAL)
    !     CLON     - COSINE OF LATITUDE (REAL)
    !     FILL     - FILL VALUE FOR UNDEFINED POINTS (REAL)
    !     RLON     - LONGITUDE IN DEGREES (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(KIND=kd),    INTENT(IN   ) :: clatr, clat, slat, clon
    REAL         ,    INTENT(IN   ) :: FILL, RLON
    REAL         ,    INTENT(  OUT) :: XLON, XLAT, YLON, YLAT
 
    REAL(KIND=kd)                   :: slon, term1, term2
 
    IF(clatr.LE.0._kd) THEN
       xlon=fill
       xlat=fill
       ylon=fill
       ylat=fill
    ELSE
       slon=sin((rlon-rlon0)/dpr)
       term1=(clat0*clat+slat0*slat*clon)/clatr
       term2=slat0*slon/clatr
       xlon=term1*clat/(dlons*clatr)
       xlat=-term2/(dlons*clatr)
       ylon=term2*clat/dlats
       ylat=term1/dlats
    ENDIF
 
  END SUBROUTINE rot_equid_cylind_map_jacob
  !
  SUBROUTINE rot_equid_cylind_grid_area(CLATR, FILL, AREA)
    !$$$  SUBPROGRAM DOCUMENTATION BLOCK
    !
    ! SUBPROGRAM:  ROT_EQUID_CYLIND_GRID_AREA  GRID BOX AREA FOR
    !                                          ROTATED EQUIDISTANT CYLINDRICAL
    !                                          GRIDS - NON "E" STAGGER.
    !
    !   PRGMMR: GAYNO     ORG: W/NMC23       DATE: 2015-01-21
    !
    ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE GRID BOX AREA FOR
    !           A ROTATED EQUIDISTANT CYLINDRICAL GRID -
    !           NON "E" STAGGER.
    !
    ! PROGRAM HISTORY LOG:
    ! 2015-01-21  GAYNO    INITIAL VERSION
    ! 2015-07-19  GAYNO    RENAME AS "ROT_EQUID_CYLIND_GRID_AREA."
    ! 2018-07-20  WESLEY   PASS IN CLATR FOR THREADING
    !
    ! USAGE:  CALL ROT_EQUID_CYLIND_GRID_AREA(FILL,AREA,CLATR)
    !
    !   INPUT ARGUMENT LIST:
    !     CLATR    - COSINE OF UNROTATED LATITUDE (REAL)
    !     FILL     - FILL VALUE FOR UNDEFINED POINTS (REAL)
    !
    !   OUTPUT ARGUMENT LIST:
    !     AREA     - AREA WEIGHTS IN M**2 (REAL)
    !
    ! ATTRIBUTES:
    !   LANGUAGE: FORTRAN 90
    !
    !$$$
    !
    IMPLICIT NONE
 
    REAL(KIND=kd),    INTENT(IN   ) :: clatr
    REAL,             INTENT(IN   ) :: FILL
    REAL,             INTENT(  OUT) :: AREA
 
    IF(clatr.LE.0._kd) THEN
       area=fill
    ELSE
       area=2._kd*(rerth**2)*clatr*(dlons/dpr)*sin(0.5_kd*dlats/dpr)
    ENDIF
 
  END SUBROUTINE rot_equid_cylind_grid_area
 
end module ip_rot_equid_cylind_grid_mod