! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Auxiliary Routines File
! 
! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
! 
! File                 : CB05_CSIRO_Util.f90
! Time                 : Thu Oct  8 08:21:51 2020
! Working directory    : /Users/lechriso/git/boxmox/boxmox/tmp_CB05_CSIRO_26998_BOXMOX
! Equation file        : CB05_CSIRO.kpp
! Output root filename : CB05_CSIRO
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



MODULE CB05_CSIRO_Util

  USE CB05_CSIRO_Parameters
  USE CB05_CSIRO_Global
  USE CB05_CSIRO_Monitor
  IMPLICIT NONE

CONTAINS



! User INLINED Utility Functions


! ****************************************************************************
! ****************************************************************************
! BOXMOX utilities
! ****************************************************************************
! ****************************************************************************


! ****************************************************************************
! toUPPER - convert ASCII characters to all upper case
!           from http://stackoverflow.com/questions/10759375/
!           how-can-i-write-a-to-upper-or-to-lower-function-in-f90
! ****************************************************************************
      FUNCTION toUPPER (str) RESULT (string)

      CHARACTER(*), INTENT(In) :: str
      CHARACTER(LEN(str))      :: string

      INTEGER                  :: ic, i

      CHARACTER(26), PARAMETER :: cap = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
      CHARACTER(26), PARAMETER :: low = 'abcdefghijklmnopqrstuvwxyz'

      string = str
      DO i = 1, LEN_TRIM(str)
        ic = INDEX(low, str(i:i))
        IF (ic > 0) string(i:i) = cap(ic:ic)
      END DO

      END FUNCTION toUPPER
! ****************************************************************************


! ****************************************************************************
! hash - a string hashing function
!        from flibs, http://flibs.sourceforge.net
! ****************************************************************************
      INTEGER FUNCTION hash( key )

      CHARACTER(len=*), INTENT(in) :: key

      INTEGER                      :: i

      INTEGER, PARAMETER           :: hash_size  = 49931
      INTEGER, PARAMETER           :: multiplier = 312

      hash = 0

      DO i = 1,LEN(key)
          hash = MODULO( multiplier * hash + ICHAR(key(i:i)), hash_size)
      ENDDO

      hash = 1 + MODULO( hash-1, hash_size )
      END FUNCTION hash
! ****************************************************************************


! ****************************************************************************
! strEqual - comparison function matching two strings
! ****************************************************************************
      LOGICAL FUNCTION strEqual(a, b)

      CHARACTER(LEN=*), INTENT(IN)   :: a, b

      LOGICAL, PARAMETER             :: ignoreCase = .TRUE.

      CHARACTER(LEN=100)             :: nice_a, nice_b
      INTEGER                        :: max_len

      nice_a = ADJUSTL(a)
      nice_b = ADJUSTL(b)

      IF (ignoreCase) THEN
        nice_a = toUPPER(nice_a)
        nice_b = toUPPER(nice_b)
      END IF

      max_len = MAX( LEN_TRIM(nice_a), LEN_TRIM(nice_b) )

      strEqual = nice_a(1:max_len) .EQ. nice_b(1:max_len)

      END FUNCTION strEqual
! ****************************************************************************


! ****************************************************************************
! writespecstring - Utility function that plots a (potentially)
!                   large number of species
! ****************************************************************************
      SUBROUTINE writespecstring(nspecies, txt, species, offset)

      INTEGER             :: nspecies, offset
      CHARACTER(LEN=*)    :: txt
      CHARACTER(LEN=15)   :: species(nspecies), fmt
      CHARACTER(LEN=2000) :: specstring

      INTEGER, PARAMETER  :: cperline=60

      INTEGER             :: i, l

      fmt = ""
      WRITE(fmt, "(A1,I3,A4)") "(", offset, "X,A)"

      i = 1; l = 1
      specstring = ""

      DO WHILE (i <= nspecies)

        specstring = TRIM(ADJUSTL(specstring))//TRIM(ADJUSTL(species(i)))
        IF (i < nspecies) specstring = TRIM(ADJUSTL(specstring))//", "

        i = i + 1

        IF (LEN_TRIM(specstring) > cperline) THEN
          IF (l == 1) specstring = txt//specstring
          l = l + 1
          WRITE(*,fmt) TRIM(specstring)
          specstring = ""
        ENDIF
      ENDDO

      IF (LEN_TRIM(specstring) > 0) THEN
        IF (l == 1) specstring = txt//specstring
        WRITE(*,fmt) TRIM(specstring)
      ENDIF

      END SUBROUTINE writespecstring
! ****************************************************************************


! ****************************************************************************
! calcCFACTOR - Calculates CFACTOR variable for unit conversion between
!               input/output files and values used in calculations / rates
! ****************************************************************************
      SUBROUTINE calcCFACTOR ( )

      ! from ppmv to mlc/cm3:
      ! 1e-6 ppmv -> vmr
      ! 1e-3 kg/m3 -> g/cm3
      ! 1/28.9645 g/cm3 -> mole/cm3
      ! 6.022e23 -> molec/cm3
      ! then multiply with air density rho
      ! 1.0E-6 * 1E-03 * (1.0/28.97) * 6.022e23 * rho
      ! CFACTOR = 2.598378e+13; (with rho == 1.25)

      CFACTOR = 1.0E-6 * 1E-03 * (1.0/28.9645) * 6.022e23 * &
                pres / (287.058 * temp)

      END SUBROUTINE calcCFACTOR
! ****************************************************************************


