ip_gaussian_grid_mod.F90

Go to the documentation of this file.

!> @file
!! @brief Gaussian grid coordinate transformations.
!! @author Mark Iredell, George Gayno, Kyle Gerheiser
!! @date July 2021
 
!> @brief Gaussian grid coordinate transformations.
!!
!! Octet numbers refer to [GRIB2 - GRID DEFINITION TEMPLATE 3.40
!! Gaussian
!! Latitude/Longitude](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp3-40.shtml).
!!
!! @author George Gayno, Mark Iredell, Kyle Gerheiser
!! @date July 2021
module ip_gaussian_grid_mod
  use ip_grid_descriptor_mod
  use ip_grid_mod
  use earth_radius_mod
  use ip_constants_mod
  use sp_mod
  implicit none
 
  private
  public :: ip_gaussian_grid
 
  type, extends(ip_grid) :: ip_gaussian_grid
     integer :: jh !< Scan mode flag in 'j' direction. When '1' points scan from N to S. When "-1" points scan from S to N.
     real :: dlon !< "i"-direction increment. GRIB2 Section 3, octets 64-67.
     real :: rlat1 !<  Latitude of first grid point. GRIB2 Section 3, octets 47-50.
     real :: rlon1 !< Longitude of first grid point. GRIB2 Section 3, octets 51-54.
     real :: rlon2 !< Longitude of last grid point. GRIB2 Section 3, octets 60-63.
     real :: hi !< Scan mode flag in 'i' direction. When '1' points scan from W to E. When "-1" points scan from E to W.
     integer :: jg !< Number of parallels between the equator and pole times 2. GRIB2 Section 3, octets 68-71.
     integer :: jscan !< Scanning mode in the 'j' direction. GRIB2 Section 3, octet 72.
   contains
     !> Initializes a gaussian grid given a grib1_descriptor object.
     !> @return N/A @copydoc ip_gaussian_grid_mod::init_grib1
     procedure :: init_grib1
     !> Initializes a gaussian grid given a grib2_descriptor object.
     !> @return N/A @copydoc ip_gaussian_grid_mod::init_grib2
     procedure :: init_grib2
     !> Calculates Earth coordinates (iopt = 1) or grid coorindates (iopt = -1)
     !> for Gaussian grids.
     !> @return N/A @copydoc ip_gaussian_grid_mod::gdswzd_gaussian
     procedure :: gdswzd => gdswzd_gaussian
  end type ip_gaussian_grid
  type, extends(ip_grid) :: ip_gaussian_grid …
 
  INTEGER :: j1 !< 'j' index of first grid point within the global array of latitudes.
  INTEGER :: jh !< Scan mode flag in 'j' direction. When '1' points scan from N to S. When "-1" points scan from S to N.
  REAL, ALLOCATABLE :: blat(:) !< Gaussian latitude for each parallel.
  REAL :: dlon !< "i"-direction increment. GRIB2 Section 3, octets 64-67.
  REAL :: rerth !< Radius of the earth. GRIB2 Section 3, octets 15-30.
  REAL, ALLOCATABLE :: ylat_row(:) !< dy/dlat for each row in 1/degrees.
 
contains
 
  !> Initializes a gaussian grid given a grib1_descriptor object.
  !>
  !> @param[inout] self The grid to initialize
  !> @param[in] g1_desc A grib1_descriptor
  !>
  !> @author Kyle Gerheiser
  !> @date July 2021
  subroutine init_grib1(self, g1_desc)
    class(ip_gaussian_grid), intent(inout) :: self
    type(grib1_descriptor), intent(in) :: g1_desc
 
    integer :: iscan, jg
 
    associate(kgds => g1_desc%gds)
      self%rerth = 6.3712e6
      self%eccen_squared = 0.0
 
      self%IM=kgds(2)
      self%JM=kgds(3)
      self%RLAT1=kgds(4)*1.e-3
      self%RLON1=kgds(5)*1.e-3
      self%RLON2=kgds(8)*1.e-3
      self%JG=kgds(10)*2
      iscan=mod(kgds(11)/128,2)
      self%JSCAN=mod(kgds(11)/64,2)
      self%HI=(-1.)**iscan
      self%JH=(-1)**self%JSCAN
      self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(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
 
      self%iwrap=nint(360 / abs(self%dlon))
      if(self%im < self%iwrap) self%iwrap = 0
 
      if(self%iwrap > 0 .and. mod(self%iwrap, 2) == 0) then
         jg=kgds(10)*2
         if(self%jm == self%jg) then
            self%jwrap1 = 1
            self%jwrap2 = 2 * self%jm + 1
         endif
      endif
 
    end associate
  subroutine init_grib1(self, g1_desc) …
  end subroutine init_grib1
 
  !> Initializes a gaussian grid given a grib2_descriptor object.
  !> @param[inout] self The grid to initialize
  !> @param[in] g2_desc A grib2_descriptor
  !>
  !> @author Kyle Gerheiser
  !> @date July 2021
  subroutine init_grib2(self, g2_desc)
    class(ip_gaussian_grid), intent(inout) :: self
    type(grib2_descriptor), intent(in) :: g2_desc
 
    integer :: iscale, iscan, jg
 
    associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
      call earth_radius(igdtmpl, igdtlen, self%rerth, self%eccen_squared)
 
      self%IM=igdtmpl(8)
      self%JM=igdtmpl(9)
      iscale=igdtmpl(10)*igdtmpl(11)
      IF(iscale==0) iscale=10**6
      self%RLAT1=float(igdtmpl(12))/float(iscale)
      self%RLON1=float(igdtmpl(13))/float(iscale)
      self%RLON2=float(igdtmpl(16))/float(iscale)
      self%JG=igdtmpl(18)*2
      iscan=mod(igdtmpl(19)/128,2)
      self%JSCAN=mod(igdtmpl(19)/64,2)
      self%HI=(-1.)**iscan
      self%JH=(-1)**self%JSCAN
      self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
 
 
      self%iwrap = nint(360 / abs(self%dlon))
      if(self%im < self%iwrap) self%iwrap = 0
      self%jwrap1 = 0
      self%jwrap2 = 0
      if(self%iwrap > 0 .and. mod(self%iwrap, 2) == 0) then
         jg = igdtmpl(18) * 2
         if(self%jm == jg) then
            self%jwrap1=1
            self%jwrap2 = 2 * self%jm + 1
         endif
      endif
      self%nscan = mod(igdtmpl(19) / 32, 2)
      self%nscan_field_pos = self%nscan
      self%kscan = 0
    end associate
 
  subroutine init_grib2(self, g2_desc) …
  end subroutine init_grib2
 
  !> Calculates Earth coordinates (iopt = 1) or grid coorindates (iopt = -1)
  !> for Gaussian grids.
  !>
  !> 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.
  !>
  !> @param[in] self The grid object gdswzd was called on.
  !> @param[in] iopt option flag
  !>            - +1 to compute earth coords of selected grid coords.
  !>            - -1 o compute grid coords of selected earth coords.
  !> @param[in] npts Maximum number of coordinates.
  !> @param[in] fill Fill value to set invalid output data.
  !>            Must be impossible value; suggested value: -9999.
  !> @param[inout] xpts Grid x point coordinates if iopt>0.
  !> @param[inout] ypts Grid y point coordinates if iopt>0.
  !> @param[inout] rlon Earth longitudes in degrees e if iopt<0
  !>                   (Acceptable range: -360. to 360.)
  !> @param[inout] rlat Earth latitudes in degrees n if iopt<0
  !>                (Acceptable range: -90. to 90.)
  !> @param[out] nret Number of valid points computed.
  !> @param[out] crot Optional clockwise vector rotation cosines.
  !> @param[out] srot Optional clockwise vector rotation sines.
  !> @param[out] xlon Optional dx/dlon in 1/degrees.
  !> @param[out] xlat Optional dx/dlat in 1/degrees.
  !> @param[out] ylon Optional dy/dlon in 1/degrees.
  !> @param[out] ylat Optional dy/dlat in 1/degrees.
  !> @param[out] area Optional area weights in m**2.
  !>
  !> @author Mark Iredell, George Gayno, Kyle Gerheiser
  !> @date July 2021
  SUBROUTINE gdswzd_gaussian(self,IOPT,NPTS,FILL, &
       XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    IMPLICIT NONE
    !
    class(ip_gaussian_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                        :: JSCAN, IM, JM
    INTEGER                        :: J, JA, JG
    INTEGER                        :: N
    !
    LOGICAL                        :: LROT, LMAP, LAREA
    !
    REAL,            ALLOCATABLE   :: ALAT(:), ALAT_JSCAN(:)
    REAL,            ALLOCATABLE   :: ALAT_TEMP(:),BLAT_TEMP(:)
    REAL                           :: HI, RLATA, RLATB, RLAT1, RLON1, RLON2
    REAL                           :: XMAX, XMIN, YMAX, YMIN, YPTSA, YPTSB
    REAL                           :: WB
 
    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
 
    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
 
    im=self%im
    jm=self%jm
 
    rlat1=self%rlat1
    rlon1=self%rlon1
    rlon2=self%rlon2
 
    jg=self%jg
    jscan=self%jscan
    hi=self%hi
 
    jh=self%jh
    dlon=self%dlon
    rerth = self%rerth
 
    ALLOCATE(alat_temp(jg))
    ALLOCATE(blat_temp(jg))
    CALL splat(4,jg,alat_temp,blat_temp)
    ALLOCATE(alat(0:jg+1))
    ALLOCATE(blat(0:jg+1))
    !$OMP PARALLEL DO PRIVATE(JA) SCHEDULE(STATIC)
    DO ja=1,jg
       alat(ja)=real(dpr*asin(alat_temp(ja)))
       blat(ja)=blat_temp(ja)
    ENDDO
    !$OMP END PARALLEL DO
    DEALLOCATE(alat_temp,blat_temp)
    alat(0)=180.-alat(1)
    alat(jg+1)=-alat(0)
    blat(0)=-blat(1)
    blat(jg+1)=blat(0)
    j1=1
    DO WHILE(j1.LT.jg.AND.rlat1.LT.(alat(j1)+alat(j1+1))/2)
       j1=j1+1
    ENDDO
    IF(lmap)THEN
       ALLOCATE(alat_jscan(jg))
       DO ja=1,jg
          alat_jscan(j1+jh*(ja-1))=alat(ja)
       ENDDO
       ALLOCATE(ylat_row(0:jg+1))
       DO ja=2,(jg-1)
          ylat_row(ja)=2.0/(alat_jscan(ja+1)-alat_jscan(ja-1))
       ENDDO
       ylat_row(1)=1.0/(alat_jscan(2)-alat_jscan(1))
       ylat_row(0)=ylat_row(1)
       ylat_row(jg)=1.0/(alat_jscan(jg)-alat_jscan(jg-1))
       ylat_row(jg+1)=ylat_row(jg)
       DEALLOCATE(alat_jscan)
    ENDIF
    xmin=0
    xmax=im+1
    IF(im.EQ.nint(360/abs(dlon))) xmax=im+2
    ymin=0.5
    ymax=jm+0.5
    nret=0
 
    !  TRANSLATE GRID COORDINATES TO EARTH COORDINATES
    IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
       !$OMP PARALLEL DO PRIVATE(N,J,WB,RLATA,RLATB) 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
             rlon(n)=mod(rlon1+dlon*(xpts(n)-1)+3600,360.)
             j=int(ypts(n))
             wb=ypts(n)-j
             rlata=alat(j1+jh*(j-1))
             rlatb=alat(j1+jh*j)
             rlat(n)=rlata+wb*(rlatb-rlata)
             nret=nret+1
             IF(lrot) CALL gaussian_vect_rot(crot(n),srot(n))
             IF(lmap) CALL gaussian_map_jacob(ypts(n),&
                  xlon(n),xlat(n),ylon(n),ylat(n))
             IF(larea) CALL gaussian_grid_area(ypts(n),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,JA,YPTSA, YPTSB, WB) REDUCTION(+:NRET) SCHEDULE(STATIC)
       DO n=1,npts
          xpts(n)=fill
          ypts(n)=fill
          IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90) THEN
             xpts(n)=1+hi*mod(hi*(rlon(n)-rlon1)+3600,360.)/dlon
             ja=min(int((jg+1)/180.*(90-rlat(n))),jg)
             IF(rlat(n).GT.alat(ja)) ja=max(ja-2,0)
             IF(rlat(n).LT.alat(ja+1)) ja=min(ja+2,jg)
             IF(rlat(n).GT.alat(ja)) ja=ja-1
             IF(rlat(n).LT.alat(ja+1)) ja=ja+1
             yptsa=1+jh*(ja-j1)
             yptsb=1+jh*(ja+1-j1)
             wb=(alat(ja)-rlat(n))/(alat(ja)-alat(ja+1))
             ypts(n)=yptsa+wb*(yptsb-yptsa)
             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 gaussian_vect_rot(crot(n),srot(n))
                IF(lmap) CALL gaussian_map_jacob(ypts(n), &
                     xlon(n),xlat(n),ylon(n),ylat(n))
                IF(larea) CALL gaussian_grid_area(ypts(n),area(n))
             ELSE
                xpts(n)=fill
                ypts(n)=fill
             ENDIF
          ENDIF
       ENDDO
       !$OMP END PARALLEL DO
    ENDIF
    DEALLOCATE(alat, blat)
    IF (ALLOCATED(ylat_row)) DEALLOCATE(ylat_row)
 
  SUBROUTINE gdswzd_gaussian(self,IOPT,NPTS,FILL, & …
  END SUBROUTINE gdswzd_gaussian
 
  !> Computes the vector rotation sines and cosines for a gaussian
  !> cylindrical grid.
  !>
  !> @param[out] crot Clockwise vector rotation cosines.
  !> @param[out] srot Clockwise vector rotation sines.
  !>
  !> @note
  !> ugrid=crot*uearth-srot*vearth;
  !> vgrid=srot*uearth+crot*vearth)
  !>
  !> @author George Gayno
  !> @date July 2021
  SUBROUTINE gaussian_vect_rot(CROT,SROT)
    IMPLICIT NONE
 
    REAL,                INTENT(  OUT) :: CROT, SROT
 
    crot=1.0
    srot=0.0
 
  SUBROUTINE gaussian_vect_rot(CROT,SROT) …
  END SUBROUTINE gaussian_vect_rot
 
  !> Computes the map jacobians for a gaussian cylindrical grid.
  !>
  !> @param[in] ypts y-index of grid point.
  !> @param[out] xlon dx/dlon in 1/degrees.
  !> @param[out] xlat dx/dlat in 1/degrees.
  !> @param[out] ylon dy/dlon in 1/degrees.
  !> @param[out] ylat dy/dlat in 1/degrees.
  !>
  !> @author George Gayno
  !> @date July 2021
  SUBROUTINE gaussian_map_jacob(YPTS, XLON, XLAT, YLON, YLAT)
    IMPLICIT NONE
 
    REAL,                INTENT(IN   ) :: YPTS
    REAL,                INTENT(  OUT) :: XLON, XLAT, YLON, YLAT
 
    xlon=1/dlon
    xlat=0.
    ylon=0.
    ylat=ylat_row(nint(ypts))
 
  SUBROUTINE gaussian_map_jacob(YPTS, XLON, XLAT, YLON, YLAT) …
  END SUBROUTINE gaussian_map_jacob
 
  !> Computes the grid box area for a gaussian cylindrical grid.
  !>
  !> @param[in] ypts y-index of grid point.
  !> @param[out] area Area weights in m^2
  !>
  !> @author Mark Iredell, George Gayno
  !> @date July 2021
  SUBROUTINE gaussian_grid_area(YPTS,AREA)
    IMPLICIT NONE
 
    REAL,            INTENT(IN   ) :: YPTS
    REAL,            INTENT(  OUT) :: AREA
 
    INTEGER                        :: J
 
    REAL                           :: WB, WLAT, WLATA, WLATB
 
    j = int(ypts)
    wb=ypts-j
    wlata=blat(j1+jh*(j-1))
    wlatb=blat(j1+jh*j)
    wlat=wlata+wb*(wlatb-wlata)
    area=real(rerth**2*wlat*dlon/dpr)
 
  SUBROUTINE gaussian_grid_area(YPTS,AREA) …
  END SUBROUTINE gaussian_grid_area
end module ip_gaussian_grid_mod