SUBROUTINE check_isotopes(q, ijb, ije, err_msg)
   USE parallel_lmdz
   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
   USE iso_params_mod, ONLY: tnat_H216O, tnat_H217O, tnat_H218O, tnat_HDO, tnat_HTO


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

   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=maxlen) :: modname, msg1, nm(2)
   INTEGER :: ixt, ipha, k, i, iq, iiso, izon, ieau, iqeau, iqpar
   INTEGER, ALLOCATABLE       ::   ix(:)
   REAL,    ALLOCATABLE, SAVE :: tnat(:)         !--- OpenMP shared variable
   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
!$OMP THREADPRIVATE(iso_eau, iso_HDO, iso_O18, iso_O17, iso_HTO)
   LOGICAL       :: ltnat1
   LOGICAL, SAVE :: first=.TRUE.
!$OMP THREADPRIVATE(first)

   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
!$OMP MASTER
      ltnat1 = .TRUE.; CALL getin('tnateq1', ltnat1)
      ALLOCATE(tnat(niso))
      iso_eau = strIdx(isoName,'H216O'); IF(iso_eau /= 0) tnat(iso_eau) = tnat_H216O
      iso_O17 = strIdx(isoName,'H217O'); IF(iso_O17 /= 0) tnat(iso_O17) = tnat_H217O
      iso_O18 = strIdx(isoName,'H218O'); IF(iso_O18 /= 0) tnat(iso_O18) = tnat_H218O
      iso_HDO = strIdx(isoName,'HDO');   IF(iso_HDO /= 0) tnat(iso_HDO) = tnat_HDO
      iso_HTO = strIdx(isoName,'HTO');   IF(iso_HTO /= 0) tnat(iso_HTO) = tnat_HTO
      IF(ltnat1) tnat(:) = 1.0
!$OMP END MASTER
!$OMP BARRIER
      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)
!$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)')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
!$OMP END DO NOWAIT
      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
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
         DO k = 1, llm
            DO i = ijb, ije
               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
!$OMP END DO NOWAIT
      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
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
         DO k = 1, llm
            DO i = ijb, ije
               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
!$OMP END DO NOWAIT
      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)
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
            DO k = 1, llm
               DO i = ijb, ije
                  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
!$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, nphas
         iq = iqIsoPha(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, 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
!$OMP END DO NOWAIT
      END DO
   END DO

END SUBROUTINE check_isotopes