! ****************************************************************************
! readConfig - Sets variables using values read in from "BOXMOX.nml"
! ****************************************************************************
      SUBROUTINE readConfig ( )

      INTEGER, PARAMETER          :: config_pid   = 22
      CHARACTER(LEN=*), PARAMETER :: config_name  = "BOXMOX.nml"

      NAMELIST /BOXMOX/ TSTART, TEND, DT, lverbose, iturb, lfix_nox

      LOGICAL            :: lnml_exists

      ! *** default values section ***

      TSTART   = 0.0D0              ! model start time, seconds since 0 LT
      TEND     = TSTART + 3600.0D0  ! model end time, seconds since 0 LT
      DT       = 10.D0              ! model time step, seconds

      lverbose = .FALSE.
      iturb    = TURB_OPT_NONE

      lfix_nox = .FALSE.

      ! *** namelist read-in section ***

      INQUIRE(FILE=config_name, EXIST=lnml_exists)

      IF (lnml_exists) THEN

        WRITE(*,*) "BOXMOX configuration file ", TRIM(config_name)," found."

        OPEN (config_pid, FILE=config_name)
        READ (config_pid, nml=BOXMOX)
        CLOSE(config_pid)

      ELSE
        WRITE(*,*) "No BOXMOX configuration file (", &
                         TRIM(config_name),") exists."
      ENDIF

      ! *** configuration printout ***

      WRITE(*,*) " "
      WRITE(*,"(A20,2X,F12.3)") "TSTART", TSTART
      WRITE(*,"(A20,2X,F12.3)") "TEND", TEND
      WRITE(*,"(A20,2X,F12.3)") "DT", DT
      WRITE(*,*) ""
      WRITE(*,"(A20,2X,L10)")   "lverbose", lverbose
      WRITE(*,*) ""
      WRITE(*,"(A20,11X,I1)")   "iturb", iturb
      WRITE(*,*) ""
      WRITE(*,"(A20,2X,L10)")   "lfix_nox", lfix_nox

      WRITE(*,*) " "

      END SUBROUTINE readConfig
! ****************************************************************************


! ****************************************************************************
! checkConfig - check configuration variables
! ****************************************************************************
      SUBROUTINE checkConfig ( )

      INTEGER                      :: status

      status = 0

      IF (TSTART < 0.0) THEN
        WRITE(*,*) "start time TSTART cannot be < 0.0"
        status = 1
      ENDIF
      IF (TEND < 0.0) THEN
        WRITE(*,*) "end time TEND cannot be < 0.0"
        status = 2
      ENDIF
      IF (TSTART >= TEND) THEN
        WRITE(*,*) "start time TSTART has to be < end time TEND"
        status = 3
      ENDIF
      IF (DT <= 0.0) THEN
        WRITE(*,*) "time step DT has to be > 0.0"
        status = 4
      ENDIF

      IF (status .NE. 0) THEN
        IF (.NOT. llibrary) THEN
          STOP
        ENDIF
      ENDIF

      END SUBROUTINE checkConfig
! ****************************************************************************


! ****************************************************************************
! writeVersion - Write "VERSION" file to run directory
! ****************************************************************************
      SUBROUTINE writeVersion ( )

      INTEGER, PARAMETER          :: version_pid  = 23
      CHARACTER(LEN=*), PARAMETER :: version_name = "VERSION"

      ! *** write BOXMOX version to output directory ***

      OPEN(version_pid, file=version_name)
      WRITE(version_pid, "(A)") "__BOXMOX_VERSION__"
      CLOSE(version_pid)

      END SUBROUTINE writeVersion
! ****************************************************************************


! ****************************************************************************
! readVals   - main input routine for the BOXMOX extensions.
!              reads values from input files with different time
!              formats, and interpolates to get the values
!              appropriate for the current time step.
! ****************************************************************************
      SUBROUTINE readVals(fname, nvar, species, values, stat, &
                          nanc, ancspecies, ancvalues, lquiet)

      CHARACTER(LEN=*), INTENT(IN)             :: fname
      INTEGER, INTENT(OUT)                     :: nvar
      CHARACTER(LEN=15), INTENT(OUT)           :: species(max_species_boxmox)
      REAL(KIND=dp), INTENT(OUT)               :: values(max_species_boxmox)
      INTEGER, INTENT(OUT)                     :: stat
      INTEGER, INTENT(OUT), OPTIONAL           :: nanc
      CHARACTER(LEN=15), INTENT(OUT),OPTIONAL :: ancspecies(max_species_boxmox)
      REAL(KIND=dp), INTENT(OUT), OPTIONAL     :: ancvalues(max_species_boxmox)
      LOGICAL, INTENT(IN), OPTIONAL            :: lquiet

      INTEGER, PARAMETER :: input_pid = 21
      REAL(KIND=dp)      :: time_cur, weight_t1, &
                            values_t0(max_species_boxmox), &
                            values_t1(max_species_boxmox), &
                            values_cur(max_species_boxmox)
      CHARACTER(LEN=15)  :: names(max_species_boxmox)
      INTEGER            :: iostat, time_format, ntoread, ntime, nanct, pid, i
      LOGICAL            :: lfile_exists, ltalkative

      names(:)     =  ""
      values_t0(:) =  0.0
      values_t1(:) =  0.0

      stat         = 0 ! everything is allright

      ltalkative = ( TIME .EQ. TSTART ) .OR. lverbose
      IF (PRESENT(lquiet)) THEN
        IF (lquiet) THEN
          ltalkative = .FALSE.
        ENDIF
      ENDIF

      ! check whether file exists
      INQUIRE(FILE=fname, EXIST=lfile_exists)
      IF (.NOT. lfile_exists) THEN
        stat = 1
        RETURN
      ELSE
        IF (ltalkative) THEN
          WRITE(*,*) ""
          WRITE(*,*) " * ", fname
        ENDIF
      END IF

      OPEN(input_pid, file=fname)

      ! line 1: number of variables
      READ(input_pid, *) nvar

      ! if no variable should be read - bail
      IF (nvar == 0) THEN
        IF (ltalkative) WRITE(*,*) fname, " declared 0 variables to be read"
        stat = 2
        CLOSE(input_pid)
        RETURN
      END IF

      ! line 2: number of ancillary variables
      READ(input_pid, *) nanct

      IF (nanct .LT. 0) THEN
        WRITE(*,"(A)") "Number of ancillary variables must be >= 0"
        CLOSE(input_pid)
        STOP
      ENDIF

      IF (nanct .GT. 1) THEN
        WRITE(*,"(A)") "Number of ancillary variables must be < 2"
        CLOSE(input_pid)
        STOP
      ENDIF

      ! line 3: time format
      READ(input_pid, *) time_format

      IF (time_format .LT. 0 .OR. time_format .GT. 2) THEN
        WRITE(*,"(3A)") "Time format has to be 0 (constant), ", &
                                              "1 (s since TSTART) or ", &
                                              "2 (hour of day)"
        CLOSE(input_pid)
        STOP
      END IF

      ntime = 0
      IF (time_format .GT. 0) ntime = 1

      ntoread = ntime + nanct + nvar

      ! if more variables than max_species_boxmox should be read - bail
      IF (ntoread > max_species_boxmox) THEN
        IF (ltalkative) THEN
          WRITE(*,*) fname, " declared more variables to be read"
          WRITE(*,*) "than allocated. Increase max_species_boxmox,"
          WRITE(*,*) "recompile, and run again."
        ENDIF
        stat = 3
        CLOSE(input_pid)
        RETURN
      END IF

      ! line 4: header line
      READ(input_pid, *) names(1:ntoread)

      ! line 4ff: values
      READ(input_pid, *) values_t0(1:ntoread)

      ! if its not constant in time values
      weight_t1 = 0
      IF (.NOT. time_format .EQ. 0) THEN
        READ(input_pid, *, IOSTAT=iostat) values_t1(1:ntoread)

        ! TIME is the KPP provided current model time
        time_cur = TIME
        IF (time_format .EQ. 2) THEN
          time_cur = TIME/3600.0_dp - (INT(TIME/3600.0_dp)/24)*24
        END IF

        ! ASSUMPTION: time is the first value read
        DO WHILE(values_t1(1) .LT. time_cur .AND. iostat .EQ. 0)
          values_t0 = values_t1

          READ(input_pid, *, IOSTAT=iostat) values_t1(1:ntoread)
        END DO

        ! take the last value in the file in case its time values don't cover
        ! the full simulation period
        IF (iostat .NE. 0) THEN
          weight_t1 = 1.0
        ELSE
          weight_t1 = (time_cur - values_t0(1)) / (values_t1(1) - values_t0(1))
        ENDIF

      ENDIF

      CLOSE(input_pid)

      values_cur = weight_t1         * values_t1 + &
                   (1.0 - weight_t1) * values_t0

      species(1:nvar) = names(ntime+nanct+1:ntime+nanct+nvar)
      values(1:nvar)  = values_cur(ntime+nanct+1:ntime+nanct+nvar)

      IF (PRESENT(nanc)) nanc = nanct

      IF (nanct .GT. 0) THEN
        IF (PRESENT(ancvalues)) THEN
          ancvalues(1:nanct) = values_cur(ntime+1:ntime+nanct)
        ENDIF
        IF (PRESENT(ancspecies)) THEN
          ancspecies(1:nanct) = names(ntime+1:ntime+nanct)
        ENDIF
      ENDIF

      END SUBROUTINE readVals
