Skip to content

Commit

Permalink
Finish updates to model_grid.F90 for GRIB2 grid definition
Browse files Browse the repository at this point in the history
template 1.

Fixes ufs-community#660.
  • Loading branch information
GeorgeGayno-NOAA committed Jun 27, 2022
1 parent 02f14d8 commit ed77aca
Showing 1 changed file with 45 additions and 9 deletions.
54 changes: 45 additions & 9 deletions sorc/chgres_cube.fd/model_grid.F90
Original file line number Diff line number Diff line change
Expand Up @@ -688,6 +688,7 @@ subroutine define_input_grid_grib2(localpet,npets)
endif

kgds = 0

call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, i_input, j_input, res)

ip1_input = i_input + 1
Expand Down Expand Up @@ -720,6 +721,16 @@ subroutine define_input_grid_grib2(localpet,npets)

deallocate(xpts, ypts)

if (localpet == 0) then
print*,'cell center 11 ',rlat(1,1), rlon(1,1)
print*,'cell center im/1 ',rlat(i_input,1), rlon(i_input,1)
print*,'cell center 1/jm ',rlat(1,j_input), rlon(1,j_input)
print*,'cell center im/jm ',rlat(i_input,j_input), rlon(i_input,j_input)
print*,'cell center 2400/1800 ',rlat(2400,1800),rlon(2400,1800)
endif



do j = 1, jp1_input
do i = 1, ip1_input
xpts_corner(i,j) = float(i) - 0.5
Expand All @@ -734,6 +745,14 @@ subroutine define_input_grid_grib2(localpet,npets)
call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2)
endif

if (localpet == 0) then
print*,'cell corner 11 ',rlat_corner(1,1), rlon_corner(1,1)
print*,'cell corner im/1 ',rlat_corner(i_input,1), rlon_corner(i_input,1)
print*,'cell corner 1/jm ',rlat_corner(1,j_input), rlon_corner(1,j_input)
print*,'cell corner im/jm ',rlat_corner(i_input,j_input), rlon_corner(i_input,j_input)
print*,'cell corner 2400 1800 ',rlat_corner(2400,1800),rlon_corner(2400,1800)
endif

deallocate(xpts_corner, ypts_corner)

if (gfld%igdtnum == 0) then ! gfs lat/lon data
Expand Down Expand Up @@ -1452,7 +1471,7 @@ 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

Expand Down Expand Up @@ -1535,11 +1554,7 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
! Lon of cent of rotation
kgds(7) = kgds(7) + 90000.0


print*,'got here ',iscale,ni,nj,kgds(4),kgds(5),kgds(6)
print*,'got here2 ',kgds(7),kgds(8)

DPR = 180.0/3.1415
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)
Expand All @@ -1563,8 +1578,8 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)

RLON = MOD(RLON0+HS*DPR*ACOS(CLON)+3600,360.0)

print*,'got here3 clatr ',clatr, slatr, clonr, slat0, clat0, rlon0
print*,'got here4 rlat ', rlat, rlon
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
Expand Down Expand Up @@ -1592,7 +1607,28 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)

print*,'got here last point ',kgds(12), kgds(13)
print*,'got here last point rotated ', rlat, rlon
stop

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

Expand Down

0 comments on commit ed77aca

Please sign in to comment.