From 0d46f836d573e4fd10991f3dcd649b76c1c74c68 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 16:05:19 +0200 Subject: [PATCH 01/20] Updated change log. --- CHANGELOG | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index 78cdbe4..41efa1d 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -600,3 +600,10 @@ Convert multiple files if orbit is stored in patches. Check on duplicate output src/cso/cso_file.py src/cso/cso_s5p.py config/tutorial/tutorial.rc + + +vx.y +---- + +Updated Fortran operator. + oper/ -- GitLab From 7c80f1ab1f538fb9fe8344ab7eec354ca053f381 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Fri, 20 Feb 2026 09:09:44 +0100 Subject: [PATCH 02/20] Updated Fortran operator. --- oper/CHANGELOG | 21 +++++ oper/src/cso_ncfile.F90 | 194 +++++++++++++++++++++++++++------------- oper/src/cso_pixels.F90 | 87 ++++++++++++++++-- oper/src/cso_tools.F90 | 106 +++++++++++++++++++++- 4 files changed, 338 insertions(+), 70 deletions(-) diff --git a/oper/CHANGELOG b/oper/CHANGELOG index d9ff3a9..7e51972 100644 --- a/oper/CHANGELOG +++ b/oper/CHANGELOG @@ -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 + +! diff --git a/oper/src/cso_ncfile.F90 b/oper/src/cso_ncfile.F90 index 1644fb4..393e296 100644 --- a/oper/src/cso_ncfile.F90 +++ b/oper/src/cso_ncfile.F90 @@ -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,7 +3445,10 @@ contains integer :: varid real :: add_offset, scale_factor - + logical :: packed + real :: fill_value_in + logical :: filled + ! --- begin ---------------------------------- ! get variable id: @@ -3393,23 +3456,32 @@ 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: - values = add_offset + scale_factor * values + 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: call self%Inq_VarUnits( varid, description, units, status ) IF_ERROR_RETURN(status=1) @@ -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 diff --git a/oper/src/cso_pixels.F90 b/oper/src/cso_pixels.F90 index 435d5a8..6c5e230 100644 --- a/oper/src/cso_pixels.F90 +++ b/oper/src/cso_pixels.F90 @@ -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 )'. ! !############################################################################### ! @@ -592,7 +598,7 @@ contains real, allocatable :: data2(:,:,:) ! (mv,nv,npix) ! --- begin ---------------------------------- - + ! by default perform broadcasts in case of reading on root, ! but skip if only this pe should read: with_bcast = .true. @@ -648,7 +654,7 @@ contains call csoc%BCast( csoc%root_id, xshp, status ) IF_NOT_OK_RETURN(status=1) end if ! broadcast - + ! selection? if ( present(slc) ) then @@ -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' ) diff --git a/oper/src/cso_tools.F90 b/oper/src/cso_tools.F90 index 46a7536..691e7f2 100644 --- a/oper/src/cso_tools.F90 +++ b/oper/src/cso_tools.F90 @@ -12,6 +12,9 @@ ! Added 'PointInsideDomain' to check if pixel center is inside domain. ! Added 'PolygonOverlapsDomain' to check if pixel footprint overlaps with domain. ! +! 2025-10, Arjo Segers +! Added 'LinInterpWeights' routine to calculate interpolation weights. +! ! !############################################################################### ! @@ -54,6 +57,7 @@ module CSO_Tools public :: LonLatRectangleArea public :: GetPolygonPoints public :: LinInterp + public :: LinInterpWeights public :: PointInsideDomain public :: PolygonOverlapsDomain @@ -600,7 +604,8 @@ contains ! Coordinate x should be strictly increasing or decreasing. ! - subroutine LinInterp( x, y, xp, yp, status, extrapolate ) + subroutine LinInterp( x, y, xp, yp, status, & + extrapolate ) ! --- in/out --------------------------------- @@ -705,6 +710,105 @@ contains status = 0 end subroutine LinInterp + ! * + + ! Return 'ii(1:2)' and 'ww(1:2)' with indices and weights of interpolation points. + + subroutine LinInterpWeights( x, xp, ii, ww, status, & + extrapolate ) + + ! --- in/out --------------------------------- + + real, intent(in) :: x(:) ! (n) + real, intent(in) :: xp + integer, intent(out) :: ii(2) + real, intent(out) :: ww(2) + integer, intent(out) :: status + logical, intent(in), optional :: extrapolate + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/LinInterpWeights' + + ! --- 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 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 ( xdir*xp < xdir*x(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 ( xdir*xp > xdir*x(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 ( xdir*xp <= xdir*x(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) ) + + ! return indices: + ii = (/ ip, ip+1 /) + ww = (/ x(ip+1) - xp, xp - x(ip) /) / ( x(ip+1) - x(ip) ) + + ! ok + status = 0 + + end subroutine LinInterpWeights ! *** -- GitLab From 4b6d249945504d53005318b676bc53307046b1b7 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 15:54:59 +0200 Subject: [PATCH 03/20] Support latest servers. --- oper/launcher | 26 ++++--- oper/launcher-mpi | 67 ++++++++++++------- oper/src/Makefile_config_metno-fahrenheit | 34 ++++++++++ oper/src/Makefile_config_metno-nebula | 40 ----------- ...nfig_tno-hpc3 => Makefile_config_tno-hpc4} | 0 5 files changed, 90 insertions(+), 77 deletions(-) create mode 100644 oper/src/Makefile_config_metno-fahrenheit delete mode 100644 oper/src/Makefile_config_metno-nebula rename oper/src/{Makefile_config_tno-hpc3 => Makefile_config_tno-hpc4} (100%) diff --git a/oper/launcher b/oper/launcher index d473310..4330489 100755 --- a/oper/launcher +++ b/oper/launcher @@ -20,7 +20,7 @@ exe='tutorial_oper_S5p.x' case `hostname` in #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # TNO / HPC3: + # TNO / HPC4: app-hpc* ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -64,20 +64,18 @@ case `hostname` in ;; #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # MET Norway / Nebula: - nebula.* ) + # MET Norway / Fahrenheit: + fahrenheit* ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # module commands: - module --force purge - module load "mpprun/4.0" - module load "nsc/.1.0" - module load "buildenv-intel/2018a-eb" - module load "netCDF/4.4.1.1-HDF5-1.8.19-nsc1-intel-2018a-eb" - #module load "Python/3.7.0-anaconda-5.3.0-extras-nsc1" + module purge + module load "buildenv-intel/2023.1.0-hpc1" + module load "netCDF-HDF5/4.9.2-1.12.2-hpc1" + module load "makedepf90/2.8.8" # machine specific settings: - MACHINE="metno-nebula" + MACHINE="metno-fahrenheit" ;; @@ -102,7 +100,7 @@ esac cd src # remove old object files etc: -make clean +#make clean # replace include file if necessary: if ( grep '#define _MPI' 'cso.inc' > /dev/null ) ; then @@ -112,8 +110,8 @@ if ( grep '#define _MPI' 'cso.inc' > /dev/null ) ; then EOF fi -## update dependencies, this requires 'makedepf90' -#make deps +# update dependencies, this requires 'makedepf90' +make deps # make executable: make ${exe} MACHINE=${MACHINE} @@ -131,7 +129,7 @@ echo "" echo "run ..." # cleanup: -make clean +#make clean # run: ./src/${exe} diff --git a/oper/launcher-mpi b/oper/launcher-mpi index fefb48e..6e9bd3b 100755 --- a/oper/launcher-mpi +++ b/oper/launcher-mpi @@ -24,20 +24,16 @@ NPE=4 case `hostname` in #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # TNO / HPC3: + # TNO / HPC4: app-hpc* ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # module commands: module purge - module load "slurm/18.08.8" - module load "gcc/8.2.0" - module load "gcc-suite/8.2.0" - module load "openmpi/4.0.5" - module load "curl/default" - module load "openssl/default" - module load "netcdf-c/4.7.4" - module load "netcdf-fortran/4.5.3" + module load "gcc-suite/9.5.0" + module load "openmpi/4.1.3" + module load "netcdf-c/4.9.2" + module load "netcdf-fortran/4.6.1" module load "makedepf90" # machine specific settings: @@ -46,20 +42,45 @@ case `hostname` in ;; #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # MET Norway / Nebula: - nebula.* ) + # MET Norway / PPI: + ppi-* ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - + + # compiler version and key used in module names: + export COMPILER_VERSION="intelPE2018" ; compkey="intel2018" + # modules do not export installation paths ... + #~ HDF5 libraray: + export HDF5_VERSION="1.10.5-${compkey}" + export HDF5_HOME="/modules/centos7/hdf5/${HDF5_VERSION}" + #~ NetCDF libraray: + export NETCDF_VERSION="4.7.0-${compkey}" + export NETCDF_HOME="/modules/centos7/netcdf/${NETCDF_VERSION}" + # module commands: module --force purge - module load "mpprun/4.0" - module load "nsc/.1.0" - module load "buildenv-intel/2018a-eb" - module load "netCDF/4.4.1.1-HDF5-1.8.19-nsc1-intel-2018a-eb" - #module load "Python/3.7.0-anaconda-5.3.0-extras-nsc1" + module load "compiler/${COMPILER_VERSION}" + module load "hdf5/${HDF5_VERSION}" + module load "netcdf/${NETCDF_VERSION}" + module load "makedepf90/default" + + # machine specific settings: + MACHINE="metno-ppi" + + ;; + + #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # MET Norway / Fahrenheit: + fahrenheit* ) + #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + # module commands: + module purge + module load "buildenv-intel/2023.1.0-hpc1" + module load "netCDF-HDF5/4.9.2-1.12.2-hpc1" + module load "makedepf90/2.8.8" # machine specific settings: - MACHINE="metno-nebula" + MACHINE="metno-fahrenheit" ;; @@ -95,8 +116,8 @@ if ! ( grep '#define _MPI' 'cso.inc' > /dev/null ) ; then EOF fi -## update dependencies: -#make deps +# update dependencies, this requires 'makedepf90' +make deps # make executable: make ${exe} MACHINE=${MACHINE} MPI=yes @@ -122,7 +143,7 @@ jbfile="${jbname}.jb" case `hostname` in #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # TNO / HPC3: + # TNO / HPC4: app-hpc* ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -161,8 +182,8 @@ EOF ;; #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # MET Norway / Nebula: - nebula.* ) + # MET Norway / Fahrenheit: + fahrenheit* ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # fill: diff --git a/oper/src/Makefile_config_metno-fahrenheit b/oper/src/Makefile_config_metno-fahrenheit new file mode 100644 index 0000000..86cc469 --- /dev/null +++ b/oper/src/Makefile_config_metno-fahrenheit @@ -0,0 +1,34 @@ +# +# Makefile configuration. +# +# Environment: +# MPI = yes|no # if "yes" use MPI wrappers +# + +# compiler: +ifeq ($(MPI),yes) + #~ parallel: + FC = mpiifort +else + #~ serial: + FC = ifort +endif + +# default flags, similar to what is used in EMEP: +FFLAGS_DEFAULT = -extend-source -recursive -r8 + +# optimization flags: +FFLAGS_OPTIM = -O0 -g -traceback -check all -check noarg_temp_created -debug-parameters all -ftrapuv -fpe0 +#FFLAGS_OPTIM = -O3 + +# combine: +FFLAGS = $(FFLAGS_DEFAULT) $(FFLAGS_OPTIM) +LDFLAGS = + +# single source for NETCDF libs: +NETCDF_INCS = -I${NETCDF_DIR}/include +NETCDF_LIBS = -L${NETCDF_DIR}/lib -lnetcdf -lnetcdff -Wl,-rpath -Wl,${NETCDF_DIR}/lib + +# collect: +INCS = $(NETCDF_INCS) +LIBS = $(NETCDF_LIBS) diff --git a/oper/src/Makefile_config_metno-nebula b/oper/src/Makefile_config_metno-nebula deleted file mode 100644 index dd237f3..0000000 --- a/oper/src/Makefile_config_metno-nebula +++ /dev/null @@ -1,40 +0,0 @@ -# -# Makefile configuration. -# -# Environment: -# MPI = yes|no # if "yes" use MPI wrappers -# - -# compiler: -ifeq ($(MPI),yes) - #~ parallel: - FC = mpiifort -else - #~ serial: - FC = ifort -endif - -# default flags, similar to what is used in EMEP: -FFLAGS_DEFAULT = -extend-source -fpp -recursive -r8 -convert big_endian -shared-intel -assume noold_maxminloc -ftz - -# optimization flags: -FFLAGS_OPTIM = -O0 -g -traceback -check all -check noarg_temp_created -debug-parameters all -ftrapuv -fpe0 -#FFLAGS_OPTIM = -O3 - -# combine: -FFLAGS = $(FFLAGS_DEFAULT) $(FFLAGS_OPTIM) -LDFLAGS = - -# library: -NETCDF_FORTRAN_HOME = ${EBROOTNETCDF} -NETCDF_FORTRAN_INCS = -I${NETCDF_FORTRAN_HOME}/include -NETCDF_FORTRAN_LIBS = -L${NETCDF_FORTRAN_HOME}/lib -lnetcdff -Wl,-rpath -Wl,${NETCDF_FORTRAN_HOME}/lib - -# library: -NETCDF_C_HOME = ${EBROOTNETCDF} -NETCDF_C_INCS = -NETCDF_C_LIBS = -L${NETCDF_C_HOME}/lib -lnetcdf -Wl,-rpath -Wl,${NETCDF_C_HOME}/lib - -# collect: -INCS = $(NETCDF_FORTRAN_INCS) $(NETCDF_C_INCS) -LIBS = $(NETCDF_FORTRAN_LIBS) $(NETCDF_C_LIBS) diff --git a/oper/src/Makefile_config_tno-hpc3 b/oper/src/Makefile_config_tno-hpc4 similarity index 100% rename from oper/src/Makefile_config_tno-hpc3 rename to oper/src/Makefile_config_tno-hpc4 -- GitLab From a2f39ed4cda227dc47c0f72337b6f93837571fe3 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 15:55:36 +0200 Subject: [PATCH 04/20] Added routines to broadcast/scatter 1D integer or character arrays. --- oper/src/cso_comm.F90 | 192 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 168 insertions(+), 24 deletions(-) diff --git a/oper/src/cso_comm.F90 b/oper/src/cso_comm.F90 index b4c63f1..6261769 100644 --- a/oper/src/cso_comm.F90 +++ b/oper/src/cso_comm.F90 @@ -13,13 +13,16 @@ ! https://computing.llnl.gov/tutorials/mpi/ ! ! -! HISTORY +! CHANGES ! -! 2023-01, Arjo Segers -! Support integer(1) and character variables. +! 2023-01, Arjo Segers +! Support integer(1) and character variables. ! -! 2025-02, Arjo Segers -! Removed comment from 'use MPI' statements after compiler complained. +! 2025-02, Arjo Segers +! Removed comment from 'use MPI' statements after compiler complained. +! +! 2026-07, Arjo Segers +! Added routine to broadcast 1D character arrays. ! !### macro's ########################################################### ! @@ -145,6 +148,7 @@ module CSO_Comm procedure :: CSO_Comm_BCast_r8_3d procedure :: CSO_Comm_BCast_r8_4d procedure :: CSO_Comm_BCast_c + procedure :: CSO_Comm_BCast_c_1d generic :: BCast => CSO_Comm_BCast_i, & CSO_Comm_BCast_i1_1d, & CSO_Comm_BCast_i_1d, & @@ -158,7 +162,8 @@ module CSO_Comm CSO_Comm_BCast_r8_2d, & CSO_Comm_BCast_r8_3d, & CSO_Comm_BCast_r8_4d, & - CSO_Comm_BCast_c + CSO_Comm_BCast_c, & + CSO_Comm_BCast_c_1d ! procedure :: CSO_Comm_Gather_i procedure :: CSO_Comm_Gather_i_1d @@ -175,8 +180,8 @@ module CSO_Comm CSO_Comm_AllGather_r1 ! procedure CSO_Comm_GatherV_i1_1d - procedure CSO_Comm_GatherV_i_1d - procedure CSO_Comm_GatherV_i_2d + procedure CSO_Comm_GatherV_i4_1d + procedure CSO_Comm_GatherV_i4_2d procedure CSO_Comm_GatherV_r4_1d procedure CSO_Comm_GatherV_r8_1d procedure CSO_Comm_GatherV_c1_2d @@ -185,8 +190,8 @@ module CSO_Comm procedure CSO_Comm_GatherV_r4_3d procedure CSO_Comm_GatherV_r8_3d generic :: GatherV => CSO_Comm_GatherV_i1_1d, & - CSO_Comm_GatherV_i_1d, & - CSO_Comm_GatherV_i_2d, & + CSO_Comm_GatherV_i4_1d, & + CSO_Comm_GatherV_i4_2d, & CSO_Comm_GatherV_r4_1d, & CSO_Comm_GatherV_r8_1d, & CSO_Comm_GatherV_c1_2d, & @@ -201,6 +206,7 @@ module CSO_Comm CSO_Comm_Gather2D_r8 ! procedure CSO_Comm_ScatterV_i1_1d + procedure CSO_Comm_ScatterV_i4_1d procedure CSO_Comm_ScatterV_r4_1d procedure CSO_Comm_ScatterV_c1_2d procedure CSO_Comm_ScatterV_r4_2d @@ -209,6 +215,7 @@ module CSO_Comm procedure CSO_Comm_ScatterV_r8_2d procedure CSO_Comm_ScatterV_r8_3d generic :: ScatterV => CSO_Comm_ScatterV_i1_1d, & + CSO_Comm_ScatterV_i4_1d, & CSO_Comm_ScatterV_r4_1d, & CSO_Comm_ScatterV_c1_2d, & CSO_Comm_ScatterV_r4_2d, & @@ -2287,6 +2294,50 @@ contains end subroutine CSO_Comm_BCast_c + ! * + + + subroutine CSO_Comm_BCast_c_1d( self, rootid, values, status ) + +#ifdef _MPI + use MPI_F08, only : MPI_DataType + use MPI_F08, only : MPI_BCast +#endif + + ! --- in/out --------------------------------- + + class(T_CSO_Comm), intent(in) :: self + integer, intent(in) :: rootid + character(len=*), intent(inout) :: values(:) + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Comm_BCast_c_1d' + + ! --- local ---------------------------------- + +#ifdef _MPI + type(MPI_DataType) :: dtype +#endif + + ! --- begin ---------------------------------- + +#ifdef _MPI + ! data type: + call self%GetDataType( 'char', 1, dtype, status ) + IF_NOT_OK_RETURN(status=1) + ! send values from root to all other pe's: + call MPI_BCast( values, len(values(1))*size(values), dtype, rootid, self%comm, ierror=status ) + IF_MPI_NOT_OK_RETURN(status=1) +#endif + + ! ok + status = 0 + + end subroutine CSO_Comm_BCast_c_1d + + ! ******************************************************************** ! *** ! *** gather @@ -2717,7 +2768,7 @@ contains ! * - subroutine CSO_Comm_GatherV_i_1d( self, send, recv, status, & + subroutine CSO_Comm_GatherV_i4_1d( self, send, recv, status, & nloc ) #ifdef _MPI @@ -2728,15 +2779,15 @@ contains ! --- in/out --------------------------------- class(T_CSO_Comm), intent(in) :: self - integer, intent(in) :: send(:) ! (max(1,nloc)) - integer, intent(out) :: recv(:) ! (sum nloc) + integer(4), intent(in) :: send(:) ! (max(1,nloc)) + integer(4), intent(out) :: recv(:) ! (sum nloc) integer, intent(out) :: status integer, intent(in), optional :: nloc ! --- const ---------------------------------- - character(len=*), parameter :: rname = mname//'/CSO_Comm_GatherV_i_1d' + character(len=*), parameter :: rname = mname//'/CSO_Comm_GatherV_i4_1d' ! --- local ---------------------------------- @@ -2744,8 +2795,8 @@ contains integer :: ntot #ifdef _MPI type(MPI_DataType) :: dtype - integer, allocatable :: recvcounts(:) ! (npes) - integer, allocatable :: displs(:) ! (npes) + integer(4), allocatable :: recvcounts(:) ! (npes) + integer(4), allocatable :: displs(:) ! (npes) #endif ! --- begin ---------------------------------- @@ -2803,11 +2854,11 @@ contains ! ok status = 0 - end subroutine CSO_Comm_GatherV_i_1d + end subroutine CSO_Comm_GatherV_i4_1d ! * - subroutine CSO_Comm_GatherV_i_2d( self, send, recv, status, & + subroutine CSO_Comm_GatherV_i4_2d( self, send, recv, status, & nloc ) #ifdef _MPI @@ -2818,15 +2869,15 @@ contains ! --- in/out --------------------------------- class(T_CSO_Comm), intent(in) :: self - integer, intent(in) :: send(:,:) ! (m,max(1,nloc)) - integer, intent(out) :: recv(:,:) ! (m,sum nloc) + integer(4), intent(in) :: send(:,:) ! (m,max(1,nloc)) + integer(4), intent(out) :: recv(:,:) ! (m,sum nloc) integer, intent(out) :: status integer, intent(in), optional :: nloc ! --- const ---------------------------------- - character(len=*), parameter :: rname = mname//'/CSO_Comm_GatherV_i_2d' + character(len=*), parameter :: rname = mname//'/CSO_Comm_GatherV_i4_2d' ! --- local ---------------------------------- @@ -2835,8 +2886,8 @@ contains #ifdef _MPI integer :: ntot type(MPI_DataType) :: dtype - integer, allocatable :: recvcounts(:) ! (npes) - integer, allocatable :: displs(:) ! (npes) + integer(4), allocatable :: recvcounts(:) ! (npes) + integer(4), allocatable :: displs(:) ! (npes) #endif ! --- begin ---------------------------------- @@ -2899,7 +2950,7 @@ contains ! ok status = 0 - end subroutine CSO_Comm_GatherV_i_2d + end subroutine CSO_Comm_GatherV_i4_2d ! * @@ -3944,6 +3995,99 @@ contains ! * + + subroutine CSO_Comm_ScatterV_i4_1d( self, send, recv, status, & + nloc ) + +#ifdef _MPI + use MPI_F08, only : MPI_DataType + use MPI_F08, only : MPI_ScatterV +#endif + + ! --- in/out --------------------------------- + + integer, parameter :: wp = 4 + + class(T_CSO_Comm), intent(in) :: self + integer(wp), intent(in) :: send(:) ! (sum nloc) + integer(wp), intent(out) :: recv(:) ! (max(1,nloc)) + integer, intent(out) :: status + + integer, intent(in), optional :: nloc + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Comm_ScatterV_i4_1d' + + ! --- local ---------------------------------- + + integer :: n +#ifdef _MPI + type(MPI_DataType) :: dtype + integer :: ntot + integer, allocatable :: sendcounts(:) ! (npes) + integer, allocatable :: displs(:) ! (npes) +#endif + + ! --- begin ---------------------------------- + + ! local size, take from optional argument if present (value is probably zero ..) + if ( present(nloc) ) then + n = nloc + else + n = size(recv) + end if + +#ifdef _MPI + + ! data type: + call self%GetDataType( 'integer', wp, dtype, status ) + IF_NOT_OK_RETURN(status=1) + + ! storage: + allocate( sendcounts(0:self%npes-1), stat=status ) + IF_NOT_OK_RETURN(status=1) + allocate( displs(0:self%npes-1), stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! collect numbers: + call self%ParInfo( n, status, ntot=ntot, sendcounts=sendcounts, displs=displs ) + IF_NOT_OK_RETURN(status=1) + + ! check ... + if ( self%root ) then + if ( size(send) /= ntot ) then + write (csol,'("send buffer has size ",i0," while ntot is ",i0)') size(send), ntot; call csoErr + TRACEBACK; status=1; return + end if + end if + + ! collect values from all pe's on root: + call MPI_ScatterV( send, sendcounts, displs, dtype, & + recv, n , dtype, & + self%root_id, self%comm, ierror=status ) + IF_MPI_NOT_OK_RETURN(status=1) + + ! clear: + deallocate( sendcounts, stat=status ) + IF_NOT_OK_RETURN(status=1) + deallocate( displs, stat=status ) + IF_NOT_OK_RETURN(status=1) + +#else + + ! just copy ... + if ( n > 0 ) recv(1:n) = send(1:n) + +#endif + + ! ok + status = 0 + + end subroutine CSO_Comm_ScatterV_i4_1d + + ! * + subroutine CSO_Comm_ScatterV_r4_1d( self, send, recv, status, & nloc ) -- GitLab From eb8476fc87c5f62bcc25ab4d7a4a33b1ef25ce86 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 15:55:47 +0200 Subject: [PATCH 05/20] Added `CompressTime` routine. --- oper/src/cso_datetime.F90 | 110 +++++++++++++++++++++++++++++++------- 1 file changed, 92 insertions(+), 18 deletions(-) diff --git a/oper/src/cso_datetime.F90 b/oper/src/cso_datetime.F90 index 291f456..2018cb3 100644 --- a/oper/src/cso_datetime.F90 +++ b/oper/src/cso_datetime.F90 @@ -232,7 +232,7 @@ ! %M : 2 digit minute ! %S : 2 digit second ! -! t = TDate( year=2021, month=9, day=6, hour=0, min=0 ) +! t = T_CSO_DateTime( year=2021, month=9, day=6, hour=0, min=0 ) ! line = 'file_%Y%m%d_%H%M.txt' ! call CSO_Format( line, t, status ) ! @@ -245,20 +245,27 @@ ! call CSO_DateTimeDefaults( [calendar='gregorian'] ) ! ! -! HISTORY +!############################################################################## +! +! CHANGES ! ! 2022-09, Arjo Segers ! Initialize some variables to avoid compiler warnings. ! ! 2023-09, Arjo Segers -! Changed "mili" to the more correct "milli". +! Changed `mili` to the more correct `milli`. +! +! 2023-10-31, Jacob van Peet +! Added date_Get to public module members to get compilation with nvfortran working. ! ! 2022-09, Arjo Segers ! Added `CSO_Format` routine. ! ! 2025-02, Arjo Segers -! Support labels 'hours', 'minutes', 'seconds', etc. +! Support labels `hours`, `minutes`, `seconds`, etc. ! +! 2026-07, Arjo Segers +! Added `CompressTime` routine. ! !### macro's ################################################################## ! @@ -294,6 +301,10 @@ module CSO_DateTimes public :: Check public :: Normalize + ! JvP 20231031: Added date_Get (which is already a member of "Get" interface) + ! to get compilation with nvfortran working. + public :: date_Get + public :: LeapYear public :: LeapDay public :: days_in_month, calc_days_in_month @@ -325,6 +336,7 @@ module CSO_DateTimes public :: CSO_ReadFromLine public :: ExpandTime + public :: CompressTime public :: Compare_Date_Num public :: Pretty @@ -2642,40 +2654,40 @@ contains ! count fractional months: incr_rTotal = dt%day + & dt%hour / 24.0 + & - dt%minute / 24.0 / 60.0 + & - dt%second / 24.0 / 60.0 / 60.0 + & - dt%millisecond / 24.0 / 60.0 / 60.0 / 1000.0 + dt%minute / 24.0 / 60.0 + & + dt%second / 24.0 / 60.0 / 60.0 + & + dt%millisecond / 24.0 / 60.0 / 60.0 / 1000.0 case ( 'hour', 'hours' ) ! 'nday' has been set to the total number of days from 0 to t; ! count fractional hours: incr_rTotal = dt%day * 24.0 + & dt%hour + & - dt%minute / 60.0 + & - dt%second / 60.0 / 60.0 + & - dt%millisecond / 60.0 / 60.0 / 1000.0 + dt%minute / 60.0 + & + dt%second / 60.0 / 60.0 + & + dt%millisecond / 60.0 / 60.0 / 1000.0 case ( 'min', 'minute', 'minutes' ) ! 'nday' has been set to the total number of days from 0 to t; ! count fractional minutes: incr_rTotal = dt%day * 24.0 * 60.0 + & dt%hour * 60.0 + & - dt%minute + & - dt%second / 60.0 + & - dt%millisecond / 60.0 / 1000.0 + dt%minute + & + dt%second / 60.0 + & + dt%millisecond / 60.0 / 1000.0 case ( 'sec', 'second', 'seconds' ) ! 'nday' has been set to the total number of days from 0 to t; ! count fractional seconds: incr_rTotal = dt%day * 24.0 * 60.0 * 60.0 + & dt%hour * 60.0 * 60.0 + & - dt%minute * 60.0 + & - dt%second + & - dt%millisecond / 1000.0 + dt%minute * 60.0 + & + dt%second + & + dt%millisecond / 1000.0 case ( 'milli', 'millisecond', 'milliseconds' ) ! 'nday' has been set to the total number of days from 0 to t; ! count fractional milli seconds: incr_rTotal = dt%day * 24.0 * 60.0 * 60.0 * 1000.0 + & dt%hour * 60.0 * 60.0 * 1000.0 + & - dt%minute * 60.0 * 1000.0 + & - dt%second * 1000.0 + & + dt%minute * 60.0 * 1000.0 + & + dt%second * 1000.0 + & dt%millisecond case default write (csol,'("do not know how to count time in unit : ",a)') trim(unit); call csoErr @@ -3370,6 +3382,68 @@ contains end subroutine ExpandTime + ! * + + + ! + ! Compress time to real value using units: + ! + ! hours since 1900-01-01 00:00:0.0 + ! + + subroutine CompressTime( t, units, value, status ) + + use CSO_String, only : CSO_ReadFromLine + + ! --- in/out --------------------------------- + + type(T_CSO_DateTime), intent(in) :: t + real, intent(out) :: value + character(len=*), intent(in) :: units + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CompressTime' + + ! --- local ---------------------------------- + + character(len=512) :: line + character(len=32) :: steps + character(len=32) :: dummy + type(T_CSO_DateTime) :: tref + + ! --- begin ---------------------------------- + + ! copy: + line = trim(units) + + ! extract step units: + call CSO_ReadFromLine( line, steps, status, sep=' ' ) + IF_NOT_OK_RETURN(status=1) + + ! extract 'since' + call CSO_ReadFromLine( line, dummy, status, sep=' ' ) + IF_NOT_OK_RETURN(status=1) + ! check ... + if ( trim(dummy) /= 'since' ) then + write (csol,'("second field should be `since`, not `",a,"`")') trim(dummy); call csoPr + TRACEBACK; status=1; return + end if + + ! extract ref time values: + call CSO_ReadFromLine( line, tref, status ) + IF_NOT_OK_RETURN(status=1) + + ! compute difference: + value = rTotal( t - tref, steps ) + + ! ok + status = 0 + + end subroutine CompressTime + + ! *** -- GitLab From b86cbfe703492684870ccf9a283e4170d5b6ec94 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 15:56:07 +0200 Subject: [PATCH 06/20] Added `CSO_GetBasename` and `CSO_SplitExt` routines. --- oper/src/cso_file.F90 | 110 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 104 insertions(+), 6 deletions(-) diff --git a/oper/src/cso_file.F90 b/oper/src/cso_file.F90 index e542526..d582be3 100644 --- a/oper/src/cso_file.F90 +++ b/oper/src/cso_file.F90 @@ -2,13 +2,16 @@ ! ! File tools. ! -! HISTORY +! CHANGES ! -! 2022-09, Arjo Segers -! Added `CSO_CheckDir` routine. +! 2022-09, Arjo Segers +! Added `CSO_CheckDir` routine. ! -! 2023-01, Arjo Segers -! Fixed problem with creation of output directories when using Intel compiler. +! 2023-01, Arjo Segers +! Fixed problem with creation of output directories when using Intel compiler. +! +! 2026-07, Arjo Segers +! Added `CSO_GetBasename` and `CSO_SplitExt` routines. ! !################################################################# ! @@ -32,6 +35,8 @@ module CSO_File public :: CSO_GetFU public :: CSO_GetDirname + public :: CSO_GetBasename + public :: CSO_SplitExt public :: CSO_CheckDir public :: T_CSO_TextFile @@ -157,7 +162,7 @@ contains ! --- const --------------------------- - character(len=*), parameter :: rname = mname//'/CSO_GetFU' + character(len=*), parameter :: rname = mname//'/CSO_GetDirname' ! --- local -------------------------- @@ -180,6 +185,99 @@ contains status = 0 end subroutine CSO_GetDirname + ! * + + + ! return filename part of path + + subroutine CSO_GetBasename( filename, basename, status ) + + ! --- in/out -------------------------- + + character(len=*), intent(in) :: filename + character(len=*), intent(out) :: basename + integer, intent(out) :: status + + ! --- const --------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_GetBasename' + + ! --- local -------------------------- + + integer :: n + integer :: i + + ! --- local --------------------------- + + ! size: + n = len_trim(filename) + ! search backward for path sep: + i = index( filename, '/', back=.true. ) + ! found? + if ( i > 0 ) then + ! part from seperator onwards: + if ( i < n ) then + basename = filename(i+1:n) + else + basename = '' + end if + else + ! no path ... + basename = trim(filename) + end if + + ! ok + status = 0 + + end subroutine CSO_GetBasename + + + ! * + + + ! split filename in basename and extension + + subroutine CSO_SplitExt( filename, rootname, ext, status ) + + ! --- in/out -------------------------- + + character(len=*), intent(in) :: filename + character(len=*), intent(out) :: rootname + character(len=*), intent(out) :: ext + integer, intent(out) :: status + + ! --- const --------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_SplitExt' + + ! --- local -------------------------- + + integer :: n + integer :: i + + ! --- local --------------------------- + + ! size: + n = len_trim(filename) + ! search backward for dot: + i = index( filename, '.', back=.true. ) + ! found? + if ( i > 0 ) then + ext = filename(i:n) + if ( i > 1 ) then + rootname = filename(1:i-1) + else + rootname = '' + end if + else + rootname = trim(filename) + ext = '' + end if + + ! ok + status = 0 + + end subroutine CSO_SplitExt ! * -- GitLab From 53e46ce3a770dc235140dfaeb945ba509ca56ff9 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 15:56:45 +0200 Subject: [PATCH 07/20] Support listings with second `filename`, this is for example used to specify both data and state files. Added method `Create` to initialize an empty listing that will be filled and saved. Added method `Select` to select records based on time interval. --- oper/src/cso_listing.F90 | 455 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 443 insertions(+), 12 deletions(-) diff --git a/oper/src/cso_listing.F90 b/oper/src/cso_listing.F90 index ac5a6b0..aa5d132 100644 --- a/oper/src/cso_listing.F90 +++ b/oper/src/cso_listing.F90 @@ -10,11 +10,15 @@ ! 2018/06/S5p_RPRO_NO2_03274.nc;2018-06-01T04:52:43.586000000;2018-06-01T04:59:12.377000000;03274 ! ! -! HISTORY +! CHANGES ! -! 2022-09, Arjo Segers -! Extended error messages. +! 2022-09, Arjo Segers +! Extended error messages. ! +! 2026-07, Arjo Segers +! Support listings with second `filename`, this is for example used to specify both data and state files. +! Added method `Create` to initialize an empty listing that will be filled and saved. +! Added method `Select` to select records based on time interval. ! !################################################################# ! @@ -50,6 +54,7 @@ module CSO_Listing type T_CSO_Record ! filename: character(len=256) :: filename + character(len=256) :: filename2 ! timerange: type(T_CSO_DateTime) :: t1, t2 ! average: @@ -57,6 +62,7 @@ module CSO_Listing contains procedure :: Init => Record_Init procedure :: Done => Record_Done + procedure :: Get => Record_Get procedure :: Match => Record_Match end type T_CSO_Record @@ -72,12 +78,23 @@ module CSO_Listing ! number of records: integer :: nrec ! storage: - type(T_CSO_Record), allocatable :: rec(:) ! (n) + type(T_CSO_Record), allocatable :: rec(:) ! (nrec) + ! selection flags: + logical, pointer :: selected(:) ! (nrec) + integer, pointer :: selection(:) ! (nrec) + integer :: nselect + ! contains procedure :: Init => Listing_Init procedure :: Done => Listing_Done procedure :: Show => Listing_Show + procedure :: Create => Listing_Create + procedure :: AddRecord => Listing_AddRecord + procedure :: PutOut => Listing_PutOut + procedure :: Get => Listing_Get + procedure :: GetRecord => Listing_GetRecord procedure :: SearchFile => Listing_SearchFile + procedure :: Select => Listing_Select end type T_CSO_Listing @@ -93,7 +110,8 @@ contains ! ============================================================== - subroutine Record_Init( self, filename, t1, t2, status ) + subroutine Record_Init( self, filename, t1, t2, status, & + filename2 ) use CSO_DateTimes, only : operator(+), operator(-), operator(*) @@ -103,6 +121,7 @@ contains character(len=*), intent(in) :: filename type(T_CSO_DateTime), intent(in) :: t1, t2 integer, intent(out) :: status + character(len=*), intent(in), optional :: filename2 ! --- const --------------------------- @@ -119,6 +138,12 @@ contains ! fill average: self%taver = self%t1 + 0.5*( self%t2 - self%t1 ) + ! extra? + if ( present(filename2) ) then + self%filename2 = trim(filename2) + else + self%filename2 = '' + end if ! ok status = 0 @@ -130,12 +155,12 @@ contains subroutine Record_Done( self, status ) - ! --- in/out ----------------- + ! --- in/out ----------------------- class(T_CSO_Record), intent(inout) :: self integer, intent(out) :: status - ! --- const ---------------------- + ! --- const ------------------------ character(len=*), parameter :: rname = mname//'/Record_Done' @@ -145,6 +170,35 @@ contains status = 0 end subroutine Record_Done + ! *** + + + subroutine Record_Get( self, status, filename, filename2, t1, t2 ) + + ! --- in/out ----------------------- + + class(T_CSO_Record), intent(in) :: self + integer, intent(out) :: status + character(len=*), intent(out), optional :: filename + character(len=*), intent(out), optional :: filename2 + type(T_CSO_DateTime), intent(out), optional :: t1, t2 + + ! --- const ------------------------ + + character(len=*), parameter :: rname = mname//'/Record_Get' + + ! --- begin ------------------------ + + ! copy: + if ( present(filename ) ) filename = trim(self%filename) + if ( present(filename2) ) filename2 = trim(self%filename2) + if ( present(t1 ) ) t1 = self%t1 + if ( present(t2 ) ) t2 = self%t2 + + ! ok + status = 0 + + end subroutine Record_Get ! *** @@ -153,7 +207,8 @@ contains ! Try if record can be assigned to target interval [t1,t2]. ! ! How the record is matched depends on the timekey: - ! 'aver' : record time average in (t1,t2] + ! 'aver' : record time average in (t1,t2] + ! 'overlap' : record time range overlaps with (t1,t2] ! ! Return value: ! match : .true. if matched @@ -183,10 +238,17 @@ contains ! switch: select case ( trim(timekey) ) - !~ average? + !~ check if record average time is in (t1,t2]: case ( 'aver' ) ! average should be in interval: match = (t1 < self%taver) .and. (self%taver <= t2) + !~ check if record range overlaps with (t1,t2] + case ( 'overlap' ) + ! compare boundaries: + ! record self%t1 self%t2 + ! o-----------------o + ! target (t1,t2] (t1,t2] (t1,t2] + match = (self%t1 <= t2) .and. (t1 < self%t2) !~ unknown ... case default @@ -210,6 +272,9 @@ contains ! === ! ============================================================== + ! + ! Create lising object from table file. + ! subroutine Listing_Init( self, filename, status ) @@ -236,6 +301,7 @@ contains character(len=1024) :: line0 integer :: irec character(len=1024) :: rec_filename + character(len=1024) :: rec_filename2 type(T_CSO_DateTime) :: rec_t1, rec_t2 ! --- begin ---------------------------- @@ -373,9 +439,17 @@ contains write (csol,'(" ",a)') trim(line); call csoErr TRACEBACK; status=1; return end if + !~ extra? + if ( len_trim(line) > 0 ) then + call CSO_ReadFromLine( line, rec_filename2, status, sep=self%sep ) + IF_NOT_OK_RETURN(status=1) + else + rec_filename2 = '' + end if ! store: - call self%rec(irec)%Init( rec_filename, rec_t1, rec_t2, status ) + call self%rec(irec)%Init( rec_filename, rec_t1, rec_t2, status, & + filename2=rec_filename2 ) IF_NOT_OK_RETURN(status=1) end do ! records @@ -388,6 +462,12 @@ contains TRACEBACK; status=1; return end if + ! storage for selection flags: + allocate( self%selected(self%nrec), source=.false., stat=status ) + IF_NOT_OK_RETURN(status=1) + allocate( self%selection(self%nrec), source=-999, stat=status ) + IF_NOT_OK_RETURN(status=1) + end if ! n>0 ! ok @@ -417,7 +497,7 @@ contains ! --- begin ------------------------ ! any? - if ( self%nrec > 0 ) then + if ( allocated(self%rec) ) then ! loop: do irec = 1, self%nrec ! clear: @@ -427,6 +507,13 @@ contains ! clear: deallocate( self%rec, stat=status ) IF_NOT_OK_RETURN(status=1) + ! clear: + if ( associated(self%selected) ) then + deallocate( self%selected, stat=status ) + IF_NOT_OK_RETURN(status=1) + deallocate( self%selection, stat=status ) + IF_NOT_OK_RETURN(status=1) + end if ! empty: self%nrec = 0 end if ! n>0 @@ -437,6 +524,167 @@ contains end subroutine Listing_Done + ! + ! Create empty lising object. + ! Specify target file to which content will be written. + ! Specify maximum number of records. + ! + + subroutine Listing_Create( self, filename, maxrec, status ) + + use CSO_File , only : CSO_GetDirname + + ! --- in/out ------------------------ + + class(T_CSO_Listing), intent(out) :: self + character(len=*), intent(in) :: filename + integer, intent(in) :: maxrec + integer, intent(out) :: status + + ! --- const --------------------------- + + character(len=*), parameter :: rname = mname//'/Listing_Create' + + ! --- local ---------------------------- + + ! --- begin ---------------------------- + + ! store: + self%filename = trim(filename) + + ! directory part: + call CSO_GetDirname( self%filename, self%dirname, status ) + IF_NOT_OK_RETURN(status=1) + + ! formatting: + self%sep = ';' + + ! storage: + allocate( self%rec(maxrec), stat=status ) + IF_NOT_OK_RETURN(status=1) + ! no records yet: + self%nrec = 0 + + ! ok + status = 0 + + end subroutine Listing_Create + + + ! *** + + + subroutine Listing_AddRecord( self, filename, t1, t2, status, & + filename2 ) + + ! --- in/out ------------------------ + + class(T_CSO_Listing), intent(inout) :: self + character(len=*), intent(in) :: filename + type(T_CSO_DateTime), intent(in) :: t1, t2 + integer, intent(out) :: status + character(len=*), intent(in), optional :: filename2 + + ! --- const --------------------------- + + character(len=*), parameter :: rname = mname//'/Listing_AddRecord' + + ! --- local ---------------------------- + + ! --- begin ---------------------------- + + ! increase counter: + self%nrec = self%nrec + 1 + ! check .. + if ( self%nrec > size(self%rec) ) then + write (csol,'("could not add record after maximum number of ",i0)') size(self%rec); call csoErr + TRACEBACK; status=1; return + end if + + ! fill: + call self%rec(self%nrec)%Init( filename, t1, t2, status, filename2=filename2 ) + IF_NOT_OK_RETURN(status=1) + + ! ok + status = 0 + + end subroutine Listing_AddRecord + + + ! *** + + + ! + ! Write listing file. + ! + + subroutine Listing_PutOut( self, status ) + + use CSO_Comm , only : csoc + use CSO_File , only : CSO_GetFU + use CSO_File , only : CSO_CheckDir + use CSO_DateTimes, only : Pretty + + ! --- in/out ------------------------ + + class(T_CSO_Listing), intent(in) :: self + integer, intent(out) :: status + + ! --- const --------------------------- + + character(len=*), parameter :: rname = mname//'/Listing_PutOut' + + ! --- local ---------------------------- + + integer :: fu + integer :: irec + + ! --- begin ---------------------------- + + ! write from root only .. + if ( csoc%root ) then + + ! create directory if necessary: + call CSO_CheckDir( self%filename, status ) + IF_NOT_OK_RETURN(status=1) + + ! select free file unit: + Call CSO_GetFU( fu, status ) + IF_NOT_OK_RETURN(status=1) + + ! create file: + open( unit=fu, file=trim(self%filename), iostat=status, form='formatted' ) + if ( status /= 0 ) then + write (csol,'("from file creation :")'); call csoErr + write (csol,'(" file name : ",a)') trim(self%filename); call csoErr + TRACEBACK; status=1; return + end if + + ! write header: + write (fu,'("filename;start_time;end_time;filename_state")') + ! loop over records: + do irec = 1, self%nrec + ! write records: + write (fu,'(7a)') trim(self%rec(irec)%filename), self%sep, & + trim(Pretty(self%rec(irec)%t1)), self%sep, & + trim(Pretty(self%rec(irec)%t2)), self%sep, & + trim(self%rec(irec)%filename2) + end do ! records + + ! close file: + close( unit=fu, iostat=status ) + if ( status /= 0 ) then + write (csol,'("from closing file:")'); call csoErr + write (csol,'(" ",a)') trim(self%filename); call csoErr + TRACEBACK; status=1; return + end if + + end if ! root + + ! ok + status = 0 + + end subroutine Listing_PutOut ! *** @@ -476,13 +724,132 @@ contains end subroutine Listing_Show + ! *** + + + subroutine Listing_Get( self, status, nrec, filename ) + + ! --- in/out ----------------- + + class(T_CSO_Listing), intent(in) :: self + integer, intent(out) :: status + integer, intent(out), optional :: nrec + character(len=*), intent(out), optional :: filename + + ! --- const ---------------------- + + character(len=*), parameter :: rname = mname//'/Listing_Get' + + ! --- local ------------------------ + + ! --- begin ------------------------ + + ! number of records: + if ( present(nrec) ) nrec = self%nrec + ! listing file: + if ( present(filename) ) filename = trim(self%filename) + + ! ok + status = 0 + + end subroutine Listing_Get + + + ! *** + + ! + ! Return record valus based on either: + ! irec : record index + ! iselect : selection index + ! + + subroutine Listing_GetRecord( self, status, irec, iselect, & + filename, filename2, t1, t2 ) + + ! --- in/out ----------------- + + class(T_CSO_Listing), intent(in) :: self + integer, intent(out) :: status + integer, intent(in), optional :: irec + integer, intent(in), optional :: iselect + character(len=*), intent(out), optional :: filename + character(len=*), intent(out), optional :: filename2 + type(T_CSO_DateTime), intent(out), optional :: t1, t2 + + ! --- const ---------------------- + + character(len=*), parameter :: rname = mname//'/Listing_GetRecord' + + ! --- local ------------------------ + + integer :: narg + integer :: indx + character(len=1024) :: fname + + ! --- begin ------------------------ + + ! check ... + narg = count( (/present(irec),present(iselect)/) ) + if ( narg /= 1 ) then + write (csol,'("exactly one of arguments `irec` and `iselect` sould be present, found: ",i0)') narg; call csoErr + TRACEBACK; status=1; return + end if + + ! record index: + if ( present(irec) ) then + ! check ... + if ( (irec < 1) .or. (irec > self%nrec) ) then + write (csol,'("record index ",i0," out of range 1,..,",i0)') irec, self%nrec; call csoErr + TRACEBACK; status=1; return + end if + ! store: + indx = irec + else + ! check ... + if ( (iselect < 1) .or. (iselect > self%nselect) ) then + write (csol,'("selection index ",i0," out of range 1,..,",i0)') iselect, self%nselect; call csoErr + TRACEBACK; status=1; return + end if + ! store: + indx = self%selection(iselect) + end if + + ! full path .. + if ( present(filename) ) then + ! copy: + call self%rec(indx)%Get( status, filename=fname ) + IF_NOT_OK_RETURN(status=1) + ! include path to listing file: + write (filename,'(a,"/",a)') trim(self%dirname), trim(fname) + end if + + ! full path .. + if ( present(filename2) ) then + ! copy: + call self%rec(indx)%Get( status, filename2=fname ) + IF_NOT_OK_RETURN(status=1) + ! include path to listing file: + write (filename2,'(a,"/",a)') trim(self%dirname), trim(fname) + end if + + ! timerange: + call self%rec(indx)%Get( status, t1=t1, t2=t2 ) + IF_NOT_OK_RETURN(status=1) + + ! ok + status = 0 + + end subroutine Listing_GetRecord + + ! *** ! ! Search record that assigned to time range [t1,t2] . ! ! How the record is matched depends on the timekey: - ! 'aver' : record time average in [t1,t2] + ! 'aver' : record time average in [t1,t2] + ! 'overlap' : record time range overlaps with (t1,t2] ! ! Return value: ! filename: full path to file, empty if not file found @@ -533,5 +900,69 @@ contains end subroutine Listing_SearchFile + ! *** + + ! + ! Search records that assigned to time range [t1,t2] . + ! + ! How the record is matched depends on the timekey: + ! 'aver' : record time average in [t1,t2] + ! 'overlap' : record time range overlaps with (t1,t2] + ! + ! Return value: + ! nselect : number of selected records + ! status: 0 if ok, >0 error + ! + + subroutine Listing_Select( self, t1, t2, timekey, nselect, selection, status ) + + use CSO_DateTimes, only : Pretty + + ! --- in/out ----------------- + + class(T_CSO_Listing), intent(inout) :: self + type(T_CSO_Datetime), intent(in) :: t1, t2 + character(len=*), intent(in) :: timekey + integer, intent(out) :: nselect + integer, pointer :: selection(:) + integer, intent(out) :: status + + ! --- const ---------------------- + + character(len=*), parameter :: rname = mname//'/Listing_Select' + + ! --- local ------------------------ + + integer :: irec + + ! --- begin ------------------------ + + ! init counter: + self%nselect = 0 + ! loop until first match: + do irec = 1, self%nrec + ! match? + call self%rec(irec)%Match( t1,t2, timekey, self%selected(irec), status ) + IF_NOT_OK_RETURN(status=1) + ! selected? + if ( self%selected(irec) ) then + ! increase counter: + self%nselect = self%nselect + 1 + ! store record index: + self%selection(self%nselect) = irec + end if + end do ! records + + ! return number of selected records: + nselect = self%nselect + ! pointer to selection array: + selection => self%selection + + ! ok + status = 0 + + end subroutine Listing_Select + + end module CSO_Listing -- GitLab From 0fa3f0da4f9d8488da081ab8cc7bc56fcaae5902 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 15:57:17 +0200 Subject: [PATCH 08/20] Updated logging messages. --- oper/src/cso_pixels.F90 | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/oper/src/cso_pixels.F90 b/oper/src/cso_pixels.F90 index 6c5e230..c0f2df7 100644 --- a/oper/src/cso_pixels.F90 +++ b/oper/src/cso_pixels.F90 @@ -32,6 +32,9 @@ ! 2025-10, Arjo Segers ! Support kernel convolution 'xa + A ( x - xa )'. ! +! 2026-07, Arjo Segers +! Updated logging messages. +! !############################################################################### ! #define TRACEBACK write (csol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call csoErr @@ -5744,6 +5747,8 @@ contains ! loop over data sets: do i = 1, self%n + ! info .. + write (csol,'(a,": put out `",a,"` ...")') rname, trim(self%values(i)%p%name); call csoPr ! selection? if ( associated(vars) ) then ! init flag: @@ -5756,10 +5761,20 @@ contains if ( found ) exit end do ! skip if not in list .. - if ( .not. found ) cycle + if ( .not. found ) then + ! info ... + write (csol,'(a,": not in output selection ...")') rname; call csoPr + ! skip: + cycle + end if end if ! selection - ! filter: - if ( self%values(i)%p%glb ) cycle + ! no need to put out data that is global on each processor: + if ( self%values(i)%p%glb ) then + ! info ... + write (csol,'(a,": no need to put out global data ...")') rname; call csoPr + ! skip: + cycle + end if ! put out: call self%values(i)%p%NcPutGather( ncf, mapping, add, status ) IF_NOT_OK_RETURN(status=1) -- GitLab From a851bbb3af21b19ea3e7c965fd9000992e0ff92b Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 16:08:09 +0200 Subject: [PATCH 09/20] Added objects for time series of input files. --- oper/src/cso_series.F90 | 1362 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 1362 insertions(+) create mode 100644 oper/src/cso_series.F90 diff --git a/oper/src/cso_series.F90 b/oper/src/cso_series.F90 new file mode 100644 index 0000000..48ad3a4 --- /dev/null +++ b/oper/src/cso_series.F90 @@ -0,0 +1,1362 @@ +!############################################################################### +! +! CSO_Series - collection of data and states +! +! CHANGES +! +! 2023-10 +! Do not write record to output listing if no file was written. +! +! 2024-10, Arjo Segers +! Support time stamps in output directory. +! +!############################################################################### +! +#define TRACEBACK write (csol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call csoErr +#define IF_NOT_OK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if +#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if +! +#define IF_NF90_NOT_OK_RETURN(action) if (status/=NF90_NOERR) then; csol=nf90_strerror(status); call csoErr; TRACEBACK; action; return; end if +! +#include "cso.inc" +! +!############################################################################### + +module CSO_Series + + use CSO_Logging , only : csol, csoPr, csoErr + use CSO_Listing , only : T_CSO_Listing + use CSO_Sat , only : T_CSO_Data_Set + use CSO_Sat , only : T_CSO_State_Set + use CSO_DateTimes , only : T_CSO_DateTime + + implicit none + + + ! --- in/out ---------------------------- + + private + + public :: T_CSO_Data_Series + public :: T_CSO_State_Series + + + ! --- const ------------------------------ + + character(len=*), parameter :: mname = 'CSO_Series' + + ! status parameters: + integer, parameter :: STAT_UNDEF = 0 + integer, parameter :: STAT_NEW = 1 + integer, parameter :: STAT_OPEN = 2 + integer, parameter :: STAT_FINISHED = 3 + integer, parameter :: STAT_DONE = 4 + integer, parameter :: STAT_UNKNOWN = 5 + ! names: + character(len=*), parameter :: STAT_NAME(0:5) = & + (/ 'UNDEF ', & + 'NEW ', & + 'OPEN ', & + 'FINISHED', & + 'DONE ', & + 'UNKNOWN ' /) + + + ! --- types -------------------------------- + + ! ~ pointer structure + + type P_CSO_Sat_Data + ! pointer to data: + type(T_CSO_Data_Set), pointer :: p + ! orbit filename: + character(len=1024) :: filename + character(len=1024) :: basename + ! timestamp, used for output dir: + type(T_CSO_DateTime) :: tstamp + contains + procedure :: InitCopy => CSO_Sat_Data_InitCopy + end type P_CSO_Sat_Data + + ! ~ + + type P_CSO_Sat_State + ! pointer to state: + type(T_CSO_State_Set), pointer :: p + ! part of filename: + character(len=32) :: name + end type P_CSO_Sat_State + + ! ~ series + + type T_CSO_Data_Series + ! name that might be used in output files etc: + character(len=64) :: name + ! storage for list of orbit file names: + type(T_CSO_Listing) :: listing + ! output settings: + logical :: output_enabled + character(len=1024) :: output_dir + character(len=64) :: output_subdir + ! pointers to data: + integer :: nrec + type(P_CSO_Sat_Data), allocatable :: psdata(:) ! (nrec) + ! status flag: + integer, allocatable :: stat(:) ! (nrec) + ! index list for currently open data sets: + integer :: nset + integer, allocatable :: irset(:) ! (nrec) + ! index list for previously opened data: + integer :: nprev + integer, allocatable :: irprev(:) ! (nrec) + integer, allocatable :: stat_prev(:) ! (nrec) + ! index list for newly opened data: + integer :: nnew + integer, allocatable :: irnew(:) ! (nrec) + ! index list for finished data: + integer :: nfin + integer, allocatable :: irfin(:) ! (nrec) + ! + contains + procedure :: Init => CSO_Series_Data_Init + procedure :: Done => CSO_Series_Data_Done + procedure :: Get => CSO_Series_Data_Get + procedure :: Setup => CSO_Series_Data_Setup + procedure :: Close => CSO_Series_Data_Close + procedure :: PostProc => CSO_Series_Data_PostProc + procedure :: PutOut => CSO_Series_Data_PutOut + procedure :: Finish => CSO_Series_Data_Finish + end type T_CSO_Data_Series + + !~ + + type T_CSO_State_Series + ! name, used for output files: + character(len=64) :: name + ! output dir: + character(len=1024) :: output_dir + character(len=64) :: output_subdir + ! listing of created files: + type(T_CSO_Listing) :: output_listing + ! attribute: + character(len=512) :: description + ! pointers to open states: + type(P_CSO_Sat_State), pointer :: psstate(:) ! (n) + ! + contains + procedure :: Init => CSO_Series_State_Init + procedure :: Done => CSO_Series_State_Done + procedure :: Get => CSO_Series_State_Get + procedure :: Setup => CSO_Series_State_Setup + procedure :: Close => CSO_Series_State_Close + procedure :: PostProc => CSO_Series_State_PostProc + procedure :: PutOut => CSO_Series_State_PutOut + end type T_CSO_State_Series + + + +contains + + + ! ==================================================================== + ! === + ! === pointer object + ! === + ! ==================================================================== + + + subroutine CSO_Sat_Data_InitCopy( self, source, status ) + + ! --- in/out --------------------------------- + + class(P_CSO_Sat_Data), intent(out) :: self + class(P_CSO_Sat_Data), intent(in) :: source + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Sat_Data_InitCopy' + + ! --- local ---------------------------------- + + ! --- begin ---------------------------------- + + ! assign: + self%p => source%p + ! copy: + self%filename = trim(source%filename) + self%basename = trim(source%basename) + self%tstamp = source%tstamp + + ! ok + status = 0 + + end subroutine CSO_Sat_Data_InitCopy + + + ! ==================================================================== + ! === + ! === Series of Data objects + ! === + ! ==================================================================== + + + subroutine CSO_Series_Data_Init( self, listing_filename, status, & + output_enabled, output_dir, output_subdir, name ) + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Series), intent(out) :: self + character(len=*), intent(in) :: listing_filename + integer, intent(out) :: status + logical, intent(in), optional :: output_enabled + character(len=*), intent(in), optional :: output_dir + character(len=*), intent(in), optional :: output_subdir + character(len=*), intent(in), optional :: name + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_Data_Init' + + ! --- local ---------------------------------- + + ! --- begin ---------------------------------- + + ! info ... + write (csol,'(a,": read listing file: ",a)') rname, trim(listing_filename); call csoPr + ! read listing file: + call self%listing%Init( listing_filename, status ) + IF_NOT_OK_RETURN(status=1) + !! show content: + !call self%listing%Show( status ) + !IF_NOT_OK_RETURN(status=1) + + ! number of records: + call self%listing%Get( status, nrec=self%nrec ) + IF_NOT_OK_RETURN(status=1) + + ! storage for pointers: + allocate( self%psdata(self%nrec), stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! storage, all undefined yet: + allocate( self%stat(self%nrec), source=STAT_UNDEF, stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! no lists yet: + self%nset = 0 + self%nnew = 0 + self%nfin = 0 + ! storage: + allocate( self%irset(self%nrec), source=-999, stat=status ) + IF_NOT_OK_RETURN(status=1) + allocate( self%irnew(self%nrec), source=-999, stat=status ) + IF_NOT_OK_RETURN(status=1) + allocate( self%irfin(self%nrec), source=-999, stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! storage: + allocate( self%irprev(self%nrec), source=-999, stat=status ) + IF_NOT_OK_RETURN(status=1) + allocate( self%stat_prev(self%nrec), source=-999, stat=status ) + IF_NOT_OK_RETURN(status=1) + ! no lists yet: + self%nprev = 0 + + ! by default put out .. + self%output_enabled = .true. + if ( present(output_enabled) ) self%output_enabled = output_enabled + + ! output directory, default here: + if ( present(output_dir) ) then + self%output_dir = trim(output_dir) + else + self%output_dir = '.' + end if + + ! output subdir, typically including time templates '%Y/%m' etc: + if ( present(output_subdir) ) then + self%output_subdir = trim(output_subdir) + else + self%output_subdir = '' + end if + ! add seperation character if necessary: + if ( len_trim(self%output_subdir) > 0 ) then + self%output_subdir = trim(self%output_subdir)//'/' + end if + + ! adhoc name .. + self%name = '' + if ( present(name) ) self%name = trim(name) + + ! ok + status = 0 + + end subroutine CSO_Series_Data_Init + + + ! *** + + + subroutine CSO_Series_Data_Done( self, status ) + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Series), intent(inout) :: self + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_Data_Done' + + ! --- local ---------------------------------- + + ! --- begin ---------------------------------- + + ! clear data with status 'finished': + call self%Finish( status ) + IF_NOT_OK_RETURN(status=1) + + ! clear: + deallocate( self%stat, stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! clear: + deallocate( self%irset, stat=status ) + IF_NOT_OK_RETURN(status=1) + deallocate( self%irnew, stat=status ) + IF_NOT_OK_RETURN(status=1) + deallocate( self%irfin, stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! clear: + deallocate( self%irprev, stat=status ) + IF_NOT_OK_RETURN(status=1) + deallocate( self%stat_prev, stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! done with listing: + call self%listing%Done( status ) + IF_NOT_OK_RETURN(status=1) + + ! ok + status = 0 + + end subroutine CSO_Series_Data_Done + + + ! *** + + + subroutine CSO_Series_Data_Get( self, status, & + nset, nnew, nfin, & + iset, inew, ifin, & + sdata, isnew ) + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Series), intent(inout) :: self + integer, intent(out) :: status + integer, intent(out), optional :: nset + integer, intent(out), optional :: nnew + integer, intent(out), optional :: nfin + integer, intent(in), optional :: iset + integer, intent(in), optional :: inew + integer, intent(in), optional :: ifin + type(T_CSO_Data_Set), pointer, optional :: sdata + logical, intent(out), optional :: isnew + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_Data_Get' + + ! --- local ---------------------------------- + + integer :: irec + + ! --- begin ---------------------------------- + + ! copy: + if ( present(nset) ) nset = self%nset + if ( present(nfin) ) nfin = self%nfin + if ( present(nnew) ) nnew = self%nnew + + ! check .. + if ( count((/present(iset),present(inew),present(ifin)/)) > 1 ) then + write (csol,'("provide max one of arguments: `iset`, `inew`")'); call csoErr + TRACEBACK; status=1; return + end if + ! init record index: + irec = -999 + ! switch: + !~ open set: + if ( present(iset) ) then + ! check ... + if ( (iset < 1) .or. (iset > self%nset) ) then + write (csol,'("argument `i` is ",i0," but should be in range 1,..,",i0)') iset, self%nset; call csoErr + TRACEBACK; status=1; return + end if + ! set record index: + irec = self%irset(iset) + !~ new set: + else if ( present(inew) ) then + ! check ... + if ( (inew < 1) .or. (inew > self%nnew) ) then + write (csol,'("argument `inew` is ",i0," but should be in range 1,..,",i0)') inew, self%nnew; call csoErr + TRACEBACK; status=1; return + end if + ! set record index: + irec = self%irnew(inew) + !~ finished set: + else if ( present(ifin) ) then + ! check ... + if ( (ifin < 1) .or. (ifin > self%nfin) ) then + write (csol,'("argument `ifin` is ",i0," but should be in range 1,..,",i0)') ifin, self%nfin; call csoErr + TRACEBACK; status=1; return + end if + ! set record index: + irec = self%irfin(ifin) + end if + + ! output: + if ( present(sdata) ) then + ! check .. + if ( irec < 0 ) then + write (csol,'("could not set `sdata` without one of: `iset`, `inew`, or `ifin`")'); call csoErr + TRACEBACK; status=1; return + end if + ! assign pointer: + sdata => self%psdata(irec)%p + end if + + ! output: + if ( present(isnew) ) then + ! check .. + if ( irec < 0 ) then + write (csol,'("could not set `isnew` without one of: `iset`, `inew`, or `ifin`")'); call csoErr + TRACEBACK; status=1; return + end if + ! set flag: + isnew = self%stat(irec) == STAT_NEW + end if + + ! ok + status = 0 + + end subroutine CSO_Series_Data_Get + + + ! *** + + + ! + ! Search orbit files with pixel times overlapping (t1,t2] . + ! + + subroutine CSO_Series_Data_Setup( self, rcF, rcbase, t1,t2, status ) + + use CSO_Rc , only : T_CSO_RcFile + use CSO_DateTimes, only : T_CSO_DateTime, operator(<) + use CSO_File , only : CSO_GetBasename + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Series), intent(inout) :: self + type(T_CSO_RcFile), intent(in) :: rcF + character(len=*), intent(in) :: rcbase + type(T_CSO_DateTime), intent(in) :: t1, t2 + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_Data_Setup' + + ! --- local ---------------------------------- + + type(T_CSO_DateTime) :: t1s, t2s + integer, pointer :: selection(:) + integer :: iset + integer :: irec + integer :: iprev + + ! --- begin ---------------------------------- + + ! info ... + write (csol,'(a,": setup series of data objects ...")') rname; call csoPr + + ! sorted time range: + if ( t1 < t2 ) then + t1s = t1 + t2s = t2 + else + t1s = t2 + t2s = t1 + end if + + ! clear finished data if necessary: + call self%Finish( status ) + IF_NOT_OK_RETURN(status=1) + + ! save current: + self%nprev = self%nset + self%irprev = self%irset + self%stat_prev = self%stat + ! reset status of currently opened data: + do iset = 1, self%nset + ! record index: + irec = self%irset(iset) + ! check status: + select case ( self%stat(irec) ) + !~ new or open: + case ( STAT_NEW, STAT_OPEN ) + ! reset status, might be still in use of not anymore: + self%stat(irec) = STAT_UNKNOWN + !~ done after call to Finish + case ( STAT_DONE ) + ! nothing + !~ strange .. + case ( STAT_UNDEF, STAT_UNKNOWN ) + write (csol,'("strange: record ",i0," has status `",a,"`")') irec, STAT_NAME(self%stat(irec)); call csoErr + TRACEBACK; status=1; return + !~ + case default + write (csol,'("unsupported status `",i0,"`")') self%stat(irec); call csoErr + TRACEBACK; status=1; return + end select + end do ! i + + ! reset counter: + self%nnew = 0 + self%irnew = -999 + + ! select records in listing file; + ! return number of selected records, + ! and pointer to array 'selection' with elements 1:nset the record numbers: + call self%listing%Select( t1s,t2s, 'overlap', self%nset, selection, status ) + IF_NOT_OK_RETURN(status=1) + ! info ... + write (csol,'(a,": number of selected records: ",i0)') rname, self%nset; call csoPr + ! any? + if ( self%nset > 0 ) then + + ! loop over selection: + do iset = 1, self%nset + ! record: + irec = selection(iset) + ! store record index: + self%irset(iset) = irec + + ! get orbit filename: + call self%listing%GetRecord( status, irec=irec, filename=self%psdata(irec)%filename ) + IF_NOT_OK_RETURN(status=1) + ! filename only: + call CSO_GetBasename( self%psdata(irec)%filename, self%psdata(irec)%basename, status ) + IF_NOT_OK_RETURN(status=1) + ! store start time as timestamp for output: + self%psdata(irec)%tstamp = t1s + ! info ... + !write (csol,'(a,": selected record ",i0," (",a,")")') rname, irec, trim(self%psdata(irec)%basename); call csoPr + write (csol,'(a,": selected record ",i0," (",a,")")') rname, irec, trim(self%psdata(irec)%filename); call csoPr + + ! check current status: + select case ( self%stat(irec) ) + + !~ never used before: + case ( STAT_UNDEF ) + + ! info ... + write (csol,'(a,": new ...")') rname; call csoPr + + ! reset status: + self%stat(irec) = STAT_NEW + ! store record index: + self%nnew = self%nnew + 1 + self%irnew(self%nnew) = irec + ! new object: + allocate( self%psdata(irec)%p, stat=status ) + IF_NOT_OK_RETURN(status=1) + ! inititalize orbit data, + ! read all footprints available in file: + call self%psdata(irec)%p%Init( rcF, rcbase, self%psdata(irec)%filename, status ) + IF_NOT_OK_RETURN(status=1) + + !~ already open: + case ( STAT_NEW, STAT_OPEN, STAT_UNKNOWN ) + ! info ... + write (csol,'(a,": already opened ...")') rname; call csoPr + ! reset state: + self%stat(irec) = STAT_OPEN + + !~ strange ... + case ( STAT_FINISHED ) + write (csol,'("strange: record ",i0," has status `",a,"`")') irec, STAT_NAME(self%stat(irec)); call csoErr + TRACEBACK; status=1; return + + !~ + case default + write (csol,'("unsupported status `",i0,"`")') self%stat(irec); call csoErr + TRACEBACK; status=1; return + end select + + end do ! i + + end if ! nselected > 0 + + ! reset counter: + self%nfin = 0 + self%irfin = -999 + ! search for finished records using save list: + do iprev = 1, self%nprev + ! record index: + irec = self%irprev(iprev) + ! testing ... + write (csol,'(a,": record ",i0," (",a,") previous status `",a,"`, now `",a,"`")') & + rname, irec, trim(self%psdata(irec)%basename), & + trim(STAT_NAME(self%stat_prev(irec))), trim(STAT_NAME(self%stat(irec))); call csoPr + ! check previous status: + select case ( self%stat_prev(irec) ) + !~ was open + case ( STAT_NEW, STAT_OPEN ) + ! check current status: + select case ( self%stat(irec) ) + !~ still open: + case ( STAT_OPEN ) + ! nothing to be done + !~ not open anymore .. + case ( STAT_UNKNOWN ) + ! info ... + write (csol,'(a,": finish data ",i0," (",a,")")') rname, & + irec, trim(self%psdata(irec)%basename); call csoPr + ! reset status: + self%stat(irec) = STAT_FINISHED + ! store index: + self%nfin = self%nfin + 1 + self%irfin(self%nfin) = irec + + !~ unknown .. + case default + write (csol,'("unsupported current status `",i0,"`")') self%stat(irec); call csoErr + TRACEBACK; status=1; return + end select ! stat + !~ other + case default + write (csol,'("unsupported previous status `",i0,"`")') self%stat_prev(irec); call csoErr + TRACEBACK; status=1; return + end select ! stat_prev + end do ! iprev + + ! post-process finished data: + call self%PostProc( status ) + IF_NOT_OK_RETURN(status=1) + + ! put out? + if ( self%output_enabled ) then + ! put out finished data: + call self%PutOut( status ) + IF_NOT_OK_RETURN(status=1) + end if + + ! ok + status = 0 + + end subroutine CSO_Series_Data_Setup + + + ! *** + + + ! + ! PostProcproeces currently open records, and mark all as 'finished'. + ! This should be called before "Done" to close open records properly. + ! + + subroutine CSO_Series_Data_Close( self, status ) + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Series), intent(inout) :: self + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_Data_Close' + + ! --- local ---------------------------------- + + integer :: iset + integer :: irec + + ! --- begin ---------------------------------- + + ! clear data with status 'finished': + call self%Finish( status ) + IF_NOT_OK_RETURN(status=1) + + ! mark currently open records as 'finished'; + ! reset counter: + self%nfin = 0 + self%irfin = -999 + ! loop over all currently open records: + do iset = 1, self%nset + ! record index: + irec = self%irset(iset) + ! testing ... + write (csol,'(a,": record ",i0," (",a,") status `",a,"`")') & + rname, irec, trim(self%psdata(irec)%basename), & + trim(STAT_NAME(self%stat_prev(irec))); call csoPr + ! check status: + select case ( self%stat(irec) ) + !~ is open + case ( STAT_NEW, STAT_OPEN ) + ! info ... + write (csol,'(a,": finish data ",i0," (",a,")")') rname, & + irec, trim(self%psdata(irec)%basename); call csoPr + ! reset status: + self%stat(irec) = STAT_FINISHED + ! store index: + self%nfin = self%nfin + 1 + self%irfin(self%nfin) = irec + !~ alread done: + case ( STAT_DONE ) + ! nothing to be changed + !~ other + case default + write (csol,'("unsupported status `",i0,"`")') self%stat(irec); call csoErr + TRACEBACK; status=1; return + end select ! stat + end do ! i + + ! postprocess data just marked as finished: + call self%PostProc( status ) + IF_NOT_OK_RETURN(status=1) + + ! put out? + if ( self%output_enabled ) then + ! put out finished data: + call self%PutOut( status ) + IF_NOT_OK_RETURN(status=1) + end if + + ! ok + status = 0 + + end subroutine CSO_Series_Data_Close + + + ! * + + + ! + ! postprocessing of finished data: + ! - setup exchange between domains + ! + + subroutine CSO_Series_Data_PostProc( self, status ) + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Series), intent(inout) :: self + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_Data_PostProc' + + ! --- local ---------------------------------- + + integer :: ifin + integer :: irec + + ! --- begin ---------------------------------- + + ! loop over finished data: + do ifin = 1, self%nfin + ! current: + irec = self%irfin(ifin) + ! info ... + write (csol,'(a,": post data ",i0," (",a,")")') & + rname, irec, trim(self%psdata(irec)%basename); call csoPr + + ! inquire info on pixels covering multiple domains, + ! setup exchange parameters: + call self%psdata(irec)%p%SetupExchange( status ) + IF_NOT_OK_RETURN(status=1) + + end do ! ifin + + ! ok + status = 0 + + end subroutine CSO_Series_Data_PostProc + + + ! * + + + ! + ! put out finished data + ! + + subroutine CSO_Series_Data_PutOut( self, status ) + + use CSO_File , only : CSO_SplitExt + use CSO_DateTimes, only : CSO_Format + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Series), intent(inout) :: self + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_Data_PutOut' + + ! --- local ---------------------------------- + + integer :: ifin + integer :: irec + character(len=8) :: ext + character(len=1024) :: rootname + character(len=1024) :: output_filename + + ! --- begin ---------------------------------- + + ! loop over finished data: + do ifin = 1, self%nfin + ! current: + irec = self%irfin(ifin) + + ! split in base and extension: + call CSO_SplitExt( self%psdata(irec)%basename, rootname, ext, status ) + IF_NOT_OK_RETURN(status=1) + ! target file: + write (output_filename,'(a,"/",a,a,"_data",a)') & + trim(self%output_dir), trim(self%output_subdir), & + trim(rootname), trim(ext) + ! evaluate time value templates: + call CSO_Format( output_filename, self%psdata(irec)%tstamp, status ) + IF_NOT_OK_RETURN(status=1) + ! setup gathering information, gather output arrays, write: + call self%psdata(irec)%p%PutOut( output_filename, status ) + IF_NOT_OK_RETURN(status=1) + + end do ! ifin + + ! ok + status = 0 + + end subroutine CSO_Series_Data_PutOut + + + ! * + + ! clear finished data + + subroutine CSO_Series_Data_Finish( self, status ) + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Series), intent(inout) :: self + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_Data_Finish' + + ! --- local ---------------------------------- + + integer :: ifin + integer :: irec + + ! --- begin ---------------------------------- + + ! loop over finished data: + do ifin = 1, self%nfin + ! current: + irec = self%irfin(ifin) + ! info ... + write (csol,'(a,": clear data ",i0," (",a,")")') & + rname, irec, trim(self%psdata(irec)%basename); call csoPr + ! done: + call self%psdata(irec)%p%Done( status ) + IF_NOT_OK_RETURN(status=1) + ! clear: + deallocate( self%psdata(irec)%p, stat=status ) + IF_NOT_OK_RETURN(status=1) + ! safety ... + nullify( self%psdata(irec)%p ) + ! reset state: + self%stat(irec) = STAT_DONE + ! safety ... + self%irfin(ifin) = -999 + end do ! ifin + ! reset counter: + self%nfin = 0 + + ! ok + status = 0 + + end subroutine CSO_Series_Data_Finish + + + ! ==================================================================== + ! === + ! === Series of Data objects + ! === + ! ==================================================================== + + + subroutine CSO_Series_State_Init( self, dser, name, status, & + description ) + + use CSO_File , only : CSO_GetBasename + + ! --- in/out --------------------------------- + + class(T_CSO_State_Series), intent(out) :: self + class(T_CSO_Data_Series), intent(in) :: dser + character(len=*), intent(in) :: name + integer, intent(out) :: status + character(len=*), intent(in), optional :: description + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_State_Init' + + ! --- local ---------------------------------- + + character(len=1024) :: listing_filename + character(len=1024) :: basename + character(len=1024) :: outfile + + ! --- begin ---------------------------------- + + ! store: + self%name = trim(name) + + ! output directory: + self%output_dir = trim(dser%output_dir) + self%output_subdir = trim(dser%output_subdir) + + ! description? + if ( present(description) ) then + self%description = trim(description) + else + self%description = trim(self%name) + end if + + ! storage: + allocate( self%psstate(dser%nrec), stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! put out? + if ( dser%output_enabled ) then + ! from name? + if ( len_trim(dser%name) > 0 ) then + ! target file for listing of written files, without subdir ("%Y/%m") ! + write (outfile,'(a,"/",a)') trim(self%output_dir), 'cso_listing_'//trim(dser%name)//'.csv' + else + ! name of listing file: + call dser%listing%Get( status, filename=listing_filename ) + IF_NOT_OK_RETURN(status=1) + ! basename of listing file: + call CSO_GetBasename( listing_filename, basename, status ) + IF_NOT_OK_RETURN(status=1) + ! target file for listing of written files: + write (outfile,'(a,"/",a)') trim(self%output_dir), trim(basename) + end if + ! init target listing: + call self%output_listing%Create( outfile, dser%nrec, status ) + IF_NOT_OK_RETURN(status=1) + end if ! output enabled + + ! ok + status = 0 + + end subroutine CSO_Series_State_Init + + + ! *** + + + subroutine CSO_Series_State_Done( self, dser, status ) + + ! --- in/out --------------------------------- + + class(T_CSO_State_Series), intent(inout) :: self + class(T_CSO_Data_Series), intent(in) :: dser + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_State_Done' + + ! --- local ---------------------------------- + + ! --- begin ---------------------------------- + + ! put out? + if ( dser%output_enabled ) then + ! write output listing: + call self%output_listing%PutOut( status ) + IF_NOT_OK_RETURN(status=1) + ! clear: + call self%output_listing%Done( status ) + IF_NOT_OK_RETURN(status=1) + end if + + ! clear: + deallocate( self%psstate, stat=status ) + IF_NOT_OK_RETURN(status=1) + + ! ok + status = 0 + + end subroutine CSO_Series_State_Done + + + ! *** + + + subroutine CSO_Series_State_Get( self, dser, status, iset, sstate ) + + ! --- in/out --------------------------------- + + class(T_CSO_State_Series), intent(inout) :: self + class(T_CSO_Data_Series), intent(in) :: dser + integer, intent(out) :: status + integer, intent(in), optional :: iset + type(T_CSO_State_Set), pointer, optional :: sstate + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_State_Get' + + ! --- local ---------------------------------- + + ! --- begin ---------------------------------- + + ! check ... + if ( present(iset) ) then + if ( (iset < 1) .or. (iset > dser%nset) ) then + write (csol,'("argument `iset` is ",i0," but should be in range 1,..,",i0)') iset, dser%nset; call csoErr + TRACEBACK; status=1; return + end if + end if + + ! assign pointer? + if ( present(sstate) ) then + if ( present(iset) ) then + sstate => self%psstate(dser%irset(iset))%p + else + write (csol,'("could not set `sstate` without `iset`")'); call csoErr + TRACEBACK; status=1; return + end if + end if + + ! ok + status = 0 + + end subroutine CSO_Series_State_Get + + + ! *** + + + ! + ! Setup states similar to data. + ! + ! In forward run: + ! - finish previously used states + ! - keep still used states + ! - init new states + ! + ! In reverse run: + ! - clear previously used states + ! - keep still used states + ! - read new states + ! + + subroutine CSO_Series_State_Setup( self, dser, rcF, rcbase, status, & + reverse ) + + use CSO_Rc , only : T_CSO_RcFile + use CSO_File , only : CSO_SplitExt + + ! --- in/out --------------------------------- + + class(T_CSO_State_Series), intent(inout) :: self + class(T_CSO_Data_Series), intent(in) :: dser + type(T_CSO_RcFile), intent(in) :: rcF + character(len=*), intent(in) :: rcbase + integer, intent(out) :: status + logical, intent(in), optional :: reverse + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_State_Setup' + + ! --- local ---------------------------------- + + logical :: is_reverse + integer :: ifin + integer :: irec + integer :: iset + !character(len=1024) :: filename2 + + ! --- begin ---------------------------------- + + ! info ... + write (csol,'(a,": setup series of state objects ...")') rname; call csoPr + + ! reverse run? + is_reverse = .false. + if ( present(reverse) ) is_reverse = reverse + + ! switch: + if ( is_reverse ) then + + ! loop over finished data: + do ifin = 1, dser%nfin + ! current: + irec = dser%irfin(ifin) + ! clear: + call self%psstate(irec)%p%Done( dser%psdata(irec)%p, status ) + IF_NOT_OK_RETURN(status=1) + end do ! ifin + + else + + ! postprocess states for data marked as 'finished' + ! (exchange contributions from domains, apply formula) + call self%PostProc( dser, status ) + IF_NOT_OK_RETURN(status=1) + + end if ! reverse + + ! put out? + if ( dser%output_enabled ) then + ! put out finished data: + call self%PutOut( dser, status ) + IF_NOT_OK_RETURN(status=1) + end if + + ! loop over open sets: + do iset = 1, dser%nset + ! record: + irec = dser%irset(iset) + ! check status: + select case ( dser%stat(irec) ) + !~ new: + case ( STAT_NEW ) + ! info ... + write (csol,'(a,": new state ",i0," (",a,") ",a)') & + rname, irec, trim(dser%psdata(irec)%basename), trim(self%name); call csoPr + ! allocate state: + allocate( self%psstate(irec)%p, stat=status ) + IF_NOT_OK_RETURN(status=1) + + !! get extra filename, might be used to read previously saved state: + !call dser%listing%GetRecord( status, irec=irec, filename2=filename2 ) + !IF_NOT_OK_RETURN(status=1) + + ! initialize simulation state, eventually read from specified file: + call self%psstate(irec)%p%Init( dser%psdata(irec)%p, rcF, rcbase, status, & + description=trim(self%description) ) + !filename=filename2 ) + IF_NOT_OK_RETURN(status=1) + + !~ already in use: + case ( STAT_OPEN, STAT_FINISHED ) + ! nothing to be done + !~ other + case default + write (csol,'("unsupported status `",i0,"`")') dser%stat(irec); call csoErr + TRACEBACK; status=1; return + end select ! stat + end do ! iset + + ! ok + status = 0 + + end subroutine CSO_Series_State_Setup + + + ! *** + + + ! + ! Close open states (postprocess) + ! + + subroutine CSO_Series_State_Close( self, dser, status ) + + ! --- in/out --------------------------------- + + class(T_CSO_State_Series), intent(inout) :: self + class(T_CSO_Data_Series), intent(in) :: dser + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_State_Close' + + ! --- local ---------------------------------- + + ! --- begin ---------------------------------- + + ! open data records have been marked by "CSO_Series_Data_Close" as 'finished', + ! call postprocessing to apply exchange, formula, and put out: + call self%PostProc( dser, status ) + IF_NOT_OK_RETURN(status=1) + + ! put out? + if ( dser%output_enabled ) then + ! put out: + call self%PutOut( dser, status ) + IF_NOT_OK_RETURN(status=1) + end if + + ! ok + status = 0 + + end subroutine CSO_Series_State_Close + + + ! *** + + + ! + ! PostProcprocess records that are marked as 'finished' + ! (exchange, formulas). + ! + + subroutine CSO_Series_State_PostProc( self, dser, status ) + + ! --- in/out --------------------------------- + + class(T_CSO_State_Series), intent(inout) :: self + class(T_CSO_Data_Series), intent(in) :: dser + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_State_PostProc' + + ! --- local ---------------------------------- + + integer :: irec + integer :: ifin + + ! --- begin ---------------------------------- + + ! loop over finished data: + do ifin = 1, dser%nfin + ! current: + irec = dser%irfin(ifin) + ! info ... + write (csol,'(a,": post state ",i0," (",a,") ",a)') & + rname, irec, trim(dser%psdata(irec)%basename), trim(self%name); call csoPr + + ! exchange simulations: + call self%psstate(irec)%p%Exchange( dser%psdata(irec)%p, status ) + IF_NOT_OK_RETURN(status=1) + ! apply formula: kernel convolution etc: + call self%psstate(irec)%p%ApplyFormulas( dser%psdata(irec)%p, status ) + IF_NOT_OK_RETURN(status=1) + + end do ! ifin + + ! ok + status = 0 + + end subroutine CSO_Series_State_PostProc + + + ! *** + + + ! + ! Put out records that are marked as 'finished'. + ! + + subroutine CSO_Series_State_PutOut( self, dser, status ) + + use CSO_File , only : CSO_SplitExt + use CSO_DateTimes, only : T_CSO_DateTime + use CSO_DateTimes, only : CSO_Format + + ! --- in/out --------------------------------- + + class(T_CSO_State_Series), intent(inout) :: self + class(T_CSO_Data_Series), intent(in) :: dser + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Series_State_PutOut' + + ! --- local ---------------------------------- + + integer :: irec + integer :: ifin + character(len=8) :: ext + character(len=1024) :: rootname + character(len=1024) :: output_filename + character(len=1024) :: data_filename + character(len=1024) :: state_filename + type(T_CSO_DateTime) :: t1, t2 + logical :: written + + ! --- begin ---------------------------------- + + ! loop over finished data: + do ifin = 1, dser%nfin + ! current: + irec = dser%irfin(ifin) + + ! split in base and extension: + call CSO_SplitExt( dser%psdata(irec)%basename, rootname, ext, status ) + IF_NOT_OK_RETURN(status=1) + ! target file: + write (state_filename,'(a,a,"_",a,a)') trim(self%output_subdir), trim(rootname), trim(self%name), trim(ext) + ! evaluate time value templates: + call CSO_Format( state_filename, dser%psdata(irec)%tstamp, status ) + IF_NOT_OK_RETURN(status=1) + ! full path: + write (output_filename,'(a,"/",a)') trim(self%output_dir), trim(state_filename) + ! setup gathering information, gather output arrays, write; + ! return flag if file is actually written, might be skipped becasuse no data: + call self%psstate(irec)%p%PutOut( dser%psdata(irec)%p, output_filename, status, & + written=written ) + IF_NOT_OK_RETURN(status=1) + + ! clear: + call self%psstate(irec)%p%Done( dser%psdata(irec)%p, status ) + IF_NOT_OK_RETURN(status=1) + + ! file written? then add record to listing: + if ( written ) then + ! name of corresponding data file: + write (data_filename,'(a,a,"_",a,a)') trim(self%output_subdir), trim(rootname), 'data', trim(ext) + ! evaluate time value templates: + call CSO_Format( data_filename, dser%psdata(irec)%tstamp, status ) + IF_NOT_OK_RETURN(status=1) + ! time range from listing: + call dser%listing%GetRecord( status, irec=irec, t1=t1, t2=t2 ) + IF_NOT_OK_RETURN(status=1) + ! add record to output listing: + call self%output_listing%AddRecord( data_filename, t1, t2, status, & + filename2=state_filename ) + IF_NOT_OK_RETURN(status=1) + end if + + end do ! ifin + + ! ok + status = 0 + + end subroutine CSO_Series_State_PutOut + + +end module CSO_Series -- GitLab From a8fc1b8b33bb26b3e9eee1307bd4aa7d9a60654b Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 15:58:48 +0200 Subject: [PATCH 10/20] Introduced `Series` object to open multiple files from a time series for a single interval. Pixel selection based on time values. --- oper/src/Makefile_deps | 43 +- oper/src/cso.F90 | 5 + oper/src/cso_sat.F90 | 707 ++++++++++++--------------------- oper/src/tutorial_oper_S5p.F90 | 600 ++++++++++++++-------------- 4 files changed, 585 insertions(+), 770 deletions(-) diff --git a/oper/src/Makefile_deps b/oper/src/Makefile_deps index a9554a5..9c24f8a 100644 --- a/oper/src/Makefile_deps +++ b/oper/src/Makefile_deps @@ -1,21 +1,22 @@ -cso.o : cso.F90 cso.inc cso_ncfile.o cso_profile.o cso_grid.o cso_tools.o cso_sat.o cso_listing.o cso_rc.o cso_string.o cso_datetime.o cso_domains.o cso_logging.o cso_comm.o -cso_comm.o : cso_comm.F90 cso.inc cso_logging.o -cso_datetime.o : cso_datetime.F90 cso.inc cso_string.o cso_logging.o -cso_domains.o : cso_domains.F90 cso.inc cso_comm.o cso_logging.o -cso_exchange.o : cso_exchange.F90 cso.inc cso_logging.o -cso_file.o : cso_file.F90 cso.inc cso_logging.o -cso_grid.o : cso_grid.F90 cso.inc cso_ncfile.o cso_tools.o cso_parray.o cso_logging.o -cso_listing.o : cso_listing.F90 cso.inc cso_string.o cso_file.o cso_datetime.o cso_logging.o -cso_logging.o : cso_logging.F90 cso.inc -cso_mapping.o : cso_mapping.F90 cso.inc cso_parray.o cso_swapping.o cso_comm.o cso_logging.o -cso_ncfile.o : cso_ncfile.F90 cso.inc cso_string.o cso_file.o cso_comm.o cso_logging.o -cso_parray.o : cso_parray.F90 cso.inc cso_logging.o -cso_pixels.o : cso_pixels.F90 cso.inc cso_profile.o cso_exchange.o cso_string.o cso_swapping.o cso_parray.o cso_domains.o cso_comm.o cso_ncfile.o cso_logging.o -cso_profile.o : cso_profile.F90 cso.inc cso_logging.o -cso_rc.o : cso_rc.F90 cso.inc cso_datetime.o cso_file.o cso_string.o cso_logging.o -cso_sat.o : cso_sat.F90 cso.inc cso_string.o cso_rc.o cso_comm.o cso_domains.o cso_swapping.o cso_exchange.o cso_mapping.o cso_pixels.o cso_ncfile.o cso_logging.o -cso_string.o : cso_string.F90 cso.inc cso_logging.o -cso_swapping.o : cso_swapping.F90 cso.inc cso_parray.o cso_domains.o cso_comm.o cso_logging.o -cso_tools.o : cso_tools.F90 cso.inc cso_logging.o -tutorial_oper_S5p.o : tutorial_oper_S5p.F90 cso.inc cso.o -tutorial_oper_adj-test.o : tutorial_oper_adj-test.F90 cso.inc cso.o +cso_comm.o : cso_comm.F90 cso.inc cso_logging.o +cso_datetime.o : cso_datetime.F90 cso.inc cso_string.o cso_logging.o +cso_domains.o : cso_domains.F90 cso.inc cso_comm.o cso_logging.o +cso_exchange.o : cso_exchange.F90 cso.inc cso_logging.o +cso.o : cso.F90 cso.inc cso_ncfile.o cso_profile.o cso_grid.o cso_tools.o cso_series.o cso_sat.o cso_listing.o cso_rc.o cso_string.o cso_datetime.o cso_domains.o cso_logging.o cso_comm.o +cso_file.o : cso_file.F90 cso.inc cso_logging.o +cso_grid.o : cso_grid.F90 cso.inc cso_ncfile.o cso_tools.o cso_parray.o cso_logging.o +cso_listing.o : cso_listing.F90 cso.inc cso_comm.o cso_string.o cso_file.o cso_datetime.o cso_logging.o +cso_logging.o : cso_logging.F90 cso.inc +cso_mapping.o : cso_mapping.F90 cso.inc cso_parray.o cso_swapping.o cso_comm.o cso_logging.o +cso_ncfile.o : cso_ncfile.F90 cso.inc cso_string.o cso_file.o cso_comm.o cso_logging.o +cso_parray.o : cso_parray.F90 cso.inc cso_logging.o +cso_pixels.o : cso_pixels.F90 cso.inc cso_tools.o cso_profile.o cso_exchange.o cso_string.o cso_swapping.o cso_parray.o cso_domains.o cso_comm.o cso_ncfile.o cso_logging.o +cso_profile.o : cso_profile.F90 cso.inc cso_logging.o +cso_rc.o : cso_rc.F90 cso.inc cso_datetime.o cso_file.o cso_string.o cso_logging.o +cso_sat.o : cso_sat.F90 cso.inc cso_file.o cso_datetime.o cso_string.o cso_rc.o cso_comm.o cso_domains.o cso_exchange.o cso_mapping.o cso_pixels.o cso_ncfile.o cso_logging.o +cso_series.o : cso_series.F90 cso.inc cso_file.o cso_rc.o cso_datetime.o cso_sat.o cso_listing.o cso_logging.o +cso_string.o : cso_string.F90 cso.inc cso_logging.o +cso_swapping.o : cso_swapping.F90 cso.inc cso_parray.o cso_domains.o cso_comm.o cso_logging.o +cso_tools.o : cso_tools.F90 cso.inc cso_logging.o +tutorial_oper_adj-test.o : tutorial_oper_adj-test.F90 cso.inc cso.o +tutorial_oper_S5p.o : tutorial_oper_S5p.F90 cso.inc cso.o diff --git a/oper/src/cso.F90 b/oper/src/cso.F90 index 046cf2b..ed566cd 100644 --- a/oper/src/cso.F90 +++ b/oper/src/cso.F90 @@ -7,6 +7,10 @@ ! ! For usage, see "tutorial.F90". ! +! CHANGES +! +! 2026-07, Arjo Segers +! Include "CSO_Series" module. ! !### macro's ########################################################### ! @@ -28,6 +32,7 @@ module CSO use CSO_Rc use CSO_Listing use CSO_Sat + use CSO_Series use CSO_Tools use CSO_Grid use CSO_Profile diff --git a/oper/src/cso_sat.F90 b/oper/src/cso_sat.F90 index 1544399..50a9c75 100644 --- a/oper/src/cso_sat.F90 +++ b/oper/src/cso_sat.F90 @@ -3,24 +3,30 @@ ! CSO_Sat - simulations of satellite retrievals ! ! -! HISTORY +! CHANGES ! -! 2022-09, Arjo Segers -! Support input and output of packed variables. +! 2022-09, Arjo Segers +! Support input and output of packed variables. ! -! 2023-01, Arjo Segers -! Support files with "sample" coordinate instead of "pixel", -! and without "corner" data (also implies no "area"). +! 2023-01, Arjo Segers +! Support files with "sample" coordinate instead of "pixel", +! and without "corner" data (also implies no "area"). ! -! 2023-01, Arjo Segers -! Added option to define datatype of variable; -! mainly used to define dtype "i1" (integer(1)) or "c" (character) -! that are used for ground observation data. +! 2023-01, Arjo Segers +! Added option to define datatype of variable; +! mainly used to define dtype "i1" (integer(1)) or "c" (character) +! that are used for ground observation data. ! -! 2023-01, Arjo Segers -! Return null pointers for corner locations if these are not available. +! 2023-01, Arjo Segers +! Return null pointers for corner locations if these are not available. ! +! 2024-03, Arjo Segers +! Removed swap routines. ! +! 2026-07, Arjo Segers +! Renamed `T_CSO_Sat_Data` and `T_CSO_Sat_State` to `T_CSO_Data_Set` and `T_CSO_State_Set` +! to facilitate `Series` object. +! !############################################################################### ! #define TRACEBACK write (csol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call csoErr @@ -43,7 +49,6 @@ module CSO_Sat use CSO_Pixels , only : T_Track_0D, T_Track_1D use CSO_Mapping , only : T_Mapping use CSO_Exchange, only : T_Exchange - use CSO_Swapping, only : T_Swapping use CSO_Domains , only : T_CSO_Selection1D implicit none @@ -53,8 +58,8 @@ module CSO_Sat private - public :: T_CSO_Sat_Data - public :: T_CSO_Sat_State + public :: T_CSO_Data_Set + public :: T_CSO_State_Set ! --- const ------------------------------ @@ -101,7 +106,7 @@ module CSO_Sat ! timestamps ! kernels ! - type T_CSO_Sat_Data + type T_CSO_Data_Set ! base for rcfile keys: character(len=32) :: rcbase ! filename: @@ -112,6 +117,8 @@ module CSO_Sat ! how to read? logical :: read_on_root logical :: read_by_me + ! has data been read? + logical :: local_data_is_read ! !~ netcdf id's type(T_NcFile) :: ncf @@ -150,6 +157,11 @@ module CSO_Sat ! ! number of pixels loaded: integer :: npix + ! + ! local time selection: + integer :: ntsel + integer, pointer :: ipix__tsel(:) + ! ! index in global arrays: integer, allocatable :: iglb(:) ! (npix) @@ -173,32 +185,27 @@ module CSO_Sat ! flag: logical :: exchange_is_setup ! - ! info on how to swap pixels to other domain decomposition: - type(T_Swapping), pointer :: swp - ! flag: - logical :: swap_is_setup - ! ! pixel data: type(T_PixelDatas) :: pd ! contains - procedure :: Init => CSO_Sat_Data_Init - procedure :: InitSwap => CSO_Sat_Data_InitSwap - procedure :: Done => CSO_Sat_Data_Done - procedure :: Get => CSO_Sat_Data_Get - procedure :: GetData => CSO_Sat_Data_GetData - procedure :: Read => CSO_Sat_Data_Read - procedure :: GetPixel => CSO_Sat_Data_GetPixel - procedure :: SetPixel => CSO_Sat_Data_SetPixel - procedure :: SetMapping => CSO_Sat_Data_SetMapping - procedure :: SetupExchange => CSO_Sat_Data_SetupExchange - procedure :: PutOut => CSO_Sat_Data_PutOut - end type T_CSO_Sat_Data + procedure :: Init => CSO_Sat_Data_Init + procedure :: Done => CSO_Sat_Data_Done + procedure :: Get => CSO_Sat_Data_Get + procedure :: GetData => CSO_Sat_Data_GetData + procedure :: GetTimeSelection => CSO_Sat_Data_GetTimeSelection + procedure :: Read => CSO_Sat_Data_Read + procedure :: GetPixel => CSO_Sat_Data_GetPixel + procedure :: SetPixel => CSO_Sat_Data_SetPixel + procedure :: SetMapping => CSO_Sat_Data_SetMapping + procedure :: SetupExchange => CSO_Sat_Data_SetupExchange + procedure :: PutOut => CSO_Sat_Data_PutOut + end type T_CSO_Data_Set ! ! state simulations ! - type T_CSO_Sat_State + type T_CSO_State_Set ! annote: character(len=256) :: description ! copied from data: @@ -215,7 +222,6 @@ module CSO_Sat ! contains procedure :: Init => CSO_Sat_State_Init - procedure :: InitSwap => CSO_Sat_State_InitSwap procedure :: Done => CSO_Sat_State_Done procedure :: GetData => CSO_Sat_State_GetData procedure :: GatherData => CSO_Sat_State_GatherData @@ -230,7 +236,7 @@ module CSO_Sat procedure :: GetForcing => CSO_Sat_State_GetForcing procedure :: PutOut => CSO_Sat_State_PutOut procedure :: ReadForcing => CSO_Sat_State_ReadForcing - end type T_CSO_Sat_State + end type T_CSO_State_Set @@ -592,7 +598,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(out) :: self + class(T_CSO_Data_Set), intent(out) :: self type(T_CSO_RcFile), intent(in) :: rcF character(len=*), intent(in) :: rcbase character(len=*), intent(in) :: filename @@ -805,158 +811,21 @@ contains call self%exch(pid)%Init( pid, status ) IF_NOT_OK_RETURN(status=1) end do ! pid - ! not setup yet: - self%exchange_is_setup = .false. end if ! npes > 0 + ! not setup yet: + self%exchange_is_setup = .false. - ! no swapping yet: - self%swap_is_setup = .false. - - ! ok - status = 0 - - end subroutine CSO_Sat_Data_Init - - - ! *** - - - subroutine CSO_Sat_Data_InitSwap( self, sdata, doms, doms_f, status ) - - use CSO_Comm , only : csoc - use CSO_Domains , only : T_CSO_Domains - use CSO_Swapping, only : T_Swapping - - ! --- in/out --------------------------------- - - class(T_CSO_Sat_Data), intent(out) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata - type(T_CSO_Domains), intent(in) :: doms - type(T_CSO_Domains), intent(in) :: doms_f - integer, intent(out) :: status - - ! --- const ---------------------------------- - - character(len=*), parameter :: rname = mname//'/CSO_Sat_Data_InitSwap' - - ! --- local ---------------------------------- - - type(T_Swapping) :: mswp - integer :: pid - - ! --- begin ---------------------------------- - - ! info .. - write (csol,'(a,": swap satellite data ...")') rname; call csoPr - - ! copy settings info: - self%rcbase = sdata%rcbase - self%filename = sdata%filename - - ! copy flags: - self%with_track = sdata%with_track - - ! copy dimensions: - self%nglb = sdata%nglb - self%ncorner = sdata%ncorner - self%nlayer = sdata%nlayer - self%nretr = sdata%nretr - - ! copy track info (needed on root only): - if ( csoc%root .and. (self%nglb > 0) ) then - ! copy dims: - self%ntx = sdata%ntx - self%nty = sdata%nty - ! copy tracks: - call self%glb_track_lon%InitCopy( sdata%glb_track_lon, status ) - IF_NOT_OK_RETURN(status=1) - call self%glb_track_lat%InitCopy( sdata%glb_track_lat, status ) - IF_NOT_OK_RETURN(status=1) - call self%glb_track_clons%InitCopy( sdata%glb_track_clons, status ) - IF_NOT_OK_RETURN(status=1) - call self%glb_track_clats%InitCopy( sdata%glb_track_clats, status ) - IF_NOT_OK_RETURN(status=1) - end if - - ! copy flags: - self%putout_mapping = sdata%putout_mapping - - ! ~ swap pixel array - - ! any pixels on some domain? - if ( sdata%slc%nsel_tot > 0 ) then - - ! reset flag: - self%swap_is_setup = .true. - ! storage: - allocate( self%swp, stat=status ) - IF_NOT_OK_RETURN(status=1) - ! info on swapping from current decomposition 'doms' to decomposition 'doms_f'; - ! use the mapping info to know with which of the new domains a pixel overplaps ; - ! define for swapping of (npix) arrays: - call self%swp%Init( 'pix', sdata%mapping%map_n, sdata%mapping%map_ii, sdata%mapping%map_jj, & - doms, doms_f, status ) - IF_NOT_OK_RETURN(status=1) - ! idem for (nmap) arrays: - call mswp%Init( 'map', sdata%mapping%map_n, sdata%mapping%map_ii, sdata%mapping%map_jj, & - doms, doms_f, status ) - IF_NOT_OK_RETURN(status=1) - - ! new number of local pixels: - self%npix = self%swp%nrecv - - ! storage per local pixel for global index: - allocate( self%iglb(max(1,self%npix)), stat=status ) - IF_NOT_OK_RETURN(status=1) - ! swap global indices: - call self%swp%Swap( sdata%iglb, self%iglb, status ) - IF_NOT_OK_RETURN(status=1) - - else - ! no pixels: - self%npix = 0 - end if - - ! store selection info, collect all info on root for: - ! - scatter after reading from root (optional) - ! - gather before writing from root (always needed): - call self%slc%Init( self%nglb, self%npix, self%iglb, .true., status ) - IF_NOT_OK_RETURN(status=1) - - ! any pixels on some domain? - if ( sdata%slc%nsel_tot > 0 ) then - - ! swap pixel data to target decomposition: - call self%pd%InitSwap( sdata%pd, self%swp, status ) - IF_NOT_OK_RETURN(status=1) - - ! swap mapping info, requires swapping of both pixel and mapping arrays: - call self%mapping%InitSwap( sdata%mapping, self%swp, mswp, status ) - IF_NOT_OK_RETURN(status=1) - - ! clear: - call mswp%Done( status ) - IF_NOT_OK_RETURN(status=1) - - end if ! any pixels on some domain + ! no time selection yet: + self%ntsel = 0 + nullify( self%ipix__tsel ) - ! exchange info needed? - if ( csoc%npes > 1 ) then - ! storage: - allocate( self%exch(0:csoc%npes-1), stat=status ) - IF_NOT_OK_RETURN(status=1) - ! loop: - do pid = 0, csoc%npes-1 - ! init empty: - call self%exch(pid)%Init( pid, status ) - IF_NOT_OK_RETURN(status=1) - end do ! pid - end if ! npes > 0 + ! set flag: + self%local_data_is_read = .false. ! ok status = 0 - end subroutine CSO_Sat_Data_InitSwap + end subroutine CSO_Sat_Data_Init ! *** @@ -968,7 +837,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(inout) :: self + class(T_CSO_Data_Set), intent(inout) :: self integer, intent(out) :: status ! --- const ---------------------------------- @@ -984,46 +853,56 @@ contains !! testing .. !write (csol,'(a,": done with sat data")') rname; call csoPr - ! * done with arrays allocated by Read - - ! local pixels? - if ( self%npix > 0 ) then + ! any pixels in global domain? + if ( self%nglb > 0 ) then - deallocate( self%loc_lons, stat=status ) - IF_NOT_OK_RETURN(status=1) - deallocate( self%loc_lats, stat=status ) + ! storage for selection flags: + deallocate( self%glb_select, stat=status ) IF_NOT_OK_RETURN(status=1) - ! corners defined? - if ( self%ncorner > 0 ) then - deallocate( self%loc_clons, stat=status ) + + ! track defined? + if ( self%with_track .and. csoc%root ) then + call self%glb_track_lon%Done( status ) IF_NOT_OK_RETURN(status=1) - deallocate( self%loc_clats, stat=status ) + call self%glb_track_lat%Done( status ) IF_NOT_OK_RETURN(status=1) - end if + call self%glb_track_clons%Done( status ) + IF_NOT_OK_RETURN(status=1) + call self%glb_track_clats%Done( status ) + IF_NOT_OK_RETURN(status=1) + end if ! track - end if ! npix > 0 + ! variables allocated by "Read": + if ( self%local_data_is_read ) then - ! output collection: - call self%slc%Done( status ) - IF_NOT_OK_RETURN(status=1) + ! global index: + deallocate( self%iglb, stat=status ) + IF_NOT_OK_RETURN(status=1) - ! global index: - deallocate( self%iglb, stat=status ) - IF_NOT_OK_RETURN(status=1) + ! output collection: + call self%slc%Done( status ) + IF_NOT_OK_RETURN(status=1) - ! * done with arrays allocated by Init + end if ! local data read - ! swap is used? - if ( self%swap_is_setup ) then - ! clear: - call self%swp%Done( status ) - IF_NOT_OK_RETURN(status=1) - ! clear: - deallocate( self%swp, stat=status ) - IF_NOT_OK_RETURN(status=1) - ! reset flag: - self%swap_is_setup = .false. - end if ! swap is used + ! local pixels? + if ( self%npix > 0 ) then + deallocate( self%loc_lons, stat=status ) + IF_NOT_OK_RETURN(status=1) + deallocate( self%loc_lats, stat=status ) + IF_NOT_OK_RETURN(status=1) + ! corners defined? + if ( self%ncorner > 0 ) then + deallocate( self%loc_clons, stat=status ) + IF_NOT_OK_RETURN(status=1) + deallocate( self%loc_clats, stat=status ) + IF_NOT_OK_RETURN(status=1) + end if + end if + + end if ! nglb > 0 + + ! * done with arrays allocated by Init ! exchange info needed? if ( csoc%npes > 1 ) then @@ -1042,22 +921,6 @@ contains call self%mapping%Done( status ) IF_NOT_OK_RETURN(status=1) - ! track defined? - if ( self%with_track .and. csoc%root ) then - call self%glb_track_lon%Done( status ) - IF_NOT_OK_RETURN(status=1) - call self%glb_track_lat%Done( status ) - IF_NOT_OK_RETURN(status=1) - call self%glb_track_clons%Done( status ) - IF_NOT_OK_RETURN(status=1) - call self%glb_track_clats%Done( status ) - IF_NOT_OK_RETURN(status=1) - end if ! track - - ! storage for selection flags: - deallocate( self%glb_select, stat=status ) - IF_NOT_OK_RETURN(status=1) - ! done with pixel data: call self%pd%Done( status ) IF_NOT_OK_RETURN(status=1) @@ -1082,7 +945,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(inout) :: self + class(T_CSO_Data_Set), intent(inout) :: self type(T_CSO_RcFile), intent(in) :: rcF character(len=*), intent(in) :: rcbase integer, intent(out) :: status @@ -1193,6 +1056,13 @@ contains end do ! dvars + ! also add time values: + ! - assume that 'time' is available in the file + ! - no extra dimensions after ('pixel',) + call self%pd%NcInit( 'time', '', self%ncf%ncid, 'time', status, & + slc=self%slc, xtype='real' ) + IF_NOT_OK_RETURN(status=1) + ! file opened? if ( self%read_by_me ) then ! close: @@ -1300,7 +1170,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(in) :: self + class(T_CSO_Data_Set), intent(in) :: self integer, intent(out) :: status integer, intent(out), optional :: nglb @@ -1462,7 +1332,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(in) :: self + class(T_CSO_Data_Set), intent(in) :: self integer, intent(out) :: status integer, intent(in), optional :: id @@ -1496,6 +1366,109 @@ contains ! *** + ! + ! Search pixels for which 'time' data is in (t1,t2] . + ! In case of and ajdjoint run, t2 might be before t1, and that case (t2,t1] is used. + ! Return values: + ! - ntsel : number of selected pixels + ! - ipix__tsel : pointer to integer array of length at least (ntsel) with pixel indices + ! + + subroutine CSO_Sat_Data_GetTimeSelection( self, t1,t2, ntsel, ipix__tsel, status ) + + use CSO_DateTimes, only : T_CSO_DateTime, operator(<), wrtcsol + use CSO_DateTimes, only : CompressTime + + ! --- in/out --------------------------------- + + class(T_CSO_Data_Set), intent(inout) :: self + type(T_CSO_DateTime), intent(in) :: t1, t2 + integer, intent(out) :: ntsel + integer, pointer :: ipix__tsel(:) + integer, intent(out) :: status + + ! --- const ---------------------------------- + + character(len=*), parameter :: rname = mname//'/CSO_Sat_Data_GetTimeSelection' + + ! --- local ---------------------------------- + + real, pointer :: tvalues(:) ! (npix) + character(len=64) :: tunits + real :: tv1, tv2 + integer :: ipix + integer :: itsel + + ! --- begin ---------------------------------- + + !! testing .. + !call wrtcsol( rname//': time range ', (/t1,t2/) ); call csoPr + + ! any local pixels? + if ( self%npix > 0 ) then + + ! get pointer to real values and units: + call self%pd%GetData( status, name='time', & + data0=tvalues, units=tunits ) + IF_NOT_OK_RETURN(status=1) + + ! convert time interval bounds to values: + if ( t1 < t2 ) then + ! "forward" with t1 self%ipix__tsel + + else + + ! no selection .. + ntsel = 0 + nullify( ipix__tsel ) + + end if ! npix > 0 + + ! ok + status = 0 + + end subroutine CSO_Sat_Data_GetTimeSelection + + + ! *** + + subroutine CSO_Sat_Data_GetPixel( self, ipix, status, & glbid, lon, lat, clons, clats, & name, value, profile, & @@ -1503,7 +1476,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(in) :: self + class(T_CSO_Data_Set), intent(in) :: self integer, intent(in) :: ipix integer, intent(out) :: status @@ -1627,7 +1600,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(inout) :: self + class(T_CSO_Data_Set), intent(inout) :: self integer, intent(in) :: ipix integer, intent(out) :: status @@ -1695,7 +1668,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(inout) :: self + class(T_CSO_Data_Set), intent(inout) :: self real, intent(in) :: area(:) ! (npix) integer, intent(in) :: nw(:) ! (npix) integer, intent(in) :: ii(:) ! (nmap) @@ -1750,7 +1723,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(inout) :: self + class(T_CSO_Data_Set), intent(inout) :: self character(len=*), intent(in) :: filename integer, intent(out) :: status @@ -1942,7 +1915,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_Data), intent(inout) :: self + class(T_CSO_Data_Set), intent(inout) :: self integer, intent(out) :: status ! --- const ---------------------------------- @@ -2243,7 +2216,7 @@ contains ! ! ! --- in/out --------------------------------- ! -! class(T_CSO_Sat_Data), intent(inout) :: self +! class(T_CSO_Data_Set), intent(inout) :: self ! type(T_Random), intent(inout) :: rnd ! real, intent(out) :: v(:) ! (npix) ! integer, intent(out) :: status @@ -2303,20 +2276,21 @@ contains subroutine CSO_Sat_State_Init( self, sdata, rcF, rcbase, status, & - description ) + description, filename ) use CSO_Rc , only : T_CSO_RcFile use CSO_String , only : CSO_ReadFromLine ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(out) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(out) :: self + type(T_CSO_Data_Set), intent(in) :: sdata type(T_CSO_RcFile), intent(in) :: rcF character(len=*), intent(in) :: rcbase integer, intent(out) :: status character(len=*), intent(in), optional :: description + character(len=*), intent(in), optional :: filename ! --- const ---------------------------------- @@ -2469,204 +2443,12 @@ contains ! *** - ! - ! Init state (defined on domain "doms_f") as swapped version of "sstate" (defined on "doms"), - ! use the info in data "sdata_f" (already defined on "doms_f") - ! - - subroutine CSO_Sat_State_InitSwap( self, sdata_f, sstate, status ) - - use CSO_Comm , only : csoc -! use CSO_Domains , only : T_CSO_Domains -! use CSO_Swapping, only : T_Swapping - -! use CSO_String , only : CSO_ReadFromLine - - ! --- in/out --------------------------------- - - class(T_CSO_Sat_State), intent(out) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata_f - type(T_CSO_Sat_State), intent(in) :: sstate -! type(T_CSO_Domains), intent(in) :: doms -! type(T_CSO_Domains), intent(in) :: doms_f - integer, intent(out) :: status - - ! --- const ---------------------------------- - - character(len=*), parameter :: rname = mname//'/CSO_Sat_State_InitSwap' - - ! --- local ---------------------------------- - -! type(T_Swapping) :: swp -! character(len=1024) :: line -! integer :: ivar -! integer :: iudim -! integer :: mxdim -! integer :: i -! character(len=32) :: vname -! character(len=256) :: xnames -! character(len=256) :: aname -! character(len=1024) :: avalue -! character(len=256) :: units - - ! --- begin ---------------------------------- - - ! copy: - self%description = trim(sstate%description) - - ! copy: - self%npix = sdata_f%npix - - ! any pixels in global domain? - if ( sdata_f%nglb > 0 ) then - - ! info ... - write (csol,'(a,": define state variables ...")') rname; call csoPr - -! ! init extra output: -! call self%pd%Init( status ) -! IF_NOT_OK_RETURN(status=1) -! -! ! line with variable names: -! call rcF%Get( trim(rcbase)//'.vars', line, status ) -! IF_NOT_OK_RETURN(status=1) -! ! loop over elements: -! do -! ! empty? -! if ( len_trim(line) == 0 ) exit -! ! next part: -! call CSO_ReadFromLine( line, vname, status, sep=' ' ) -! IF_NOT_OK_RETURN(status=1) -! -! ! dimension names: -! call rcF%Get( trim(rcbase)//'.var.'//trim(vname)//'.dims', xnames, status ) -! IF_NOT_OK_RETURN(status=1) -! ! info ... -! write (csol,'(a,": variable `",a,"` ...")') rname, trim(vname); call csoPr -! !! testing ... -! !write (csol,'(a,": dimensions : ",a)') rname, trim(xnames); call csoPr -! -! ! define new variable, return index: -! call self%pd%Def( vname, xnames, ivar, status ) -! IF_NOT_OK_RETURN(status=1) -! -! ! attribute names: -! call rcF%Get( trim(rcbase)//'.var.'//trim(vname)//'.attrs', xnames, status ) -! IF_NOT_OK_RETURN(status=1) -! ! loop: -! do -! ! empty? -! if ( len_trim(xnames) == 0 ) exit -! ! next name: -! call CSO_ReadFromLine( xnames, aname, status, sep=' ' ) -! IF_NOT_OK_RETURN(status=1) -! ! attribute value: -! call rcF%Get( trim(rcbase)//'.var.'//trim(vname)//'.attr.'//trim(aname), avalue, status ) -! IF_NOT_OK_RETURN(status=1) -! ! store: -! call self%pd%SetAttr( ivar, aname, avalue, status ) -! IF_NOT_OK_RETURN(status=1) -! end do ! attributes -! -! ! formula defined? -! call rcF%Get( trim(rcbase)//'.var.'//trim(vname)//'.formula', aname, status, default='' ) -! IF_ERROR_RETURN(status=1) -! ! defined? -! if ( len_trim(aname) > 0 ) then -! ! read formula terms: -! call rcF%Get( trim(rcbase)//'.var.'//trim(vname)//'.formula_terms', avalue, status ) -! IF_NOT_OK_RETURN(status=1) -! !! info ... -! !write (csol,'(a,": formula : ",a)') rname, trim(aname); call csoPr -! !write (csol,'(a,": terms : ",a)') rname, trim(avalue); call csoPr -! ! store: -! call self%pd%SetFormula( ivar, aname, avalue, status ) -! IF_NOT_OK_RETURN(status=1) -! end if ! formula -! -! end do ! elements - - ! any pixels on some domain? - if ( sdata_f%slc%nsel_tot > 0 ) then - - ! check .. - if ( .not. sdata_f%swap_is_setup ) then - write (csol,'("swap of data is not setup yet; sdata structure not swapped yet?")'); call csoErr - TRACEBACK; status=1; return - end if - - ! swap pixel data to target decomposition, - ! use the swapping info that was setup in "sdata_f": - call self%pd%InitSwap( sstate%pd, sdata_f%swp, status ) - IF_NOT_OK_RETURN(status=1) - - end if ! any pixels on some domain - -! ! extract counter for user defined variables: -! call self%pd%Get( status, n=self%nvar ) -! IF_NOT_OK_RETURN(status=1) -! -! ! extract total number of dimensions: -! call self%pd%Get( status, ndim=mxdim ) -! IF_NOT_OK_RETURN(status=1) -! ! maximum storage for user defined dim names: -! allocate( self%udimnames(mxdim), stat=status ) -! IF_NOT_OK_RETURN(status=1) -! ! init counter: -! iudim = 0 -! ! loop over all dims: -! do i = 1, mxdim -! ! dimension name: -! call self%pd%GetDim( i, status, name=vname ) -! IF_NOT_OK_RETURN(status=1) -! ! default or user defined? -! select case ( trim(vname) ) -! !~ number of apriori layers: -! case ( 'layer' ) -! ! set length: -! call self%pd%SetDim( vname, sdata%nlayer, status ) -! IF_NOT_OK_RETURN(status=1) -! !~ number of retrieval layers: -! case ( 'retr' ) -! ! set length: -! call self%pd%SetDim( vname, sdata%nretr, status ) -! IF_NOT_OK_RETURN(status=1) -! !~ user defined: -! case default -! ! increase counter: -! iudim = iudim + 1 -! ! store: -! self%udimnames(iudim) = trim(vname) -! end select -! end do ! dims -! ! store number of user defined dims: -! self%nudim = iudim -! -! ! output selections: -! call self%outkeys%Init( rcF, rcbase, status ) -! IF_NOT_OK_RETURN(status=1) - - ! output selections: - call self%outkeys%InitCopy( sstate%outkeys, status ) - IF_NOT_OK_RETURN(status=1) - - end if ! any global pixels - - ! ok - status = 0 - - end subroutine CSO_Sat_State_InitSwap - - - ! *** - - subroutine CSO_Sat_State_Done( self, sdata, status ) ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(inout) :: self + type(T_CSO_Data_Set), intent(in) :: sdata integer, intent(out) :: status ! --- const ---------------------------------- @@ -2714,7 +2496,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(in) :: self + class(T_CSO_State_Set), intent(in) :: self integer, intent(out) :: status integer, intent(out), optional :: nudim @@ -2763,7 +2545,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(in) :: self + class(T_CSO_State_Set), intent(in) :: self integer, intent(in) :: id integer, intent(out) :: status @@ -2799,7 +2581,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self + class(T_CSO_State_Set), intent(inout) :: self character(len=*), intent(in) :: name integer, intent(in) :: length integer, intent(out) :: status @@ -2829,7 +2611,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self + class(T_CSO_State_Set), intent(inout) :: self integer, intent(out) :: status ! --- const ---------------------------------- @@ -2863,7 +2645,7 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(in) :: self + class(T_CSO_State_Set), intent(in) :: self integer, intent(out) :: status integer, intent(in), optional :: id @@ -2906,8 +2688,8 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(in) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(in) :: self + type(T_CSO_Data_Set), intent(in) :: sdata integer, intent(out) :: status integer, intent(in), optional :: id @@ -2942,8 +2724,8 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(in) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(in) :: self + type(T_CSO_Data_Set), intent(in) :: sdata integer, intent(in) :: ipix integer, intent(out) :: status @@ -2989,8 +2771,8 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(inout) :: self + type(T_CSO_Data_Set), intent(in) :: sdata integer, intent(in) :: ipix integer, intent(out) :: status @@ -3067,8 +2849,8 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(inout) :: self + type(T_CSO_Data_Set), intent(in) :: sdata integer, intent(out) :: status character(len=*), intent(in), optional :: outkey @@ -3123,8 +2905,8 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(inout) :: self + type(T_CSO_Data_Set), intent(in) :: sdata integer, intent(out) :: status character(len=*), intent(in), optional :: outkey @@ -3174,8 +2956,8 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self - type(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(inout) :: self + type(T_CSO_Data_Set), intent(in) :: sdata integer, intent(in) :: ipix real, intent(out) :: hf(:) ! (nlayer) integer, intent(out) :: status @@ -3279,23 +3061,25 @@ contains subroutine CSO_Sat_State_PutOut( self, sdata, filename, status, & - outkey ) + outkey, written ) use NetCDF , only : NF90_Create, NF90_Close use NetCDF , only : NF90_Def_Dim use NetCDF , only : NF90_EndDef use CSO_Comm , only : csoc + use CSO_File , only : CSO_CheckDir use CSO_NcFile, only : T_NcFile ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self - class(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(inout) :: self + class(T_CSO_Data_Set), intent(in) :: sdata character(len=*), intent(in) :: filename integer, intent(out) :: status character(len=*), intent(in), optional :: outkey + logical, intent(out), optional :: written ! --- const ---------------------------------- @@ -3325,12 +3109,19 @@ contains nullify( vars ) end if + ! copy flag? + if ( present(written) ) written = sdata%slc%nout > 0 + ! any pixels to be put out? if ( sdata%slc%nout > 0 ) then ! collect on root, this is much faster than parallel write ... if ( csoc%root ) then + ! create directory if needed: + call CSO_CheckDir( filename, status ) + IF_NOT_OK_RETURN(status=1) + ! create file: call ncf%Init( filename, 'w', status ) IF_NOT_OK_RETURN(status=1) @@ -3394,8 +3185,8 @@ contains ! --- in/out --------------------------------- - class(T_CSO_Sat_State), intent(inout) :: self - class(T_CSO_Sat_Data), intent(in) :: sdata + class(T_CSO_State_Set), intent(inout) :: self + class(T_CSO_Data_Set), intent(in) :: sdata character(len=*), intent(in) :: filename integer, intent(out) :: status diff --git a/oper/src/tutorial_oper_S5p.F90 b/oper/src/tutorial_oper_S5p.F90 index 85b09b6..a91bf5b 100644 --- a/oper/src/tutorial_oper_S5p.F90 +++ b/oper/src/tutorial_oper_S5p.F90 @@ -5,12 +5,15 @@ ! ! CHANGES ! -! 2023-08, Arjo Segers -! Check on `pix_all` value to know if pixels are assigned -! to any of the sub-domains. +! 2023-08, Arjo Segers +! Check on `pix_all` value to know if pixels are assigned +! to any of the sub-domains. ! -! 2024-09, Arjo Segers -! Ensure that pixels crossing the dateline are simulated correctly. +! 2024-09, Arjo Segers +! Ensure that pixels crossing the dateline are simulated correctly. +! +! 2026-07, Arjo Segers +! Pixel selection based on time values rather than orbit file. ! ! !### macro's ##################################################### @@ -37,9 +40,10 @@ program Tutorial_Oper_S5p use CSO, only : T_CSO_RcFile use CSO, only : T_CSO_DateTime, T_CSO_TimeDelta, CSO_DateTime, CSO_TimeDelta, Pretty, & operator(+), operator(-), operator(*), operator(>=) - use CSO, only : T_CSO_Listing - use CSO, only : T_CSO_Sat_Data - use CSO, only : T_CSO_Sat_State + use CSO, only : T_CSO_Data_Series + use CSO, only : T_CSO_State_Series + use CSO, only : T_CSO_Data_Set + use CSO, only : T_CSO_State_Set use CSO, only : T_GridMapping use CSO, only : PolygonOverlapsDomain @@ -112,15 +116,17 @@ program Tutorial_Oper_S5p ! orbit listing: character(len=1024) :: listing_filename - type(T_CSO_Listing) :: listing - ! orbit files: - character(len=1024) :: orbit_filename ! orbit data and state: - type(T_CSO_Sat_Data) :: sdata - type(T_CSO_Sat_State) :: sstate - ! output filename: - character(len=1024) :: output_filename - + type(T_CSO_Data_Set), pointer :: sdata + type(T_CSO_State_Set), pointer :: sstate + ! series of data and state objects: + type(T_CSO_Data_Series) :: dSeries + type(T_CSO_State_Series) :: sSeries + ! dataset loop: + integer :: nset + integer :: iset + logical :: isnew + ! global pixel arrays: integer :: nglb real, pointer :: glb_lon(:) ! (nglb) @@ -146,6 +152,11 @@ program Tutorial_Oper_S5p integer :: iudim character(len=32) :: udimname + ! time selection: + integer :: ntsel + integer :: itsel + integer, pointer :: ipix__tsel(:) + ! user defined variables: integer :: nuvar integer :: iuvar @@ -389,13 +400,15 @@ program Tutorial_Oper_S5p !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! info ... - write (logfu,'(a,": read listing file: ",a)') rname, trim(listing_filename) - ! read listing file: - call listing%Init( listing_filename, status ) + write (logfu,'(a,": initialize series of data files;")') rname + write (logfu,'(a,": use listing file: ",a)') rname, trim(listing_filename) + ! init data series; + ! use listing to search for files with overlapping timerange; + call dSeries%Init( listing_filename, status ) + IF_NOT_OK_RETURN(status=1) + ! init state series: + call sSeries%Init( dSeries, 'state', status ) IF_NOT_OK_RETURN(status=1) - !! show content: - !call listing%Show( status ) - !IF_NOT_OK_RETURN(status=1) ! model concentration array: allocate( conc(nlon,nlat,nz), stat=status ) @@ -425,12 +438,12 @@ program Tutorial_Oper_S5p ! pressure units .. hp_units = 'Pa' - ! cloud cover with single layer + ! cloud cover with single layer over half of the domain allocate( cc(nlon,nlat,nz), stat=status ) IF_NOT_OK_RETURN(status=1) cc = 0.0 - cc(:,:,4) = 1.0 - cc(:,:,5) = 0.2 + cc(:,1:int(nlat/2),4) = 1.0 + cc(:,1:int(nlat/2),5) = 0.2 ! total cloud cover: allocate( tcc(nlon,nlat), stat=status ) IF_NOT_OK_RETURN(status=1) @@ -444,80 +457,121 @@ program Tutorial_Oper_S5p t = t1 do ! info ... + write (logfu,'(a,":")') rname write (logfu,'(a,": time : ",a)') rname, trim(Pretty(t)) + write (logfu,'(a,":")') rname - ! get name of CSO file for orbit with time average - ! in interval around current time: - call listing%SearchFile( t-0.5*dt, t+0.5*dt, 'aver', orbit_filename, status ) + ! select orbit filenames for time interval around current time; + ! mark objects that are not neeeded anynmore + ! keep objects that are still needed + ! initialize new objects + call dSeries%Setup( rcF, rcbase, t-0.5*dt, t+0.5*dt, status ) IF_NOT_OK_RETURN(status=1) - !! testing ... - !if ( t%hour /= 11 ) orbit_filename = '' - - ! defined? - if ( len_trim(orbit_filename) > 0 ) then - - ! info .. - write (logfu,'(a,": orbit file : ",a)') rname, trim(orbit_filename) + ! number of open data sets: + call dSeries%Get( status, nset=nset ) + IF_NOT_OK_RETURN(status=1) + ! info ... + write (logfu,'(a,": number of open datasets: ",i0)') rname, nset - ! ~ orbit data + ! loop over open data sets: + do iset = 1, nset + ! info ... + write (logfu,'(a,": dataset ",i0," / ",i0)') rname, iset, nset - ! inititalize orbit data, - ! read all footprints available in file: - call sdata%Init( rcF, rcbase, orbit_filename, status ) + ! get pointer to data object, obtain flag whether this data is new: + call dSeries%Get( status, iset=iset, sdata=sdata, isnew=isnew ) IF_NOT_OK_RETURN(status=1) + ! new dataset? + if ( isnew ) then + + ! obtain from the orbit: + ! - number of footprints (global, all pixels in file) + ! - pointers to arrays with footprint centers + ! - pointers to arrays with footprint corners + ! - pointer to to selection flags, should be used to select pixels overlapping with local domain + call sdata%Get( status, nglb=nglb, & + glb_lon=glb_lon, glb_lat=glb_lat, & + glb_clons=glb_clons, glb_clats=glb_clats, & + glb_select=glb_select ) + IF_NOT_OK_RETURN(status=1) - ! obtain from the orbit: - ! - number of footprints (global, all pixels in file) - ! - pointers to arrays with footprint centers - ! - pointers to arrays with footprint corners - ! - pointer to to selection flags, should be used to select pixels overlapping with local domain - call sdata%Get( status, nglb=nglb, & - glb_lon=glb_lon, glb_lat=glb_lat, & - glb_clons=glb_clons, glb_clats=glb_clats, & - glb_select=glb_select ) - IF_NOT_OK_RETURN(status=1) + ! select pixels that overlap with local domain; + ! loop over global pixels: + do iglb = 1, nglb + !!~ testing, select all: + !glb_select(1:nglb) = .true. + !!~ deceide if pixel center is part of local domain: + !call PointInsideDomain( glb_lon(iglb), glb_lat(iglb), & + ! (/dom_west,dom_east,dom_south,dom_north/), & + ! glb_select(iglb), status ) + !~ check if footprint overlaps (partly) with local domain, + ! and if it is entirely in global domain: + call PolygonOverlapsDomain( glb_clons(:,iglb), glb_clats(:,iglb), & + (/dom_west,dom_east,dom_south,dom_north/), & + glb_select(iglb), status, & + inside_domain=(/west,east,south,north/) ) + IF_NOT_OK_RETURN(status=1) + end do ! iglb - ! select pixels that overlap with local domain; - ! loop over global pixels: - do iglb = 1, nglb - !!~ testing, select all: - !glb_select(1:nglb) = .true. - !!~ deceide if pixel center is part of local domain: - !call PointInsideDomain( glb_lon(iglb), glb_lat(iglb), & - ! (/dom_west,dom_east,dom_south,dom_north/), & - ! glb_select(iglb), status ) - !~ check if footprint overlaps (partly) with local domain, - ! and if it is entirely in global domain: - call PolygonOverlapsDomain( glb_clons(:,iglb), glb_clats(:,iglb), & - (/dom_west,dom_east,dom_south,dom_north/), & - glb_select(iglb), status, & - inside_domain=(/west,east,south,north/) ) + ! read orbit, locally store only pixels that are flagged in 'glb_select' + ! as having overlapping with this domain: + call sdata%Read( rcF, rcbase, status ) IF_NOT_OK_RETURN(status=1) - end do ! iglb - ! read orbit, locally store only pixels that are flagged in 'glb_select' - ! as having overlapping with this domain: - call sdata%Read( rcF, rcbase, status ) - IF_NOT_OK_RETURN(status=1) + ! number of local pixels read: + call sdata%Get( status, npix=npix ) + IF_NOT_OK_RETURN(status=1) + ! any local pixels? + if ( npix > 0 ) then - ! ~ simulation + ! pointers to (*,pixel) arrays: + ! - footprint centers (not used, for inspiration ..) + ! - footprint corners + call sdata%Get( status, lons=lons, lats=lats, clons=clons, clats=clats ) + IF_NOT_OK_RETURN(status=1) - ! obtain info on track: - ! - number of local pixels - ! - total number of local pixels over all sub-domains - call sdata%Get( status, npix=npix, npix_all=npix_all ) + ! info ... + write (logfu,'(a,": compute mapping weights ...")') rname + ! get pointers to mapping arrays for all pixels: + ! - areas(1:npix) : pixel area [m2] + ! - iw0(1:npix), nw(1:npix) : offset and number of elements in ii/jj/ww + ! - ii(:), jj(:), ww(:) : cell and weight arrays for mapping to footprint, + call GridMapping%GetWeights( clons, clats, & + areas, iw0, nw, ii, jj, ww, status ) + IF_NOT_OK_RETURN(status=1) + + ! info ... + write (logfu,'(a,": store ...")') rname + ! store mapping weights, might be saved and used to create gridded averages; + ! cell indices ii/jj need to be the global index numbers! + call sdata%SetMapping( areas, nw, dom_ilon0+ii, dom_ilat0+jj, ww, status ) + IF_NOT_OK_RETURN(status=1) + + end if ! npix > 0 + + end if ! isnew + + end do ! iset + + ! setup state objects for open data files: + call sSeries%Setup( dSeries, rcF, rcbase, status ) + IF_NOT_OK_RETURN(status=1) + + ! loop over open data sets: + do iset = 1, nset + + ! get pointer to data object, + ! and flag to know whether this is new data: + call dSeries%Get( status, iset=iset, sdata=sdata, isnew=isnew ) IF_NOT_OK_RETURN(status=1) - ! any pixels at all? - if ( npix_all > 0 ) then + ! get pointer to state object: + call sSeries%Get( dSeries, status, iset=iset, sstate=sstate ) + IF_NOT_OK_RETURN(status=1) - ! initialize simulation state; - ! optional arguments: - ! description='long name' : used for output attributes - call sstate%Init( sdata, rcF, rcbase, status, & - description='simulated retrievals' ) - IF_NOT_OK_RETURN(status=1) + ! new state? + if ( isnew ) then ! number of user defined dimensions: call sstate%Get( status, nudim=nudim ) @@ -550,227 +604,181 @@ program Tutorial_Oper_S5p call sstate%EndDef( status ) IF_NOT_OK_RETURN(status=1) - ! any local pixels? - if ( npix > 0 ) then - - ! pointers to (*,pixel) arrays: - ! - footprint centers (not used, for inspiration ..) - ! - footprint corners - call sdata%Get( status, lons=lons, lats=lats, clons=clons, clats=clats ) - IF_NOT_OK_RETURN(status=1) - - ! info ... - write (logfu,'(a,": compute mapping weights ...")') rname - ! get pointers to mapping arrays for all pixels: - ! - areas(1:npix) : pixel area [m2] - ! - iw0(1:npix), nw(1:npix) : offset and number of elements in ii/jj/ww - ! - ii(:), jj(:), ww(:) : cell and weight arrays for mapping to footprint, - call GridMapping%GetWeights( clons, clats, & - areas, iw0, nw, ii, jj, ww, status ) - IF_NOT_OK_RETURN(status=1) - - ! info ... - write (logfu,'(a,": store ...")') rname - ! store mapping weights, might be saved and used to create gridded averages; - ! cell indices ii/jj need to be the global index numbers! - call sdata%SetMapping( areas, nw, dom_ilon0+ii, dom_ilat0+jj, ww, status ) - IF_NOT_OK_RETURN(status=1) - - ! number of user defined variables: - call sstate%Get( status, nuvar=nuvar ) - IF_NOT_OK_RETURN(status=1) - ! any user defined? - if ( nuvar > 0 ) then - - ! storage for variable names and units: - allocate( uvarnames(nuvar), stat=status ) - IF_NOT_OK_RETURN(status=1) - allocate( uvarunits(nuvar), stat=status ) - IF_NOT_OK_RETURN(status=1) - ! fill: - call sstate%Get( status, uvarnames=uvarnames, uvarunits=uvarunits ) - IF_NOT_OK_RETURN(status=1) - - ! loop over variables: - do iuvar = 1, nuvar - ! info .. - write (logfu,'(a,": user defined variable: ",a)') rname, trim(uvarnames(iuvar)) - ! switch: - select case ( trim(uvarnames(iuvar)) ) - - !~ model concentrations - case ( 'mod_conc' ) - ! check units: - if ( trim(uvarunits(iuvar)) /= trim(conc_units) ) then - write (logfu,'("ERROR - variable `",a,"` requires conversion to `",a,"` from `",a,"`")') & - trim(uvarnames(iuvar)), trim(uvarunits(iuvar)), trim(conc_units) - TRACEBACK; status=1; call Exit(status) - end if - ! get pointer to target array with shape (nz,npix): - call sstate%GetData( status, name=uvarnames(iuvar), data1=data1 ) - IF_NOT_OK_RETURN(status=1) - ! loop over pixels: - do ipix = 1, npix - ! any source contributions? - if ( nw(ipix) > 0 ) then - ! init sum: - data1(:,ipix) = 0.0 - ! loop over source contributions: - do iw = iw0(ipix)+1, iw0(ipix)+nw(ipix) - ! add contribution: - data1(:,ipix) = data1(:,ipix) + conc(ii(iw),jj(iw),:) * ww(iw)/areas(ipix) - end do ! iw - end if ! nw > 0 - end do ! ipix - - !~ model half-level pressures: - case ( 'mod_hp' ) - ! check units: - if ( trim(uvarunits(iuvar)) /= trim(hp_units) ) then - write (logfu,'("ERROR - variable `",a,"` requires conversion to `",a,"` from `",a,"`")') & - trim(uvarnames(iuvar)), trim(uvarunits(iuvar)), trim(hp_units) - TRACEBACK; status=1; call Exit(status) - end if - ! get pointer to target array with shape (nz,npix): - call sstate%GetData( status, name=uvarnames(iuvar), data1=data1 ) - IF_NOT_OK_RETURN(status=1) - ! loop over pixels: - do ipix = 1, npix - ! any source contributions? - if ( nw(ipix) > 0 ) then - ! init sum: - data1(:,ipix) = 0.0 - ! loop over source contributions: - do iw = iw0(ipix)+1, iw0(ipix)+nw(ipix) - ! add contribution: - data1(:,ipix) = data1(:,ipix) + hp(ii(iw),jj(iw),:) * ww(iw)/areas(ipix) - end do ! iw - end if ! nw > 0 - end do ! ipix - - !~ cloud cover: total: - case ( 'mod_tcc' ) - ! check units: - if ( trim(uvarunits(iuvar)) /= trim(cc_units) ) then - write (logfu,'("ERROR - variable `",a,"` requires conversion to `",a,"` from `",a,"`")') & - trim(uvarnames(iuvar)), trim(uvarunits(iuvar)), trim(cc_units) - TRACEBACK; status=1; call Exit(status) - end if - ! get pointer to target array with shape (npix): - call sstate%GetData( status, name=uvarnames(iuvar), data0=data0 ) - IF_NOT_OK_RETURN(status=1) - ! loop over pixels: - do ipix = 1, npix - ! any source contributions? - if ( nw(ipix) > 0 ) then - ! init sum: - data0(ipix) = 0.0 - ! loop over source contributions: - do iw = iw0(ipix)+1, iw0(ipix)+nw(ipix) - ! add contribution: - data0(ipix) = data0(ipix) + tcc(ii(iw),jj(iw)) * ww(iw)/areas(ipix) - end do ! iw - end if ! nw > 0 - end do ! ipix - - !~ cloud cover: profile: - case ( 'mod_cc' ) - ! check units: - if ( trim(uvarunits(iuvar)) /= trim(cc_units) ) then - write (logfu,'("ERROR - variable `",a,"` requires conversion to `",a,"` from `",a,"`")') & - trim(uvarnames(iuvar)), trim(uvarunits(iuvar)), trim(cc_units) - TRACEBACK; status=1; call Exit(status) - end if - ! get pointer to target array with shape (nz,npix): - call sstate%GetData( status, name=uvarnames(iuvar), data1=data1 ) - IF_NOT_OK_RETURN(status=1) - ! loop over pixels: - do ipix = 1, npix - ! any source contributions? - if ( nw(ipix) > 0 ) then - ! init sum: - data1(:,ipix) = 0.0 - ! loop over source contributions: - do iw = iw0(ipix)+1, iw0(ipix)+nw(ipix) - ! add contribution: - data1(:,ipix) = data1(:,ipix) + cc(ii(iw),jj(iw),:) * ww(iw)/areas(ipix) - end do ! iw - end if ! nw > 0 - end do ! ipix - - !~ - case default - write (logfu,'("unsupported variable name `",a,"`")') trim(uvarnames(iuvar)) - TRACEBACK; status=1; call Exit(status) - end select - - end do ! user defined variables - - ! info ... - write (logfu,'(a,": end ...")') rname - - ! clear: - deallocate( uvarnames, stat=status ) - IF_NOT_OK_RETURN(status=1) - deallocate( uvarunits, stat=status ) - IF_NOT_OK_RETURN(status=1) + end if ! isnew - end if ! nuvar > 0 - - end if ! npix > 0 - - ! ~ formula + ! ~ simulation - ! inquire info on pixels covering multiple domains, - ! setup exchange parameters: - call sdata%SetupExchange( status ) - IF_NOT_OK_RETURN(status=1) - ! exchange simulations: - call sstate%Exchange( sdata, status ) - IF_NOT_OK_RETURN(status=1) + ! number of local pixels read: + call sdata%Get( status, npix=npix ) + IF_NOT_OK_RETURN(status=1) - ! apply formula: kernel convolution etc: - call sstate%ApplyFormulas( sdata, status ) - IF_NOT_OK_RETURN(status=1) + ! setup selection of pixels to be simulated for this time; + ! integer array ipix__tsel(1:ntsel) has the selected pixel indices: + call sdata%GetTimeSelection( t1, t2, ntsel, ipix__tsel, status ) + IF_NOT_OK_RETURN(status=1) + ! info .. + write (logfu,'(a,": selected ",i6," / ",i6," pixels for current time")') rname, ntsel, npix - ! ~ put out + ! any pixels selected? + if ( ntsel > 0 ) then - ! target file for selected data: - write (output_filename,'("CSO_output_",i4.4,2i2.2,"_",2i2.2,"_data.nc")') & - t%year, t%month, t%day, t%hour, t%minute - ! setup gathering information, gather output arrays, write: - call sdata%PutOut( output_filename, status ) + ! pointers to mapping arrays that were stored when this data was new; + ! note that ii and jj are global indices! + call sdata%Get( status, nw=nw, iw0=iw0, ii=ii, jj=jj, ww=ww, areas=areas ) IF_NOT_OK_RETURN(status=1) - ! target file for simulation state: - write (output_filename,'("CSO_output_",i4.4,2i2.2,"_",2i2.2,"_state.nc")') & - t%year, t%month, t%day, t%hour, t%minute - ! gether output arrays, write: - call sstate%PutOut( sdata, output_filename, status ) + ! number of user defined variables: + call sstate%Get( status, nuvar=nuvar ) IF_NOT_OK_RETURN(status=1) + ! any user defined? + if ( nuvar > 0 ) then - ! ~ clear + ! storage for variable names and units: + allocate( uvarnames(nuvar), stat=status ) + IF_NOT_OK_RETURN(status=1) + allocate( uvarunits(nuvar), stat=status ) + IF_NOT_OK_RETURN(status=1) + ! fill: + call sstate%Get( status, uvarnames=uvarnames, uvarunits=uvarunits ) + IF_NOT_OK_RETURN(status=1) - ! clear: - nullify( glb_lon ) - nullify( glb_lat ) - nullify( glb_select ) + ! loop over variables: + do iuvar = 1, nuvar + ! info .. + write (logfu,'(a,": user defined variable: ",a)') rname, trim(uvarnames(iuvar)) + ! switch: + select case ( trim(uvarnames(iuvar)) ) + + !~ model concentrations + case ( 'mod_conc' ) + ! check units: + if ( trim(uvarunits(iuvar)) /= trim(conc_units) ) then + write (logfu,'("ERROR - variable `",a,"` requires conversion to `",a,"` from `",a,"`")') & + trim(uvarnames(iuvar)), trim(uvarunits(iuvar)), trim(conc_units) + TRACEBACK; status=1; call Exit(status) + end if + ! get pointer to target array with shape (nz,npix): + call sstate%GetData( status, name=uvarnames(iuvar), data1=data1 ) + IF_NOT_OK_RETURN(status=1) + ! loop over time selection: + do itsel = 1, ntsel + ! pixel index: + ipix = ipix__tsel(itsel) + ! any source contributions? + if ( nw(ipix) > 0 ) then + ! init sum: + data1(:,ipix) = 0.0 + ! loop over source contributions: + do iw = iw0(ipix)+1, iw0(ipix)+nw(ipix) + ! add contribution: + data1(:,ipix) = data1(:,ipix) + conc(ii(iw),jj(iw),:) * ww(iw)/areas(ipix) + end do ! iw + end if ! nw > 0 + end do ! itsel + + !~ model half-level pressures: + case ( 'mod_hp' ) + ! check units: + if ( trim(uvarunits(iuvar)) /= trim(hp_units) ) then + write (logfu,'("ERROR - variable `",a,"` requires conversion to `",a,"` from `",a,"`")') & + trim(uvarnames(iuvar)), trim(uvarunits(iuvar)), trim(hp_units) + TRACEBACK; status=1; call Exit(status) + end if + ! get pointer to target array with shape (nz,npix): + call sstate%GetData( status, name=uvarnames(iuvar), data1=data1 ) + IF_NOT_OK_RETURN(status=1) + ! loop over time selection: + do itsel = 1, ntsel + ! pixel index: + ipix = ipix__tsel(itsel) + ! any source contributions? + if ( nw(ipix) > 0 ) then + ! init sum: + data1(:,ipix) = 0.0 + ! loop over source contributions: + do iw = iw0(ipix)+1, iw0(ipix)+nw(ipix) + ! add contribution: + data1(:,ipix) = data1(:,ipix) + hp(ii(iw),jj(iw),:) * ww(iw)/areas(ipix) + end do ! iw + end if ! nw > 0 + end do ! itsel + + !~ cloud cover: total: + case ( 'mod_tcc' ) + ! check units: + if ( trim(uvarunits(iuvar)) /= trim(cc_units) ) then + write (logfu,'("ERROR - variable `",a,"` requires conversion to `",a,"` from `",a,"`")') & + trim(uvarnames(iuvar)), trim(uvarunits(iuvar)), trim(cc_units) + TRACEBACK; status=1; call Exit(status) + end if + ! get pointer to target array with shape (npix): + call sstate%GetData( status, name=uvarnames(iuvar), data0=data0 ) + IF_NOT_OK_RETURN(status=1) + ! loop over time selection: + do itsel = 1, ntsel + ! pixel index: + ipix = ipix__tsel(itsel) + ! any source contributions? + if ( nw(ipix) > 0 ) then + ! init sum: + data0(ipix) = 0.0 + ! loop over source contributions: + do iw = iw0(ipix)+1, iw0(ipix)+nw(ipix) + ! add contribution: + data0(ipix) = data0(ipix) + tcc(ii(iw),jj(iw)) * ww(iw)/areas(ipix) + end do ! iw + end if ! nw > 0 + end do ! itsel + + !~ cloud cover: profile: + case ( 'mod_cc' ) + ! check units: + if ( trim(uvarunits(iuvar)) /= trim(cc_units) ) then + write (logfu,'("ERROR - variable `",a,"` requires conversion to `",a,"` from `",a,"`")') & + trim(uvarnames(iuvar)), trim(uvarunits(iuvar)), trim(cc_units) + TRACEBACK; status=1; call Exit(status) + end if + ! get pointer to target array with shape (nz,npix): + call sstate%GetData( status, name=uvarnames(iuvar), data1=data1 ) + IF_NOT_OK_RETURN(status=1) + ! loop over time selection: + do itsel = 1, ntsel + ! pixel index: + ipix = ipix__tsel(itsel) + ! any source contributions? + if ( nw(ipix) > 0 ) then + ! init sum: + data1(:,ipix) = 0.0 + ! loop over source contributions: + do iw = iw0(ipix)+1, iw0(ipix)+nw(ipix) + ! add contribution: + data1(:,ipix) = data1(:,ipix) + cc(ii(iw),jj(iw),:) * ww(iw)/areas(ipix) + end do ! iw + end if ! nw > 0 + end do ! itsel + + !~ + case default + write (logfu,'("unsupported variable name `",a,"`")') trim(uvarnames(iuvar)) + TRACEBACK; status=1; call Exit(status) + end select + + end do ! user defined variables - ! done with state: - call sstate%Done( sdata, status ) - IF_NOT_OK_RETURN(status=1) + ! info ... + write (logfu,'(a,": end ...")') rname - end if ! npix_all > 0 + ! clear: + deallocate( uvarnames, stat=status ) + IF_NOT_OK_RETURN(status=1) + deallocate( uvarunits, stat=status ) + IF_NOT_OK_RETURN(status=1) - ! done with data: - call sdata%Done( status ) - IF_NOT_OK_RETURN(status=1) + end if ! nuvar > 0 - !! testing ... - !write (logfu,'("break after first orbit")') - !exit + end if ! ntsel > 0 - end if ! orbit found in listing + end do ! data sets ! end? if ( t >= t2 ) exit @@ -778,8 +786,18 @@ program Tutorial_Oper_S5p t = t + dt end do ! time loop - ! done with listing: - call listing%Done( status ) + ! postprocess currently open files: + call dSeries%Close( status ) + IF_NOT_OK_RETURN(status=1) + ! close states too: + call sSeries%Close( dSeries, status ) + IF_NOT_OK_RETURN(status=1) + + ! done with series of states: + call sSeries%Done( dSeries, status ) + IF_NOT_OK_RETURN(status=1) + ! done with series: + call dSeries%Done( status ) IF_NOT_OK_RETURN(status=1) ! done with grid mapping: -- GitLab From 765c50cf5d8e3aae61814d625fdc73cdd62fcfe6 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 16:00:02 +0200 Subject: [PATCH 11/20] Updated cleaning of temporary files. --- oper/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/oper/Makefile b/oper/Makefile index 1a178a8..25584bb 100644 --- a/oper/Makefile +++ b/oper/Makefile @@ -24,7 +24,7 @@ clean: rm -f tutorial*.err* rm -f tutorial*.log* rm -f cso.log* - rm -f CSO_*.nc + rm -f *.nc clean-all: clean (cd src; make clean) -- GitLab From bec3b3842be3fdb8314af5d96efc726fe6679d68 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Tue, 21 Jul 2026 16:00:16 +0200 Subject: [PATCH 12/20] Updated change log. --- oper/CHANGELOG | 46 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/oper/CHANGELOG b/oper/CHANGELOG index 7e51972..d09b2fa 100644 --- a/oper/CHANGELOG +++ b/oper/CHANGELOG @@ -8,7 +8,7 @@ v2.6 Improved packing of constant variables. src/cso_ncfile.F90 -Do not pack 'time' variables. +Do not pack `time` variables. src/cso_pixels.F90 Support files with "sample" coordinate instead of "pixel", @@ -39,7 +39,7 @@ Added "GetWeightsCell" methodes to assign points to grid cells. Trap undefined values in kernel application. src/cso_pixels.F90 -Added formula 'TropoPressures' to create variable with half-level pressures +Added formula `TropoPressures` to create variable with half-level pressures at surface, top of boundary layer, and top of troposhere. src/cso_pixels.F90 src/cso_tools.F90 @@ -57,7 +57,7 @@ Ensure that pixels crossing the dateline are simulated correctly. v2.11 ----- -Removed comment from 'use MPI' statements after compiler complained. +Removed comment from `use MPI` statements after compiler complained. src/cso_comm.F90 Corrected spell error. @@ -66,7 +66,7 @@ Corrected spell error. Support 2D real arrays. src/cso_parray.F90 -Support labels 'hours', 'minutes', 'seconds', etc. +Support labels `hours`, `minutes`, `seconds`, etc. src/cso_datetime.F90 @@ -80,14 +80,46 @@ Changes to avoid compiler warnings. Updated comment. src/cso_mapping.F90 -Support 'fill_value' argument when reading 1D real array. +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. +Added `LinInterpWeights` routine to calculate interpolation weights. src/cso_tools.F90 -! +Support latest servers. + src/Makefile_config_metno-nebula # deleted + src/Makefile_config_metno-fahrenheit + src/Makefile_config_tno-hpc3 # deleted + src/Makefile_config_tno-hpc4 + launcher + launcher-mpi + +Added routines to broadcast/scatter 1D integer or character arrays. + src/cso_comm.F90 + +Added `CompressTime` routine. + src/cso_datetime.F90 + +Added `CSO_GetBasename` and `CSO_SplitExt` routines. + src/cso_file.F90 + +Support listings with second `filename`, this is for example used to specify both data and state files. +Added method `Create` to initialize an empty listing that will be filled and saved. +Added method `Select` to select records based on time interval. + src/cso_listing.F90 + +Updated logging messages. + src/cso_pixels.F90 + +Introduced `Series` object to open multiple files from a time series for a single interval. +Pixel selection based on time values. + src/cso_sat.F90 + src/cso_series.F90 + src/cso.F90 + src/tutorial_oper_S5p.F90 + src/Makefile_deps + -- GitLab From acbc7df7eaad1699793c2c02d24bec41ec2c525d Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Wed, 22 Jul 2026 11:14:46 +0200 Subject: [PATCH 13/20] Moved listing class to seperate module. --- doc/source/pymod-cso_listing.rst | 5 + doc/source/pymods.rst | 1 + src/cso/__init__.py | 5 + src/cso/cso_catalogue.py | 6 +- src/cso/cso_colhub.py | 21 +- src/cso/cso_earthaccess.py | 15 +- src/cso/cso_file.py | 567 ----------------------------- src/cso/cso_gridded.py | 6 +- src/cso/cso_inquire.py | 6 +- src/cso/cso_listing.py | 604 +++++++++++++++++++++++++++++++ src/cso/cso_s5p.py | 14 +- src/cso/cso_superobs.py | 6 +- src/cso/cso_viirs.py | 6 +- 13 files changed, 672 insertions(+), 590 deletions(-) create mode 100644 doc/source/pymod-cso_listing.rst create mode 100644 src/cso/cso_listing.py diff --git a/doc/source/pymod-cso_listing.rst b/doc/source/pymod-cso_listing.rst new file mode 100644 index 0000000..49ad7f7 --- /dev/null +++ b/doc/source/pymod-cso_listing.rst @@ -0,0 +1,5 @@ +.. Documentation for module. + +.. Import documentation from ".py" file: +.. automodule:: cso.cso_listing + diff --git a/doc/source/pymods.rst b/doc/source/pymods.rst index ce62f4a..b156b79 100644 --- a/doc/source/pymods.rst +++ b/doc/source/pymods.rst @@ -35,6 +35,7 @@ Overview of the Python module(s) and classes. pymod-cso_s5p pymod-cso_viirs pymod-cso_file + pymod-cso_listing pymod-cso_gridded pymod-cso_superobs pymod-cso_colocate diff --git a/src/cso/__init__.py b/src/cso/__init__.py index bdeab0b..a181724 100644 --- a/src/cso/__init__.py +++ b/src/cso/__init__.py @@ -14,6 +14,9 @@ # 2025-02, Arjo Segers # Support VIIRS AOD data. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# # ------------------------------------------------- @@ -35,6 +38,7 @@ The ``cso`` module is the top level access to the various Actual implementations can be found in submodules: * :py:mod:`cso.cso_file` +* :py:mod:`cso.cso_listing` * :py:mod:`cso.cso_inquire` * :py:mod:`cso.cso_dataspace` * :py:mod:`cso.cso_pal` @@ -138,6 +142,7 @@ except: # import entities from sub-modules: from .cso_file import * +from .cso_listing import * from .cso_inquire import * from .cso_dataspace import * from .cso_colhub import * diff --git a/src/cso/cso_catalogue.py b/src/cso/cso_catalogue.py index 916a565..732adcd 100644 --- a/src/cso/cso_catalogue.py +++ b/src/cso/cso_catalogue.py @@ -40,6 +40,9 @@ # Extended unit conversions. # Fixed string formatting. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# ######################################################################## @@ -319,6 +322,7 @@ class CSO_Catalogue(CSO_CatalogueBase): # tools: from . import cso_file + from . import cso_listing from . import cso_plot # setup graphics back-end: @@ -353,7 +357,7 @@ class CSO_Catalogue(CSO_CatalogueBase): raise Exception # endif # read listing: - lst = cso_file.CSO_Listing(lstfile, indent=indent) + lst = cso_listing.CSO_Listing(lstfile, indent=indent) # subset of files in timerange: xlst = lst.Select(tr=(t1, t2)) # number of records: diff --git a/src/cso/cso_colhub.py b/src/cso/cso_colhub.py index 89decab..f3761f7 100644 --- a/src/cso/cso_colhub.py +++ b/src/cso/cso_colhub.py @@ -35,6 +35,9 @@ # 2026-04, Arjo Segers # Skip files that cannot be accessed, for example a broken link. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# ######################################################################## ### @@ -162,7 +165,7 @@ class CSO_ColHubMirror_Inquire(utopya.UtopyaRc): import fnmatch # tools: - from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -190,7 +193,7 @@ class CSO_ColHubMirror_Inquire(utopya.UtopyaRc): # initiallize for (re)creation, # define extra columns to ensure correct content even for empty files: - listing = cso_file.CSO_Listing( columns=["orbit","processing","collection","processor_version","production_time","href"] ) + listing = cso_listing.CSO_Listing( columns=["orbit","processing","collection","processor_version","production_time","href"] ) # archive directories archive_dirs = self.GetSetting("dir").split() @@ -402,7 +405,7 @@ class CSO_ColHubMirror_Missing(utopya.UtopyaRc): import fnmatch # tools: - from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -430,7 +433,7 @@ class CSO_ColHubMirror_Missing(utopya.UtopyaRc): # initiallize for (re)creation, # define extra columns to ensure correct content even for empty files: - listing = cso_file.CSO_Listing( columns=["orbit","processing","collection","processor_version","production_time","href"] ) + listing = cso_listing.CSO_Listing( columns=["orbit","processing","collection","processor_version","production_time","href"] ) # table with all available files: @@ -441,7 +444,7 @@ class CSO_ColHubMirror_Missing(utopya.UtopyaRc): ) listfile_all = filedate.strftime(listfile_all) # read: - listing_all = cso_file.CSO_Listing(listfile_all) + listing_all = cso_listing.CSO_Listing(listfile_all) # table with currently already available files: listfiles = self.GetSetting("curr.file").split(";") @@ -456,7 +459,7 @@ class CSO_ColHubMirror_Missing(utopya.UtopyaRc): listfile = filedate.strftime(listfile) listfile = listfile.strip() # read: - listings_curr.append(cso_file.CSO_Listing(listfile)) + listings_curr.append(cso_listing.CSO_Listing(listfile)) # endfor # extract orbits: @@ -651,7 +654,7 @@ class CSO_ColHubMirror_Cleanup(utopya.UtopyaRc): import fnmatch # tools: - from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -669,7 +672,7 @@ class CSO_ColHubMirror_Cleanup(utopya.UtopyaRc): ) listfile1 = filedate1.strftime(listfile1) # read: - listing1 = cso_file.CSO_Listing(listfile1) + listing1 = cso_listing.CSO_Listing(listfile1) # table files in extra archive: listfile2= self.GetSetting("mirror2.file") @@ -679,7 +682,7 @@ class CSO_ColHubMirror_Cleanup(utopya.UtopyaRc): ) listfile2 = filedate2.strftime(listfile2) # read: - listing2 = cso_file.CSO_Listing(listfile2) + listing2 = cso_listing.CSO_Listing(listfile2) # remove duplicate files? if not, rename: remove_duplicates = self.GetSetting("mirror2.remove_duplicates", totype="bool") diff --git a/src/cso/cso_earthaccess.py b/src/cso/cso_earthaccess.py index 21c8c67..d7345ce 100644 --- a/src/cso/cso_earthaccess.py +++ b/src/cso/cso_earthaccess.py @@ -16,6 +16,9 @@ # 2025-09, Arjo Segers # Added 'blacklist' for problematic URL's. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# ######################################################################## ### @@ -164,7 +167,7 @@ class CSO_EarthAccess_Inquire(utopya.UtopyaRc): import pandas # tools: - from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -225,7 +228,7 @@ class CSO_EarthAccess_Inquire(utopya.UtopyaRc): logging.info(f"{indent}create %s ..." % lst_file) # initiallize for (re)creation: - listing = cso_file.CSO_Listing() + listing = cso_listing.CSO_Listing() # start of first month: tm = datetime.datetime(t1.year, t1.month, 1) @@ -512,6 +515,7 @@ class CSO_EarthAccess_Download(utopya.UtopyaRc): # tools: from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -549,7 +553,7 @@ class CSO_EarthAccess_Download(utopya.UtopyaRc): # endif # new empty listing: - listing = cso_file.CSO_Listing() + listing = cso_listing.CSO_Listing() # loop: for ifile in range(len(filename__templates)): # current: @@ -571,7 +575,7 @@ class CSO_EarthAccess_Download(utopya.UtopyaRc): # info .. logging.info(f"{indent}read inquire table: %s" % filename) # read: - lst = cso_file.CSO_Listing(filename) + lst = cso_listing.CSO_Listing(filename) # extend: listing.Append(lst) # endfor # filenames @@ -775,6 +779,7 @@ class CSO_EarthAccess_Download_Listing(utopya.UtopyaRc): # tools: from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -817,7 +822,7 @@ class CSO_EarthAccess_Download_Listing(utopya.UtopyaRc): cso_file.CheckDir(lst_file, dmode=dmode) # initiallize for (re)creation: - listing = cso_file.CSO_Listing(indent=f"{indent} ") + listing = cso_listing.CSO_Listing(indent=f"{indent} ") # keep scanned roots for progress info: subdirs = [] diff --git a/src/cso/cso_file.py b/src/cso/cso_file.py index 1cb4e51..9f3e872 100644 --- a/src/cso/cso_file.py +++ b/src/cso/cso_file.py @@ -8,9 +8,6 @@ # 2022-09, Arjo Segers # Fixed computation of packing parameters for array with constant value. # -# 2023-08, Arjo Segers -# Added "method" argument to CSO_Listing.Select(). -# # 2023, Lewis Blake, Arjo Segers # Formatted using "black". # @@ -23,27 +20,14 @@ # 2024-08, Arjo Segers # Updated formatting. # -# 2025-01, Arjo Segers -# Extened CSO_Listing class with selection and other methods. -# # 2025-02, Arjo Segers # Option to CheckDir to set directory creation mode. # Accept dmode=None as default. # # 2025-04, Arjo Segers -# Convert datetime time values in listing files after reading to keep default str values in table. -# Updated methods of CSO_Listing class: -# - GetValues now returns the index column by defaul. -# - added GetValue method -# Removed usage of undefined filename in error messages. -# -# 2025-04, Arjo Segers # Enable zlib compression only for numerical data. # Avoid warnings from packing in case of all-nan values. # Open a file rather than loading it. -# Extended sort options for listing files. -# Support creation of listing file objects without filename. -# Support selection of multiple records from listing file. # # 2025-04, Arjo Segers # Changed imports for python packaging. @@ -55,21 +39,8 @@ # 2025-09, Arjo Segers # Let 'GetData' return all attributes rather than just 'units'. # -# 2025-09, Arjo Segers -# Improved speed of CSO_Listing.Selection method. -# -# 2025-10, Arjo Segers -# Replace ' and ' by ' & ' when selecting records. -# # 2026-01, Arjo Segers # Updated use of Dataset.dims following deprication warning. -# Removed whitespace from template replacement. -# -# 2026-04, Arjo Segers -# Added optional "columns" argument to initialize empty listing object. -# -# 2026-05, Arjo Segers -# Support `name` argument for GetRecord method. # ######################################################################## @@ -85,23 +56,12 @@ Create and access file with satellite data extract. - -Methods -======= - -The following methods are defined: - -* :py:meth:`CheckDir` -* :py:meth:`Pack_DataArray` - - Class hierchy ============= The classes are defined according to the following hierchy: * :py:class:`.CSO_File` -* :py:class:`.CSO_Listing` Methods and classes @@ -1286,533 +1246,6 @@ class CSO_File(object): # endclass CSO_File -######################################################################## -### -### table with extract files -### -######################################################################## - - -class CSO_Listing(object): - - """ - Storage for table with file properties. - The records describe properties of :py:class:`CSO_File` files, - for example the time range of the pixels included. - This could be used to quickly scan an archive for filenames - with pixels within a requested interval. - - The content of the table file looks like:: - - filename;start_time;end_time - 2007/06/RAL-IASI-CH4_20070601T003855.nc;2007-06-01 00:45:36.827000000;2007-06-01 02:23:57.567000000 - 2007/06/RAL-IASI-CH4_20070601T022359.nc;2007-06-01 02:23:59.512000000;2007-06-01 04:05:57.328000000 - ... - - Optional arguments: - - * ``filename`` : listing file read into table - * ``columns`` : extra columns that will be created if an empty listing is initialized - - """ - - def __init__(self, filename=None, columns=[], indent=""): - """ - Initialize empty table or read existing. - """ - - # modules: - import os - import pandas - - # csv style: - self.sep = ";" - # head for index column: - self.index_label = "filename" - - # store filename: - self.filename = filename - - # directory name: - if self.filename is not None: - # check .. - if not os.path.isfile(filename): - logging.error("listing file not found: %s" % filename) - raise Exception - # endif - - # base directory: - self.dirname = os.path.dirname(self.filename) - # could be empty .. - if len(self.dirname) == 0: - self.dirname = os.curdir - # endif - - # info ... - logging.info(f"{indent}read listing {filename} ...") - # read: - self.df = pandas.read_csv( - filename, - sep=self.sep, - index_col=self.index_label, - dtype="str", - ) - # convert datetime columns, as the default is "str" - # the "parse_dates" argument could not be used: - self.df["start_time"] = pandas.to_datetime(self.df["start_time"]) - self.df["end_time"] = pandas.to_datetime(self.df["end_time"]) - - else: - # not defined yet, assume current location: - self.dirname = os.curdir - - # new empty table: - self.df = pandas.DataFrame(columns=["start_time", "end_time"] + columns) - - # endif - - # enddef __init__ - - # * - - def Save(self, filename, dmode=None, indent=""): - """ - Write table to file. - """ - - # info ... - logging.info(f"{indent}save listing {filename} ...") - - # create directory if needed: - CheckDir(filename, dmode=dmode) - - # write all columns: - columns = list(self.df.keys()) - # index only once .. - if self.index_label in columns: - columns.remove(self.index_label) - - # save, also write the index column: - self.df.to_csv(filename, sep=self.sep, columns=columns, index_label=self.index_label) - - # enddef Save - - # * - - def __len__(self): - """ - Returns number of records in table. - """ - - return len(self.df) - - # enddef __len__ - - # * - - def GetRecord(self, irec=None, name=None): - """ - Return record for corresponding 0-based row index or index name. - """ - - # modules: - import numpy - - # by name? - if name is not None: - indx, = numpy.where( self.df.index == name ) - if len(indx) == 0: - logging.error(f"value `{name}` not found in index") - raise Exception - #endif - irec = indx[0] - #endif - - # extract: - rec = self.df.iloc[irec].copy() - # add index values: - rec[self.index_label] = self.df.index[irec] - - # ok - return rec - - # enddef GetRecord - - # * - - def GetValue(self, filename, key): - """ - Return ``key`` value for record of ``filename``. - """ - - # copy: - return self.df.at[filename, key] - - # enddef GetValue - - # * - - def GetValues(self, name=None): - """ - Return :py:class:`pandas.Series` object with all values for provided column ``name``, - or the filenames if ``name`` is not provided. - """ - - # switch .. - if name is None: - return self.df.index - elif name in self.df.keys(): - return self.df[name] - else: - logging.error(f"column '{name}' not found in listing") - raise Exception - # endif - - # enddef GetValues - - # * - - def __contains__(self, key): - """ - Returns True if key is in the table index. - """ - - return key in self.df.index - - # enddef __contains__ - - # * - - def keys(self): - """ - Returns list of column names. - """ - - return self.df.keys() - - # enddef keys - - # * - - def Cleanup(self, indent=""): - """ - Remove records from table if filename is not present anymore. - """ - - # modules: - import os - - # loop over records: - for fname in self.df.index: - # file missing? - if not os.path.isfile(os.path.join(self.dirname, fname)): - # info ... - logging.info(indent + 'drop record "%s", file not present anymore ..' % fname) - # remove record: - self.df.drop(fname, inplace=True) - # endif - # endfor - - # enddef Cleanup - - # * - - def Check(self, filename, indent=""): - """ - Return exception if no record for ``filename`` is present. - """ - - # modules: - import os - - # relative: - fname = os.path.relpath(filename, self.dirname) - - # check .. - if fname not in self.df.index: - logging.error(f"file '{fname}' is not a record in table: {self.filename}") - raise Exception - # endif - - # enddef Check - - # * - - def Update(self, filename, csf, tr=None, xcolumns=[], indent=""): - """ - Add record (or replace) for ``filename`` with information - from the ``csf`` :py:class:`CSO_File` object. - - The time range of pixels is read from the ``csf`` object, - unless the :py:class:`datetime.datetime` tupple ``tr=(t1,t2)`` is supplied. - - The optional ``xcolumns`` specify a list of global attributes - to be added as extra columns. - """ - - # modules: - import os - import numpy - - # timerange: - if tr is None: - # get from file: - tr = csf.GetTimeRange() - # endif - # check, in some products the time values were found to be undefined ... - if numpy.isnan(tr[0]) or numpy.isnan(tr[1]): - logging.error("found undefined time range while updating listing record") - raise Exception - # endif - # store: - self.df.at[filename, "start_time"] = tr[0] - self.df.at[filename, "end_time"] = tr[1] - - # extra: - for name in xcolumns: - # get value: - avalue = csf.GetAttr(name) - # copy: - self.df.at[filename, name] = avalue - # endfor # xcolumns - - # enddef Update - - # * - - def UpdateRecord(self, filename, data, indent=""): - """ - Add record (or replace) for ``filename`` with columns and values defined in the ``data`` dict. - """ - - # modules: - import os - - # loop over columns: - for key in data.keys(): - self.df.at[filename, key] = data[key] - # endfor - - # enddef UpdateRecord - - # * - - def Select( - self, tr=None, method="overlap", expr=None, blacklist=[], verbose=True, indent="", **kwargs - ): - """ - Return :py:class:`CSO_Listing` object with selection of records. - - Optional arguments: - - * ``tr=(t1,t2)`` : select records that could be assigned to this timerange; - requires the ``method`` keyword: - - * ``method='overlap'`` (default) : select records that have pixels overlapping with interval - * ``method='middle'`` : select pixels with the middle of ``start_time`` and ``end_time`` - within the inerval ``(t1,t2]`` - - * ``name=value`` arguments select records based on column values - - * The exprssion ``expr`` provides a list of ';'-seperated selection expressions - with templates ``%{..}`` for the column names:: - - (%{processor_version} == '020400') & (%{processing} == 'RPRO') ; ... - - This is evaluted after the previously described selections. - Eventually skip files in ``blacklist``. - - """ - - # modules: - import os - import pandas - - # init result with entire dataset: - df = self.df - - # * - - # time range specified? - if tr is not None: - # convert if needed: - t1, t2 = pandas.Timestamp(tr[0]), pandas.Timestamp(tr[1]) - # switch: - if method == "overlap": - # select: - df = df[(df["start_time"] <= t2) & (df["end_time"] >= t1)] - # - elif method == "middle": - # mid times: - mid_time = df["start_time"] + (df["end_time"] - df["start_time"]) / 2 - # select: - df = df[(mid_time > t1) & (mid_time <= t2)] - # - else: - logging.error('unsupported method "%s" with argument ``tr``' % method) - raise Exception - # endif - # endif - - # * - - # other selections: - for key, value in kwargs.items(): - # check .. - if key not in df.keys(): - logging.error(f"key '{key}' not defined in listing: {self.filename}") - raise Exception - # endif - # select: - df = df[df[key] == value] - # endfor # key/value pairs - - # * - - # evaluate selection expression? - if expr is not None: - # split into sequnece of critera; first expression with macthing record will be used: - selections = expr.split(";") - - # initialize empty list of selected record indices (filenames): - selected = [] - # storage for status label: "selected", "blacklisted", ... - filestatus = {} - # loop over selection criteria, - # this should give either none or a single file: - for selection in selections: - # replace templates: - # %{orbit} == '12345' - # to: - # df['orbit'] == '12345' - for key in self.df.keys(): - selection = selection.replace(f"%{{{key}}}", f"df['{key}']").strip() - # endfor - ## testing ... - #logging.info(f"{indent}selection `{selection}` ...") - # adhoc: replace some operators that are not allowded anymore ... - selection = selection.replace( " and ", " & " ) - # evaluate: - xdf = df[ eval(selection) ] - # any? - if len(xdf) > 0: - # selected indices (filenames): - selected = [] - for indx, xrec in xdf.iterrows(): - # skip? - if os.path.basename(indx) in blacklist: - filestatus[indx] = "blacklisted" - continue - # endif - # store: - selected.append(indx) - # set status: - filestatus[indx] = "selected" - # endif - # endfor # records - - # any selected? - if len(selected) > 0: - # leave: - break - # endif - - # endfor # selection criteria - - # show selection? - if verbose: - # info ... - logging.info(f"{indent}available record(s):") - # loop: - for fname, row in df.iterrows(): - line = fname - if fname in filestatus.keys(): - line = line + f" [{filestatus[fname]}]" - logging.info(f"{indent} {line}") - # endfor - # endif # verbose - - # no match? - if len(selected) == 0: - # info ... - logging.warning(" no match with any selection criterium ...") - ## degug: - # for selection in selections : - # logging.warning( ' %s' % selection.strip() ) - # logging.warning( ' record:' ) - # for key in rec.keys() : - # logging.warning( ' "%s" : "%s"' % (key,rec[key]) ) - # create empty dataframe as result: - df = pandas.DataFrame(columns=df.columns) - else: - # extract selected record(s): - df = df.loc[selected] - # endif - - # endif # expr defined - - # * - - # return as listing: - lst = CSO_Listing() - lst.df = df - - # ok - return lst - - # enddef Select - - # * - - def Append(self, lst, path=""): - """ - Append content of other listing. - """ - - # modules: - import os - import pandas - - # prefix paths? - if len(path) > 0: - # copy, make "filename" a normal column again: - xdf = lst.df.reset_index() - # add path: - for key in ["filename", "filename_state"]: - xdf[key] = os.path.join(path, xdf[key]) - # endfor - # add rows to current rows, use the original index: - self.df = pandas.concat((self.df, xdf.set_index(self.index_label))) - else: - # append: - if len(self.df) == 0: - self.df = lst.df - else: - self.df = pandas.concat((self.df, lst.df)) - # endif - # endif - - # enddef Append - - # * - - def Sort(self, by=None): - """ - Sort listing table by index (default, this is the "filename") or by a named column. - """ - - # sort index or values: - if by is None: - self.df.sort_index(inplace=True) - else: - self.df.sort_values(by, inplace=True) - # endif - - # endef Sort - - -# endclass CSO_Listing - ######################################################################## ### diff --git a/src/cso/cso_gridded.py b/src/cso/cso_gridded.py index 7dcdd3c..2d84f3d 100644 --- a/src/cso/cso_gridded.py +++ b/src/cso/cso_gridded.py @@ -46,6 +46,9 @@ # 2026-01, Arjo Segers # Fixed timerange setup. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# ######################################################################## @@ -281,6 +284,7 @@ class CSO_GriddedAverage(utopya.UtopyaRc): # tools: from . import cso_file from . import cso_mapping + from . import cso_listing # init base object: utopya.UtopyaRc.__init__(self, rcfile, rcbase=rcbase, env=env) @@ -662,7 +666,7 @@ class CSO_GriddedAverage(utopya.UtopyaRc): # new file? if listing_curr != listing_latest: # read listing: - listing = cso_file.CSO_Listing(listing_curr, indent=f"{indent} ") + listing = cso_listing.CSO_Listing(listing_curr, indent=f"{indent} ") # store name: listing_latest = listing_curr # endif diff --git a/src/cso/cso_inquire.py b/src/cso/cso_inquire.py index 772271e..e50b18c 100644 --- a/src/cso/cso_inquire.py +++ b/src/cso/cso_inquire.py @@ -37,6 +37,9 @@ # 2026-03, Arjo Segers # Added `CSO_Inquire_Plot_Avail` class. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# ######################################################################## ### @@ -581,6 +584,7 @@ class CSO_Inquire_Plot_Avail(utopya.UtopyaRc): # tools: from . import cso_file + from . import cso_listing # setup graphics back-end: matplotlib.use("Agg") @@ -643,7 +647,7 @@ class CSO_Inquire_Plot_Avail(utopya.UtopyaRc): raise Exception #endif # read: - lst = cso_file.CSO_Listing( filenames[0] ) + lst = cso_listing.CSO_Listing( filenames[0] ) # time range: t1 = lst.df["start_time"].min() diff --git a/src/cso/cso_listing.py b/src/cso/cso_listing.py new file mode 100644 index 0000000..37d2fe7 --- /dev/null +++ b/src/cso/cso_listing.py @@ -0,0 +1,604 @@ +# +# CHANGES +# +# 2023-08, Arjo Segers +# Added "method" argument to CSO_Listing.Select(). +# +# 2025-01, Arjo Segers +# Extened CSO_Listing class with selection and other methods. +# +# 2025-04, Arjo Segers +# Convert datetime time values in listing files after reading to keep default str values in table. +# Updated methods of CSO_Listing class: +# - GetValues now returns the index column by defaul. +# - added GetValue method +# Removed usage of undefined filename in error messages. +# +# 2025-04, Arjo Segers +# Extended sort options for listing files. +# Support creation of listing file objects without filename. +# Support selection of multiple records from listing file. +# +# 2025-04, Arjo Segers +# Changed imports for python packaging. +# +# 2025-09, Arjo Segers +# Improved speed of CSO_Listing.Selection method. +# +# 2025-10, Arjo Segers +# Replace ' and ' by ' & ' when selecting records. +# +# 2026-01, Arjo Segers +# Removed whitespace from template replacement. +# +# 2026-04, Arjo Segers +# Added optional "columns" argument to initialize empty listing object. +# +# 2026-05, Arjo Segers +# Support `name` argument for GetRecord method. +# + +######################################################################## +### +### help +### +######################################################################## + +""" +********************** +``cso_listing`` module +********************** + +Create and access listing files of data sets. + + +Class hierchy +============= + +The classes are defined according to the following hierchy: + +* :py:class:`.CSO_Listing` + + +Classes +======= +""" + +# modules: +import logging + + +######################################################################## +### +### listing of data files +### +######################################################################## + + +class CSO_Listing(object): + + """ + Storage for table with file properties. + The records describe properties of :py:class:`CSO_File` files, + for example the time range of the pixels included. + This could be used to quickly scan an archive for filenames + with pixels within a requested interval. + + The content of the table file looks like:: + + filename;start_time;end_time + 2007/06/RAL-IASI-CH4_20070601T003855.nc;2007-06-01 00:45:36.827000000;2007-06-01 02:23:57.567000000 + 2007/06/RAL-IASI-CH4_20070601T022359.nc;2007-06-01 02:23:59.512000000;2007-06-01 04:05:57.328000000 + ... + + Optional arguments: + + * ``filename`` : listing file read into table + * ``columns`` : extra columns that will be created if an empty listing is initialized + + """ + + def __init__(self, filename=None, columns=[], indent=""): + """ + Initialize empty table or read existing. + """ + + # modules: + import os + import pandas + + # csv style: + self.sep = ";" + # head for index column: + self.index_label = "filename" + + # store filename: + self.filename = filename + + # directory name: + if self.filename is not None: + # check .. + if not os.path.isfile(filename): + logging.error("listing file not found: %s" % filename) + raise Exception + # endif + + # base directory: + self.dirname = os.path.dirname(self.filename) + # could be empty .. + if len(self.dirname) == 0: + self.dirname = os.curdir + # endif + + # info ... + logging.info(f"{indent}read listing {filename} ...") + # read: + self.df = pandas.read_csv( + filename, + sep=self.sep, + index_col=self.index_label, + dtype="str", + ) + # convert datetime columns, as the default is "str" + # the "parse_dates" argument could not be used: + self.df["start_time"] = pandas.to_datetime(self.df["start_time"]) + self.df["end_time"] = pandas.to_datetime(self.df["end_time"]) + + else: + # not defined yet, assume current location: + self.dirname = os.curdir + + # new empty table: + self.df = pandas.DataFrame(columns=["start_time", "end_time"] + columns) + + # endif + + # enddef __init__ + + # * + + def Save(self, filename, dmode=None, indent=""): + """ + Write table to file. + """ + + # info ... + logging.info(f"{indent}save listing {filename} ...") + + # create directory if needed: + CheckDir(filename, dmode=dmode) + + # write all columns: + columns = list(self.df.keys()) + # index only once .. + if self.index_label in columns: + columns.remove(self.index_label) + + # save, also write the index column: + self.df.to_csv(filename, sep=self.sep, columns=columns, index_label=self.index_label) + + # enddef Save + + # * + + def __len__(self): + """ + Returns number of records in table. + """ + + return len(self.df) + + # enddef __len__ + + # * + + def GetRecord(self, irec=None, name=None): + """ + Return record for corresponding 0-based row index or index name. + """ + + # modules: + import numpy + + # by name? + if name is not None: + indx, = numpy.where( self.df.index == name ) + if len(indx) == 0: + logging.error(f"value `{name}` not found in index") + raise Exception + #endif + irec = indx[0] + #endif + + # extract: + rec = self.df.iloc[irec].copy() + # add index values: + rec[self.index_label] = self.df.index[irec] + + # ok + return rec + + # enddef GetRecord + + # * + + def GetValue(self, filename, key): + """ + Return ``key`` value for record of ``filename``. + """ + + # copy: + return self.df.at[filename, key] + + # enddef GetValue + + # * + + def GetValues(self, name=None): + """ + Return :py:class:`pandas.Series` object with all values for provided column ``name``, + or the filenames if ``name`` is not provided. + """ + + # switch .. + if name is None: + return self.df.index + elif name in self.df.keys(): + return self.df[name] + else: + logging.error(f"column '{name}' not found in listing") + raise Exception + # endif + + # enddef GetValues + + # * + + def __contains__(self, key): + """ + Returns True if key is in the table index. + """ + + return key in self.df.index + + # enddef __contains__ + + # * + + def keys(self): + """ + Returns list of column names. + """ + + return self.df.keys() + + # enddef keys + + # * + + def Cleanup(self, indent=""): + """ + Remove records from table if filename is not present anymore. + """ + + # modules: + import os + + # loop over records: + for fname in self.df.index: + # file missing? + if not os.path.isfile(os.path.join(self.dirname, fname)): + # info ... + logging.info(indent + 'drop record "%s", file not present anymore ..' % fname) + # remove record: + self.df.drop(fname, inplace=True) + # endif + # endfor + + # enddef Cleanup + + # * + + def Check(self, filename, indent=""): + """ + Return exception if no record for ``filename`` is present. + """ + + # modules: + import os + + # relative: + fname = os.path.relpath(filename, self.dirname) + + # check .. + if fname not in self.df.index: + logging.error(f"file '{fname}' is not a record in table: {self.filename}") + raise Exception + # endif + + # enddef Check + + # * + + def Update(self, filename, csf, tr=None, xcolumns=[], indent=""): + """ + Add record (or replace) for ``filename`` with information + from the ``csf`` :py:class:`CSO_File` object. + + The time range of pixels is read from the ``csf`` object, + unless the :py:class:`datetime.datetime` tupple ``tr=(t1,t2)`` is supplied. + + The optional ``xcolumns`` specify a list of global attributes + to be added as extra columns. + """ + + # modules: + import os + import numpy + + # timerange: + if tr is None: + # get from file: + tr = csf.GetTimeRange() + # endif + # check, in some products the time values were found to be undefined ... + if numpy.isnan(tr[0]) or numpy.isnan(tr[1]): + logging.error("found undefined time range while updating listing record") + raise Exception + # endif + # store: + self.df.at[filename, "start_time"] = tr[0] + self.df.at[filename, "end_time"] = tr[1] + + # extra: + for name in xcolumns: + # get value: + avalue = csf.GetAttr(name) + # copy: + self.df.at[filename, name] = avalue + # endfor # xcolumns + + # enddef Update + + # * + + def UpdateRecord(self, filename, data, indent=""): + """ + Add record (or replace) for ``filename`` with columns and values defined in the ``data`` dict. + """ + + # modules: + import os + + # loop over columns: + for key in data.keys(): + self.df.at[filename, key] = data[key] + # endfor + + # enddef UpdateRecord + + # * + + def Select( + self, tr=None, method="overlap", expr=None, blacklist=[], verbose=True, indent="", **kwargs + ): + """ + Return :py:class:`CSO_Listing` object with selection of records. + + Optional arguments: + + * ``tr=(t1,t2)`` : select records that could be assigned to this timerange; + requires the ``method`` keyword: + + * ``method='overlap'`` (default) : select records that have pixels overlapping with interval + * ``method='middle'`` : select pixels with the middle of ``start_time`` and ``end_time`` + within the inerval ``(t1,t2]`` + + * ``name=value`` arguments select records based on column values + + * The exprssion ``expr`` provides a list of ';'-seperated selection expressions + with templates ``%{..}`` for the column names:: + + (%{processor_version} == '020400') & (%{processing} == 'RPRO') ; ... + + This is evaluted after the previously described selections. + Eventually skip files in ``blacklist``. + + """ + + # modules: + import os + import pandas + + # init result with entire dataset: + df = self.df + + # * + + # time range specified? + if tr is not None: + # convert if needed: + t1, t2 = pandas.Timestamp(tr[0]), pandas.Timestamp(tr[1]) + # switch: + if method == "overlap": + # select: + df = df[(df["start_time"] <= t2) & (df["end_time"] >= t1)] + # + elif method == "middle": + # mid times: + mid_time = df["start_time"] + (df["end_time"] - df["start_time"]) / 2 + # select: + df = df[(mid_time > t1) & (mid_time <= t2)] + # + else: + logging.error('unsupported method "%s" with argument ``tr``' % method) + raise Exception + # endif + # endif + + # * + + # other selections: + for key, value in kwargs.items(): + # check .. + if key not in df.keys(): + logging.error(f"key '{key}' not defined in listing: {self.filename}") + raise Exception + # endif + # select: + df = df[df[key] == value] + # endfor # key/value pairs + + # * + + # evaluate selection expression? + if expr is not None: + # split into sequnece of critera; first expression with macthing record will be used: + selections = expr.split(";") + + # initialize empty list of selected record indices (filenames): + selected = [] + # storage for status label: "selected", "blacklisted", ... + filestatus = {} + # loop over selection criteria, + # this should give either none or a single file: + for selection in selections: + # replace templates: + # %{orbit} == '12345' + # to: + # df['orbit'] == '12345' + for key in self.df.keys(): + selection = selection.replace(f"%{{{key}}}", f"df['{key}']").strip() + # endfor + ## testing ... + #logging.info(f"{indent}selection `{selection}` ...") + # adhoc: replace some operators that are not allowded anymore ... + selection = selection.replace( " and ", " & " ) + # evaluate: + xdf = df[ eval(selection) ] + # any? + if len(xdf) > 0: + # selected indices (filenames): + selected = [] + for indx, xrec in xdf.iterrows(): + # skip? + if os.path.basename(indx) in blacklist: + filestatus[indx] = "blacklisted" + continue + # endif + # store: + selected.append(indx) + # set status: + filestatus[indx] = "selected" + # endif + # endfor # records + + # any selected? + if len(selected) > 0: + # leave: + break + # endif + + # endfor # selection criteria + + # show selection? + if verbose: + # info ... + logging.info(f"{indent}available record(s):") + # loop: + for fname, row in df.iterrows(): + line = fname + if fname in filestatus.keys(): + line = line + f" [{filestatus[fname]}]" + logging.info(f"{indent} {line}") + # endfor + # endif # verbose + + # no match? + if len(selected) == 0: + # info ... + logging.warning(" no match with any selection criterium ...") + ## degug: + # for selection in selections : + # logging.warning( ' %s' % selection.strip() ) + # logging.warning( ' record:' ) + # for key in rec.keys() : + # logging.warning( ' "%s" : "%s"' % (key,rec[key]) ) + # create empty dataframe as result: + df = pandas.DataFrame(columns=df.columns) + else: + # extract selected record(s): + df = df.loc[selected] + # endif + + # endif # expr defined + + # * + + # return as listing: + lst = CSO_Listing() + lst.df = df + + # ok + return lst + + # enddef Select + + # * + + def Append(self, lst, path=""): + """ + Append content of other listing. + """ + + # modules: + import os + import pandas + + # prefix paths? + if len(path) > 0: + # copy, make "filename" a normal column again: + xdf = lst.df.reset_index() + # add path: + for key in ["filename", "filename_state"]: + xdf[key] = os.path.join(path, xdf[key]) + # endfor + # add rows to current rows, use the original index: + self.df = pandas.concat((self.df, xdf.set_index(self.index_label))) + else: + # append: + if len(self.df) == 0: + self.df = lst.df + else: + self.df = pandas.concat((self.df, lst.df)) + # endif + # endif + + # enddef Append + + # * + + def Sort(self, by=None): + """ + Sort listing table by index (default, this is the "filename") or by a named column. + """ + + # sort index or values: + if by is None: + self.df.sort_index(inplace=True) + else: + self.df.sort_values(by, inplace=True) + # endif + + # endef Sort + + +# endclass CSO_Listing + + +######################################################################## +### +### end +### +######################################################################## + diff --git a/src/cso/cso_s5p.py b/src/cso/cso_s5p.py index 838867e..d9bd459 100644 --- a/src/cso/cso_s5p.py +++ b/src/cso/cso_s5p.py @@ -88,6 +88,9 @@ # 2026-05, Arjo Segers # Convert multiple files if orbit is stored in patches. Check on duplicate output files. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# ######################################################################## ### @@ -2745,6 +2748,7 @@ class CSO_S5p_Convert(utopya.UtopyaRc): import pandas # tools: + from . import cso_listing from . import cso_file from . import cso_dataspace from . import cso_pal @@ -2789,7 +2793,7 @@ class CSO_S5p_Convert(utopya.UtopyaRc): raise Exception # endif # init storage: - listing = cso_file.CSO_Listing() + listing = cso_listing.CSO_Listing() # loop: for ifile in range(len(filename__templates)): # current: @@ -2804,7 +2808,7 @@ class CSO_S5p_Convert(utopya.UtopyaRc): # expand time templates filename = t0.strftime(filename__template) # read: - lst = cso_file.CSO_Listing(filename) + lst = cso_listing.CSO_Listing(filename) # add: listing.Append(lst) # endfor # filenames @@ -3627,6 +3631,7 @@ class CSO_S5p_Listing(utopya.UtopyaRc): # tools: from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -3670,7 +3675,7 @@ class CSO_S5p_Listing(utopya.UtopyaRc): logging.info(f"{indent} base directory: %s ..." % bdir) # initiallize for (re)creation: - listing = cso_file.CSO_Listing() + listing = cso_listing.CSO_Listing() # info ... logging.info(f"{indent} cleanup records if necessary ...") @@ -3854,6 +3859,7 @@ class CSO_S5p_Download_Listing(utopya.UtopyaRc): # tools: from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -3897,7 +3903,7 @@ class CSO_S5p_Download_Listing(utopya.UtopyaRc): logging.info(f"{indent} base directory: %s ..." % bdir) # initiallize for (re)creation: - listing = cso_file.CSO_Listing(indent=f"{indent} ") + listing = cso_listing.CSO_Listing(indent=f"{indent} ") # info ... logging.info(f"{indent} cleanup records if necessary ...") diff --git a/src/cso/cso_superobs.py b/src/cso/cso_superobs.py index f33b64c..914900b 100644 --- a/src/cso/cso_superobs.py +++ b/src/cso/cso_superobs.py @@ -33,6 +33,9 @@ # 2026-01, Arjo Segers # Fixed escape characters in (document) strings. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# ######################################################################## @@ -205,6 +208,7 @@ class CSO_SuperObs(utopya.UtopyaRc): # tools: from . import cso_file from . import cso_mapping + from . import cso_listing # init base object: utopya.UtopyaRc.__init__(self, rcfile, rcbase=rcbase, env=env) @@ -306,7 +310,7 @@ class CSO_SuperObs(utopya.UtopyaRc): # info .. logging.info("input listing: %s" % listing_file) # read: - listing = cso_file.CSO_Listing(listing_file) + listing = cso_listing.CSO_Listing(listing_file) # target file: outfile_template = self.GetSetting("output.file") diff --git a/src/cso/cso_viirs.py b/src/cso/cso_viirs.py index 24137bf..de9d0bd 100644 --- a/src/cso/cso_viirs.py +++ b/src/cso/cso_viirs.py @@ -7,6 +7,9 @@ # 2025-09, Arjo Segers # Fix longitude/latitude arrays with no-data values, interpolate between surrounding locations. # +# 2026-07, Arjo Segers +# Use new `cso_listing` module. +# ######################################################################## @@ -1261,6 +1264,7 @@ class CSO_VIIRS_Convert(utopya.UtopyaRc): # tools: from . import cso_file + from . import cso_listing # info ... logging.info(f"{indent}") @@ -1287,7 +1291,7 @@ class CSO_VIIRS_Convert(utopya.UtopyaRc): # listing table: listing_file = self.GetSetting("input.listing") # read: - listing = cso_file.CSO_Listing(listing_file) + listing = cso_listing.CSO_Listing(listing_file) # info ... logging.info(f"{indent}number of files : {len(listing)}") -- GitLab From 145579fc3e9fb55474c9f5d57dc002cc52e654dd Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Wed, 22 Jul 2026 11:16:06 +0200 Subject: [PATCH 14/20] Use listing files to create plots of simulated orbits. --- config/tutorial/tutorial.rc | 49 +++--- src/cso/cso_catalogue.py | 290 +++++++++++++++++++++--------------- 2 files changed, 202 insertions(+), 137 deletions(-) diff --git a/config/tutorial/tutorial.rc b/config/tutorial/tutorial.rc index b4e689d..4bfff9c 100644 --- a/config/tutorial/tutorial.rc +++ b/config/tutorial/tutorial.rc @@ -26,7 +26,7 @@ cso.tutorial.elements : inquire \ ! ... preprocessor steps one by one: !cso.tutorial.elements : inquire -!cso.tutorial.elements : convert +!cso.tutorial.elements : convert listing !cso.tutorial.elements : listing !cso.tutorial.elements : catalogue ! @@ -260,7 +260,8 @@ my.attr.email : Your.Name@cso.org ! base location for work directories: !my.work : /work/${USER}/CSO-Tutorial -my.work : /Scratch/${USER}/CSO-Tutorial +!my.work : /Scratch/${USER}/CSO-Tutorial +my.work : ${WORK}/CSO-Tutorial !---------------------------------------------------------- @@ -298,9 +299,12 @@ VIRTUAL_ENV : ! Obtain names of all available S5p files. ! Stored as csv with processing type, processor version, filenames, etc. -! full time range: -cso.tutorial.inquire-table-dataspace.timerange.start : 2018-04-30 00:00:00 -cso.tutorial.inquire-table-dataspace.timerange.end : 2025-01-01 00:00:00 +!! full time range: +!cso.tutorial.inquire-table-dataspace.timerange.start : 2018-04-30 00:00:00 +!cso.tutorial.inquire-table-dataspace.timerange.end : 2025-01-01 00:00:00 +! testing ... +cso.tutorial.inquire-table-dataspace.timerange.start : 2018-06-01 00:00:00 +cso.tutorial.inquire-table-dataspace.timerange.end : 2018-07-01 00:00:00 ! API url: cso.tutorial.inquire-table-dataspace.url : https://stac.dataspace.copernicus.eu/v1 @@ -347,8 +351,8 @@ cso.tutorial.inquire-plot.output.file : ${my.work}/Copernicus/Cope ! renew existing files (True|False) ? -cso.tutorial.convert.renew : True -!cso.tutorial.convert.renew : False +!cso.tutorial.convert.renew : True +cso.tutorial.convert.renew : False ! time range: cso.tutorial.convert.timerange.start : ${my.timerange.start} @@ -1148,13 +1152,18 @@ cso.tutorial.sim-catalogue.timerange.start : ${my.timerange.start} cso.tutorial.sim-catalogue.timerange.end : ${my.timerange.end} cso.tutorial.sim-catalogue.timerange.step : hour -! input files: -cso.tutorial.sim-catalogue.input.data.file : ${my.work}/CSO-oper/CSO_output_%Y%m%d_%H%M_data.nc -cso.tutorial.sim-catalogue.input.state.file : ${my.work}/CSO-oper/CSO_output_%Y%m%d_%H%M_state.nc +! listing with "data" and "state" files: +cso.tutorial.sim-catalogue.input.listing : ${my.work}/CSO-oper/listing.csv -! target files, time tempates and variable name are replaced: -cso.tutorial.sim-catalogue.output.file : ${my.work}/CSO-sim-catalogue/S5p/NO2/${my.region}/%Y/%m/%d/S5p_NO2_%Y%m%d_%H%M_%{var}.png +! target directory for figures: +! - time values : average time of pixels in file +! full filenames includes "basename" of simulation files (without "_data.nc" or "_state.nc") +! and the variable name: +! /_.png +cso.tutorial.sim-catalogue.output.dir : ${my.work}/CSO-sim-catalogue/S5p/NO2/${my.region}/%Y/%m/%d +!! domain of converted data: +!cso.tutorial.sim-catalogue.domain : ${my.region.west} ${my.region.east} ${my.region.south} ${my.region.north} ! map domain used for simulations (west east south north): cso.tutorial.sim-catalogue.domain : -10 30 35 65 !!~ globe: @@ -1238,20 +1247,20 @@ cso.tutorial.sim-catalogue-index.date.newpage : True ! content type: cso.tutorial.sim-catalogue-index.date.type : table-row ! define row values: -cso.tutorial.sim-catalogue-index.date.name : time -cso.tutorial.sim-catalogue-index.date.values : Range( 0, 23, 1, '%2.2i00' ) +cso.tutorial.sim-catalogue-index.date.name : basename +cso.tutorial.sim-catalogue-index.date.values : CsvFile( '%{date[0:4]}/%{date[4:6]}/%{date[6:8]}/basenames.csv' ) ! content type: -cso.tutorial.sim-catalogue-index.date.time.type : table-col +cso.tutorial.sim-catalogue-index.date.basename.type : table-col ! define column values: -cso.tutorial.sim-catalogue-index.date.time.name : var -cso.tutorial.sim-catalogue-index.date.time.values : ${cso.tutorial.sim-catalogue.vars} +cso.tutorial.sim-catalogue-index.date.basename.name : var +cso.tutorial.sim-catalogue-index.date.basename.values : ${cso.tutorial.sim-catalogue.vars} ! content type: -cso.tutorial.sim-catalogue-index.date.time.var.type : img +cso.tutorial.sim-catalogue-index.date.basename.var.type : img ! define image: -cso.tutorial.sim-catalogue-index.date.time.var.img : %{date[0:4]}/%{date[4:6]}/%{date[6:8]}/S5p_NO2_%{date}_%{time}_%{var}.png -cso.tutorial.sim-catalogue-index.date.time.var.kwargs : height=300 +cso.tutorial.sim-catalogue-index.date.basename.var.img : %{date[0:4]}/%{date[4:6]}/%{date[6:8]}/%{basename}_%{var}.png +cso.tutorial.sim-catalogue-index.date.basename.var.kwargs : height=300 diff --git a/src/cso/cso_catalogue.py b/src/cso/cso_catalogue.py index 732adcd..e4d47d8 100644 --- a/src/cso/cso_catalogue.py +++ b/src/cso/cso_catalogue.py @@ -42,6 +42,7 @@ # # 2026-07, Arjo Segers # Use new `cso_listing` module. +# Updated `CSO_SimCatalogue` class to use lisint files created by observation operator code. # @@ -332,16 +333,16 @@ class CSO_Catalogue(CSO_CatalogueBase): CSO_CatalogueBase.__init__(self, rcfile, rcbase=rcbase, env=env) # info ... - logging.info(indent + "") - logging.info(indent + "** plot CSO data") - logging.info(indent + "") + logging.info(f"{indent}") + logging.info(f"{indent}** plot CSO data") + logging.info(f"{indent}") # time range: t1 = self.GetSetting("timerange.start", totype="datetime") t2 = self.GetSetting("timerange.end", totype="datetime") # info ... tfmt = "%Y-%m-%d %H:%M" - logging.info(indent + "timerange: [%s,%s]" % (t1.strftime(tfmt), t2.strftime(tfmt))) + logging.info(f"{indent}timerange: [{t1.strftime(tfmt)},{t2.strftime(tfmt)}]") # target files: figfile_template = self.GetSetting("output.file") @@ -364,7 +365,7 @@ class CSO_Catalogue(CSO_CatalogueBase): nrec = len(xlst) # info ... - logging.info(" records within timerange: %i" % nrec) + logging.info(f"{indent} records within timerange: {nrec}") # any data? if nrec > 0: @@ -400,18 +401,18 @@ class CSO_Catalogue(CSO_CatalogueBase): # extra map properties: bmp_kwargs = self.GetSetting("map", totype="dict", default=dict()) - # storage for csv files with orbits per subdir: + # storage for csv files with orbits per subdir, used for index pages: orbits = {} # base path, filenames are relative to listing file: inputdir = os.path.dirname(lstfile) # info .. - logging.info(indent + "loop over orbit files ...") + logging.info(f"{indent}loop over orbit files ...") # loop: for fname in xlst.GetValues(): # info .. - logging.info(indent + " %s ..." % fname) + logging.info(f"{indent} {fname} ...") # full path: filename = os.path.join(inputdir, fname) @@ -459,12 +460,12 @@ class CSO_Catalogue(CSO_CatalogueBase): # renew? if (not os.path.isfile(figfile)) or renew: # info .. - logging.info(indent + " create %s ..." % figfile) + logging.info(f"{indent} create %s ..." % figfile) # read? if orb is None: # info ... - logging.info(indent + " read ...") + logging.info(f"{indent} read ...") # init storage and read: orb = cso_file.CSO_File(filename=filename) # endif @@ -472,11 +473,8 @@ class CSO_Catalogue(CSO_CatalogueBase): # display as local standard time: taver_lst = taver + tzone_dt # annote: - title = "%s (%s) %s" % ( - taver_lst.strftime("%Y-%m-%d %H:%M"), - tzone_label, - orbit, - ) + tstamp = taver_lst.strftime("%Y-%m-%d %H:%M") + title = f"{tstamp} ({tzone_label}) {orbit}" # settings for this variable: vkey = "var.%s" % varname @@ -489,9 +487,6 @@ class CSO_Catalogue(CSO_CatalogueBase): # long name used in labels: long_name = self.GetSetting(f"var.{varname}.long_name", default=varname) - # track defined in files? - on_track = "track_longitude" in orb.ds - # switch: if ptype == "map": # style: @@ -499,6 +494,9 @@ class CSO_Catalogue(CSO_CatalogueBase): vmax = eval(self.GetSetting(f"var.{varname}.vmax", default="None")) colors = eval(self.GetSetting(f"var.{varname}.colors", default="None")) + # track defined in files? + on_track = "track_longitude" in orb.ds + # extract corner grids and values: if on_track: xx, yy, values, attrs = orb.GetTrack(sstate.ds[vname]) @@ -612,7 +610,7 @@ class CSO_Catalogue(CSO_CatalogueBase): # endfor # orbit files - # write orbit lists: + # write orbit lists, used for index pages: for figdir in orbits.keys(): # new: df = pandas.DataFrame({"orbit": orbits[figdir]}) @@ -645,24 +643,26 @@ class CSO_SimCatalogue(CSO_CatalogueBase): :align: center :alt: image catalogue - The catalogue is created in an output directory specified with:: - - .output.dir : /Scratch/CSO/catalogue - A time range for which images should be created is specified with:: ! specifiy timerange: .timerange.start : 2012-01-01 00:00 .timerange.end : 2012-12-31 23:59 - ! step is one of: hour | day | month - .timerange.step : hour The time range is used to scan for output files from the satellite observation operator. - Both a *data* file with for example footprints and observations, - as well as a *state* file with simulations is needed:: - - .input.data.file : /Scratch/model/output/CSO_output_%Y%m%d_%H%M_data.nc - .input.state.file : /Scratch/model/output/CSO_output_%Y%m%d_%H%M_state.nc + The operator has written a *listing* csv file that has a list of the output files; + specify this file using:: + + .input.listing : /Scratch/model/output/listing.csv + + Each record in the *listing* provides the name of *data* file (with information on the pixels, retrievals, kernels, etc), + the time time range for which this file has data, + and the name of the corresponding *state* file with simulations:: + + + filename ;start_time ;end_time ;filename_state + S5p_RPRO_NO2_03273_20180601T020251_data.nc;2018-06-01 3:12:53;2018-06-01 3:17:42;S5p_RPRO_NO2_03273_20180601T020251_state.nc + : Specify a list of variables to be plotted, for example the retrieved and simulated column:: @@ -677,7 +677,7 @@ class CSO_SimCatalogue(CSO_CatalogueBase): Optionally specify target units that are different from the input:: ! convert units: - .var.vcd.units : 1e15 mlc/cm2 + .var.vcd.units : umol/m2 The value range of the colorbar could be tunes using (default limits are based on data values):: @@ -693,13 +693,23 @@ class CSO_SimCatalogue(CSO_CatalogueBase): .var.vcd.long_name : retrieved vertical column density - The name of the created image files is read from:: + The figures are created in an output directory specified in the settings; + the path could contain time templates that are replaced by the average of the time of + the plotted data:: - ! target files, time tempates are replaced: - .output.file : %Y/%m/%d/S5p_RPRO_NO2_%Y%m%d_%H%M_%{var}.png + ! target directory for figures: + ! - time values : average time of pixels in file + .output.dir : /Scratch/model/catalogue/%Y/%m/%d - 2018/06/01/S5p_RPRO_NO2_03278__vcd.png - S5p_RPRO_NO2_03278__qa_value.png + The full path the created image files is formed from the basename of the *data* (and *state*) files + plus the variable name:: + + /Scratch/model/catalogue/2018/06/01/S5p_RPRO_NO2_03273_20180601T020251_yr.png + S5p_RPRO_NO2_03273_20180601T020251_ys.png + basenames.csv + + The output direcotry also contains a file ``basenames.csv`` with the basenames of the image files; + this is used to create index pages to browse through the figures. Enable the following flag to re-create existing files, by default only non-existing files are created:: @@ -733,6 +743,7 @@ class CSO_SimCatalogue(CSO_CatalogueBase): import matplotlib.pyplot as plt # tools: + from . import cso_listing from . import cso_file from . import cso_plot @@ -743,110 +754,147 @@ class CSO_SimCatalogue(CSO_CatalogueBase): CSO_CatalogueBase.__init__(self, rcfile, rcbase=rcbase, env=env) # info ... - logging.info(indent + "") - logging.info(indent + "** plot CSO simulations") - logging.info(indent + "") + logging.info(f"{indent}") + logging.info(f"{indent}** plot CSO simulations") + logging.info(f"{indent}") # time range: t1 = self.GetSetting("timerange.start", totype="datetime") t2 = self.GetSetting("timerange.end", totype="datetime") - tstep = self.GetSetting("timerange.step") # info ... tfmt = "%Y-%m-%d %H:%M" - logging.info( - indent + "timerange: [%s,%s] by %s" % (t1.strftime(tfmt), t2.strftime(tfmt), tstep) - ) - - # target file template: - figfile_template = self.GetSetting("output.file") - - # template for data files: - datafile_template = self.GetSetting("input.data.file") - statefile_template = self.GetSetting("input.state.file") + logging.info(f"{indent}timerange: [{t1.strftime(tfmt)},{t2.strftime(tfmt)}]") - # variables to be plotted: - varnames = self.GetSetting("vars").split() + # listing file: + lstfile = self.GetSetting("input.listing") + # info ... + logging.info("listing file:") + logging.info(" %s" % lstfile) + # should exist ... + if not os.path.isfile(lstfile): + logging.error("file not found: %s" % lstfile) + raise Exception + # endif + # read listing: + lst = cso_listing.CSO_Listing(lstfile, indent=indent) + # subset of files in timerange: + xlst = lst.Select(tr=(t1, t2)) + # number of records: + nrec = len(xlst) + # info ... + logging.info(f"{indent} records within timerange: {nrec}") - # renew existing figures? - renew = self.GetSetting("renew", totype="bool") + # base path, filenames are relative to listing file: + inputdir = os.path.dirname(lstfile) - # plot domain: - domain = list(map(float, self.GetSetting("domain").split())) + # any data? + if nrec > 0: - # figure size: - figsize = eval(self.GetSetting("figsize")) + # variables to be plotted: + varnames = self.GetSetting("vars").split() - # no-data color: - color_nan = self.GetSetting("color_nan", default="0.90") - # extra map properties: - bmp_kwargs = self.GetSetting("map", totype="dict", default=dict()) + # renew existing figures? + renew = self.GetSetting("renew", totype="bool") - # info .. - logging.info(indent + "time loop ...") - # time loop: - t = t1 - while t <= t2: - # info .. - logging.info(indent + " %s ..." % t.strftime("%Y-%m-%d %H:%M")) + # plot domain: + domain = list(map(float, self.GetSetting("domain").split())) - # timestep: - if tstep in ["hour", "hourly"]: - long_timestamp = t.strftime("%Y-%m-%d %H:%M") - dt = datetime.timedelta(0, 3600) - elif tstep in ["day", "daily"]: - long_timestamp = t.strftime("%Y-%m-%d") - dt = datetime.timedelta(1) - elif tstep in ["month", "monthly"]: - t = datetime.datetime(t.year, t.month, 1, 0, 0) - long_timestamp = t.strftime("%Y-%m") - weekday, nday = calendar.monthrange(t.year, t.month) - dt = datetime.timedelta(nday) + # mid of domain, used for local time: + domain_lon_mid = 0.5 * (domain[0] + domain[1]) + # time zone shift based on bands of 360/24=15 deg + tzone_dhour = int(numpy.round(domain_lon_mid / 15.0)) + # increment: + tzone_dt = pandas.Timedelta(value=tzone_dhour, unit="hour") + # annote: + if tzone_dhour < 0: + tzone_label = "UTC%i" % tzone_dhour + elif tzone_dhour > 0: + tzone_label = "UTC+%i" % tzone_dhour else: - logging.error('unsupported time step "%s"' % tstep) - raise Exception + tzone_label = "UTC" # endif - # input files: - datafile = t.strftime(datafile_template) - statefile = t.strftime(statefile_template) + # figure size: + figsize = eval(self.GetSetting("figsize")) + + # no-data color: + color_nan = self.GetSetting("color_nan", default="0.90") + # extra map properties: + bmp_kwargs = self.GetSetting("map", totype="dict", default=dict()) + + # target file template: + figdir_template = self.GetSetting("output.dir") - # state file present? - if os.path.isfile(statefile): - # info ... - logging.info(indent + " found state file: %s" % statefile) + # storage for csv files with basenames per subdir, used for index pages: + basenames = {} + + # info .. + logging.info(f"{indent}loop over orbit files ...") + # loop: + for fname in xlst.GetValues(): + # info .. + logging.info(f"{indent} {fname} ...") + # full path: + datafile = os.path.join(inputdir, fname) + # corresponding state file: + statefile = os.path.join(inputdir, xlst.GetValue(fname,"filename_state")) # check ... - if not os.path.isfile(datafile): - logging.error("data file not found: %s" % datafile) + if not os.path.isfile(statefile): + logging.error(f"state file not found: {statefile}") raise Exception # endif - # no data file read yet .. + # basename of file without extension: + basename,ext = os.path.splitext( os.path.basename(datafile) ) + basename = basename.replace("_data","") + + # time range: + tr1 = xlst.GetValue(fname, "start_time") + tr2 = xlst.GetValue(fname, "end_time") + # mid time: + taver = tr1 + 0.5 * (tr2 - tr1) + + # target directory for figures: + figdir = taver.strftime( figdir_template ) + + # update storage for basenames: + if figdir not in basenames.keys(): + basenames[figdir] = [] + #endif + # store current: + basenames[figdir].append(basename) + + # no data files read yet .. sdata = None + sstate = None # loop over variables ... for varname in varnames: - # target file: - figfile = t.strftime(figfile_template) - figfile = figfile.replace("%{var}", varname) + + # full path to target file: + figfile = os.path.join( figdir, f"{basename}_{varname}.png" ) # renew? if (not os.path.isfile(figfile)) or renew: # info .. - logging.info(indent + " create %s ..." % figfile) + logging.info(f"{indent} create {figfile} ...") # read? if sdata is None: # info ... - logging.info(indent + " read ...") + logging.info(f"{indent} read ...") # init storage and read: sdata = cso_file.CSO_File(filename=datafile) # init storage and read: sstate = cso_file.CSO_File(filename=statefile) # endif + # display as local standard time: + taver_lst = taver + tzone_dt # annote: - title = long_timestamp + tstamp = taver_lst.strftime("%Y-%m-%d %H:%M") + title = f"{tstamp} ({tzone_label})" # settings for this variable: vkey = "var.%s" % varname @@ -874,8 +922,10 @@ class CSO_SimCatalogue(CSO_CatalogueBase): else: vfile, vname = "data", vsource # endif + # track defined in files? on_track = "track_longitude" in sdata.ds + # extract corner grids and values: if vfile == "state": # check ... @@ -937,7 +987,7 @@ class CSO_SimCatalogue(CSO_CatalogueBase): bmp.update(bmp_kwargs) # info ... - logging.info(indent + " plot pattern ...") + logging.info(f"{indent} plot pattern ...") # create map figure: if on_track: fig = cso_plot.QuickPat( @@ -981,22 +1031,28 @@ class CSO_SimCatalogue(CSO_CatalogueBase): else: # info .. - logging.info(indent + " keep %s ..." % figfile) + logging.info(f"{indent} keep %s ..." % figfile) # endif # renew # endfor # variables - ## testing ... - # break - - else: - logging.info(f"{indent} no state file: {statefile}") - # endif # state file present + # endfor # records in listing file - # next: - t = t + dt + # info .. + logging.info(f"{indent}write basename tables ..") + # write basename lists, used for index pages: + for figdir in basenames.keys(): + # new: + df = pandas.DataFrame({"basename": basenames[figdir]}) + # target file: + fname = os.path.join(figdir, "basenames.csv") + # info .. + logging.info(f"{indent} {fname} ..") + # write: + df.to_csv(fname, index=False, sep=";") + # endfor - # endfor # time loop + # endif # nrec > 0 # enddef __init__ @@ -1123,9 +1179,9 @@ class CSO_GriddedCatalogue(CSO_CatalogueBase): CSO_CatalogueBase.__init__(self, rcfile, rcbase=rcbase, env=env) # info ... - logging.info(indent + "") - logging.info(indent + "** plot CSO simulations") - logging.info(indent + "") + logging.info(f"{indent}") + logging.info(f"{indent}** plot CSO simulations") + logging.info(f"{indent}") # time range: t1 = self.GetSetting("timerange.start", totype="datetime") @@ -1134,7 +1190,7 @@ class CSO_GriddedCatalogue(CSO_CatalogueBase): # info ... tfmt = "%Y-%m-%d %H:%M" logging.info( - indent + "timerange: [%s,%s] by %s" % (t1.strftime(tfmt), t2.strftime(tfmt), tstep) + f"{indent}timerange: [%s,%s] by %s" % (t1.strftime(tfmt), t2.strftime(tfmt), tstep) ) # target file template: @@ -1154,12 +1210,12 @@ class CSO_GriddedCatalogue(CSO_CatalogueBase): bmp_kwargs = self.GetSetting("map", totype="dict", default=dict()) # info .. - logging.info(indent + "time loop ...") + logging.info(f"{indent}time loop ...") # time loop: t = t1 while t <= t2: # info .. - logging.info(indent + " %s ..." % t.strftime("%Y-%m-%d %H:%M")) + logging.info(f"{indent} %s ..." % t.strftime("%Y-%m-%d %H:%M")) # timestep: if tstep in ["hour", "hourly"]: @@ -1191,7 +1247,7 @@ class CSO_GriddedCatalogue(CSO_CatalogueBase): # renew? if (not os.path.isfile(figfile)) or renew: # info .. - logging.info(indent + " create %s ..." % figfile) + logging.info(f"{indent} create %s ..." % figfile) # input file(s): if "-" in varname: -- GitLab From aa1a1e23737ef9d37075c5b1749dbc54541e869c Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Wed, 22 Jul 2026 11:16:34 +0200 Subject: [PATCH 15/20] Updated introduction text. --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index 5f3c04f..ec7b9ce 100644 --- a/README.md +++ b/README.md @@ -121,6 +121,9 @@ or creates it if not present yet: # install local codes: pip install --editable . + # optionally also install dependencies to generate documenation locally: + pip install --editable .[docs] + fi -- GitLab From bef61a3e5dc64885d4cc8f24869f2cb2acdd2e7e Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Wed, 22 Jul 2026 11:18:13 +0200 Subject: [PATCH 16/20] Updated cleanup action. --- oper/Makefile | 1 + 1 file changed, 1 insertion(+) diff --git a/oper/Makefile b/oper/Makefile index 25584bb..b838b28 100644 --- a/oper/Makefile +++ b/oper/Makefile @@ -25,6 +25,7 @@ clean: rm -f tutorial*.log* rm -f cso.log* rm -f *.nc + rm -f listing.csv clean-all: clean (cd src; make clean) -- GitLab From f8eaa16dd928c6cce2689270666bca971f24faa0 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Wed, 22 Jul 2026 11:18:31 +0200 Subject: [PATCH 17/20] Added debug templates. --- oper/src/cso_sat.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/oper/src/cso_sat.F90 b/oper/src/cso_sat.F90 index 50a9c75..2b5a8f1 100644 --- a/oper/src/cso_sat.F90 +++ b/oper/src/cso_sat.F90 @@ -1412,7 +1412,7 @@ contains data0=tvalues, units=tunits ) IF_NOT_OK_RETURN(status=1) - ! convert time interval bounds to values: + ! convert time interval bounds to values with similar units: if ( t1 < t2 ) then ! "forward" with t1 Date: Wed, 22 Jul 2026 11:22:11 +0200 Subject: [PATCH 18/20] Updated template settings for obseration operator. --- oper/tutorial_oper_S5p.rc | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/oper/tutorial_oper_S5p.rc b/oper/tutorial_oper_S5p.rc index a194197..4aefbc5 100644 --- a/oper/tutorial_oper_S5p.rc +++ b/oper/tutorial_oper_S5p.rc @@ -44,11 +44,7 @@ tutorial.mapping.levels : 5 !------------------------------------------------- ! template for listing with converted files: -tutorial.S5p.no2.listing : ../CSO-data/CAMS/S5p/NO2/C03/listing.csv -!!~ testing global orbits: -!tutorial.S5p.no2.listing : ../CSO-data/globe/S5p/NO2/C03/listing.csv -!!~ testing super observations: -!tutorial.S5p.no2.listing : ../CSO-data/globe/S5p/NO2/C03__sup1x1__qa08/listing.csv +tutorial.S5p.no2.listing : /work/CSO-Tutorial/CSO-data/CAMS/S5p/NO2/C03/listing.csv ! also read info on original track (T|F)? ! if enabled, this will be stored in the output too: @@ -89,7 +85,7 @@ tutorial.S5p.no2.dvar.yr.source : vcd ! retrieval error covariance: !~ dimensions, copied from data file: -tutorial.S5p.no2.dvar.vr.dims : retr retr +tutorial.S5p.no2.dvar.vr.dims : retr retr0 !~ source variable: tutorial.S5p.no2.dvar.vr.source : vcd_errvar -- GitLab From d60228419a880157f6b562fd3f14a6e3cc62dd87 Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Wed, 22 Jul 2026 11:23:06 +0200 Subject: [PATCH 19/20] Updated change log. --- CHANGELOG | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/CHANGELOG b/CHANGELOG index 41efa1d..1354e31 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -605,5 +605,13 @@ Convert multiple files if orbit is stored in patches. Check on duplicate output vx.y ---- -Updated Fortran operator. +Updated observation operator to use time selections of pixels in data files. oper/ + +Created seperatue module `cso_listing` for the listing objects. + src/cso/cso*.py + doc/source/ + +Use listing files created by observation operator. + src/cso/cso_catalogue.py + -- GitLab From da871cbc9e1810ffb3062e36d2dfceb52e244c5c Mon Sep 17 00:00:00 2001 From: Arjo Segers Date: Wed, 22 Jul 2026 11:32:38 +0200 Subject: [PATCH 20/20] Defined version 2.18. --- CHANGELOG | 4 ++-- doc/source/history.rst | 5 +++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 1354e31..d2cd83e 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -602,8 +602,8 @@ Convert multiple files if orbit is stored in patches. Check on duplicate output config/tutorial/tutorial.rc -vx.y ----- +v2.18 +----- Updated observation operator to use time selections of pixels in data files. oper/ diff --git a/doc/source/history.rst b/doc/source/history.rst index 994a0ba..6292af9 100644 --- a/doc/source/history.rst +++ b/doc/source/history.rst @@ -138,6 +138,11 @@ A summary of the versions and changes. | Convert multiple files if orbit is stored in patches. Check on duplicate output files. | *(2026-05)* +* | *v2.18* + | Updated oservation operator code: select multiple data files that overlap with simulation time, use time selection for pixels. + | *(2026-07)* + + To be done ========== -- GitLab