! ****************************************************************************


! ****************************************************************************
! readInput - Reads values from .csv, applies function
!
! "t"<something> is for "target" (i.e. the vector to apply "action" over)
! "r"<something> is for "read" from file
! "u"<something> is for "useable" (i.e., exists in the "target" vector)
! "m"<something> is for "missing" (i.e., asked for, but not found in file)
!
! ****************************************************************************
      SUBROUTINE readInput (filename, nt, tnames, tvalues, &
                            action, lzero, quiet)

      CHARACTER(LEN=*), INTENT(IN)  :: filename
      INTEGER, INTENT(IN)           :: nt
      REAL(KIND=dp)                 :: tvalues(nt)
      CHARACTER(LEN=15), INTENT(IN) :: tnames(nt)
      REAL(KIND=dp), OPTIONAL       :: action  ! function to be called
      LOGICAL, INTENT(IN), OPTIONAL :: lzero   ! call action with 0.0
                                               ! for 'missing' species
      LOGICAL, INTENT(IN), OPTIONAL :: quiet

      INTEGER                       :: nr, nu, nm, i, k, m, &
                                       iostat, ierrstat
      CHARACTER(LEN=15)             :: rnames(max_species_boxmox), &
                                       unames(max_species_boxmox), &
                                       mnames(max_species_boxmox)
      REAL(KIND=dp)                 :: rvalues(max_species_boxmox)
      LOGICAL                       :: is_used, first, &
                                       talkative

      CALL readVals(filename, nr, rnames, rvalues, iostat)

      IF (iostat .NE. 0) RETURN

      talkative = .TRUE.
      IF (PRESENT(quiet)) talkative = .NOT. quiet

      unames(:) = "none"
      nu        = 0

      first  = lverbose

      DO i = 1, nr
        DO k = 1, nt
          IF ( strEqual(rnames(i), tnames(k)) ) THEN

            IF (PRESENT(action)) THEN
              tvalues(k) = action( tvalues(k), k, rvalues(i), first)
            ELSE
              tvalues(k) = rvalues(i)
              IF ( (lverbose .OR. first) .AND. talkative ) THEN
                WRITE(*,"(6X,A15,1X,E12.6)") tnames(k), tvalues(k)/CFACTOR
              ENDIF
            ENDIF

            first = .FALSE.

            nu = nu + 1
            unames(nu) = TRIM(ADJUSTL(rnames(i)))
          ENDIF
        ENDDO
      ENDDO

      IF (PRESENT(lzero)) THEN
        IF (lzero) THEN
          DO k = 1, nt
            is_used = .FALSE.
            DO i = 1, nu
              IF ( strEqual(unames(i), tnames(k)) ) is_used = .TRUE.
            ENDDO
            IF (.NOT. is_used) THEN
              tvalues(k) = action( tvalues(k), k, 0.0, first)
              first = .FALSE.
            ENDIF
          ENDDO
        ENDIF
      ENDIF

      nm = nr - nu
      mnames(:) = ""

      IF (nm > 0) THEN
        m = 1
        DO i = 1, nr
          is_used = .FALSE.
          DO k = 1, nu
            IF ( strEqual(rnames(i), unames(k)) ) is_used = .TRUE.
          ENDDO

          IF (.NOT. is_used) THEN
            mnames(m) = TRIM(ADJUSTL(rnames(i)))
            m = m + 1
          ENDIF
        ENDDO
      ENDIF

      IF ( ((TIME .EQ. TSTART) .OR. lverbose) .AND. talkative) THEN
        IF (nu > 0) CALL writespecstring(nu, "useable:", unames, 6)
        IF (nm > 0) CALL writespecstring(nm, "unknown: ", mnames, 6)
      ENDIF

      END SUBROUTINE readInput
! ****************************************************************************


! ****************************************************************************
! InitialCondition - Sets initial values
! ****************************************************************************
      REAL(KIND=dp) FUNCTION InitialCondition(x, k, y, lfirst)

      REAL(KIND=dp), INTENT(IN) :: x, y
      INTEGER, INTENT(IN)       :: k
      LOGICAL, INTENT(IN)       :: lfirst

      InitialCondition = y * CFACTOR

      IF (lverbose) THEN
        IF (lfirst) WRITE(*,"(6X,A15,1X,A12,1X)") "species        ", "IC"
        WRITE(*,"(6X,A15,1X,E12.6,1X)") SPC_NAMES(k), InitialCondition/CFACTOR
      ENDIF

      END FUNCTION InitialCondition
