chgres_cube: Implement WMO grib2 template 1 (rotated lat-lon) read ca…

…pability (#902)

Add ability to read GRIB2 data that uses the WMO-standard rotated lat-lon grid
template (GRIB2 grid template 1). Update 'readthedocs'.

Fixes #901.

docs/source/chgres_cube.rst

@@ -369,7 +369,7 @@ Namelist variables with “input” in their name refer to data input to chgres_
               * Set to 2 to create a boundary condition file. Use this option for all but the initialization time.
       * halo_blend - Integer number of row/columns to apply halo blending into the domain, where model and lateral boundary tendencies are applied.
       * halo_bndy - Integer number of rows/columns that exist within the halo, where pure lateral boundary conditions are applied.
-      * external_model - Name of source model for input data. Valid options: 'GFS', 'NAM', 'RAP', 'HRRR'. (Default: 'GFS')
+      * external_model - Name of source model for input data. Valid options: 'GFS', 'NAM', 'RAP', 'HRRR', 'RRFS'. (Default: 'GFS')
 
 **Optional Entries**
 

sorc/chgres_cube.fd/atm_input_data.F90

@@ -3018,6 +3018,14 @@ subroutine read_winds(u,v,localpet,octet_23,rlevs,lugb,pdt_num)
      print*, "- CALL CALCALPHA_ROTLATLON with center lat,lon = ",latin1,lov
      call calcalpha_rotlatlon(lat,lon,latin1,lov,alpha)
 
+   elseif (gfld%igdtnum == 1) then ! grid definition template number - non-E stagger rotated lat/lon grid
+
+     latin1 = real(float(gfld%igdtmpl(20))/1.0E6, kind=esmf_kind_r4) + 90.0_esmf_kind_r4
+     lov = real(float(gfld%igdtmpl(21))/1.0E6, kind=esmf_kind_r4)
+
+     print*, "- CALL CALCALPHA_ROTLATLON with center lat,lon = ",latin1,lov
+     call calcalpha_rotlatlon(lat,lon,latin1,lov,alpha) 
+
    elseif (gfld%igdtnum == 30) then ! grid definition template number - lambert conformal grid.
 
      lov = real(float(gfld%igdtmpl(14))/1.0E6, kind=esmf_kind_r4)
@@ -3087,7 +3095,7 @@ subroutine read_winds(u,v,localpet,octet_23,rlevs,lugb,pdt_num)
           u(:,:,vlev) = u_tmp
           v(:,:,vlev) = v_tmp
         endif
-      else if (gfld%igdtnum == 32769) then ! grid definition template number - rotated lat/lon grid
+      else if (gfld%igdtnum == 32769 .or. gfld%igdtnum == 1) then ! grid definition template number - rotated lat/lon grid
         ws = sqrt(u_tmp**2 + v_tmp**2)
         wd = real((atan2(-u_tmp,-v_tmp) / d2r), kind=esmf_kind_r4) ! calculate grid-relative wind direction
         wd = real((wd + alpha + 180.0), kind=esmf_kind_r4) ! Rotate from grid- to earth-relative direction

sorc/chgres_cube.fd/model_grid.F90

@@ -675,7 +675,7 @@ subroutine define_input_grid_grib2(npets)
  elseif (gfld%igdtnum == 30) then
    print*,"- INPUT DATA ON LAMBERT CONFORMAL GRID."
    input_grid_type = 'lambert'
- elseif (gfld%igdtnum == 32769) then
+ elseif (gfld%igdtnum == 32769 .or. gfld%igdtnum == 1) then
    print*,"- INPUT DATA ON ROTATED LAT/LON GRID."
    input_grid_type = 'rotated_latlon'
  else
@@ -1477,10 +1477,14 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
 
  integer, intent(in   )  :: igdtnum, igdtlen, igdstmpl(igdtlen)
  integer, intent(  out)  :: kgds(200), ni, nj
- integer                 :: iscale
+ integer                 :: iscale, i
 
  real,    intent(  out)  :: res
 
+ real                    :: clatr, slatr, clonr, dpr, slat
+ real                    :: slat0, clat0, clat, clon, rlat
+ real                    :: rlon0, rlon, hs
+
  kgds=0
 
  if (igdtnum.eq.32769) then        ! rot lat/lon b grid
@@ -1528,6 +1532,112 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
      res = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) &
              * 0.5 * 111.0
 
