SUBROUTINE check_isotopes(q, ijb, ije, err_msg)
  USE strings_mod, ONLY: strIdx, msg
  USE infotrac,    ONLY: isotope, isoSelect, iH2O, isoCheck, isoName, nqtot, niso, nitr, nzon, npha, iTraPha, iZonIso, tnat
  USE parallel_lmdz
  IMPLICIT NONE
#include "dimensions.h"
  REAL,             INTENT(INOUT) :: q(ijb_u:ije_u,llm,nqtot)
  INTEGER,          INTENT(IN)    :: ijb, ije    !--- Can be local and different from ijb_u,ije_u, for example in qminimum
  CHARACTER(LEN=*), INTENT(IN)    :: err_msg     !--- Error message to display
  CHARACTER(LEN=256) :: msg1, modname
  INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau
  INTEGER, ALLOCATABLE :: ix(:)
  REAL    :: xtractot, xiiso, deltaD, q1, q2
  REAL, PARAMETER :: borne     = 1e19,  &
                     errmax    = 1e-8,  &        !--- Max. absolute error
                     errmaxrel = 1e-8,  &        !--- Max. relative error
                     qmin      = 1e-11, &
                     deltaDmax = 200.0, &
                     deltaDmin =-999.9, &
                     ridicule  = 1e-12
  INTEGER, SAVE :: ixH2O, ixHDO, ixO18
!OMP THREADPRIVATE(ixH2O, ixHDO, ixO18)
  LOGICAL, SAVE :: first=.TRUE.
