Updated Fortran operator. (a6db1fef) · Commits · CAMS / CSO

CHANGELOG

+6 −0

Original line number	Diff line number	Diff line
		@@ -549,3 +549,9 @@ Plot pixels als polygons in case track information is not present.
		src/cso/cso_plot.py


		vx.y
		----

		Updated Fortran operator.
		oper/

oper/CHANGELOG

+21 −0

Original line number	Diff line number	Diff line
		@@ -70,3 +70,24 @@ Support labels 'hours', 'minutes', 'seconds', etc.
		src/cso_datetime.F90


		vx.y
		----

		Changes to avoid compiler warnings.
		src/cso_comm.F90
		src/cso_ncfile.F90

		Updated comment.
		src/cso_mapping.F90

		Support 'fill_value' argument when reading 1D real array.
		src/cso_ncfile.F90
		src/cso_pixels.F90

		Support kernel convolution including prior state.
		src/cso_pixels.F90

		Added 'LinInterpWeights' routine to calculate interpolation weights.
		src/cso_tools.F90

		!

oper/src/cso_ncfile.F90

+133 −61

Original line number	Diff line number	Diff line
		@@ -25,6 +25,14 @@
		! 2023-11, Arjo Segers
		! Close files also in 'read' mode ...
		!
		! 2025-09, Arjo Segers
		! Support 'fill_value' argument when reading 1D real array.
		!
		! 2025-09, Arjo Segers
		! Renamed 'start' and 'count' arguments to 'vstart' and 'vcount'
		! to avoid that compiler gets confused with 'count()' function ...
		!
		!
		!###############################################################################
		!
		#define TRACEBACK write (csol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call csoErr
		@@ -266,6 +274,7 @@ module CSO_NcFile
		procedure :: Inq_Variable => NcFile_Inq_Variable
		procedure :: Inq_VarPacking => NcFile_Inq_VarPacking
		procedure :: Inq_VarMissing => NcFile_Inq_VarMissing
		procedure :: Inq_VarFillValue => NcFile_Inq_VarFillValue
		procedure :: Inq_VarUnits => NcFile_Inq_VarUnits
		!
		procedure :: NcFile_Get_Var_i_1d
		@@ -2935,8 +2944,8 @@ contains

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

		! try to get first packing variable:
		status = NF90_Get_Att( self%ncid, varid, 'missing_value', missing_value )
		! try to get fill value:
		status = NF90_Get_Att( self%ncid, varid, '_FillValue', missing_value )
		! no error, missing value available
		if ( status == NF90_NOERR ) then
		! missing value available,
		@@ -2964,6 +2973,57 @@ contains
		! *


		! Return missing value, status -1 if not found.

		subroutine NcFile_Inq_VarFillValue( self, varid, fill_value, status )

		use NetCDF, only : NF90_Get_Att
		use NetCDF, only : NF90_ENOTATT

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

		class(T_NcFile), intent(in) :: self
		integer, intent(in) :: varid
		real, intent(out) :: fill_value
		integer, intent(out) :: status

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

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

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

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

		! try to get fill value:
		status = NF90_Get_Att( self%ncid, varid, '_FillValue', fill_value )
		! no error, missing value available
		if ( status == NF90_NOERR ) then
		! missing value available,
		! remain in output argument

		! attribute not found ?
		else if ( status == NF90_ENOTATT ) then
		! no missing value, set dummy values:
		fill_value = -huge(1.0)
		! warning status:
		status = -1; return
		!
		else
		! some error ...
		csol=NF90_StrError(status); call csoErr
		TRACEBACK; status=1; return
		end if

		! ok
		status = 0

		end subroutine NcFile_Inq_VarFillValue


		! *


		!
		! Return units from attribute of variabled opened with 'varid'.
		! If attribute is not present, try to extract from 'description'
		@@ -3042,7 +3102,7 @@ contains


		subroutine NcFile_Get_Var_i1_1d( self, description, values, units, status, &
		start, count, missing_value )
		vstart, vcount, missing_value )

		use NetCDF, only : NF90_Get_Var
		use NetCDF, only : NF90_Get_Att
		@@ -3055,7 +3115,7 @@ contains
		integer(1), intent(out) :: values(:)
		character(len=*), intent(out) :: units
		integer, intent(out) :: status
		integer, intent(in), optional :: start(:), count(:)
		integer, intent(in), optional :: vstart(:), vcount(:)
		real, intent(out), optional :: missing_value

		! --- const --------------------------------------
		@@ -3074,7 +3134,7 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! read:
		status = NF90_Get_Var( self%ncid, varid, values, start=start, count=count )
		status = NF90_Get_Var( self%ncid, varid, values, start=vstart, count=vcount )
		IF_NF90_NOT_OK_RETURN(status=1)

		! packed?
		@@ -3105,7 +3165,7 @@ contains


		subroutine NcFile_Get_Var_i_1d( self, description, values, units, status, &
		start, count, missing_value )
		vstart, vcount, missing_value )

		use NetCDF, only : NF90_Get_Var
		use NetCDF, only : NF90_Get_Att
		@@ -3118,7 +3178,7 @@ contains
		integer, intent(out) :: values(:)
		character(len=*), intent(out) :: units
		integer, intent(out) :: status
		integer, intent(in), optional :: start(:), count(:)
		integer, intent(in), optional :: vstart(:), vcount(:)
		real, intent(out), optional :: missing_value

		! --- const --------------------------------------
		@@ -3137,7 +3197,7 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! read:
		status = NF90_Get_Var( self%ncid, varid, values, start=start, count=count )
		status = NF90_Get_Var( self%ncid, varid, values, start=vstart, count=vcount )
		IF_NF90_NOT_OK_RETURN(status=1)

		! packed?
		@@ -3166,7 +3226,7 @@ contains
		! *

		subroutine NcFile_Get_Var_c_2d( self, description, values, units, status, &
		start, count )
		vstart, vcount )

		use NetCDF, only : NF90_Get_Var
		use NetCDF, only : NF90_Get_Att
		@@ -3179,7 +3239,7 @@ contains
		character(len=1), intent(out) :: values(:,:)
		character(len=*), intent(out) :: units
		integer, intent(out) :: status
		integer, intent(in), optional :: start(:), count(:)
		integer, intent(in), optional :: vstart(:), vcount(:)

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

		@@ -3203,7 +3263,7 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! safety ..
		if ( present(start) .or. present(count) ) then
		if ( present(vstart) .or. present(vcount) ) then
		write (csol,'("optional arguments `start` or `count` not supported yet for char arrays")'); call csoErr
		TRACEBACK; status=1; return
		end if
		@@ -3235,7 +3295,7 @@ contains
		! *

		subroutine NcFile_Get_Var_i_2d( self, description, values, units, status, &
		start, count, missing_value )
		vstart, vcount, missing_value )

		use NetCDF, only : NF90_Get_Var
		use NetCDF, only : NF90_Get_Att
		@@ -3248,7 +3308,7 @@ contains
		integer, intent(out) :: values(:,:)
		character(len=*), intent(out) :: units
		integer, intent(out) :: status
		integer, intent(in), optional :: start(:), count(:)
		integer, intent(in), optional :: vstart(:), vcount(:)
		real, intent(out), optional :: missing_value

		! --- const --------------------------------------
		@@ -3267,7 +3327,7 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! read:
		status = NF90_Get_Var( self%ncid, varid, values, start=start, count=count )
		status = NF90_Get_Var( self%ncid, varid, values, start=vstart, count=vcount )
		IF_NF90_NOT_OK_RETURN(status=1)

		! packed?
		@@ -3298,7 +3358,7 @@ contains


		subroutine NcFile_Get_Var_i_3d( self, description, values, units, status, &
		start, count, missing_value )
		vstart, vcount, missing_value )

		use NetCDF, only : NF90_Get_Var
		use NetCDF, only : NF90_Get_Att
		@@ -3311,7 +3371,7 @@ contains
		integer, intent(out) :: values(:,:,:)
		character(len=*), intent(out) :: units
		integer, intent(out) :: status
		integer, intent(in), optional :: start(:), count(:)
		integer, intent(in), optional :: vstart(:), vcount(:)
		real, intent(out), optional :: missing_value

		! --- const --------------------------------------
		@@ -3330,7 +3390,7 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! read:
		status = NF90_Get_Var( self%ncid, varid, values, start=start, count=count )
		status = NF90_Get_Var( self%ncid, varid, values, start=vstart, count=vcount )
		IF_NF90_NOT_OK_RETURN(status=1)

		! packed?
		@@ -3361,7 +3421,7 @@ contains


		subroutine NcFile_Get_Var_r_1d( self, description, values, units, status, &
		start, count, missing_value )
		vstart, vcount, fill_value )

		use NetCDF, only : NF90_Get_Var
		use NetCDF, only : NF90_Get_Att
		@@ -3374,8 +3434,8 @@ contains
		real, intent(out) :: values(:)
		character(len=*), intent(out) :: units
		integer, intent(out) :: status
		integer, intent(in), optional :: start(:), count(:)
		real, intent(out), optional :: missing_value
		integer, intent(in), optional :: vstart(:), vcount(:)
		real, intent(in), optional :: fill_value

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

		@@ -3385,6 +3445,9 @@ contains

		integer :: varid
		real :: add_offset, scale_factor
		logical :: packed
		real :: fill_value_in
		logical :: filled

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

		@@ -3393,21 +3456,30 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! read:
		status = NF90_Get_Var( self%ncid, varid, values, start=start, count=count )
		status = NF90_Get_Var( self%ncid, varid, values, start=vstart, count=vcount )
		IF_NF90_NOT_OK_RETURN(status=1)

		! packed?
		! obtain fill value of input;
		! status -1 if not defined, fill_value_in=-huge(1.0) then:
		call self%Inq_VarFillValue( varid, fill_value_in, status )
		IF_ERROR_RETURN(status=1)

		! packed? status -1 if no packing attributes defined:
		call self%Inq_VarPacking( varid, add_offset, scale_factor, status )
		IF_ERROR_RETURN(status=1)
		if ( status == 0 ) then
		! unpack:
		where ( values /= fill_value_in )
		values = add_offset + scale_factor * values
		end where
		end if

		! Missing value?
		if ( present( missing_value ) ) then
		call self%Inq_VarMissing( varid, missing_value, status )
		IF_ERROR_RETURN(status=1)
		! reset to fill value?
		if ( present(fill_value) ) then
		! reset input values that have their fill_value:
		where ( values == fill_value_in )
		values = fill_value
		endwhere
		end if

		! get units:
		@@ -3422,7 +3494,7 @@ contains
		! *

		subroutine NcFile_Get_Var_r_2d( self, description, values, units, status, &
		start, count, missing_value )
		vstart, vcount, missing_value )

		use NetCDF, only : NF90_Inquire_Dimension
		use NetCDF, only : NF90_Inquire_Variable, NF90_Get_Var
		@@ -3436,7 +3508,7 @@ contains
		real, intent(out) :: values(:,:)
		character(len=*), intent(out) :: units
		integer, intent(out) :: status
		integer, intent(in), optional :: start(:), count(:)
		integer, intent(in), optional :: vstart(:), vcount(:)
		real, intent(out), optional :: missing_value

		! --- const --------------------------------------
		@@ -3460,28 +3532,28 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! check start:
		if ( any((/present(start),present(count)/)) .and. &
		(.not. all((/present(start),present(count)/))) ) then
		if ( any((/present(vstart),present(vcount)/)) .and. &
		(.not. all((/present(vstart),present(vcount)/))) ) then
		write (csol,'("specify both start and count")'); call csoErr
		TRACEBACK; status=1; return
		end if

		! combine slabs?
		combine = .false.
		if ( present(start) ) combine = start(1) < 1
		if ( present(vstart) ) combine = vstart(1) < 1
		! switch:
		if ( combine ) then

		! storage:
		allocate( xstart(size(start)), stat=status )
		allocate( xstart(size(vstart)), stat=status )
		IF_NOT_OK_RETURN(status=1)
		allocate( xcount(size(count)), stat=status )
		allocate( xcount(size(vcount)), stat=status )
		IF_NOT_OK_RETURN(status=1)
		allocate( dimids(size(count)), stat=status )
		allocate( dimids(size(vcount)), stat=status )
		IF_NOT_OK_RETURN(status=1)
		! copy:
		xstart = start
		xcount = count
		xstart = vstart
		xcount = vcount

		! start index:
		x1 = xstart(1)
		@@ -3536,17 +3608,17 @@ contains
		else

		! read:
		status = NF90_Get_Var( self%ncid, varid, values, start=start, count=count )
		status = NF90_Get_Var( self%ncid, varid, values, start=vstart, count=vcount )
		if ( status /= NF90_NOERR ) then
		csol=NF90_StrError(status); call csoErr
		write (csol,'("while reading:")'); call csoErr
		write (csol,'(" file name : ",a)') trim(self%filename); call csoErr
		write (csol,'(" variable description : ",a)') trim(description); call csoErr
		if ( present(start) ) then
		write (csol,*) ' start : ', start; call csoErr
		if ( present(vstart) ) then
		write (csol,*) ' start : ', vstart; call csoErr
		end if
		if ( present(count) ) then
		write (csol,*) ' count : ', count; call csoErr
		if ( present(vcount) ) then
		write (csol,*) ' count : ', vcount; call csoErr
		end if
		TRACEBACK; status=1; return
		end if
		@@ -3579,7 +3651,7 @@ contains
		! *

		subroutine NcFile_Get_Var_r_3d( self, description, values, units, status, &
		start, count, missing_value )
		vstart, vcount, missing_value )

		use NetCDF, only : NF90_Inquire_Dimension
		use NetCDF, only : NF90_Inquire_Variable, NF90_Get_Var
		@@ -3593,7 +3665,7 @@ contains
		real, intent(out) :: values(:,:,:)
		character(len=*), intent(out) :: units
		integer, intent(out) :: status
		integer, intent(in), optional :: start(:), count(:)
		integer, intent(in), optional :: vstart(:), vcount(:)
		real, intent(out), optional :: missing_value

		! --- const --------------------------------------
		@@ -3617,28 +3689,28 @@ contains
		IF_NOT_OK_RETURN(status=1)

		! check start:
		if ( any((/present(start),present(count)/)) .and. &
		(.not. all((/present(start),present(count)/))) ) then
		if ( any((/present(vstart),present(vcount)/)) .and. &
		(.not. all((/present(vstart),present(vcount)/))) ) then
		write (csol,'("specify both start and count")'); call csoErr
		TRACEBACK; status=1; return
		end if

		! combine slabs?
		combine = .false.
		if ( present(start) ) combine = start(1) < 1
		if ( present(vstart) ) combine = vstart(1) < 1
		! switch:
		if ( combine ) then

		! storage:
		allocate( xstart(size(start)), stat=status )
		allocate( xstart(size(vstart)), stat=status )
		IF_NOT_OK_RETURN(status=1)
		allocate( xcount(size(count)), stat=status )
		allocate( xcount(size(vcount)), stat=status )
		IF_NOT_OK_RETURN(status=1)
		allocate( dimids(size(count)), stat=status )
		allocate( dimids(size(vcount)), stat=status )
		IF_NOT_OK_RETURN(status=1)
		! copy:
		xstart = start
		xcount = count
		xstart = vstart
		xcount = vcount

		! start index:
		x1 = xstart(1)
		@@ -3693,17 +3765,17 @@ contains
		else

		! read:
		status = NF90_Get_Var( self%ncid, varid, values, start=start, count=count )
		status = NF90_Get_Var( self%ncid, varid, values, start=vstart, count=vcount )
		if ( status /= NF90_NOERR ) then
		csol=NF90_StrError(status); call csoErr
		write (csol,'("while reading:")'); call csoErr
		write (csol,'(" file name : ",a)') trim(self%filename); call csoErr
		write (csol,'(" variable description : ",a)') trim(description); call csoErr
		if ( present(start) ) then
		write (csol,*) ' start : ', start; call csoErr
		if ( present(vstart) ) then
		write (csol,*) ' start : ', vstart; call csoErr
		end if
		if ( present(count) ) then
		write (csol,*) ' count : ', count; call csoErr
		if ( present(vcount) ) then
		write (csol,*) ' count : ', vcount; call csoErr
		end if
		TRACEBACK; status=1; return
		end if

