open TELEMAC-MASCARET The mathematically superior suite of solvers

Hi there,
I apologise for posting this thread under TELEMAC-3D but it kind of relates it to.
I have been trying to incorporate TELEMAC-3D variables into the water quality file of SUBIEF. However, when compiling the module using the command wq2subief it ignores TELEMAC-3D variables as it seems not to recognise them.
Any hints on how to include temperature and salinity in a subief water quality model?

Thanks alot

Zeinab

Hi
Subief is no longer maintained so a lot of problems could exist with old version.
Depending your interest, you could try to include the water quality equations directly in Telemac-3D

Regards

Thank you very much for your response.
I would like to try to model the water quality varaible using TELEMAC-3D, only problem is that I don't know how. I haven't come across a single example that demonstrates how to do that!!!
Any clues, direction or perhaps an example will be very much appreciated.

Thanks alot

Zeinab

Hi again,
What I am trying to do is to model the decay of e.coli as a function of temperature, salinity and solar radiation using SUBIEF-3D (I managed to do a constant decay with it) but I don't mind doing that in TELEMAC-3D if this doesn't involve too much hassle. So maybe if I was given directions as to how to accomplish this I will be very grateful.
Best Regards,
Zeinab

Hello Zeinab,

To code your own decay rate for a tracer, you can implement SOURCE_TRAC subroutine in TELEMAC-3D.

The arrays you can play with are S1TA and S0TA.
- S0TA acts explicitely on the tracer value at the begining of the time step as follows:
TA = TA + S0TA ... + ADVECTION+DIFFUSION +...
- S1TA acts implicitely on the tracer value at the end of the time step as follows:
TA *( 1.-S1TA ) = TA + ... + ADVECTION+DIFFUSION

Please note the sign in front of S0TA and S1TA, and that S1TA is a multiplier of TA, the tracer concentration.

So, say you have
INTEGER IPOIN3, IPOIN2, IPLAN
DOUBLE PRECISION DECAYDAY,DECAYRATE
Here, for instance, the formulation use a T90 ( Loss of 'DECAYRATE' % IN 'DECAYDAY' time ), or how many days it takes for a tracer to lose so much % of it smass
DECAYDAY = 0.125D0
= 1/8 of a day, just to make the simulation quick
DECAYRATE = 0.9D0
=> T90, a T50 would be with DECAYRATE = 0.5D0

And we assume NUMBER OF TRACERS = 3 ...

You can write for instance for TRACER #2:
- For TRACER #2, here I used a balanced formulation with both S0TA/2 + TA*S1TA/2 ...

! UNIFORM DECAY RATE FOR TRACER NUMBER 2 (%ADR(2)) -
IF(NTRAC.GE.2) THEN
!
S0TA%ADR(2)%P%TYPR='Q'
S1TA%ADR(2)%P%TYPR='Q'
DO IPOIN3 = 1,NPOIN3
S0TA%ADR(2)%P%R(IPOIN3) = TA%ADR(2)%P%R(IPOIN3) *
& DLOG( 1.D0-DECAYRATE ) / ( DECAYDAY*24*3600.D0 ) / 2.D0
S1TA%ADR(2)%P%R(IPOIN3) =
& - DLOG( 1.D0-DECAYRATE ) / ( DECAYDAY*24*3600.D0 ) / 2.D0
ENDDO
!
ENDIF

You can write for instance for TRACER #3:
- For TRACER #3, here I used an implicit formulation with only TA*S1TA/2 ...
but ( DECAYDAY / IPLAN ) for each plan to see vertical variation


