Skip to content

Commit

Permalink
chgres_cube: Implement WMO grib2 template 1 (rotated lat-lon) read ca…
Browse files Browse the repository at this point in the history 
…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.
  • Loading branch information
@LarissaReames-NOAA
LarissaReames-NOAA committed Feb 27, 2024
1 parent 415f616 commit 1dac855
Show file tree
Hide file tree
Showing 4 changed files with 125 additions and 7 deletions.
2 changes: 1 addition & 1 deletion docs/source/chgres_cube.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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**

Expand Down
10 changes: 9 additions & 1 deletion sorc/chgres_cube.fd/atm_input_data.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
114 changes: 112 additions & 2 deletions sorc/chgres_cube.fd/model_grid.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
6 changes: 3 additions & 3 deletions sorc/chgres_cube.fd/program_setup.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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 " // &
Expand Down

0 comments on commit 1dac855

Please sign in to comment.