From fc9befe4c6f26cd1b163ed6abe19e7a598fccadf Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 18 Jan 2022 21:43:04 +0000 Subject: [PATCH 01/14] Begin update of model_grid.F90 to use g2 library. Fixes #591 --- sorc/chgres_cube.fd/CMakeLists.txt | 2 + sorc/chgres_cube.fd/model_grid.F90 | 348 ++++++++++++++++++++++++++--- 2 files changed, 323 insertions(+), 27 deletions(-) diff --git a/sorc/chgres_cube.fd/CMakeLists.txt b/sorc/chgres_cube.fd/CMakeLists.txt index 7b6ad9388..1b086207d 100644 --- a/sorc/chgres_cube.fd/CMakeLists.txt +++ b/sorc/chgres_cube.fd/CMakeLists.txt @@ -42,10 +42,12 @@ target_include_directories(chgres_cube_lib INTERFACE ${mod_dir}) target_link_libraries( chgres_cube_lib PUBLIC + g2::g2_d nemsio::nemsio sfcio::sfcio sigio::sigio bacio::bacio_4 + ip::ip_d sp::sp_d w3nco::w3nco_d esmf diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 028693f35..e358ffa13 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -618,7 +618,7 @@ end subroutine define_input_grid_mosaic subroutine define_input_grid_gfs_grib2(localpet, npets) use wgrib2api - use mpi + use grib_mod use program_setup, only : data_dir_input_grid, & grib2_file_input_grid @@ -629,36 +629,92 @@ subroutine define_input_grid_gfs_grib2(localpet, npets) character(len=250) :: the_file integer :: i, j, rc, clb(2), cub(2) - integer :: ierr + integer :: lugb, ierr + integer :: k, jdisc, jgdtn, jpdtn, lugi + integer :: jids(200), jgdt(200), jpdt(200) + logical :: unpack real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) - real(esmf_kind_r4), allocatable :: lat4(:,:), lon4(:,:) real(esmf_kind_r8), pointer :: lat_src_ptr(:,:) real(esmf_kind_r8), pointer :: lon_src_ptr(:,:) real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:) real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:) - real(esmf_kind_r8) :: deltalon + real(esmf_kind_r8) :: deltalat, deltalon + + real :: lat11, lon11 type(esmf_polekind_flag) :: polekindflag(2) + type(gribfield) :: gfld + print*,"- DEFINE INPUT GRID OBJECT FOR INPUT GRIB2 DATA." + print*,'model grid - gfs latlon' + num_tiles_input_grid = 1 the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid + if (localpet == 0) then print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) rc=grb2_mk_inv(the_file,inv_file) - if (rc /=0) call error_handler("OPENING GRIB2 FILE",rc) + if (rc /=0) call error_handler("OPENING GRIB2 FILE using wgrib2",rc) endif -! Wait for localpet 0 to create inventory - call mpi_barrier(mpi_comm_world, ierr) + lugb=12 + call baopenr(lugb,the_file,rc) + print*,'after g2 open ',rc + + j = 0 ! search at beginning of file + lugi = 0 ! no grib index file + jdisc = 0 ! search for discipline - meterological products + jpdtn = 0 ! search for product definition template number + jgdtn = 0 ! search for grid definition template number; 0 - lat/lon grid + jids = -9999 ! array of values in identification section, set to wildcard + jgdt = -9999 ! array of values in grid definition template 3.m + jpdt = -9999 ! array of values in product definition template 4.n + unpack = .false. ! unpack data + + call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + unpack, k, gfld, rc) + + print*,'after getgb2 temp num ', gfld%igdtnum + print*,'after getgb2 template ', gfld%igdtmpl + + call baclose(lugb, rc) + + i_input = gfld%igdtmpl(8) ! oct 31-34 + j_input = gfld%igdtmpl(9) ! oct 35-38 + + deltalon = float( gfld%igdtmpl(17) )/ 1.e6 + deltalat = float( gfld%igdtmpl(18) )/ 1.e6 + +!print*,'after getgb2 ', deltalon, deltalat - rc = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, & - lat=lat4, lon=lon4) - if (rc /= 1) call error_handler("READING GRIB2 FILE", rc) + allocate(longitude(i_input,j_input)) + allocate(latitude(i_input,j_input)) + +! I think wgrib2 flips the n/s poles by default. +! So for now, follow that convention. +! I don't think the g2 library does that. +! lat11 = float(gfld%igdtmpl(12))/1.e6 +! do j = 1, j_input +! print*,'lat ', j, lat11 - (float(j-1)*dy) +! enddo + +! Follow wgrib2 convention. + lat11 = -(float(gfld%igdtmpl(12))/1.e6) + do j = 1, j_input + print*,'lat ', j, lat11 + (float(j-1)*deltalat) + latitude(:,j) = lat11 + (float(j-1)*deltalat) + enddo + + lon11 = float(gfld%igdtmpl(13))/1.e6 + do i = 1, i_input + print*,'lon ', i, lon11 + (float(i-1)*deltalon) + longitude(i,:) = lon11 + (float(i-1)*deltalon) + enddo ip1_input = i_input + 1 jp1_input = j_input + 1 @@ -693,21 +749,6 @@ subroutine define_input_grid_gfs_grib2(localpet, npets) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldCreate", rc) - allocate(longitude(i_input,j_input)) - allocate(latitude(i_input,j_input)) - - do i = 1, i_input - longitude(i,:) = real(lon4(i,:),kind=esmf_kind_r8) - enddo - - do i = 1, j_input - latitude(:,i) = real(lat4(:,i),kind=esmf_kind_r8) - enddo - - deallocate(lat4, lon4) - - deltalon = abs(longitude(2,1)-longitude(1,1)) - print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE." call ESMF_FieldScatter(longitude_input_grid, longitude, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -807,6 +848,8 @@ end subroutine define_input_grid_gfs_grib2 subroutine define_input_grid_grib2(localpet, npets) use mpi + use grib_mod + use gdswzd_mod use netcdf use wgrib2api use program_setup, only : grib2_file_input_grid, data_dir_input_grid, & @@ -829,12 +872,97 @@ subroutine define_input_grid_grib2(localpet, npets) character(len=10000) :: temp_msg character(len=10) :: temp_num = 'NA' + integer :: k, jdisc, jgdtn, jpdtn, lugb, lugi + integer :: jids(200), jgdt(200), jpdt(200) + integer :: kgds(200), ni, nj, nret + real :: res + real, allocatable :: rlat(:,:), rlon(:,:),xpts(:,:),ypts(:,:) + real, allocatable :: rlatc(:,:), rlonc(:,:) + logical :: unpack + + type(gribfield) :: gfld + num_tiles_input_grid = 1 !inv_file = "chgres.inv" the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid temp_file = trim(fix_dir_input_grid)//"/latlon_grid3.32769.nc" +!!!!!! + + lugb=12 + call baopenr(lugb,the_file,error) + print*,'after g2 open ',error + + j = 0 ! search at beginning of file + lugi = 0 ! no grib index file + jdisc = 0 ! search for discipline - meterological products + jpdtn = 0 ! search for product definition template number + jgdtn = -1 ! search for grid definition template number + jids = -9999 ! array of values in identification section, set to wildcard + jgdt = -9999 ! array of values in grid definition template 3.m + jpdt = -9999 ! array of values in product definition template 4.n + unpack = .true. ! unpack data + + call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + unpack, k, gfld, error) + + print*,'after getgb2 temp num ', gfld%igdtnum + print*,'after getgb2 template ', gfld%igdtmpl + + call baclose(lugb,error) + + kgds = 0 + call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, ni, nj, res) + + print*,'after conversion kgds 1 ',ni, nj + print*,'after conversion kgds 2 ',kgds(1:25) + + allocate(rlat(ni,nj)) + allocate(rlon(ni,nj)) + allocate(xpts(ni,nj)) + allocate(ypts(ni,nj)) + + call gdswzd(kgds,0,(ni*nj),-9999.,xpts,ypts,rlon,rlat,nret) + + if (kgds(1) == 205) then + where (rlon > 180.0) rlon = rlon - 360.0 + endif + + print*,'after gdswzd nret ',nret + print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) + print*,'after gdswzd lat/lon ni/11', rlat(ni,1),rlon(ni,1) + print*,'after gdswzd lat/lon 11/nj', rlat(1,nj),rlon(1,nj) + print*,'after gdswzd lat/lon ni/nj', rlat(ni,nj),rlon(ni,nj) + print*,'after gdswzd lat/lon mid ', rlat(ni/2,nj/2),rlon(ni/2,nj/2) + + deallocate(xpts,ypts) + + allocate(rlatc(ni+1,nj+1)) + allocate(rlonc(ni+1,nj+1)) + allocate(xpts(ni+1,nj+1)) + allocate(ypts(ni+1,nj+1)) + + do j = 1, nj+1 + do i = 1, ni+1 + xpts(i,j) = float(i) - 0.5 + ypts(i,j) = float(j) - 0.5 + enddo + enddo + + call gdswzd(kgds,1,((ni+1)*(nj+1)),-9999.,xpts,ypts,rlonc,rlatc,nret) + if (kgds(1) == 205) then + where (rlonc > 180.0) rlonc = rlonc - 360.0 + endif + print*,'after gdswzd nret ',nret + print*,'after gdswzd corner lat/lon 11', rlatc(1,1),rlonc(1,1) + print*,'after gdswzd corner lat/lon ni/11', rlatc(ni+1,1),rlonc(ni+1,1) + print*,'after gdswzd corner lat/lon 11/nj', rlatc(1,nj+1),rlonc(1,nj+1) + print*,'after gdswzd corner lat/lon ni/nj', rlatc(ni+1,nj+1),rlonc(ni+1,nj+1) + print*,'after gdswzd corner lat/lon mid ', rlatc(ni/2,nj/2),rlonc(ni/2,nj/2) + +!!!!! + call ESMF_FieldGather(latitude_target_grid, lat_target, rootPet=0, tile=1, rc=error) if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & call error_handler("IN FieldGather", error) @@ -860,6 +988,8 @@ subroutine define_input_grid_grib2(localpet, npets) input_grid_type = "rotated_latlon" + print*,'model grid - rotated lat lon' + error=nf90_open(trim(temp_file),nf90_nowrite,ncid) call netcdf_err(error, 'opening: '//trim(temp_file)) @@ -886,10 +1016,26 @@ subroutine define_input_grid_grib2(localpet, npets) error=nf90_get_var(ncid, id_var, longitude_one_tile) call netcdf_err(error, 'reading field' ) + print*,'after netcdf lat/lon 11', latitude_one_tile(1,1),longitude_one_tile(1,1) + print*,'after netcdf lat/lon ni/11', & + latitude_one_tile(ni,1),longitude_one_tile(ni,1) + print*,'after netcdf lat/lon 11/nj', & + latitude_one_tile(1,nj),longitude_one_tile(1,nj) + print*,'after netcdf lat/lon ni/nj', & + latitude_one_tile(ni,nj),longitude_one_tile(ni,nj) + print*,'after netcdf lat/lon mid ', & + latitude_one_tile(ni/2,nj/2),longitude_one_tile(ni/2,nj/2) + elseif (temp_num == "3.0" .or. temp_num == "3.30" .or. temp_num=="30" .or. temp_num == "0") then - if (temp_num =="3.0" .or. temp_num == "0") input_grid_type = "latlon" - if (temp_num =="3.30" .or. temp_num=='30') input_grid_type = "lambert" + if (temp_num =="3.0" .or. temp_num == "0") then + input_grid_type = "latlon" + print*,'model grid - latlon' + endif + if (temp_num =="3.30" .or. temp_num=='30')then + input_grid_type = "lambert" + print*,'model grid - lambert' + endif error = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, & lat=latitude_one_tile, lon=longitude_one_tile) @@ -898,6 +1044,16 @@ subroutine define_input_grid_grib2(localpet, npets) if (localpet==0) print*, "from file lon(1:10,1) = ", longitude_one_tile(1:10,1) if (localpet==0) print*, "from file lat(1,1:10) = ", latitude_one_tile(1,1:10) + print*,'after wgrib2 lat/lon 11', latitude_one_tile(1,1),longitude_one_tile(1,1) + print*,'after wgrib2 lat/lon ni/11', & + latitude_one_tile(ni,1),longitude_one_tile(ni,1) + print*,'after wgrib2 lat/lon 11/nj', & + latitude_one_tile(1,nj),longitude_one_tile(1,nj) + print*,'after wgrib2 lat/lon ni/nj', & + latitude_one_tile(ni,nj),longitude_one_tile(ni,nj) + print*,'after wgrib2 lat/lon mid ', & + latitude_one_tile(ni/2,nj/2),longitude_one_tile(ni/2,nj/2) + elseif (temp_num=="NA") then error = 0 call error_handler("Grid template number cannot be read from the input file. Please " //& @@ -991,6 +1147,13 @@ subroutine define_input_grid_grib2(localpet, npets) do i = clb(1), cub(1) lon_src_ptr(i,j)=real(longitude_one_tile(i,j),esmf_kind_r8) lat_src_ptr(i,j)=real(latitude_one_tile(i,j),esmf_kind_r8) +! lon_src_ptr(i,j)=rlon(i,j) +! lat_src_ptr(i,j)=rlat(i,j) + if (abs(rlon(i,j)-longitude_one_tile(i,j)) > 0.1 .or. & + abs(rlat(i,j)-latitude_one_tile(i,j)) > 0.1) then + print*,'lat/lon diff ', rlon(i,j), longitude_one_tile(i,j), & + rlat(i,j),latitude_one_tile(i,j) + endif enddo enddo @@ -1063,6 +1226,19 @@ subroutine define_input_grid_grib2(localpet, npets) call get_cell_corners(real(latitude_one_tile,esmf_kind_r8), & real(longitude_one_tile, esmf_kind_r8), & lat_src_ptr, lon_src_ptr, dx, clb, cub) + + do j = clb(2),cub(2) + do i = clb(1), cub(1) + if (abs(rlonc(i,j)-lon_src_ptr(i,j)) > 0.1 .or. & + abs(rlatc(i,j)-lat_src_ptr(i,j)) > 0.1) then + print*,'lat/lon corn diff ', rlonc(i,j), lon_src_ptr(i,j), & + rlatc(i,j),lat_src_ptr(i,j) + endif +! lon_src_ptr(i,j)=rlonc(i,j) +! lat_src_ptr(i,j)=rlatc(i,j) + enddo + enddo + endif elseif (trim(input_grid_type) == "rotated_latlon") then !Read the corner coords from file @@ -1083,9 +1259,17 @@ subroutine define_input_grid_grib2(localpet, npets) do i = clb(1), cub(1) lon_src_ptr(i,j)=real(lon_corners(i,j),esmf_kind_r8) lat_src_ptr(i,j)=real(lat_corners(i,j),esmf_kind_r8) +! lon_src_ptr(i,j)=rlonc(i,j) +! lat_src_ptr(i,j)=rlatc(i,j) + if (abs(rlonc(i,j)-lon_corners(i,j)) > 0.1 .or. & + abs(rlatc(i,j)-lat_corners(i,j)) > 0.1) then + print*,'lat/lon corn diff ', rlonc(i,j), lon_corners(i,j), & + rlatc(i,j),lat_corners(i,j) + endif enddo enddo + error= nf90_close(ncid) endif @@ -1564,6 +1748,8 @@ subroutine get_cell_corners( latitude, longitude, latitude_sw, longitude_sw, dx, integer :: i, j + print*,'in routine get_cell_corners' + d = sqrt((dx**2.0_esmf_kind_r8)/2.0_esmf_kind_r8) do j = clb(2),cub(2) @@ -1737,4 +1923,112 @@ subroutine cleanup_input_target_grid_data end subroutine cleanup_input_target_grid_data +!> Convert the grib2 grid description template to +!! to the grib2 grid description sections. +!! +!! @param [in] igdtnum grib2 grid description template number +!! @param [in] igdstmpl length of grib2 grid description template +!! @param [in] igdtlen array of grib2 grid description template octets. +!! @param [out] kgds array of grib1 grid description octets. +!! @param [out] ni i-dimension of grid +!! @param [out] nj j-dimension of grid +!! @param [out] res resolution of grid in km +!! @author George Gayno NCEP/EMC + subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) + + implicit none + + integer, intent(in ) :: igdtnum, igdtlen, igdstmpl(igdtlen) + integer, intent( out) :: kgds(200), ni, nj + integer :: iscale + + real, intent( out) :: res + + kgds=0 + + if (igdtnum.eq.32769) then ! rot lat/lon b grid + + 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(15))/float(iscale)*1000.) ! octs 18-20, + ! Lat of cent of rotation + kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, + ! Lon of cent of rotation + kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, + ! Di + kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, + ! Dj + + 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(12)=nint(float(igdstmpl(20))/float(iscale)*1000.) ! octs 29-31, Lat of + ! last grid point + kgds(13)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 32-34, Lon of + ! last grid point + + kgds(19)=0 ! oct 4, # vert coordinate parameters + kgds(20)=255 ! oct 5, used for thinned grids, set to 255 + + res = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) & + * 0.5 * 111.0 + + elseif(igdtnum==30) then + + kgds(1)=3 ! oct 6, lambert conformal + kgds(2)=igdstmpl(8) ! octs 7-8, Ni + ni = kgds(2) + kgds(3)=igdstmpl(9) ! octs 9-10, Nj + nj = kgds(3) + + iscale = 1e6 + kgds(4) = nint(float(igdstmpl(10))/1000.0) + kgds(5) = nint(float(igdstmpl(11))/1000.0) + + kgds(6)=0 ! oct 17, resolution and component flags + if (igdstmpl(1)==2 ) kgds(6)=64 + if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) ) kgds(6)=kgds(6)+128 + if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8 + + kgds(7) = nint(float(igdstmpl(14))/1000.0) + kgds(8) = nint(float(igdstmpl(15))/1000.0) + kgds(9) = nint(float(igdstmpl(16))/1000.0) + kgds(10) = 0 + + kgds(11) = 0 ! oct 28, scan mode + if (btest(igdstmpl(18),7)) kgds(11) = 128 + if (btest(igdstmpl(18),6)) kgds(11) = kgds(11) + 64 + if (btest(igdstmpl(18),5)) kgds(11) = kgds(11) + 32 + + kgds(12) = nint(float(igdstmpl(19))/1000.0) + kgds(13) = nint(float(igdstmpl(20))/1000.0) + kgds(14) = -90 + kgds(15) = 0 + + else + print*,'grid not defined', igdtnum + + + endif + + end subroutine gdt_to_gds + end module model_grid From 1ebb2f74edcdefcf390658420ee99263cb2b8a9d Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 25 Jan 2022 13:17:47 +0000 Subject: [PATCH 02/14] Start new routine to set the esmf grid object for grib2 data. Fixes #591. --- sorc/chgres_cube.fd/model_grid.F90 | 129 ++++++++++++++++++++++++++++- 1 file changed, 126 insertions(+), 3 deletions(-) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index e358ffa13..98e5ca293 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -127,14 +127,17 @@ subroutine define_input_grid(localpet, npets) trim(input_type) == "gfs_sigio" .or. & trim(input_type) == "gaussian_netcdf") then call define_input_grid_gaussian(localpet, npets) - elseif (trim(external_model) == "GFS" .and. trim(input_type) == "grib2") then - call define_input_grid_gfs_grib2(localpet,npets) +!elseif (trim(external_model) == "GFS" .and. trim(input_type) == "grib2") then +! call define_input_grid_gfs_grib2(localpet,npets) elseif (trim(input_type) == "grib2") then - call define_input_grid_grib2(localpet,npets) + call define_input_grid_grib2_gg(localpet,npets) +! call define_input_grid_grib2(localpet,npets) else call define_input_grid_mosaic(localpet, npets) endif + call abort + end subroutine define_input_grid !> Define grid object for input data on global gaussian grids. @@ -839,6 +842,84 @@ subroutine define_input_grid_gfs_grib2(localpet, npets) end subroutine define_input_grid_gfs_grib2 + subroutine define_input_grid_grib2_gg(localpet,npets) + + use mpi + use grib_mod + use gdswzd_mod + use program_setup, only : grib2_file_input_grid, data_dir_input_grid + + implicit none + + integer, intent(in) :: localpet, npets + + character(len=500) :: the_file + + integer :: j, k, jdisc, jgdtn, jpdtn, lugb, lugi + integer :: jids(200), jgdt(200), jpdt(200), error + integer :: ni, nj, kgds(200), nret + + logical :: unpack + + real :: res + real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:) + + type(gribfield) :: gfld + + the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid + + lugb=12 + + if (localpet == 0) then + call baopenr(lugb,the_file,error) + print*,'after g2 open ',error + + j = 0 ! search at beginning of file + lugi = 0 ! no grib index file + jdisc = 0 ! search for discipline - meterological products + jpdtn = 0 ! search for product definition template number + jgdtn = -1 ! search for grid definition template number + jids = -9999 ! array of values in identification section, set to wildcard + jgdt = -9999 ! array of values in grid definition template 3.m + jpdt = -9999 ! array of values in product definition template 4.n + unpack = .true. ! unpack data + + call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + unpack, k, gfld, error) + + print*,'after getgb2 temp num ', gfld%igdtnum + print*,'after getgb2 template ', gfld%igdtmpl + + call baclose(lugb,error) + + kgds = 0 + call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, ni, nj, res) + + print*,'after conversion kgds 1 ',ni, nj + print*,'after conversion kgds 2 ',kgds(1:25) + + allocate(rlat(ni,nj)) + allocate(rlon(ni,nj)) + allocate(xpts(ni,nj)) + allocate(ypts(ni,nj)) + + call gdswzd(kgds,0,(ni*nj),-9999.,xpts,ypts,rlon,rlat,nret) + + print*,'after gdswzd nret ',nret + + print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) + print*,'after gdswzd lat/lon ni/11', rlat(ni,1),rlon(ni,1) + print*,'after gdswzd lat/lon 11/nj', rlat(1,nj),rlon(1,nj) + print*,'after gdswzd lat/lon ni/nj', rlat(ni,nj),rlon(ni,nj) + print*,'after gdswzd lat/lon mid ', rlat(ni/2,nj/2),rlon(ni/2,nj/2) + + endif + + call MPI_BARRIER(MPI_COMM_WORLD, error) + stop + + end subroutine define_input_grid_grib2_gg + !> Define input grid object for non-GFS grib2 data. !! !! @param [in] localpet ESMF local persistent execution thread @@ -930,12 +1011,16 @@ subroutine define_input_grid_grib2(localpet, npets) endif print*,'after gdswzd nret ',nret + + if (localpet) then print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) print*,'after gdswzd lat/lon ni/11', rlat(ni,1),rlon(ni,1) print*,'after gdswzd lat/lon 11/nj', rlat(1,nj),rlon(1,nj) print*,'after gdswzd lat/lon ni/nj', rlat(ni,nj),rlon(ni,nj) print*,'after gdswzd lat/lon mid ', rlat(ni/2,nj/2),rlon(ni/2,nj/2) + endif + deallocate(xpts,ypts) allocate(rlatc(ni+1,nj+1)) @@ -954,13 +1039,18 @@ subroutine define_input_grid_grib2(localpet, npets) if (kgds(1) == 205) then where (rlonc > 180.0) rlonc = rlonc - 360.0 endif + + if (localpet) then print*,'after gdswzd nret ',nret print*,'after gdswzd corner lat/lon 11', rlatc(1,1),rlonc(1,1) print*,'after gdswzd corner lat/lon ni/11', rlatc(ni+1,1),rlonc(ni+1,1) print*,'after gdswzd corner lat/lon 11/nj', rlatc(1,nj+1),rlonc(1,nj+1) print*,'after gdswzd corner lat/lon ni/nj', rlatc(ni+1,nj+1),rlonc(ni+1,nj+1) print*,'after gdswzd corner lat/lon mid ', rlatc(ni/2,nj/2),rlonc(ni/2,nj/2) + endif + + call abort !!!!! call ESMF_FieldGather(latitude_target_grid, lat_target, rootPet=0, tile=1, rc=error) @@ -2023,9 +2113,42 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) kgds(14) = -90 kgds(15) = 0 + elseif(igdtnum==0) then ! lat/lon grid + + iscale=igdstmpl(10)*igdstmpl(11) + if (iscale == 0) iscale = 1e6 + kgds(1)=0 ! oct 6 + 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(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point + kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point + kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, di + kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, dj + + 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(12)=0 ! octs 29-32, reserved + kgds(19)=0 ! oct 4, # vert coordinate parameters + kgds(20)=255 ! oct 5, used for thinned grids, set to 255 + else print*,'grid not defined', igdtnum + call abort + endif From d348db9a157c6f48dad9108a8abb63feb7bff7f9 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 25 Jan 2022 16:00:56 -0600 Subject: [PATCH 03/14] Begin updates to add non-latlon grids. Fixes #591 --- sorc/chgres_cube.fd/model_grid.F90 | 253 +++++++++++++++++++++++++---- 1 file changed, 224 insertions(+), 29 deletions(-) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 98e5ca293..f227ebd6b 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -136,8 +136,6 @@ subroutine define_input_grid(localpet, npets) call define_input_grid_mosaic(localpet, npets) endif - call abort - end subroutine define_input_grid !> Define grid object for input data on global gaussian grids. @@ -844,6 +842,7 @@ end subroutine define_input_grid_gfs_grib2 subroutine define_input_grid_grib2_gg(localpet,npets) + use wgrib2api use mpi use grib_mod use gdswzd_mod @@ -855,14 +854,25 @@ subroutine define_input_grid_grib2_gg(localpet,npets) character(len=500) :: the_file - integer :: j, k, jdisc, jgdtn, jpdtn, lugb, lugi - integer :: jids(200), jgdt(200), jpdt(200), error - integer :: ni, nj, kgds(200), nret + integer :: i, j, k, jdisc, jgdtn, jpdtn, lugb, lugi + integer :: jids(200), jgdt(200), jpdt(200), rc + integer :: kgds(200), nret, clb(2), cub(2) logical :: unpack real :: res real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:) + real, allocatable :: rlon_corner(:,:),rlat_corner(:,:) + real(esmf_kind_r8), allocatable :: latitude(:,:) + real(esmf_kind_r8), allocatable :: longitude(:,:) + real(esmf_kind_r8), allocatable :: latitude_corner(:,:) + real(esmf_kind_r8), allocatable :: longitude_corner(:,:) + real(esmf_kind_r8), pointer :: lat_src_ptr(:,:) + real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:) + real(esmf_kind_r8), pointer :: lon_src_ptr(:,:) + real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:) + + type(esmf_polekind_flag) :: polekindflag(2) type(gribfield) :: gfld @@ -870,9 +880,8 @@ subroutine define_input_grid_grib2_gg(localpet,npets) lugb=12 - if (localpet == 0) then - call baopenr(lugb,the_file,error) - print*,'after g2 open ',error + call baopenr(lugb,the_file,rc) + print*,'after g2 open ',rc j = 0 ! search at beginning of file lugi = 0 ! no grib index file @@ -885,38 +894,222 @@ subroutine define_input_grid_grib2_gg(localpet,npets) unpack = .true. ! unpack data call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & - unpack, k, gfld, error) + unpack, k, gfld, rc) - print*,'after getgb2 temp num ', gfld%igdtnum - print*,'after getgb2 template ', gfld%igdtmpl + if (localpet == 0)then + print*,'after getgb2 temp num ', gfld%igdtnum + print*,'after getgb2 template ', gfld%igdtmpl + endif - call baclose(lugb,error) + call baclose(lugb,rc) kgds = 0 - call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, ni, nj, res) + call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, i_input, j_input, res) - print*,'after conversion kgds 1 ',ni, nj - print*,'after conversion kgds 2 ',kgds(1:25) + if (localpet == 0) then + print*,'after conversion kgds ni/nj ',i_input, j_input + print*,'after conversion kgds ',kgds(1:25) + endif - allocate(rlat(ni,nj)) - allocate(rlon(ni,nj)) - allocate(xpts(ni,nj)) - allocate(ypts(ni,nj)) + allocate(rlat(i_input,j_input)) + allocate(rlon(i_input,j_input)) + allocate(rlat_corner(i_input,j_input)) + allocate(rlon_corner(i_input,j_input)) + allocate(xpts(i_input,j_input)) + allocate(ypts(i_input,j_input)) - call gdswzd(kgds,0,(ni*nj),-9999.,xpts,ypts,rlon,rlat,nret) + call gdswzd(kgds,0,(i_input*j_input),-9999.,xpts,ypts,rlon,rlat,nret) print*,'after gdswzd nret ',nret - print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) - print*,'after gdswzd lat/lon ni/11', rlat(ni,1),rlon(ni,1) - print*,'after gdswzd lat/lon 11/nj', rlat(1,nj),rlon(1,nj) - print*,'after gdswzd lat/lon ni/nj', rlat(ni,nj),rlon(ni,nj) - print*,'after gdswzd lat/lon mid ', rlat(ni/2,nj/2),rlon(ni/2,nj/2) + if (localpet == 0) then + print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) + print*,'after gdswzd lat/lon ni/11', rlat(i_input,1),rlon(i_input,1) + print*,'after gdswzd lat/lon 11/nj', rlat(1,j_input),rlon(1,j_input) + print*,'after gdswzd lat/lon ni/nj', rlat(i_input,j_input),rlon(i_input,j_input) + print*,'after gdswzd lat/lon mid ', rlat(i_input/2,j_input/2),rlon(i_input/2,j_input/2) + endif + + xpts = xpts - 0.5 + ypts = ypts - 0.5 + + call gdswzd(kgds,1,(i_input*j_input),-9999.,xpts,ypts,rlon_corner,rlat_corner,nret) + + print*,'after gdswzd corner nret ',nret + + if (localpet == 0) then + print*,'after gdswzd lat/lon corner 11', rlat_corner(1,1),rlon_corner(1,1) + print*,'after gdswzd lat/lon corner ni/11', rlat_corner(i_input,1),rlon_corner(i_input,1) + print*,'after gdswzd lat/lon corner 11/nj', rlat_corner(1,j_input),rlon_corner(1,j_input) + print*,'after gdswzd lat/lon corner ni/nj', rlat_corner(i_input,j_input),rlon_corner(i_input,j_input) + print*,'after gdswzd lat/lon corner mid ', rlat_corner(i_input/2,j_input/2),rlon_corner(i_input/2,j_input/2) + endif + + + if (gfld%igdtnum == 0) then ! gfs lat/lon data + + print*,"- CALL GridCreate1PeriDim FOR INPUT GRID." + + polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE + + input_grid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & + maxIndex=(/i_input,j_input/), & + polekindflag=polekindflag, & + periodicDim=1, & + poleDim=2, & + coordSys=ESMF_COORDSYS_SPH_DEG, & + regDecomp=(/1,npets/), & + indexflag=ESMF_INDEX_GLOBAL, rc=rc) + + else + + print*,"- CALL GridCreateNoPeriDim FOR INPUT GRID." + + input_grid = ESMF_GridCreateNoPeriDim(maxIndex=(/i_input,j_input/), & + indexflag=ESMF_INDEX_GLOBAL, & + rc=rc) + + endif + + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridCreate1PeriDim", rc) + + print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE." + latitude_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="input_grid_latitude", rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE." + longitude_input_grid = ESMF_FieldCreate(input_grid, & + typekind=ESMF_TYPEKIND_R8, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + name="input_grid_longitude", rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldCreate", rc) + + allocate(latitude(i_input,j_input)) + allocate(longitude(i_input,j_input)) + + latitude = rlat + longitude = rlon + + print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE." + call ESMF_FieldScatter(longitude_input_grid, longitude, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE." + call ESMF_FieldScatter(latitude_input_grid, latitude, rootpet=0, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN FieldScatter", rc) + + print*,"- CALL GridAddCoord FOR INPUT GRID." + call ESMF_GridAddCoord(input_grid, & + staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord", rc) + + print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." + nullify(lon_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=1, & + farrayPtr=lon_src_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", rc) + + print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." + nullify(lat_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CENTER, & + coordDim=2, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=lat_src_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", rc) + + do j = clb(2), cub(2) + do i = clb(1), cub(1) + lon_src_ptr(i,j) = longitude(i,j) + if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8 + lat_src_ptr(i,j) = latitude(i,j) + enddo + enddo + + ip1_input = i_input + 1 + jp1_input = j_input + 1 + + if (gfld%igdtnum == 0) then ! gfs lat/lon data + + allocate(latitude_corner(i_input,jp1_input)) + allocate(longitude_corner(i_input,jp1_input)) + + latitude_corner(:,1:j_input) = rlat_corner + latitude_corner(:,jp1_input) = -(latitude_corner(:,1)) + longitude_corner(:,1:j_input) = rlon_corner + longitude_corner(:,jp1_input) = longitude_corner(:,j_input) + if (localpet == 0) then + do j = 1, jp1_input + print*,'lat corner ',j, latitude_corner(1,j) + enddo + endif + + else + + allocate(latitude_corner(i_input,j_input)) + allocate(longitude_corner(i_input,j_input)) + + latitude_corner = rlat_corner + longitude_corner = rlon_corner + + endif + + print*,"- CALL GridAddCoord FOR INPUT GRID." + call ESMF_GridAddCoord(input_grid, & + staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridAddCoord", rc) + + print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." + nullify(lon_corner_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CORNER, & + coordDim=1, & + farrayPtr=lon_corner_src_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", rc) + + print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." + nullify(lat_corner_src_ptr) + call ESMF_GridGetCoord(input_grid, & + staggerLoc=ESMF_STAGGERLOC_CORNER, & + coordDim=2, & + computationalLBound=clb, & + computationalUBound=cub, & + farrayPtr=lat_corner_src_ptr, rc=rc) + if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & + call error_handler("IN GridGetCoord", rc) + + do j = clb(2), cub(2) + do i = clb(1), cub(1) + lon_corner_src_ptr(i,j) = longitude_corner(i,j) + if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8 + lat_corner_src_ptr(i,j) = latitude_corner(i,j) + enddo + enddo + + if (localpet == 0) then + print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) + rc=grb2_mk_inv(the_file,inv_file) + if (rc /=0) call error_handler("OPENING GRIB2 FILE using wgrib2",rc) endif - call MPI_BARRIER(MPI_COMM_WORLD, error) - stop + call MPI_BARRIER(MPI_COMM_WORLD, rc) end subroutine define_input_grid_grib2_gg @@ -2122,7 +2315,8 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) 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(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point + kgds(7)=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 @@ -2130,7 +2324,8 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) 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(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point +! kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point + kgds(4)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, di kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, dj From 2425c5d99493a941ce823fec4eab36511526f54a Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 27 Jan 2022 15:03:02 -0600 Subject: [PATCH 04/14] Call gdswzd twice. Once for the grid point centers and once for the corner points. Fixes #591. --- sorc/chgres_cube.fd/model_grid.F90 | 79 ++++++++++++++++-------------- 1 file changed, 42 insertions(+), 37 deletions(-) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index f227ebd6b..3bc8a00d8 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -863,6 +863,7 @@ subroutine define_input_grid_grib2_gg(localpet,npets) real :: res real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:) real, allocatable :: rlon_corner(:,:),rlat_corner(:,:) + real, allocatable :: xpts_corner(:,:),ypts_corner(:,:) real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) real(esmf_kind_r8), allocatable :: latitude_corner(:,:) @@ -903,6 +904,17 @@ subroutine define_input_grid_grib2_gg(localpet,npets) call baclose(lugb,rc) + if (gfld%igdtnum == 0) then + input_grid_type = 'latlon' + elseif (gfld%igdtnum == 30) then + input_grid_type = 'lambert' + elseif (gfld%igdtnum == 32769) then + input_grid_type = 'rotated_latlon' + else + print*,'grid not supported' + call abort + endif + kgds = 0 call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, i_input, j_input, res) @@ -911,14 +923,26 @@ subroutine define_input_grid_grib2_gg(localpet,npets) print*,'after conversion kgds ',kgds(1:25) endif + ip1_input = i_input + 1 + jp1_input = j_input + 1 + allocate(rlat(i_input,j_input)) allocate(rlon(i_input,j_input)) - allocate(rlat_corner(i_input,j_input)) - allocate(rlon_corner(i_input,j_input)) allocate(xpts(i_input,j_input)) allocate(ypts(i_input,j_input)) + allocate(rlat_corner(ip1_input,jp1_input)) + allocate(rlon_corner(ip1_input,jp1_input)) + allocate(xpts_corner(ip1_input,jp1_input)) + allocate(ypts_corner(ip1_input,jp1_input)) + + do j = 1, j_input + do i = 1, i_input + xpts(i,j) = float(i) + ypts(i,j) = float(j) + enddo + enddo - call gdswzd(kgds,0,(i_input*j_input),-9999.,xpts,ypts,rlon,rlat,nret) + call gdswzd(kgds,1,(i_input*j_input),-9999.,xpts,ypts,rlon,rlat,nret) print*,'after gdswzd nret ',nret @@ -930,22 +954,26 @@ subroutine define_input_grid_grib2_gg(localpet,npets) print*,'after gdswzd lat/lon mid ', rlat(i_input/2,j_input/2),rlon(i_input/2,j_input/2) endif - xpts = xpts - 0.5 - ypts = ypts - 0.5 + do j = 1, jp1_input + do i = 1, ip1_input + xpts_corner(i,j) = float(i) - 0.5 + ypts_corner(i,j) = float(j) - 0.5 + enddo + enddo - call gdswzd(kgds,1,(i_input*j_input),-9999.,xpts,ypts,rlon_corner,rlat_corner,nret) + call gdswzd(kgds,1,(ip1_input*jp1_input),-9999.,xpts_corner,ypts_corner,rlon_corner,rlat_corner,nret) print*,'after gdswzd corner nret ',nret if (localpet == 0) then print*,'after gdswzd lat/lon corner 11', rlat_corner(1,1),rlon_corner(1,1) - print*,'after gdswzd lat/lon corner ni/11', rlat_corner(i_input,1),rlon_corner(i_input,1) - print*,'after gdswzd lat/lon corner 11/nj', rlat_corner(1,j_input),rlon_corner(1,j_input) - print*,'after gdswzd lat/lon corner ni/nj', rlat_corner(i_input,j_input),rlon_corner(i_input,j_input) + print*,'after gdswzd lat/lon corner ni/11', rlat_corner(ip1_input,1),rlon_corner(ip1_input,1) + print*,'after gdswzd lat/lon corner 11/nj', rlat_corner(1,jp1_input),rlon_corner(1,jp1_input) + print*,'after gdswzd lat/lon corner ni/nj', rlat_corner(ip1_input,jp1_input),rlon_corner(ip1_input,jp1_input) print*,'after gdswzd lat/lon corner mid ', rlat_corner(i_input/2,j_input/2),rlon_corner(i_input/2,j_input/2) + print*,'after gdswzd max/min corner ',maxval(rlat_corner),minval(rlat_corner),maxval(rlon_corner),minval(rlon_corner) endif - if (gfld%igdtnum == 0) then ! gfs lat/lon data print*,"- CALL GridCreate1PeriDim FOR INPUT GRID." @@ -1040,34 +1068,11 @@ subroutine define_input_grid_grib2_gg(localpet,npets) enddo enddo - ip1_input = i_input + 1 - jp1_input = j_input + 1 - - if (gfld%igdtnum == 0) then ! gfs lat/lon data - - allocate(latitude_corner(i_input,jp1_input)) - allocate(longitude_corner(i_input,jp1_input)) + allocate(latitude_corner(ip1_input,jp1_input)) + allocate(longitude_corner(ip1_input,jp1_input)) - latitude_corner(:,1:j_input) = rlat_corner - latitude_corner(:,jp1_input) = -(latitude_corner(:,1)) - longitude_corner(:,1:j_input) = rlon_corner - longitude_corner(:,jp1_input) = longitude_corner(:,j_input) - - if (localpet == 0) then - do j = 1, jp1_input - print*,'lat corner ',j, latitude_corner(1,j) - enddo - endif - - else - - allocate(latitude_corner(i_input,j_input)) - allocate(longitude_corner(i_input,j_input)) - - latitude_corner = rlat_corner - longitude_corner = rlon_corner - - endif + latitude_corner = rlat_corner + longitude_corner = rlon_corner print*,"- CALL GridAddCoord FOR INPUT GRID." call ESMF_GridAddCoord(input_grid, & From d00235e360131b6e821147016bbb77127694fb8f Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 1 Feb 2022 15:26:23 -0600 Subject: [PATCH 05/14] Output lat/lons to a netcdf file. --- sorc/chgres_cube.fd/model_grid.F90 | 66 +++++++++++++++++++++++++++++- 1 file changed, 65 insertions(+), 1 deletion(-) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 3bc8a00d8..cdaf6b673 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -842,6 +842,7 @@ end subroutine define_input_grid_gfs_grib2 subroutine define_input_grid_grib2_gg(localpet,npets) + use netcdf use wgrib2api use mpi use grib_mod @@ -857,12 +858,19 @@ subroutine define_input_grid_grib2_gg(localpet,npets) integer :: i, j, k, jdisc, jgdtn, jpdtn, lugb, lugi integer :: jids(200), jgdt(200), jpdt(200), rc integer :: kgds(200), nret, clb(2), cub(2) + integer :: error, ncid, dim_i, dim_ip1, dim_j, dim_jp1 + integer :: fsize=65536 + integer :: header_buffer_val = 16384 + integer :: id_gridlat, id_gridlon, id_gridrot + integer :: id_gridlat_corners, id_gridlon_corners + integer :: i_super, j_super, dim_i_super, dim_j_super logical :: unpack real :: res real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:) real, allocatable :: rlon_corner(:,:),rlat_corner(:,:) + real, allocatable :: rlon_super(:,:),rlat_super(:,:) real, allocatable :: xpts_corner(:,:),ypts_corner(:,:) real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) @@ -965,6 +973,12 @@ subroutine define_input_grid_grib2_gg(localpet,npets) print*,'after gdswzd corner nret ',nret + i_super = 2 * i_input + 1 + j_super = 2 * j_input + 1 + + allocate(rlon_super(i_super,j_super)) + allocate(rlat_super(i_super,j_super)) + if (localpet == 0) then print*,'after gdswzd lat/lon corner 11', rlat_corner(1,1),rlon_corner(1,1) print*,'after gdswzd lat/lon corner ni/11', rlat_corner(ip1_input,1),rlon_corner(ip1_input,1) @@ -972,6 +986,55 @@ subroutine define_input_grid_grib2_gg(localpet,npets) print*,'after gdswzd lat/lon corner ni/nj', rlat_corner(ip1_input,jp1_input),rlon_corner(ip1_input,jp1_input) print*,'after gdswzd lat/lon corner mid ', rlat_corner(i_input/2,j_input/2),rlon_corner(i_input/2,j_input/2) print*,'after gdswzd max/min corner ',maxval(rlat_corner),minval(rlat_corner),maxval(rlon_corner),minval(rlon_corner) + ncid = 34 + error = nf90_create("./latlon.nc", IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), & + ncid, initialsize=0, chunksize=fsize) + call netcdf_err(error, 'CREATING FILE=latlon.nc' ) + + error = nf90_def_dim(ncid, 'ny', j_input, dim_j) + call netcdf_err(error, 'DEFINING j DIMENSION' ) + + error = nf90_def_dim(ncid, 'nx', i_input, dim_i) + call netcdf_err(error, 'DEFINING i DIMENSION' ) + + error = nf90_def_dim(ncid, 'ny_stag', jp1_input, dim_jp1) + call netcdf_err(error, 'DEFINING jp1 DIMENSION' ) + + error = nf90_def_dim(ncid, 'nx_stag', ip1_input, dim_ip1) + call netcdf_err(error, 'DEFINING ip1 DIMENSION' ) + + error = nf90_def_var(ncid, 'gridlat', NF90_FLOAT, (/dim_i,dim_j/), id_gridlat) + call netcdf_err(error, 'DEFINING gridlat' ) + + error = nf90_def_var(ncid, 'gridlon', NF90_FLOAT, (/dim_i,dim_j/), id_gridlon) + call netcdf_err(error, 'DEFINING gridlon' ) + + error = nf90_def_var(ncid, 'gridrot', NF90_FLOAT, (/dim_i,dim_j/), id_gridrot) + call netcdf_err(error, 'DEFINING gridrot' ) + + error = nf90_def_var(ncid, 'gridlat_corners', NF90_FLOAT, (/dim_ip1,dim_jp1/), id_gridlat_corners) + call netcdf_err(error, 'DEFINING gridlat_corners' ) + + error = nf90_def_var(ncid, 'gridlon_corners', NF90_FLOAT, (/dim_ip1,dim_jp1/), id_gridlon_corners) + call netcdf_err(error, 'DEFINING gridlon_corners' ) + + error = nf90_enddef(ncid, header_buffer_val,4,0,4) + call netcdf_err(error, 'DEFINING HEADER' ) + + error = nf90_put_var(ncid, id_gridlat, real(rlat,4) ) + call netcdf_err(error, 'writing gridlat' ) + + error = nf90_put_var(ncid, id_gridlon, real(rlon,4) ) + call netcdf_err(error, 'writing gridlon' ) + + error = nf90_put_var(ncid, id_gridlat_corners, real(rlat_corner,4) ) + call netcdf_err(error, 'writing gridlat_corner' ) + + error = nf90_put_var(ncid, id_gridlon_corners, real(rlon_corner,4) ) + call netcdf_err(error, 'writing gridlon_corner' ) + + error = nf90_close(ncid) + endif if (gfld%igdtnum == 0) then ! gfs lat/lon data @@ -1210,7 +1273,7 @@ subroutine define_input_grid_grib2(localpet, npets) print*,'after gdswzd nret ',nret - if (localpet) then + if (localpet == 0) then print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) print*,'after gdswzd lat/lon ni/11', rlat(ni,1),rlon(ni,1) print*,'after gdswzd lat/lon 11/nj', rlat(1,nj),rlon(1,nj) @@ -1248,6 +1311,7 @@ subroutine define_input_grid_grib2(localpet, npets) endif + call abort !!!!! From 474401f680e976938a8684bdc9b559009d9e328d Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 4 Feb 2022 14:59:36 +0000 Subject: [PATCH 06/14] Adjust longitude convention of rotated lat/lon grids in output netcdf file. --- sorc/chgres_cube.fd/model_grid.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index cdaf6b673..5997f5192 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -986,6 +986,12 @@ subroutine define_input_grid_grib2_gg(localpet,npets) print*,'after gdswzd lat/lon corner ni/nj', rlat_corner(ip1_input,jp1_input),rlon_corner(ip1_input,jp1_input) print*,'after gdswzd lat/lon corner mid ', rlat_corner(i_input/2,j_input/2),rlon_corner(i_input/2,j_input/2) print*,'after gdswzd max/min corner ',maxval(rlat_corner),minval(rlat_corner),maxval(rlon_corner),minval(rlon_corner) + + if (gfld%igdtnum == 32769) then + where(rlon > 180.0) rlon = rlon - 360.0 + where(rlon_corner > 180.0) rlon_corner = rlon_corner - 360.0 + endif + ncid = 34 error = nf90_create("./latlon.nc", IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), & ncid, initialsize=0, chunksize=fsize) From c66b5c67c3a2b5037c57ce29d599ec81675e7b30 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 4 Feb 2022 19:29:05 +0000 Subject: [PATCH 07/14] Remove unused routines from model_grid.F90 Part of #591. --- sorc/chgres_cube.fd/model_grid.F90 | 789 +---------------------------- 1 file changed, 8 insertions(+), 781 deletions(-) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 5997f5192..015cf9bde 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -116,7 +116,7 @@ module model_grid !! @author George Gayno NCEP/EMC subroutine define_input_grid(localpet, npets) - use program_setup, only : input_type, external_model + use program_setup, only : input_type implicit none @@ -127,11 +127,8 @@ subroutine define_input_grid(localpet, npets) trim(input_type) == "gfs_sigio" .or. & trim(input_type) == "gaussian_netcdf") then call define_input_grid_gaussian(localpet, npets) -!elseif (trim(external_model) == "GFS" .and. trim(input_type) == "grib2") then -! call define_input_grid_gfs_grib2(localpet,npets) elseif (trim(input_type) == "grib2") then - call define_input_grid_grib2_gg(localpet,npets) -! call define_input_grid_grib2(localpet,npets) + call define_input_grid_grib2(localpet,npets) else call define_input_grid_mosaic(localpet, npets) endif @@ -610,237 +607,14 @@ subroutine define_input_grid_mosaic(localpet, npets) end subroutine define_input_grid_mosaic -!> Define input grid object for GFS grib2 data. Only works for data on -!! global lat/lon or gaussian grids. +!> Define input grid object for non-GFS grib2 data. !! !! @param [in] localpet ESMF local persistent execution thread !! @param [in] npets Number of persistent execution threads -!! @author George Gayno NCEP/EMC - subroutine define_input_grid_gfs_grib2(localpet, npets) - - use wgrib2api - use grib_mod - use program_setup, only : data_dir_input_grid, & - grib2_file_input_grid - - implicit none - - integer, intent(in) :: localpet, npets - - character(len=250) :: the_file - - integer :: i, j, rc, clb(2), cub(2) - integer :: lugb, ierr - integer :: k, jdisc, jgdtn, jpdtn, lugi - integer :: jids(200), jgdt(200), jpdt(200) - logical :: unpack - - real(esmf_kind_r8), allocatable :: latitude(:,:) - real(esmf_kind_r8), allocatable :: longitude(:,:) - real(esmf_kind_r8), pointer :: lat_src_ptr(:,:) - real(esmf_kind_r8), pointer :: lon_src_ptr(:,:) - real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:) - real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:) - real(esmf_kind_r8) :: deltalat, deltalon - - real :: lat11, lon11 - - type(esmf_polekind_flag) :: polekindflag(2) - - type(gribfield) :: gfld - - print*,"- DEFINE INPUT GRID OBJECT FOR INPUT GRIB2 DATA." - - print*,'model grid - gfs latlon' - - num_tiles_input_grid = 1 - - the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid - - if (localpet == 0) then - print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) - rc=grb2_mk_inv(the_file,inv_file) - if (rc /=0) call error_handler("OPENING GRIB2 FILE using wgrib2",rc) - endif - - lugb=12 - call baopenr(lugb,the_file,rc) - print*,'after g2 open ',rc - - j = 0 ! search at beginning of file - lugi = 0 ! no grib index file - jdisc = 0 ! search for discipline - meterological products - jpdtn = 0 ! search for product definition template number - jgdtn = 0 ! search for grid definition template number; 0 - lat/lon grid - jids = -9999 ! array of values in identification section, set to wildcard - jgdt = -9999 ! array of values in grid definition template 3.m - jpdt = -9999 ! array of values in product definition template 4.n - unpack = .false. ! unpack data - - call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & - unpack, k, gfld, rc) - - print*,'after getgb2 temp num ', gfld%igdtnum - print*,'after getgb2 template ', gfld%igdtmpl - - call baclose(lugb, rc) - - i_input = gfld%igdtmpl(8) ! oct 31-34 - j_input = gfld%igdtmpl(9) ! oct 35-38 - - deltalon = float( gfld%igdtmpl(17) )/ 1.e6 - deltalat = float( gfld%igdtmpl(18) )/ 1.e6 - -!print*,'after getgb2 ', deltalon, deltalat - - allocate(longitude(i_input,j_input)) - allocate(latitude(i_input,j_input)) - -! I think wgrib2 flips the n/s poles by default. -! So for now, follow that convention. -! I don't think the g2 library does that. -! lat11 = float(gfld%igdtmpl(12))/1.e6 -! do j = 1, j_input -! print*,'lat ', j, lat11 - (float(j-1)*dy) -! enddo - -! Follow wgrib2 convention. - lat11 = -(float(gfld%igdtmpl(12))/1.e6) - do j = 1, j_input - print*,'lat ', j, lat11 + (float(j-1)*deltalat) - latitude(:,j) = lat11 + (float(j-1)*deltalat) - enddo - - lon11 = float(gfld%igdtmpl(13))/1.e6 - do i = 1, i_input - print*,'lon ', i, lon11 + (float(i-1)*deltalon) - longitude(i,:) = lon11 + (float(i-1)*deltalon) - enddo - - ip1_input = i_input + 1 - jp1_input = j_input + 1 - - polekindflag(1:2) = ESMF_POLEKIND_MONOPOLE - - print*,"- CALL GridCreate1PeriDim FOR INPUT GRID." - input_grid = ESMF_GridCreate1PeriDim(minIndex=(/1,1/), & - maxIndex=(/i_input,j_input/), & - polekindflag=polekindflag, & - periodicDim=1, & - poleDim=2, & - coordSys=ESMF_COORDSYS_SPH_DEG, & - regDecomp=(/1,npets/), & - indexflag=ESMF_INDEX_GLOBAL, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridCreate1PeriDim", rc) - - print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE." - latitude_input_grid = ESMF_FieldCreate(input_grid, & - typekind=ESMF_TYPEKIND_R8, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - name="input_grid_latitude", rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldCreate", rc) - - print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE." - longitude_input_grid = ESMF_FieldCreate(input_grid, & - typekind=ESMF_TYPEKIND_R8, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - name="input_grid_longitude", rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldCreate", rc) - - print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE." - call ESMF_FieldScatter(longitude_input_grid, longitude, rootpet=0, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldScatter", rc) - - print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE." - call ESMF_FieldScatter(latitude_input_grid, latitude, rootpet=0, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldScatter", rc) - - print*,"- CALL GridAddCoord FOR INPUT GRID." - call ESMF_GridAddCoord(input_grid, & - staggerloc=ESMF_STAGGERLOC_CENTER, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridAddCoord", rc) - - print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." - nullify(lon_src_ptr) - call ESMF_GridGetCoord(input_grid, & - staggerLoc=ESMF_STAGGERLOC_CENTER, & - coordDim=1, & - farrayPtr=lon_src_ptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord", rc) - - print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." - nullify(lat_src_ptr) - call ESMF_GridGetCoord(input_grid, & - staggerLoc=ESMF_STAGGERLOC_CENTER, & - coordDim=2, & - computationalLBound=clb, & - computationalUBound=cub, & - farrayPtr=lat_src_ptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord", rc) - - do j = clb(2), cub(2) - do i = clb(1), cub(1) - lon_src_ptr(i,j) = longitude(i,j) - if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8 - lat_src_ptr(i,j) = latitude(i,j) - enddo - enddo - - print*,"- CALL GridAddCoord FOR INPUT GRID." - call ESMF_GridAddCoord(input_grid, & - staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridAddCoord", rc) - - print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." - nullify(lon_corner_src_ptr) - call ESMF_GridGetCoord(input_grid, & - staggerLoc=ESMF_STAGGERLOC_CORNER, & - coordDim=1, & - farrayPtr=lon_corner_src_ptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord", rc) - - print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." - nullify(lat_corner_src_ptr) - call ESMF_GridGetCoord(input_grid, & - staggerLoc=ESMF_STAGGERLOC_CORNER, & - coordDim=2, & - computationalLBound=clb, & - computationalUBound=cub, & - farrayPtr=lat_corner_src_ptr, rc=rc) - if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord", rc) - - do j = clb(2), cub(2) - do i = clb(1), cub(1) - lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon) - if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8 - if (j == 1) then - lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8 - cycle - endif - if (j == jp1_input) then - lat_corner_src_ptr(i,j) = +90.0_esmf_kind_r8 - cycle - endif - lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j)) - enddo - enddo - - deallocate(latitude,longitude) - - end subroutine define_input_grid_gfs_grib2 - - subroutine define_input_grid_grib2_gg(localpet,npets) +!! @author Larissa Reames +!! @author Jeff Beck +!! @author George Gayno + subroutine define_input_grid_grib2(localpet,npets) use netcdf use wgrib2api @@ -863,14 +637,12 @@ subroutine define_input_grid_grib2_gg(localpet,npets) integer :: header_buffer_val = 16384 integer :: id_gridlat, id_gridlon, id_gridrot integer :: id_gridlat_corners, id_gridlon_corners - integer :: i_super, j_super, dim_i_super, dim_j_super logical :: unpack real :: res real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:) real, allocatable :: rlon_corner(:,:),rlat_corner(:,:) - real, allocatable :: rlon_super(:,:),rlat_super(:,:) real, allocatable :: xpts_corner(:,:),ypts_corner(:,:) real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) @@ -973,12 +745,6 @@ subroutine define_input_grid_grib2_gg(localpet,npets) print*,'after gdswzd corner nret ',nret - i_super = 2 * i_input + 1 - j_super = 2 * j_input + 1 - - allocate(rlon_super(i_super,j_super)) - allocate(rlat_super(i_super,j_super)) - if (localpet == 0) then print*,'after gdswzd lat/lon corner 11', rlat_corner(1,1),rlon_corner(1,1) print*,'after gdswzd lat/lon corner ni/11', rlat_corner(ip1_input,1),rlon_corner(ip1_input,1) @@ -1185,460 +951,8 @@ subroutine define_input_grid_grib2_gg(localpet,npets) call MPI_BARRIER(MPI_COMM_WORLD, rc) - end subroutine define_input_grid_grib2_gg - -!> Define input grid object for non-GFS grib2 data. -!! -!! @param [in] localpet ESMF local persistent execution thread -!! @param [in] npets Number of persistent execution threads -!! @author Larissa Reames -!! @author Jeff Beck - subroutine define_input_grid_grib2(localpet, npets) - - use mpi - use grib_mod - use gdswzd_mod - use netcdf - use wgrib2api - use program_setup, only : grib2_file_input_grid, data_dir_input_grid, & - fix_dir_input_grid - implicit none - - character(len=500) :: the_file, temp_file - - integer, intent(in) :: localpet, npets - - integer :: error, extra, i, j, clb(2), cub(2) - - real(esmf_kind_r4), allocatable :: latitude_one_tile(:,:), lat_corners(:,:) - real(esmf_kind_r4), allocatable :: longitude_one_tile(:,:), lon_corners(:,:) - real(esmf_kind_r8) :: lat_target(i_target,j_target), & - lon_target(i_target,j_target) - real(esmf_kind_r8) :: deltalon, dx - integer :: ncid,id_var, id_dim - real(esmf_kind_r8), pointer :: lat_src_ptr(:,:), lon_src_ptr(:,:) - character(len=10000) :: temp_msg - character(len=10) :: temp_num = 'NA' - - integer :: k, jdisc, jgdtn, jpdtn, lugb, lugi - integer :: jids(200), jgdt(200), jpdt(200) - integer :: kgds(200), ni, nj, nret - real :: res - real, allocatable :: rlat(:,:), rlon(:,:),xpts(:,:),ypts(:,:) - real, allocatable :: rlatc(:,:), rlonc(:,:) - logical :: unpack - - type(gribfield) :: gfld - - num_tiles_input_grid = 1 - - !inv_file = "chgres.inv" - the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid - temp_file = trim(fix_dir_input_grid)//"/latlon_grid3.32769.nc" - -!!!!!! - - lugb=12 - call baopenr(lugb,the_file,error) - print*,'after g2 open ',error - - j = 0 ! search at beginning of file - lugi = 0 ! no grib index file - jdisc = 0 ! search for discipline - meterological products - jpdtn = 0 ! search for product definition template number - jgdtn = -1 ! search for grid definition template number - jids = -9999 ! array of values in identification section, set to wildcard - jgdt = -9999 ! array of values in grid definition template 3.m - jpdt = -9999 ! array of values in product definition template 4.n - unpack = .true. ! unpack data - - call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & - unpack, k, gfld, error) - - print*,'after getgb2 temp num ', gfld%igdtnum - print*,'after getgb2 template ', gfld%igdtmpl - - call baclose(lugb,error) - - kgds = 0 - call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, ni, nj, res) - - print*,'after conversion kgds 1 ',ni, nj - print*,'after conversion kgds 2 ',kgds(1:25) - - allocate(rlat(ni,nj)) - allocate(rlon(ni,nj)) - allocate(xpts(ni,nj)) - allocate(ypts(ni,nj)) - - call gdswzd(kgds,0,(ni*nj),-9999.,xpts,ypts,rlon,rlat,nret) - - if (kgds(1) == 205) then - where (rlon > 180.0) rlon = rlon - 360.0 - endif - - print*,'after gdswzd nret ',nret - - if (localpet == 0) then - print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) - print*,'after gdswzd lat/lon ni/11', rlat(ni,1),rlon(ni,1) - print*,'after gdswzd lat/lon 11/nj', rlat(1,nj),rlon(1,nj) - print*,'after gdswzd lat/lon ni/nj', rlat(ni,nj),rlon(ni,nj) - print*,'after gdswzd lat/lon mid ', rlat(ni/2,nj/2),rlon(ni/2,nj/2) - endif - - - deallocate(xpts,ypts) - - allocate(rlatc(ni+1,nj+1)) - allocate(rlonc(ni+1,nj+1)) - allocate(xpts(ni+1,nj+1)) - allocate(ypts(ni+1,nj+1)) - - do j = 1, nj+1 - do i = 1, ni+1 - xpts(i,j) = float(i) - 0.5 - ypts(i,j) = float(j) - 0.5 - enddo - enddo - - call gdswzd(kgds,1,((ni+1)*(nj+1)),-9999.,xpts,ypts,rlonc,rlatc,nret) - if (kgds(1) == 205) then - where (rlonc > 180.0) rlonc = rlonc - 360.0 - endif - - if (localpet) then - print*,'after gdswzd nret ',nret - print*,'after gdswzd corner lat/lon 11', rlatc(1,1),rlonc(1,1) - print*,'after gdswzd corner lat/lon ni/11', rlatc(ni+1,1),rlonc(ni+1,1) - print*,'after gdswzd corner lat/lon 11/nj', rlatc(1,nj+1),rlonc(1,nj+1) - print*,'after gdswzd corner lat/lon ni/nj', rlatc(ni+1,nj+1),rlonc(ni+1,nj+1) - print*,'after gdswzd corner lat/lon mid ', rlatc(ni/2,nj/2),rlonc(ni/2,nj/2) - endif - - - - call abort -!!!!! - - call ESMF_FieldGather(latitude_target_grid, lat_target, rootPet=0, tile=1, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldGather", error) - call ESMF_FieldGather(longitude_target_grid, lon_target, rootPet=0, tile=1, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldGather", error) - - if (localpet==0) then - print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) - error=grb2_mk_inv(the_file,inv_file) - if (error /=0) call error_handler("OPENING GRIB2 FILE",error) - error = grb2_inq(the_file, inv_file,grid_desc=temp_msg) - i = index(temp_msg, "grid_template=") + len("grid_template=") - j = index(temp_msg,":winds(") - temp_num=temp_msg(i:j-1) - endif - call MPI_BARRIER(MPI_COMM_WORLD, error) - call MPI_BCAST(temp_num,10,MPI_CHAR,0,MPI_COMM_WORLD,error) - - ! Wgrib2 can't properly read the lat/lon arrays of data on NCEP rotated lat/lon - ! grids, so read in lat/lon from fixed coordinate file - if (trim(temp_num)=="3.32769" .or. trim(temp_num)=="32769") then - - input_grid_type = "rotated_latlon" - - print*,'model grid - rotated lat lon' - - error=nf90_open(trim(temp_file),nf90_nowrite,ncid) - call netcdf_err(error, 'opening: '//trim(temp_file)) - - error=nf90_inq_dimid(ncid, 'nx', id_dim) - call netcdf_err(error, 'reading nx id' ) - error=nf90_inquire_dimension(ncid,id_dim,len=i_input) - call netcdf_err(error, 'reading nx value' ) - - error=nf90_inq_dimid(ncid, 'ny', id_dim) - call netcdf_err(error, 'reading ny id' ) - error=nf90_inquire_dimension(ncid,id_dim,len=j_input) - call netcdf_err(error, 'reading ny value' ) - - allocate(longitude_one_tile(i_input,j_input)) - allocate(latitude_one_tile(i_input,j_input)) - - error=nf90_inq_varid(ncid, 'gridlat', id_var) - call netcdf_err(error, 'reading field id' ) - error=nf90_get_var(ncid, id_var, latitude_one_tile) - call netcdf_err(error, 'reading field' ) - - error=nf90_inq_varid(ncid, 'gridlon', id_var) - call netcdf_err(error, 'reading field id' ) - error=nf90_get_var(ncid, id_var, longitude_one_tile) - call netcdf_err(error, 'reading field' ) - - print*,'after netcdf lat/lon 11', latitude_one_tile(1,1),longitude_one_tile(1,1) - print*,'after netcdf lat/lon ni/11', & - latitude_one_tile(ni,1),longitude_one_tile(ni,1) - print*,'after netcdf lat/lon 11/nj', & - latitude_one_tile(1,nj),longitude_one_tile(1,nj) - print*,'after netcdf lat/lon ni/nj', & - latitude_one_tile(ni,nj),longitude_one_tile(ni,nj) - print*,'after netcdf lat/lon mid ', & - latitude_one_tile(ni/2,nj/2),longitude_one_tile(ni/2,nj/2) - - elseif (temp_num == "3.0" .or. temp_num == "3.30" .or. temp_num=="30" .or. temp_num == "0") then - - if (temp_num =="3.0" .or. temp_num == "0") then - input_grid_type = "latlon" - print*,'model grid - latlon' - endif - if (temp_num =="3.30" .or. temp_num=='30')then - input_grid_type = "lambert" - print*,'model grid - lambert' - endif - - error = grb2_inq(the_file,inv_file,':PRES:',':surface:',nx=i_input, ny=j_input, & - lat=latitude_one_tile, lon=longitude_one_tile) - if (error /= 1) call error_handler("READING FILE", error) - - - if (localpet==0) print*, "from file lon(1:10,1) = ", longitude_one_tile(1:10,1) - if (localpet==0) print*, "from file lat(1,1:10) = ", latitude_one_tile(1,1:10) - print*,'after wgrib2 lat/lon 11', latitude_one_tile(1,1),longitude_one_tile(1,1) - print*,'after wgrib2 lat/lon ni/11', & - latitude_one_tile(ni,1),longitude_one_tile(ni,1) - print*,'after wgrib2 lat/lon 11/nj', & - latitude_one_tile(1,nj),longitude_one_tile(1,nj) - print*,'after wgrib2 lat/lon ni/nj', & - latitude_one_tile(ni,nj),longitude_one_tile(ni,nj) - print*,'after wgrib2 lat/lon mid ', & - latitude_one_tile(ni/2,nj/2),longitude_one_tile(ni/2,nj/2) - - elseif (temp_num=="NA") then - error = 0 - call error_handler("Grid template number cannot be read from the input file. Please " //& - "check that the wgrib2 executable is in your path.", error) - else - error = 0 - call error_handler("Unknown input file grid template number. Must be one of: " //& - "3, 3.30, 3.32769", error) - endif - - print*,"- I/J DIMENSIONS OF THE INPUT GRID TILES ", i_input, j_input - - ip1_input = i_input + 1 - jp1_input = j_input + 1 - -!----------------------------------------------------------------------- -! Create ESMF grid object for the model grid. -!----------------------------------------------------------------------- - - extra = npets / num_tiles_input_grid - - print*,"- CALL GridCreateNoPeriDim FOR INPUT MODEL GRID" - input_grid = ESMF_GridCreateNoPeriDim(maxIndex=(/i_input,j_input/), & - indexflag=ESMF_INDEX_GLOBAL, & - rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridCreateNoPeriDim", error) - - -!----------------------------------------------------------------------- -! Read the mask and lat/lons. -!----------------------------------------------------------------------- - - print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE." - latitude_input_grid = ESMF_FieldCreate(input_grid, & - typekind=ESMF_TYPEKIND_R8, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - name="input_grid_latitude", & - rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldCreate", error) - - print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE." - longitude_input_grid = ESMF_FieldCreate(input_grid, & - typekind=ESMF_TYPEKIND_R8, & - staggerloc=ESMF_STAGGERLOC_CENTER, & - name="input_grid_longitude", & - rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldCreate", error) - - print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE. " - call ESMF_FieldScatter(latitude_input_grid, real(latitude_one_tile,esmf_kind_r8), rootpet=0, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldScatter", error) - -print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE." - call ESMF_FieldScatter(longitude_input_grid, real(longitude_one_tile,esmf_kind_r8), rootpet=0, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN FieldScatter", error) - - - print*,"- CALL GridAddCoord FOR INPUT GRID." - call ESMF_GridAddCoord(input_grid, & - staggerloc=ESMF_STAGGERLOC_CENTER, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridAddCoord", error) - - - print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." - nullify(lon_src_ptr) - call ESMF_GridGetCoord(input_grid, & - staggerLoc=ESMF_STAGGERLOC_CENTER, & - coordDim=1, & - farrayPtr=lon_src_ptr, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord", error) - - print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." - nullify(lat_src_ptr) - call ESMF_GridGetCoord(input_grid, & - staggerLoc=ESMF_STAGGERLOC_CENTER, & - coordDim=2, & - computationalLBound=clb, & - computationalUBound=cub, & - farrayPtr=lat_src_ptr, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord", error) - - do j = clb(2),cub(2) - do i = clb(1), cub(1) - lon_src_ptr(i,j)=real(longitude_one_tile(i,j),esmf_kind_r8) - lat_src_ptr(i,j)=real(latitude_one_tile(i,j),esmf_kind_r8) -! lon_src_ptr(i,j)=rlon(i,j) -! lat_src_ptr(i,j)=rlat(i,j) - if (abs(rlon(i,j)-longitude_one_tile(i,j)) > 0.1 .or. & - abs(rlat(i,j)-latitude_one_tile(i,j)) > 0.1) then - print*,'lat/lon diff ', rlon(i,j), longitude_one_tile(i,j), & - rlat(i,j),latitude_one_tile(i,j) - endif - enddo - enddo - - print*,"- CALL GridAddCoord FOR INPUT GRID." - call ESMF_GridAddCoord(input_grid, & - staggerloc=ESMF_STAGGERLOC_CORNER, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridAddCoord", error) - - print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD." - nullify(lon_src_ptr) - call ESMF_GridGetCoord(input_grid, & - staggerLoc=ESMF_STAGGERLOC_CORNER, & - coordDim=1, & - farrayPtr=lon_src_ptr, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord", error) - - print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD." - nullify(lat_src_ptr) - call ESMF_GridGetCoord(input_grid, & - staggerLoc=ESMF_STAGGERLOC_CORNER, & - coordDim=2, & - computationalLBound=clb, & - computationalUBound=cub, & - farrayPtr=lat_src_ptr, rc=error) - if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & - call error_handler("IN GridGetCoord", error) - - - ! If we have data on a lat/lon or lambert grid, create staggered coordinates - if(trim(input_grid_type)=="latlon" .or. trim(input_grid_type) == "lambert") then - if (trim(input_grid_type) == "latlon") then - - deltalon = abs(longitude_one_tile(2,1)-longitude_one_tile(1,1)) - do j = clb(2), cub(2) - do i = clb(1), cub(1) - - if (i == ip1_input) then - lon_src_ptr(i,j) = longitude_one_tile(i_input,1) + (0.5_esmf_kind_r8*deltalon) - else - lon_src_ptr(i,j) = longitude_one_tile(i,1) - (0.5_esmf_kind_r8*deltalon) - endif - - if (j == jp1_input) then - lat_src_ptr(i,j) = latitude_one_tile(1,j_input) + (0.5_esmf_kind_r8*deltalon) - else - lat_src_ptr(i,j) = latitude_one_tile(1,j) - (0.5_esmf_kind_r8*deltalon) - endif - - enddo - enddo - else - if (localpet==0) then - !cmdline_msg = "wgrib2 "//trim(the_file)//" -d 1 -grid &> temp2.out" - !call system(cmdline_msg) - !open(4,file="temp2.out") - !do i = 1,6 - ! read(4,"(A)") temp_msg2 - !enddo - !close(4) - print*, trim(temp_msg) - i = index(temp_msg, "Dx ") + len("Dx ") - j = index(temp_msg," m Dy") - read(temp_msg(i:j-1),"(F9.6)") dx - endif - call MPI_BARRIER(MPI_COMM_WORLD,error) - call MPI_BCAST(dx,1,MPI_REAL8,0,MPI_COMM_WORLD,error) - - call get_cell_corners(real(latitude_one_tile,esmf_kind_r8), & - real(longitude_one_tile, esmf_kind_r8), & - lat_src_ptr, lon_src_ptr, dx, clb, cub) - - do j = clb(2),cub(2) - do i = clb(1), cub(1) - if (abs(rlonc(i,j)-lon_src_ptr(i,j)) > 0.1 .or. & - abs(rlatc(i,j)-lat_src_ptr(i,j)) > 0.1) then - print*,'lat/lon corn diff ', rlonc(i,j), lon_src_ptr(i,j), & - rlatc(i,j),lat_src_ptr(i,j) - endif -! lon_src_ptr(i,j)=rlonc(i,j) -! lat_src_ptr(i,j)=rlatc(i,j) - enddo - enddo - - endif - elseif (trim(input_grid_type) == "rotated_latlon") then !Read the corner coords from file - - allocate(lon_corners(ip1_input,jp1_input)) - allocate(lat_corners(ip1_input,jp1_input)) - - error=nf90_inq_varid(ncid, 'gridlon_corners', id_var) - call netcdf_err(error, 'reading field id' ) - error=nf90_get_var(ncid, id_var, lon_corners) - call netcdf_err(error, 'reading field' ) - - error=nf90_inq_varid(ncid, 'gridlat_corners', id_var) - call netcdf_err(error, 'reading field id' ) - error=nf90_get_var(ncid, id_var, lat_corners) - call netcdf_err(error, 'reading field' ) - - do j = clb(2),cub(2) - do i = clb(1), cub(1) - lon_src_ptr(i,j)=real(lon_corners(i,j),esmf_kind_r8) - lat_src_ptr(i,j)=real(lat_corners(i,j),esmf_kind_r8) -! lon_src_ptr(i,j)=rlonc(i,j) -! lat_src_ptr(i,j)=rlatc(i,j) - if (abs(rlonc(i,j)-lon_corners(i,j)) > 0.1 .or. & - abs(rlatc(i,j)-lat_corners(i,j)) > 0.1) then - print*,'lat/lon corn diff ', rlonc(i,j), lon_corners(i,j), & - rlatc(i,j),lat_corners(i,j) - endif - enddo - enddo - - - error= nf90_close(ncid) - endif - - nullify(lon_src_ptr) - nullify(lat_src_ptr) - - deallocate(longitude_one_tile) - deallocate(latitude_one_tile) - end subroutine define_input_grid_grib2 - + !> Setup the esmf grid object for the target grid. !! !! @param [in] localpet ESMF local persistent execution thread @@ -2074,93 +1388,6 @@ subroutine get_model_latlons(mosaic_file, orog_dir, num_tiles, tile, & end subroutine get_model_latlons -!> For grids with equal cell sizes (e.g., lambert conformal), get -!! latitude and longitude of the grid cell corners. -!! -!! @param [in] latitude grid box center latitude -!! @param [in] longitude grid box center longitude -!! @param [inout] latitude_sw latitude of the 'southwest' corner of grid box -!! @param [inout] longitude_sw longitude of the 'southwest' corner of grid box -!! @param [in] dx grid cell side size in meters -!! @param [in] clb lower bounds of indices processed by this mpi task -!! @param [in] cub upper bounds of indices processed by this mpi task -!! @author Larissa Reames -!! @author Jeff Beck - subroutine get_cell_corners( latitude, longitude, latitude_sw, longitude_sw, dx,clb,cub) - implicit none - - real(esmf_kind_r8), intent(in) :: latitude(i_input,j_input) - real(esmf_kind_r8), intent(inout), pointer :: latitude_sw(:,:) - real(esmf_kind_r8), intent(in) :: longitude(i_input, j_input) - real(esmf_kind_r8), intent(inout), pointer :: longitude_sw(:,:) - real(esmf_kind_r8), intent(in) :: dx - - integer, intent(in) :: clb(2), cub(2) - - real(esmf_kind_r8) :: lat1, lon1, lat2, lon2, d, brng - - - real(esmf_kind_r8), parameter :: pi = 3.14159265359 - real(esmf_kind_r8), parameter :: R = 6371200.0 - real(esmf_kind_r8), parameter :: bearingInDegrees = 135.0 - - integer :: i, j - - print*,'in routine get_cell_corners' - - d = sqrt((dx**2.0_esmf_kind_r8)/2.0_esmf_kind_r8) - - do j = clb(2),cub(2) - do i = clb(1), cub(1) - - if (j == jp1_input .and. i == ip1_input) then - lat1 = latitude(i_input,j_input) * ( pi / 180.0_esmf_kind_r8 ) - lon1 = longitude(i_input,j_input) * ( pi / 180.0_esmf_kind_r8 ) - brng = 315.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8 - lat2 = asin( sin( lat1 ) * cos( d / R ) + cos( lat1 ) * sin( d / R ) * cos( brng ) ); - lon2= lon1 + atan2( sin( brng ) * sin( d / R ) * cos( lat1 ), cos( d / R ) - sin( lat1 ) * sin( lat2 ) ); - latitude_sw(ip1_input,jp1_input) = lat2 * 180.0_esmf_kind_r8 / pi - longitude_sw(ip1_input,jp1_input) = lon2 * 180.0_esmf_kind_r8 / pi - cycle - endif - - if (i == ip1_input) then - brng = 225.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8 - lat1 = latitude(i_input,j) * ( pi / 180.0_esmf_kind_r8 ) - lon1 = longitude(i_input,j) * ( pi / 180.0_esmf_kind_r8 ) - lat2 = asin( sin( lat1 ) * cos( d / R ) + cos( lat1 ) * sin( d / R ) * cos( brng ) ); - lon2= lon1 + atan2( sin( brng ) * sin( d / R ) * cos( lat1 ), cos( d / R ) - sin( lat1 ) * sin( lat2 ) ); - latitude_sw(ip1_input,j) = lat2 * 180.0_esmf_kind_r8 / pi - longitude_sw(ip1_input,j) = lon2 * 180.0_esmf_kind_r8 / pi - cycle - endif - - if (j == jp1_input) then - brng = 45.0_esmf_kind_r8 * pi / 180.0_esmf_kind_r8 - lat1 = latitude(i,j_input) * ( pi / 180.0_esmf_kind_r8 ) - lon1 = longitude(i,j_input) * ( pi / 180.0_esmf_kind_r8 ) - lat2 = asin( sin( lat1 ) * cos( d / R ) + cos( lat1 ) * sin( d / R ) * cos( brng ) ); - lon2= lon1 + atan2( sin( brng ) * sin( d / R ) * cos( lat1 ), cos( d / R ) - sin( lat1 ) * sin( lat2 ) ); - latitude_sw(i,jp1_input) = lat2 * 180.0_esmf_kind_r8 / pi - longitude_sw(i,jp1_input) = lon2 * 180.0_esmf_kind_r8 / pi - cycle - endif - - lat1 = latitude(i,j) * ( pi / 180.0_esmf_kind_r8 ) - lon1 = longitude(i,j) * ( pi / 180.0_esmf_kind_r8 ) - - brng = bearingInDegrees * ( pi / 180.0_esmf_kind_r8 ); - lat2 = asin( sin( lat1 ) * cos( d / R ) + cos( lat1 ) * sin( d / R ) * cos( brng ) ); - lon2= lon1 + atan2( sin( brng ) * sin( d / R ) * cos( lat1 ), cos( d / R ) - sin( lat1 ) * sin( lat2 ) ); - - latitude_sw(i,j) = lat2 * 180.0_esmf_kind_r8 / pi - longitude_sw(i,j) = lon2 * 180.0_esmf_kind_r8 / pi - - enddo - enddo - - end subroutine get_cell_corners - !> Read the model land mask and terrain for a single tile !! from the orography file. !! From a4d9699b076fd4bdb392ccb393c4aa3070c0d6bf Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 4 Feb 2022 20:40:21 +0000 Subject: [PATCH 08/14] Some cleanup to model_grid.F90 --- sorc/chgres_cube.fd/model_grid.F90 | 66 ++++++++++++++++++------------ 1 file changed, 40 insertions(+), 26 deletions(-) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 015cf9bde..95bb5debe 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -607,7 +607,7 @@ subroutine define_input_grid_mosaic(localpet, npets) end subroutine define_input_grid_mosaic -!> Define input grid object for non-GFS grib2 data. +!> Define input grid object for grib2 input data. !! !! @param [in] localpet ESMF local persistent execution thread !! @param [in] npets Number of persistent execution threads @@ -661,48 +661,44 @@ subroutine define_input_grid_grib2(localpet,npets) lugb=12 + print*,"- OPEN AND READ INPUT DATA GRIB2 FILE: ", trim(the_file) call baopenr(lugb,the_file,rc) - print*,'after g2 open ',rc - - j = 0 ! search at beginning of file - lugi = 0 ! no grib index file - jdisc = 0 ! search for discipline - meterological products - jpdtn = 0 ! search for product definition template number - jgdtn = -1 ! search for grid definition template number - jids = -9999 ! array of values in identification section, set to wildcard - jgdt = -9999 ! array of values in grid definition template 3.m - jpdt = -9999 ! array of values in product definition template 4.n - unpack = .true. ! unpack data + if (rc /= 0) call error_handler("OPENING FILE", rc) + +! Read the first record and get the grid definition template. + + j = 0 ! Search at beginning of file + lugi = 0 ! No grib index file + jdisc = -1 ! Search for any discipline + jpdtn = -1 ! Search for any product definition template number + jgdtn = -1 ! Search for any grid definition template number + jids = -9999 ! Array of values in identification section, set to wildcard. + jgdt = -9999 ! Array of values in grid definition template, set to wildcard. + jpdt = -9999 ! Array of values in product definition template, set to wildcard. + unpack = .false. ! unpack data call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & unpack, k, gfld, rc) - - if (localpet == 0)then - print*,'after getgb2 temp num ', gfld%igdtnum - print*,'after getgb2 template ', gfld%igdtmpl - endif + if (rc /= 0) call error_handler("DEGRIBBING INPUT FILE.", rc) call baclose(lugb,rc) if (gfld%igdtnum == 0) then + print*,"- INPUT DATA ON LAT/LON GRID." input_grid_type = 'latlon' elseif (gfld%igdtnum == 30) then + print*,"- INPUT DATA ON LAMBERT CONFORMAL GRID." input_grid_type = 'lambert' elseif (gfld%igdtnum == 32769) then + print*,"- INPUT DATA ON ROTATED LAT/LON GRID." input_grid_type = 'rotated_latlon' else - print*,'grid not supported' - call abort + call error_handler("INPUT GRID TEMPLATE NOT SUPPORTED.", 2) endif kgds = 0 call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, i_input, j_input, res) - if (localpet == 0) then - print*,'after conversion kgds ni/nj ',i_input, j_input - print*,'after conversion kgds ',kgds(1:25) - endif - ip1_input = i_input + 1 jp1_input = j_input + 1 @@ -722,9 +718,14 @@ subroutine define_input_grid_grib2(localpet,npets) enddo enddo + print*,"- COMPUTE GRID CELL CENTER COORDINATES." call gdswzd(kgds,1,(i_input*j_input),-9999.,xpts,ypts,rlon,rlat,nret) - print*,'after gdswzd nret ',nret + if (nret /= (i_input*j_input)) then + call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2) + endif + + deallocate(xpts, ypts) if (localpet == 0) then print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) @@ -741,9 +742,14 @@ subroutine define_input_grid_grib2(localpet,npets) enddo enddo + print*,"- COMPUTE GRID CELL CORNER COORDINATES." call gdswzd(kgds,1,(ip1_input*jp1_input),-9999.,xpts_corner,ypts_corner,rlon_corner,rlat_corner,nret) - print*,'after gdswzd corner nret ',nret + if (nret /= (ip1_input*jp1_input)) then + call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2) + endif + + deallocate(xpts_corner, ypts_corner) if (localpet == 0) then print*,'after gdswzd lat/lon corner 11', rlat_corner(1,1),rlon_corner(1,1) @@ -859,6 +865,8 @@ subroutine define_input_grid_grib2(localpet,npets) latitude = rlat longitude = rlon + deallocate (rlat, rlon) + print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE." call ESMF_FieldScatter(longitude_input_grid, longitude, rootpet=0, rc=rc) if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) & @@ -903,12 +911,16 @@ subroutine define_input_grid_grib2(localpet,npets) enddo enddo + deallocate(latitude, longitude) + allocate(latitude_corner(ip1_input,jp1_input)) allocate(longitude_corner(ip1_input,jp1_input)) latitude_corner = rlat_corner longitude_corner = rlon_corner + deallocate (rlat_corner, rlon_corner) + print*,"- CALL GridAddCoord FOR INPUT GRID." call ESMF_GridAddCoord(input_grid, & staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc) @@ -943,6 +955,8 @@ subroutine define_input_grid_grib2(localpet,npets) enddo enddo + deallocate(latitude_corner, longitude_corner) + if (localpet == 0) then print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) rc=grb2_mk_inv(the_file,inv_file) From b498e5d48eae57504a2a4a029dca45df5cdf5e71 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 4 Feb 2022 21:09:27 +0000 Subject: [PATCH 09/14] Remove unused variable used to store the path and name of the RAP grid coordinate file. Part of #591. --- sorc/chgres_cube.fd/input_data.F90 | 5 +---- sorc/chgres_cube.fd/program_setup.F90 | 5 ----- ush/chgres_cube.sh | 4 +--- 3 files changed, 2 insertions(+), 12 deletions(-) diff --git a/sorc/chgres_cube.fd/input_data.F90 b/sorc/chgres_cube.fd/input_data.F90 index dab617f35..03a5850ea 100644 --- a/sorc/chgres_cube.fd/input_data.F90 +++ b/sorc/chgres_cube.fd/input_data.F90 @@ -6144,7 +6144,7 @@ subroutine read_winds(file,inv,u,v,localpet) use wgrib2api use netcdf - use program_setup, only : get_var_cond, fix_dir_input_grid + use program_setup, only : get_var_cond use model_grid, only : input_grid_type implicit none @@ -6165,7 +6165,6 @@ subroutine read_winds(file,inv,u,v,localpet) character(len=20) :: vname character(len=50) :: method_u, method_v - character(len=250) :: file_coord character(len=10000) :: temp_msg d2r=acos(-1.0_esmf_kind_r8) / 180.0_esmf_kind_r8 @@ -6176,8 +6175,6 @@ subroutine read_winds(file,inv,u,v,localpet) allocate(u(0,0,0)) allocate(v(0,0,0)) endif - - file_coord = trim(fix_dir_input_grid)//"/latlon_grid3.32769.nc" vname = "u" call get_var_cond(vname,this_miss_var_method=method_u, this_miss_var_value=value_u, & diff --git a/sorc/chgres_cube.fd/program_setup.F90 b/sorc/chgres_cube.fd/program_setup.F90 index 5e854a2c3..722e701a7 100644 --- a/sorc/chgres_cube.fd/program_setup.F90 +++ b/sorc/chgres_cube.fd/program_setup.F90 @@ -53,10 +53,6 @@ module program_setup !! 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=500), public :: fix_dir_input_grid = "NULL" !< Directory containing files of latitude and - !! and longitude for certain GRIB2 input data. - - integer, parameter, public :: max_tracers=100 !< Maximum number of atmospheric tracers processed. integer, public :: num_tracers !< Number of atmospheric tracers to be processed. integer, public :: num_tracers_input !< Number of atmospheric tracers in input file. @@ -193,7 +189,6 @@ subroutine read_setup_namelist(filename) tracers_input, & halo_bndy, & halo_blend, & - fix_dir_input_grid, & nsoill_out, & thomp_mp_climo_file diff --git a/ush/chgres_cube.sh b/ush/chgres_cube.sh index c60eba605..bcba5f5d1 100755 --- a/ush/chgres_cube.sh +++ b/ush/chgres_cube.sh @@ -24,8 +24,7 @@ CRES=${CRES:-96} # FIXufs - Location of ufs_utils root fixed data directory. # FIXfv3 - Location of target grid orography and 'grid' files. # FIXsfc - Location of target grid surface climatological files. -# FIXam - Location of vertical coordinate definition file for target grid -# and RAP grib2 input grid lat/lon definition file. +# FIXam - Location of vertical coordinate definition file for target grid. #---------------------------------------------------------------------------- ufs_ver=${ufs_ver:-v1.0.0} @@ -249,7 +248,6 @@ cat << EOF > ./fort.41 mosaic_file_target_grid="${MOSAIC_FILE_TARGET_GRID}" fix_dir_target_grid="${FIXsfc}" orog_dir_target_grid="${FIXfv3}" - fix_dir_input_grid="${FIXam}" orog_files_target_grid="${OROG_FILES_TARGET_GRID}" vcoord_file_target_grid="${VCOORD_FILE}" mosaic_file_input_grid="${MOSAIC_FILE_INPUT_GRID}" From fbffe3c906f0a32dd42c620e84eb38b9e3ad39fe Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 24 Mar 2022 15:41:02 +0000 Subject: [PATCH 10/14] Remove all dependencies on wgrib2 from the build. Fixes #591. --- CMakeLists.txt | 1 - README.md | 1 - cmake/Findwgrib2.cmake | 41 ----------------------- modulefiles/build.cheyenne.intel | 1 - modulefiles/build.hera.gnu.lua | 3 -- modulefiles/build.hera.intel.lua | 3 -- modulefiles/build.jet.intel.lua | 3 -- modulefiles/build.odin.intel | 3 -- modulefiles/build.orion.intel.lua | 3 -- modulefiles/build.s4.intel | 1 - modulefiles/build.stampede.intel | 1 - modulefiles/build.wcoss_cray.intel | 1 - modulefiles/build.wcoss_dell_p3.intel.lua | 3 -- sorc/chgres_cube.fd/CMakeLists.txt | 1 - sorc/chgres_cube.fd/input_data.F90 | 3 +- sorc/chgres_cube.fd/model_grid.F90 | 12 ------- 16 files changed, 1 insertion(+), 80 deletions(-) delete mode 100644 cmake/Findwgrib2.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index 38bcfb5ed..68869d5c2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -70,7 +70,6 @@ find_package(sigio 2.3.0 REQUIRED) find_package(sp 2.3.3 REQUIRED) find_package(ip 3.3.3 REQUIRED) find_package(g2 3.4.3 REQUIRED) -find_package(wgrib2 2.0.8 REQUIRED) find_package(sigio 2.3.0 REQUIRED) # EMC requires executables in ./exec diff --git a/README.md b/README.md index 3a192662c..1241c9b6d 100644 --- a/README.md +++ b/README.md @@ -40,7 +40,6 @@ This package uses the [hpc-stack](https://github.com/NOAA-EMC/hpc-stack) for the - [NCEPLIBS-sp](https://github.com/NOAA-EMC/NCEPLIBS-sp) - [NCEPLIBS-ip](https://github.com/NOAA-EMC/NCEPLIBS-ip) - [NCEPLIBS-g2](https://github.com/NOAA-EMC/NCEPLIBS-g2) - - [NCEPLIBS-wgrib2](https://github.com/NOAA-EMC/NCEPLIBS-wgrib2) And for the following third party libraries: diff --git a/cmake/Findwgrib2.cmake b/cmake/Findwgrib2.cmake deleted file mode 100644 index 68b8833c4..000000000 --- a/cmake/Findwgrib2.cmake +++ /dev/null @@ -1,41 +0,0 @@ -# This module produces the target wgrib2::wgrib2 - -find_path(WGRIB2_INCLUDES wgrib2api.mod) -find_library(WGRIB2_LIB libwgrib2.a) -find_library(WGRIB2_API_LIB libwgrib2_api.a) - -add_library(wgrib2::wgrib2 UNKNOWN IMPORTED) - -# Library builds are different between CMake build and the make build. -# libwgrib2_api.a is only necessary in the CMake build and must come first when linking -if(WGRIB2_API_LIB) - # CMake build. Need both. - set(first_lib ${WGRIB2_API_LIB}) - set(second_lib ${WGRIB2_LIB}) -else() - # Makefile build. Only need libwgrib2.a - set(first_lib ${WGRIB2_LIB}) - set(second_lib "") -endif() - -set_target_properties(wgrib2::wgrib2 PROPERTIES - IMPORTED_LOCATION "${first_lib}" - INTERFACE_INCLUDE_DIRECTORIES "${WGRIB2_INCLUDES}" - INTERFACE_LINK_LIBRARIES "${second_lib}") - -set(WGRIB2_LIBRARIES "${first_lib}" "${second_lib}") - -find_program(WGRIB2_EXE wgrib2) -execute_process(COMMAND ${WGRIB2_EXE} --version OUTPUT_VARIABLE version_str) - -# Wgrib2 changed how it output --version from "v0.x.y.z" to "vx.y.z" starting in wgrib2 3.0 -if(version_str MATCHES "^v0.*") - string(SUBSTRING "${version_str}" 3 5 version) -else() - string(SUBSTRING "${version_str}" 1 5 version) -endif() - -find_package_handle_standard_args(wgrib2 - REQUIRED_VARS WGRIB2_LIBRARIES WGRIB2_INCLUDES WGRIB2_EXE - VERSION_VAR version - ) diff --git a/modulefiles/build.cheyenne.intel b/modulefiles/build.cheyenne.intel index 570a424f4..c1971328f 100644 --- a/modulefiles/build.cheyenne.intel +++ b/modulefiles/build.cheyenne.intel @@ -20,7 +20,6 @@ module load w3nco/2.4.1 module load sigio/2.3.2 module load sfcio/1.4.1 -module load wgrib2/2.0.8 module load netcdf/4.7.4 setenv ESMFMKFILE /glade/p/ral/jntp/GMTB/tools/NCEPLIBS-ufs-v2.0.0/intel-19.1.1/mpt-2.19/lib64/esmf.mk diff --git a/modulefiles/build.hera.gnu.lua b/modulefiles/build.hera.gnu.lua index 90bca87ae..6b57b13c3 100644 --- a/modulefiles/build.hera.gnu.lua +++ b/modulefiles/build.hera.gnu.lua @@ -49,9 +49,6 @@ load(pathJoin("sfcio", sfcio_ver)) sigio_ver=os.getenv("sigio_ver") or "2.3.2" load(pathJoin("sigio", sigio_ver)) -wgrib2_ver=os.getenv("wgrib2_ver") or "2.0.8" -load(pathJoin("wgrib2", wgrib2_ver)) - nccmp_ver=os.getenv("nccmp_ver") or "1.8.7.0" load(pathJoin("nccmp", nccmp_ver)) diff --git a/modulefiles/build.hera.intel.lua b/modulefiles/build.hera.intel.lua index 48335e748..036b47fc6 100644 --- a/modulefiles/build.hera.intel.lua +++ b/modulefiles/build.hera.intel.lua @@ -43,9 +43,6 @@ load(pathJoin("sfcio", sfcio_ver)) sigio_ver=os.getenv("sigio_ver") or "2.3.2" load(pathJoin("sigio", sigio_ver)) -wgrib2_ver=os.getenv("wgrib2_ver") or "2.0.8" -load(pathJoin("wgrib2", wgrib2_ver)) - zlib_ver=os.getenv("zlib_ver") or "1.2.11" load(pathJoin("zlib", zlib_ver)) diff --git a/modulefiles/build.jet.intel.lua b/modulefiles/build.jet.intel.lua index 50f23bb96..2a18d626d 100644 --- a/modulefiles/build.jet.intel.lua +++ b/modulefiles/build.jet.intel.lua @@ -55,9 +55,6 @@ load(pathJoin("nemsio", nemsio_ver)) g2_ver=os.getenv("g2_ver") or "3.4.3" load(pathJoin("g2", g2_ver)) -wgrib2_ver=os.getenv("wgrib2_ver") or "2.0.8" -load(pathJoin("wgrib2", wgrib2_ver)) - prod_util_ver=os.getenv("prod_util_ver") or "1.2.2" load(pathJoin("prod_util", prod_util_ver)) diff --git a/modulefiles/build.odin.intel b/modulefiles/build.odin.intel index 362aa3a84..7904f665c 100644 --- a/modulefiles/build.odin.intel +++ b/modulefiles/build.odin.intel @@ -25,7 +25,6 @@ module load sigio module load sfcio module load nemsio module load g2 -module load wgrib2 #module load esmf/8.0.0 #setenv ESMFMKFILE /oldscratch/ywang/external/NCEPLIBS_SRW/lib64/esmf.mk @@ -33,5 +32,3 @@ setenv ESMFMKFILE /oldscratch/ywang/external/NCEPLIBS_SRWv2.0/lib64/esmf.mk setenv CMAKE_Fortran_COMPILER ftn setenv CMAKE_C_COMPILER cc - -#setenv WGRIB2_ROOT /oldscratch/ywang/external/NCEPLIBS_SRWv2.0/wgrib2-2.0.8 diff --git a/modulefiles/build.orion.intel.lua b/modulefiles/build.orion.intel.lua index d641cdaf5..fb1ec0569 100644 --- a/modulefiles/build.orion.intel.lua +++ b/modulefiles/build.orion.intel.lua @@ -40,9 +40,6 @@ load(pathJoin("sfcio", sfcio_ver)) sigio_ver=os.getenv("sigio_ver") or "2.3.2" load(pathJoin("sigio", sigio_ver)) -wgrib2_ver=os.getenv("wgrib2_ver") or "2.0.8" -load(pathJoin("wgrib2", wgrib2_ver)) - zlib_ver=os.getenv("zlib_ver") or "1.2.11" load(pathJoin("zlib", zlib_ver)) diff --git a/modulefiles/build.s4.intel b/modulefiles/build.s4.intel index 8a9a0b0eb..37cdf2a0c 100644 --- a/modulefiles/build.s4.intel +++ b/modulefiles/build.s4.intel @@ -16,7 +16,6 @@ module load sp/2.3.3 module load w3nco/2.4.1 module load sfcio/1.4.1 module load sigio/2.3.2 -module load wgrib2/2.0.8 module load jasper/2.0.22 module load zlib/1.2.11 diff --git a/modulefiles/build.stampede.intel b/modulefiles/build.stampede.intel index 5d5829507..e91e253ca 100644 --- a/modulefiles/build.stampede.intel +++ b/modulefiles/build.stampede.intel @@ -32,7 +32,6 @@ module load sigio module load sfcio module load nemsio module load g2 -module load wgrib2 #setenv ESMFMKFILE /work/00315/tg455890/stampede2/regional_fv3/NCEPLIBS_SRWv2.0/lib64/esmf.mk diff --git a/modulefiles/build.wcoss_cray.intel b/modulefiles/build.wcoss_cray.intel index 7ae0c107f..d2c7c438e 100644 --- a/modulefiles/build.wcoss_cray.intel +++ b/modulefiles/build.wcoss_cray.intel @@ -24,7 +24,6 @@ module load sp/2.3.3 module load w3nco/2.4.1 module load sfcio/1.4.1 module load sigio/2.3.2 -module load wgrib2/2.0.8 setenv ZLIB_ROOT /usrx/local/prod/zlib/1.2.7/intel/haswell setenv PNG_ROOT /usrx/local/prod/png/1.2.49/intel/haswell diff --git a/modulefiles/build.wcoss_dell_p3.intel.lua b/modulefiles/build.wcoss_dell_p3.intel.lua index 0c0f3955e..9d7408bc8 100644 --- a/modulefiles/build.wcoss_dell_p3.intel.lua +++ b/modulefiles/build.wcoss_dell_p3.intel.lua @@ -64,9 +64,6 @@ load(pathJoin("sfcio", sfcio_ver)) sigio_ver=os.getenv("sigio_ver") or "2.3.2" load(pathJoin("sigio", sigio_ver)) -wgrib2_ver=os.getenv("wgrib2_ver") or "2.0.8" -load(pathJoin("wgrib2", wgrib2_ver)) - prepend_path("MODULEPATH", "/usrx/local/dev/modulefiles") prod_util_ver=os.getenv("prod_util_ver") or "1.1.3" diff --git a/sorc/chgres_cube.fd/CMakeLists.txt b/sorc/chgres_cube.fd/CMakeLists.txt index 1b086207d..81300a6a4 100644 --- a/sorc/chgres_cube.fd/CMakeLists.txt +++ b/sorc/chgres_cube.fd/CMakeLists.txt @@ -51,7 +51,6 @@ target_link_libraries( sp::sp_d w3nco::w3nco_d esmf - wgrib2::wgrib2 MPI::MPI_Fortran NetCDF::NetCDF_Fortran) diff --git a/sorc/chgres_cube.fd/input_data.F90 b/sorc/chgres_cube.fd/input_data.F90 index 420b92955..ed7f60cd0 100644 --- a/sorc/chgres_cube.fd/input_data.F90 +++ b/sorc/chgres_cube.fd/input_data.F90 @@ -41,8 +41,7 @@ module input_data ip1_input, jp1_input, & num_tiles_input_grid, & latitude_input_grid, & - longitude_input_grid, & - inv_file + longitude_input_grid implicit none diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 06dafc4ca..8eb5a3cf7 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -17,8 +17,6 @@ module model_grid character(len=5), allocatable, public :: tiles_target_grid(:) !< Tile names of target grid. - character(len=10), public :: inv_file = "chgres.inv" - !< wgrib2 inventory file character(len=50), public :: input_grid_type = "latlon" !< map projection of input grid @@ -617,7 +615,6 @@ end subroutine define_input_grid_mosaic subroutine define_input_grid_grib2(localpet,npets) use netcdf - use wgrib2api use mpi use grib_mod use gdswzd_mod @@ -957,15 +954,6 @@ subroutine define_input_grid_grib2(localpet,npets) deallocate(latitude_corner, longitude_corner) -! get rid of this? - if (localpet == 0) then - print*,'- OPEN AND INVENTORY GRIB2 FILE: ',trim(the_file) - rc=grb2_mk_inv(the_file,inv_file) - if (rc /=0) call error_handler("OPENING GRIB2 FILE using wgrib2",rc) - endif - - call MPI_BARRIER(MPI_COMM_WORLD, rc) - end subroutine define_input_grid_grib2 !> Setup the esmf grid object for the target grid. From 73248cfe3230c20e29616c8bf6d9c28ff0514a48 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 24 Mar 2022 16:01:00 +0000 Subject: [PATCH 11/14] Remove wgrib2 from workflow 'yml' files. Fixes #591. --- .github/workflows/debug-docs-test_coverage.yml | 3 --- .github/workflows/intel.yml | 3 --- .github/workflows/linux-mac-nceplibs-mpi.yml | 3 --- .github/workflows/netcdf-versions.yml | 3 --- 4 files changed, 12 deletions(-) diff --git a/.github/workflows/debug-docs-test_coverage.yml b/.github/workflows/debug-docs-test_coverage.yml index 96face7b7..6d953cd95 100644 --- a/.github/workflows/debug-docs-test_coverage.yml +++ b/.github/workflows/debug-docs-test_coverage.yml @@ -97,10 +97,7 @@ jobs: export CC=mpicc export CXX=mpicxx export FC=mpifort - # Findwgrib2 in module form does not search -version - # as NCEPLIBS installs it export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:~/jasper/lib;~/jasper/lib64" - export wgrib2_ROOT=`find ~/nceplibs -type d -maxdepth 1 -iname "wgrib2*"` cmake .. -DCMAKE_PREFIX_PATH='~/jasper;~/nceplibs' -DCMAKE_BUILD_TYPE=Debug -DENABLE_DOCS=On -DCMAKE_Fortran_FLAGS="-g -fprofile-arcs -ftest-coverage -O0" make -j2 diff --git a/.github/workflows/intel.yml b/.github/workflows/intel.yml index e1743a8a1..240003849 100644 --- a/.github/workflows/intel.yml +++ b/.github/workflows/intel.yml @@ -160,10 +160,7 @@ jobs: export ESMFMKFILE=~/esmf/lib/esmf.mk cd ufs_utils mkdir build && cd build - # Findwgrib2 in module form does not search -version - # as NCEPLIBS installs it export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:~/jasper/lib;~/jasper/lib64" - export wgrib2_ROOT=`find ~/nceplibs -type d -maxdepth 1 -iname "wgrib2*"` cmake .. -DCMAKE_BUILD_TYPE=Debug -DCMAKE_PREFIX_PATH='~;~/jasper;~/nceplibs;~/netcdf' make -j2 diff --git a/.github/workflows/linux-mac-nceplibs-mpi.yml b/.github/workflows/linux-mac-nceplibs-mpi.yml index 9250598fa..68d9da6ae 100644 --- a/.github/workflows/linux-mac-nceplibs-mpi.yml +++ b/.github/workflows/linux-mac-nceplibs-mpi.yml @@ -214,11 +214,8 @@ jobs: export CC=mpicc export CXX=mpicxx export FC=mpifort - # Findwgrib2 in module form does not search -version - # as NCEPLIBS installs it export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:~/jasper/lib;~/jasper/lib64" export DYLD_LIBRARY_PATH="${LD_LIBRARY_PATH}:~/jasper/lib;~/jasper/lib64" - export wgrib2_ROOT=`find ~/nceplibs -type d -maxdepth 1 -iname "wgrib2*"` cmake .. -DCMAKE_PREFIX_PATH='~/jasper;~/nceplibs;~/netcdf' make -j2 diff --git a/.github/workflows/netcdf-versions.yml b/.github/workflows/netcdf-versions.yml index 2040e6b72..98ecd746a 100644 --- a/.github/workflows/netcdf-versions.yml +++ b/.github/workflows/netcdf-versions.yml @@ -148,10 +148,7 @@ jobs: export CC=mpicc export CXX=mpicxx export FC=mpifort - # Findwgrib2 in module form does not search -version - # as NCEPLIBS installs it export LD_LIBRARY_PATH="${LD_LIBRARY_PATH}:~/jasper/lib;~/jasper/lib64" - export wgrib2_ROOT=`find ~/nceplibs -type d -maxdepth 1 -iname "wgrib2*"` cmake .. -DCMAKE_PREFIX_PATH='~/jasper;~/nceplibs;~/netcdf' -DCMAKE_BUILD_TYPE=Debug make -j2 From 7604cf81723f9185df00b1aae8997b9289e3149c Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 28 Mar 2022 19:56:23 +0000 Subject: [PATCH 12/14] Fix pole flip bug in gdt_to_gds for GFS GRIB2 data. Fixes #591. --- sorc/chgres_cube.fd/model_grid.F90 | 36 +++++++++++++----------------- 1 file changed, 16 insertions(+), 20 deletions(-) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index 8eb5a3cf7..b790266a6 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -1511,16 +1511,16 @@ subroutine cleanup_input_target_grid_data end subroutine cleanup_input_target_grid_data -!> Convert the grib2 grid description template to -!! to the grib2 grid description sections. +!> Convert the GRIB2 grid description template to +!! to the GRIB1 grid description section. !! -!! @param [in] igdtnum grib2 grid description template number -!! @param [in] igdstmpl length of grib2 grid description template -!! @param [in] igdtlen array of grib2 grid description template octets. -!! @param [out] kgds array of grib1 grid description octets. -!! @param [out] ni i-dimension of grid -!! @param [out] nj j-dimension of grid -!! @param [out] res resolution of grid in km +!! @param [in] igdtnum GRIB2 grid description template number. +!! @param [in] igdstmpl Length of grib2 grid description template. +!! @param [in] igdtlen Array of GRIB2 grid description template octets. +!! @param [out] kgds Array of GRIB1 grid description octets. +!! @param [out] ni I-dimension of grid. +!! @param [out] nj J-dimension of grid. +!! @param [out] res Resolution of grid in km. !! @author George Gayno NCEP/EMC subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) @@ -1615,13 +1615,12 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) iscale=igdstmpl(10)*igdstmpl(11) if (iscale == 0) iscale = 1e6 - kgds(1)=0 ! oct 6 + kgds(1)=0 ! oct 6, data representation type. 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(7)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point + 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 @@ -1629,11 +1628,10 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) 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(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point - kgds(4)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point + kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point - kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, di - kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, dj + kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, "i" resolution. + kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, "j" resolution. kgds(11) = 0 ! oct 28, scan mode if (btest(igdstmpl(19),7)) kgds(11) = 128 @@ -1645,10 +1643,8 @@ subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res) kgds(20)=255 ! oct 5, used for thinned grids, set to 255 else - print*,'grid not defined', igdtnum - - call abort - + + call error_handler("UNRECOGNIZED INPUT GRID TYPE ", 1) endif From 9f66174294686efcfd79569b05b199ef10fb8f5a Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 29 Mar 2022 15:02:29 +0000 Subject: [PATCH 13/14] Add diagnostic print of lat/lon corner/center differences. Fixes #591. --- sorc/chgres_cube.fd/model_grid.F90 | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index b790266a6..f4d92a72e 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -634,12 +634,14 @@ subroutine define_input_grid_grib2(localpet,npets) integer :: header_buffer_val = 16384 integer :: id_gridlat, id_gridlon, id_gridrot integer :: id_gridlat_corners, id_gridlon_corners + integer :: id_gridlat_diff, id_gridlon_diff logical :: unpack real :: res real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:) real, allocatable :: rlon_corner(:,:),rlat_corner(:,:) + real, allocatable :: rlon_diff(:,:),rlat_diff(:,:) real, allocatable :: xpts_corner(:,:),ypts_corner(:,:) real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) @@ -701,6 +703,8 @@ subroutine define_input_grid_grib2(localpet,npets) allocate(rlat(i_input,j_input)) allocate(rlon(i_input,j_input)) + allocate(rlat_diff(i_input,j_input)) + allocate(rlon_diff(i_input,j_input)) allocate(xpts(i_input,j_input)) allocate(ypts(i_input,j_input)) allocate(rlat_corner(ip1_input,jp1_input)) @@ -761,6 +765,13 @@ subroutine define_input_grid_grib2(localpet,npets) where(rlon_corner > 180.0) rlon_corner = rlon_corner - 360.0 endif + do j = 1, j_input + do i = 1, i_input + rlat_diff(i,j) = rlat_corner(i,j) - rlat(i,j) + rlon_diff(i,j) = rlon_corner(i,j) - rlon(i,j) + enddo + enddo + ncid = 34 error = nf90_create("./latlon.nc", IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), & ncid, initialsize=0, chunksize=fsize) @@ -793,6 +804,12 @@ subroutine define_input_grid_grib2(localpet,npets) error = nf90_def_var(ncid, 'gridlon_corners', NF90_FLOAT, (/dim_ip1,dim_jp1/), id_gridlon_corners) call netcdf_err(error, 'DEFINING gridlon_corners' ) + error = nf90_def_var(ncid, 'gridlat_corner_minus_cent', NF90_FLOAT, (/dim_i,dim_j/), id_gridlat_diff) + call netcdf_err(error, 'DEFINING gridlat_corner_minus_cent' ) + + error = nf90_def_var(ncid, 'gridlon_corner_minus_cent', NF90_FLOAT, (/dim_i,dim_j/), id_gridlon_diff) + call netcdf_err(error, 'DEFINING gridlon_corner_minus_cent' ) + error = nf90_enddef(ncid, header_buffer_val,4,0,4) call netcdf_err(error, 'DEFINING HEADER' ) @@ -802,6 +819,12 @@ subroutine define_input_grid_grib2(localpet,npets) error = nf90_put_var(ncid, id_gridlon, real(rlon,4) ) call netcdf_err(error, 'writing gridlon' ) + error = nf90_put_var(ncid, id_gridlat_diff, real(rlat_diff,4) ) + call netcdf_err(error, 'writing gridlat_diff' ) + + error = nf90_put_var(ncid, id_gridlon_diff, real(rlon_diff,4) ) + call netcdf_err(error, 'writing gridlon_diff' ) + error = nf90_put_var(ncid, id_gridlat_corners, real(rlat_corner,4) ) call netcdf_err(error, 'writing gridlat_corner' ) From 4be3bfd59899c194a0e9984574538eacd0163768 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 19 Apr 2022 15:28:59 +0000 Subject: [PATCH 14/14] Remove some unused variables and the write of the diagnostic file containing grid point lat/lons. Fixes #591. --- sorc/chgres_cube.fd/model_grid.F90 | 125 +++------------------------- sorc/chgres_cube.fd/static_data.F90 | 3 +- sorc/chgres_cube.fd/surface.F90 | 2 +- 3 files changed, 15 insertions(+), 115 deletions(-) diff --git a/sorc/chgres_cube.fd/model_grid.F90 b/sorc/chgres_cube.fd/model_grid.F90 index f4d92a72e..cdd91336b 100644 --- a/sorc/chgres_cube.fd/model_grid.F90 +++ b/sorc/chgres_cube.fd/model_grid.F90 @@ -614,35 +614,27 @@ end subroutine define_input_grid_mosaic !! @author George Gayno subroutine define_input_grid_grib2(localpet,npets) - use netcdf - use mpi use grib_mod use gdswzd_mod use program_setup, only : grib2_file_input_grid, data_dir_input_grid implicit none - integer, intent(in) :: localpet, npets + integer, intent(in) :: localpet, npets - character(len=500) :: the_file + character(len=500) :: the_file + + integer :: i, j, k, jdisc, jgdtn, jpdtn, lugb, lugi + integer :: jids(200), jgdt(200), jpdt(200), rc + integer :: kgds(200), nret, clb(2), cub(2) - integer :: i, j, k, jdisc, jgdtn, jpdtn, lugb, lugi - integer :: jids(200), jgdt(200), jpdt(200), rc - integer :: kgds(200), nret, clb(2), cub(2) - integer :: error, ncid, dim_i, dim_ip1, dim_j, dim_jp1 - integer :: fsize=65536 - integer :: header_buffer_val = 16384 - integer :: id_gridlat, id_gridlon, id_gridrot - integer :: id_gridlat_corners, id_gridlon_corners - integer :: id_gridlat_diff, id_gridlon_diff - - logical :: unpack - - real :: res - real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:) - real, allocatable :: rlon_corner(:,:),rlat_corner(:,:) - real, allocatable :: rlon_diff(:,:),rlat_diff(:,:) - real, allocatable :: xpts_corner(:,:),ypts_corner(:,:) + logical :: unpack + + real :: res + real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:) + real, allocatable :: rlon_corner(:,:),rlat_corner(:,:) + real, allocatable :: rlon_diff(:,:),rlat_diff(:,:) + real, allocatable :: xpts_corner(:,:),ypts_corner(:,:) real(esmf_kind_r8), allocatable :: latitude(:,:) real(esmf_kind_r8), allocatable :: longitude(:,:) real(esmf_kind_r8), allocatable :: latitude_corner(:,:) @@ -728,14 +720,6 @@ subroutine define_input_grid_grib2(localpet,npets) deallocate(xpts, ypts) - if (localpet == 0) then - print*,'after gdswzd lat/lon 11', rlat(1,1),rlon(1,1) - print*,'after gdswzd lat/lon ni/11', rlat(i_input,1),rlon(i_input,1) - print*,'after gdswzd lat/lon 11/nj', rlat(1,j_input),rlon(1,j_input) - print*,'after gdswzd lat/lon ni/nj', rlat(i_input,j_input),rlon(i_input,j_input) - print*,'after gdswzd lat/lon mid ', rlat(i_input/2,j_input/2),rlon(i_input/2,j_input/2) - endif - do j = 1, jp1_input do i = 1, ip1_input xpts_corner(i,j) = float(i) - 0.5 @@ -752,89 +736,6 @@ subroutine define_input_grid_grib2(localpet,npets) deallocate(xpts_corner, ypts_corner) - if (localpet == 0) then - print*,'after gdswzd lat/lon corner 11', rlat_corner(1,1),rlon_corner(1,1) - print*,'after gdswzd lat/lon corner ni/11', rlat_corner(ip1_input,1),rlon_corner(ip1_input,1) - print*,'after gdswzd lat/lon corner 11/nj', rlat_corner(1,jp1_input),rlon_corner(1,jp1_input) - print*,'after gdswzd lat/lon corner ni/nj', rlat_corner(ip1_input,jp1_input),rlon_corner(ip1_input,jp1_input) - print*,'after gdswzd lat/lon corner mid ', rlat_corner(i_input/2,j_input/2),rlon_corner(i_input/2,j_input/2) - print*,'after gdswzd max/min corner ',maxval(rlat_corner),minval(rlat_corner),maxval(rlon_corner),minval(rlon_corner) - - if (gfld%igdtnum == 32769) then - where(rlon > 180.0) rlon = rlon - 360.0 - where(rlon_corner > 180.0) rlon_corner = rlon_corner - 360.0 - endif - - do j = 1, j_input - do i = 1, i_input - rlat_diff(i,j) = rlat_corner(i,j) - rlat(i,j) - rlon_diff(i,j) = rlon_corner(i,j) - rlon(i,j) - enddo - enddo - - ncid = 34 - error = nf90_create("./latlon.nc", IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), & - ncid, initialsize=0, chunksize=fsize) - call netcdf_err(error, 'CREATING FILE=latlon.nc' ) - - error = nf90_def_dim(ncid, 'ny', j_input, dim_j) - call netcdf_err(error, 'DEFINING j DIMENSION' ) - - error = nf90_def_dim(ncid, 'nx', i_input, dim_i) - call netcdf_err(error, 'DEFINING i DIMENSION' ) - - error = nf90_def_dim(ncid, 'ny_stag', jp1_input, dim_jp1) - call netcdf_err(error, 'DEFINING jp1 DIMENSION' ) - - error = nf90_def_dim(ncid, 'nx_stag', ip1_input, dim_ip1) - call netcdf_err(error, 'DEFINING ip1 DIMENSION' ) - - error = nf90_def_var(ncid, 'gridlat', NF90_FLOAT, (/dim_i,dim_j/), id_gridlat) - call netcdf_err(error, 'DEFINING gridlat' ) - - error = nf90_def_var(ncid, 'gridlon', NF90_FLOAT, (/dim_i,dim_j/), id_gridlon) - call netcdf_err(error, 'DEFINING gridlon' ) - - error = nf90_def_var(ncid, 'gridrot', NF90_FLOAT, (/dim_i,dim_j/), id_gridrot) - call netcdf_err(error, 'DEFINING gridrot' ) - - error = nf90_def_var(ncid, 'gridlat_corners', NF90_FLOAT, (/dim_ip1,dim_jp1/), id_gridlat_corners) - call netcdf_err(error, 'DEFINING gridlat_corners' ) - - error = nf90_def_var(ncid, 'gridlon_corners', NF90_FLOAT, (/dim_ip1,dim_jp1/), id_gridlon_corners) - call netcdf_err(error, 'DEFINING gridlon_corners' ) - - error = nf90_def_var(ncid, 'gridlat_corner_minus_cent', NF90_FLOAT, (/dim_i,dim_j/), id_gridlat_diff) - call netcdf_err(error, 'DEFINING gridlat_corner_minus_cent' ) - - error = nf90_def_var(ncid, 'gridlon_corner_minus_cent', NF90_FLOAT, (/dim_i,dim_j/), id_gridlon_diff) - call netcdf_err(error, 'DEFINING gridlon_corner_minus_cent' ) - - error = nf90_enddef(ncid, header_buffer_val,4,0,4) - call netcdf_err(error, 'DEFINING HEADER' ) - - error = nf90_put_var(ncid, id_gridlat, real(rlat,4) ) - call netcdf_err(error, 'writing gridlat' ) - - error = nf90_put_var(ncid, id_gridlon, real(rlon,4) ) - call netcdf_err(error, 'writing gridlon' ) - - error = nf90_put_var(ncid, id_gridlat_diff, real(rlat_diff,4) ) - call netcdf_err(error, 'writing gridlat_diff' ) - - error = nf90_put_var(ncid, id_gridlon_diff, real(rlon_diff,4) ) - call netcdf_err(error, 'writing gridlon_diff' ) - - error = nf90_put_var(ncid, id_gridlat_corners, real(rlat_corner,4) ) - call netcdf_err(error, 'writing gridlat_corner' ) - - error = nf90_put_var(ncid, id_gridlon_corners, real(rlon_corner,4) ) - call netcdf_err(error, 'writing gridlon_corner' ) - - error = nf90_close(ncid) - - endif - if (gfld%igdtnum == 0) then ! gfs lat/lon data print*,"- CALL GridCreate1PeriDim FOR INPUT GRID." diff --git a/sorc/chgres_cube.fd/static_data.F90 b/sorc/chgres_cube.fd/static_data.F90 index ad9634da0..081d6beac 100644 --- a/sorc/chgres_cube.fd/static_data.F90 +++ b/sorc/chgres_cube.fd/static_data.F90 @@ -46,8 +46,7 @@ module static_data !! @author George Gayno NCEP/EMC subroutine get_static_fields(localpet) - use model_grid, only : target_grid, & - num_tiles_target_grid, & + use model_grid, only : num_tiles_target_grid, & i_target, j_target implicit none diff --git a/sorc/chgres_cube.fd/surface.F90 b/sorc/chgres_cube.fd/surface.F90 index 7df645cfa..8213d17ba 100644 --- a/sorc/chgres_cube.fd/surface.F90 +++ b/sorc/chgres_cube.fd/surface.F90 @@ -323,7 +323,7 @@ subroutine interp(localpet) integer :: clb_target(2), cub_target(2) integer :: isrctermprocessing integer :: num_fields - integer :: sotyp_ind, vgfrc_ind, mmvg_ind, lai_ind + integer :: vgfrc_ind, mmvg_ind, lai_ind integer, allocatable :: search_nums(:) integer(esmf_kind_i4), pointer :: unmapped_ptr(:) integer(esmf_kind_i4), pointer :: mask_input_ptr(:,:)