Ensure that pixels crossing the dateline are simulated correctly. (c175ecde) · Commits · CAMS / CSO

oper/src/cso_grid.F90

+344 −91

Original line number	Diff line number	Diff line
		@@ -9,6 +9,12 @@
		! Fixed wrong gather of lon/lat arrays.
		! Support input and output of packed variables.
		!
		! 2023-01, Arjo Segers
		! Added "GetWeightsCell" methodes to assign points to grid cells.
		!
		! 2024-09, Arjo Segers
		! Ensure that pixels crossing the dateline are simulated correctly.
		!
		!
		!###############################################################################
		!
		@@ -44,9 +50,13 @@ module CSO_Grid

		type :: T_GridMapping
		! grid description:
		real :: west, south
		!~ global bounding box
		real :: west, east, south, north
		!~ spacing:
		real :: dlon, dlat
		!~ local domain offset in global index space:
		integer :: ilon0, ilat0
		!~ local shape:
		integer :: nlon, nlat
		! mapping type:
		integer :: levels
		@@ -74,6 +84,10 @@ module CSO_Grid
		procedure :: GridMapping_GetWeights_1d
		generic :: GetWeights => GridMapping_GetWeights_0d, &
		GridMapping_GetWeights_1d
		procedure :: GridMapping_GetWeightsCell_0d
		procedure :: GridMapping_GetWeightsCell_1d
		generic :: GetWeightsCell => GridMapping_GetWeightsCell_0d, &
		GridMapping_GetWeightsCell_1d
		end type T_GridMapping


		@@ -95,8 +109,8 @@ contains
		! - cell offset in global domain
		!

		subroutine GridMapping_Init( self, west , dlon, ilon0, nlon, &
		south, dlat, ilat0, nlat, &
		subroutine GridMapping_Init( self, west , east , dlon, ilon0, nlon, &
		south, north, dlat, ilat0, nlat, &
		levels, status )

		use CSO_PArray, only : CSO_PArray_Init
		@@ -104,7 +118,7 @@ contains
		! --- in/out ---------------------------------

		class(T_GridMapping), intent(out) :: self
		real, intent(in) :: west, south
		real, intent(in) :: west, east, south, north
		real, intent(in) :: dlon, dlat
		integer, intent(in) :: ilon0, ilat0
		integer, intent(in) :: nlon, nlat
		@@ -121,10 +135,12 @@ contains

		! store:
		self%west = west
		self%east = east
		self%dlon = dlon
		self%ilon0 = ilon0
		self%nlon = nlon
		self%south = south
		self%north = north
		self%dlat = dlat
		self%ilat0 = ilat0
		self%nlat = nlat
		@@ -281,7 +297,47 @@ contains
		! total area of polygon:
		area = sum(self%wwp)

		! index box:
		! Assuming that pixels were alredy filtered on overlap,
		! so (xxp(:),yyp(:)) should all be inside domain.
		! However, if pixels were overlapping the dateline,
		! the longitudes might be outside [-180,180] and therefore outside the domain;
		! reset their values:
		do
		! points left from dateline?
		if ( any( self%xxp < self%west ) ) then
		! reset:
		where ( self%xxp < self%west )
		self%xxp = self%xxp + 360.0
		endwhere
		else
		! leave:
		exit
		end if
		end do
		do
		! points right from dateline?
		if ( any( self%xxp > self%east ) ) then
		! reset:
		where ( self%xxp > self%east )
		self%xxp = self%xxp - 360.0
		endwhere
		else
		! leave:
		exit
		end if
		end do
		! safety check:
		do ip = 1, np
		if ( (self%xxp(ip) < self%west ) .or. (self%xxp(ip) > self%east ) .or. &
		(self%yyp(ip) < self%south) .or. (self%yyp(ip) > self%north) ) then
		write (csol,'("found pixel centroid (",f12.5,",",f12.5,")")') self%xxp(ip), self%yyp(ip); call csoErr
		write (csol,'("outside global domain [",f12.5,",",f12.5,"] x [",f12.5,",",f12.5,"]")') &
		self%west, self%east, self%south, self%north; call csoErr
		TRACEBACK; status=1; return
		end if
		end do

		! subset of global space, used to speedup loops:
		i1 = max( 1, ceiling( (minval(self%xxp) - self%west )/self%dlon ) - self%ilon0 )
		i2 = min( self%nlon, ceiling( (maxval(self%xxp) - self%west )/self%dlon ) - self%ilon0 )
		j1 = max( 1, ceiling( (minval(self%yyp) - self%south)/self%dlat ) - self%ilat0 )
		@@ -465,6 +521,203 @@ contains

		end subroutine GridMapping_GetWeights_1d

		!
		! ***
		!

		!
		! Assign points (x,y) to grid cell (i,j).
		! The interpolation weight w is set to 1.0.
		! The area is for footprints used to normalize,
		! here set to 1.0 to avoid division by zero.
		!

		subroutine GridMapping_GetWeightsCell_0d( self, x, y, area, i, j, w, status )

		! --- in/out ---------------------------------

		class(T_GridMapping), intent(inout) :: self
		real, intent(in) :: x, y ! [degree]
		real, intent(out) :: area
		integer, intent(out) :: i, j
		real , intent(out) :: w
		integer, intent(out) :: status

		! --- const ----------------------------------

		character(len=*), parameter :: rname = mname//'/GridMapping_GetWeightsCell_0d'

		! --- local ----------------------------------

		! --- begin ----------------------------------

		! target cell, first cell if exactly on first edge:
		if ( x == self%west ) then
		i = 1 - self%ilon0
		else
		i = ceiling( (x-self%west )/self%dlon ) - self%ilon0
		end if
		if ( y == self%south ) then
		j = 1 - self%ilat0
		else
		j = ceiling( (y-self%south)/self%dlat ) - self%ilat0
		end if
		! in range?
		if ( (i >= 1) .and. (i <= self%nlon) .and. (j >= 1) .and. (j <= self%nlat) ) then
		! full weight:
		w = 1.0
		else
		! outside map, no weight:
		w = 0.0
		end if

		! dummy "pixel" area; used for normization with weights, so set to unity:
		area = 1.0

		! ok
		status = 0

		end subroutine GridMapping_GetWeightsCell_0d

		! *

		!
		! Array of pixel centers, select corresponding cell.
		! Input:
		! x, y : pixel centers
		! Ouptut:
		! area : pixel area (here 1.0, corners are unknown)
		! iw0, nw : per pixel the offset and number of elements in ii/jj/ww arrays ;
		! here nw==1 for all pixels
		! ii, jj : source cell indices
		! ww : source cell weights (here always 1.0)
		!

		subroutine GridMapping_GetWeightsCell_1d( self, x, y, &
		area, iw0, nw, ii, jj, ww, status )

		use CSO_PArray, only : CSO_PArray_Reshape

		! --- in/out ---------------------------------

		class(T_GridMapping), intent(inout) :: self
		real, intent(in) :: x(:), y(:) ! (npix) [degree]
		real, pointer :: area(:) ! (npix) [m2]
		integer, pointer :: iw0(:) ! (npix)
		integer, pointer :: nw(:) ! (npix)
		integer, pointer :: ii(:), jj(:) ! (nw)
		real, pointer :: ww(:) ! (nw) [m2]
		integer, intent(out) :: status

		! --- const ----------------------------------

		character(len=*), parameter :: rname = mname//'/GridMapping_GetWeightsCell_1d'

		! --- local ----------------------------------

		integer :: npix
		integer :: ipix
		real :: pix_area
		integer :: pix_i, pix_j
		real :: pix_w
		integer :: pix_nw
		integer :: nnew

		! --- begin ----------------------------------

		! number of pixels:
		npix = size(x)
		! check ...
		if ( size(y) /= npix ) then
		write (csol,'("arrays x (",i0,") and y (",i0,") should have same shape")') &
		size(x), size(y); call csoErr
		TRACEBACK; status=1; return
		end if

		! storage for pixel area; check current storage:
		call CSO_PArray_Reshape( self%all_area, npix, status )
		IF_NOT_OK_RETURN(status=1)
		call CSO_PArray_Reshape( self%all_iw0 , npix, status )
		IF_NOT_OK_RETURN(status=1)
		call CSO_PArray_Reshape( self%all_nw , npix, status )
		IF_NOT_OK_RETURN(status=1)

		! reset counter:
		self%nall = 0
		self%all_nw = 0

		! loop over pixels:
		do ipix = 1, npix

		! pixel mapping weights:
		call self%GetWeightsCell( x(ipix), y(ipix), &
		pix_area, pix_i, pix_j, pix_w, status )
		if ( status /= 0 ) then
		write (csol,'("could not compute mapping weights for pixel ",i0)') ipix; call csoErr
		write (csol,*) ' xx = ', x(ipix); call csoErr
		write (csol,*) ' yy = ', y(ipix); call csoErr
		TRACEBACK; status=ipix; return
		end if
		! single source cell:
		pix_nw = 1

		! store pixel area:
		self%all_area(ipix) = pix_area

		! offset and number of overlapping cells (might be 0 ...):
		self%all_iw0 (ipix) = self%nall
		self%all_nw (ipix) = pix_nw

		! any overlap? some pixels might not overlap with domain ...
		if ( pix_nw > 0 ) then
		! exceeds maximum storage?
		if ( self%nall + pix_nw > self%mxall ) then
		! new size, extend with 1 value extra per cell until it fits ...
		do
		self%mxall = self%mxall + self%nlon*self%nlat
		if ( self%nall + pix_nw <= self%mxall ) exit
		end do
		! extend arrays, copy current:
		call CSO_PArray_Reshape( self%all_ii , self%mxall, status )
		IF_NOT_OK_RETURN(status=1)
		call CSO_PArray_Reshape( self%all_jj , self%mxall, status )
		IF_NOT_OK_RETURN(status=1)
		call CSO_PArray_Reshape( self%all_ww , self%mxall, status )
		IF_NOT_OK_RETURN(status=1)
		end if
		! store pixel mapping:
		self%all_ii (self%nall+1:self%nall+pix_nw) = pix_i
		self%all_jj (self%nall+1:self%nall+pix_nw) = pix_j
		self%all_ww (self%nall+1:self%nall+pix_nw) = pix_w
		! increase counter:
		self%nall = self%nall + pix_nw
		end if ! nw > 0

		end do ! ipix

		! truncate to length that is actually used:
		call CSO_PArray_Reshape( self%all_ii , self%nall, status )
		IF_NOT_OK_RETURN(status=1)
		call CSO_PArray_Reshape( self%all_jj , self%nall, status )
		IF_NOT_OK_RETURN(status=1)
		call CSO_PArray_Reshape( self%all_ww , self%nall, status )
		IF_NOT_OK_RETURN(status=1)
		! reset maximum size:
		self%mxall = self%nall

		! set pointers:
		area => self%all_area
		iw0 => self%all_iw0
		nw => self%all_nw
		ii => self%all_ii
		jj => self%all_jj
		ww => self%all_ww

		! ok
		status = 0

		end subroutine GridMapping_GetWeightsCell_1d


		! ====================================================================
		! ===

oper/src/cso_tools.F90

+503 −84

Original line number	Diff line number	Diff line
		@@ -2,6 +2,17 @@
		!
		! CSO_Tools
		!
		!
		! HISTORY
		!
		! 2022-11, Arjo Segers
		! Added 'LinInterp' routine to perform linear interpolation.
		!
		! 2024-01, Arjo Segers
		! Added 'PointInsideDomain' to check if pixel center is inside domain.
		! Added 'PolygonOverlapsDomain' to check if pixel footprint overlaps with domain.
		!
		!
		!###############################################################################
		!
		#define TRACEBACK write (csol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call csoErr
		@@ -42,6 +53,9 @@ module CSO_Tools
		public :: LonLatTriangleArea
		public :: LonLatRectangleArea
		public :: GetPolygonPoints
		public :: LinInterp
		public :: PointInsideDomain
		public :: PolygonOverlapsDomain


		! --- const ------------------------------
		@@ -56,6 +70,12 @@ module CSO_Tools
		! conversion from degrees to radians:
		real, parameter :: deg2rad = pi/180.0 ! rad/deg

		! bounding box indices:
		integer, parameter :: ibb_west = 1
		integer, parameter :: ibb_east = 2
		integer, parameter :: ibb_south = 3
		integer, parameter :: ibb_north = 4


		contains

		@@ -572,6 +592,405 @@ contains
		end subroutine GetPolygonPoints


		! ***

		!
		! Linear interpolation values y(1:n) defined on x(1:n) to target location xp.
		! Use 'extrapolate=.true.' to allow extrapolation outside bounds of x.
		! Coordinate x should be strictly increasing or decreasing.
		!

		subroutine LinInterp( x, y, xp, yp, status, extrapolate )

		! --- in/out ---------------------------------

		real, intent(in) :: x(:) ! (n)
		real, intent(in) :: y(:) ! (n)
		real, intent(in) :: xp
		real, intent(out) :: yp
		integer, intent(out) :: status
		logical, intent(in), optional :: extrapolate

		! --- const ----------------------------------

		character(len=*), parameter :: rname = mname//'/LinInterp'

		! --- local ----------------------------------

		integer :: n
		integer :: xdir
		integer :: i
		integer :: ip
		logical :: extrapol

		! --- begin ----------------------------------

		! flag
		extrapol = .false.
		if ( present(extrapolate) ) extrapol = extrapolate

		! sizxe:
		n = size(x)

		! check ...
		if ( size(y) /= n ) then
		write (csol,'("size of y (",i0,") should match size of x (",i0,")")') size(y), n; call csoErr
		TRACEBACK; status=1; return
		end if

		! check direction:
		if ( x(n) > x(1) ) then
		! direction factor:
		xdir = 1
		else
		! direction factor:
		xdir = -1
		end if
		! coordinate should be strictly increasing if multiplied with 'xdir':
		if ( any( xdir*( x(2:n) - x(1:n-1) ) <= 0.0 ) ) then
		write (csol,'("coordinate should be strictly increasing or decreasing:")'); call csoErr
		do i = 1, n
		write (csol,'(" ",i6," ",es16.6)') i, x(i); call csoErr
		end do
		TRACEBACK; status=1; return
		end if

		! search interval 1,..,n-1 holding xp:
		!~ left from coordinates?
		if ( xdirxp < xdirx(1) ) then
		! allowed?
		if ( extrapol ) then
		! extrapolate from first cell:
		ip = 1
		else
		write (csol,'("xp ",es16.6," outside x range [",es16.6,",",es16.6,"]")') xp, x(1), x(n); call csoErr
		TRACEBACK; status=1; return
		end if
		!~ right from coordinates?
		else if ( xdirxp > xdirx(n) ) then
		! allowed?
		if ( extrapol ) then
		! extrapolate from last cell:
		ip = n-1
		else
		write (csol,'("xp ",es16.6," outside x range [",es16.6,",",es16.6,"]")') xp, x(1), x(n); call csoErr
		TRACEBACK; status=1; return
		end if
		!~ within coordinates:
		else
		! loop over end points of intervals:
		do i = 2, n
		! below end of interval? then this holds the target point:
		if ( xdirxp <= xdirx(i) ) then
		! interpolate from cell with this end-point:
		ip = i - 1
		exit
		end if
		end do ! i
		end if ! interpolate or extrapolate:

		! linear interpolation (or extrapolation):
		yp = ( ( x(ip+1) - xp ) * y(ip) + ( xp - x(ip) ) * y(ip+1) )/( x(ip+1) - x(ip) )

		!! testing ..
		!write (csol,'("interpolation on x of y to xp:")'); call csoPr
		!do i = 1, n
		! write (csol,'(" ",i6," ",2f16.6)') i, x(i), y(i); call csoPr
		! if ( i == ip ) then
		! write (csol,'(" ---> ",2f16.6)') xp, yp; call csoPr
		! end if
		!end do

		! ok
		status = 0

		end subroutine LinInterp


		! ***


		! *


		!
		! Return .true. if point (x,y) in degrees is inside domain [west,east,south,north]:
		!
		! +----------+
		! \| * \|
		! \| \|
		! +----------+
		!

		subroutine PointInsideDomain( x, y, domain, inside, status )

		! --- in/out ---------------------------------

		real, intent(in) :: x, y
		real, intent(in) :: domain(4) ! [west,east,south,north]
		logical, intent(out) :: inside
		integer, intent(out) :: status

		! --- const ----------------------------------

		character(len=*), parameter :: rname = mname//'/PointInsideDomain'

		! --- local ----------------------------------

		! --- begin ----------------------------------

		! point inside box?
		inside = ( x >= domain(ibb_west ) ) .and. ( x <= domain(ibb_east ) ) .and. &
		( y >= domain(ibb_south) ) .and. ( y <= domain(ibb_north) )

		! ok
		status = 0

		end subroutine PointInsideDomain


		!
		! Return bounding box "[west,east,south,north]" around polygon with corners:
		! [ (x(1),y(1)), .., (x(n),y(n)) ]
		! Coordinates are assumed to be in degrees, result is in [-180,180]
		! (thus west might be larger than east!)
		!

		subroutine BoundingBoxFromPolygon( xx, yy, bbox, status )

		! --- in/out ---------------------------------

		real, intent(in) :: xx(:) ! (n)
		real, intent(in) :: yy(:) ! (n)
		real, intent(out) :: bbox(4)
		integer, intent(out) :: status

		! --- const ----------------------------------

		character(len=*), parameter :: rname = mname//'/BoundingBoxFromPolygon'

		! --- local ----------------------------------

		! --- begin ----------------------------------

		! set to min/max range:
		bbox(ibb_west) = minval(xx)
		bbox(ibb_east) = maxval(xx)
		bbox(ibb_south) = minval(yy)
		bbox(ibb_north) = maxval(yy)

		! ok
		status = 0

		end subroutine BoundingBoxFromPolygon


		! *


		!
		! Return .true. if bounding box 1 overlaps with box 2.
		!
		! +--------+
		! \| +-------+
		! \| \|xxx\| \|
		! +--------+ \|
		! +-------+
		!

		subroutine BoundingBoxOverlaps( bbox1, bbox2, overlap, status )!, debug )

		! --- in/out ---------------------------------

		real, intent(in) :: bbox1(4)
		real, intent(in) :: bbox2(4)
		logical, intent(out) :: overlap
		integer, intent(out) :: status
		!logical, intent(in), optional :: debug

		! --- const ----------------------------------

		character(len=*), parameter :: rname = mname//'/BoundingBoxOverlaps'

		! --- local ----------------------------------

		!logical :: dbg
		logical :: overlap_we

		! --- begin ----------------------------------

		!! flag:
		!dbg = .false.
		!if ( present(debug) ) dbg = debug

		! safety check ..
		if ( (abs(bbox1(ibb_west)) > 180.0+360.0) .or. &
		(abs(bbox1(ibb_east)) > 180.0+360.0) ) then
		write (csol,'("bounding box with unsupported longitudes ",2f16.8)') &
		bbox1(ibb_west), bbox1(ibb_east); call csoErr
		TRACEBACK; status=1; return
		end if

		! overlap in west/east direction, also check shifts over 360 deg:
		overlap_we = ( ( bbox1(ibb_west ) <= bbox2(ibb_east ) ) .and. &
		( bbox1(ibb_east ) >= bbox2(ibb_west ) ) ) .or. &
		( ( bbox1(ibb_west ) - 360.0 <= bbox2(ibb_east ) ) .and. &
		( bbox1(ibb_east ) - 360.0 >= bbox2(ibb_west ) ) ) .or. &
		( ( bbox1(ibb_west ) + 360.0 <= bbox2(ibb_east ) ) .and. &
		( bbox1(ibb_east ) + 360.0 >= bbox2(ibb_west ) ) )

		! combine with overlap in south/north direction:
		overlap = overlap_we .and. &
		( bbox1(ibb_south) <= bbox2(ibb_north) ) .and. &
		( bbox1(ibb_north) >= bbox2(ibb_south) )

		! ok
		status = 0

		end subroutine BoundingBoxOverlaps


		! *


		!
		! Return .true. if bounding box 1 is inside box 2.
		!
		! 2 +----------+
		! \| 1 +---+ \|
		! \| \|xxx\| \|
		! \| +---+ \|
		! +----------+
		!

		subroutine BoundingBoxInside( bbox1, bbox2, inside, status )

		! --- in/out ---------------------------------

		real, intent(in) :: bbox1(4)
		real, intent(in) :: bbox2(4)
		logical, intent(out) :: inside
		integer, intent(out) :: status

		! --- const ----------------------------------

		character(len=*), parameter :: rname = mname//'/BoundingBoxInside'

		! --- local ----------------------------------

		! --- begin ----------------------------------

		! west or east outside [-180,+180] ?
		if ( bbox1(ibb_west ) < -180.0 ) then
		! check on 360-degree shift:
		inside = ( bbox1(ibb_west ) + 360.0 >= bbox2(ibb_west ) ) .and. &
		( bbox1(ibb_east ) <= bbox2(ibb_east ) ) .and. &
		( bbox1(ibb_south) >= bbox2(ibb_south) ) .and. &
		( bbox1(ibb_north) <= bbox2(ibb_north) )

		else if ( bbox1(ibb_east ) > 180.0 ) then
		! check on 360-degree shift:
		inside = ( bbox1(ibb_west ) >= bbox2(ibb_west ) ) .and. &
		( bbox1(ibb_east ) -360.0 <= bbox2(ibb_east ) ) .and. &
		( bbox1(ibb_south) >= bbox2(ibb_south) ) .and. &
		( bbox1(ibb_north) <= bbox2(ibb_north) )
		else
		! check using original values:
		inside = ( bbox1(ibb_west ) >= bbox2(ibb_west ) ) .and. &
		( bbox1(ibb_east ) <= bbox2(ibb_east ) ) .and. &
		( bbox1(ibb_south) >= bbox2(ibb_south) ) .and. &
		( bbox1(ibb_north) <= bbox2(ibb_north) )
		end if

		! ok
		status = 0

		end subroutine BoundingBoxInside


		! *


		!
		! Return .true. if polygon with corners in degrees:
		! [ (x(1),y(1)), .., (x(n),y(n)) ]
		! overlaps with domain bounding box:
		! [west,east,south,north].
		! If "inside_domain" is present it defines a boundig box
		! that should have the polygion inside, for example a global domain.
		!
		! +-------------------+ inside_domain
		! \| +---+ domain \|
		! \| \| \| \|
		! \| [xx] \| \|
		! \| +---+ \|
		! +-------------------+
		!

		subroutine PolygonOverlapsDomain( xx, yy, domain, overlaps, status, &
		inside_domain )!, debug )

		! --- in/out ---------------------------------

		real, intent(in) :: xx(:) ! (n) degrees
		real, intent(in) :: yy(:) ! (n) degrees
		real, intent(in) :: domain(4) ! [west,east,south,north]
		logical, intent(out) :: overlaps
		integer, intent(out) :: status
		real, intent(in), optional :: inside_domain(4) ! [west,east,south,north]
		!logical, intent(in), optional :: debug

		! --- const ----------------------------------

		character(len=*), parameter :: rname = mname//'/PolygonOverlapsDomain'

		! --- local ----------------------------------

		!logical :: dbg
		real :: bbox(4)
		logical :: inside

		! --- begin ----------------------------------

		!! flag:
		!dbg = .false.
		!if ( present(debug) ) dbg = debug

		! bounding box (west,east,south,north),
		! longitudes remain outside [-180,180] if already so on input:
		call BoundingBoxFromPolygon( xx, yy, bbox, status )
		IF_NOT_OK_RETURN(status=1)

		!! testing ..
		!if ( dbg ) then
		! write (csol,*) 'xxx1 xx range ', minval(xx), maxval(xx); call csoPr
		! write (csol,*) '..x1 bbox ', bbox; call csoPr
		!end if

		! check if box overlaps with local domain:
		call BoundingBoxOverlaps( bbox, domain, overlaps, status )!, debug=debug )
		IF_NOT_OK_RETURN(status=1)

		!! testing ..
		!if ( dbg ) then
		! write (csol,*) '..x1 overlaps ', overlaps; call csoPr
		!end if

		! check if inside specified (global) domain?
		! no need to check if does not overlap with local domain ..
		if ( present(inside_domain) .and. overlaps ) then
		! check if box is entirely in (global) domain:
		call BoundingBoxInside( bbox, inside_domain, inside, status )
		IF_NOT_OK_RETURN(status=1)
		! update:
		overlaps = overlaps .and. inside
		end if

		! ok
		status = 0

		end subroutine PolygonOverlapsDomain

		end module CSO_Tools

oper/src/tutorial_oper_S5p.F90

+71 −65

File changed.

Preview size limit exceeded, changes collapsed.