! ****************************************************************************


! ****************************************************************************
! Deposition - Applies deposition velocities
! ****************************************************************************
      REAL(KIND=dp) FUNCTION Deposition(x, k, y, lfirst)

      REAL(KIND=dp), INTENT(IN) :: x, y
      INTEGER, INTENT(IN)       :: k
      LOGICAL, INTENT(IN)       :: lfirst

      REAL(KIND=dp)             :: dCdt

      Deposition = x

      ! we don't remove fixed variables
      IF (k .LT. NFIXST) THEN
        ! we assume deposition velocities in comparable units
        ! (i.e. cm s-1)
        ! dxdydz(3) == DZ
        dCdt = - y * x / (dxdydz(3) * 100.0)
        Deposition = MAX( 0.0, x + DT * dCdt )

        IF (lverbose) THEN
          IF (lfirst) WRITE(*,"(6X,A15,1X,3(A12,1X))") "species        ", &
            "before", "v_dep", "after"
          WRITE(*,"(6X,A15,1X,3(E12.6,1X))") SPC_NAMES(k), x/CFACTOR, y, &
                                            Deposition/CFACTOR
        ENDIF

      ENDIF

      END FUNCTION Deposition
! ****************************************************************************


! ****************************************************************************
! Emission - Adds emissions
! ****************************************************************************
      REAL(KIND=dp) FUNCTION Emission(x, k, y, lfirst)

      REAL(KIND=dp), INTENT(IN) :: x, y
      INTEGER, INTENT(IN)       :: k
      LOGICAL, INTENT(IN)       :: lfirst

      REAL(KIND=dp)             :: cur

      ! we assume emissions to be a flux in units comparable to the
      ! concentration values used in KPP
      ! e.g. molecules / m2 / s vs molecules / m3, or
      !      moles / cm2 / s vs moles / cm3, ...
      ! the time dimensions unit has to be seconds!

      ! dxdydz(3) is PBL height in m
      cur = y * DT / (dxdydz(3) * 100.0)

      Emission      = x + cur

      IF (lverbose) THEN
        IF (lfirst) WRITE(*,"(6X,A15,1X,3(A12,1X))") "species        ", &
          "before", "tendency", "after"
        WRITE(*,"(6X,A15,1X,3(E12.6,1X))") SPC_NAMES(k), x/CFACTOR, &
                                          cur/DT/CFACTOR, Emission/CFACTOR
      ENDIF

      END FUNCTION Emission
! ****************************************************************************


! ****************************************************************************
! Mix - calculate concentration after turbulent mixing
! ****************************************************************************
      REAL(KIND=dp) FUNCTION Mix(x, k, y, lfirst)

      REAL(KIND=dp), INTENT(IN) :: x, y
      REAL(KIND=dp)             :: frac_bg
      INTEGER, INTENT(IN)       :: k
      LOGICAL, INTENT(IN)       :: lfirst

      Mix = x

      frac_bg = Kturb * DT

      IF (lfirst .AND. lverbose) THEN
        WRITE(*,"(6X,A,1X,F6.3)") "fractional amount of background air", &
                                   frac_bg
        WRITE(*,"(6X,A15,1X,3(A12,1X))") "species        ", &
        "before", "background", "after"
      ENDIF

      ! we don't mix fixed variables
      IF (k .LT. NFIXST) THEN

        ! seems like we found something
        Mix =        frac_bg  * (y * CFACTOR) + &
              (1.0 - frac_bg) * x

        IF (lverbose) THEN
          WRITE(*,"(6X,A15,1X,3(E12.6,1X))") SPC_NAMES(k), x/CFACTOR, &
                                            y, Mix/CFACTOR
        ENDIF

      ENDIF

      END FUNCTION Mix
! ****************************************************************************


! ****************************************************************************
! UpdateKturb - calculate Kturb for mixing
! ****************************************************************************
      SUBROUTINE UpdateKturb()

      Kturb = 0.0
      SELECT CASE (iturb)
        CASE (TURB_OPT_NONE)
          Kturb = 0.0
        CASE (TURB_OPT_SPECIFIED)
          CALL UpdateKturbSpecified()
        CASE (TURB_OPT_VOLUME)
          CALL UpdateKturbVolume()
        CASE (TURB_OPT_TRACER)
          CALL UpdateKturbTracer()
        CASE DEFAULT
          WRITE(*,"(A,1X,I1)") "Unknown mixing option", iturb
      END SELECT

      END SUBROUTINE UpdateKturb
! ****************************************************************************


! ****************************************************************************
! UpdateKturbSpecified - calc. Kturb for mixing using specified Kturb
! ****************************************************************************
      SUBROUTINE UpdateKturbSpecified()

      INTEGER                       :: nr, na, i, iostat, ierrstat
      CHARACTER(LEN=15)             :: rnames(max_species_boxmox), &
                                       anames(max_species_boxmox)
      REAL(KIND=dp)                 :: rvalues(max_species_boxmox), &
                                       avalues(max_species_boxmox)
      LOGICAL                       :: lfound

      Kturb = 0.0

      CALL readVals('Background.csv', &
                    nr, rnames, rvalues, iostat, &
                    na, anames, avalues, lquiet=.TRUE.)

      IF (iostat .NE. 0) THEN
        RETURN
      ENDIF

      lfound = .FALSE.
      DO i = 1, na
        IF ( strEqual(anames(i), 'Kturb') ) THEN
          Kturb         = avalues(i)
          lfound = .TRUE.
        ENDIF
      ENDDO

      IF (.NOT. lfound) THEN
        WRITE(*,*) "Kturb variable not found in Background.csv"
        STOP
      ENDIF

      END SUBROUTINE UpdateKturbSpecified
! ****************************************************************************


! ****************************************************************************
! UpdateKturbVolume - calc. Kturb for mixing using box volume changes
! ****************************************************************************
      SUBROUTINE UpdateKturbVolume()

      REAL(KIND=dp)             :: vol, vol_old

      vol     = dxdydz(1)     * dxdydz(2)     * dxdydz(3)
      vol_old = dxdydz_old(1) * dxdydz_old(2) * dxdydz_old(3)

      ! maybe we have a mixing case where the box volume changes?
      ! increasing volume - dilute against background
      ! decreasing volume - ignore (?) for now

      ! in the PBL case, analogue:
      ! maybe we changed the PBLH? update concentration vector
      ! rising BL - mixing against zero air or bg values
      ! shrinking BL - creates residual layer, decouples

      IF (vol > vol_old) THEN
        Kturb = ( (vol - vol_old) / vol_old ) / DT
      ENDIF

      END SUBROUTINE UpdateKturbVolume