oper/src/cso_pixels.F90

+79 −8

Original line number	Diff line number	Diff line
		@@ -25,6 +25,12 @@
		! 2023-01, Arjo Segers
		! Support integer(1) and character variables.
		!
		! 2025-09, Arjo Segers
		! Support fill values when reading 1D pixel array,
		! to support for example files where variables have pixels either over land or water.
		!
		! 2025-10, Arjo Segers
		! Support kernel convolution 'xa + A ( x - xa )'.
		!
		!###############################################################################
		!
		@@ -721,8 +727,10 @@ contains
		! storage for all data:
		allocate( data0(xshp(1)), stat=status )
		IF_NOT_OK_RETURN(status=1)
		! read data and units from variable defined by description:
		call ncf%Get_Var( description, data0, dunits, status )
		! read data and units from variable defined by description;
		! set no-data values to fill value:
		call ncf%Get_Var( description, data0, dunits, status, &
		fill_value=self%fill_value )
		IF_NOT_OK_RETURN(status=1)
		else
		! dummy ...
		@@ -743,8 +751,10 @@ contains
		! storage for all data:
		allocate( data0(xshp(1)), stat=status )
		IF_NOT_OK_RETURN(status=1)
		! read data and units from variable defined by description:
		call ncf%Get_Var( description, data0, dunits, status )
		! read data and units from variable defined by description;
		! set no-data values to fill value:
		call ncf%Get_Var( description, data0, dunits, status, &
		fill_value=self%fill_value )
		IF_NOT_OK_RETURN(status=1)
		! copy selections:
		call slc%Copy( data0, self%data0, status )
		@@ -938,8 +948,10 @@ contains
		case ( '0D' )
		! read here?
		if ( self%read_by_me ) then
		! read data and units from variable defined by description:
		call ncf%Get_Var( description, self%data0, dunits, status )
		! read data and units from variable defined by description;
		! reset no-data values to fill_value:
		call ncf%Get_Var( description, self%data0, dunits, status, &
		fill_value=self%fill_value )
		IF_NOT_OK_RETURN(status=1)
		end if ! read by me
		! need to broadcast?
		@@ -4686,6 +4698,8 @@ contains
		character(len=64) :: A_units
		real, pointer :: A_m_data(:,:,:) ! (nretr,nlayer,npix)
		character(len=64) :: A_m_units
		real, pointer :: xa_data(:,:) ! (nretr,npix)
		character(len=64) :: xa_units
		real, pointer :: x_data(:,:) ! (nlayer,npix)
		character(len=64) :: x_units
		real, pointer :: Sx_data(:,:) ! (nlayer,npix)
		@@ -4978,6 +4992,63 @@ contains
		hp_data(hp_itop,ipix) = 200.0e2 ! Pa
		end do ! ipix

		!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		!~ kernel convolution:
		case ( 'xa + A ( x - xa )' )
		!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		! pointer to target array:
		call self%GetData( status, id=id, data1=y_data, units=y_units )
		IF_NOT_OK_RETURN(status=1)
		! get pointers to source arrays:
		call self%GetFormulaData( p%formula_terms, 'A', status, pd=pd, data2=A_data, units=A_units )
		IF_NOT_OK_RETURN(status=1)
		call self%GetFormulaData( p%formula_terms, 'xa', status, pd=pd, data1=xa_data, units=xa_units, fill_value=fill_value )
		IF_NOT_OK_RETURN(status=1)
		call self%GetFormulaData( p%formula_terms, 'x', status, pd=pd, data1=x_data, units=x_units, fill_value=fill_value )
		IF_NOT_OK_RETURN(status=1)
		! check units:
		if ( trim(A_units) /= '1' ) then
		write (csol,'("A units `",a,"` should be `1`")') trim(A_units); call csoErr
		write (csol,'(" formula : ",a)') trim(p%formula); call csoErr
		write (csol,'(" formula_terms : ",a)') trim(p%formula_terms); call csoErr
		write (csol,'(" variable : ",a)') trim(p%name); call csoErr
		TRACEBACK; status=1; return
		end if
		if ( trim(xa_units) /= trim(y_units) ) then
		write (csol,'("output units `",a,"` should be equal to xa units `",a,"`")') trim(y_units), trim(xa_units); call csoErr
		write (csol,'(" formula : ",a)') trim(p%formula); call csoErr
		write (csol,'(" formula_terms : ",a)') trim(p%formula_terms); call csoErr
		write (csol,'(" variable : ",a)') trim(p%name); call csoErr
		TRACEBACK; status=1; return
		end if
		if ( trim(x_units) /= trim(y_units) ) then
		write (csol,'("output units `",a,"` should be equal to x units `",a,"`")') trim(y_units), trim(x_units); call csoErr
		write (csol,'(" formula : ",a)') trim(p%formula); call csoErr
		write (csol,'(" formula_terms : ",a)') trim(p%formula_terms); call csoErr
		write (csol,'(" variable : ",a)') trim(p%name); call csoErr
		TRACEBACK; status=1; return
		end if
		! apply:
		do ipix = 1, npix
		! filter on no-data ..
		if ( x_data(1,ipix) == fill_value ) cycle
		! with some compilers problem using "matmul"; instead,
		! loop over target layers:
		do iretr = 1, size(y_data,1)
		! fill with dot procuct:
		y_data(iretr,ipix) = xa_data(iretr,ipix) + &
		sum( A_data(iretr,:,ipix) * ( x_data(:,ipix) - xa_data(:,ipix) ) )
		!! testing ...
		!write (csol,*) 'xxx ipix = ', ipix, '; iretr = ', iretr; call csoPr
		!write (csol,*) '..x A = ', A_data(iretr,1:5,ipix); call csoPr
		!write (csol,*) '..x x = ', x_data(1:5,ipix); call csoPr
		!write (csol,*) '..x xnan = ', fill_value, x_data(1,ipix) == fill_value; call csoPr
		!write (csol,*) '..x y = ', y_data(:,ipix); call csoPr
		!TRACEBACK; status=1; return
		end do ! iretr
		end do ! ipix

		!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		!~ kernel convolution, no apriori
		case ( 'A x' )

oper/src/cso_tools.F90

+105 −1

File changed.

Preview size limit exceeded, changes collapsed.