Menu ▾ ▴

Agregar_Modelos

Anonymous

Nuevos modelos en NChem

Es posible agregar nuevos modelos a NChem, creando el código fuente necesario en las siguientes unidades de programa:

  • Subrutina con los datos de configuración del modelo (sub1);
  • Subrutina de cálculo de los elementos de la matriz de coeficientes estequiométricos del modelo (sub2);
  • Subrutina de cálculo de los procesos cinéticos del modelo (sub3);
  • Módulo de declaración de parámetros del modelo (mod1) (OPCIONAL);
  • Subrutina de inicialización de "alias" de parámetros y variables de estado del modelo (sub4) (OPCIONAL);
  • Factores de conversión para las variables de estado (fnc1) (OPCIONAL);
  • Subrutina de cálculo del pH en función de la carga eléctrica neta (sub5) (OPCIONAL);
  • Función de cálculo de la carga eléctrica neta en función del pH (fcn2) (OPCIONAL)
  • Función que proporciona la leyenda del fichero Stoichiometry.txt (fcn3).

y modificando el código de la subrutina ASPMOD de asociación de punteros a las subrutinas y funciones propias del modelo. La subrutina ASPMOD se encuentra en el fichero pointer.f

Modificación de la subrutina de asociación de punteros ASPMOD

ASPMOD(MODELO,IERR)

Esta subrutina pertenece al módulo PUNTEROS. El objeto de esta subrutina es asociar diversos punteros a las subrutinas y funciones propias del modelo identificado por la etiqueta alfanumérica MODELO.

Argumentos

CHARACTER*8 MODELO ! Etiqueta que identifica el modelo a utilizar en los cálculos (INPUT)
INTEGER IERR !Identificador de error devuelto por la subrutina (OUTPUT)

Declaración de variables

  • Las seis subrutinas propias del modelo (sub1, sub2, sub3, sub4, sub5 y sub6) deben declararse como EXTERNAL.

  • La funciones fnc1 y fcn2 deben declararse como REAL(DOUBLE) y EXTERNAL, siendo DOUBLE=KIND(1.d0).

    REAL(DOUBLE) fnc1, fcn2
    EXTERNAL sub1, sub2, sub3, sub4, sub5, sub6, fnc1, fcn2

Código de ejemplo

En la estuctura IF...ELSE IF...ELSE...END IF, debe crearse una opción para el valor de la etiqueta alfanumérica del nuevo modelo, dentro del cual se hará la asociación de los punteros: 

ELSE IF (INDEX(MODELO,constante).EQ.1 ) THEN
  DATMODEL => sub1 ! sub1 = nombre de la subrutina con los datos de configuración del modelo 
  LEYENDA  => fnc3 ! fcn3 = nombre de la función con la leyenda del fichero 
                   !        Stoichiometry.txt 
  SUBESTEQ => sub2 ! sub2 = nombre de la subrutina de cálculo de los elementos de la matriz de coeficientes
                   !        estequiométricos del modelo
  PROCESOS => sub3 ! sub3 = nombre de la subrutina de cálculo de los procesos cinéticos del modelo
  INITMOD  => sub4 ! sub4 = nombre de la subrutina de inicialización de "alias" de parámetros y variables de estado del
                   !        modelo. Si el modelo no necesita la función sub5, INITMOD debe apuntar a NULL()
  SELCONT  => fnc1 ! fnc1 = nombre de la función REAL(DOUBLE) de cálculo del balance estequiométrico de materia y carga
                   !        eléctrica. Si el modelo no necesita la función fnc1, SELCONT debe apuntar a NULL()
  CALCPH   => sub5 ! sub5 = nombre de la subrutina de cálculo del pH en función de la carga neta. Si el modelo no calcula
                   !        el pH, CALCPH debe apuntar a NULL().
  BALZ     => fnc2 ! fnc2 = nombre de la función REAL(DOUBLE) de cálculo de la carga eléctrica neta en función del pH; 
                   !        Si el modelo no calcula el pH, CALCPH debe apuntar a NULL()
ELSE

donde "constante" es una constante alfanumérica que se comparará con el valor del argumento MODELO.

Subrutina con los datos de configuración del modelo (sub1)

SUBROUTINE sub1()