! ****************************************************************************


! ****************************************************************************
! UpdateKturbTracer - calc. Kturb for mixing using tracer (t, tname) conc.
!                         changes
! ****************************************************************************
      SUBROUTINE UpdateKturbTracer()

      INTEGER                       :: nr, na, i, iostat, ierrstat
      CHARACTER(LEN=15)             :: rnames(max_species_boxmox), &
                                       anames(max_species_boxmox), &
                                       tracer_name
      REAL(KIND=dp)                 :: rvalues(max_species_boxmox), &
                                       avalues(max_species_boxmox), &
                                       prev, obs, bg, frac_bg
      LOGICAL                       :: lfound_prev, lfound_bg

      Kturb = 0.0

      CALL readVals('Background.csv', &
                    nr, rnames, rvalues, iostat, &
                    na, anames, avalues, lquiet=.TRUE.)

      IF (iostat .NE. 0) THEN
        RETURN
      ENDIF

      ! ASSUMPTION: first ancillary species is mixing tracer
      tracer_name = anames(1)
      obs         = avalues(1) * CFACTOR

      prev        = 0.0
      lfound_prev = .FALSE.
      DO i = 1, NSPEC
        IF ( strEqual(SPC_NAMES(i), tracer_name) ) THEN
          prev        = C(i)
          lfound_prev = .TRUE.
        ENDIF
      ENDDO

      bg        = 0.0
      lfound_bg = .FALSE.
      DO i = 1, nr
        IF ( strEqual(rnames(i), tracer_name) ) THEN
          bg        = rvalues(i) * CFACTOR
          lfound_bg = .TRUE.
        ENDIF
      ENDDO

      IF (.NOT. lfound_bg) THEN
        WRITE(*,*) "Tracer ", TRIM(tracer_name), &
                   " not found in Background.csv"
        STOP
      ENDIF
      IF (.NOT. lfound_prev) THEN
        WRITE(*,*) "Tracer ", TRIM(tracer_name), &
                   " not found in mechanism"
        STOP
      ENDIF

      IF ( ABS( bg - prev ) < 1e-10 ) THEN
        Kturb = 0.0
      ELSE
        frac_bg = ( obs - prev ) / ( bg - prev )
        frac_bg = MIN( frac_bg, 1.0 )
        frac_bg = MAX( frac_bg, 0.0 )

        Kturb       = frac_bg / DT
      ENDIF

      END SUBROUTINE UpdateKturbTracer
! ****************************************************************************


! ****************************************************************************
! UpdatePhotolysisRates - well, it does what it seems to do
! ****************************************************************************
      SUBROUTINE UpdatePhotolysisRates()

      INTEGER             :: i, iostat, nr, jidx
      CHARACTER(LEN=15)   :: rnames(max_species_boxmox)
      REAL(KIND=dp)       :: rvalues(max_species_boxmox)

      INTEGER             :: j, jidx2

      rnames(:)  = "               "
      rvalues(:) = 0.0

      jvalues(:) = 0.0

      CALL readVals('PhotolysisRates.csv', nr, rnames, rvalues, iostat)

      IF (iostat .NE. 0) THEN
        RETURN
      ENDIF

      DO i = 1, nr
        IF (lverbose) THEN
          WRITE(*,"(6X,A15,1X,E12.6)")   rnames(i), rvalues(i)
        ENDIF
        jidx          = hash(rnames(i))
        jvalues(jidx) = rvalues(i)
      ENDDO

      ! first time step only: make sure hashing actually made unique hashes
      IF (TIME .EQ. TSTART) THEN
        DO i = 1, nr
          jidx          = hash(rnames(i))
          DO j = 1, nr
            jidx2          = hash(rnames(j))
            IF ((jidx .EQ. jidx2) .AND. .NOT. (i .EQ. j)) THEN
              WRITE(*,*) "Photolysis rate label hash collision:"
              WRITE(*,*) rnames(i), "(idx ",i,", hash ", jidx, ") == ", &
                         rnames(j), "(idx ",j,", hash ", jidx2,")"
              IF ( strEqual( rnames(i), rnames(j) ) ) THEN
                WRITE(*,*) "Remove duplicate label in PhotolysisRates.csv!"
              ELSE
                WRITE(*,*) "Please contact the BOXMOX developers!"
              ENDIF
              STOP
            ENDIF
          ENDDO
        ENDDO

        CALL writespecstring(nr, "known:", rnames, 6)

      ENDIF

      END SUBROUTINE UpdatePhotolysisRates
! ****************************************************************************


! ****************************************************************************
! UpdateAerosol - update aerosol properties
! ****************************************************************************
      SUBROUTINE UpdateAerosol()

      INTEGER, PARAMETER :: naer = 2, sad_idx = 1, r_idx = 2
      REAL(KIND=dp)      :: aer(naer), errval = -987654321.0
      CHARACTER(LEN=15)  :: aernames(naer)
      INTEGER            :: iostat, i

      aernames(sad_idx) = "sad"; aer(sad_idx) = errval
      aernames(r_idx)   = "r";   aer(r_idx)   = errval

      CALL readInput('Aerosol.csv', naer, aernames, aer, quiet=lverbose)

      IF (aer(sad_idx) .NE. errval) AER_SAD  = aer(sad_idx)
      IF (aer(r_idx)   .NE. errval) AER_R    = aer(r_idx)

      IF (lverbose) THEN
        WRITE(*,"(6X,A15,1X,F8.2)")   "SAD            ", AER_SAD
        WRITE(*,"(6X,A15,1X,F8.2)")   "R              ", AER_R
      ENDIF

      END SUBROUTINE UpdateAerosol
! ****************************************************************************


! ****************************************************************************
! UpdateEnvironment - Read environmental values
! ****************************************************************************
      SUBROUTINE UpdateEnvironment()

      INTEGER, PARAMETER :: nenv = 4, temp_idx = 1, pres_idx = 2, &
                            pblh_idx = 3, vol_idx = 4
      REAL(KIND=dp)      :: env(nenv), errval = -987654321.0
      CHARACTER(LEN=15)  :: envnames(nenv)
      INTEGER            :: iostat, i

      REAL(KIND=dp)      :: temp_old, pres_old, pblh, vol

      envnames(temp_idx)  = "TEMP";  env(temp_idx)  = errval
      envnames(pres_idx)  = "PRES";  env(pres_idx)  = errval
      envnames(pblh_idx)  = "PBLH";  env(pblh_idx)  = errval
      envnames(vol_idx)   = "VOL";   env(vol_idx)   = errval

      CALL readInput('Environment.csv', nenv, envnames, env, quiet=lverbose)

      temp_old       = temp
      pres_old       = pres

      dxdydz_old(:)  = dxdydz(:)

      ! those are just pseudo variables for dxdydz properties
      pblh           = dxdydz(3)
      vol            = dxdydz(1) * dxdydz(2) * dxdydz(3)

      IF (env(pblh_idx)  .NE. errval .AND. env(vol_idx)   .NE. errval) &
      THEN
        WRITE(*,*) "Cannot work with changing PBLH and VOL at the same time."
        STOP
      ENDIF

      IF (env(temp_idx)  .NE. errval) temp  = env(temp_idx)
      IF (env(pres_idx)  .NE. errval) pres  = env(pres_idx)
      IF (env(vol_idx)   .NE. errval) vol   = env(vol_idx)

      IF (env(pblh_idx)  .NE. errval) THEN
        pblh  = env(pblh_idx)
        dxdydz(3) = pblh
      ENDIF

      ! Adjust species vector
      DO i = 1, NSPEC
        C(i) = C(i) * (temp_old * pres) / (temp * pres_old)
      ENDDO

      ! recalculate unit conversion
      CALL calcCFACTOR()

      IF (lverbose) THEN
        WRITE(*,"(6X,A15,1X,F8.2)")   "TEMP           ", temp
        WRITE(*,"(6X,A13,1X,F10.2)")  "PRES           ", pres
        WRITE(*,"(6X,A15,1X,F8.2)")   "PBLH           ", dxdydz(3)
        WRITE(*,"(6X,A15,1X,F8.2)")   "VOL            ", &
                                              dxdydz(1)*dxdydz(2)*dxdydz(3)
        WRITE(*,*)                    ""
        WRITE(*,"(6X,A7,1X,E16.8)")   "CFACTOR", CFACTOR
      ENDIF

      END SUBROUTINE UpdateEnvironment
! ****************************************************************************


! ****************************************************************************
! InitFixNOx - Initialize the NOx fixing data (for steady state calculations)
! ****************************************************************************
      SUBROUTINE InitFixNOx(n, names, values)

      INTEGER, INTENT(IN)           :: n
      REAL(KIND=dp)                 :: values(n)
      CHARACTER(LEN=15), INTENT(IN) :: names(n)

      INTEGER                       :: i, nfound_tracers

      initial_total_nox = 0.0
      nfound_tracers    = 0
      DO i = 1, n
        IF ( strEqual(names(i), 'NO             ') ) THEN
          initial_total_nox = initial_total_nox + values(i)
          nfound_tracers    = nfound_tracers + 1
        ENDIF
        IF ( strEqual(names(i), 'NO2            ') ) THEN
          initial_total_nox = initial_total_nox + values(i)
          nfound_tracers    = nfound_tracers + 1
        ENDIF
      ENDDO

      IF (.NOT. nfound_tracers == 2) THEN
        WRITE(*,*) "NOx fixing failed, found (only) ", &
                   nfound_tracers, " species."
        STOP
      ENDIF

      END SUBROUTINE InitFixNOx
! ****************************************************************************


! ****************************************************************************
! FixNOx - keep total NOx at initial levels (for steady state calculations)
! ****************************************************************************
      SUBROUTINE FixNOx(n, names, values)

      INTEGER, INTENT(IN)           :: n
      REAL(KIND=dp), INTENT(INOUT)  :: values(n)
      CHARACTER(LEN=15), INTENT(IN) :: names(n)

      INTEGER                       :: i, nfound_tracers
      REAL(KIND=dp)                 :: current_total_nox, scaleFactor

      current_total_nox = 0.0
      nfound_tracers = 0
      DO i = 1, n
        IF ( strEqual(names(i), 'NO             ') ) THEN
          current_total_nox      = current_total_nox + values(i)
          nfound_tracers = nfound_tracers + 1
        ENDIF
        IF ( strEqual(names(i), 'NO2            ') ) THEN
          current_total_nox      = current_total_nox + values(i)
          nfound_tracers = nfound_tracers + 1
        ENDIF
      ENDDO

      IF (.NOT. nfound_tracers == 2) THEN
        WRITE(*,*) "NOx fixing failed, found (only) ", &
                   nfound_tracers, " species."
        STOP
      ENDIF

      scaleFactor = initial_total_nox
      IF (current_total_nox > 0.0) THEN
        scaleFactor = initial_total_nox / current_total_nox
      ENDIF

      DO i = 1, n
        IF ( strEqual(names(i), 'NO             ') ) THEN
          values(i) = values(i) * scaleFactor
        ENDIF
        IF ( strEqual(names(i), 'NO2            ') ) THEN
          values(i) = values(i) * scaleFactor
        ENDIF
      ENDDO

      END SUBROUTINE FixNOx
! ****************************************************************************


! End INLINED Utility Functions

! Utility Functions from KPP_HOME/util/util
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! UTIL - Utility functions
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! ****************************************************************
!
! InitSaveData - Opens the data file for writing
!   Parameters :
!
! ****************************************************************

      SUBROUTINE InitSaveData ()

      USE CB05_CSIRO_Parameters
!* ----------------------------------------------------
!* BOXMOX extension - Jerome Barre - 20150410
!*                    Christoph Knote - 20151030
!* add header to output species and sparse jacobian indices
!* and Hessian
      USE CB05_CSIRO_JacobianSP
#if 0 > 0
      USE CB05_CSIRO_HessianSP
#endif
!* ----------------------------------------------------

!* ----------------------------------------------------
!* BOXMOX extension - Christoph Knote - 20140610
!* add header to output files and open rates file
      INTEGER     :: i
!* ----------------------------------------------------

!* ----------------------------------------------------
!* BOXMOX extension - Christoph Knote - 20161129
!* renamed concentration output
!      open(10, file='CB05_CSIRO.dat')
      open(10, file='CB05_CSIRO.conc')
!* ----------------------------------------------------

!* ----------------------------------------------------
!* BOXMOX extension - Christoph Knote - 20140610
!*                    Jerome Barre - 20150410
!*                    Christoph Knote - 20151030
!* add header to output files and open rates file

      ! add header to output file
      WRITE(10, "(A24,10000(1X,A24))") "time", ADJUSTR(SPC_NAMES(LOOKAT))

      open(11, file='CB05_CSIRO.rates')

      WRITE(11, "(I5)") NREACT

      DO i = 1, NREACT
        WRITE(11, "(A)") EQN_NAMES(i)
      ENDDO

      ! add header to output file
      WRITE(11, "(A24,10000(1X,I24))") "time", (/(i, i=1,NREACT, 1)/)

!* write header to Jacobian output
      open(12, file='CB05_CSIRO.jacobian')

!* species names
      WRITE(12, "(10000(1X,A24))") ADJUSTR(SPC_NAMES(LOOKAT))
!* species index
      WRITE(12, "(10000(1X,I5))") LOOKAT
!* Jacobian non zero row index
      WRITE(12, "(100000(1X,I5))") (LU_IROW(i), i=1,LU_NONZERO)
!* Jacobian non zero column index
      WRITE(12, "(100000(1X,I5))") (LU_ICOL(i), i=1,LU_NONZERO)

!* write header to Hessian output
#if 0 > 0
      open(13, file='CB05_CSIRO.hessian')

!* species names
      WRITE(13, "(10000(1X,A24))") ADJUSTR(SPC_NAMES(LOOKAT))
!* species index
      WRITE(13, "(10000(1X,I5))") LOOKAT
!* Hessian i, j, k species
      WRITE(13, "(100000(1X,I5))") (IHESS_I(i), i=1,NHESS)
      WRITE(13, "(100000(1X,I5))") (IHESS_J(i), i=1,NHESS)
      WRITE(13, "(100000(1X,I5))") (IHESS_K(i), i=1,NHESS)
#endif
!* ----------------------------------------------------

      END SUBROUTINE InitSaveData

! End of InitSaveData function
! ****************************************************************

! ****************************************************************
!
! SaveData - Write LOOKAT species in the data file
!   Parameters :
!
! ****************************************************************

      SUBROUTINE SaveData ()

      USE CB05_CSIRO_Global
      USE CB05_CSIRO_Monitor
!* ----------------------------------------------------
!* BOXMOX extension - Jerome Barre - 20150410
!*                    Christoph Knote - 20151030
      USE CB05_CSIRO_Jacobian
#if 0 > 0
      USE CB05_CSIRO_Hessian
#endif
!* ----------------------------------------------------

!* ----------------------------------------------------
!* BOXMOX extension - Christoph Knote - 20151030
      REAL(KIND=dp)      :: JVS(LU_NONZERO)           ! d(dc/dt)i/dcj
#if 0 > 0
      REAL(KIND=dp)      :: HESS(NHESS)               !
#endif
      INTEGER i

      WRITE(10,999) (TIME-TSTART)/3600.D0,  &
                   (C(LOOKAT(i))/CFACTOR, i=1,NLOOKAT)

!* ----------------------------------------------------
!* BOXMOX extension - Christoph Knote - 20140610
!*                    Jerome Barre - 20150410
!*                    Christoph Knote - 20151030
!* write out reaction rates
      WRITE(11,999) (TIME-TSTART)/3600.D0,  &
                   (RCONST(i), i=1,NREACT)
!* call Jacobian function to output Jacobian values
      CALL Jac_SP( VAR, FIX, RCONST, JVS )
      WRITE(12,9999) (JVS(i), i=1,LU_NONZERO)
!* call Hessian function to output Hessian values
#if 0 > 0
      CALL Hessian( VAR, FIX, RCONST, HESS )
      WRITE(13,9999) (HESS(i), i=1,NHESS)
#endif
!* ---------------------------------------------------

!* ----------------------------------------------------
!* BOXMOX extension - Christoph Knote - 20140610
!* might have more than 100 output variables
! 999   FORMAT(E24.16,100(1X,E24.16))
999    FORMAT(E24.16,10000(1X,E24.16E3))
9999   FORMAT(E24.16,100000(1X,E24.16E3))
!* ----------------------------------------------------

      END SUBROUTINE SaveData

! End of SaveData function
! ****************************************************************

! ****************************************************************
!
! CloseSaveData - Close the data file
!   Parameters :
!
! ****************************************************************

      SUBROUTINE CloseSaveData ()

      USE CB05_CSIRO_Parameters

      CLOSE(10)

!* ----------------------------------------------------
!* BOXMOX extension - Christoph Knote - 20140610
!*                    Jerome Barre - 20150410
!                     Christoph Knote - 20153010
!* close rates file
      CLOSE(11)
!* close Jacobian file
      CLOSE(12)
!* close Hessian file
#if 0 > 0
      CLOSE(13)
#endif
!* ---------------------------------------------------

      END SUBROUTINE CloseSaveData

! End of CloseSaveData function
! ****************************************************************

! ****************************************************************
!
! GenerateMatlab - Generates MATLAB file to load the data file
!   Parameters :
!                It will have a character string to prefix each
!                species name with.
!
! ****************************************************************

      SUBROUTINE GenerateMatlab ( PREFIX )

      USE CB05_CSIRO_Parameters
      USE CB05_CSIRO_Global
      USE CB05_CSIRO_Monitor


      CHARACTER(LEN=8) PREFIX
      INTEGER i

      open(20, file='CB05_CSIRO.m')
      write(20,*) 'load CB05_CSIRO.dat;'
      write(20,990) PREFIX
990   FORMAT(A1,'c = CB05_CSIRO;')
      write(20,*) 'clear CB05_CSIRO;'
      write(20,991) PREFIX, PREFIX
991   FORMAT(A1,'t=',A1,'c(:,1);')
      write(20,992) PREFIX
992   FORMAT(A1,'c(:,1)=[];')

      do i=1,NLOOKAT
        write(20,993) PREFIX, SPC_NAMES(LOOKAT(i)), PREFIX, i
993     FORMAT(A1,A6,' = ',A1,'c(:,',I2,');')
      end do

      CLOSE(20)

      END SUBROUTINE GenerateMatlab

! End of GenerateMatlab function
! ****************************************************************


! ****************************************************************
!                            
! tag2num - convert equation tags to kpp reaction number
!   Arguments :
!      id        - string with the equation tag
!
! ****************************************************************

ELEMENTAL INTEGER FUNCTION tag2num ( id )

  USE CB05_CSIRO_Monitor, ONLY: EQN_TAGS

  CHARACTER(LEN=*), INTENT(IN) :: id
  INTEGER i

  tag2num = 0 ! mz_rs_20050115
  DO i = 1, SIZE(EQN_TAGS)
    IF (TRIM(EQN_TAGS(i)) == TRIM(id)) THEN
      tag2num = i ! mz_rs_20050115
      EXIT
    ENDIF
  END DO

END FUNCTION tag2num

! End of tag2num function
! ****************************************************************

! End Utility Functions from KPP_HOME/util/util
! End of UTIL function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Shuffle_user2kpp - function to copy concentrations from USER to KPP
!   Arguments :
!      V_USER    - Concentration of variable species in USER's order
!      V         - Concentrations of variable species (local)
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Shuffle_user2kpp ( V_USER, V )

! V_USER - Concentration of variable species in USER's order
  REAL(kind=dp) :: V_USER(NVAR)
! V - Concentrations of variable species (local)
  REAL(kind=dp) :: V(NVAR)

  V(61) = V_USER(1)
  V(62) = V_USER(2)
  V(63) = V_USER(3)
  V(53) = V_USER(4)
  V(56) = V_USER(5)
  V(14) = V_USER(6)
  V(60) = V_USER(7)
  V(54) = V_USER(8)
  V(12) = V_USER(9)
  V(34) = V_USER(10)
  V(16) = V_USER(11)
  V(23) = V_USER(12)
  V(17) = V_USER(13)
  V(59) = V_USER(14)
  V(48) = V_USER(15)
  V(51) = V_USER(16)
  V(27) = V_USER(17)
  V(52) = V_USER(18)
  V(50) = V_USER(19)
  V(49) = V_USER(20)
  V(44) = V_USER(21)
  V(36) = V_USER(22)
  V(57) = V_USER(23)
  V(30) = V_USER(24)
  V(10) = V_USER(25)
  V(26) = V_USER(26)
  V(19) = V_USER(27)
  V(58) = V_USER(28)
  V(6) = V_USER(29)
  V(20) = V_USER(30)
  V(21) = V_USER(31)
  V(55) = V_USER(32)
  V(22) = V_USER(33)
  V(28) = V_USER(34)
  V(41) = V_USER(35)
  V(38) = V_USER(36)
  V(40) = V_USER(37)
  V(8) = V_USER(38)
  V(43) = V_USER(39)
  V(29) = V_USER(40)
  V(45) = V_USER(41)
  V(33) = V_USER(42)
  V(37) = V_USER(43)
  V(15) = V_USER(44)
  V(46) = V_USER(45)
  V(47) = V_USER(46)
  V(4) = V_USER(47)
  V(1) = V_USER(48)
  V(2) = V_USER(49)
  V(11) = V_USER(50)
  V(9) = V_USER(51)
  V(42) = V_USER(52)
  V(13) = V_USER(56)
  V(5) = V_USER(57)
  V(24) = V_USER(58)
  V(31) = V_USER(59)
  V(32) = V_USER(60)
  V(18) = V_USER(61)
  V(35) = V_USER(62)
  V(25) = V_USER(63)
      
END SUBROUTINE Shuffle_user2kpp

! End of Shuffle_user2kpp function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Shuffle_kpp2user - function to restore concentrations from KPP to USER
!   Arguments :
!      V         - Concentrations of variable species (local)
!      V_USER    - Concentration of variable species in USER's order
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE Shuffle_kpp2user ( V, V_USER )

! V - Concentrations of variable species (local)
  REAL(kind=dp) :: V(NVAR)
! V_USER - Concentration of variable species in USER's order
  REAL(kind=dp) :: V_USER(NVAR)

  V_USER(1) = V(61)
  V_USER(2) = V(62)
  V_USER(3) = V(63)
  V_USER(4) = V(53)
  V_USER(5) = V(56)
  V_USER(6) = V(14)
  V_USER(7) = V(60)
  V_USER(8) = V(54)
  V_USER(9) = V(12)
  V_USER(10) = V(34)
  V_USER(11) = V(16)
  V_USER(12) = V(23)
  V_USER(13) = V(17)
  V_USER(14) = V(59)
  V_USER(15) = V(48)
  V_USER(16) = V(51)
  V_USER(17) = V(27)
  V_USER(18) = V(52)
  V_USER(19) = V(50)
  V_USER(20) = V(49)
  V_USER(21) = V(44)
  V_USER(22) = V(36)
  V_USER(23) = V(57)
  V_USER(24) = V(30)
  V_USER(25) = V(10)
  V_USER(26) = V(26)
  V_USER(27) = V(19)
  V_USER(28) = V(58)
  V_USER(29) = V(6)
  V_USER(30) = V(20)
  V_USER(31) = V(21)
  V_USER(32) = V(55)
  V_USER(33) = V(22)
  V_USER(34) = V(28)
  V_USER(35) = V(41)
  V_USER(36) = V(38)
  V_USER(37) = V(40)
  V_USER(38) = V(8)
  V_USER(39) = V(43)
  V_USER(40) = V(29)
  V_USER(41) = V(45)
  V_USER(42) = V(33)
  V_USER(43) = V(37)
  V_USER(44) = V(15)
  V_USER(45) = V(46)
  V_USER(46) = V(47)
  V_USER(47) = V(4)
  V_USER(48) = V(1)
  V_USER(49) = V(2)
  V_USER(50) = V(11)
  V_USER(51) = V(9)
  V_USER(52) = V(42)
  V_USER(56) = V(13)
  V_USER(57) = V(5)
  V_USER(58) = V(24)
  V_USER(59) = V(31)
  V_USER(60) = V(32)
  V_USER(61) = V(18)
  V_USER(62) = V(35)
  V_USER(63) = V(25)
      
END SUBROUTINE Shuffle_kpp2user

! End of Shuffle_kpp2user function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! GetMass - compute total mass of selected atoms
!   Arguments :
!      CL        - Concentration of all species (local)
!      Mass      - value of mass balance
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

SUBROUTINE GetMass ( CL, Mass )

! CL - Concentration of all species (local)
  REAL(kind=dp) :: CL(NSPEC)
! Mass - value of mass balance
  REAL(kind=dp) :: Mass(1)

      
END SUBROUTINE GetMass

! End of GetMass function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



END MODULE CB05_CSIRO_Util