+   elseif (igdtnum.eq.1) then        ! rot lat/lon b grid using standard wmo
+                                    ! template.
+
+     iscale=igdstmpl(10)*igdstmpl(11)
+     if (iscale == 0) iscale = 1e6
+     kgds(1)=205                    ! oct 6,     rotated lat/lon for Non-E
+                                    !            Stagger grid
+     kgds(2)=igdstmpl(8)            ! octs 7-8,  Ni
+     ni = kgds(2)
+     kgds(3)=igdstmpl(9)            ! octs 9-10, Nj
+     nj = kgds(3)
+
+     kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.)  ! octs 11-13, Lat of
+                                                            ! 1st grid point
+     kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.)  ! octs 14-16, Lon of
+                                                            ! 1st grid point
+
+     kgds(6)=0                      ! oct 17, resolution and component flags
+     if (igdstmpl(1)==2 ) kgds(6)=64
+     if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) )  kgds(6)=kgds(6)+128
+     if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
+
+     kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.)       ! octs 18-20,
+                                                                 ! Lat of cent of rotation
+     kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.)       ! octs 21-23,
+                                                                 ! Lon of cent of rotation
+     kgds(7) = kgds(7) + 90000.0
+     print*, "INPUT LAT, LON CENTER ", kgds(7), kgds(8)
+
+     DPR = 180.0/3.1415926
+     CLATR=COS((float(kgds(4))/1000.0)/DPR)
+     SLATR=SIN((float(kgds(4))/1000.0)/DPR)
+     CLONR=COS((float(kgds(5))/1000.0)/DPR)
+     SLAT0=SIN((float(kgds(7))/1000.0)/DPR)
+     CLAT0=COS((float(kgds(7))/1000.0)/DPR)
+
+     SLAT=CLAT0*SLATR+SLAT0*CLATR*CLONR
+     CLAT=SQRT(1-SLAT**2)
+     CLON=(CLAT0*CLATR*CLONR-SLAT0*SLATR)/CLAT
+     CLON=MIN(MAX(CLON,-1.0),1.0)
+
+     RLAT=DPR*ASIN(SLAT)
+
+     RLON0=float(kgds(8))/1000.0
+
+     if ((kgds(5)-kgds(8)) > 0) then
+       HS = -1.0
+     else
+       HS = 1.0
+     endif
+
+     RLON = MOD(RLON0+HS*DPR*ACOS(CLON)+3600,360.0)
+
+     kgds(4)=nint(rlat*1000.)  ! octs 11-13, Lat of
+     kgds(5)=nint(rlon*1000.)  ! octs 14-16, Lon of
+
+     kgds(12)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 29-31, Lat of
+                                                            ! last grid point
+     kgds(13)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 32-34, Lon of
+                                                            ! last grid point
+
+     CLATR=COS((float(kgds(12))/1000.0)/DPR)
+     SLATR=SIN((float(kgds(12))/1000.0)/DPR)
+     CLONR=COS((float(kgds(13))/1000.0)/DPR)
+
+     SLAT=CLAT0*SLATR+SLAT0*CLATR*CLONR
+     RLAT=DPR*ASIN(SLAT)
+
+     CLAT=SQRT(1-SLAT**2)
+     CLON=(CLAT0*CLATR*CLONR-SLAT0*SLATR)/CLAT
+     CLON=MIN(MAX(CLON,-1.0),1.0)
+
+     if ((kgds(13)-kgds(8)) > 0) then
+       HS = -1.0
+     else
+       HS = 1.0
+     endif
+
+     RLON = MOD(RLON0+HS*DPR*ACOS(CLON)+3600,360.0)
+
+     print*,'got here last point ',kgds(12), kgds(13)
+     print*,'got here last point rotated ', rlat, rlon
+
+     kgds(12)=nint(rlat*1000.)  ! octs 11-13, Lat of
+     kgds(13)=nint(rlon*1000.)  ! octs 14-16, Lon of
+
+     kgds(9)=igdstmpl(17)
+     kgds(10)=igdstmpl(18)
+
+     kgds(11) = 0                   ! oct 28, scan mode
+     if (btest(igdstmpl(19),7)) kgds(11) = 128
+     if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) +  64
+     if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) +  32
+
+     kgds(19)=0    ! oct 4, # vert coordinate parameters
+     kgds(20)=255  ! oct 5, used for thinned grids, set to 255
+
+     res = ((float(kgds(9)) / 1.e6) + (float(kgds(10)) / 1.e6)) &
+             * 0.5 * 111.0
+
+
+     do i = 1, 25
+       print*,'final kgds ',i,kgds(i)
+     enddo
+
+
    elseif(igdtnum==30) then
 
      kgds(1)=3                      ! oct 6,     lambert conformal

sorc/chgres_cube.fd/program_setup.F90

@@ -54,7 +54,7 @@ module program_setup
 !!                                    gaussian nemsio files
 !!                                 - "gfs_sigio" for spectral gfs
 !!                                    gfs sigio/sfcio files.
- character(len=20),  public      :: external_model="GFS"  !< The model that the input data is derived from. Current supported options are: "GFS", "HRRR", "NAM", "RAP". Default: "GFS"
+ character(len=20),  public      :: external_model="GFS"  !< The model that the input data is derived from. Current supported options are: "GFS", "HRRR", "NAM", "RAP", "RRFS". Default: "GFS"
 
  integer, parameter, public      :: max_tracers=100 !< Maximum number of atmospheric tracers processed.
  integer, public                 :: num_tracers !< Number of atmospheric tracers to be processed.
@@ -318,8 +318,8 @@ subroutine read_setup_namelist(filename)
 !-------------------------------------------------------------------------
 
  if (trim(input_type) == "grib2") then
-   if (.not. any((/character(4)::"GFS","NAM","RAP","HRRR"/)==trim(external_model))) then
-     call error_handler( "KNOWN SUPPORTED external_model INPUTS ARE GFS, NAM, RAP, AND HRRR. " // &
+   if (.not. any((/character(4)::"GFS","NAM","RAP","HRRR","RRFS"/)==trim(external_model))) then
+     call error_handler( "KNOWN SUPPORTED external_model INPUTS ARE GFS, NAM, RAP, HRRR, AND RRFS. " // &
     "IF YOU WISH TO PROCESS GRIB2 DATA FROM ANOTHER MODEL, YOU MAY ATTEMPT TO DO SO AT YOUR OWN RISK. " // &
     "ONE WAY TO DO THIS IS PROVIDE NAM FOR external_model AS IT IS A RELATIVELY STRAIGHT-" // &
     "FORWARD REGIONAL GRIB2 FILE. YOU MAY ALSO COMMENT OUT THIS ERROR MESSAGE IN " // &