From 1dac85544dd4134de7cb899214c7db67bbe8c0b4 Mon Sep 17 00:00:00 2001 From: Larissa Reames <52886575+LarissaReames-NOAA@users.noreply.github.com> Date: Tue, 27 Feb 2024 10:28:59 -0600 Subject: [PATCH] chgres_cube: Implement WMO grib2 template 1 (rotated lat-lon) read capability (#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 | 2 +- sorc/chgres_cube.fd/atm_input_data.F90 | 10 ++- sorc/chgres_cube.fd/model_grid.F90 | 114 ++++++++++++++++++++++++- sorc/chgres_cube.fd/program_setup.F90 | 6 +- 4 files changed, 125 insertions(+), 7 deletions(-) diff --git a/docs/source/chgres_cube.rst b/docs/source/chgres_cube.rst index 71f1274c3..9af27656e 100644 --- a/docs/source/chgres_cube.rst +++ b/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** diff --git a/sorc/chgres_cube.fd/atm_input_data.F90 b/sorc/chgres_cube.fd/atm_input_data.F90 index 7d995ee0a..fd7d254c1 100644 --- a/sorc/chgres_cube.fd/atm_input_data.F90 +++ b/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 diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 44431f59c..af3d2952f 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/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 diff --git a/sorc/chgres_cube.fd/program_setup.F90 b/sorc/chgres_cube.fd/program_setup.F90 index efec5e3a4..895f71c28 100644 --- a/sorc/chgres_cube.fd/program_setup.F90 +++ b/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 " // &