Skip to content

Fm/feature/unst 8910 nesting read 3d vectors3#1019

Open
mihai-pricope wants to merge 32 commits into
mainfrom
fm/feature/UNST-8910_nesting_read_3d_vectors3
Open

Fm/feature/unst 8910 nesting read 3d vectors3#1019
mihai-pricope wants to merge 32 commits into
mainfrom
fm/feature/UNST-8910_nesting_read_3d_vectors3

Conversation

@mihai-pricope

Copy link
Copy Markdown
Contributor

What was done

Read quantities with time interpolated Z-Dimension from _his.nc files (_his.nc files are produced by dflowfm, so this will enable nesting).

Evidence of the work done

  • Video/figures
    <add video/figures if applicable>
  • Clear from the issue description
  • Not applicable

Tests

  • Tests updated
    <add testcase numbers if applicable, Issue number>
  • Not applicable

Documentation

  • Documentation updated
    <add description of changes if applicable, Issue number>
  • Not applicable

Issue link

https://issuetracker.deltares.nl/browse/UNST-8910

@arthurvd arthurvd left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi Mihai, a few small remarks as well as a few high level questions. Have a look at them, and we can have a chat to discuss some of them if you want.


do idIdx = 1, size(sourceItemIds)
sourceItemId = sourceItemIds(idIdx)
success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

catch error, such that subsequent iterations cannot hide any errors:

Suggested change
success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)
success = ecAddConnectionSourceItem(ecInstancePtr, connectionId, sourceItemId)
if (.not. success) then
exit
end if

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

integer :: ncid !< unique NetCDF ncid
character(len=maxFileNameLen) :: ncfilename !< netCDF filename
integer, allocatable, dimension(:) :: dimlen !< lengths of dimensions
integer, allocatable, dimension(:) :: variable_dimension !< lengths of dimensions

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rename to variable_ndims (recognized NetCDF concept), or alternatively variable_rank.

(Current name is ambiguous whether rank/shape/dimension ids are stored.)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

integer, intent(in) :: l_id
integer, dimension(:), intent(in) :: dims
integer, intent(in) :: timelevel
integer, intent(in) :: func

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Move this to last dummy arg declaration line (should be in order of appearance).
And nice to have docstring (should be a BC_FUNC_* enum value that specifies the type of timeseries.)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

if (.not. ecItemSetType(instancePtr, zTargetItemId, itemPT%accessType)) then
return
end if
if (.not. ecItemSetQuantity(instancePtr, zTargetItemId, itemPT%quantityPtr%id)) then

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does the zTargetItem get the same id as the actual data item (e.g. 'salinity'). Is this required for correct workings? (I'd expect each item to have a descriptive id, so something describe the vertical z coordinate_.

Is bcBlockPtr%ncptr%variable_names(bcBlockPtr%ncptr%vertical_coordinate_id) more appropriate?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I made some modification but it seems that if the item is not called polytim_item, is not gonna work.
zQuantityId = ecInstanceCreateQuantity(instancePtr)
if (.not. (ecQuantitySet(instancePtr, zQuantityId, name='polytim_item'))) then
return
end if
need to debug further to figure out why

Comment on lines +1890 to +1907
! ToDo, more generic approach to determine veriabel name from type of boundary
call str_lower(quantityname)
if (index(trim(bctfilename)//'|', '_his.nc|') > 0) then
! History file
if (strcmpi(quantityname,'waterlevelbnd' )) quantityname = 'waterlevel'
if (strcmpi(quantityname,'salinitybnd' )) quantityname = 'salinity'
if (strcmpi(quantityname,'temperaturebnd' )) quantityname = 'temperature'
if (strcmpi(quantityname,'uxuyadvectionvelocitybnd')) quantityname = 'x_velocity'
else
! Old existing nc files
if (strcmpi(quantityname,'waterlevelbnd' )) quantityname = 'waterlevelbnd'
if (strcmpi(quantityname,'salinitybnd' )) quantityname = 'so'
if (strcmpi(quantityname,'temperaturebnd' )) quantityname = 'thetao'
if (strcmpi(quantityname,'uxuyadvectionvelocitybnd')) quantityname = 'ux'
! Old existing nc files (temporary fix for Tom, durban)
! if (strcmpi(quantityname,'so' )) quantityname = 'salinity'
! if (strcmpi(quantityname,'thetao' )) quantityname = 'temperature'
end if

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you not use Thomas's utility function here? --> m_ec_support::ecSupportNetcdfGetQuantityCandidateNames()

  • extend the waterlevelbnd case to include a fallback name
  • add salinity and other bnds
  • Where needed, call ecProviderSearchStdOrVarnames(fileReaderPtr, i, idvar, ncstdnames, ncvarnames, uservarnames=nccustomnames) and same for fallback names.

(I'm ok if you create a separate issue for this to do it in a future PR.)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added this. https://issuetracker.deltares.nl/browse/UNST-10104
I have also a few ideas around it, but I think at this point this can easily scoped out of this PR

Comment on lines +1227 to 1230
! Only interpolate if both T0 and T1 contain at least 1 valid value, else set valuesT to dmiss
if (all(valuesT0 == dmiss) .or. all(valuesT1 == dmiss)) then
valuesT = dmiss
else

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This does not handle the situation where most values are not dmiss, but on a particular index i both valuesT0 and valuesT1 have a dmiss. So then we fall through to the do-loop on line 1246, which starts to do math on line 1255. Should you not catch the double dmiss there (at the index i level).

Comment on lines +1718 to +1719
! TK_Temp: Deal with one sided interpolation, set kL or kR to 0 if arr1D does not contain valid values
! TK_Temp: Left side

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
! TK_Temp: Deal with one sided interpolation, set kL or kR to 0 if arr1D does not contain valid values
! TK_Temp: Left side
! Deal with one sided interpolation, set kL or kR to 0 if arr1D does not contain valid values at the current time.
! Left side:

Comment on lines 1749 to 1750
select case (connection%converterPtr%operandType)
case (operand_replace_element, operand_replace, operand_add)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Warning: you will probably encounter merge conflicts here with changes from Thomas from PR #959. Try to adopt his new apply_operand() approach in your code fragment here.

if (kL > 0) then
kbeginL = vectormax * maxlay_src * (kL - 1) + 1! refers to source right column
kendL = kbeginL + vectormax*maxlay_src - 1
if (all(connection%sourceItemsPtr(1)%ptr%targetFieldPtr%arr1Dptr(kbeginL:kendL) == dmiss)) then

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The code here (and below too for kR) enables one-sided interpolation if ALL layer values are dmiss in the kL point. But what if only some layers contain dmiss on point kL: will those dmiss values end up in a linear interpolation below, or not?

block
integer :: zItemId

if (present(zTargetItemId) .and. present(zFileReaderId)) then

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

High-level question about the code below: here you're creating a subconverter + connection for the extra z quantity. How does that relate to the new code in ec_addtimespacerelation(), where also some extra association is made between multiple source items and a single target item (code over there changed ecFindItemInFileReader() --> ecFindItemsInFileReader() and runs over do idIdx = 1, size(sourceItemIds).

@mihai-pricope mihai-pricope force-pushed the fm/feature/UNST-8910_nesting_read_3d_vectors3 branch from 70aa590 to f4139ea Compare July 3, 2026 06:37
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants