TNO Intern

Commit a8fc1b8b authored by Arjo Segers's avatar Arjo Segers
Browse files

Introduced `Series` object to open multiple files from a time series for a... 

Introduced `Series` object to open multiple files from a time series for a single interval. Pixel selection based on time values.
parent a851bbb3

oper/src/Makefile_deps

+22 −21
Original line number Diff line number Diff line 
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.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_string.o cso_file.o cso_datetime.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_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_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_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_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_S5p.o : tutorial_oper_S5p.F90 cso.inc cso.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

oper/src/cso.F90

+5 −0
Original line number Diff line number Diff line
@@ -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

oper/src/cso_sat.F90

+249 −458

File changed.

Preview size limit exceeded, changes collapsed.

oper/src/tutorial_oper_S5p.F90

+309 −291
Original line number Diff line number Diff line
@@ -12,6 +12,9 @@ 
! 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,14 +116,16 @@ 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
@@ -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,28 +457,33 @@ 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)
@@ -501,24 +519,60 @@ program Tutorial_Oper_S5p 
        call sdata%Read( rcF, rcbase, status )

        IF_NOT_OK_RETURN(status=1)



      ! ~ simulation

        ! number of local pixels read:

        call sdata%Get( status, npix=npix )

        IF_NOT_OK_RETURN(status=1)

        ! any local pixels?

        if ( npix > 0 ) then



      ! 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 )

          ! 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)



      ! any pixels at all?

      if ( npix_all > 0 ) then

          ! 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)



        ! initialize simulation state;

        ! optional arguments:

        !   description='long name'    : used for output attributes

        call sstate%Init( sdata, rcF, rcbase, status, &

                             description='simulated retrievals' )

          ! 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)



      ! get pointer to state object:

      call sSeries%Get( dSeries, status, iset=iset, sstate=sstate )

      IF_NOT_OK_RETURN(status=1)



      ! new state?

      if ( isnew ) then



        ! number of user defined dimensions:

        call sstate%Get( status, nudim=nudim )

        IF_NOT_OK_RETURN(status=1)
@@ -550,30 +604,27 @@ program Tutorial_Oper_S5p 
        call sstate%EndDef( status )

        IF_NOT_OK_RETURN(status=1)



        ! any local pixels?

        if ( npix > 0 ) then

      end if ! isnew



          ! pointers to (*,pixel) arrays:

          ! - footprint centers  (not used, for inspiration ..)

          ! - footprint corners

          call sdata%Get( status, lons=lons, lats=lats, clons=clons, clats=clats )

      ! ~ simulation



      ! number of local pixels read:

      call sdata%Get( status, npix=npix )

      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 )

      ! 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



          ! 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 )

      ! any pixels selected?

      if ( ntsel > 0 ) then



        ! 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)



        ! number of user defined variables:
@@ -609,8 +660,10 @@ program Tutorial_Oper_S5p 
                ! 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

                ! loop over time selection:

                do itsel = 1, ntsel

                  ! pixel index:

                  ipix = ipix__tsel(itsel)

                  ! any source contributions?

                  if ( nw(ipix) > 0 ) then

                    ! init sum:
@@ -621,7 +674,7 @@ program Tutorial_Oper_S5p 
                      data1(:,ipix) = data1(:,ipix) + conc(ii(iw),jj(iw),:) * ww(iw)/areas(ipix)

                    end do ! iw

                  end if ! nw > 0

                  end do ! ipix

                end do ! itsel



              !~ model half-level pressures:

              case ( 'mod_hp' )
@@ -634,8 +687,10 @@ program Tutorial_Oper_S5p 
                ! 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

                ! loop over time selection:

                do itsel = 1, ntsel

                  ! pixel index:

                  ipix = ipix__tsel(itsel)

                  ! any source contributions?

                  if ( nw(ipix) > 0 ) then

                    ! init sum:
@@ -646,7 +701,7 @@ program Tutorial_Oper_S5p 
                      data1(:,ipix) = data1(:,ipix) + hp(ii(iw),jj(iw),:) * ww(iw)/areas(ipix)

                    end do ! iw

                  end if ! nw > 0

                  end do ! ipix

                end do ! itsel



              !~ cloud cover: total:

              case ( 'mod_tcc' )
@@ -659,8 +714,10 @@ program Tutorial_Oper_S5p 
                ! 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

                ! loop over time selection:

                do itsel = 1, ntsel

                  ! pixel index:

                  ipix = ipix__tsel(itsel)

                  ! any source contributions?

                  if ( nw(ipix) > 0 ) then

                    ! init sum:
@@ -671,7 +728,7 @@ program Tutorial_Oper_S5p 
                      data0(ipix) = data0(ipix) + tcc(ii(iw),jj(iw)) * ww(iw)/areas(ipix)

                    end do ! iw

                  end if ! nw > 0

                  end do ! ipix

                end do ! itsel



              !~ cloud cover: profile:

              case ( 'mod_cc' )
@@ -684,8 +741,10 @@ program Tutorial_Oper_S5p 
                ! 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

                ! loop over time selection:

                do itsel = 1, ntsel

                  ! pixel index:

                  ipix = ipix__tsel(itsel)

                  ! any source contributions?

                  if ( nw(ipix) > 0 ) then

                    ! init sum:
@@ -696,7 +755,7 @@ program Tutorial_Oper_S5p 
                      data1(:,ipix) = data1(:,ipix) + cc(ii(iw),jj(iw),:) * ww(iw)/areas(ipix)

                    end do ! iw

                  end if ! nw > 0

                  end do ! ipix

                end do ! itsel



              !~

              case default
@@ -717,60 +776,9 @@ program Tutorial_Oper_S5p 


        end if ! nuvar > 0



        end if ! npix > 0

      end if ! ntsel > 0



        ! ~ formula



        ! 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)



        ! apply formula: kernel convolution etc:

        call sstate%ApplyFormulas( sdata, status )

        IF_NOT_OK_RETURN(status=1)



        ! ~ put out



        ! 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 )

        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 )

        IF_NOT_OK_RETURN(status=1)



        ! ~ clear



        ! clear:

        nullify( glb_lon )

        nullify( glb_lat )

        nullify( glb_select )



        ! done with state:

        call sstate%Done( sdata, status )

        IF_NOT_OK_RETURN(status=1)



      end if ! npix_all > 0



      ! done with data:

      call sdata%Done( status )

      IF_NOT_OK_RETURN(status=1)



      !! testing ...

      !write (logfu,'("break after first orbit")')

      !exit



    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: