MODULE lmdz_check_isotopes
  IMPLICIT NONE; PRIVATE
  PUBLIC check_isotopes_seq

CONTAINS


  SUBROUTINE check_isotopes_seq(q, ip1jmp1, err_msg)
    USE lmdz_strings, ONLY: maxlen, msg, strIdx, strStack, int2str, real2str
    USE lmdz_infotrac, ONLY: nqtot, niso, nphas, isotope, isoCheck, iqIsoPha, isoSelect, &
            ntiso, iH2O, nzone, tracers, isoName, itZonIso, getKey

    USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
    IMPLICIT NONE

    REAL, INTENT(INOUT) :: q(ip1jmp1, llm, nqtot)
    INTEGER, INTENT(IN) :: ip1jmp1
    CHARACTER(LEN = *), INTENT(IN) :: err_msg    !--- Error message to display
    CHARACTER(LEN = maxlen) :: modname, msg1, nm(2)
    INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
    INTEGER, ALLOCATABLE :: ix(:)
    REAL, ALLOCATABLE, SAVE :: tnat(:)
    REAL :: xtractot, xiiso, deltaD, q1, q2
    REAL, PARAMETER :: borne = 1e19, &
            errmax = 1e-8, &       !--- Max. absolute error
            errmaxrel = 1e-3, &       !--- Max. relative error
            qmin = 1e-11, &
            deltaDmax = 1000.0, &
            deltaDmin = -999.0, &
            ridicule = 1e-12
    INTEGER, SAVE :: iso_eau, iso_HDO, iso_O18, &
            iso_O17, iso_HTO
    LOGICAL, SAVE :: first = .TRUE.
    LOGICAL, PARAMETER :: tnat1 = .TRUE.

    modname = 'check_isotopes'
    IF(.NOT.isoCheck)    RETURN                   !--- No need to check => finished
    IF(isoSelect('H2O')) RETURN                   !--- No H2O isotopes group found
    IF(niso == 0)        RETURN                   !--- No isotopes => finished
    IF(first) THEN
      iso_eau = strIdx(isoName, 'H216O')
      iso_HDO = strIdx(isoName, 'HDO')
      iso_O18 = strIdx(isoName, 'H218O')
      iso_O17 = strIdx(isoName, 'H217O')
      iso_HTO = strIdx(isoName, 'HTO')
      IF (tnat1) THEN
        tnat(:) = 1.0
      else
        IF(getKey('tnat', tnat)) CALL abort_gcm(modname, 'missing isotopic parameter', 1)
      endif
      first = .FALSE.
    END IF
    CALL msg('31: err_msg=' // TRIM(err_msg), modname)

    !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
    modname = 'check_isotopes:iso_verif_noNaN'
    DO ixt = 1, ntiso
      DO ipha = 1, nphas
        iq = iqIsoPha(ixt, ipha)
        DO k = 1, llm
          DO i = 1, ip1jmp1
            IF(ABS(q(i, k, iq)) < borne) CYCLE
            WRITE(msg1, '(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')TRIM(isoName(ixt)), i, k, iq, q(i, k, iq)
            CALL msg(msg1, modname)
            CALL abort_gcm(modname, 'Error with isotopes: ' // TRIM(err_msg), 1)
          END DO
        END DO
      END DO
    END DO

    !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
    modname = 'check_isotopes:iso_verif_egalite'
    ixt = iso_eau
    IF(ixt /= 0) THEN
      DO ipha = 1, nphas
        iq = iqIsoPha(ixt, ipha)
        iqpar = tracers(iq)%iqParent
        DO k = 1, llm
          DO i = 1, ip1jmp1
            q1 = q(i, k, iqpar)
            q2 = q(i, k, iq)
            !--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
            !    This would be probably required to sum from smallest to highest concentrations ; the corresponding
            !    indices vector can be computed once only (in the initializations stage), using mean concentrations.
            !              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
            IF(ABS(q1 - q2) <= errmax .OR. ABS(q1 - q2) / MAX(MAX(ABS(q1), ABS(q2)), 1e-18) <= errmaxrel) THEN
              q(i, k, iq) = q1                 !--- Bidouille pour convergence
              !                 q(i,k,tracers(iqPar)%iqDesc) = q(i,k,tracers(iqPar)%iqDesc) * q1 / q2
              CYCLE
            END IF
            CALL msg('ixt, iq = ' // TRIM(strStack(int2str([ixt, iq]))), modname)
            msg1 = '(' // TRIM(strStack(int2str([i, k]))) // ')'
            CALL msg(TRIM(tracers(iqpar)%name) // TRIM(msg1) // ' = ' // TRIM(real2str(q1)), modname)
            CALL msg(TRIM(tracers(iq)%name) // TRIM(msg1) // ' = ' // TRIM(real2str(q2)), modname)
            CALL abort_gcm(modname, 'Error with isotopes: ' // TRIM(err_msg), 1)
          END DO
        END DO
      END DO
    END IF

    !--- CHECK DELTA ANOMALIES
    modname = 'check_isotopes:iso_verif_aberrant'
    ix = [ iso_HDO, iso_O18 ]
    nm = ['deltaD  ', 'deltaO18']
    DO iiso = 1, SIZE(ix)
      ixt = ix(iiso)
      IF(ixt  == 0) CYCLE
      DO ipha = 1, nphas
        iq = iqIsoPha(ixt, ipha)
        iqpar = tracers(iq)%iqParent
        DO k = 1, llm
          DO i = 1, ip1jmp1
            q1 = q(i, k, iqpar)
            q2 = q(i, k, iq)
            !--- IMPROVEMENT in case at least one isotope is not negligible compared to the main isotopic form.
            !    This would be probably required to sum from smallest to highest concentrations ; the corresponding
            !    indices vector can be computed once only (in the initializations stage), using mean concentrations.
            !              q2 = SUM(q(i,k,tracers(iqPar)%iqDesc), DIM=3)
            IF(q2 <= qmin) CYCLE
            deltaD = (q2 / q1 / tnat(ixt) - 1.) * 1000.
            IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
            CALL msg('ixt, iq = ' // TRIM(strStack(int2str([ixt, iq]))), modname)
            msg1 = '(' // TRIM(strStack(int2str([i, k]))) // ')'
            CALL msg(TRIM(tracers(iqpar)%name) // TRIM(msg1) // ' = ' // TRIM(real2str(q1)), modname)
            CALL msg(TRIM(tracers(iq)%name) // TRIM(msg1) // ' = ' // TRIM(real2str(q2)), modname)
            CALL msg(TRIM(nm(iiso)) // TRIM(real2str(deltaD)), modname)
            CALL abort_gcm(modname, 'Error with isotopes: ' // TRIM(err_msg), 1)
          END DO
        END DO
      END DO
    END DO

    IF(nzone == 0) RETURN

    !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
    modname = 'check_isotopes:iso_verif_aberrant'
    IF(iso_eau /= 0 .AND. iso_HDO /= 0) THEN
      DO izon = 1, nzone
        ixt = itZonIso(izon, iso_HDO)
        ieau = itZonIso(izon, iso_eau)
        DO ipha = 1, nphas
          iq = iqIsoPha(ixt, ipha)
          iqeau = iqIsoPha(ieau, ipha)
          DO k = 1, llm
            DO i = 1, ip1jmp1
              q1 = q(i, k, iqeau)
              q2 = q(i, k, iq)
              IF(q2<=qmin) CYCLE
              deltaD = (q2 / q1 / tnat(iso_HDO) - 1.) * 1000.
              IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
              CALL msg('izon, ipha = ' // TRIM(strStack(int2str([izon, ipha]))), modname)
              CALL msg('ixt, ieau = ' // TRIM(strStack(int2str([ ixt, ieau]))), modname)
              msg1 = '(' // TRIM(strStack(int2str([i, k]))) // ')'
              CALL msg(TRIM(tracers(iqeau)%name) // TRIM(msg1) // ' = ' // TRIM(real2str(q1)), modname)
              CALL msg(TRIM(tracers(iq)%name) // TRIM(msg1) // ' = ' // TRIM(real2str(q2)), modname)
              CALL msg('deltaD = ' // TRIM(real2str(deltaD)), modname)
              CALL abort_gcm(modname, 'Error with isotopes: ' // TRIM(err_msg), 1)
            END DO
          END DO
        END DO
      END DO
    END IF

    !--- CHECK FOR TAGGING TRACERS CONSERVATION (PARENT AND TAGGING TRACERS SUM OVER ALL REGIONS MUST BE EQUAL)
    DO iiso = 1, niso
      DO ipha = 1, nphas
        iq = iqIsoPha(iiso, ipha)
        DO k = 1, llm
          DO i = 1, ip1jmp1
            xiiso = q(i, k, iq)
            xtractot = SUM(q(i, k, iqIsoPha(itZonIso(1:nzone, iiso), ipha)))
            IF(ABS(xtractot - xiiso) > errmax .AND. ABS(xtractot - xiiso) / MAX(MAX(ABS(xtractot), ABS(xiiso)), 1e-18) > errmaxrel) THEN
              CALL msg('Error in iso_verif_aberrant trac: ' // TRIM(err_msg))
              CALL msg('iiso, ipha = ' // TRIM(strStack(int2str([iiso, ipha]))), modname)
              CALL msg('q(' // TRIM(strStack(int2str([i, k]))) // ',:) = ' // TRIM(strStack(real2str(q(i, k, :)))), modname)
              CALL abort_gcm(modname, 'Error with isotopes: ' // TRIM(err_msg), 1)
            END IF
            IF(ABS(xtractot) <= ridicule) CYCLE
            DO izon = 1, nzone
              q(i, k, iq) = q(i, k, iq) / xtractot * xiiso !--- Bidouille pour convergence
            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE check_isotopes_seq


END MODULE lmdz_check_isotopes