add switch for two different init_grib2 approaches for non-E arakawa grids

src/ip_grid_mod.F90

@@ -28,10 +28,14 @@ module ip_grid_mod
   integer, public, parameter :: ROT_EQUID_CYLIND_E_GRID_ID_GRIB2 = 32768 !< Integer grid number for rotated equidistant cylindrical E-stagger grid (grib2)
   integer, public, parameter :: ROT_EQUID_CYLIND_B_GRID_ID_GRIB2 = 32769 !< Integer grid number for rotated equidistant cylindrical B-stagger grid (grib2)
 
+  logical, public, save :: ncep_post_arakawa=.false. !< Use ncep_post/wgrib2-compatible version of init_grib2() for non-E Arakawa grids (enable with use_ncep_post_arakawa())
+
   private
   public :: ip_grid
   public :: gdswzd_interface
   public :: operator(==)
+  public :: use_ncep_post_arakawa
+  public :: unuse_ncep_post_arakawa
 
   !> Abstract grid that holds fields and methods common to all grids.
   !! ip_grid is meant to be subclassed when implementing a new grid.
@@ -172,6 +176,26 @@ end subroutine init_grib2_interface
 
 contains
 
+  !> Enables ncep_post/wgrib2-compatible non-E Arakawa grib2 grids
+  !> by setting 'ncep_post_arakawa=.true.'.
+  !> This subroutine should be called prior to init_grib2().
+  !>
+  !> @author Alex Richert
+  !> @date May 2024
+  subroutine use_ncep_post_arakawa() bind(c)
+    ncep_post_arakawa = .true.
+  end subroutine use_ncep_post_arakawa
+
+  !> Disables ncep_post/wgrib2-compatible non-E Arakawa grib2 grids
+  !> by setting 'ncep_post_arakawa=.false.'.
+  !> This subroutine should be called prior to init_grib2().
+  !>
+  !> @author Alex Richert
+  !> @date May 2024
+  subroutine unuse_ncep_post_arakawa() bind(c)
+    ncep_post_arakawa = .false.
+  end subroutine unuse_ncep_post_arakawa
+
   !> Compares two grids.
   !>
   !> @param[in] grid1 An ip_grid

src/ip_rot_equid_cylind_grid_mod.F90

@@ -130,16 +130,37 @@ subroutine init_grib1(self, g1_desc)
   end subroutine init_grib1
 
   !> Initializes a Rotated equidistant cylindrical grid given a
-  !> grib2_descriptor object.
+  !> grib2_descriptor object. Call 'use_ncep_post_arakawa()' before
+  !> using this subroutine to use ncep_post-compatible grid definition.
   !>
   !> @param[inout] self The grid to initialize
   !> @param[in] g2_desc A grib2_descriptor
   !>
-  !> @author Gayno @date 2007-NOV-15
+  !> @author Alex Richert @date 2024-MAY-20
   subroutine init_grib2(self, g2_desc)
     class(ip_rot_equid_cylind_grid), intent(inout) :: self
     type(grib2_descriptor), intent(in) :: g2_desc
 
+    if (ncep_post_arakawa.and.(g2_desc%gdt_num.eq.32769)) then
+      call init_grib2_ncep_post(self, g2_desc)
+    else
+      call init_grib2_default(self, g2_desc)
+    endif
+
+  end subroutine init_grib2
+
+  !> Initializes a Rotated equidistant cylindrical grid given a
+  !> grib2_descriptor object. Uses template definitions consistent
+  !> with WMO standards.
+  !>
+  !> @param[inout] self The grid to initialize
+  !> @param[in] g2_desc A grib2_descriptor
+  !>
+  !> @author Gayno @date 2007-NOV-15
+  subroutine init_grib2_default(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
@@ -192,7 +213,92 @@ subroutine init_grib2(self, g2_desc)
       self%nscan = mod(igdtmpl(19) / 32, 2)
       self%nscan_field_pos = self%nscan
     end associate
-  end subroutine init_grib2
+  end subroutine init_grib2_default
+
+  !> Initializes a Rotated equidistant cylindrical grid given a
+  !> grib2_descriptor object. Uses template number definitions in
+  !> a manner compatible with wgrib2 and ncep_post, which is based
+  !> on grib1.
+  !>
+  !> @param[inout] self The grid to initialize
+  !> @param[in] g2_desc A grib2_descriptor
+  !>
+  !> @author Alex Richert, George Gayno @date 2024-MAY-20
+  subroutine init_grib2_ncep_post(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
+    real(kd) :: hs, hs2, slat1, slat2, slatr, clon1, clon2, clat1, clat2, clatr, clonr, rlonr, rlatr
+
+    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(15))/FLOAT(ISCALE)
+
+      self%RLON0=FLOAT(IGDTMPL(16))/FLOAT(ISCALE)
+
+      RLAT2=FLOAT(IGDTMPL(20))/FLOAT(ISCALE)
+      RLON2=FLOAT(IGDTMPL(21))/FLOAT(ISCALE)
+
+      self%IROT=MOD(IGDTMPL(14)/8,2)
+      self%IM=IGDTMPL(8)
+      self%JM=IGDTMPL(9)
+
+      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
+      IF (self%WBD > 180.0) self%WBD = self%WBD - 360.0
+      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)
+
+      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_ncep_post
 
   !> GDS wizard for rotated equidistant cylindrical.
   !>

tests/CMakeLists.txt

@@ -100,6 +100,7 @@ foreach(kind ${kinds})
   add_test(test_rotatedB_spectral_vector_grib2_${kind} test_vector_grib2_${kind} 205 4)
   add_test(test_rotatedE_budget_vector_grib2_${kind} test_vector_grib2_${kind} 203 3)
   add_test(test_rotatedB_direct_spectral_vector_grib2_${kind} test_vector_grib2_${kind} 32769 4)
+  add_test(test_rotatedB_direct_ncep_post_spectral_vector_grib2_${kind} test_vector_grib2_${kind} 32769b 4)
   add_test(test_rotatedE_direct_budget_vector_grib2_${kind} test_vector_grib2_${kind} 32768 3)
 
   # vector station point tests

tests/interp_mod_grib2.F90

@@ -3,6 +3,7 @@
 ! Kyle Gerheiser June, 2021
 
   use ip_mod
+  use ip_grid_mod
   implicit none
 
 contains
@@ -547,13 +548,14 @@ subroutine interp_vector(grid, interp_opt)
        output_gdtmpl=gdtmpl203
        i_output = output_gdtmpl(8)
        j_output = output_gdtmpl(9)
-    case ('32769')
+    case ('32769','32769b')
        output_gdtnum = 32769
        output_gdtlen = gdtlen205
        allocate(output_gdtmpl(output_gdtlen))
        output_gdtmpl=gdtmpl205
        i_output = output_gdtmpl(8)
        j_output = output_gdtmpl(9)
+       if (trim(grid).eq.'32769b') call use_ncep_post_arakawa()
     case default
        print*,"ERROR: ENTER VALID GRID NUMBER."
        stop 55