
#ifdef ISOVERIF
! $Id: $

MODULE isotopes_verif_mod
  USE isotopes_mod, ONLY: iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO, small => ridicule, tnat
  USE infotrac_phy, ONLY: isoName, isoPhas, isoZone, iqIsoPha, tracers, nqtot, isoSelect, ntraciso => ntiso, &
            nbIso,  niso,   ntiso,   nphas,   nzone, itZonIso, isotope, ixIso, isoFamilies, delPhase
  USE  strings_mod, ONLY: cat, dispOutliers, lunout, num2str, msg, strStack, strIdx, maxlen
  USE phys_local_var_mod, ONLY: qx => qx_seri
  USE ioipsl_getin_p_mod, ONLY : getin_p

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: delta2r, delta,  checkNaN, checkEqu, checkMass !, deltaR
  PUBLIC :: r2delta, excess, checkPos, checkBnd, checkTrac, abs_err, rel_err

  PUBLIC :: errmax, errmaxrel, deltalim, faccond, o17_verif, tmin_verif, faible_evap, deltaDfaible, deltalim_snow, tmax_verif, nlevmaxo17, &
            errmax_sol

  PUBLIC :: iso_verif_init

  PUBLIC :: iso_verif_aberrant,            iso_verif_aberrant_nostop
  PUBLIC :: iso_verif_aberrant_choix,      iso_verif_aberrant_choix_nostop
  PUBLIC :: iso_verif_aberrant_vect2D
  PUBLIC :: iso_verif_aberrant_vect2Dch
  PUBLIC :: iso_verif_aberrant_encadre,    iso_verif_aberrant_enc_nostop
  PUBLIC :: iso_verif_aberrant_enc_vect2D, iso_verif_aberrant_enc_vect2D_ns
  PUBLIC :: iso_verif_aberrant_enc_par2D
  PUBLIC :: iso_verif_aberrant_enc_choix_nostop

  PUBLIC :: iso_verif_aberrant_o17, iso_verif_aberrant_o17_nostop
  PUBLIC :: iso_verif_O18_aberrant, iso_verif_O18_aberrant_nostop
  PUBLIC :: iso_verif_O18_aberrant_enc_vect2D

  PUBLIC :: iso_verif_egalite,        iso_verif_egalite_nostop
  PUBLIC :: iso_verif_egalite_choix,  iso_verif_egalite_choix_nostop
  PUBLIC :: iso_verif_egalite_par2D
  PUBLIC :: iso_verif_egalite_vect1D, iso_verif_egalite_vect2D
  PUBLIC :: iso_verif_egalite_std_vect

  PUBLIC :: iso_verif_noNaN, iso_verif_noNaN_nostop
  PUBLIC :: iso_verif_noNaN_par2D
  PUBLIC :: iso_verif_noNaN_vect1D
  PUBLIC :: iso_verif_noNaN_vect2D

  PUBLIC :: iso_verif_positif,        iso_verif_positif_nostop
  PUBLIC :: iso_verif_positif_vect
  PUBLIC :: iso_verif_positif_choix,  iso_verif_positif_choix_nostop
  PUBLIC :: iso_verif_positif_choix_vect
  PUBLIC :: iso_verif_positif_strict, iso_verif_positif_strict_nostop

#ifdef ISOTRAC
  PUBLIC :: iso_verif_tag17_q_deltaD_vect
  PUBLIC :: iso_verif_tag17_q_deltaD_vect_ret3D
  PUBLIC :: iso_verif_tag17_q_deltaD_chns
  PUBLIC :: iso_verif_trac17_q_deltaD

  PUBLIC :: iso_verif_tracdd_vect
  PUBLIC :: iso_verif_tracdD_choix_nostop

  PUBLIC :: iso_verif_traceur, iso_verif_traceur_nostop
  PUBLIC :: iso_verif_traceur_vect
  PUBLIC :: iso_verif_traceur_choix, iso_verif_traceur_choix_nostop

  PUBLIC :: iso_verif_tracnps
  PUBLIC :: iso_verif_tracnps_vect
  PUBLIC :: iso_verif_tracnps_choix_nostop
  PUBLIC :: iso_verif_tracpos_vect
  PUBLIC :: iso_verif_tracpos_choix, iso_verif_tracpos_choix_nostop

  PUBLIC :: iso_verif_trac_masse_vect
  PUBLIC :: iso_verif_traceur_justmass, iso_verif_traceur_jm_nostop
  PUBLIC :: iso_verif_traceur_noNaN_nostop
  PUBLIC :: iso_verif_traceur_noNaN_vect
  PUBLIC :: iso_verif_tracm_choix_nostop
#endif

  PUBLIC :: deltaD, deltaO, dexcess, delta_all, delta_to_R, o17excess


! SUPPRIMEE: deltaDfaible_lax=-180.0


!===============================================================================================================================
  INTERFACE delta2r;   MODULE PROCEDURE    delta2r_0D,    delta2r_1D,    delta2r_2D;                 END INTERFACE delta2r
  INTERFACE r2delta;   MODULE PROCEDURE    r2delta_0D,    r2delta_1D,    r2delta_2D;                 END INTERFACE r2delta
  INTERFACE checkNaN;  MODULE PROCEDURE   checkNaN_0D,   checkNaN_1D,   checkNaN_2D,    checkNaN_3D; END INTERFACE checkNaN
  INTERFACE checkPos;  MODULE PROCEDURE   checkPos_0D,   checkPos_1D,   checkPos_2D,    checkPos_3D; END INTERFACE checkPos
  INTERFACE checkBnd;  MODULE PROCEDURE   checkBnd_0D,   checkBnd_1D,   checkBnd_2D,    checkBnd_3D; END INTERFACE checkBnd
  INTERFACE checkEqu;  MODULE PROCEDURE   checkEqu_0D,   checkEqu_1D,   checkEqu_2D,    checkEqu_3D; END INTERFACE checkEqu
  INTERFACE checkMass; MODULE PROCEDURE checkMassQ_0D, checkMassQ_1D, checkMassQ_2D, checkMassQx_2D; END INTERFACE checkMass
  INTERFACE delta;     MODULE PROCEDURE     deltaQ_0D,     deltaQ_1D,     deltaQ_2D,     deltaQx_2D; END INTERFACE delta
  INTERFACE deltaR;    MODULE PROCEDURE     deltaR_0D,     deltaR_1D,     deltaR_2D;                 END INTERFACE deltaR
  INTERFACE excess;    MODULE PROCEDURE    excessQ_0D,    excessQ_1D,    excessQ_2D, &
                                           excessR_0D,    excessR_1D,    excessR_2D,    excessQx_2D; END INTERFACE excess
!===============================================================================================================================
  TYPE :: r; REAL,    ALLOCATABLE :: r(:); END TYPE r
  TYPE :: i; INTEGER, ALLOCATABLE :: i(:); END TYPE i
!===============================================================================================================================
                                                 !    DESCRIPTION                                     PREVIOUS NAME
  REAL, PARAMETER :: abs_err      = 1.e-8,     & !--- Max. absolute error allowed                     errmax
                     rel_err      = 1.e-3,     & !--- Max. relative error allowed                     errmax_rel
                     max_val      = 1.e19,     & !--- Max. value allowed      (NaN     if greater)    borne
                     faccond      = 1.1,       & !--- In liquids, R(max_delD)*faccond is allowed
                     min_T        = 5.0,       & !--- Min. realistic temperature value (in K)         Tmin_verif
                     max_T        = 1000.0,    & !--- Max. realistic temperature value (in K)         Tmax_verif
                     low_delD     = -300.0,    & !--- Max quantity of vapour making outliers from     deltaDfaible
                                                 !    balance with ground acceptable
                     min_delDevap = 3.0,       & !--- Less demanding with ORCHIDEE evap outliers      faible_evap
                     slope_HDO    = 8.0,       & !--- H[2]HO  excess law: slope
                     slope_O17    = 0.528        !--- H2[17]O excess law: slope
  CHARACTER(LEN=3), PARAMETER ::               &
                     progr_HDO    = 'lin',     & !--- H[2]HO  excess law: progression type
                     progr_O17    = 'log'        !--- H2[17]O excess law: progression type
  REAL,      SAVE :: min_delD,                 & !--- Min. deltaD allowed     (defined in iso.def)    deltaDmin
                     max_delD,                 & !--- Max. deltaD for vapour  (outlier if greater)    deltalim
                     max_delDtrac,             & !--- Max. deltaD for tracers (defined in iso.def)
                     max_delDsnow,             & !--- Max. deltaD for snow    (outlier if greater)
                     min_Dexc,                 & !--- Min. value of   D-excess (excess of H2[18]O     dexcess_min
                     max_Dexc,                 & !--- Max. value    "    "     w/resp. to H[2]HO)     dexcess_max
                     min_O17exc,               & !--- Min. value of O17-excess (excess of H2[17]O     O17excess_bas
                     max_O17exc                  !--- Max. value    "    "     w/resp. to H2[18]O)    O17excess_haut

  LOGICAL,   SAVE :: checkO17                    !--- Flag to trigger delta(H2[17]O) checking         O17_verif
  INTEGER,   SAVE :: max_O17nlev                 !---                                                 nlevmaxO17
!$OMP THREADPRIVATE(min_delD, max_delDtrac, min_Dexc, min_O17exc, checkO17)
!$OMP THREADPRIVATE(max_delD, max_delDsnow, max_Dexc, max_O17exc, max_O17nlev)
!      logical bidouille_anti_divergence ! .TRUE. => xt relaxed to q to avoid errors accumulation leading to divergence
                                         ! Moved to wateriso2.h, rzad at execution

  !=== QUANTITIES DEPENDING ON THE ISOTOPES FAMILIES (LENGTH: nbIso) ; SET UP FIRST TIME ONLY IN "isoCheck_init"
  TYPE(r), ALLOCATABLE, SAVE :: dmin(:),dmax(:),&!--- Thresholds for correct delta  for each isotope (len: isotopes families nb)
                                emin(:),emax(:)  !--- Thresholds for correct excess for each isotope (len: isotopes families nb)
  TYPE(i), ALLOCATABLE, SAVE :: iMajr(:)         !--- Indices of the major isotopes (to be compared with parent concentration)
  CHARACTER(LEN=6), ALLOCATABLE, SAVE ::       &
                       sIsoCheck(:)              !--- Activated checking routines keywords list ('npbem' ; 'x' if not activated)

  !=== VALUES TO BE SET UP BEFORE CALLING ELEMENTAL FUNCTIONS
  REAL,    SAVE :: epsRel, epsAbs, epsPos        !--- NEEDED BY "isEq???".               Set before each checkEqual call
  REAL,    SAVE :: min_bnd, max_bnd              !--- NEEDED BY "isBounded".             Set before each checkBound call
  REAL,    SAVE :: tnat0                         !--- NEEDED BY "r2delta" and "delta2r". Set before each delta*     call
  REAL,             SAVE :: slope                !--- NEEDED BY "isoExcess".             Set using "set_isoParams"
  CHARACTER(LEN=3), SAVE :: lawType              !---  "  "
  INTEGER,          SAVE :: iso_ref              !---  "  "

! variables de vérifications
real errmax ! erreur maximale en absolu.
parameter (errmax=abs_err)
real errmax_sol ! erreur maximale en absolu.
parameter (errmax_sol=5e-7)
      real errmaxrel ! erreur maximale en relatif autorisée
!      parameter (errmaxrel=1e10)
      parameter (errmaxrel=rel_err)
      real borne ! valeur maximale que n'importe quelle variable peut
                 ! atteindre, en val abs; utile pour vérif que pas NAN
      parameter (borne=max_val)
      real, save ::  deltalim ! deltalim est le maximum de deltaD qu'on
                             ! autorise dans la vapeur, au-dela, en suppose que c'est abérrant.
                             ! dans le liquide, on autorise deltalim*faccond.
!$OMP THREADPRIVATE(deltalim)
!      parameter (deltalim=1e10)
!      parameter (deltalim=300.0)
       ! maintenant défini dans iso.def

       real, save :: deltalimtrac ! max de deltaD dans les traceurs, défini dans iso.def
!$OMP THREADPRIVATE(deltalimtrac)

      real, save ::  deltalim_snow ! deltalim est le maximum de deltaD qu'on
                             ! autorise dans la neige, au-dela, en suppose que c'est abérrant.
!$OMP THREADPRIVATE(deltalim_snow)
!      parameter (deltalim_snow=500.0) ! initalisé dans iso_init

	real, save ::  deltaDmin
!$OMP THREADPRIVATE(deltaDmin)
!	parameter (deltaDmin=-900.0)
    ! maintentant, défini dans iso.def

      real, save ::  o17excess_bas,o17excess_haut ! bornes inf et sup de l'O17-excess
!      parameter(o17excess_bas=-200.0,o17excess_haut=120)
!$OMP THREADPRIVATE(o17excess_bas,o17excess_haut)
      integer nlevmaxO17
!$OMP THREADPRIVATE(nlevmaxO17)

      logical, save ::  O17_verif
!$OMP THREADPRIVATE(O17_verif)
!      parameter (O17_verif=.true.)

	real, save ::  dexcess_min,dexcess_max
!$OMP THREADPRIVATE(dexcess_min,dexcess_max)

!      real faccond  ! dans le liquide, on autorise R(deltalim)*faccond.
!      parameter (faccond=1.1)
      
!      logical bidouille_anti_divergence ! si true, alors on fait un
!                                        ! rappel régulier des xt4 vers les q pour 
!                                        !éviter accumulation d'érreurs et  divergences
!      parameter (bidouille_anti_divergence=.true.)
!      parameter (bidouille_anti_divergence=.false.)
    ! bidouille_anti_divergence a été déplacé dans wateriso2.h et est lue à l'execution

        
	real deltaDfaible ! deltaD en dessous duquel la vapeur est tellement faible 
        !que on peut accepter que la remise à l'équilibre du sol avec 
        ! cette vapeur donne des deltaDevap aberrants.
        parameter (deltaDfaible=low_delD)

    	real faible_evap ! faible évaporation: on est plus laxiste 
        !pour les deltaD aberrant dans le cas de l'évap venant d'orchidee
        parameter (faible_evap=min_delDevap)

    	real Tmin_verif
        parameter (Tmin_verif=min_T) ! en K
    	real Tmax_verif
        parameter (Tmax_verif=max_T) ! en K


! subroutines de vérifications génériques, à ne pas modifier

 
CONTAINS


