README.md

adists_sediment_cores

Python script to extract numerical sediment cores from adists binary results

données:

• h[0:tmax] : liste des hauteurs totales de sédiment à chaque pas de temps
• masse[0:nb_sediment*nb_affluents][0:tmax] : la masse de chaqun des sédiments qui constituent le fond a chaque pas de temps adis-ts n'a qu'une seule couche et ne mélange pas les sédiments : on a besoin d'un algo pour recréer ces couches note : on part du principe qu'on n'a rien pour relier la masse à la hauteur

sorties :

• la hauteur de chaque couche
• sa proportion en chaque sédiment

algorythm : (pseudo-code, ne fonctionne pas)

layers = [] # liste de couches, vide pour l'instant
l = 0 # numéro de la couche courante
for t in range(nb_timesteps):
    # on crée une nouvelle couche en fin de liste
    # une couche contient:
    # [0] : la hauteur de sédiments
    # [1] : la masse des sédiments dans la couche
    # [2:nb_sediments+2] : la proportion de chaque classe de sédiments dans la couche
    layers.append(np.zeros(nb_sediments+2,float))
    layers[l][0] = h[t]-h[t-1] # hauteur potentiellement négative
    layers[l][1] = 0.0
    for s in range(nb_sediments):
        # pour l'instant on met la masse de chaque sédiment dans layers[l][2:nb_sediments+2]
        # la masse d'une classe de sédiment est la somme des masses des classes correspondantes de chaque affluent
        # il s'agit d'une diférence avec le pas de temps précédent, ces masses peuvent donc être négatives si il y a érosion
        layers[l][s+2] = sum(masse[s+j][t]-masse[s+j][t-1] for j in range(0,nb_afluent*nb_sediments,nb_sediments)
        residual = layers[l][s+2]
        if residual < 0.0:
            # il y a eu érosion pour cette classe
            # dans notre carotte on ne peut pas avoir une masse négative
            layers[l][s+2] = 0.0
        # on ajoute la masse de cette classe à la masse totale de la couche
        layers[l][1] += layers[l][s+2]
        l2=l
        # si residual est négatif alors on doit aller chercher du sédiment dans les couches du dessous
        while residual < 0.0:
            l2 -= 1
            # calcul préliminaire
            diff = min(-1*residual, layers[l2][s+2]) # diff > 0 : ce que l'on doit prélever
            ratio = (layers[l2][1]-diff)/layers[l2][1] # ratio < 1 : pour réduire l'épaisseur de la couche.
            dh = layers[l2][0]-(ratio*layers[l2][0]) # dh > 0
            # modif de la couche
            residual += layers[l2][s+2]
            layers[l2][s+2] -= diff
            layers[l2][1] -= diff
            layers[l2][0] -= dh
            layers[l][0] += dh # doit finir par compenser les hauteurs négatives
        # fin du while sur residual
    # fin de la boucle sur les classes de sédiments
    l += 1
# on a fini, on va refaire une boucle sur les couches pour obtenir les proportions
for layer in layers:
    for s in range(nb_sediments):
        if layers[l][1] > 0.0: # normalement jamais négatif
            # on transforme les masses en proportions
            layer[s+2] /= layer[1]