En esta subrutina se deberá dar valores a algunas de las variables declaradas en el módulo GLOBAL:

  • NCONT, número de ecuaciones de continuidad
  • NVARS, número de variables de estado del modelo
  • NPARS, número de parámetros del modelo
  • NPROC, número de procesos del modelo
  • MODELO, etiqueta identificativa del modelo
  • VAR_PH, número de variable de estado que corresponde al pH
  • CNAME(NCONT), nombre de los componentes de las ecuaciones continuidad
  • ROTVAR(NVARS), nombre de las variables de estado
  • YREF(NVARS), valores de referencia de variables de estado
  • FASE(NVARS), indicador de la fase de cada variable de estado
  • ROTPAR(NPARS), nombre de los parámetros del modelo
  • PREF(NPARS), valores de referencia de los parámetros
  • ROTPRC(NPROC), nombre de los procesos

Igualmente debe reservarse memoria para estas variables mediante la llamada a la subrutina MEMMODEL.

La llamada a la subrutina MEMMODEL debe efectuarse después de asignados los valores a las variables NCONT, NVARS, NPARS y NPROC pues son éstas las que determinan el espacio de memoria necesario para almacenar el resto de variables.

Declaración de variables

      USE GLOBAL

Código de ejemplo

Un ejemplo de código de esta subrutina sería: 

C 0 Dimensiones del modelo
C 0.1 Número de ecuaciones de continuidad (INTEGER)
      NCONT = 2 
C 0.2 Número de variables de estado (INTEGER)
      NVARS = 10 ! INTEGER
C 0.3 Número de parámetros (INTEGER)
      NPARS = 27
C 0.4 Número de procesos (INTEGER)
      NPROC = 18
C 0.5 Reserva de memoria para las variables del modelo
      CALL MEMMODEL
C 1 Identificación del modelo
C 1.1 Nombre del modelo (CHARACTER*8)
      MODELO = 'GaP'
C 1.2 Posición del pH entre las variables de estado (INTEGER) 
      VAR_PH = 0 !0 indica que el modelo no utiliza el valor del pH
C 1.3 Nombre de los componentes de las ecuaciones de continuidad (CHARACTER*10)
      CNAME(1) = 'Nitrogen'
      CNAME(2) = 'COD'
C 2 Variables de estado
C 2.1 Nombre de las variables de estado (CHARACTER*10, DIMENSION(NVARS))
      ROTVAR(1)  = 'So2'
      ... 
      ROTVAR(10) = 'Xi'
C 2.2 Valores de referencia (REAL(DOUBLE), DIMENSION(NVARS))
      YREF(1)  =   0.0  ! So2
      ...
      YREF(10) =  25.0  ! Xi
C 2.3 Fases de las variables de estado (CHARACTER, DIMENSION(NVARS))
      FASE(1)='s'      
      ...
      FASE(10)='x' 
C 3 Parámetros
C 3.1 Nombre de los parámetros (CHARACTER*10, DIMENSION(NPARS))
      ROTPAR(1)  = 'ynh4'
      ...
      ROTPAR(27) = 'bpred'
C 3.2 Valores de referencia (REAL(DOUBLE), DIMENSION(NPARS))
      PREF(1)  =  0.18    ! ynh4
      ...
      PREF(27) =  3.60    ! bpred
C 4 Procesos
C 4.1 Nombre de los procesos (CHARACTER*10, DIMENSION(NPROC))
      ROTPRC(1)  = 'GrXnh4'
      ...
      ROTPRC(18) = 'Aeration'

Función con leyenda del fichero Stoichiometry.txt (fcn3)

CHARACTER(80) FUNCTION FCN3(TIPO,I)

Esta función ha de proporcionar el valor de la leyenda para cada elemento I de los componentes (TIPO='C') o de los procesos del modelo. La escritura se realiza en el fichero Stoichiometry.txt

Argumentos

INTEGER I
CHARACTER TIPO

I indica el número de componente o de proceso para el que se requiere la leyenda.

TIPO indica si se debe producir la escritura de la información relativa a los componentes (TIPO='C') o a los procesos.

Código de ejemplo

IF (TIPO=='C') THEN
  SELECT CASE(I) ! Para los NVAR componentes
  CASE(1)
    FCN3='Dissolved oxygen (mg O2/l)'
  CASE(2)
    FCN3='Ammonium plus ammonia nitrogen (mg N/l)'
    ...
  END SELECT
ELSE ! Para los NPROC procesos
  SELECT CASE(I)
  CASE(1)
    FCN3='Growth of ammonia oxidiziers'
  CASE(2)
    FCN3='Maintenance of ammonia oxidiziers'
    ...
  END SELECT 
END IF

Subrutina de cálculo de los elementos de la matriz de coeficientes estequiométricos del modelo (sub2)

SUBROUTINE SUB2(N,NU)

Esta subrutina debe calcular los valores de cada uno de los elementos de la matriz de coeficientes estequimétricos del modelo NU(1:NPROC,1:NVARS).

Argumentos

INTEGER N
REAL(DOUBLE) NU(N,*)

Declaración de variables

Los valores de las constantes cinéticas deberán obtenerse del vector PARAM declarado en el módulo GLOBAL o bien a través de las variables declaradas en el módulo mod1: 

USE GLOBAL, ONLY: NVARS,PARAM

Código de ejemplo

C Inicializacion de los elementos de la matriz NU()
      NU(:,:NVARS) = 0.d0
C NU(proceso,variable de estado):
C Growth on Xnh4
      NU(1,so2)=ONE - 3.43d0/ynh4
      NU(1,snh4)=-ONE/ynh4 - inbm
      NU(1,sno2)=ONE/ynh4
      NU(1,xnh4)=ONE
C Maintenance of Xnh4
      NU(2,so2)=-3.43d0
...

Subrutina de cálculo de los procesos cinéticos del modelo (sub3)

SUBROUTINE SUB3(TEMP,KLA,S,PH)

Calcula los valores de los procesos cinéticos en un reactor a partir de los valores de los prámetros del modelo y las variables de estado en ese reactor.

Argumentos

REAL(DOUBLE) TEMP,KLA,S(*)
REAL(DOUBLE), OPTIONAL :: PH

Declaración de variables

A los valores de las constantes biocinéticas del modelo y de las variables de estado en el reactor que está siendo calculado se accede, respectivamente, mediante las variables PARAM(NPARS) y VAR(NVARS) declaradas en el módulo GLOBAL, así que será necesario incluir la instrucción: 

USE GLOBAL, ONLY: PARAMS, VAR

El programa de NChem hace, que previamente a la llamada a esta subrutina, el vector VAR contenga los valores de las variables de estado del reactor que está siendo calculado.

Código de ejemplo

Se deben calcular los valores de cada elemento del vector S(NPROCS)

A fin de mejorar la legibilidad del código de la subrutina, pueden utilizarse las siguientes opciones:

  • Alternativamente al empleo del vector PARAM, se puede utilizar una lista de variables cuyos nombres identifiquen a cada parámetro. Estas variables estarán de declaradas en el módulo mod1 y por tanto deberá declararse su uso al comenzo de la subrutina. P ej.: 
     USE PARAMS_GAP, ONLY: munh4, ko2nh4, knh4nh4, mnh4, bnh4, enh4,
&  muno2, ko2no2, kno2no2, mno2, bno2, eno2, muh, ko2h, kssh,
&  kno3h, mh, eh, bh, mupre, ko2pred, bpred

,siendo PARAMS_GAP el nombre particular del módulo mod1 para este modelo.

Previamente, deberá asignárse los valores de los parámetros en la subrutina sub5.

  • Alternativamente al empleo del vector VAR(NVARS), se puede utilizar una lista de variables POINTER cuyo nombre identifique a la variable de estado y alineadas en un bloque COMMON rotulado: 
    REAL(DOUBLE), POINTER :: So2, Snh4, Sno2, Sno3, Ss, Xnh4, Xno2, Xh, Xpre, Xi
COMMON/VARS_GAP/ So2, Snh4, Sno2, Sno3, Ss, Xnh4, Xno2, Xh, Xpre, Xi

,siendo VARS_GAP el nombre particular del rótulo de este bloque COMMON para este modelo.

Previamente, la subrutina sub5 deberá asociar cada puntero a cada elemento del vector VAR

El resto de la subrutina deberá establecer el cálculo de los diferentes elementos del vector S, p. ej.: 

C Growth on Xnh4
      s(1)  = munh4 * (so2/(ko2nh4+so2)) *
     &                (snh4/(knh4nh4+snh4)) * xnh4
C Maintenance of Xnh4
      s(2)  = mnh4 * (so2/(ko2nh4+so2))*
     &               (snh4/(knh4nh4+snh4)) * xnh4
...

Subrutina de inicialización de "alias" de parámetros y variables de estado del modelo (sub4)

SUBROUTINE SUB4()

Esta subrutina deberá establecer:

  • una lista de variables REAL(DOUBLE), POINTER como alias de las variables de estado a cuyos elementos se les asocian los elementos del vector VAR(NVARS) declarados en el módulo GLOBAL. Estas variables deberán estara lineadas en un bloque COMMON rotulado;
  • una lista de variables REAL(DOUBLE) como alias de los parámetros biocinéticos del modelo y asigna valores almacenados en el vector PARAM(NPARS) declarado en el módulo GLOBAL

En ambos casos: DOUBLE=KIND(1.d0)