LOGICAL FUNCTION checkiIso(iIso, sub, iFamily) RESULT(lerr)
  INTEGER,                    INTENT(IN) :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
  INTEGER,          OPTIONAL, INTENT(IN) :: iFamily
  CHARACTER(LEN=maxlen) :: subn
  INTEGER :: isav
  LOGICAL :: ll
  subn = 'checkiIso'; IF(PRESENT(sub)) subn = TRIM(sub)//TRIM(subn)
  IF(PRESENT(iFamily)) THEN
     iSav = ixIso; lerr = isoSelect(iFamily); IF(lerr) RETURN
     lerr = iIso <= 0 .OR. iIso > isotope%niso
     ll = isoSelect(iSav)
  ELSE
     lerr = iIso <= 0 .OR. iIso > isotope%niso
  END IF
  CALL msg('wrong isotopic index iIso = '//TRIM(num2str(iIso))//' (must be >0, <= '//TRIM(num2str(iIso))//')', subn, lerr)
END FUNCTION checkiIso


LOGICAL FUNCTION isoCheck_init() RESULT(lerr)
  CHARACTER(LEN=maxlen) :: isoFamily, modname='isoCheck_init'
  INTEGER :: iIso
  LOGICAL :: isoCheck

  CALL msg('Starting...', modname)
  ALLOCATE(dmin(nbIso), dmax(nbIso), emin(nbIso), emax(nbIso), sIsoCheck(nbIso), iMajr(nbIso))
  CALL getin_p('isoCheck', isoCheck, .FALSE.)
  sIsoCheck = 'xxxxxx'; IF(isoCheck) sIsoCheck = 'npmdeb'  !=== Single keyword for all isotopes families -> CAN BE GENERALIZED !

  DO iIso = 1, nbIso
     isoFamily = isotope%parent
     SELECT CASE(isoFamily)
        !-----------------------------------------------------------------------------------------------------------------------
        CASE('H2O')
        !-----------------------------------------------------------------------------------------------------------------------
           iMajr(iIso)%i = [iso_eau]
           WRITE(*,*)'iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO = ', iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO
           iso_ref     = iso_eau
           CALL getin_p('min_delD',     min_delD,    -900.0)
           CALL getin_p('max_delD',     max_delD,     300.0)
           CALL getin_p('max_delDtrac', max_delDtrac, 300.0)
           max_delDsnow =  200.0 + max_delD
           WRITE(*,*)'max_delDsnow = ', max_delDsnow
           IF(iso_O17 > 0 .AND. iso_O18 > 0) THEN
              CALL getin_p('checkO17',    checkO17, .FALSE.)
              CALL getin_p('min_O17exc',  min_O17exc, -200.0)
              CALL getin_p('max_O17exc',  max_O17exc,  120.0)
              CALL getin_p('max_O17nlev', max_O17nlev, 1000)
           END IF
           IF(iso_HDO > 0 .AND. iso_O18 > 0) THEN
              CALL getin_p('min_Dexc', min_Dexc, -100.0)
              CALL getin_p('max_Dexc', max_Dexc, 1000.0)
           END IF

           !=== HANDY VECTORS OF ACCEPTABLE LIMITS FOR delta
           ALLOCATE(dmin(iIso)%r(niso)); dmin(iIso)%r(:)=-max_val
           ALLOCATE(dmax(iIso)%r(niso)); dmax(iIso)%r(:)=+max_val
           IF(iso_HDO /= 0) THEN; dmin(iIso)%r(iso_HDO) = min_delD;     dmax(iIso)%r(iso_HDO) = max_delD;     END IF
           IF(iso_O17 /= 0) THEN; dmin(iIso)%r(iso_O17) =-max_val;      dmax(iIso)%r(iso_O17) = max_val;      END IF
           IF(iso_O18 /= 0) THEN; dmin(iIso)%r(iso_O18) = min_delD/2.0; dmax(iIso)%r(iso_O18) = max_delD/8.0; END IF

           !=== HANDY VECTORS OF ACCEPTABLE LIMITS FOR excess
           ALLOCATE(emin(iIso)%r(niso)); emin(iIso)%r(:)=-max_val
           ALLOCATE(emax(iIso)%r(niso)); emax(iIso)%r(:)=+max_val
           IF(iso_O17 /= 0) THEN; emin(iIso)%r(iso_O17) = min_O17exc;   emax(iIso)%r(iso_O17) = max_O17exc;   END IF
           IF(iso_O18 /= 0) THEN; emin(iIso)%r(iso_O18) = min_Dexc;     emax(iIso)%r(iso_O18) = max_Dexc;     END IF
        !-----------------------------------------------------------------------------------------------------------------------
        CASE DEFAULT
           lerr = .TRUE.
           CALL msg('routine not yet set up for isotopes family "'//TRIM(isoFamily)//'"', modname, lerr)
        !-----------------------------------------------------------------------------------------------------------------------
     END SELECT
  END DO

END FUNCTION isoCheck_init


!===============================================================================================================================
!=== ANCILLARY ROUTINES
!===============================================================================================================================
!=== CHECK ISOTOPIC INDEX "iIso" AND DETERMINE FEW USEFULL QUANTITIES, FOR ELEMENTAL FUNCTIONS AND FOR ERROR MESSAGES:
!===      lerr = set_isoParams(iIso, err_tag, sub[, mss][, phas][, iqIso][, iqPar])
!===
!=== ARGUMENTS:                                                                         Intent  Optional?  Default value
!===    * iIso     Index to be checked                                                    IN    MANDATORY
!===    * err_tag  Tag for error messages, appended at the beginning of "mss"             IN    OPTIONAL
!===    * sub      Calling sequence, appended with current routine name (":" separator)   IN    OPTIONAL
!===    * mss      Message in case absurd values are found                                OUT   OPTIONAL
!===    * phas     Phase number, needed for the computation of "iqIso" and "iqPar"        IN    OPTIONAL
!===    * iqIso    Index in "qx" of the "iIso"th isotope                                  OUT   OPTIONAL
!===    * iqPar    Index in "qx" of the "iIso"th isotope parent                           OUT   OPTIONAL
!===
!=== NOTES:          * SET CORRECT "slope", "iRef" AND "iso_ref" FOR ISOTOPIC EXCESS COMPUTATION.
!===                      SO FAR: H[2]HO-excess, H2[17]O-excess (with respect to H2[18]O)
!===                 * THE FOLLOWING DELTAs NEED TO BE COMPUTED:  delta-H[2]HO,  delta-H2[17]O,  delta-H2[18]O
!===============================================================================================================================
LOGICAL FUNCTION set_isoParams(iIso, err_tag, sub, mss, phas, iqIso, iqPar) RESULT(lerr)
  INTEGER,                         INTENT(IN)  :: iIso
  CHARACTER(LEN=*),                INTENT(IN)  :: err_tag
  CHARACTER(LEN=maxlen), OPTIONAL, INTENT(OUT) :: sub, mss
  CHARACTER(LEN=*),      OPTIONAL, INTENT(IN)  :: phas
  INTEGER,               OPTIONAL, INTENT(OUT) :: iqIso, iqPar
  INTEGER :: iq, iPha, ix
  LOGICAL :: lExc
  CHARACTER(LEN=maxlen) :: iname, subn
  subn = 'set_isoParams'; IF(PRESENT(sub)) subn = TRIM(sub)//TRIM(subn)
  lerr = checkiIso(iIso, subn); IF(lerr) RETURN
  IF(PRESENT(phas)) THEN
     iPha = INDEX(isoPhas, phas); lerr = iPha == 0
     CALL msg('isotopes family "'//TRIM(isotope%parent)//'" has no "'//TRIM(phas)//'" phase', subn, lerr)
  ELSE
     iPha = 0; lerr = PRESENT(iqIso) .OR. PRESENT(iqPar)
     CALL msg('missing phase, needed to compute iqIso or iqPar', subn, lerr)
  END IF
  IF(lerr) RETURN
  ix = INDEX(sub, ':')
  SELECT CASE(sub(ix+1:LEN_TRIM(sub)))
     CASE('delta' ); iname = vname(iIso, 'd', iPha)
        IF(iPha /= 0)   iq = iqIsoPha(iIso, iPha)
        IF(PRESENT(iqIso)) iqIso = iq
        IF(PRESENT(iqPar)) iqPar = tracers(iq)%iqParent
     CASE('excess'); iname = vname(iIso, 'e', iPha)
        IF     (iIso == iso_HDO) THEN
           slope = slope_HDO; lawType = progr_HDO; iso_ref = iso_O18
        ELSE IF(iIso == iso_O17) THEN
          slope = slope_O17; lawType = progr_O17; iso_ref = iso_O18
        ELSE
           CALL msg('missing parameters for "'//TRIM(iname)//'" (iIso = '//TRIM(num2str(iIso))//'")', subn)
           lerr = .TRUE.; RETURN
        END IF
  END SELECT
  IF(PRESENT(mss)) THEN
     mss = 'absurd '//TRIM(iname)//' values found'
     IF(err_tag /= '') mss = TRIM(err_tag)//':'//TRIM(mss)
  END IF
END FUNCTION set_isoParams
!===============================================================================================================================
!=== DETERMINE THE NAME OF THE QUANTITY COMPUTED DEPENDING ON "typ": "d"elta, "e"xcess OR "r"atio.
!===============================================================================================================================
CHARACTER(LEN=maxlen) FUNCTION vName(iIso, typ, iPha, iZon) RESULT(nam)
  INTEGER,           INTENT(IN) :: iIso
  CHARACTER(LEN=*),  INTENT(IN) :: typ
  INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon
  INTEGER :: ix, ip
  IF(iIso == 0) THEN
     IF(typ == 'd') nam = 'delta'
     IF(typ == 'e') nam = 'excess'
     IF(typ == 'r') nam = 'ratio'
  ELSE
     ip = 0;    IF(PRESENT(iPha)) ip = iPha
     ix = iIso; IF(PRESENT(iZon)) ix = itZonIso(iZon, iIso)
     IF(ip  /=  0 ) nam = tracers(iqIsoPha(ix, iPha))%name
     IF(ip  ==  0 ) nam = isoName(ix)
     IF(typ == 'd') nam = 'delta-'//TRIM(nam)
     IF(typ == 'e') nam = TRIM(nam)//'-excess'
     IF(typ == 'r') nam = TRIM(nam)//'-ratio'
  END IF
END FUNCTION vName
!===============================================================================================================================
!=== SAME AS vName, FOR SEVERAL INDICES "iIso"
!===============================================================================================================================
FUNCTION vNames(iIso, typs, iPha, iZon) RESULT(nam)
  INTEGER,           INTENT(IN) :: iIso(:)
  CHARACTER(LEN=*),  INTENT(IN) :: typs
  INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon
  CHARACTER(LEN=maxlen)         :: nam(SIZE(iIso))
  INTEGER :: it
  DO it = 1, SIZE(nam)
     nam(it) = vName(iIso(it), typs(it:it), iPha, iZon)
  END DO
END FUNCTION vNames
!===============================================================================================================================
!=== LIST OF TWO VARIABLES FOR DELTA CASE: ratio(iIso), delta(iIso)
!===============================================================================================================================
FUNCTION vNamDe(iIso, iPha, iZon) RESULT(nam)
  INTEGER,           INTENT(IN) :: iIso
  INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon
  CHARACTER(LEN=maxlen) :: nam(2)
  nam = vNames([iIso, iso_ref], 'rd', iPha, iZon)
END FUNCTION vNamDe
!===============================================================================================================================
!=== LIST OF THREE VARIABLES FOR EXCESS CASE: delta(iIso), delta(iParent), excess(iIso)
!===============================================================================================================================
FUNCTION vNamEx(iIso, iPha, iZon) RESULT(nam)
  INTEGER,           INTENT(IN) :: iIso
  INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon
  CHARACTER(LEN=maxlen) :: nam(3)
  nam = vNames([iIso, iso_ref, iIso], 'dde', iPha, iZon)
END FUNCTION vNamEx
!===============================================================================================================================
!=== GET THE VARIABLES NAMES LIST (DEFAULT VALUE OR FROM tracers(:)%name, isoName(:) OR OPTIONAL ARGUMENT nam(:))
!===============================================================================================================================
LOGICAL FUNCTION gettNam(nq, sub, tnam, nam) RESULT(lerr)
  INTEGER,                       INTENT(IN)  ::   nq(:)
  CHARACTER(LEN=*),              INTENT(IN)  ::  sub
  CHARACTER(LEN=*), ALLOCATABLE, INTENT(OUT) :: tnam(:)
  CHARACTER(LEN=*), OPTIONAL,    INTENT(IN)  ::  nam(:)
  INTEGER :: lq, ln, rk
  rk = SIZE(nq); lq = SIZE(nq, DIM=rk)
  IF(PRESENT(nam)) THEN; ln = SIZE(nam)
     lerr = ALL([1, lq] /= ln)
     CALL msg('SIZE(q,'//TRIM(num2str(rk))//')='//TRIM(num2str(lq))//' /= SIZE(nam)='//num2str(ln), sub, lerr); IF(lerr) RETURN
     tnam = nam
  ELSE
     tnam = ['q']
     IF(lq == ntiso) tnam = isoName(:)
     IF(lq == nqtot) tnam = tracers(:)%name
  END IF
END FUNCTION gettNam
!===============================================================================================================================




!===============================================================================================================================
!===        BASIC ELEMENTAL FUNCTIONS
!===   FUNCTION isNaN  (a)                         Not a Number (NaNs) detection: TRUE if |a|   > max_val
!===   FUNCTION isNeg  (a)                         Positiveness detection:        TRUE IF  a    < epsPos
!===   FUNCTION isEqAbs(a,b)                (1)    Absolute equality checking:    TRUE if |a-b| < epsAbs
!===   FUNCTION isEqRel(a,b)                (1)    Relative equality checking:    TRUE if |a-b|/max(|a|,|b|,1e-18) < epsRel
!===   FUNCTION isEqual(a,b)                (1)    Equality checking (both criteria)
!===   FUNCTION delta2ratio(ratioIso)       (2)    Delta  function from a   ratio value
!===   FUNCTION ratio2delta(ratioIso)       (2)    Ratio  function from a   delta value
!===   FUNCTION isoExcess(delIso, delRef)   (3)    Excess function from two delta values
!===
!===   NOTES:     (1): epsPos and epsAbs must be defined before calling
!===              (2): tnat0             must be defined before calling
!===              (3): slope and lawType must be defined before calling
!=== 
!===============================================================================================================================
ELEMENTAL LOGICAL FUNCTION isNaN(a) RESULT(lout)                     !--- IS "a" AN OUTLIER,   ie: a < -max_val or a > max_val ?
  REAL, INTENT(IN) :: a
  lout = a < -max_val .OR. max_val < a
END FUNCTION isNaN
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL LOGICAL FUNCTION isBounded(a) RESULT(lout)                 !--- IS "a" BOUNDED,      ie:    min_bnd < a < max_bnd ?
  REAL, INTENT(IN) :: a
  lout= min_bnd < a .AND. a < max_bnd
END FUNCTION isBounded
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL LOGICAL FUNCTION isNeg(a) RESULT(lout)                     !--- IS "a" NEGATIVE OR NOT, ie: a ≤ epsPos ?
  REAL, INTENT(IN) :: a
  lout = a <= epsPos
END FUNCTION isNeg
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL LOGICAL FUNCTION isPos(a) RESULT(lout)                     !--- IS "a" POSITIVE OR NOT, ie: a ≥ epsPos ?
  REAL, INTENT(IN) :: a
  lout = a >= epsPos
END FUNCTION isPos
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL LOGICAL FUNCTION isEqAbs(a,b) RESULT(lout)                 !--- IS "a" ABSOLUTELY EQ TO "b", ie: |a-b|<eps ?
  REAL, INTENT(IN) :: a, b
  lout = abs(a-b) < epsAbs
END FUNCTION isEqAbs
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL LOGICAL FUNCTION isEqRel(a,b) RESULT(lout)                 !--- IS "a" RELATIVELY EQ TO "b", ie: |a-b|/max(|a|,|b|,1e-18)<eps ?
  REAL, INTENT(IN) :: a, b
  lout = ABS(a-b)/MAX(ABS(a),ABS(b),1e-18) < epsRel
END FUNCTION isEqRel
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL LOGICAL FUNCTION isEqual(a,b) RESULT(lout)                 !--- IS "a" EQUAL TO "b", DUAL REL. AND ABS. CRITERIA
  REAL, INTENT(IN) :: a, b
  lout = isEqAbs(a,b); IF(lout) lout=isEqRel(a,b)
END FUNCTION isEqual
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL REAL FUNCTION ratio2delta(ratioIso) RESULT(delta)          !--- CONVERT AN ISOTOPIC RATIO INTO A DELTA
  REAL, INTENT(IN) :: ratioIso
  delta = (ratioIso / tnat0 - 1.0) * 1000.0                          !  /!\ tnat0 must be defined !!!
END FUNCTION ratio2delta
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL REAL FUNCTION delta2ratio(delta) RESULT(ratioIso)          !--- CONVERT A DELTA INTO AN ISOTOPIC RATIO
  REAL, INTENT(IN) :: delta
  ratioIso = (delta / 1000.0 + 1.0) * tnat0                          !  /!\ tnat0 must be defined !!!
END FUNCTION delta2ratio
!-------------------------------------------------------------------------------------------------------------------------------
ELEMENTAL REAL FUNCTION isoExcess(deltaIso, deltaRef) RESULT(excess) !--- COMPUTE EXCESS USING TWO ISOTOPIC RATIOS
  REAL, INTENT(IN) :: deltaIso, deltaRef                             !  /!\ slope and lawType must be defined !!!
  SELECT CASE(lawType)
    CASE('lin'); excess =              deltaIso        - slope *       deltaRef
    CASE('log'); excess = 1.e6*(LOG(1.+deltaIso/1000.) - slope *LOG(1.+deltaRef/1000.))
  END SELECT
END FUNCTION isoExcess
!-------------------------------------------------------------------------------------------------------------------------------


!===============================================================================================================================
!===   R = delta2r(iIso, delta[, lerr]):  convert an anomaly "delta" into a ratio "R"
!=== ARGUMENTS:                                                                     Intent  Optional?
!===   * iIso     Isotope index in "isoName"                                        INPUT   MANDATORY
!===   * delta    delta ratio for isotope of index "iIso"                           INPUT   MANDATORY
!===   * lerr     Error code                                                        OUTPUT  OPTIONAL
!===   * R        Isotopic ratio (output value, rank 1 to 3)                        OUTPUT  MANDATORY
!===============================================================================================================================
REAL FUNCTION delta2r_0D(iIso, delta, lerr) RESULT(R)
  INTEGER,           INTENT(IN)  :: iIso
  REAL,              INTENT(IN)  :: delta
  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
  LOGICAL :: ler
  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
  tnat0 = tnat(iIso); R = delta2ratio(delta)
END FUNCTION delta2r_0D
!-------------------------------------------------------------------------------------------------------------------------------
FUNCTION delta2r_1D(iIso, delta, lerr) RESULT(R)
  REAL,              ALLOCATABLE :: R(:)
  INTEGER,           INTENT(IN)  :: iIso
  REAL,              INTENT(IN)  :: delta(:)
  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
  LOGICAL :: ler
  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
  tnat0 = tnat(iIso); R = delta2ratio(delta)
END FUNCTION delta2r_1D
!-------------------------------------------------------------------------------------------------------------------------------
FUNCTION delta2r_2D(iIso, delta, lerr) RESULT(R)
  REAL,              ALLOCATABLE :: R(:,:)
  INTEGER,           INTENT(IN)  :: iIso
  REAL,              INTENT(IN)  :: delta(:,:)
  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
  LOGICAL :: ler
  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
  tnat0 = tnat(iIso); R = delta2ratio(delta)
END FUNCTION delta2r_2D
!===============================================================================================================================
!===   delta = r2delta(iIso, R[, lerr]):  convert a ratio "R" into an anomaly "delta"
!=== ARGUMENTS:                                                                     Intent  Optional?
!===   * iIso      Isotope index in "isoName"                                       INPUT   MANDATORY
!===   * R         Isotopic ratio for isotope of index "iIso"                       INPUT   MANDATORY
!===   * lerr      Error code                                                       OUTPUT  OPTIONAL
!===   * delta     delta ratio (output value, rank 1 to 3)                          OUTPUT  MANDATORY
!===============================================================================================================================
REAL FUNCTION r2delta_0D(iIso, R, lerr) RESULT(delta)
  INTEGER,           INTENT(IN)  :: iIso
  REAL,              INTENT(IN)  :: R
  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
  LOGICAL :: ler
  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
  tnat0 = tnat(iIso); delta = ratio2delta(R)
END FUNCTION r2delta_0D
!-------------------------------------------------------------------------------------------------------------------------------
FUNCTION r2delta_1D(iIso, R, lerr) RESULT(delta)
  REAL,              ALLOCATABLE :: delta(:)
  INTEGER,           INTENT(IN)  :: iIso
  REAL,              INTENT(IN)  :: R(:)
  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
  LOGICAL :: ler
  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
  tnat0 = tnat(iIso); delta = ratio2delta(R)
END FUNCTION r2delta_1D
!-------------------------------------------------------------------------------------------------------------------------------
FUNCTION r2delta_2D(iIso, R, lerr) RESULT(delta)
  REAL,              ALLOCATABLE :: delta(:,:)
  INTEGER,           INTENT(IN)  :: iIso
  REAL,              INTENT(IN)  :: R(:,:)
  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
  LOGICAL :: ler
  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
  tnat0 = tnat(iIso); delta = ratio2delta(R)
END FUNCTION r2delta_2D
!===============================================================================================================================



!===============================================================================================================================
!===============================================================================================================================
!=== CHECK WHETHER INPUT FIELD "q" CONTAINS Not A Number (NaN) VALUES OR NOT
!===
!===   lerr = checkNaN(q[, err_tag][, subn][, nam][, nmax])
!===
!=== ARGUMENTS:                                                                     intent  optional?  default value
!===   * q         Concentration, scalar or array of rank 1 to 3    (checkNaN)      INPUT   MANDATORY
!===   * err_tag   Error message to display if NaNs are found                       INPUT   OPTIONAL   "NaNs found"
!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkNaN"
!===   * nam       Name(s) of the tracers (last index)                              INPUT   OPTIONAL   "q"
!===   * nmax      Maximum number of outliers values to compute for each tracer     INPUT   OPTIONAL   All values
!===============================================================================================================================
LOGICAL FUNCTION checkNaN_0D(q, err_tag, nam, subn) RESULT(lerr)
  REAL,                       INTENT(IN) :: q
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub ='checkNaN';   IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  lerr = gettNam([1], sub, tnm, [nam]); IF(lerr) RETURN
  lerr = dispOutliers( [isNaN(q)], [q], [1], tag, tnm, sub)
END FUNCTION checkNaN_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkNaN_1D(q, err_tag, nam, subn, nmax) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub ='checkNaN';   IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkNaN_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkNaN_2D(q, err_tag, nam, subn, nmax) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub ='checkNaN';   IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkNaN_2D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkNaN_3D(q, err_tag, nam, subn, nmax) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:,:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub ='checkNaN';   IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkNaN_3D
!===============================================================================================================================


!===============================================================================================================================
!===============================================================================================================================
!=== CHECK WHETHER INPUT FIELD "q" HAVE ONLY POSITIVE VALUES OR NOT
!===
!=== lerr = checkPos(q[, err_tag][, nam][, subn][, nmax][, tny])                    Check for negative values in "q" .
!===
!=== ARGUMENTS:                                                                     Intent  Optional?  Default value
!===   * q         Concentration, scalar or array of rank 1 to 3    (checkNaN)      INPUT   MANDATORY
!===   * err_tag   Error message to display if NaNs are found                       INPUT   OPTIONAL   "Negative values found"
!===   * nam       Name(s) of the tracers (last index)                              INPUT   OPTIONAL   "q"
!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkPositive"
!===   * nmax      Maximum number of outliers values to compute for each tracer     INPUT   OPTIONAL   All values
!===   * tny       Value (-) under which a "q" is considered to be negative         INPUT   OPTIONAL   small
!===============================================================================================================================
LOGICAL FUNCTION checkPos_0D(q, err_tag, nam, subn, tny) RESULT(lerr)
  REAL,                       INTENT(IN) :: q
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn
  REAL,             OPTIONAL, INTENT(IN) :: tny
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub = 'checkPos';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  epsPos = -small;              IF(PRESENT(tny )) epsPos = -tny
  lerr = gettNam(SHAPE(q), sub, tnm, [nam]); IF(lerr) RETURN
  lerr = dispOutliers( [isNeg(q)], [q], [1], tag, tnm, sub)
END FUNCTION checkPos_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkPos_1D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: tny
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub = 'checkPos';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  epsPos = -small;              IF(PRESENT(tny )) epsPos = -tny
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkPos_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkPos_2D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: tny
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub = 'checkPos';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  epsPos = -small;              IF(PRESENT(tny )) epsPos = -tny
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkPos_2D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkPos_3D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:,:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: tny
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub = 'checkPos';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  epsPos = -small;              IF(PRESENT(tny )) epsPos = -tny
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkPos_3D
!===============================================================================================================================


!===============================================================================================================================
!===============================================================================================================================
!=== CHECK WHETHER INPUT FIELD "q" HAVE VALUES INSIDE GIVEN BOUNDS OR NOT (IE "REASONNABLE" VALUES)
!===
!=== lerr = checkBnd(q[, err_tag][, nam][, subn][, nmax][, qmin][, qmax])
!===
!=== ARGUMENTS:                                                                Intent  Optional?  Default value
!===   * q       Tested field, scalar or array of rank 1 to 3                  INPUT   MANDATORY
!===   * err_tag Error message to display if out-of-bounds are found           INPUT   OPTIONAL   "Negative values found"
!===   * nam     Name(s) of the tracers (last index)                           INPUT   OPTIONAL   "q"
!===   * subn    Calling subroutine name                                       INPUT   OPTIONAL   "checkBnd"
!===   * nmax    Maximum number of outliers values to compute for each tracer  INPUT   OPTIONAL   All values
!===   * qmin    Lower bound. Correct values must be > lBnd                    INPUT   OPTIONAL    0
!===   * qmax    Upper bound. Correct values must be < uBnd                    INPUT   OPTIONAL   max_val
!===============================================================================================================================
LOGICAL FUNCTION checkBnd_0D(q, err_tag, nam, subn, qmin, qmax) RESULT(lerr)
  REAL,                       INTENT(IN) :: q
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn
  REAL,             OPTIONAL, INTENT(IN) :: qmin, qmax
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub = 'checkBnd';           IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  min_bnd = 0.;               IF(PRESENT(qmin)) min_bnd = qmin
  max_bnd = max_val;          IF(PRESENT(qmax)) max_bnd = qmax
  lerr = gettNam(SHAPE(q), sub, tnm, [nam]); IF(lerr) RETURN
  tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)'
  lerr = dispOutliers( [isBounded(q)], [q], [1], tag, tnm, sub)
END FUNCTION checkBnd_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkBnd_1D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  REAL,             OPTIONAL, INTENT(IN) :: qmin, qmax
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub = 'checkBnd';           IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  min_bnd = 0.;               IF(PRESENT(qmin)) min_bnd = qmin
  max_bnd = max_val;          IF(PRESENT(qmax)) max_bnd = qmax
  tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)'
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkBnd_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkBnd_2D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  REAL,             OPTIONAL, INTENT(IN) :: qmin, qmax
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub = 'checkBnd';           IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  min_bnd = 0.;               IF(PRESENT(qmin)) min_bnd = qmin
  max_bnd = max_val;          IF(PRESENT(qmax)) max_bnd = qmax
  tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)'
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkBnd_2D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkBnd_3D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr)
  REAL,                       INTENT(IN) :: q(:,:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
  REAL,             OPTIONAL, INTENT(IN) :: qmin, qmax
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  sub = 'checkBnd';           IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
  min_bnd = 0.;               IF(PRESENT(qmin)) min_bnd = qmin
  max_bnd = max_val;          IF(PRESENT(qmax)) max_bnd = qmax
  tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)'
  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
  lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
END FUNCTION checkBnd_3D
!===============================================================================================================================


!===============================================================================================================================
!=== CHECK WHETHER FIELDS "a" AND "b" ARE EQUAL OR NOT (ABSOLUTELY, THE RELATIVELY)
!===
!=== lerr = checkEqu(a, b[, err_tag][, nam][, subn][, nmax][, absErr][, relErr])
!===
!=== ARGUMENTS:                                                                     Intent  Optional?  Default value
!===   * a, b      Fields,  rank 0 to 3                                             INPUT   MANDATORY
!===   * err_tag   Error message to display if fields are not matching              INPUT   OPTIONAL   "mismatching values"
!===   * nam       Name(s) of the tracers (last index)                              INPUT   OPTIONAL   "a","b" or from isoNames
!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkEqu"
!===   * nmax      Maximum number of outliers values to compute for each tracer     INPUT   OPTIONAL   All values
!===   * absErr    Maximum value for |q-r| to consider "q" and "q" as equal         INPUT   OPTIONAL    abs_err
!===   * relErr    Maximum value for |q-r|/max(|q|,|r|,1e-18)      " "              INPUT   OPTIONAL    rel_err
!===============================================================================================================================
LOGICAL FUNCTION checkEqu_0D(a, b, err_tag, nam, subn, absErr, relErr) RESULT(lerr)
  REAL,                       INTENT(IN) :: a, b
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn
  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
  CHARACTER(LEN=maxlen)                  :: tag, sub
  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
  REAL   :: aErr, rErr
  tag = 'mismatching value:'; IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  tnm = ['a  ','b  ','a-b'];  IF(PRESENT(nam    )) tnm = nam                   !--- Variables names
  sub = 'checkEqu';           IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)
  epsAbs = abs_err;           IF(PRESENT(absErr )) epsAbs = absErr             !--- Absolute error
  epsRel = rel_err;           IF(PRESENT(relErr )) epsRel = relErr             !--- Relative error
  lerr = dispOutliers([isEqual(a,b)], cat(a, b, a-b), [1], tag,tnm,sub)
