ipxwafs3.f90

 SUBROUTINE ipxwafs3(IDIR, NUMPTS_THIN, NUMPTS_FULL, KM, NUM_OPT, OPT_PTS, &
                    IGDTLEN, IGDTMPL_THIN, DATA_THIN, IB_THIN, BITMAP_THIN,  &
                    IGDTMPL_FULL, DATA_FULL, IB_FULL, BITMAP_FULL, IRET)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM:  IPXWAFS3   EXPAND OR CONTRACT WAFS GRIDS
!   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
!
! ABSTRACT: THIS SUBPROGRAM TRANSFORMS BETWEEN THE THINNED WAFS GRIDS
!           USED FOR TRANSMITTING TO THE AVIATION COMMUNITY
!           AND THEIR FULL EXPANSION AS USED FOR GENERAL INTERPOLATION
!           AND GRAPHICS.  THE THINNED WAFS GRIDS ARE LATITUDE-LONGITUDE
!           GRIDS WHERE THE NUMBER OF POINTS IN EACH ROW DECREASE
!           TOWARD THE POLE.  THIS INFORMATION IS STORED
!           IN THE GRIB 2 GRID DEFINITION TEMPLATE (SECTION 3) 
!           STARTING AT OCTET 73. THE FULL GRID COUNTERPARTS
!           HAVE AN EQUAL NUMBER OF POINTS PER ROW.
!           THE TRANSFORM BETWEEN THE FULL AND THINNED WAFS
!           GRID IS DONE BY NEAREST NEIGHBOR AND IS
!           NOT REVERSIBLE.  THIS ROUTINE WORKS WITH
!           BITMAPPED DATA.
!
! PROGRAM HISTORY LOG:
!   07-07-13  TROJAN   - INITIAL VERSION BASED ON IPXWAFS2
!   2015-JUL  GAYNO    - CONVERT TO GRIB 2
!
! USAGE:   CALL IPXWAFS3(IDIR, NUMPTS_THIN, NUMPTS_FULL, KM, NUM_OPT, OPT_PTS, &
!                   IGDTLEN, IGDTMPL_THIN, DATA_THIN, IB_THIN, BITMAP_THIN,  &
!                   IGDTMPL_FULL, DATA_FULL, IB_FULL, BITMAP_FULL, IRET)
!
!   INPUT ARGUMENT LIST:
!     IDIR         - INTEGER TRANSFORM OPTION
!                   (+1 TO EXPAND THINNED FIELDS TO FULL FIELDS)
!                   (-1 TO CONTRACT FULL FIELDS TO THINNED FIELDS)
!     NUMPTS_THIN  - INTEGER NUMBER OF GRID POINTS - THINNED GRID.  MUST BE
!                    3447.
!     NUMPTS_FULL  - INTEGER NUMBER OF GRID POINTS - FULL GRID. MUST
!                    BE 5329.
!     KM           - INTEGER NUMBER OF FIELDS TO TRANSFORM
!     NUM_OPT      - INTEGER NUMBER OF VALUES TO DESCRIBE THE THINNED
!                    GRID.  MUST BE 73.  DIMENSION OF ARRAY OPT_PTS.
!     OPT_PTS      - INTEGER (NUM_OPT) NUMBER OF GRID POINTS PER ROW -
!                    THINNED GRID - IF IDIR=+1
!     IGDTLEN      - INTEGER GRID DEFINTION TEMPLATE ARRAY LENGTH.  MUST BE
!                    19 FOR LAT/LON GRIDS. CORRESPONDS TO THE GFLD%IGDTLEN
!                    COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                    SAME FOR THIN AND FULL GRIDS WHICH ARE BOTH LAT/LON.
!     IGDTMPL_THIN - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY -
!                    THINNED GRID - IF IDIR=+1. CORRESPONDS TO THE
!                    GFLD%IGDTMPL COMPONENT OF THE NCEP G2 LIBRARY
!                    GRIDMOD DATA STRUCTURE (SECTION 3 INFO):
!                    (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):  SET TO MISSING FOR THINNED GRID., 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, OCTETS 47-50
!                    (13): LONGITUDE OF FIRST GRID POINT, OCTETS 51-54
!                    (14): RESOLUTION AND COMPONENT FLAGS, OCTET 55
!                    (15): LATITUDE OF LAST GRID POINT, OCTETS 56-59
!                    (16): LONGITUDE OF LAST GRID POINT, OCTETS 60-63
!                    (17): SET TO MISSING FOR THINNED GRID, OCTETS 64-67
!                    (18): J-DIRECTION INCREMENT, OCTETS 68-71
!                    (19): SCANNING MODE, OCTET 72
!     DATA_THIN    - REAL (NUMPTS_THIN,KM) THINNED GRID FIELDS IF IDIR=+1
!     IB_THIN      - INTEGER (KM) BITMAP FLAGS THINNED GRID - IF IDIR=+1
!     BITMAP_THIN  - LOGICAL (NUMPTS_THIN,KM) BITMAP FIELDS THIN GRID - IF IDIR=+1
!     IGDTMPL_FULL - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY -
!                    FULL GRID - IF IDIR=-1. CORRESPONDS TO THE
!                    GFLD%IGDTMPL COMPONENT OF THE NCEP G2 LIBRARY
!                    GRIDMOD DATA STRUCTURE. SAME AS IGDTMPL_THIN
!                    EXCEPT:
!                    (8):  NUMBER OF POINTS ALONG A PARALLEL, OCTS 31-34
!                    (17): I-DIRECTION INCREMENT, OCTETS 64-67
!     DATA_FULL    - REAL (NUMPTS_FULL,KM) FULL GRID FIELDS IF IDIR=-1
!     IB_FULL      - INTEGER (KM) BITMAP FLAGS FULL GRID - IF IDIR=-1
!     BITMAP_FULL  - LOGICAL (NUMPTS_FULL,KM) BITMAP FIELDS FULL GRID - IF IDIR=-1
!
!   OUTPUT ARGUMENT LIST:
!     OPT_PTS      - INTEGER (NUM_OPT) NUMBER OF GRID POINTS PER ROW -
!                    THINNED GRID - IF IDIR=-1
!     IGDTMPL_THIN - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY -
!                    THINNED GRID - IF IDIR=-1. CORRESPONDS TO THE
!                    GFLD%IGDTMPL COMPONENT OF THE NCEP G2 LIBRARY
!                    GRIDMOD DATA STRUCTURE.  DEFINED ABOVE.
!     DATA_THIN    - REAL (NUMPTS_THIN,KM) THINNED GRID FIELDS IF IDIR=-1
!     IB_THIN      - INTEGER (KM) BITMAP FLAGS THINNED GRID - IF IDIR=-1
!     BITMAP_THIN  - LOGICAL (NUMPTS_THIN,KM) BITMAP FIELDS THIN GRID - IF IDIR=-1
!     IGDTMPL_FULL - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY -
!                    FULL GRID - IF IDIR=+1. CORRESPONDS TO THE
!                    GFLD%IGDTMPL COMPONENT OF THE NCEP G2 LIBRARY
!                    GRIDMOD DATA STRUCTURE.  DEFINED ABOVE.
!     DATA_FULL    - REAL (NUMPTS_FULL,KM) FULL GRID FIELDS IF IDIR=+1
!     IB_FULL      - INTEGER (KM) BITMAP FLAGS FULL GRID - IF IDIR=+1
!     BITMAP_FULL  - LOGICAL (NUMPTS_FULL,KM) BITMAP FIELDS FULL GRID - IF IDIR=+1
!     IRET         - INTEGER RETURN CODE
!                     0    SUCCESSFUL TRANSFORMATION
!                     1    IMPROPER GRID SPECIFICATION
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 IMPLICIT NONE
!
 INTEGER,               INTENT(IN   ) :: NUM_OPT
 INTEGER,               INTENT(INOUT) :: OPT_PTS(NUM_OPT)
 INTEGER,               INTENT(IN   ) :: IDIR, KM, NUMPTS_THIN, NUMPTS_FULL
 INTEGER,               INTENT(IN   ) :: IGDTLEN
 INTEGER,               INTENT(INOUT) :: IGDTMPL_THIN(IGDTLEN)
 INTEGER,               INTENT(INOUT) :: IGDTMPL_FULL(IGDTLEN)
 INTEGER,               INTENT(INOUT) :: IB_THIN(KM), IB_FULL(KM)
 INTEGER,               INTENT(  OUT) :: IRET
!
 LOGICAL(KIND=1),       INTENT(INOUT) :: BITMAP_THIN(NUMPTS_THIN,KM)
 LOGICAL(KIND=1),       INTENT(INOUT) :: BITMAP_FULL(NUMPTS_FULL,KM)
!
 REAL,                  INTENT(INOUT) :: DATA_THIN(NUMPTS_THIN,KM)
 REAL,                  INTENT(INOUT) :: DATA_FULL(NUMPTS_FULL,KM)
!
 INTEGER,               PARAMETER     :: MISSING=-1
!
 INTEGER                              :: SCAN_MODE, I, J, K, IDLAT, IDLON
 INTEGER                              :: IA, IB, IM, IM1, IM2, NPWAFS(73)
 INTEGER                              :: IS1, IS2, ISCAN, ISCALE
!
 LOGICAL                              :: TEST1, TEST2
!
 REAL                                 :: DLON, HI
 REAL                                 :: RAT1, RAT2, RLON1, RLON2
 REAL                                 :: WA, WB, X1, X2
!
 DATA npwafs/ &
       73, 73, 73, 73, 73, 73, 73, 73, 72, 72, 72, 71, 71, 71, 70,&
       70, 69, 69, 68, 67, 67, 66, 65, 65, 64, 63, 62, 61, 60, 60,&
       59, 58, 57, 56, 55, 54, 52, 51, 50, 49, 48, 47, 45, 44, 43,&
       42, 40, 39, 38, 36, 35, 33, 32, 30, 29, 28, 26, 25, 23, 22,&
       20, 19, 17, 16, 14, 12, 11,  9,  8,  6,  5,  3,  2/
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  TRANSFORM GDS
 iret=0
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  REG LAT/LON GRIDS HAVE 19 GDT ELEMENTS.
 IF (igdtlen /= 19 .OR. numpts_thin/=3447 .OR. numpts_full/=5329) THEN
   iret=1
   RETURN
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  EXPAND THINNED GDS TO FULL GDS
 IF(idir.GT.0) THEN
   scan_mode=igdtmpl_thin(19)
   iscale=igdtmpl_thin(10)*igdtmpl_thin(11)
   IF(iscale==0) iscale=10**6
   idlat=nint(1.25*float(iscale))
   test1=all(opt_pts==npwafs)
   test2=all(opt_pts==npwafs(73:1:-1))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  SOME CHECKS TO ENSURE THIS IS A WAFS GRID
   IF(scan_mode==64 .AND. igdtmpl_thin(9)==73 .AND. &
      idlat==igdtmpl_thin(18) .AND. (test1 .OR. test2) ) THEN
     igdtmpl_full=igdtmpl_thin
     im=73
     igdtmpl_full(8)=im
     rlon1=float(igdtmpl_full(13))/float(iscale)
     rlon2=float(igdtmpl_full(16))/float(iscale)
     iscan=mod(igdtmpl_full(19)/128,2)
     hi=(-1.)**iscan
     dlon=hi*(mod(hi*(rlon2-rlon1)-1+3600,360.)+1)/(im-1)
     igdtmpl_full(17)=nint(dlon*float(iscale))
   ELSE
     iret=1
   ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  CONTRACT FULL GDS TO THINNED GDS
 ELSEIF(idir.LT.0) THEN
   scan_mode=igdtmpl_full(19)
   iscale=igdtmpl_full(10)*igdtmpl_full(11)
   IF(iscale==0) iscale=10**6
   idlat=nint(1.25*float(iscale))
   idlon=idlat
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  SOME CHECKS TO ENSURE THIS IS A WAFS GRID
   IF(scan_mode==64 .AND. igdtmpl_full(8)==73 .AND. igdtmpl_full(9)==73 .AND. &
      num_opt==73 .AND. idlat==igdtmpl_full(18) .AND. idlon==igdtmpl_full(17))THEN
     igdtmpl_thin=igdtmpl_full
     igdtmpl_thin(8)=missing
     igdtmpl_thin(17)=missing
     IF(igdtmpl_thin(12)==0) THEN  ! IS LATITUDE OF ROW 1 THE EQUATOR?
       opt_pts=npwafs
     ELSE
       opt_pts=npwafs(73:1:-1)
     ENDIF
   ELSE
     iret=1
   ENDIF
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  TRANSFORM FIELDS
 IF(iret.EQ.0) THEN
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  EXPAND THINNED FIELDS TO FULL FIELDS
   IF(idir.EQ.1) THEN
     DO k=1,km
       is1=0
       is2=0
       ib_full(k)=0
       DO j=1,igdtmpl_full(9)
         im1=opt_pts(j)
         im2=igdtmpl_full(8)
         rat1=float(im1-1)/float(im2-1)
         DO i=1,im2
           x1=(i-1)*rat1+1
           ia=x1
           ia=min(max(ia,1),im1-1)
           ib=ia+1
           wa=ib-x1
           wb=x1-ia
           IF(wa.GE.wb) THEN
             IF(ib_thin(k).EQ.0.OR.bitmap_thin(is1+ia,k)) THEN
               data_full(is2+i,k)=data_thin(is1+ia,k)
               bitmap_full(is2+i,k)=.true.
             ELSE
               data_full(is2+i,k)=0.0
               bitmap_full(is2+i,k)=.false.
               ib_full(k)=1
             ENDIF
           ELSE
             IF(ib_thin(k).EQ.0.OR.bitmap_thin(is1+ib,k)) THEN
               data_full(is2+i,k)=data_thin(is1+ib,k)
               bitmap_full(is2+i,k)=.true.
             ELSE
               data_full(is2+i,k)=0.0
               bitmap_full(is2+i,k)=.false.
               ib_full(k)=1
             ENDIF
           ENDIF
         ENDDO
         is1=is1+im1
         is2=is2+im2
       ENDDO
     ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  CONTRACT FULL FIELDS TO THINNED FIELDS
   ELSEIF(idir.EQ.-1) THEN
     DO k=1,km
       is1=0
       is2=0
       ib_thin(k)=0
       DO j=1,igdtmpl_full(9)
         im1=opt_pts(j)
         im2=igdtmpl_full(8)
         rat2=float(im2-1)/float(im1-1)
         DO i=1,im1
           x2=(i-1)*rat2+1
           ia=x2
           ia=min(max(ia,1),im2-1)
           ib=ia+1
           wa=ib-x2
           wb=x2-ia
           IF(wa.GE.wb) THEN
             IF(ib_full(k).EQ.0.OR.bitmap_full(is2+ia,k)) THEN
               data_thin(is1+i,k)=data_full(is2+ia,k)
               bitmap_thin(is1+i,k)=.true.
             ELSE
               data_thin(is1+i,k)=0.0
               bitmap_thin(is1+i,k)=.false.
               ib_thin(k)=1
             ENDIF
           ELSE
             IF(ib_full(k).EQ.0.OR.bitmap_full(is2+ib,k)) THEN
               data_thin(is1+i,k)=data_full(is2+ib,k)
               bitmap_thin(is1+i,k)=.true.
             ELSE
               data_thin(is1+i,k)=0.0
               bitmap_thin(is1+i,k)=.false.
               ib_thin(k)=1
             ENDIF
           ENDIF
         ENDDO
         is1=is1+im1
         is2=is2+im2
       ENDDO
     ENDDO
   ENDIF
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 END SUBROUTINE ipxwafs3