Por ejemplo: 

      USE GLOBAL, ONLY: PARAM, VAR
      USE PARAMS_GAP, ONLY: ynh4, yno2, yh, ypred, fxi, munh4, ko2nh4,
     &  knh4nh4, mnh4, bnh4, enh4, muno2,  ko2no2, kno2no2, mno2, bno2,
     &  eno2, muh, ko2h, kssh, kno3h, mh, eh, bh, mupre, ko2pred, bpred

      REAL(DOUBLE), POINTER :: So2, Snh4, Sno2, Sno3, Ss, Xnh4,
     &  Xno2, Xh, Xpre, Xi
      COMMON/VARS_GAP/ So2, Snh4, Sno2, Sno3, Ss, Xnh4, Xno2, Xh,
     &  Xpre, Xi

      So2  => VAR(1)
      Snh4 => VAR(2)
...      
      ynh4  =   PARAM(1)
      yno2  =   PARAM(2)
...

donde siendo PARAMS_GAP y VARS_GAP, respectivemente, los nombres particulares del módulo mod1 y del rótulo del bloque COMMON de variables de estado para este modelo.

Módulo de declaración de parámetros del modelo (mod1)

El único proposito de este módulo es la declaración de las variables correspondientes a cada uno de los parámetros biocinéticos y constantes estequiométricas del modelo. Por ejemplo: 

      MODULE PARAMS_GAP
        USE CONSTANTES, ONLY: DOUBLE
        REAL(DOUBLE) :: ynh4, yno2 , yh, ypred, fxi, munh4, ko2nh4,
     &    knh4nh4, mnh4, bnh4, enh4, muno2, ko2no2, kno2no2, mno2,
     &    bno2, eno2, muh, ko2h, kssh, kno3h, mh, eh, bh, mupre,
     &    ko2pred, bpred
      END MODULE PARAMS_GAP

Factores de conversión para las variables de estado (fnc1)

REAL(DOUBLE) FUNCTION FCN1(IVAR,IREF)

Esta función REAL(DOUBLE) debe proporcionar el factor de conversión (relación estequiométrica) del componente de referencia IREF en la variable de estado IVAR.

Argumentos

INTEGER IVAR, IREF

Código de ejemplo

Para un modelo de diez variables de estado donde se emplean las equivalencias de éstas en nitrógeno y DQO, un posible código sería: 

 INTEGER I
 REAL(DOUBLE) REF(10,2)
C                           N         DQO
 DATA (REF(1,I),I=1,2)  /  0.00D+0, -1.00D0/    ! So2
 DATA (REF(2,I),I=1,2)  /  1.00D+0,  0.00D0/    ! Snh4
...
 DATA (REF(10,I),I=1,2) /  2.00D-2,  1.00D0/    ! Xi

 FCN1 = REF(IVAR,IREF)

Subrutina de cálculo del pH en función de la carga eléctrica neta (sub5)

SUBROUTINE SUB5(TEMP,Z,PH)

Argumentos

REAL(DOUBLE) TEMP, Z, PH

, siendo TEMP la temperatura en ºC, Z la carga eléctrica neta (cationes inertes - aniones inertes) en mol/l y PH el valor del pH calculado.

Por "inertes" entendemos las variables de estado constituidas por especies iónicas que no son alteradas por los procesos que el modelo describe.

Función de cálculo de la carga eléctrica neta en función del pH (fcn2)

REAL(DOUBLE) FUNCTION FNC2,TEMP,IREAC,NVARS,Y)

Calcula el exceso de cationes inertes sobre aniones inertes en un reactor IREAC, aplicando un balance de cargas eléctricas en el mismo:

CAT_ADM = CATIONESI - ANIONESI = ANIONESNI - CATIONESNI

, siendo:

  • CATIONESI = concentraciones de cationes inertes (no alterados por los procesos del modelo), en mol/l;
  • CATIONESNI = concentraciones de cationes no inertes (alterados por los procesos del modelo), en mol/l;
  • ANIONESI = concentraciones de aniones inertes (no alterados por los procesos del modelo), en mol/l;
  • ANIONESNI = concentraciones de aniones no inertes (alterados por los procesos del modelo), en mol/l

Argumentos

INTEGER IREAC
INTEGER NVARS
REAL(DOUBLE) Y(NVARS,*)
REAL(DOUBLE) PH
REAL(DOUBLE) TEMP

, siendo IREAC el número de reactor; NVARS el número de variables de estado; Y(NVARS,NREAC) la matriz variables de estado en los rectores y PH y TEMP, respectivamente, el valor del pH y la temperatura en ºC en el reactor IREAC.

auth0 logo