END FUNCTION checkEqu_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkEqu_1D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr)
  REAL,                       INTENT(IN) :: a(:), b(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
  CHARACTER(LEN=maxlen) :: tag,  tnm(3), sub
  REAL    :: aErr, rErr
  tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
  tnm = ['a  ','b  ','a-b'];  IF(PRESENT(nam    )) tnm = nam                   !--- Variables names
  sub = 'checkEqu';           IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)
  epsAbs = abs_err;           IF(PRESENT(absErr )) epsAbs = absErr             !--- Absolute error
  epsRel = rel_err;           IF(PRESENT(relErr )) epsRel = relErr             !--- Relative error
  lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax)
END FUNCTION checkEqu_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkEqu_2D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr)
  REAL,                       INTENT(IN) :: a(:,:), b(:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
  CHARACTER(LEN=maxlen) :: tag,  tnm(3), sub=''
  REAL    :: aErr, rErr
  tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//TRIM(tag); END IF
  tnm = ['a  ','b  ','a-b'];  IF(PRESENT(nam    )) tnm = nam                   !--- Variables names
  sub = 'checkEqu';           IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)
  epsAbs = abs_err;           IF(PRESENT(absErr )) epsAbs = absErr             !--- Absolute error
  epsRel = rel_err;           IF(PRESENT(relErr )) epsRel = relErr             !--- Relative error
  lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax)
END FUNCTION checkEqu_2D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkEqu_3D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr)
  REAL,                       INTENT(IN) :: a(:,:,:), b(:,:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
  CHARACTER(LEN=maxlen) :: tag,  tnm(3), sub=''
  REAL    :: aErr, rErr
  tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//TRIM(tag); END IF
  tnm = ['a  ','b  ','a-b'];  IF(PRESENT(nam    )) tnm = nam                   !--- Variables names
  sub = 'checkEqu';           IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)
  epsAbs = abs_err;           IF(PRESENT(absErr )) epsAbs = absErr             !--- Absolute error
  epsRel = rel_err;           IF(PRESENT(relErr )) epsRel = relErr             !--- Relative error
  lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax)
END FUNCTION checkEqu_3D
!===============================================================================================================================


!===============================================================================================================================
!=== CHECK FOR MASS CONSERVATION, IE. WHETHER THE CONCENTRATION OF A TRACER IS EQUAL TO THE SUM OF ITS CHILDREN CONCENTRATIONS.
!=== EXCEPTION (PROBABLY TO BE SUPPRESSED): GENERATION 1 TRACERS ARE COMPARED TO THE MAJOR NEXT GEBNERATION ISOTOPES ONLY
!===
!=== lerr = checkMass([qt][, err_tag][, subn][, nmax][, absErr][, relErr])
!===
!=== ARGUMENTS:    Meaning                                                          Intent  Optional?  Default value
!===   * qt        Tracers+isotopes+tags stacked along last index, rank 1 to 3      INPUT   OPTIONAL   qx (RANK 3 CASE ONLY !)
!===   * err_tag   Error message to display if fields are not matching              INPUT   OPTIONAL   "mismatching values"
!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkEqu"
!===   * nmax      Maximum number of outliers values to compute for each tracer     INPUT   OPTIONAL   All values
!===   * absErr    Maximum value for |q-r| to consider "q" and "q" as equal.        INPUT   OPTIONAL    abs_err
!===   * relErr    Maximum value for |q-r|/max(|q|,|r|,1e-18)      " "              INPUT   OPTIONAL    rel_err
!===
!=== REMARKS:
!===   * The concentration of each tracer is compared to the sum of its childrens concentrations (they must match).
!===   * In the isotopes case, only a subset of childrens is used (usually only the major contributor).
!===   * "qt" last dimension size can either be:
!===      - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes (single unknown phase)
!===      - nqtot   ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx").
!===============================================================================================================================
LOGICAL FUNCTION checkMassQ_0D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr)
  REAL,                       INTENT(IN) :: qt(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
  CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3)
  REAL                  :: qs, qp
  INTEGER, ALLOCATABLE :: iqChld(:)
  INTEGER :: iIso, iPha,  nqChld, iq, it
  tag = '';           IF(PRESENT(err_tag)) tag = err_tag                       !--- Tag for error message
  sub  = 'checkMass'; IF(PRESENT(subn   )) sub = TRIM(sub)//':'//subn          !--- Calling subroutine chain
  epsAbs = abs_err;   IF(PRESENT(absErr )) epsAbs = absErr                     !--- Absolute error (for isEqual routine)
  epsRel = rel_err;   IF(PRESENT(relErr )) epsRel = relErr                     !--- Relative error (for isEqual routine)
  IF(SIZE(qt) == ntiso+1) THEN
     iqChld = iMajr(it)%i
     pNam = isotope%parent
     DO it = 0, niso
        IF(it /= 0) iqChld = itZonIso(:,it)
        IF(it /= 0) pNam = isoName(it)
        qp = qt(it+1)
        qs = SUM(qt(1+iqChld), DIM=1)                                          !--- Tracer compared to its major isotopes
        sNam = TRIM(strStack(isoName(iqChld),'+'))                             !--- Names of contributing childrens
        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
        lerr = dispOutliers([isEqual(qp, qs)], cat(qp, qs, qp-qs), SHAPE(qt), tag, nam, sub, nmax); IF(lerr) RETURN
     END DO
  ELSE
     DO iq = 1, SIZE(qt)
        pNam   = tracers(iq)%name                                              !--- Current tracer name with phase
        nqChld = tracers(iq)%nqChildren                                        !--- Number of childrens (next generation only)
        IF(nqChld == 0) CYCLE                                                  !--- No descendants
        iIso = tracers(iq)%iso_iGroup                                          !--- Isotopes family index
        iPha = tracers(iq)%iso_iPhase                                          !--- Index of current phase in isoPhases(:)
        IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha)                   !--- Indices of major contributors in qt
        IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld)                  !--- Indices of every contributor  in qt
        qp = qt(iq)                                                            !--- Linearized (reprofiled) parent field
        qs = SUM(qt(iqChld))                                                   !--- Sum of contributing child(ren) field(s)
        sNam = TRIM(strStack(tracers(iqChld)%name ,'+'))                       !--- Names of contributing child(ren)
        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
        lerr = dispOutliers([isEqual(qp, qs)], cat(qp, qs, qp-qs), [1], tag, nam, sub, nmax); IF(lerr) RETURN
     END DO
  END IF
END FUNCTION checkMassQ_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkMassQ_1D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr)
  REAL,             TARGET,   INTENT(IN) :: qt(:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
  CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3)
  REAL, DIMENSION(SIZE(qt(:,1))) :: qs, qp
  INTEGER, ALLOCATABLE :: iqChld(:)
  INTEGER :: iIso, iPha,  nqChld, iq, it
  tag = '';           IF(PRESENT(err_tag)) tag = err_tag                       !--- Tag for error message
  sub  = 'checkMass'; IF(PRESENT(subn   )) sub = TRIM(sub)//':'//subn          !--- Calling subroutine chain
  epsAbs = abs_err;   IF(PRESENT(absErr )) epsAbs = absErr                     !--- Absolute error (for isEqual routine)
  epsRel = rel_err;   IF(PRESENT(relErr )) epsRel = relErr                     !--- Relative error (for isEqual routine)
  IF(SIZE(qt, DIM=2) == ntiso+1) THEN
     iqChld = iMajr(it)%i
     pNam = isotope%parent
     DO it = 0, niso
        IF(it /= 0) iqChld = itZonIso(:,it)
        IF(it /= 0) pNam = isoName(it)
        qp = PACK(qt(:,it+1), MASK=.TRUE.)
        qs = SUM(qt(:,1+iqChld), DIM=2)                                        !--- Tracer compared to its major isotopes
        sNam = TRIM(strStack(isoName(iqChld),'+'))                             !--- Names of contributing childrens
        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
        lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,1)), tag, nam, sub, nmax); IF(lerr) RETURN
     END DO
  ELSE
     DO iq = 1, SIZE(qt, DIM=2)
        pNam   = tracers(iq)%name                                              !--- Current tracer name with phase
        nqChld = tracers(iq)%nqChildren                                        !--- Number of childrens (next generation only)
        IF(nqChld == 0) CYCLE                                                  !--- No descendants
        iIso = tracers(iq)%iso_iGroup                                          !--- Isotopes family index
        iPha = tracers(iq)%iso_iPhase                                          !--- Index of current phase in isoPhases(:)
        IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha)                   !--- Indices of major contributors in qt
        IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld)                  !--- Indices of every contributor  in qt
        qp = PACK(qt(:,iq), MASK=.TRUE.)                                       !--- Linearized (reprofiled) parent field
        qs = PACK(SUM(qt(:,iqChld), DIM=2), MASK=.TRUE.)                       !--- Sum of contributing child(ren) field(s)
        sNam = TRIM(strStack(tracers(iqChld)%name ,'+'))                       !--- Names of contributing child(ren)
        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
        nam(1) = pnam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
       lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,1)), tag, nam, sub, nmax); IF(lerr) RETURN
     END DO
  END IF
END FUNCTION checkMassQ_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkMassQ_2D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr)
  REAL,             TARGET,   INTENT(IN) :: qt(:,:,:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
  CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3)
  REAL, DIMENSION(SIZE(qt(:,:,1))) :: qs, qp
  INTEGER, ALLOCATABLE :: iqChld(:)
  INTEGER :: iIso, iPha,  nqChld, iq, it
  tag = '';           IF(PRESENT(err_tag)) tag = err_tag                       !--- Tag for error message
  sub  = 'checkMass'; IF(PRESENT(subn   )) sub = TRIM(sub)//':'//subn          !--- Calling subroutine chain
  epsAbs = abs_err;   IF(PRESENT(absErr )) epsAbs = absErr                     !--- Absolute error (for isEqual routine)
  epsRel = rel_err;   IF(PRESENT(relErr )) epsRel = relErr                     !--- Relative error (for isEqual routine)
  IF(SIZE(qt, DIM=3) == ntiso+1) THEN
     iqChld = iMajr(it)%i
     pNam = isotope%parent
     DO it = 0, niso
        IF(it /= 0) iqChld = itZonIso(:,it)
        IF(it /= 0) pNam = isoName(it)
        qp = PACK(qt(:,:,it+1), MASK=.TRUE.)
        qs = PACK(SUM(qt(:,:,1+iqChld), DIM=3), MASK=.TRUE.)                   !--- Tracer compared to its major isotopes
        sNam = TRIM(strStack(isoName(iqChld),'+'))                             !--- Names of contributing childrens
        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
        lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,:,1)), tag, nam, sub, nmax); IF(lerr) RETURN
     END DO
  ELSE
     DO iq = 1, SIZE(qt, DIM=3)
        pNam   = tracers(iq)%name                                              !--- Current tracer name with phase
        nqChld = tracers(iq)%nqChildren                                        !--- Number of childrens (next generation only)
        IF(nqChld == 0) CYCLE                                                  !--- No descendants
        iIso = tracers(iq)%iso_iGroup                                          !--- Isotopes family index
        iPha = tracers(iq)%iso_iPhase                                          !--- Index of current phase in isoPhases(:)
        IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha)                   !--- Indices of major contributors in qt
        IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld)                  !--- Indices of every contributor  in qt
        qp = PACK(qt(:,:,iq), MASK=.TRUE.)                                     !--- Linearized (reprofiled) parent field
        qs = PACK(SUM(qt(:,:,iqChld), DIM=3), MASK=.TRUE.)                     !--- Sum of contributing child(ren) field(s)
        sNam = TRIM(strStack(tracers(iqChld)%name ,'+'))                       !--- Names of contributing child(ren)
        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
        lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,:,1)), tag, nam, sub, nmax); IF(lerr) RETURN
     END DO
  END IF
END FUNCTION checkMassQ_2D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION checkMassQx_2D(err_tag, subn, nmax, absErr, relErr) RESULT(lerr)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
  lerr = checkMassQ_2D(qx, err_tag, subn, nmax, absErr, relErr)
END FUNCTION checkMassQx_2D
!===============================================================================================================================


!===============================================================================================================================
!===============================================================================================================================
!=== COMPUTE THE ISOTOPIC ANOMALY "DELTA" AND CHECK THE VALUES
!===
!=== lerr = deltaR(rIso, iIso[, err_tag][, subn][, nmax][, delMax][, delIso]                [, lCheck])
!=== lerr = deltaQ(qt,   iIso[, err_tag][, subn][, nmax][, delMax][, delIso][, Phas][, qmin][, lCheck])
!=== lerr = deltaQ(      iIso[, err_tag][, subn][, nmax][, delMax][, delIso][, Phas][, qmin][, lCheck])
!===
!=== ARGUMENTS:    Meaning                                                          Intent  Optional?  Default      Present:
!===   * rIso      Field of rank 0 to 3                                             INPUT   MANDATORY               CASE 1 ONLY
!===   * qt        Stacked fields, rank 1 to 3 (last index for species)             INPUT   MANDATORY               CASE 2 ONLY
!===   * iIso      Index of tested isotope in "isoName"                             INPUT   MANDATORY
!===   * err_tag   Error message to display if fields are not matching              INPUT   OPTIONAL   "??????????"
!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "delta[R]"
!===   * nmax      Maximum number of lines to print for outliers                    INPUT   OPTIONAL   <all lines>  RANK>0 ONLY
!===   * delMax    Maximum value for the isotopic delta                             INPUT   OPTIONAL    dmax(iIso)
!===   * delIso    Isotopic anomaly delta for isotope iIso                          OUTPUT  OPTIONAL
!===   * phas      Phase                                                            INPUT   OPTIONAL                CASE 2 ONLY
!===   * qmin      Values lower than this threshold are not checked                 INPUT   OPTIONAL    small       CASE 2 ONLY
!===   * lCheck    Trigger checking routines with tables for strage values          INPUT   OPTIONAL   .FALSE.
!===
!=== REMARKS:
!===   * In the version 2 and 3 of the routine, values are checked only if q>=qmin
!===   * In the version 2, the size of the last dimension of "qt" can either be:
!===     - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes.
!===     - nqtot   ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx").
!===============================================================================================================================
LOGICAL FUNCTION deltaR_0D(rIso, iIso, err_tag, subn, delMax, delIso, lCheck) RESULT(lerr)
  REAL,                       INTENT(IN)  :: rIso
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn
  REAL,             OPTIONAL, INTENT(IN)  :: delMax
  REAL,             OPTIONAL, INTENT(OUT) :: delIso
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss
  REAL    :: dmn, dmx, dIso
  LOGICAL :: lc
  IF(niso == 0) RETURN
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
  dmn = dMin(ixIso)%r(iIso)
  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Check the index + define few vars
  mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
  tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                                 !=== Compute delta(isotope iIso in phase iPhas)
  IF(PRESENT(delIso)) delIso = dIso
  IF(.NOT.lc)   RETURN
  lerr = dispOutliers([dIso<dmn .OR. dIso>dmx], cat(rIso, dIso), [1], mss, vNamDe(iIso), sub)
END FUNCTION deltaR_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION deltaR_1D(rIso, iIso, err_tag, subn, nmax, delMax, delIso, lCheck) RESULT(lerr)
  REAL,                       INTENT(IN)  :: rIso(:)
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn
  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
  REAL,             OPTIONAL, INTENT(IN)  :: delMax
  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(rIso,1))
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss
  REAL    :: dmn, dmx, dIso(SIZE(rIso,1))
  LOGICAL :: lc
  IF(niso == 0) RETURN
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
  dmn = dMin(ixIso)%r(iIso)
  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Check the index + define few vars
  mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
  tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                                 !=== Compute delta(isotope iIso in phase iPhas)
  IF(PRESENT(delIso)) delIso = dIso
  IF(.NOT.lc)   RETURN
  lerr = dispOutliers(PACK(dIso<dmn .OR. dIso>dmx,.TRUE.), cat(rIso, dIso), SHAPE(dIso), mss, vNamDe(iIso), sub, nmax)