!OMP THREADPRIVATE(first)

  modname = 'check_isotopes'
  IF(first) THEN
    iH2O = -1
    IF(isoSelect('H2O')) RETURN
    ixH2O = strIdx(isoName,'H2[16]O')
    ixHDO = strIdx(isoName,'H[2]HO')
    ixO18 = strIdx(isoName,'H2[18]O')
    first = .FALSE.
  ELSE
    IF(iH2O == -1)      RETURN
    IF(isoSelect(iH2O)) RETURN
  END IF
  IF(.NOT.isoCheck .OR. niso == 0) RETURN        !--- No need to check or no isotopes => finished

  !--- CHECK FOR NaNs (FOR ALL THE ISOTOPES, INCLUDING GEOGRAPHIC TAGGING TRACERS)
  DO ixt = 1, nitr
    DO ipha = 1, npha
      iq = iTraPha(ixt,ipha)
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO k = 1, llm
        DO i = ijb, ije
          IF(ABS(q(i,k,iq))<borne) CYCLE
          WRITE(msg1,'(s,"(",i0,",",i0,",",i0,") = ",ES12.4)')isoName(ixt),i,k,iq,q(i,k,iq); CALL msg(msg1)
          CALL abort_gcm(modname, 'Error in iso_verif_noNaN: '//TRIM(err_msg), 1)
          STOP
        END DO
      END DO
!$OMP END DO NOWAIT
    END DO
  END DO

  !--- CHECK CONSERVATION (MAIN ISOTOPE AND PARENT CONCENTRATIONS MUST BE EQUAL)
  ixt = iH2O
  IF(ixt /= 0) THEN
    DO ipha = 1, npha
      iq = iTraPha(ixt,ipha)
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO k = 1, llm
        DO i = ijb, ije
          q1 = q(i,k,ipha); q2 = q(i,k,iq)
          IF(ABS(q1-q2) <= errmax .OR. ABS(q1-q2)/MAX(MAX(ABS(q1),ABS(q2)),1e-18) < errmaxrel) CYCLE
          WRITE(msg1,'("ixt = ",i0)')ixt;                                      CALL msg(msg1)
          WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q2;  CALL msg(msg1)
          WRITE(msg1,'("q(",i0,",",i0,",ipha=",i0,") = ",ES12.4)')i,k,ipha,q1; CALL msg(msg1)
          CALL abort_gcm(modname, 'Error in iso_verif_egalite: '//TRIM(err_msg), 1)
          q(i,k,iq) = q(i,k,ipha)                !--- Bidouille pour convergence
        END DO
      END DO
!$OMP END DO NOWAIT
    END DO
  END IF

  !--- CHECK DELTA ANOMALIES
  ix = [ixHDO, ixO18]
  DO iiso = 1, SIZE(ix)
    ixt = ix(iiso)
    IF(ixt  == 0) CYCLE
    DO ipha = 1, npha
      iq = iTraPha(ixt,ipha)
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO k = 1, llm
        DO i = ijb, ije
          q1 = q(i,k,ipha); q2 = q(i,k,iq)
          IF(q2 <= qmin) CYCLE
          deltaD = (q2/q1/tnat(ixt)-1)*1000
          IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
          WRITE(msg1,'("ixt = ",i0)')ixt;                                     CALL msg(msg1)
          WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q2; CALL msg(msg1)
          WRITE(msg1,'("q=",ES12.4)')q(i,k,:);                                CALL msg(msg1)
          WRITE(msg1,'("deltaD=",ES12.4)')deltaD;                             CALL msg(msg1)
          CALL abort_gcm(modname, 'Error in iso_verif_aberrant: '//TRIM(err_msg), 1)
        END DO
      END DO
!$OMP END DO NOWAIT
    END DO
  END DO

  !--- CHECK FOR TAGGING TRACERS DELTAD ANOMALIES
  IF(nitr == 0) RETURN
  IF(ixH2O /= 0 .AND. ixHDO /= 0) THEN
    DO izon = 1, nzon
      ixt  = iZonIso(izon, ixHDO)
      ieau = iZonIso(izon, ixH2O)
      DO ipha = 1, npha
        iq    = iTraPha(ixt,  ipha)
        iqeau = iTraPha(ieau, ipha)
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
        DO k = 1, llm
          DO i = ijb, ije
            IF(q(i,k,iq)<=qmin) CYCLE
            deltaD = (q(i,k,iq)/q(i,k,iqeau)/tnat(ixHDO)-1)*1000
            IF(deltaD <= deltaDmax .AND. deltaD >= deltaDmin) CYCLE
            WRITE(msg1,'("izon, ipha =",2i0)')izon, ipha;                              CALL msg(msg1)
            WRITE(msg1,'( "ixt, ieau =",2i0)') ixt, ieau;                              CALL msg(msg1)
            WRITE(msg1,'("q(",i0,",",i0,",iq=",i0,") = ",ES12.4)')i, k, iq, q(i,k,iq); CALL msg(msg1)
            WRITE(msg1,'("deltaD=",ES12.4)')deltaD;                                    CALL msg(msg1)
            CALL abort_gcm(modname, 'Error in iso_verif_aberrant trac: '//TRIM(err_msg), 1)
          END DO
        END DO
!$OMP END DO NOWAIT
      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, npha
      iq = iTraPha(iiso, ipha)
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
      DO k = 1, llm
        DO i = ijb, ije
          xiiso = q(i,k,iq)
          xtractot = SUM(q(i, k, iTraPha(iZonIso(1:nzon,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))
            WRITE(msg1,'("iiso, ipha =",2i0)')iiso, ipha;              CALL msg(msg1)
            WRITE(msg1,'("i, k =",2i0)')i, k;                          CALL msg(msg1)
            WRITE(msg1,'("q(",i0,",",i0,":) = ",ES12.4)')i,k,q(i,k,:); CALL msg(msg1)
            STOP
          END IF
          IF(ABS(xtractot) <= ridicule) CYCLE
          DO izon = 1, nzon
            ixt = iZonIso(izon, iiso) 
            q(i,k,iq) = q(i,k,iq) / xtractot * xiiso
          END DO
        END DO
      END DO
!$OMP END DO NOWAIT
    END DO
  END DO

END SUBROUTINE check_isotopes