! VERTICALLY VARIABLE DECAY RATE FOR TRACER NUMBER 3 (%ADR(3))
! ... A FUNCTION OF PLAN NUMBER TO DO IT SIMPLY ... (IPLAN=1 IS THE BOTTOM)
IF(NTRAC.GE.3) THEN
!
S1TA%ADR(3)%P%TYPR='Q'
DO IPOIN2 = 1,NPOIN2
DO IPLAN = 1,NPLAN
IPOIN3 = IPOIN2 + (IPLAN-1)*NPOIN2
S1TA%ADR(3)%P%R(IPOIN3) =
& - DLOG( 1.D0-DECAYRATE ) / ( DECAYDAY*24*3600.D0 ) * IPLAN
ENDDO
ENDDO
!
ENDIF

Finally, for TEMPERATURE dependence, if you name your TARCERS, with one at least being TEMPERATURE for instance (or SALINITY)
NAMES OF TRACERS ='TEMPERATURE';'UNIFORM DECAY';'DEPTH DECAY'

Then in your PRINCI, you can directly access the value of TRACER #1 at any IPOIN with TA%ADR(1)%P%R(IPOIN) and use it the DECAY formulation ...
or in a mor generic way, you can used IND_T (or IND_S for SALINITY), which is the indice for which TEMPERATURE (SALINITY) has been recognised automatically by TELEMAC: TA%ADR(IND_T)%P%R(IPOIN)

Hope this helps.

Sébastien.

Thank you very much Sébastien for your invaluable help with this issue.
I have incorporated the effect of temperature and salinity on the decay of the tracer by slightly modifying the code you sent me.
But do you still need a negative sign in
S1TA%ADR(3)%P%R(IPOIN3) =
& - DLOG( 1.D0-DECAYRATE ) / ( DECAYDAY*24*3600.D0 ) * IPLAN

I thought you don't since:
TA *( 1.-S1TA ) = TA + ... + ADVECTION+DIFFUSION

The modified lines are as follows:

DOUBLE PRECISION k1,k2,k3

IF(NTRAC.GE.3) THEN
!
S1TA%ADR(3)%P%TYPR='Q'
DO IPOIN2 = 1,NPOIN2
DO IPLAN = 1,NPLAN
IPOIN3 = IPOIN2 + (IPLAN-1)*NPOIN2
S1TA%ADR(3)%P%R(IPOIN3) =
& (k1+k2*TA%ADR(2)%P%R(IPOIN3))*
& k3**(TA%ADR(1)%P%R(IPOIN3)-20.D0)

I tested it on the simple example (square pond) and it produced exactly the analytical solution so that's sorted thank you so much.
I hope now to test this in the main study I am working on. Thanks alot for your help.
Zeinab

Dear All,
I tried to run the decay rate subroutine source_trac on the original model of Dublin Bay which is in Version 5P5 but requires substantial work to update it to a recent version. I get the following error message regarding TYPR in S1TA%ADR(3)%P%TYPR and S0TA%ADR(3)%P%TYPR saying : This is not a field name that is defined in the encompassing structure even though the version of source_trac is V5P5.
Any ideas what do I need use in the declarations to overcome that?
Thanks a lot,
Zeinab

Dear Zeinab,

This is simple, the component TYPR was added in version 5.6, so you will not find it in version 5.5. In this latter version, the field S1TA%ADR(3)%P%R is initialised to 0, and you must treat it in any case because you do not have the hint of TYPR saying that all is zero (TYPR just allows an optimisation). Anyway in your case S1TA%ADR(3)%P%R is not zero, so, compared to version 6.0 you just have to remove the test on TYPR.

With best regards, have a good week-end,

Jean-Michel Hervouet

Thanks alot Jean-Michel,
Apologies for the late reply. Hadn't realised that you responded to my query.
I assumed that TYPR has appeared in later versions as V5.8 seems to recongnise it with no problems.
I deleted the TYPR lines from the subroutine, assuming that it will be set to zero. But gave a value for S1TA%ADR(3)%P%TYPR and the model ran well.
Thank you so much for the response.
Your support and help is very much appreciated.
Best Regards

Zeinab

Zeinab,

I am happy to see it works.
The negative sign was because DLOG(1-d) is negative. You are right to say that you do not need it in your case.

Best regards,
Sébastien.