END FUNCTION deltaR_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION deltaR_2D(rIso, iIso, err_tag, subn, nmax, delMax, delIso, lCheck) RESULT(lerr)
  REAL,                       INTENT(IN)  :: rIso(:,:)
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn
  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
  REAL,             OPTIONAL, INTENT(IN)  :: delMax
  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(rIso,1),SIZE(rIso,2))
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss
  REAL    :: dmn, dmx, emx, dIso(SIZE(rIso,1),SIZE(rIso,2))
  LOGICAL :: lc
  IF(niso == 0) RETURN
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
  dmn = dMin(ixIso)%r(iIso)
  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Check the index + define few vars
  mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
  tnat0 = tnat(iIso); dIso = ratio2delta(rIso)
  IF(PRESENT(delIso)) delIso = dIso
  IF(.NOT.lc)   RETURN
  lerr = dispOutliers( [dIso<dmn.OR.dIso>dmx], cat(PACK(rIso,.TRUE.),PACK(dIso,.TRUE.)), SHAPE(dIso),mss,vNamDe(iIso),sub,nmax)
END FUNCTION deltaR_2D
!===============================================================================================================================
LOGICAL FUNCTION deltaQ_0D(qt, iIso, err_tag, subn, delMax, delIso, phas, qmin, lCheck) RESULT(lerr)
  REAL, DIMENSION(:),         INTENT(IN)  :: qt
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
  REAL,             OPTIONAL, INTENT(IN)  :: delMax, qmin
  REAL,             OPTIONAL, INTENT(OUT) :: delIso
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
  REAL    :: dmn, dmx, qmn, rIso, dIso
  LOGICAL :: lc, lq
  INTEGER :: iqIso, iqPar, iZon, iPha, ii
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
  lq = SIZE(qt, DIM=1) == niso+1

  !=== CHECK PHASE
  IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
  IF(lerr)       RETURN
  dmn = dMin(ixIso)%r(iIso)
  m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'

  !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
  ii = iIso
  DO iZon = 0, nzone
     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     DO iPha = 1, nphas
        p = isoPhas(iPha:iPha)
        lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN    !--- Ckeck the index + define few vars
        mss = TRIM(mss)//TRIM(m)
        IF(     lq) rIso = qt(iIso+1)/qt(1)                                    !--- "qt(1+niso)": parent, then isotopes
        IF(.NOT.lq) rIso = qt(iqIso)/qt(iqPar)                                 !--- Fields taken from a "qt" similar to "qx"
        tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                           !=== Compute delta(isotope iIso in phase "p")
        IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso       !--- Return  delta(iIso, phase=p)
        IF(.NOT.lc) CYCLE
        lerr = dispOutliers( [qt(iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], cat(rIso, dIso), [1], mss, vNamde(iIso), sub)
        IF(lerr) RETURN
     END DO
     IF(.NOT.lc) EXIT
  END DO
END FUNCTION deltaQ_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION deltaQ_1D(qt, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr)
  REAL, DIMENSION(:,:),       INTENT(IN)  :: qt
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
  REAL,             OPTIONAL, INTENT(IN)  :: delMax, qmin
  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(qt,1))
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
  REAL    :: dmn, dmx, qmn
  REAL, DIMENSION(SIZE(qt,1)) :: rIso, dIso
  LOGICAL :: lc, lq
  INTEGER :: iqIso, iqPar, iZon, iPha, ii
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
  lq = SIZE(qt, DIM=2) == niso+1

  !=== CHECK PHASE
  IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
  IF(lerr)       RETURN
  dmn = dMin(ixIso)%r(iIso)
  m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'

  !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
  ii = iIso
  DO iZon = 0, nzone
     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     DO iPha = 1, nphas
        p = isoPhas(iPha:iPha)
        lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN    !--- Ckeck the index + define few vars
        mss = TRIM(mss)//TRIM(m)
        IF(     lq) rIso = qt(:,iIso+1)/qt(:,1)                                !--- "qt(1+niso)": parent, then isotopes
        IF(.NOT.lq) rIso = qt(:,iqIso)/qt(:,iqPar)                             !--- Fields taken from a "qt" similar to "qx"
        tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                           !=== Compute delta(iIso, phase=p)
        IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso       !--- Return  delta(iIso, phase=p)
        IF(.NOT.lc) CYCLE
        lerr = dispOutliers( [qt(:,iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], &
                        cat(PACK(rIso,.TRUE.), PACK(dIso,.TRUE.)), SHAPE(dIso), mss, vNamde(iIso), sub, nmax)
        IF(lerr) RETURN
     END DO
     IF(.NOT.lc) EXIT
  END DO
END FUNCTION deltaQ_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION deltaQ_2D(qt, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr)
  REAL, DIMENSION(:,:,:),     INTENT(IN)  :: qt
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
  REAL,             OPTIONAL, INTENT(IN)  :: delMax, qmin
  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(qt,1),SIZE(qt,2))
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
  REAL    :: dmn, dmx, qmn
  REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)) :: rIso, dIso
  LOGICAL :: lc, lq
  INTEGER :: iqIso, iqPar, iZon, iPha, ii
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
  lq = SIZE(qt, DIM=3) == niso+1

  !=== CHECK PHASE
  IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
  IF(lerr)       RETURN
  dmn = dMin(ixIso)%r(iIso)
  m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'

  !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
  ii = iIso
  DO iZon = 0, nzone
     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     DO iPha = 1, nphas
        p = isoPhas(iPha:iPha)
        lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN    !--- Ckeck the index + define few vars
        mss = TRIM(mss)//' '//TRIM(m)
        IF(     lq) rIso = qt(:,:,iIso+1)/qt(:,:,1)                            !--- "qt(1+niso)": parent, then isotopes
        IF(.NOT.lq) rIso = qt(:,:,iqIso)/qt(:,:,iqPar)                         !--- Fields taken from a "qt" similar to "qx"
        tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                           !=== Compute delta(iIso, phase=p)
        IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso       !--- Return  delta(iIso, phase=p)
        IF(.NOT.lc) CYCLE
        lerr = dispOutliers( [qt(:,:,iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], &
                        cat(PACK(rIso,.TRUE.), PACK(dIso,.TRUE.)), SHAPE(dIso), mss, vNamDe(iIso), sub, nmax)
        IF(lerr) RETURN
     END DO
     IF(.NOT.lc) EXIT
  END DO
END FUNCTION deltaQ_2D
!===============================================================================================================================
LOGICAL FUNCTION deltaQx_2D(iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr)
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
  REAL,             OPTIONAL, INTENT(IN)  :: delMax, qmin
  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(qx,1),SIZE(qx,2))
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  lerr = deltaQ_2D(qx, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck)
END FUNCTION deltaQx_2D
!===============================================================================================================================


!===============================================================================================================================
!===============================================================================================================================
!=== COMPUTE THE ISOTOPIC EXCESS AND CHECK THE VALUE (D-excess and O17-excess so far)
!===
!=== lerr = excess(rIso,rRef,iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef]              [,lCheck])
!=== lerr = excess(qt,       iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef][,phas][,qmin][,lCheck])
!=== lerr = excess(          iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef][,phas][,qmin][,lCheck])
!===
!=== ARGUMENTS:                                                                     Intent  Optional?  Default value
!===   * rIso,rRef Considered isotope + reference ratios: field of rank 0 to 3      INPUT   MANDATORY              (CASE 1 ONLY)
!===   * qt        Stacked fields, rank 1 to 3 (last index for species)             INPUT   MANDATORY              (CASE 2 ONLY)
!===   * iIso      Index of tested species                                          INPUT   MANDATORY
!===   * err_tag   Error tag sent from the calling routine                          INPUT   OPTIONAL    none
!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "excess"
!===   * nmax      Maximum number of lines to print for outliers                    INPUT   OPTIONAL   <all lines> (RANK>0 ONLY)
!===   * delMax    Maximum value for the isotopic delta                             INPUT   OPTIONAL
!===   * excMax    Maximum value for the isotopic excess                            INPUT   OPTIONAL    dmax(iIso)
!===   * delIso    Isotopic anomaly delta for isotope iIso                          OUTPUT  OPTIONAL
!===   * excIso    Isotopic excess for isotope iIso                                 OUTPUT  OPTIONAL
!===   * delRef    Isotopic anomaly delta for isotope iIso reference                OUTPUT  OPTIONAL
!===   * phas      Phase name (one element of "isoPhas")                            INPUT   MANDATORY              (CASE 2 ONLY)
!===   * qmin      Values lower than this threshold are not checked                 INPUT   OPTIONAL    small      (CASE 2 ONLY)
!===
!=== REMARKS:
!===   * In the version 2 and 3 of the routine, values are checked only if q>=qmin
!===   * In the version 2, the size of the last dimension of "qt" can either be:
!===     - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes.
!===     - nqtot   ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx").
!===============================================================================================================================
LOGICAL FUNCTION excessR_0D(rIso, rRef, iIso, err_tag, subn, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr)
  REAL,                       INTENT(IN)  :: rIso, rRef
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn
  REAL,             OPTIONAL, INTENT(IN)  :: delMax, excMax
  REAL,             OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss
  REAL    :: deIso, deRef, exIso, dmx, drx, emn, emx
  LOGICAL :: lc
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
  sub = 'excess';            IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)  !--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Ckeck the index + define few vars
  emn = eMin(ixIso)%r(iIso)
  drx = dMax(ixIso)%r(iso_ref)
  mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
  lerr = deltaR_0D(rIso, iIso,    tag, sub, dmx, deIso, lCheck) &              !=== COMPUTE deltaIso AND deltaRef
    .OR. deltaR_0D(rRef, iso_ref, tag, sub, drx, deRef, lCheck)
  IF(lerr)    RETURN
  exIso = isoExcess(deIso, deRef)                                              !=== COMPUTE THE EXCESS
  IF(PRESENT(delIso)) delIso = deIso
  IF(PRESENT(excIso)) excIso = exIso
  IF(PRESENT(delRef)) delRef = deRef
  IF(.NOT.lc) RETURN
  lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), [1], mss, vNamEx(iIso), sub)
END FUNCTION excessR_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION excessR_1D(rIso, rRef, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr)
  REAL, DIMENSION(:),                      INTENT(IN)  :: rIso, rRef
  INTEGER,                                 INTENT(IN)  :: iIso
  CHARACTER(LEN=*),              OPTIONAL, INTENT(IN)  :: err_tag, subn
  INTEGER,                       OPTIONAL, INTENT(IN)  :: nmax
  REAL,                          OPTIONAL, INTENT(IN)  :: delMax, excMax
  REAL, DIMENSION(SIZE(rIso,1)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
  LOGICAL,                       OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss, m
  REAL, DIMENSION(SIZE(rIso,1)) :: deIso, deRef, exIso
  REAL    :: dmx, drx, emn, emx
  LOGICAL :: lc
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
  sub = 'excess';            IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)  !--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Ckeck the index + define few vars
  emn = eMin(ixIso)%r(iIso)
  drx = dMax(ixIso)%r(iso_ref)
  mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
  lerr = deltaR_1D(rIso, iIso,    tag, sub, nmax, dmx, deIso, lCheck) &        !=== COMPUTE deltaIso AND deltaRef
    .OR. deltaR_1D(rRef, iso_ref, tag, sub, nmax, drx, deRef, lCheck)
  IF(lerr)    RETURN
  exIso = isoExcess(deIso, deRef)                                              !=== COMPUTE THE EXCESS
  IF(PRESENT(delIso)) delIso = deIso
  IF(PRESENT(excIso)) excIso = exIso
  IF(PRESENT(delRef)) delRef = deRef
  IF(.NOT.lc) RETURN
  lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), SHAPE(exIso), mss, vNamEx(iIso), sub, nmax)
END FUNCTION excessR_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION excessR_2D(rIso, rRef, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr)
  REAL, DIMENSION(:,:),                                 INTENT(IN)  :: rIso, rRef
  INTEGER,                                              INTENT(IN)  :: iIso
  CHARACTER(LEN=*),                           OPTIONAL, INTENT(IN)  :: err_tag, subn
  INTEGER,                                    OPTIONAL, INTENT(IN)  :: nmax
  REAL,                                       OPTIONAL, INTENT(IN)  :: delMax, excMax
  REAL, DIMENSION(SIZE(rIso,1),SIZE(rIso,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
  LOGICAL,                                    OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss, m
  REAL, DIMENSION(SIZE(rIso,1),SIZE(rIso,2)) :: deIso, deRef, exIso
  REAL    :: dmx, drx, emn, emx
  LOGICAL :: lc
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
  sub = 'excess';            IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)  !--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Ckeck the index + define few vars
  emn = eMin(ixIso)%r(iIso)
  drx = dMax(ixIso)%r(iso_ref)
  mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
  lerr = deltaR_2D(rIso, iIso,    tag, sub, nmax, dmx, deIso, lCheck) &        !=== COMPUTE deltaIso AND deltaRef
    .OR. deltaR_2D(rRef, iso_ref, tag, sub, nmax, drx, deRef, lCheck)
  IF(lerr)    RETURN
  exIso = isoExcess(deIso, deRef)                                              !=== COMPUTE THE EXCESS
  IF(PRESENT(delIso)) delIso = deIso
  IF(PRESENT(excIso)) excIso = exIso
  IF(PRESENT(delRef)) delRef = deRef
  IF(.NOT.lc) RETURN
  lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(PACK(deIso,.TRUE.), PACK(deRef,.TRUE.), PACK(exIso,.TRUE.)), &
                     SHAPE(exIso), TRIM(mss)//' '//TRIM(m), vNamEx(iIso), sub, nmax)
END FUNCTION excessR_2D
!===============================================================================================================================
LOGICAL FUNCTION excessQ_0D(qt, iIso, err_tag, subn, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck) RESULT(lerr)
  REAL,                       INTENT(IN)  :: qt(:)
  INTEGER,                    INTENT(IN)  :: iIso
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
  REAL,             OPTIONAL, INTENT(IN)  :: delMax, excMax, qmin
  REAL,             OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
  REAL    :: deIso, exIso, deRef
  REAL    :: dmx, drx, emn, emx, qmn
  INTEGER :: iZon, iPha, ii
  LOGICAL :: lc
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
  sub = 'excess';            IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger

  !=== CHECK PHASE
  lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
  IF(lerr)       RETURN
  emn = eMin(ixIso)%r(iIso)
  drx = dMax(ixIso)%r(iso_ref)
  m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'

  !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
  ii = iIso
  DO iZon = 0, nzone
     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN                  !--- Check ii, set iso_ref + excess params
     mss = TRIM(mss)//TRIM(m)
     DO iPha = 1, nphas
        p = isoPhas(iPha:iPha)                                                 !--- Current phase
        lerr = deltaQ_0D(qt, ii,      tag, sub, dmx, deIso, p, qmn, lCheck) &  !=== COMPUTE deltaIso AND deltaRef
          .OR. deltaQ_0D(qt, iso_ref, tag, sub, drx, deRef, p, qmn, lCheck)
        exIso = isoExcess(deIso, deRef)                                        !=== COMPUTE THE EXCESS
        IF(iZon == 0 .AND. p == pha) THEN
           IF(PRESENT(delIso)) delIso = deIso                                  !--- Return delta(iIso,   phase=p)
           IF(PRESENT(delRef)) delRef = deRef                                  !--- Return delta(iso_ref,phase=p)
           IF(PRESENT(excIso)) excIso = exIso                                  !--- Return excess(iIso,  phase=p)
        END IF
        IF(.NOT.lc) CYCLE
        lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), [1], mss, vNamEx(ii, iPha, iZon), sub)
        IF(lerr) RETURN
     END DO
     IF(.NOT.lc) EXIT
  END DO
END FUNCTION excessQ_0D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION excessQ_1D(qt, iIso, err_tag, subn, nmax, delMax,excMax,delIso,excIso,delRef, phas, qmin, lCheck) RESULT(lerr)
  REAL,                                  INTENT(IN)  :: qt(:,:)
  INTEGER,                               INTENT(IN)  :: iIso
  CHARACTER(LEN=*),            OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
  INTEGER,                     OPTIONAL, INTENT(IN)  :: nmax
  REAL,                        OPTIONAL, INTENT(IN)  :: delMax, excMax, qmin
  REAL, DIMENSION(SIZE(qt,1)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
  LOGICAL,                     OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
  REAL, DIMENSION(SIZE(qt,1)) :: deIso, exIso, deRef
  REAL    :: dmx, drx, emn, emx, qmn
  INTEGER :: iZon, iPha, ii
  LOGICAL :: lc
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
  sub = 'excess';            IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger

  !=== CHECK PHASE
  lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
  IF(lerr)       RETURN
  emn = eMin(ixIso)%r(iIso)
  drx = dMax(ixIso)%r(iso_ref)
  m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'

  !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
  ii = iIso
  DO iZon = 0, nzone
     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN                  !--- Check ii, set iso_ref + excess params
     mss = TRIM(mss)//TRIM(m)
     DO iPha = 1, nphas
        p = isoPhas(iPha:iPha)                                                 !--- Current phase
        lerr = deltaQ_1D(qt, ii,      tag,sub,nmax, dmx, deIso,p,qmn,lCheck) & !=== COMPUTE deltaIso AND deltaRef
          .OR. deltaQ_1D(qt, iso_ref, tag,sub,nmax, drx, deRef,p,qmn,lCheck)
        exIso = isoExcess(deIso, deRef)                                        !=== COMPUTE THE EXCESS
        IF(iZon == 0 .AND. p == pha) THEN
           IF(PRESENT(delIso)) delIso = deIso                                  !--- Return delta(iIso,   phase=p)
           IF(PRESENT(delRef)) delRef = deRef                                  !--- Return delta(iso_ref,phase=p)
           IF(PRESENT(excIso)) excIso = exIso                                  !--- Return excess(iIso,  phase=p)
        END IF
        IF(.NOT.lc) CYCLE
        lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), &
                             SHAPE(exIso), mss, vNamEx(ii, iPha, iZon), sub, nmax)
        IF(lerr) RETURN
     END DO
     IF(.NOT.lc) EXIT
  END DO
END FUNCTION excessQ_1D
!-------------------------------------------------------------------------------------------------------------------------------
LOGICAL FUNCTION excessQ_2D(qt, iIso, err_tag, subn, nmax, delMax,excMax,delIso,excIso,delRef, phas, qmin, lCheck) RESULT(lerr)
  REAL,                                             INTENT(IN)  :: qt(:,:,:)
  INTEGER,                                          INTENT(IN)  :: iIso
  CHARACTER(LEN=*),                       OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
  INTEGER,                                OPTIONAL, INTENT(IN)  :: nmax
  REAL,                                   OPTIONAL, INTENT(IN)  :: delMax, excMax, qmin
  REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
  LOGICAL,                                OPTIONAL, INTENT(IN)  :: lCheck
  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
  REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)) :: deIso, exIso, deRef
  REAL    :: dmx, drx, emn, emx, qmn
  INTEGER :: iZon, iPha, ii
  LOGICAL :: lc
  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
  sub = 'excess';            IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger

  !=== CHECK PHASE
  lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
  IF(lerr)       RETURN
  emn = eMin(ixIso)%r(iIso)
  drx = dMax(ixIso)%r(iso_ref)
  m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'

  !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
  ii = iIso
  DO iZon = 0, nzone
     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN                  !--- Check ii, set iso_ref + excess params
     mss = TRIM(mss)//TRIM(m)
     DO iPha = 1, nphas
        p = isoPhas(iPha:iPha)                                                 !--- Current phase
        lerr = deltaQ_2D(qt, ii,      tag,sub,nmax, dmx, deIso,p,qmn,lCheck) & !=== COMPUTE deltaIso AND deltaRef
          .OR. deltaQ_2D(qt, iso_ref, tag,sub,nmax, drx, deRef,p,qmn,lCheck)
        exIso = isoExcess(deIso, deRef)                                        !=== COMPUTE THE EXCESS
        IF(iZon == 0 .AND. p == pha) THEN
           IF(PRESENT(delIso)) delIso = deIso                                  !--- Return delta(iIso,   phase=p)
           IF(PRESENT(delRef)) delRef = deRef                                  !--- Return delta(iso_ref,phase=p)
           IF(PRESENT(excIso)) excIso = exIso                                  !--- Return excess(iIso,  phase=p)
        END IF
        IF(.NOT.lc) CYCLE
        lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(PACK(deIso,.TRUE.), PACK(deRef,.TRUE.), PACK(exIso,.TRUE.)), &
                             SHAPE(exIso), mss, vNamEx(ii, iPha, iZon), sub, nmax)
        IF(lerr) RETURN
     END DO
     IF(.NOT.lc) EXIT
  END DO
