diff --git a/src/ip_grid_mod.F90 b/src/ip_grid_mod.F90 index 07f3b20c..c428e1a2 100644 --- a/src/ip_grid_mod.F90 +++ b/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 diff --git a/src/ip_rot_equid_cylind_grid_mod.F90 b/src/ip_rot_equid_cylind_grid_mod.F90 index 88b6cbec..e3fcb488 100644 --- a/src/ip_rot_equid_cylind_grid_mod.F90 +++ b/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. !> diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ab5236c2..93553b8f 100644 --- a/tests/CMakeLists.txt +++ b/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 diff --git a/tests/data/baseline_data/grib2/vector/4_byte_bin/grid32769b.opt4.bin_4 b/tests/data/baseline_data/grib2/vector/4_byte_bin/grid32769b.opt4.bin_4 new file mode 100644 index 00000000..47f67eff Binary files /dev/null and b/tests/data/baseline_data/grib2/vector/4_byte_bin/grid32769b.opt4.bin_4 differ diff --git a/tests/data/baseline_data/grib2/vector/8_byte_bin/grid32769b.opt4.bin_8 b/tests/data/baseline_data/grib2/vector/8_byte_bin/grid32769b.opt4.bin_8 new file mode 100644 index 00000000..f8c469ac Binary files /dev/null and b/tests/data/baseline_data/grib2/vector/8_byte_bin/grid32769b.opt4.bin_8 differ diff --git a/tests/interp_mod_grib2.F90 b/tests/interp_mod_grib2.F90 index cef88be6..15f77a35 100644 --- a/tests/interp_mod_grib2.F90 +++ b/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