open TELEMAC-MASCARET The mathematically superior suite of solvers

Hello community. I am developing a river sediment transport model within which also exist urban areas, so as you can imagine it is necessary configure the bottom with two types of behavior; first the river zone with erodible bottom and second the urban area bottom non erodible.

For this purpose it is well known that modifying the user_bed_init.f subroutine we can set erodible and non erodible zones with ESTRATUM thickness. Reading old posts and with own expertise I have realized that the best way is not using a selafin file with info of erodible and non erodible zones.

The tested way is assigning the node ID in ESTRATUM vector. But in this particular case we are talking of hundred of thousands of nodes. Under this view I am trying to optimize this setting using the user_bed_init.f subroutine.

First with the help of a selafin file read in python I got two vector with node ID's, one with info of the erodible zone; and the other one with the info of the non erodible zone. This two vector were saved in separated txt files (ero_zones.txt; non_ero_zones.txt)

After that I modified the fortran subroutine to read both txt files, and with the ID info of each zone, assign in the respective ESTRATUM the thickness of 5 meters for erodible zone, and 0 meters for non erodible zone.

All three files (user_bed_init.f; ero_zones.txt; non_ero_zones.txt) in a first try were into the USER_FORTRAN folder. During the model running appeared an error noting that there are no ero_zones.o and non_ero_zones.o files in the running folder.

After that I take out both txt files from USER_FORTRAN folder, but keeping those in the working folder. With this modification the user_bed_init.f file was compiled but apparently there is no change in the model bottom.

Is this idea possible? What was the mistake in the first try? If the subroutine is compiled with instructions of reading txt files, why it doesn´t work?

Best regards.

Hi,
I am a big friend of defining rigid bed in the selafin file.
In your case you have to ensure that your geometry file and your txt files with the IDs are fit together...

we use the following procedure:
1) definition of a "RIGID BED" variable in the gaia geometry file
2) reading this variable in temporary variable T1 in user_bed_init.f
...
USE DECLARATIONS_GAIA, ONLY: AVA0,SED_THICK,GAI_FILES,GAIGEO,T1,
& MESH,ZF
...
INTEGER IERR, RECORD, I
DOUBLE PRECISION BID
...
RECORD = 0
BID = 0.D0
CALL FIND_VARIABLE(GAI_FILES(GAIGEO)%FMT,GAI_FILES(GAIGEO)%LU,
& 'RIGID BED ',T1%R, MESH%NPOIN,IERR,
& RECORD=RECORD,TIME_RECORD=BID)
3) calculating bed thickness by bottom - rigid bed in user_bed_init.f
ESTRATUM(ISTRAT,IPOIN) = MAX(ZF%R(IPOIN)-T1%R(IPOIN), 0.D0)

This procedure works fine for our applications.

If you want to have IDs instead of the rigid bed, you can save them in the selafin file. Of course the calculation of the ESTRATUM would be different then.

Best regards,
Rebekka

Thank you so much Rebekka. In fact it looks a kind difficult.

One question, with the method that you have described, it is possible to define both bottom characteristics, i mean a erodible zone and a non erodible zone?

Could you show me a full user_bed_init.f file that yo have used?

Is necessary to define in the cas file something related with the RIGID BED info form the selafin file?

Best regards

A selafin-file can have several variables like bottom or rigid bed and it has a value for each grid node. If the values of bottom and rigid bed are equal at one node it means there is rigid bed = non-erodable node. If the bottom is higher than the rigid bed you have an erodable node. So you can define both.
I hope this answers your question.

You need a preprocessor which can define the rigid bed (or an ID) for each point. And there it is clever to use polygons=zones.

There is nothing to do in the cas file.

I have added you a full user_bed_init.f

Best regards,
Rebekka

Thank you so much. Please could you help me with more info about de BED RIGID data that you use to define erodible and non erodible zones?.

Because I have a doubt with the thickness calculation ESTRATUM(ISTRAT,IPOIN) = MAX(ZF%R(IPOIN)-T1%R(IPOIN), 0.D0). Because under my understanding or misunderstanding I see two possibles scenarios:

1. Erodible zone. Example: ZF%R(IPOIN) = 2800; T1%R(IPOIN) = 2795 (it is 5 meter of erodible zone?) MAX(5,0)=5 -- Erodible zone
2. Non erodible zone. Example: ZF%R(IPOIN) = 2750; T1%R(IPOIN) = 2750 (it is 0 meter of erodible zone?) MAX(0,0)=0 -- Non erodible zone

So in this case the BED RIGID is related with the topography?

Best regards
Ricardo

Hi Ricardo,
yes the rigid bed is related to the topography. That is how we do it. You can of course directly read a thickness or an ID. Gaia needs the thickness. So if you would put a thickness in the geometry file you would later set the thickness variable of gaia:
ESTRATUM(ISTRAT,IPOIN) = T1%R(IPOIN)

Best regards,
Rebekka

Thank you so much Rebekka, I will try with your idea and I will bring with the news.

Best regards
Ricardo