END FUNCTION excessQ_2D
!===============================================================================================================================
LOGICAL FUNCTION excessQx_2D(iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck) RESULT(lerr)
  INTEGER,                                          INTENT(IN)  :: iIso
  CHARACTER(LEN=*),                       OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
  INTEGER,                                OPTIONAL, INTENT(IN)  :: nmax
  REAL,                                   OPTIONAL, INTENT(IN)  :: delMax, excMax, qmin
  REAL, DIMENSION(SIZE(qx,1),SIZE(qx,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
  LOGICAL,                                OPTIONAL, INTENT(IN)  :: lCheck
  lerr = excessQ_2D(qx, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck)
END FUNCTION excessQx_2D
!===============================================================================================================================


!===============================================================================================================================
!=== GENERAL PURPOSE CHECKING ROUTINE:   lerr = checkTrac(err_tag[, subn][, nmax][, sCheck])
!===
!=== ARGUMENTS:                                                                     Intent  Optional?  Default value
!===   * err_tag   Error tag sent from the calling routine                          INPUT   OPTIONAL    none
!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkTrac"
!===   * nmax      Maximum number of lines to print for outliers                    INPUT   OPTIONAL    32
!===   * sCheck    Triggers for verifications                                       INPUT   OPTIONAL    sIsoChech
!===
!=== REMARKS:
!===   Tunable thresholds available in the individual routines (delMax, excMax, qmin, tny, absErr, relErr) have not been kept
!===   as optional arguments in this routine, because the adequate values are tracer or isotope-dependent.
!===   For special usages, a specific version can be written, or individual routines with special thresholds can be called.
!===============================================================================================================================
LOGICAL FUNCTION checkTrac(err_tag, subn, nmax, sCheck) RESULT(lerr)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, sCheck
  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
  INTEGER :: nmx, i, iPha, iIso, iq
  LOGICAL :: lNan, lPos, lMas, lDel, lExc, lBnd
  CHARACTER(LEN=maxlen) :: tag, sub, chk
  CHARACTER(LEN=maxlen), ALLOCATABLE :: tnam(:)

!!!!! GERER TNAM

  IF(niso == 0) RETURN
  tag = '';                 IF(PRESENT(err_tag)) tag = err_tag
  sub = 'checkTrac';        IF(PRESENT(subn   )) sub = TRIM(sub)//TRIM(subn)
  nmx = 32;                 IF(PRESENT(nmax   )) nmx = nmax
  chk =  sIsoCheck(ixIso) ; IF(PRESENT(sCheck )) chk = sCheck
  SELECT CASE(isotope%parent)
     CASE('H2O')
        DO iPha = 1, isotope%nphas
           DO iIso = 1, niso
              iq = iqIsoPha(iIso,iPha); tnam = tracers(iq)%name!; qmin = tracers(iq)%qmin; qmax = tracers(iq)%qmax
              IF(chk(1:1) == 'n') lerr = checkNaN (qx(:,:,iq),tag,tnam,sub,nmx); IF(lerr) RETURN
              IF(chk(2:2) == 'p') lerr = checkPos (qx(:,:,iq),tag,tnam,sub,nmx); IF(lerr) RETURN ! tny
              IF(chk(3:3) == 'm') lerr = checkMass(           tag,sub,     nmx); IF(lerr) RETURN ! absErr, relErr
             !IF(chk(6:6) == 'b') lerr = checkBnd (qx(:,:,iq),tag,tnam,sub,nmx,qmin,qmax); IF(lerr) RETURN
           END DO
           iIso = iso_HDO; IF(chk(4:4) == 'd') lerr =  delta(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN
           iIso = iso_O18; IF(chk(4:4) == 'd') lerr =  delta(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN
           iIso = iso_HDO; IF(chk(5:5) == 'e') lerr = excess(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN
           iIso = iso_O17; IF(chk(5:5) == 'e') lerr = excess(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN
        END DO
  END SELECT
END FUNCTION checkTrac
!===============================================================================================================================


  
        SUBROUTINE iso_verif_init()
        use ioipsl_getin_p_mod, ONLY : getin_p
        use isotopes_mod, ONLY: iso_O17, iso_O18, iso_HDO
        implicit none

        write(*,*) 'iso_verif_init 99: entree'
        deltalim=300.0
        deltalim_snow=500.0
        deltaDmin=-900.0
        deltalimtrac=2000.0
        O17_verif=.false.
        o17excess_bas=-200.0
        o17excess_haut=120.0
        dexcess_min=-100.0
        dexcess_max=1000.0

call getin_p('deltalim',deltalim)
deltalim_snow=deltalim+200.0
call getin_p('deltaDmin',deltaDmin)
call getin_p('deltalimtrac',deltalimtrac)

write(*,*) 'iso_O17,iso_O18,iso_HDO=',iso_O17,iso_O18,iso_HDO
if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
  call getin_p('O17_verif',O17_verif)  
  call getin_p('o17excess_bas',o17excess_bas)  
  call getin_p('o17excess_haut',o17excess_haut)  
  call getin_p('nlevmaxO17',nlevmaxO17)  
endif
if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
  call getin_p('dexcess_min',dexcess_min)  
  call getin_p('dexcess_max',dexcess_max)  
endif

write(*,*) 'deltalim=',deltalim
write(*,*) 'deltaDmin=',deltaDmin      
write(*,*) 'deltalimtrac=',deltalimtrac  
write(*,*) 'O17_verif=',O17_verif
write(*,*) 'o17excess_bas=',o17excess_bas
write(*,*) 'o17excess_haut=',o17excess_haut  
write(*,*) 'dexcess_min=',dexcess_min
write(*,*) 'dexcess_max=',dexcess_max

        end SUBROUTINE iso_verif_init

      subroutine iso_verif_egalite(a,b,err_msg)
        implicit none
        ! compare a et b. Si pas egal à errmax près, on affiche message
        ! d'erreur et stoppe

        ! input:
        real a, b
        character*(*) err_msg ! message d''erreur à afficher

        ! local
        !integer iso_verif_egalite_choix_nostop


        if (iso_verif_egalite_choix_nostop(a,b,err_msg, &
     &           errmax,errmaxrel).eq.1) then
                stop
        endif
        
#ifdef ISOVERIF
#else
        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
        stop
#endif           

        return
        end subroutine iso_verif_egalite

        !*****************

        function iso_verif_egalite_nostop(a,b,err_msg)
        implicit none
        ! compare a et b. Si pas egal à errmax près, on affiche message
        ! d'erreur et stoppe

        ! input:
        real a, b
        character*(*) err_msg ! message d''erreur à afficher
        !ouptut
        integer iso_verif_egalite_nostop
        ! local 
        !integer iso_verif_egalite_choix_nostop


        iso_verif_egalite_nostop=iso_verif_egalite_choix_nostop &
     &          (a,b,err_msg,errmax,errmaxrel)

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                  

        return
        end function iso_verif_egalite_nostop



        !************************************

        subroutine iso_verif_aberrant(R,err_msg)
        use isotopes_mod, ONLY: ridicule, iso_HDO
        implicit none
        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
        ! d'erreur et stoppe

        ! input:
        real R
        character*(*) err_msg ! message d''erreur à afficher
        !real deltaD
        !integer iso_verif_aberrant_choix_nostop


        if (iso_verif_aberrant_choix_nostop(R,1.0,ridicule, &
     &           deltalim,err_msg).eq.1) then
             stop
        endif 

#ifdef ISOVERIF
        if (.not.(iso_HDO.gt.0)) then    
         write(*,*) 'iso86: err_msg=',err_msg,': pourquoi verif?'
          stop
        endif
#else        
        write(*,*) 'iso 90: err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   
                
        end subroutine iso_verif_aberrant

        subroutine iso_verif_aberrant_encadre(R,err_msg)
        use isotopes_mod, ONLY: ridicule, iso_HDO
        implicit none
        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
        ! d'erreur et stoppe

        ! input:
        real R
        character*(*) err_msg ! message d''erreur à afficher
        !real deltaD

        !integer iso_verif_aberrant_enc_choix_nostop


        if (iso_verif_aberrant_enc_choix_nostop(R,1.0,ridicule, &
     &           deltalim,err_msg).eq.1) then
             write(*,*) 'deltaD=',deltaD(R)
             call abort_physic('isotopes_verif_mod > iso_verif_aberrant_encadre',err_msg,1)
             !stop             
        endif 

#ifdef ISOVERIF
        if (.not.(iso_HDO.gt.0)) then    
         write(*,*) 'iso86: err_msg=',err_msg,': pourquoi verif?'
          stop
        endif
#else        
        write(*,*) 'iso 90: err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   
        
        end subroutine iso_verif_aberrant_encadre

        !************************************

        subroutine iso_verif_aberrant_choix(xt,q,qmin,deltaDmax,err_msg)
        use isotopes_mod, ONLY: iso_HDO
        implicit none
        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
        ! d'erreur et stoppe

        ! input:
        real xt,q,qmin,deltaDmax
        character*(*) err_msg ! message d''erreur à afficher
        !real deltaD

        ! locals
        !integer iso_verif_aberrant_choix_nostop


        if (iso_verif_aberrant_choix_nostop(xt,q,qmin, &
     &           deltaDmax,err_msg).eq.1) then
             stop
        endif

#ifdef ISOVERIF
        if (.not.(iso_HDO.gt.0)) then    
         write(*,*) 'iso122: err_msg=',err_msg,': pourquoi verif?'
          stop
        endif
#else
        write(*,*) 'iso126: err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   
        
        end subroutine iso_verif_aberrant_choix

        !************************************

        function iso_verif_aberrant_nostop(R,err_msg)
        use isotopes_mod, ONLY: ridicule,iso_HDO
        implicit none
        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
        ! d'erreur et stoppe

        ! input:
        real R
        character*(*) err_msg ! message d''erreur à afficher
        integer iso_verif_aberrant_nostop ! output: 1 si erreur, 0 sinon
        !real deltaD

        ! locals
        !integer iso_verif_aberrant_choix_nostop

        iso_verif_aberrant_nostop=iso_verif_aberrant_choix_nostop &
     &           (R,1.0,ridicule,deltalim,err_msg)

#ifdef ISOVERIF
        if (.not.(iso_HDO.gt.0)) then    
         write(*,*) 'iso156: err_msg=',err_msg,': pourquoi verif?'
          stop
        endif
#else
        write(*,*) 'iso160: err_msg=',err_msg,': pourquoi verif?'
        stop
#endif           
        
        return
        end function iso_verif_aberrant_nostop

        function iso_verif_aberrant_enc_nostop(R,err_msg)
        use isotopes_mod, ONLY: ridicule,iso_HDO
        implicit none
        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
        ! d'erreur et stoppe

        ! input:
        real R
        character*(*) err_msg ! message d''erreur à afficher
        integer iso_verif_aberrant_enc_nostop ! output: 1 si erreur, 0 sinon
        !real deltaD

        ! locals
        !integer iso_verif_aberrant_enc_choix_nostop

        iso_verif_aberrant_enc_nostop= &
     &           iso_verif_aberrant_enc_choix_nostop &
     &           (R,1.0,ridicule,deltalim,err_msg)

#ifdef ISOVERIF
        if (.not.(iso_HDO.gt.0)) then    
         write(*,*) 'iso156: err_msg=',err_msg,': pourquoi verif?'
          stop
        endif
#else
        write(*,*) 'iso160: err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   
        return
        end function iso_verif_aberrant_enc_nostop

        !************************************

        function iso_verif_aberrant_choix_nostop(xt,q,   &
     &            qmin,deltaDmax,err_msg)

        use isotopes_mod, ONLY: iso_HDO
        implicit none
        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
        ! d'erreur et stoppe

        ! input:
        real xt,q,qmin,deltaDmax
        character*(*) err_msg ! message d''erreur à afficher
        ! output
        integer iso_verif_aberrant_choix_nostop
        ! locals
        !real deltaD
        !integer iso_verif_noNaN_nostop        


        iso_verif_aberrant_choix_nostop=0

#ifdef ISOVERIF        
        if (iso_verif_noNaN_nostop(q,err_msg).eq.1) then
            write(*,*) 'q=',q
            iso_verif_aberrant_choix_nostop=1
        endif      
        if (iso_verif_noNaN_nostop(xt,err_msg).eq.1) then
            write(*,*) 'xt=',xt
            iso_verif_aberrant_choix_nostop=1
        endif
#endif

        if (q.gt.qmin) then
        if ((deltaD(xt/q).gt.deltaDmax).or. &
     &          (deltaD(xt/q).lt.-borne).or. &
     &          (xt.lt.-borne).or. &
     &          (xt.gt.borne)) then                 
            write(*,*) 'erreur detectee par '// &
     &             'iso_verif_aberrant_choix_nostop:'
            write(*,*) err_msg
            write(*,*) 'q,deltaD=',q,deltaD(xt/q)
            write(*,*) 'deltaDmax=',deltaDmax
            iso_verif_aberrant_choix_nostop=1
            if (abs(xt-q)/q.lt.errmax) then
                write(*,*) 'attention, n''a-t-on pas confondu'// &
     &           ' iso_HDO et iso_eau dans l''appel à la verif?' 
            endif
        endif
        endif

#ifdef ISOVERIF
        if (.not.(iso_HDO.gt.0)) then    
         write(*,*) 'iso205: err_msg=',err_msg,': pourquoi verif?'
          stop
        endif
#else
        write(*,*) 'iso 209: err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   
        
        return
        end function iso_verif_aberrant_choix_nostop

        function iso_verif_aberrant_enc_choix_nostop(xt,q,   &
     &            qmin,deltaDmax,err_msg)
        use isotopes_mod, ONLY: iso_HDO
        implicit none
        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
        ! d'erreur et stoppe

        ! input:
        real xt,q,qmin,deltaDmax
        character*(*) err_msg ! message d''erreur à afficher
        ! output
        integer iso_verif_aberrant_enc_choix_nostop
        ! locals
        !real deltaD

        iso_verif_aberrant_enc_choix_nostop=0
        if (q.gt.qmin) then
        if ((deltaD(xt/q).gt.deltaDmax).or. &
     &          (deltaD(xt/q).lt.deltaDmin)) then                 
            write(*,*) 'erreur detectee par '// &
     &             'iso_verif_aberrant_choix_nostop:'
            write(*,*) err_msg
            write(*,*) 'q,deltaD=',q,deltaD(xt/q)
            iso_verif_aberrant_enc_choix_nostop=1
            if (abs(xt-q)/q.lt.errmax) then
                write(*,*) 'attention, n''a-t-on pas confondu'// &
     &           ' iso_HDO et iso_eau dans l''appel à la verif?' 
            endif
        endif
        endif

#ifdef ISOVERIF
        if (.not.(iso_HDO.gt.0)) then    
         write(*,*) 'iso205: err_msg=',err_msg,': pourquoi verif?'
          stop
        endif
#else
        write(*,*) 'iso 209: err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   
        
        return
        end function iso_verif_aberrant_enc_choix_nostop

        !*******************

        subroutine iso_verif_aberrant_o17(R17,R18,err_msg)
        implicit none
        ! si l'O17-excess est aberrant, on affiche un message
        ! d'erreur et stoppe

        ! input:
        real R17,R18
        character*(*) err_msg ! message d''erreur à afficher
        !real o17excess

        ! locals
        !integer iso_verif_aberrant_o17_nostop

!        write(*,*) 'O17_verif=',O17_verif
        if (O17_verif) then
            if (iso_verif_aberrant_o17_nostop(R17,R18,err_msg) &
     &           .eq.1) then
                stop
            endif
        endif !if (O17_verif) then

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   
        
        return
        end subroutine iso_verif_aberrant_o17

        !*******************

        function iso_verif_aberrant_o17_nostop(R17,R18,err_msg)
        USE isotopes_mod, ONLY: tnat,iso_O17,iso_O18
        implicit none
        ! si l'O17-excess est aberrant, on affiche un message
        ! d'erreur et renvoit 1

        ! input:
        real R17,R18
        character*(*) err_msg ! message d''erreur à afficher
        !local
        !real o17excess
        ! output
        integer iso_verif_aberrant_o17_nostop

        if (O17_verif) then

        iso_verif_aberrant_o17_nostop=0
        if ((o17excess(R17,R18).gt.o17excess_haut).or. &
     &            (o17excess(R17,R18).lt.o17excess_bas)) then  
            write(*,*) 'erreur detectee par iso_verif_aberrant_O17:'
            write(*,*) err_msg
            write(*,*) 'o17excess=',o17excess(R17,R18)
            write(*,*) 'deltaO17=',(R17/tnat(iso_o17)-1.0)*1000.0
            write(*,*) 'deltaO18=',(R18/tnat(iso_O18)-1.0)*1000.0
            ! attention, vérifier que la ligne suivante est bien activée
            iso_verif_aberrant_o17_nostop=1
        endif

        endif  !if (O17_verif) then

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   
        
        return
        end function iso_verif_aberrant_o17_nostop


        !************************************

        subroutine iso_verif_noNaN(x,err_msg)
        implicit none
        ! si x est NaN, on affiche message
        ! d'erreur et stoppe

        ! input:
        real x
        character*(*) err_msg ! message d''erreur à afficher

        ! locals
        !integer iso_verif_noNAN_nostop

        if (iso_verif_noNAN_nostop(x,err_msg).eq.1) then
            stop
        endif
#ifdef ISOVERIF
#else
        write(*,*) 'err_msg iso443=',err_msg,': pourquoi verif?'
        stop
#endif           
        
        end subroutine iso_verif_noNaN

                !************************************

        function iso_verif_noNaN_nostop(x,err_msg)
        implicit none
        ! si x est NaN, on affiche message
        ! d'erreur et return 1 si erreur

        ! input:
        real x
        character*(*) err_msg ! message d''erreur à afficher

        ! output 
        integer iso_verif_noNAN_nostop

        if ((x.gt.-borne).and.(x.lt.borne)) then
                iso_verif_noNAN_nostop=0
        else
            write(*,*) 'erreur detectee par iso_verif_nonNAN:'
            write(*,*) err_msg
            write(*,*) 'x=',x
            iso_verif_noNAN_nostop=1
        endif      

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg iso 482=',err_msg,': pourquoi verif?'
        stop
#endif           

        return
        end function iso_verif_noNaN_nostop

        subroutine iso_verif_noNaN_vect2D(x,err_msg,ni,n,m)
        implicit none
        ! si x est NaN, on affiche message
        ! d'erreur et return 1 si erreur

        ! input:
          integer n,m,ni
        real x(ni,n,m)
        character*(*) err_msg ! message d''erreur à afficher

        ! output 

        ! locals        
        integer i,j,ixt

      do i=1,n
       do j=1,m
        do ixt=1,ni
         if ((x(ixt,i,j).gt.-borne).and. &
     &            (x(ixt,i,j).lt.borne)) then
         else !if ((x(ixt,i,j).gt.-borne).and.
            write(*,*) 'erreur detectee par iso_verif_nonNAN:'
            write(*,*) err_msg
            write(*,*) 'x,ixt,i,j=',x(ixt,i,j),ixt,i,j
            stop
         endif  !if ((x(ixt,i,j).gt.-borne).and.
        enddo !do ixt=1,ni
       enddo !do j=1,m
      enddo !do i=1,n     

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg iso525=',err_msg,': pourquoi verif?'
        stop
#endif           

        end subroutine iso_verif_noNaN_vect2D


        subroutine iso_verif_noNaN_vect1D(x,err_msg,ni,n)
        implicit none
        ! si x est NaN, on affiche message
        ! d'erreur et return 1 si erreur

        ! input:
          integer n,ni
        real x(ni,n)
        character*(*) err_msg ! message d''erreur à afficher

        ! output 

        ! locals        
        integer i,ixt

      do i=1,n
        do ixt=1,ni
         if ((x(ixt,i).gt.-borne).and. &
     &            (x(ixt,i).lt.borne)) then
         else !if ((x(ixt,i,j).gt.-borne).and.
            write(*,*) 'erreur detectee par iso_verif_nonNAN:'
            write(*,*) err_msg
            write(*,*) 'x,ixt,i=',x(ixt,i),ixt,i
            stop
         endif  !if ((x(ixt,i,j).gt.-borne).and.
        enddo !do ixt=1,ni
      enddo !do i=1,n     

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg iso525=',err_msg,': pourquoi verif?'
        stop
#endif           

        end subroutine iso_verif_noNaN_vect1D



        !************************
        subroutine iso_verif_egalite_choix(a,b,err_msg,erabs,errel)
        implicit none
        ! compare a et b. Si pas egal, on affiche message
        ! d'erreur et stoppe
        ! pour egalite, on verifie erreur absolue et arreur relative

        ! input:
        real a, b
        real erabs,errel !erreur absolue et relative
        character*(*) err_msg ! message d''erreur à afficher

        ! locals
        !integer iso_verif_egalite_choix_nostop

        if (iso_verif_egalite_choix_nostop( &
     &            a,b,err_msg,erabs,errel).eq.1) then
           stop
        endif

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
        stop
#endif           

        end subroutine iso_verif_egalite_choix

!************************
        function iso_verif_egalite_choix_nostop &
     &           (a,b,err_msg,erabs,errel)
        implicit none
        ! compare a et b. Si pas egal, on affiche message
        ! d'erreur et stoppe
        ! pour egalite, on verifie erreur absolue et arreur relative

        ! input:
        real a, b
        real erabs,errel !erreur absolue et relative
        character*(*) err_msg ! message d''erreur à afficher

        ! output
        integer iso_verif_egalite_choix_nostop
        ! locals
        !integer iso_verif_noNaN_nostop

        iso_verif_egalite_choix_nostop=0

#ifdef ISOVERIF
        if (iso_verif_noNaN_nostop(a,err_msg).eq.1) then
            write(*,*) 'a=',a
            iso_verif_egalite_choix_nostop=1
        endif      
        if (iso_verif_noNaN_nostop(b,err_msg).eq.1) then
            write(*,*) 'b=',b
            iso_verif_egalite_choix_nostop=1
        endif      
#endif

        if (abs(a-b).gt.erabs) then
          if (abs((a-b)/max(max(abs(b),abs(a)),1e-18)) &
     &            .gt.errel) then
            write(*,*) 'erreur detectee par iso_verif_egalite:'
            write(*,*) err_msg
            write(*,*) 'a=',a
            write(*,*) 'b=',b
            write(*,*) 'erabs,errel=',erabs,errel
            iso_verif_egalite_choix_nostop=1
!            stop
          endif
        endif

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
        stop
#endif           
        
        return
        end function iso_verif_egalite_choix_nostop       



        !******************
        subroutine iso_verif_positif(x,err_msg)
        use isotopes_mod, ONLY: ridicule
        implicit none
        ! si x<0, on plante.
        ! si très limite, on le met à 0.

        ! input:
        real x
        character*(*) err_msg ! message d''erreur à afficher

        ! locals 
        !integer iso_verif_positif_choix_nostop

        if (iso_verif_positif_choix_nostop(x,ridicule,err_msg) &
     &           .eq.1) then
           stop
        endif
     
#ifdef ISOVERIF
#else
        write(*,*) 'iso_verif 637: err_msg=',err_msg, &
     &           ': pourquoi verif?'
        stop
#endif
        
        end subroutine iso_verif_positif

        !******************
        subroutine iso_verif_positif_vect(x,n,err_msg)
        use isotopes_mod, ONLY: ridicule
        implicit none
        ! si x<0, on plante.

        ! input:
        integer n
        real x(n)
        character*(*) err_msg ! message d''erreur à afficher

        ! locals
        integer i
        integer iso_verif_positif_nostop
        integer ifaux

        iso_verif_positif_nostop=0
        do i=1,n
          if (x(i).lt.-ridicule) then
            iso_verif_positif_nostop=1
            ifaux=i
          endif
        enddo 
        if (iso_verif_positif_nostop.eq.1) then
          write(*,*) 'erreur detectee par iso_verif_positif_vect:'
          write(*,*) err_msg
          write(*,*) 'i,x=',ifaux,x(ifaux)
          stop
        endif   
        
        end subroutine iso_verif_positif_vect

        subroutine iso_verif_positif_choix_vect(x,n,ridic,err_msg)
        implicit none
        ! si x<0, on plante.

        ! input:
        integer n
        real x(n)
        real ridic
        character*(*) err_msg ! message d''erreur à afficher
        ! locals 
        integer i
        integer iso_verif_positif_nostop
        integer ifaux

        iso_verif_positif_nostop=0
        do i=1,n
          if (x(i).lt.-ridic) then
                iso_verif_positif_nostop=1
                ifaux=i
          endif
        enddo 
        if (iso_verif_positif_nostop.eq.1) then                        
         write(*,*) 'erreur detectee par iso_verif_positif_choix_vect:'
         write(*,*) err_msg
         write(*,*) 'i,x=',ifaux,x(ifaux)
         stop
        endif   
        
        end subroutine iso_verif_positif_choix_vect


!******************
        subroutine iso_verif_positif_strict(x,err_msg)
        implicit none
        ! si x<0, on plante.
        ! si très limite, on le met à 0.

        ! input:
        real x
        character*(*) err_msg ! message d''erreur à afficher

        ! locals
        !integer iso_verif_positif_strict_nostop

        if (iso_verif_positif_strict_nostop(x,err_msg) &
     &           .eq.1) then
           stop
        endif            
        
        end subroutine iso_verif_positif_strict

!******************

        function iso_verif_positif_strict_nostop(x,err_msg)
        implicit none
        ! si x<0, on plante.
        ! si très limite, on le met à 0.

        ! input:
        real x
        character*(*) err_msg ! message d''erreur à afficher*

        ! output
        integer iso_verif_positif_strict_nostop

        if (x.gt.0.0) then
            iso_verif_positif_strict_nostop=0
        else     
            write(*,*) 'erreur detectee par iso_verif_positif_strict:'
            write(*,*) err_msg
            write(*,*) 'x=',x  
            iso_verif_positif_strict_nostop=1    
!            stop  
        endif    
        
        return
        end function iso_verif_positif_strict_nostop

!******************

        subroutine iso_verif_positif_choix(x,ridic,err_msg)
        implicit none
        ! si x<0, on plante.
        ! si très limite, on le met à 0.

        ! input:
        real x
        real ridic
        character*(*) err_msg ! message d''erreur à afficher

        ! locals
        !integer iso_verif_positif_choix_nostop

        if (iso_verif_positif_choix_nostop(x,ridic,err_msg) &
     &           .eq.1) then
           stop
        endif
     
#ifdef ISOVERIF
#else
        write(*,*) 'iso_verif 801: err_msg=',err_msg, &
     &           ': pourquoi verif?'
        stop
#endif           
        
        end subroutine iso_verif_positif_choix       

        !******************
        function iso_verif_positif_nostop(x,err_msg)
        use isotopes_mod, ONLY: ridicule
        implicit none
        ! si x<0, on plante.
        ! si très limite, on le met à 0.

        ! input:
        real x
        character*(*) err_msg ! message d''erreur à afficher

        ! output
        integer iso_verif_positif_nostop

        ! locals
        !integer iso_verif_positif_choix_nostop

        iso_verif_positif_nostop=iso_verif_positif_choix_nostop &
     &          (x,ridicule,err_msg)

#ifdef ISOVERIF
#else
        write(*,*) 'iso_verif 837: err_msg=',err_msg, &
     &           ': pourquoi verif?'
        stop

#endif         
        
        return
        end function iso_verif_positif_nostop

        !******************
        function iso_verif_positif_choix_nostop(x,ridic,err_msg)
        implicit none
        ! si x<0, on plante.
        ! si très limite, on le met à 0.

        ! input:
        real x
        real ridic
        character*(*) err_msg ! message d''erreur à afficher

        ! output
        integer iso_verif_positif_choix_nostop

        if (x.lt.-ridic) then
            write(*,*) 'erreur detectee par iso_verif_positif_nostop:'
            write(*,*) err_msg
            write(*,*) 'x=',x
            iso_verif_positif_choix_nostop=1
        else
!            x=max(x,0.0)
            iso_verif_positif_choix_nostop=0
        endif

#ifdef ISOVERIF
#else
        write(*,*) 'iso_verif 877: err_msg=',err_msg, &
     &           ': pourquoi verif?'
        stop
#endif
        
        return
        end function iso_verif_positif_choix_nostop


        !**************

        
        subroutine iso_verif_O18_aberrant(Rd,Ro,err_msg)
        implicit none

        ! vérifie que:
        ! 1) deltaD et deltaO18 dans bonne gamme
        ! 2) dexcess dans bonne gamme
        ! input:
        real Rd,Ro
        character*(*) err_msg ! message d''erreur à afficher

        ! local
        !integer iso_verif_O18_aberrant_nostop

        if (iso_verif_O18_aberrant_nostop(Rd,Ro,err_msg).eq.1) then
            stop
        endif

        end subroutine iso_verif_O18_aberrant
        
        function iso_verif_O18_aberrant_nostop(Rd,Ro,err_msg)
        USE isotopes_mod, ONLY: tnat, iso_HDO, iso_O18
        implicit none

        ! vérifie que:
        ! 1) deltaD et deltaO18 dans bonne gamme
        ! 2) dexcess dans bonne gamme

        ! input:
        real Rd,Ro
        character*(*) err_msg ! message d''erreur à afficher

        ! outputs
        integer iso_verif_O18_aberrant_nostop

        !locals
        real deltaD,deltao,dexcess

        deltaD=(Rd/tnat(iso_hdo)-1)*1000
        deltao=(Ro/tnat(iso_O18)-1)*1000
        dexcess=deltaD-8*deltao

        iso_verif_O18_aberrant_nostop=0
        if ((deltaD.lt.deltaDmin).or.(deltao.lt.deltaDmin/2.0).or. &
     &        (deltaD.gt.deltalim).or.(deltao.gt.deltalim/8.0).or. &
     &        ((deltaD.gt.-500.0).and.((dexcess.lt.dexcess_min) &
     &        .or.(dexcess.gt.dexcess_max)))) then 
            write(*,*) 'erreur detectee par iso_verif_O18_aberrant:'
            write(*,*) err_msg
            write(*,*) 'delta180=',deltao
            write(*,*) 'deltaD=',deltaD
            write(*,*) 'Dexcess=',dexcess
            write(*,*) 'tnat=',tnat
!            stop
            iso_verif_O18_aberrant_nostop=1
          endif

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
        stop
#endif                   

        return
        end function iso_verif_O18_aberrant_nostop


        ! **********
        function deltaD(R)
        USE isotopes_mod, ONLY: tnat,iso_HDO
        implicit none
        real R,deltaD

        
        if (iso_HDO.gt.0) then
           deltaD=(R/tnat(iso_HDO)-1)*1000.0
        else
            write(*,*) 'iso_verif_egalite>deltaD 260: iso_HDO.gt.0=', &
     &           iso_HDO.gt.0
        endif
        return
        end function deltaD

        ! **********
        function deltaO(R)
        USE isotopes_mod, ONLY: tnat,iso_O18
        implicit none
        real R,deltaO
        
        if (iso_O18.gt.0) then
           deltaO=(R/tnat(iso_O18)-1)*1000.0
        else
            write(*,*) 'iso_verif_egalite>deltaO18 260: iso_O18.gt.0=', &
     &           iso_O18.gt.0
        endif
        return
        end function deltaO

        ! **********
        function dexcess(RD,RO)
        USE isotopes_mod, ONLY: tnat,iso_O18,iso_HDO
        implicit none
        real RD,RO,deltaD,deltaO,dexcess
        
        if ((iso_O18.gt.0).and.(iso_HDO.gt.0)) then
           deltaD=(RD/tnat(iso_HDO)-1)*1000.0
           deltaO=(RO/tnat(iso_O18)-1)*1000.0
           dexcess=deltaD-8*deltaO
        else
            write(*,*) 'iso_verif_egalite 1109: iso_O18,iso_HDO=',iso_O18,iso_HDO
        endif
        return
        end function dexcess


        ! **********
        function delta_all(R,ixt)
        USE isotopes_mod, ONLY: tnat
        implicit none
        real R,delta_all
        integer ixt 
        
        delta_all=(R/tnat(ixt)-1)*1000.0
        return
        end function delta_all

        ! **********
        function delta_to_R(delta,ixt)
        USE isotopes_mod, ONLY: tnat
        implicit none
        real delta,delta_to_R
        integer ixt
        
        delta_to_R=(delta/1000.0+1.0)*tnat(ixt)
        return
        end function delta_to_R

         ! **********
        function o17excess(R17,R18)
        USE isotopes_mod, ONLY: tnat,iso_O18,iso_O17
        implicit none
        real R17,R18,o17excess
        
        if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
            
           o17excess=1e6*(log(R17/tnat(iso_o17)) &
     &           -0.528*log(R18/tnat(iso_O18)))
!           write(*,*) 'o17excess=',o17excess
        else
            write(*,*) 'iso_verif_egalite>deltaD 260: iso_O17.gt.0,18=', &
     &           iso_O17.gt.0,iso_O18.gt.0
        endif
        return
        end function o17excess

        !       ****************

          subroutine iso_verif_egalite_vect2D( &
     &           xt,q,err_msg,ni,n,m)
        
        USE isotopes_mod, ONLY: iso_eau
          implicit none

          ! inputs 
          integer n,m,ni
          real q(n,m)
          real xt(ni,n,m)
          character*(*) err_msg

        ! locals
        integer iso_verif_egalite_nostop_loc
        integer i,j,ixt
        integer ifaux,jfaux

        !write(*,*) 'iso_verif_egalite_vect2D 1099 tmp: q(2,1),xt(iso_eau,2,1)=',q(2,1),xt(iso_eau,2,1)
        !write(*,*) 'ni,n,m=',ni,n,m,errmax,errmaxrel
        if (iso_eau.gt.0) then
        iso_verif_egalite_nostop_loc=0
        do i=1,n
         do j=1,m
          if (abs(q(i,j)-xt(iso_eau,i,j)).gt.errmax) then
           if (abs((q(i,j)-xt(iso_eau,i,j))/ &
     &           max(max(abs(q(i,j)),abs(xt(iso_eau,i,j))),1e-18)) &
     &           .gt.errmaxrel) then
              iso_verif_egalite_nostop_loc=1
              ifaux=i
              jfaux=j
           endif
          endif 
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_egalite_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'xt,q=',xt(iso_eau,ifaux,jfaux),q(ifaux,jfaux)
          stop
        endif 
        endif
        
#ifdef ISOVERIF
        call iso_verif_noNaN_vect2D(xt,err_msg,ni,n,m)
#endif         

        return
        end subroutine iso_verif_egalite_vect2D

        subroutine iso_verif_egalite_vect1D( &
     &           xt,q,err_msg,ni,n)

        USE isotopes_mod, ONLY: iso_eau
        implicit none

        ! inputs 
        integer n,ni
        real q(n)
        real xt(ni,n)
        character*(*) err_msg

        ! locals
        integer iso_verif_egalite_nostop_loc
        integer i
        integer ifaux

        if (iso_eau.gt.0) then
        iso_verif_egalite_nostop_loc=0
        do i=1,n
          if (abs(q(i)-xt(iso_eau,i)).gt.errmax) then
           if (abs((q(i)-xt(iso_eau,i))/ &
     &           max(max(abs(q(i)),abs(xt(iso_eau,i))),1e-18)) &
     &           .gt.errmaxrel) then
              iso_verif_egalite_nostop_loc=1
              ifaux=i
           endif !if (abs((q(i)-xt(iso_eau,i))/
          endif !if (abs(q(i)-xt(iso_eau,i)).gt.errmax) then
        enddo !do i=1,n

        if (iso_verif_egalite_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:'
          write(*,*) err_msg
          write(*,*) 'i=',ifaux
          write(*,*) 'xt,q=',xt(iso_eau,ifaux),q(ifaux)
          stop
        endif  !if (iso_verif_egalite_nostop.eq.1) then
        endif !if (iso_eau.gt.0) then

        end subroutine iso_verif_egalite_vect1D       

        subroutine iso_verif_egalite_std_vect( &
     &           a,b,err_msg,n,m,errmax,errmaxrel)

          implicit none

          ! inputs 
          integer n,m,ni
          real a(n,m)
          real b(n,m)
          character*(*) err_msg
          real errmax,errmaxrel

        ! locals
        integer iso_verif_egalite_nostop_loc
        integer i,j
        integer ifaux,jfaux

        iso_verif_egalite_nostop_loc=0
        do i=1,n
         do j=1,m
          if (abs(a(i,j)-b(i,j)).gt.errmax) then
           if (abs((a(i,j)-b(i,j))/ &
     &           max(max(abs(a(i,j)),abs(b(i,j))),1e-18)) &
     &           .gt.errmaxrel) then
              iso_verif_egalite_nostop_loc=1
              ifaux=i
              jfaux=j
           endif
          endif
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_egalite_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'a,b=',a(ifaux,jfaux),b(ifaux,jfaux)
          stop
        endif 

        return
        end subroutine iso_verif_egalite_std_vect

        subroutine iso_verif_aberrant_vect2D( &
     &           xt,q,err_msg,ni,n,m)
        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
          implicit none

          ! inputs   
          integer n,m,ni
          real q(n,m)
          real xt(ni,n,m)
          character*(*) err_msg

        ! locals
        integer iso_verif_aberrant_nostop_loc
        integer i,j
        integer ifaux,jfaux
        !real deltaD

        if (iso_HDO.gt.0) then
        iso_verif_aberrant_nostop_loc=0
        do i=1,n
         do j=1,m
          if (q(i,j).gt.ridicule) then 
            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .gt.deltalim).or. &
     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .lt.-borne)) then    
              iso_verif_aberrant_nostop_loc=1
              ifaux=i
              jfaux=j
           endif
          endif
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_aberrant_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par iso_verif_aberrant_vect2D:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
     &           /q(ifaux,jfaux))
          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
          stop
        endif  
        endif !if (iso_HDO.gt.0) then

        end subroutine iso_verif_aberrant_vect2D       

        subroutine iso_verif_aberrant_enc_vect2D( &
     &           xt,q,err_msg,ni,n,m)

        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
          implicit none

          ! inputs   
          integer n,m,ni
          real q(n,m)
          real xt(ni,n,m)
          character*(*) err_msg

        ! locals
        integer iso_verif_aberrant_nostop_loc
        integer i,j
        integer ifaux,jfaux
        !real deltaD

        if (iso_HDO.gt.0) then
        iso_verif_aberrant_nostop_loc=0
        do i=1,n
         do j=1,m
          if (q(i,j).gt.ridicule) then 
            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .gt.deltalim).or. &
     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .lt.deltaDmin).or. &
     &           (xt(iso_HDO,i,j).lt.-borne).or. &
     &           (xt(iso_HDO,i,j).gt.borne)) then     
              iso_verif_aberrant_nostop_loc=1
              ifaux=i
              jfaux=j
           endif
          endif
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_aberrant_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par ', &
     &           'iso_verif_aberrant_enc_vect2D:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
     &           /q(ifaux,jfaux))
          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
          write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux)
          call abort_physic('isotopes_verif_mod','iso_verif_aberrant_enc_vect2D',1)
        endif  
        endif !if (iso_HDO.gt.0) then

        end subroutine iso_verif_aberrant_enc_vect2D       

        
        subroutine iso_verif_aberrant_enc_vect2D_ns( &
     &           xt,q,err_msg,ni,n,m)

        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
          implicit none

          ! inputs   
          integer n,m,ni
          real q(n,m)
          real xt(ni,n,m)
          character*(*) err_msg

        ! locals
        integer iso_verif_aberrant_nostop_loc
        integer i,j
        integer ifaux,jfaux
        !real deltaD

        if (iso_HDO.gt.0) then
        iso_verif_aberrant_nostop_loc=0
        do i=1,n
         do j=1,m
          if (q(i,j).gt.ridicule) then 
            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .gt.deltalim).or. &
     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .lt.deltaDmin)) then    
              iso_verif_aberrant_nostop_loc=1
              ifaux=i
              jfaux=j
           endif
          endif
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_aberrant_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par ', &
     &           'iso_verif_aberrant_vect2D_ns:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
     &           /q(ifaux,jfaux))
          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
!          stop
        endif  
        endif !if (iso_HDO.gt.0) then

        end subroutine iso_verif_aberrant_enc_vect2D_ns       


         subroutine iso_verif_aberrant_vect2Dch( &
     &           xt,q,err_msg,ni,n,m,deltaDmax)

        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
          implicit none


          ! inputs   
          integer n,m,ni
          real q(n,m)
          real xt(ni,n,m)
          character*(*) err_msg
          real deltaDmax

        ! locals
        integer iso_verif_aberrant_nostop_loc
        integer i,j
        integer ifaux,jfaux
        !real deltaD

        if (iso_HDO.gt.0) then
        iso_verif_aberrant_nostop_loc=0
        do i=1,n
         do j=1,m
          if (q(i,j).gt.ridicule) then 
            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .gt.deltaDmax).or. &
     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .lt.-borne)) then    
              iso_verif_aberrant_nostop_loc=1
              ifaux=i
              jfaux=j
           endif
          endif
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_aberrant_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par iso_verif_aberrant_vect2D:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
     &           /q(ifaux,jfaux))
          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
          stop
        endif  
        endif !if (iso_HDO.gt.0) then

        end subroutine iso_verif_aberrant_vect2Dch     

        subroutine iso_verif_O18_aberrant_enc_vect2D( &
     &           xt,q,err_msg,ni,n,m)

        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO,iso_O18
          implicit none

          ! inputs   
          integer n,m,ni
          real q(n,m)
          real xt(ni,n,m)
          character*(*) err_msg

        ! locals
        integer iso_verif_aberrant_nostop_loc
        integer i,j
        integer ifaux,jfaux
        real deltaDloc,deltaoloc,dexcessloc

        if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
        iso_verif_aberrant_nostop_loc=0
        do i=1,n
         do j=1,m
          if (q(i,j).gt.ridicule) then 

            deltaDloc=(xt(iso_HDO,i,j)/q(i,j)/tnat(iso_hdo)-1)*1000
            deltaoloc=(xt(iso_O18,i,j)/q(i,j)/tnat(iso_O18)-1)*1000
            dexcessloc=deltaDloc-8*deltaoloc
            if ((deltaDloc.lt.deltaDmin).or.(deltaoloc.lt.deltaDmin/2.0).or. &
     &        (deltaDloc.gt.deltalim).or.(deltaoloc.gt.deltalim/8.0).or. &
     &        ((deltaDloc.gt.-500.0).and.((dexcessloc.lt.dexcess_min) &
     &        .or.(dexcessloc.gt.dexcess_max)))) then
              iso_verif_aberrant_nostop_loc=1
              ifaux=i
              jfaux=j
              write(*,*) 'deltaD,deltao,dexcess=',deltaDloc,deltaoloc,dexcessloc
           endif
          endif
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_aberrant_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par ', &
     &           'iso_verif_aberrant_enc_vect2D:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
          write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux)
          call abort_physic('isotopes_verif_mod','iso_verif_aberrant_enc_vect2D',1)
        endif  
        endif !if (iso_HDO.gt.0) then

        end subroutine iso_verif_O18_aberrant_enc_vect2D    


        subroutine select_dim23_from4D(n1,n2,n3,n4, &
     &          var,vec,i1,i4)
        implicit none

        ! inputs
        integer n1,n2,n3,n4
        real var(n1,n2,n3,n4)
        integer i1,i4
        ! outputs
        real vec(n2,n3)
        ! locals
        integer i2,i3

        do i2=1,n2
         do i3=1,n3
          vec(i2,i3)=var(i1,i2,i3,i4)
         enddo
        enddo

        return
        end subroutine select_dim23_from4D

        
        subroutine select_dim4_from4D(ntime,nlev,nlat,nlon, &
     &          var,vec,itime,ilev,ilat)
        implicit none

        ! inputs
        integer ntime,nlev,nlat,nlon
        real var(ntime,nlev,nlat,nlon)
        integer itime,ilev,ilat
        ! outputs
        real vec(nlon)
        ! locals
        integer ilon

        do ilon=1,nlon
          vec(ilon)=var(itime,ilev,ilat,ilon)
        enddo

        return
        end subroutine select_dim4_from4D

        subroutine select_dim5_from5D(n1,n2,n3,n4,n5, &
     &          var,vec,i1,i2,i3,i4)
        implicit none

        ! inputs
        integer n1,n2,n3,n4,n5
        real var(n1,n2,n3,n4,n5)
        integer i1,i2,i3,i4
        ! outputs
        real vec(n5)
        ! locals
        integer i5

        do i5=1,n5
          vec(i5)=var(i1,i2,i3,i4,i5)
        enddo

        end subroutine select_dim5_from5D

        
        subroutine select_dim3_from3D(ntime,nlat,nlon, &
     &          var,vec,itime,ilat)
        implicit none

        ! inputs
        integer ntime,nlat,nlon
        real var(ntime,nlat,nlon)
        integer itime,ilat
        ! outputs
        real vec(nlon)
        ! locals
        integer ilon

        do ilon=1,nlon
          vec(ilon)=var(itime,ilat,ilon)
        enddo

        end subroutine select_dim3_from3D

        
        subroutine select_dim23_from3D(n1,n2,n3, &
     &          var,vec,i1)
        implicit none

        ! inputs
        integer n1,n2,n3
        real var(n1,n2,n3)
        integer i1
        ! outputs
        real vec(n2,n3)
        ! locals
        integer i2,i3

        do i2=1,n2
         do i3=1,n3
          vec(i2,i3)=var(i1,i2,i3)
         enddo
        enddo

        end subroutine select_dim23_from3D

        subroutine putinto_dim23_from4D(n1,n2,n3,n4, &
     &          var,vec,i1,i4)
        implicit none

        ! inputs
        integer n1,n2,n3,n4
        real vec(n2,n3)
        integer i1,i4
        ! inout
        real var(n1,n2,n3,n4)
        ! locals
        integer i2,i3

       do i2=1,n2
        do i3=1,n3
          var(i1,i2,i3,i4)=vec(i2,i3)
         enddo
        enddo

        end subroutine putinto_dim23_from4D

        
        subroutine putinto_dim12_from4D(n1,n2,n3,n4, &
     &          var,vec,i3,i4)
        implicit none

        ! inputs
        integer n1,n2,n3,n4
        real vec(n1,n2)
        integer i3,i4
        ! inout
        real var(n1,n2,n3,n4)
        ! locals
        integer i1,i2

       do i1=1,n1
        do i2=1,n2
          var(i1,i2,i3,i4)=vec(i1,i2)
         enddo
        enddo

        end subroutine putinto_dim12_from4D
        
        subroutine putinto_dim23_from3D(n1,n2,n3, &
     &          var,vec,i1)
        implicit none

        ! inputs
        integer n1,n2,n3
        real vec(n2,n3)
        integer i1
        ! inout
        real var(n1,n2,n3)
        ! locals
        integer i2,i3

       do i2=1,n2
        do i3=1,n3
          var(i1,i2,i3)=vec(i2,i3)
         enddo
        enddo

        end subroutine putinto_dim23_from3D

        

        subroutine iso_verif_noNaN_par2D(x,err_msg,ni,n,m,ib,ie)
        implicit none
        ! si x est NaN, on affiche message
        ! d'erreur et return 1 si erreur

        ! input:
          integer n,m,ni,ib,ie
        real x(ni,n,m)
        character*(*) err_msg ! message d''erreur à afficher

        ! output 

        ! locals        
        integer i,j,ixt

      do i=ib,ie
       do j=1,m
        do ixt=1,ni
         if ((x(ixt,i,j).gt.-borne).and. &
     &            (x(ixt,i,j).lt.borne)) then
         else !if ((x(ixt,i,j).gt.-borne).and.
            write(*,*) 'erreur detectee par iso_verif_nonNAN:'
            write(*,*) err_msg
            write(*,*) 'x,ixt,i,j=',x(ixt,i,j),ixt,i,j
            stop
         endif  !if ((x(ixt,i,j).gt.-borne).and.
        enddo !do ixt=1,ni
       enddo !do j=1,m
      enddo !do i=1,n     

#ifdef ISOVERIF
#else
        write(*,*) 'err_msg iso1772=',err_msg,': pourquoi verif?'
        stop
#endif           

        return
        end subroutine iso_verif_noNaN_par2D

        
        subroutine iso_verif_aberrant_enc_par2D( &
     &           xt,q,err_msg,ni,n,m,ib,ie)

        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
          implicit none

          ! inputs   
          integer n,m,ni,ib,ie
          real q(n,m)
          real xt(ni,n,m)
          character*(*) err_msg

        ! locals
        integer iso_verif_aberrant_nostop_loc
        integer i,j
        integer ifaux,jfaux
        !real deltaD

        if (iso_HDO.gt.0) then
        iso_verif_aberrant_nostop_loc=0
        do i=ib,ie
         do j=1,m
          if (q(i,j).gt.ridicule) then 
            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .gt.deltalim).or. &
     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
     &                   .lt.deltaDmin)) then    
              iso_verif_aberrant_nostop_loc=1
              ifaux=i
              jfaux=j
           endif
          endif
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_aberrant_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par iso_verif_aberrant_par2D:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
     &           /q(ifaux,jfaux))
          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
          write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux)
          stop
        endif  
        endif !if (iso_HDO.gt.0) then

        end subroutine iso_verif_aberrant_enc_par2D       

        
          subroutine iso_verif_egalite_par2D( &
     &           xt,q,err_msg,ni,n,m,ib,ie)
        
        USE isotopes_mod, ONLY: iso_eau
          implicit none

          ! inputs 
          integer n,m,ni,ib,ie
          real q(n,m)
          real xt(ni,n,m)
          character*(*) err_msg

        ! locals
        integer iso_verif_egalite_nostop_loc
        integer i,j
        integer ifaux,jfaux

        if (iso_eau.gt.0) then
        iso_verif_egalite_nostop_loc=0
        do i=ib,ie
         do j=1,m
          if (abs(q(i,j)-xt(iso_eau,i,j)).gt.errmax) then
           if (abs((q(i,j)-xt(iso_eau,i,j))/ &
     &           max(max(abs(q(i,j)),abs(xt(iso_eau,i,j))),1e-18)) &
     &           .gt.errmaxrel) then
              iso_verif_egalite_nostop_loc=1
              ifaux=i
              jfaux=j
           endif
          endif
         enddo !do j=1,m
        enddo !do i=1,n

        if (iso_verif_egalite_nostop_loc.eq.1) then
          write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:'
          write(*,*) err_msg
          write(*,*) 'i,j=',ifaux,jfaux
          write(*,*) 'xt,q=',xt(iso_eau,ifaux,jfaux),q(ifaux,jfaux)
          stop
        endif 
        endif

        end subroutine iso_verif_egalite_par2D

#ifdef ISOTRAC

      function iso_verif_traceur_choix_nostop(x,err_msg, &
     &       errmax,errmaxrel,ridicule_trac,deltalimtrac)
        use isotopes_mod, ONLY: iso_HDO
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher
       real errmax,errmaxrel,ridicule_trac,deltalimtrac

       ! output
       integer iso_verif_traceur_choix_nostop

       ! locals
       !integer iso_verif_traceur_noNaN_nostop 
       !integer iso_verif_tracm_choix_nostop
       !integer iso_verif_tracdD_choix_nostop
       !integer iso_verif_tracpos_choix_nostop

        iso_verif_traceur_choix_nostop=0  

        ! verif noNaN
        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
             iso_verif_traceur_choix_nostop=1
        endif 
        
        ! verif masse
        if (iso_verif_tracm_choix_nostop(x,err_msg, &
     &           errmax,errmaxrel).eq.1) then
             iso_verif_traceur_choix_nostop=1
        endif              

        ! verif deltaD
        if (iso_HDO.gt.0) then
        if (iso_verif_tracdD_choix_nostop(x,err_msg, &
     &           ridicule_trac,deltalimtrac).eq.1) then
             iso_verif_traceur_choix_nostop=1
        endif  
        endif !if (iso_HDO.gt.0) then  

        ! verif pas aberramment negatif
        if (iso_verif_tracpos_choix_nostop(x,err_msg, &
     &           1e-5).eq.1) then
             iso_verif_traceur_choix_nostop=1
        endif 

        end function iso_verif_traceur_choix_nostop

        function iso_verif_tracnps_choix_nostop(x,err_msg, &
     &       errmax,errmaxrel,ridicule_trac,deltalimtrac)
        USE isotopes_mod, ONLY: iso_HDO
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant
        ! on ne vérfie pas la positivité
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher
       real errmax,errmaxrel,ridicule_trac,deltalimtrac

       ! output
       integer iso_verif_tracnps_choix_nostop

       ! locals
       !integer iso_verif_traceur_noNaN_nostop 
       !integer iso_verif_tracm_choix_nostop
       !integer iso_verif_tracdD_choix_nostop

        iso_verif_tracnps_choix_nostop=0  

        ! verif noNaN
        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
             iso_verif_tracnps_choix_nostop=1
        endif 
        
        ! verif masse
        if (iso_verif_tracm_choix_nostop(x,err_msg, &
     &           errmax,errmaxrel).eq.1) then
             iso_verif_tracnps_choix_nostop=1
        endif              

        ! verif deltaD
        if (iso_HDO.gt.0) then
        if (iso_verif_tracdD_choix_nostop(x,err_msg, &
     &           ridicule_trac,deltalimtrac).eq.1) then
             iso_verif_tracnps_choix_nostop=1
        endif   
        endif ! if (iso_HDO.gt.0) then  

        return
        end function iso_verif_tracnps_choix_nostop

        function iso_verif_tracpos_choix_nostop(x,err_msg,seuil)
        use isotopes_mod, only: isoName
        implicit none

        ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher
       real seuil

       ! output
       integer iso_verif_tracpos_choix_nostop

       ! locals
       integer lnblnk
       integer iiso,ixt
       !integer iso_verif_positif_choix_nostop

       iso_verif_tracpos_choix_nostop=0

       do ixt=niso+1,ntraciso
          if (iso_verif_positif_choix_nostop(x(ixt),seuil,err_msg// &
     &           ', verif positif, iso'//TRIM(isoName(ixt))).eq.1) then
            iso_verif_tracpos_choix_nostop=1
          endif
        enddo

        end function iso_verif_tracpos_choix_nostop


        function iso_verif_traceur_noNaN_nostop(x,err_msg)
        use isotopes_mod, only: isoName
        implicit none

        ! on vérifie juste que pas NaN
        ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher

       ! output
       integer iso_verif_traceur_noNaN_nostop

       ! locals
       integer lnblnk
       integer iiso,ixt
       !integer iso_verif_nonaN_nostop

       iso_verif_traceur_noNaN_nostop=0

        do ixt=niso+1,ntraciso
!          write(*,*) 'iso_verif_traceurs 154: iiso,ixt=',iiso,ixt
          if (iso_verif_noNaN_nostop(x(ixt),err_msg// &
     &           ', verif trac no NaN, iso'//TRIM(isoName(ixt))) &
     &           .eq.1) then
            iso_verif_traceur_noNaN_nostop=1
          endif
        enddo

        end function iso_verif_traceur_noNaN_nostop

        function iso_verif_tracm_choix_nostop(x,err_msg, &
     &           errmaxin,errmaxrelin)

        use isotopes_mod, ONLY: ridicule,isoName
        ! on vérifie juste bilan de masse
        implicit none
        
        ! inputs
        real x(ntraciso)
        character*(*) err_msg ! message d''erreur à afficher
        real errmaxin,errmaxrelin

        ! output
        integer iso_verif_tracm_choix_nostop

       ! locals
       !integer iso_verif_egalite_choix_nostop
       integer iiso,izone,ixt
       real xtractot

       iso_verif_tracm_choix_nostop=0

        do iiso=1,niso

          xtractot=0.0
          do izone=1,nzone  
            ixt=itZonIso(izone,iiso) 
            xtractot=xtractot+x(ixt)
          enddo

          if (iso_verif_egalite_choix_nostop(xtractot,x(iiso), &
     &        err_msg//', verif trac egalite1, iso '// &
     &        TRIM(isoName(iiso)), &
     &        errmaxin,errmaxrelin).eq.1) then
            write(*,*) 'iso_verif_traceur 202: x=',x
!            write(*,*) 'xtractot=',xtractot
            do izone=1,nzone  
              ixt=itZonIso(izone,iiso)
              write(*,*) 'izone,iiso,ixt=',izone,iiso,ixt
            enddo
            iso_verif_tracm_choix_nostop=1
          endif

          ! verif ajoutée le 19 fev 2011
          if ((abs(xtractot).lt.ridicule**2).and. &
     &           (abs(x(iiso)).gt.ridicule)) then
            write(*,*) err_msg,', verif masse traceurs, iso ', &
     &          TRIM(isoName(iiso))
            write(*,*) 'iso_verif_traceur 209: x=',x
!            iso_verif_tracm_choix_nostop=1
          endif

        enddo !do iiso=1,ntraceurs_iso  

        end function iso_verif_tracm_choix_nostop

        function iso_verif_tracdD_choix_nostop(x,err_msg, &
     &           ridicule_trac,deltalimtrac)

        USE isotopes_mod, ONLY: iso_eau, iso_HDO
        use isotrac_mod, only: strtrac
        ! on vérifie juste deltaD
        implicit none
                
        ! inputs
        real x(ntraciso)
        character*(*) err_msg ! message d''erreur à afficher
        real ridicule_trac,deltalimtrac

        ! output
        integer iso_verif_tracdD_choix_nostop        

       ! locals
       integer izone,ieau,ixt
       !integer iso_verif_aberrant_choix_nostop

        iso_verif_tracdD_choix_nostop=0

        if ((iso_eau.gt.0).and.(iso_HDO.gt.0)) then
        do izone=1,nzone
             ieau=itZonIso(izone,iso_eau)
             ixt=itZonIso(izone,iso_HDO)

             if (iso_verif_aberrant_choix_nostop(x(ixt),x(ieau), &
     &           ridicule_trac,deltalimtrac,err_msg// &
     &           ', verif trac no aberrant zone '//strtrac(izone)) &
     &           .eq.1) then
               iso_verif_tracdD_choix_nostop=1
             endif
!             if (x(ieau).gt.ridicule) then
!              call iso_verif_aberrant(x(ixt)/x(ieau),
!     :           err_msg//', verif trac no aberrant zone '
!     :           //strtrac(izone))
!             endif
        enddo !do izone=1,nzone
       endif ! if ((iso_eau.gt.0).and.(iso_HDO.gt.0)) then

       end function iso_verif_tracdD_choix_nostop

INTEGER FUNCTION iso_verif_tag17_q_deltaD_chns(x,err_msg) RESULT(res)
  USE isotopes_mod, ONLY: iso_HDO, iso_eau, ridicule
  USE isotrac_mod,  ONLY: nzone_temp, option_traceurs
  IMPLICIT NONE
  REAL,             INTENT(IN) :: x(ntraciso)
  CHARACTER(LEN=*), INTENT(IN) :: err_msg
  INTEGER :: ieau, ixt, ieau1
  res = 0
  IF(ALL([17,18]/=option_traceurs)) RETURN
  !--- Check whether * deltaD(highest tagging layer) < 200 permil
  !                  * q <
  ieau=itZonIso(nzone_temp,iso_eau)
  ixt=itZonIso(nzone_temp,iso_HDO)
  IF(x(ieau)>ridicule) THEN
    IF(iso_verif_positif_nostop(-200.0-deltaD(x(ixt)/x(ieau)), err_msg//': deltaDt05 trop fort')==1) THEN
      res=1; write(*,*) 'x=',x
    END IF
  END IF
  IF(iso_verif_positif_nostop(2.0e-3-x(ieau),err_msg//': qt05 trop fort')==1) THEN
    res=1; write(*,*) 'x=',x
  END IF
  !--- Check whether q is small ; then, qt01 < 10%
  IF(x(iso_eau)<2.0e-3) THEN
    ieau1= itZonIso(1,iso_eau)
    IF(iso_verif_positif_nostop(0.1-(x(ieau1)/x(iso_eau)),err_msg//': qt01 trop abondant')==1) THEN
      res=1; write(*,*) 'x=',x
    END IF
  END IF
END FUNCTION iso_verif_tag17_q_deltaD_chns

SUBROUTINE iso_verif_trac17_q_deltaD(x,err_msg)
  USE isotrac_mod,  ONLY: nzone_temp, option_traceurs
  IMPLICIT NONE
  REAL,             INTENT(IN) :: x(ntraciso)
  CHARACTER(LEN=*), INTENT(IN) :: err_msg
  IF(ALL([17,18]/=option_traceurs)) RETURN
  IF(nzone_temp>=5) THEN
    IF(iso_verif_tag17_q_deltaD_chns(x,err_msg)==1) STOP
  END IF
END SUBROUTINE iso_verif_trac17_q_deltaD

      subroutine iso_verif_traceur(x,err_msg)
        use isotrac_mod, only: ridicule_trac
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher

       ! locals
       !integer iso_verif_traceur_choix_nostop  

        if (iso_verif_traceur_choix_nostop(x,err_msg, &
     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
     &       .eq.1) then
                stop
        endif

        end subroutine iso_verif_traceur

        
      subroutine iso_verif_traceur_retourne3D(x,n1,n2,n3, &
     &           i1,i2,i3,err_msg)
        use isotrac_mod, only: ridicule_trac

        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       integer n1,n2,n3
       real x(n1,n2,n3,ntraciso)
       character*(*) err_msg ! message d''erreur à afficher
       integer i1,i2,i3

       ! locals
       !integer iso_verif_traceur_choix_nostop  
       real xiso(ntraciso)

        call select_dim4_from4D(n1,n2,n3,ntraciso, &
     &      x,xiso,i1,i2,i3)
        if (iso_verif_traceur_choix_nostop(xiso,err_msg, &
     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
     &       .eq.1) then
                stop
        endif

        end subroutine iso_verif_traceur_retourne3D

        subroutine iso_verif_traceur_retourne4D(x,n1,n2,n3,n4, &
     &           i1,i2,i3,i4,err_msg)
        use isotrac_mod, only: ridicule_trac

        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       integer n1,n2,n3,n4
       real x(n1,n2,n3,n4,ntraciso)
       character*(*) err_msg ! message d''erreur à afficher
       integer i1,i2,i3,i4

       ! locals
       !integer iso_verif_traceur_choix_nostop  
       real xiso(ntraciso)

        call select_dim5_from5D(n1,n2,n3,n4,ntraciso, &
     &      x,xiso,i1,i2,i3,i4)
        if (iso_verif_traceur_choix_nostop(xiso,err_msg, &
     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
     &       .eq.1) then
                stop
        endif

        end subroutine iso_verif_traceur_retourne4D

        
      subroutine iso_verif_traceur_retourne2D(x,n1,n2, &
     &           i1,i2,err_msg)
        use isotrac_mod, only: ridicule_trac
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       integer n1,n2
       real x(n1,n2,ntraciso)
       character*(*) err_msg ! message d''erreur à afficher
       integer i1,i2

       ! locals
       !integer iso_verif_traceur_choix_nostop  
       real xiso(ntraciso)

        call select_dim3_from3D(n1,n2,ntraciso, &
     &      x,xiso,i1,i2)
        if (iso_verif_traceur_choix_nostop(xiso,err_msg, &
     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
     &       .eq.1) then
                stop
        endif

        end subroutine iso_verif_traceur_retourne2D

        subroutine iso_verif_traceur_vect(x,n,m,err_msg)
        USE isotopes_mod, ONLY: iso_HDO
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       integer n,m
       real x(ntraciso,n,m)
       character*(*) err_msg ! message d''erreur à afficher

       ! locals
       logical iso_verif_traceur_nostop
       integer i,j,ixt,iiso,izone,ieau
       integer ifaux,jfaux,ixtfaux
       
       call iso_verif_traceur_noNaN_vect(x,n,m,err_msg)        

        ! verif masse: iso_verif_tracm_choix_nostop
        call iso_verif_trac_masse_vect(x,n,m,err_msg,errmax,errmaxrel)

        ! verif deltaD: iso_verif_tracdD_choix_nostop    
        if (iso_HDO.gt.0) then  
        call iso_verif_tracdd_vect(x,n,m,err_msg)      
        endif !if (iso_HDO.gt.0) then        

        ! verif pas aberramment negatif: iso_verif_tracpos_choix_nostop
        call iso_verif_tracpos_vect(x,n,m,err_msg,1e-5)  
        
        end subroutine iso_verif_traceur_vect

        subroutine iso_verif_tracnps_vect(x,n,m,err_msg)
        USE isotopes_mod, ONLY: iso_HDO
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       integer n,m
       real x(ntraciso,n,m)
       character*(*) err_msg ! message d''erreur à afficher

       ! locals
       logical iso_verif_traceur_nostop
       integer i,j,ixt,iiso,izone,ieau
       integer ifaux,jfaux,ixtfaux
       
       call iso_verif_traceur_noNaN_vect(x,n,m,err_msg)        

        ! verif masse: iso_verif_tracm_choix_nostop
        call iso_verif_trac_masse_vect(x,n,m,err_msg,errmax,errmaxrel)

        ! verif deltaD: iso_verif_tracdD_choix_nostop    
        if (iso_HDO.gt.0) then  
        call iso_verif_tracdd_vect(x,n,m,err_msg)      
        endif !if (iso_HDO.gt.0) then        
        
        end subroutine iso_verif_tracnps_vect


        subroutine iso_verif_traceur_noNaN_vect(x,n,m,err_msg)
        implicit none
       
       ! inputs
       integer n,m
       real x(ntraciso,n,m)
       character*(*) err_msg ! message d''erreur à afficher

       ! locals
       logical iso_verif_traceur_nostop
       integer i,j,ixt,iiso
       integer ifaux,jfaux,ixtfaux


       iso_verif_traceur_nostop=.false.
        ! verif noNaN: iso_verif_traceur_noNaN_nostop        
        do j=1,m
        do i=1,n
        do ixt=niso+1,ntraciso
          if ((x(ixt,i,j).gt.-borne).and.(x(ixt,i,j).lt.borne)) then
          else !if ((x.gt.-borne).and.(x.lt.borne)) then
              iso_verif_traceur_nostop=.true.
              ifaux=i
              jfaux=i
          endif !if ((x.gt.-borne).and.(x.lt.borne)) then
        enddo !do ixt=niso+1,ntraciso
        enddo ! do i=1,n
        enddo ! do j=1,m
        

        if (iso_verif_traceur_nostop) then
            write(*,*) 'erreur detectée par iso_verif_nonNAN ', &
     &           'dans iso_verif_traceur_vect'
            write(*,*) ''
            write(*,*) err_msg
            write(*,*) 'x=',x(:,ifaux,jfaux)
            stop
        endif

        end subroutine iso_verif_traceur_noNaN_vect

        
        subroutine iso_verif_trac_masse_vect(x,n,m,err_msg, &
     &            errmax,errmaxrel)
        use isotopes_mod, only: isoName
        implicit none
       
        ! inputs
       integer n,m
       real x(ntraciso,n,m)
       character*(*) err_msg ! message d''erreur à afficher
       real errmax,errmaxrel

       ! locals
       logical iso_verif_traceur_nostop
       integer i,j,ixt,iiso,izone
       integer ifaux,jfaux,ixtfaux
       real xtractot(n,m)
       real xiiso(n,m)

        do iiso=1,niso        
        do j=1,m 
         do i=1,n       
          xtractot(i,j)=0.0
          xiiso(i,j)=x(iiso,i,j)
          do izone=1,nzone
            ixt=itZonIso(izone,iiso) 
            xtractot(i,j)=xtractot(i,j)+x(ixt,i,j)            
          enddo !do izone=1,nzone
         enddo !do i=1,n
        enddo !do j=1,m
        

        call iso_verif_egalite_std_vect( &
     &           xtractot,xiiso, &
     &           err_msg//', verif trac egalite2, iso ' &
     &           //TRIM(isoName(iiso)), &
     &           n,m,errmax,errmaxrel)
        enddo !do iiso=1,niso

        end subroutine iso_verif_trac_masse_vect

        subroutine iso_verif_tracdd_vect(x,n,m,err_msg) 
        use isotopes_mod, only: iso_HDO,iso_eau
        use isotrac_mod, only: strtrac
        implicit none
       
        ! inputs
       integer n,m
       real x(ntraciso,n,m)
       character*(*) err_msg ! message d''erreur à afficher

       ! locals
       integer i,j,iiso,izone,ieau,ixt
       real xiiso(niso,n,m)
       real xeau(n,m)
       integer lnblnk

       if (iso_HDO.gt.0) then
        do izone=1,nzone
          ieau=itZonIso(izone,iso_eau)
          do iiso=1,niso
           ixt=itZonIso(izone,iiso)
           do j=1,m
            do i=1,n
             xiiso(iiso,i,j)=x(ixt,i,j)
            enddo !do i=1,n
           enddo !do j=1,m
          enddo !do iiso=1,niso
           
          do j=1,m
           do i=1,n
            xeau(i,j)=x(ieau,i,j)
           enddo !do i=1,n
          enddo !do j=1,m
          
          call iso_verif_aberrant_vect2Dch( &
     &           xiiso,xeau,err_msg//strtrac(izone),niso,n,m, &
     &           deltalimtrac)
         enddo !do izone=1,nzone
        endif !if (iso_HDO.gt.0) then

        end subroutine iso_verif_tracdd_vect

        subroutine iso_verif_tracpos_vect(x,n,m,err_msg,seuil) 
        implicit none

       ! inputs
       integer n,m
       real x(ntraciso,n,m)
       character*(*) err_msg ! message d''erreur à afficher
       real seuil

       ! locals
       integer i,j,ixt
       logical iso_verif_traceur_nostop
       integer ifaux,jfaux,ixtfaux

        iso_verif_traceur_nostop=.false.        
        do j=1,m
        do i=1,n
        do ixt=niso+1,ntraciso
          if (x(ixt,i,j).lt.-seuil) then
              ifaux=i
              jfaux=j
              ixtfaux=ixt
              iso_verif_traceur_nostop=.true.
          endif
        enddo !do ixt=niso+1,ntraciso 
        enddo !do i=1,n
        enddo !do j=1,m        

        if (iso_verif_traceur_nostop) then
            write(*,*) 'erreur detectée par verif positif ', &
     &           'dans iso_verif_traceur_vect'
            write(*,*) ''
            write(*,*) err_msg
            write(*,*) 'x=',x(:,ifaux,jfaux)
            stop
        endif

        end subroutine iso_verif_tracpos_vect



        subroutine iso_verif_tracnps(x,err_msg)
        use isotrac_mod, only: ridicule_trac

        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher

       ! locals
       !integer iso_verif_tracnps_choix_nostop  

        if (iso_verif_tracnps_choix_nostop(x,err_msg, &
     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
     &       .eq.1) then
                stop
        endif

        end subroutine iso_verif_tracnps

        subroutine iso_verif_tracpos_choix(x,err_msg,seuil)
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher
       real seuil

       ! locals
       !integer iso_verif_tracpos_choix_nostop  

        if (iso_verif_tracpos_choix_nostop(x,err_msg,seuil) &
     &       .eq.1) then
                stop
        endif

        end subroutine iso_verif_tracpos_choix

        subroutine iso_verif_traceur_choix(x,err_msg, &
     &       errmax,errmaxrel,ridicule_trac_loc,deltalimtrac)
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher
       real errmax,errmaxrel,ridicule_trac_loc,deltalimtrac

       ! locals
       !integer iso_verif_traceur_choix_nostop  

        if (iso_verif_traceur_choix_nostop(x,err_msg, &
     &       errmax,errmaxrel,ridicule_trac_loc,deltalimtrac) &
     &       .eq.1) then
                stop
        endif

        end subroutine iso_verif_traceur_choix

        function iso_verif_traceur_nostop(x,err_msg)
        use isotrac_mod, only: ridicule_trac
        !use isotopes_verif, only: errmax,errmaxrel,deltalimtrac
        implicit none
        ! vérifier des choses sur les traceurs
        ! * toutes les zones donne t l'istope total
        ! * pas de deltaD aberrant

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher

       ! output
       integer iso_verif_traceur_nostop

       ! locals
       !integer iso_verif_traceur_choix_nostop  

        iso_verif_traceur_nostop= &
     &       iso_verif_traceur_choix_nostop(x,err_msg, &
     &       errmax,errmaxrel,ridicule_trac,deltalimtrac)

        end function iso_verif_traceur_nostop


      subroutine iso_verif_traceur_justmass(x,err_msg)
        implicit none
        ! on vérifie que noNaN et masse

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher

       ! locals
       !integer iso_verif_traceur_noNaN_nostop
       !integer iso_verif_tracm_choix_nostop

        ! verif noNaN
        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
             stop
        endif 
        
        ! verif masse
        if (iso_verif_tracm_choix_nostop(x,err_msg, &
     &           errmax,errmaxrel).eq.1) then
             stop
        endif    
        
        end subroutine iso_verif_traceur_justmass

        function iso_verif_traceur_jm_nostop(x,err_msg)
        implicit none
        ! on vérifie que noNaN et masse

        ! on prend les valeurs pas défaut pour 
        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
       
       ! inputs
       real x(ntraciso)
       character*(*) err_msg ! message d''erreur à afficher

       ! output
       integer iso_verif_traceur_jm_nostop

       ! locals
!       integer iso_verif_traceur_noNaN_nostop
       !integer iso_verif_tracm_choix_nostop

        iso_verif_traceur_jm_nostop=0
!        ! verif noNaN
!        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
!             iso_verif_traceur_jm_nostop=1
!        endif 
        
        ! verif masse
        if (iso_verif_tracm_choix_nostop(x,err_msg, &
     &           errmax,errmaxrel).eq.1) then
             iso_verif_traceur_jm_nostop=1
        endif   
        
        end function iso_verif_traceur_jm_nostop

        subroutine iso_verif_tag17_q_deltaD_vect(x,n,m,err_msg)
        USE isotopes_mod, ONLY: tnat,iso_eau, ridicule,iso_HDO
        use isotrac_mod, only: option_traceurs,nzone_temp
        implicit none

        ! inputs
        integer n,m
        real x(ntraciso,n,m)
        character*(*) err_msg

        ! locals
        !integer iso_verif_positif_nostop
        !real deltaD
        integer ieau,ixt,ieau1
        integer i,k

        if ((option_traceurs.eq.17).or. &
     &           (option_traceurs.eq.18)) then
        ! verifier que deltaD du tag de la couche la plus haute <
        ! 200 permil, et vérifier que son q est inférieur à 
        ieau=itZonIso(nzone_temp,iso_eau)
        ixt=itZonIso(nzone_temp,iso_HDO)
        ieau1=itZonIso(1,iso_eau)
        do i=1,n
         do k=1,m
           if (x(ieau,i,k).gt.ridicule) then
             if ((x(ixt,i,k)/x(ieau,i,k)/tnat(iso_HDO)-1)*1000 &
     &            .gt.-200.0) then
                write(*,*) err_msg,', vect:deltaDt05 trop fort'
                write(*,*) 'i,k=',i,k
                write(*,*) 'x(:,i,k)=',x(:,i,k)
                stop
             endif !if (iso_verif_positif_nostop(-200.0-deltaD(x(ixt),x(ieau)),
           endif !if (x(ieau).gt.ridicule) then
           if (x(ieau,i,k).gt.2.0e-3) then
                write(*,*) err_msg,', vect:qt05 trop fort'
                write(*,*) 'i,k=',i,k
                write(*,*) 'x(:,i,k)=',x(:,i,k)
                stop
           endif !if (iso_verif_positif_nostop(1.0e-3-x(ieau),
           if (x(iso_eau,i,k).lt.2.0e-3) then
                if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then
                   write(*,*) err_msg,', vect: qt01 trop abondant'
                   write(*,*) 'i,k=',i,k
                   write(*,*) 'ieau1,iso_eau,x(ieau1,i,k),', &
     &                 'x(iso_eau,i,k)=',ieau1,iso_eau, &
     &                  x(ieau1,i,k),x(iso_eau,i,k)  
                   write(*,*) 'x(:,i,k)=',x(:,i,k)
                   stop
                endif !if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then
            endif
          enddo !do k=1,m
        enddo !do i=1,n

        endif !if (option_traceurs.eq.17) then

        end subroutine iso_verif_tag17_q_deltaD_vect


        subroutine iso_verif_tag17_q_deltaD_vect_ret3D(x,n,m,nq,err_msg)
        USE isotopes_mod, ONLY: tnat,iso_eau,iso_HDO,ridicule
        use isotrac_mod, only: option_traceurs,nzone_temp
        implicit none

        ! inputs
        integer n,m,nq
        real x(n,m,nq,ntraciso)
        character*(*) err_msg

        ! locals
        !integer iso_verif_positif_nostop
        !real deltaD
        integer ieau,ixt,ieau1
        integer i,k,iq

        if ((option_traceurs.eq.17).or. &
     &           (option_traceurs.eq.18)) then
        ! verifier que deltaD du tag de la couche la plus haute <
        ! 200 permil, et vérifier que son q est inférieur à 
        ieau=itZonIso(nzone_temp,iso_eau)
        ixt=itZonIso(nzone_temp,iso_HDO)
        ieau1=itZonIso(1,iso_eau)
        do iq=1,nq
        do i=1,n
         do k=1,m
           if (x(i,k,iq,ieau).gt.ridicule) then
             if ((x(i,k,iq,ixt)/x(i,k,iq,ieau)/tnat(iso_HDO)-1)*1000 &
     &            .gt.-200.0) then
                write(*,*) err_msg,', vect:deltaDt05 trop fort'
                write(*,*) 'i,k=',i,k
                write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:)
                stop
             endif !if (iso_verif_positif_nostop(-200.0-deltaD(x(ixt),x(ieau)),
           endif !if (x(ieau).gt.ridicule) then
           if (x(i,k,iq,ieau).gt.2.0e-3) then
                write(*,*) err_msg,', vect:qt05 trop fort'
                write(*,*) 'i,k=',i,k
                write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:)
                stop
           endif !if (iso_verif_positif_nostop(1.0e-3-x(ieau),
           if (x(i,k,iq,iso_eau).lt.2.0e-3) then
                if (x(i,k,iq,ieau1)/x(i,k,iq,iso_eau).gt.0.1) then
                   write(*,*) err_msg,', vect: qt01 trop abondant'
                   write(*,*) 'i,k=',i,k
                   write(*,*) 'ieau1,iso_eau,x(i,k,iq,ieau1),', &
     &                 'x(i,k,iq,ieau)=',ieau1,iso_eau, &
     &                  x(i,k,iq,ieau1),x(i,k,iq,iso_eau)  
                   write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:)
                   stop
                endif !if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then
            endif
          enddo !do k=1,m
        enddo !do i=1,n
        enddo ! do iq=1,nq

        endif !if (option_traceurs.eq.17) then

        end subroutine iso_verif_tag17_q_deltaD_vect_ret3D


#endif
! endif ISOTRAC

END MODULE isotopes_verif_mod

#endif         
! endif ISOVERIF 

