ipxwafs.F90

Go to the documentation of this file.

!> @file
!> @brief Expand or contract wafs grids.
!> @author Iredell @date 96-04-10
 
!> Expand or contract wafs grids.
!>
!> 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 linear interpolation and is not
!> reversible.
!>
!> This routine does not work for bitmapped data.
!>
!> ### Program History Log
!> Date | Programmer | Comments
!> -----|------------|---------
!> 96-04-10 | iredell | initial version
!> 2015-july | gayno | convert to grib 2
!>
!> @param[in] idir integer transform option
!> - (+1 to expand thinned fields to full fields)
!> - (-1 to contract full fields to thinned fields)
!> @param[in] numpts_thin integer number of grid points - thinned
!> grid.  must be 3447.
!> @param[in] numpts_full integer number of grid points - full
!> grid. must be 5329 (73^2).
!> @param[in] km integer number of fields to transform
!> @param[in] num_opt integer number of values to describe the thinned
!> grid. must be 73.  dimension of array opt_pts.
!> @param[inout] opt_pts integer (num_opt) number of grid points per
!> row - thinned grid - if idir=+1
!> @param[in] 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.
!> @param[inout] 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
!> @param[inout] data_thin real (numpts_thin,km) thinned grid fields if idir=+1
 
!> @param[inout] 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
!> @param[inout] data_full real (numpts_full,km) full grid fields if idir=-1
!> @param[out] iret integer return code
!> - 0 successful transformation
!> - 1 improper grid specification
!>
!> @author Iredell @date 96-04-10
 SUBROUTINE ipxwafs(IDIR, NUMPTS_THIN, NUMPTS_FULL, KM, NUM_OPT, &
                    OPT_PTS, IGDTLEN, IGDTMPL_THIN, DATA_THIN,  &
                    IGDTMPL_FULL, DATA_FULL, IRET)
 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(  OUT) :: IRET
!
 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
       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=int(x1)
           ia=min(max(ia,1),im1-1)
           ib=ia+1
           wa=ib-x1
           wb=x1-ia
           data_full(is2+i,k)=wa*data_thin(is1+ia,k)+wb*data_thin(is1+ib,k)
         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
       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=int(x2)
           ia=min(max(ia,1),im2-1)
           ib=ia+1
           wa=ib-x2
           wb=x2-ia
           data_thin(is1+i,k)=wa*data_full(is2+ia,k)+wb*data_full(is2+ib,k)
         ENDDO
         is1=is1+im1
         is2=is2+im2
       ENDDO
     ENDDO
   ENDIF
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 SUBROUTINE ipxwafs(IDIR, NUMPTS_THIN, NUMPTS_FULL, KM, NUM_OPT, & …
 END SUBROUTINE ipxwafs