Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

chgres_cube: Implement WMO grib2 template 1 (rotated lat-lon) read capability #902

Merged
merged 4 commits into from
Feb 27, 2024
Merged

chgres_cube: Implement WMO grib2 template 1 (rotated lat-lon) read capability #902

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8df41cc
Implement chgres_cube ability to ready WMO-standard rotated lat-lon g…
LarissaReames-NOAA Feb 22, 2024
8ba94b8
Merge branch 'develop' of https://github.com/NOAA-EMC/UFS_UTILS into …
LarissaReames-NOAA Feb 23, 2024
492eed9
Added documention noting 'RRFS' as a supported option for external_mo…
LarissaReames-NOAA Feb 23, 2024
804ea27
Modified calculation of latin1 for grid template 1 in wind rotation
LarissaReames-NOAA Feb 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/source/chgres_cube.rst
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
GeorgeGayno-NOAA marked this conversation as resolved.
Show resolved Hide resolved
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
Open in desktop
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
Loading