program camse implicit none real, parameter :: REARTH = 6370. !km real, parameter :: ZERO = 0. real, parameter :: FOUR = 4. real, parameter :: NINETY = 180. real, parameter :: ONE80 = 180. real, parameter :: THREE60 = 360. integer :: astat, ncid, ncStat integer :: nPolygons, maxVtx integer :: dimid, varid integer :: nMatchLon, nMatchLat, nUnique, nMatchVtx, nVtx integer :: minUniqueVtx integer :: polyNdx, vtxNdx integer :: minNdx(1), maxNdx(1) integer :: Quadcnt(4) integer, allocatable :: Quad(:) integer, allocatable :: cntVtx(:) real :: matchVtxlon, matchVtxlat real :: PI, D2R, R2D real :: EarthSA real, allocatable :: wrkLon(:) real, allocatable :: polyarea(:), polylon(:), polylat(:) real, allocatable :: polyVtxlon(:,:), polyVtxlat(:,:) character(len=64) :: mess character(len=64) :: filenm character(len=64) :: Fullfilenm logical :: longCross logical, allocatable :: vtxMask(:) include 'netcdf.inc' PI = 4.*atan( 1.) D2R = PI/ONE80 ; R2D = ONE80/PI EarthSA = FOUR*PI*REARTH*REARTH !----------------------------------------------------------------------- ! specify and open netcdf file !----------------------------------------------------------------------- filenm = 'conus_30_x8.g_scrip.nc' Fullfilenm = '../data/cam/' // trim( filenm) mess = 'Failed to open ' // trim(Fullfilenm) call handle_ncerr( nf_open( trim(Fullfilenm), nf_noclobber, ncid ), mess ) !----------------------------------------------------------------------- ! get number polygons !----------------------------------------------------------------------- mess = 'Failed to get grid_size dimension id' call handle_ncerr( nf_inq_dimid( ncid, 'grid_size', dimid ), mess ) mess = 'Failed to get grid_size dimension' call handle_ncerr( nf_inq_dimlen( ncid, dimid, nPolygons ), mess ) write(*,*) write(*,*) 'There are ', nPolygons, 'horizontal grid elements' write(*,*) !----------------------------------------------------------------------- ! get number vertices/polygon !----------------------------------------------------------------------- mess = 'Failed to get grid_corners dimension id' call handle_ncerr( nf_inq_dimid( ncid, 'grid_corners', dimid ), mess ) mess = 'Failed to get grid_corners dimension' call handle_ncerr( nf_inq_dimlen( ncid, dimid, maxVtx ), mess ) write(*,*) write(*,*) 'There are a maximum of ', maxVtx, ' vertices per polygon' write(*,*) !----------------------------------------------------------------------- ! allocate wrk variables !----------------------------------------------------------------------- allocate( polyarea(nPolygons), & polylon(nPolygons), polylat(nPolygons), cntVtx(nPolygons), & vtxMask(maxVtx), Quad(maxVtx), wrkLon(maxVtx), stat=astat) if( astat /= 0 ) then write(*,*) 'Failed to allocate area, lon, lat; error = ',astat stop 'AllocErr' endif allocate( polyVtxlon(maxVtx,nPolygons), polyVtxlat(maxVtx,nPolygons), stat=astat) if( astat /= 0 ) then write(*,*) 'Failed to allocate poly Vtx lon, lat; error = ',astat stop 'AllocErr' endif !----------------------------------------------------------------------- ! get poly center longitude (degrees) !----------------------------------------------------------------------- mess = 'Failed to get grid_center_lon variable id' call handle_ncerr( nf_inq_varid( ncid, 'grid_center_lon', varid ), mess ) mess = 'Failed to get grid_center_lon' call handle_ncerr( nf_get_var_real( ncid, varid, polylon ), mess ) write(*,*) write(*,*) 'Min,Max poly center lon = ',minval(polylon),maxval(polylon) write(*,*) !----------------------------------------------------------------------- ! get poly center latitude (degrees) !----------------------------------------------------------------------- mess = 'Failed to get grid_center_lat variable id' call handle_ncerr( nf_inq_varid( ncid, 'grid_center_lat', varid ), mess ) mess = 'Failed to get grid_center_lat' call handle_ncerr( nf_get_var_real( ncid, varid, polylat ), mess ) write(*,*) write(*,*) 'Min,Max poly center lat = ',minval(polylat),maxval(polylat) write(*,*) !----------------------------------------------------------------------- ! get poly area (radians) !----------------------------------------------------------------------- mess = 'Failed to get grid_area variable id' call handle_ncerr( nf_inq_varid( ncid, 'grid_area', varid ), mess ) mess = 'Failed to get grid_area' call handle_ncerr( nf_get_var_real( ncid, varid, polyarea ), mess ) !----------------------------------------------------------------------- ! get poly Vtx longitude (degrees) !----------------------------------------------------------------------- mess = 'Failed to get grid_corner_lon variable id' call handle_ncerr( nf_inq_varid( ncid, 'grid_corner_lon', varid ), mess ) mess = 'Failed to get grid_corner_lon' call handle_ncerr( nf_get_var_real( ncid, varid, polyVtxlon ), mess ) write(*,*) write(*,*) 'Min,Max poly vtx lon = ',minval(polyVtxlon),maxval(polyVtxlon) write(*,*) !----------------------------------------------------------------------- ! get poly Vtx latitude (degrees) !----------------------------------------------------------------------- mess = 'Failed to get grid_corner_lat variable id' call handle_ncerr( nf_inq_varid( ncid, 'grid_corner_lat', varid ), mess ) mess = 'Failed to get grid_corner_lat' call handle_ncerr( nf_get_var_real( ncid, varid, polyVtxlat ), mess ) write(*,*) write(*,*) 'Min,Max poly vtx lat = ',minval(polyVtxlat),maxval(polyVtxlat) write(*,*) !----------------------------------------------------------------------- ! area diagnostics !----------------------------------------------------------------------- minNdx(:) = minloc( polyarea(:) ) maxNdx(:) = maxloc( polyarea(:) ) write(*,*) write(*,*) 'Min poly_area centered at ',polylon(minNdx(1)),polylat(minNdx(1)) write(*,*) 'Max poly_area centered at ',polylon(maxNdx(1)),polylat(maxNdx(1)) write(*,*) 'Min,Max poly_area = ',EarthSA*polyarea(minNdx(1)),EarthSA*polyarea(maxNdx(1)) write(*,*) 'Min,Max poly_area DX = ',sqrt(EarthSA*polyarea(minNdx(1))), & sqrt(EarthSA*polyarea(maxNdx(1))) write(*,*) 'Sum poly_area = ',sum(polyarea),FOUR*PI write(*,*) !----------------------------------------------------------------------- ! check unique vtx count !----------------------------------------------------------------------- minUniqueVtx = 1000 do polyNdx = 1,nPolygons matchVtxlon = polyVtxlon(maxVtx,polyNdx) nMatchLon = count( polyVtxlon(1:maxVtx,polyNdx) == matchVtxlon ) matchVtxlat = polyVtxlat(maxVtx,polyNdx) nMatchLat = count( polyVtxlat(1:maxVtx,polyNdx) == matchVtxlat ) nMatchVtx = min( nMatchLon,nMatchLat ) if( nMatchLon /= nMatchLat ) then write(*,*) nMatchLon,' != ',nMatchLat,' for poly #',polyNdx vtxMask(:) = polyVtxlon(:,polyNdx) == matchVtxlon write(*,*) 'Longitude match mask: ',vtxMask(:) write(*,*) 'Vtx longitudes' write(*,*) polyVtxlon(:,polyNdx) vtxMask(:) = polyVtxlat(:,polyNdx) == matchVtxlat write(*,*) 'Latitude match mask: ',vtxMask(:) write(*,*) 'Vtx latitudes' write(*,*) polyVtxlat(:,polyNdx) endif nUnique = maxVtx - nMatchvtx + 1 cntVtx(polyNdx) = nUnique write(*,*) 'Unique vtx count = ',nUnique,' for poly #',polyNdx if( mod(nUnique,2) /= 0 ) then write(*,*) 'Unique vtx count is odd' endif minUniqueVtx = min( minUniqueVtx,nUnique ) where( polyVtxlon(:,polyNdx) > ONE80 ) wrkLon(:) = polyVtxlon(:,polyNdx) - THREE60 elsewhere wrkLon(:) = polyVtxlon(:,polyNdx) endwhere !----------------------------------------------------------------------- ! place vtx in "quadrant" by longitude !----------------------------------------------------------------------- Quadcnt(:) = 0 do vtxNdx = 1,nUnique if( ZERO <= wrkLon(vtxNdx) .and. wrkLon(vtxNdx) < NINETY ) then Quad(vtxNdx) = 1 Quadcnt(1) = Quadcnt(1) + 1 elseif( NINETY <= wrkLon(vtxNdx) .and. wrkLon(vtxNdx) < ONE80 ) then Quad(vtxNdx) = 2 Quadcnt(2) = Quadcnt(2) + 1 elseif( -ONE80 <= wrkLon(vtxNdx) .and. wrkLon(vtxNdx) < -NINETY ) then Quad(vtxNdx) = 3 Quadcnt(3) = Quadcnt(3) + 1 elseif( -NINETY <= wrkLon(vtxNdx) .and. wrkLon(vtxNdx) < ZERO ) then Quad(vtxNdx) = 4 Quadcnt(4) = Quadcnt(4) + 1 endif enddo if( minval(wrkLon(:)) * maxval(wrkLon(:)) < ZERO ) then if( count( Quadcnt(:) /= 0 ) > 2 ) then write(*,*) 'Polygon #',polyNdx,' vertices cross both 0 and 180' elseif( any(Quad(:nUnique) == 1) .and. any(Quad(:nUnique) == 4) ) then write(*,*) 'Polygon #',polyNdx,' vertices cross 0' elseif( any(Quad(:nUnique) == 2) .and. any(Quad(:nUnique) == 3) ) then write(*,*) 'Polygon #',polyNdx,' vertices cross 180' endif write(*,'(''Vtx cross longitudes: '',10f10.4)') wrkLon(1:nUnique) write(*,'(''Bounding polygon lats: '',2f10.4)') & minval(polyVtxlat(:,polyNdx)), maxval(polyVtxlat(:,polyNdx)) endif enddo write(*,*) ' ' write(*,*) 'Min Unique vtx count = ',minval(cntVtx(:)) write(*,*) 'Max Unique vtx count = ',maxval(cntVtx(:)) polyNdx = minNdx(1) nVtx = cntVtx(polyNdx) write(*,*) ' ' write(*,*) 'Debug check; nVtx,polyNdx = ',nVtx,polyNdx write(*,*) 'Min avg sub poly area = ',2.*EarthSA*polyarea(polyNdx)/real(nVtx) deallocate( polyarea, polylon, polylat, polyVtxlon, polyVtxlat ) deallocate( vtxMask, wrkLon, Quad ) ncStat = nf_close( ncid ) CONTAINS subroutine handle_ncerr( ret, mes ) !--------------------------------------------------------------------- ! ... netcdf error handling routine !--------------------------------------------------------------------- !--------------------------------------------------------------------- ! ... dummy arguments !--------------------------------------------------------------------- integer, intent(in) :: ret character(len=*), intent(in) :: mes if( ret /= nf_noerr ) then write(*,*) nf_strerror( ret ) stop 'netcdf error' endif end subroutine handle_ncerr end program camse