Changeset 5774


Ignore:
Timestamp:
Jul 11, 2025, 5:44:41 PM (3 days ago)
Author:
dcugnet
Message:

Add new isotopes checking routines.
To be tested together with the S. Nguyen's concatenation/decontatenation routine.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmdiso/isotopes_verif_mod.F90

    r4982 r5774  
    44
    55MODULE isotopes_verif_mod
    6 !use isotopes_mod, ONLY:
    7 !#ifdef ISOTRAC
    8 !   USE isotrac_mod, ONLY: nzone
    9 !#endif
    10 USE infotrac_phy, ONLY: ntraciso=>ntiso, niso, itZonIso, nzone
    11 implicit none
    12 save
     6  USE isotopes_mod, ONLY: iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO, small => ridicule, tnat
     7  USE infotrac_phy, ONLY: isoName, isoPhas, isoZone, iqIsoPha, tracers, nqtot, isoSelect, ntraciso => ntiso, &
     8            nbIso,  niso,   ntiso,   nphas,   nzone, itZonIso, isotope, ixIso, isoFamilies, delPhase
     9  USE  strings_mod, ONLY: cat, dispOutliers, lunout, num2str, msg, strStack, strIdx, maxlen
     10  USE phys_local_var_mod, ONLY: qx => qx_seri
     11  USE ioipsl_getin_p_mod, ONLY : getin_p
     12
     13  IMPLICIT NONE
     14
     15  PRIVATE
     16
     17  PUBLIC :: delta2r, delta,  checkNaN, checkEqu, checkMass !, deltaR
     18  PUBLIC :: r2delta, excess, checkPos, checkBnd, checkTrac, abs_err, rel_err
     19
     20  PUBLIC :: errmax, errmaxrel, deltalim, faccond, o17_verif, tmin_verif, faible_evap, deltaDfaible, deltalim_snow, tmax_verif, nlevmaxo17, &
     21            errmax_sol
     22
     23  PUBLIC :: iso_verif_init
     24
     25  PUBLIC :: iso_verif_aberrant,            iso_verif_aberrant_nostop
     26  PUBLIC :: iso_verif_aberrant_choix,      iso_verif_aberrant_choix_nostop
     27  PUBLIC :: iso_verif_aberrant_vect2D
     28  PUBLIC :: iso_verif_aberrant_vect2Dch
     29  PUBLIC :: iso_verif_aberrant_encadre,    iso_verif_aberrant_enc_nostop
     30  PUBLIC :: iso_verif_aberrant_enc_vect2D, iso_verif_aberrant_enc_vect2D_ns
     31  PUBLIC :: iso_verif_aberrant_enc_par2D
     32  PUBLIC :: iso_verif_aberrant_enc_choix_nostop
     33
     34  PUBLIC :: iso_verif_aberrant_o17, iso_verif_aberrant_o17_nostop
     35  PUBLIC :: iso_verif_O18_aberrant, iso_verif_O18_aberrant_nostop
     36  PUBLIC :: iso_verif_O18_aberrant_enc_vect2D
     37
     38  PUBLIC :: iso_verif_egalite,        iso_verif_egalite_nostop
     39  PUBLIC :: iso_verif_egalite_choix,  iso_verif_egalite_choix_nostop
     40  PUBLIC :: iso_verif_egalite_par2D
     41  PUBLIC :: iso_verif_egalite_vect1D, iso_verif_egalite_vect2D
     42  PUBLIC :: iso_verif_egalite_std_vect
     43
     44  PUBLIC :: iso_verif_noNaN, iso_verif_noNaN_nostop
     45  PUBLIC :: iso_verif_noNaN_par2D
     46  PUBLIC :: iso_verif_noNaN_vect1D
     47  PUBLIC :: iso_verif_noNaN_vect2D
     48
     49  PUBLIC :: iso_verif_positif,        iso_verif_positif_nostop
     50  PUBLIC :: iso_verif_positif_vect
     51  PUBLIC :: iso_verif_positif_choix,  iso_verif_positif_choix_nostop
     52  PUBLIC :: iso_verif_positif_choix_vect
     53  PUBLIC :: iso_verif_positif_strict, iso_verif_positif_strict_nostop
     54
     55#ifdef ISOTRAC
     56  PUBLIC :: iso_verif_tag17_q_deltaD_vect
     57  PUBLIC :: iso_verif_tag17_q_deltaD_vect_ret3D
     58  PUBLIC :: iso_verif_tag17_q_deltaD_chns
     59  PUBLIC :: iso_verif_trac17_q_deltaD
     60
     61  PUBLIC :: iso_verif_tracdd_vect
     62  PUBLIC :: iso_verif_tracdD_choix_nostop
     63
     64  PUBLIC :: iso_verif_traceur, iso_verif_traceur_nostop
     65  PUBLIC :: iso_verif_traceur_vect
     66  PUBLIC :: iso_verif_traceur_choix, iso_verif_traceur_choix_nostop
     67
     68  PUBLIC :: iso_verif_tracnps
     69  PUBLIC :: iso_verif_tracnps_vect
     70  PUBLIC :: iso_verif_tracnps_choix_nostop
     71  PUBLIC :: iso_verif_tracpos_vect
     72  PUBLIC :: iso_verif_tracpos_choix, iso_verif_tracpos_choix_nostop
     73
     74  PUBLIC :: iso_verif_trac_masse_vect
     75  PUBLIC :: iso_verif_traceur_justmass, iso_verif_traceur_jm_nostop
     76  PUBLIC :: iso_verif_traceur_noNaN_nostop
     77  PUBLIC :: iso_verif_traceur_noNaN_vect
     78  PUBLIC :: iso_verif_tracm_choix_nostop
     79#endif
     80
     81  PUBLIC :: deltaD, deltaO, dexcess, delta_all, delta_to_R, o17excess
     82
     83
     84! SUPPRIMEE: deltaDfaible_lax=-180.0
     85
     86
     87!===============================================================================================================================
     88  INTERFACE delta2r;   MODULE PROCEDURE    delta2r_0D,    delta2r_1D,    delta2r_2D;                 END INTERFACE delta2r
     89  INTERFACE r2delta;   MODULE PROCEDURE    r2delta_0D,    r2delta_1D,    r2delta_2D;                 END INTERFACE r2delta
     90  INTERFACE checkNaN;  MODULE PROCEDURE   checkNaN_0D,   checkNaN_1D,   checkNaN_2D,    checkNaN_3D; END INTERFACE checkNaN
     91  INTERFACE checkPos;  MODULE PROCEDURE   checkPos_0D,   checkPos_1D,   checkPos_2D,    checkPos_3D; END INTERFACE checkPos
     92  INTERFACE checkBnd;  MODULE PROCEDURE   checkBnd_0D,   checkBnd_1D,   checkBnd_2D,    checkBnd_3D; END INTERFACE checkBnd
     93  INTERFACE checkEqu;  MODULE PROCEDURE   checkEqu_0D,   checkEqu_1D,   checkEqu_2D,    checkEqu_3D; END INTERFACE checkEqu
     94  INTERFACE checkMass; MODULE PROCEDURE checkMassQ_0D, checkMassQ_1D, checkMassQ_2D, checkMassQx_2D; END INTERFACE checkMass
     95  INTERFACE delta;     MODULE PROCEDURE     deltaQ_0D,     deltaQ_1D,     deltaQ_2D,     deltaQx_2D; END INTERFACE delta
     96  INTERFACE deltaR;    MODULE PROCEDURE     deltaR_0D,     deltaR_1D,     deltaR_2D;                 END INTERFACE deltaR
     97  INTERFACE excess;    MODULE PROCEDURE    excessQ_0D,    excessQ_1D,    excessQ_2D, &
     98                                           excessR_0D,    excessR_1D,    excessR_2D,    excessQx_2D; END INTERFACE excess
     99!===============================================================================================================================
     100  TYPE :: r; REAL,    ALLOCATABLE :: r(:); END TYPE r
     101  TYPE :: i; INTEGER, ALLOCATABLE :: i(:); END TYPE i
     102!===============================================================================================================================
     103                                                 !    DESCRIPTION                                     PREVIOUS NAME
     104  REAL, PARAMETER :: abs_err      = 1.e-8,     & !--- Max. absolute error allowed                     errmax
     105                     rel_err      = 1.e-3,     & !--- Max. relative error allowed                     errmax_rel
     106                     max_val      = 1.e19,     & !--- Max. value allowed      (NaN     if greater)    borne
     107                     faccond      = 1.1,       & !--- In liquids, R(max_delD)*faccond is allowed
     108                     min_T        = 5.0,       & !--- Min. realistic temperature value (in K)         Tmin_verif
     109                     max_T        = 1000.0,    & !--- Max. realistic temperature value (in K)         Tmax_verif
     110                     low_delD     = -300.0,    & !--- Max quantity of vapour making outliers from     deltaDfaible
     111                                                 !    balance with ground acceptable
     112                     min_delDevap = 3.0,       & !--- Less demanding with ORCHIDEE evap outliers      faible_evap
     113                     slope_HDO    = 8.0,       & !--- H[2]HO  excess law: slope
     114                     slope_O17    = 0.528        !--- H2[17]O excess law: slope
     115  CHARACTER(LEN=3), PARAMETER ::               &
     116                     progr_HDO    = 'lin',     & !--- H[2]HO  excess law: progression type
     117                     progr_O17    = 'log'        !--- H2[17]O excess law: progression type
     118  REAL,      SAVE :: min_delD,                 & !--- Min. deltaD allowed     (defined in iso.def)    deltaDmin
     119                     max_delD,                 & !--- Max. deltaD for vapour  (outlier if greater)    deltalim
     120                     max_delDtrac,             & !--- Max. deltaD for tracers (defined in iso.def)
     121                     max_delDsnow,             & !--- Max. deltaD for snow    (outlier if greater)
     122                     min_Dexc,                 & !--- Min. value of   D-excess (excess of H2[18]O     dexcess_min
     123                     max_Dexc,                 & !--- Max. value    "    "     w/resp. to H[2]HO)     dexcess_max
     124                     min_O17exc,               & !--- Min. value of O17-excess (excess of H2[17]O     O17excess_bas
     125                     max_O17exc                  !--- Max. value    "    "     w/resp. to H2[18]O)    O17excess_haut
     126
     127  LOGICAL,   SAVE :: checkO17                    !--- Flag to trigger delta(H2[17]O) checking         O17_verif
     128  INTEGER,   SAVE :: max_O17nlev                 !---                                                 nlevmaxO17
     129!$OMP THREADPRIVATE(min_delD, max_delDtrac, min_Dexc, min_O17exc, checkO17)
     130!$OMP THREADPRIVATE(max_delD, max_delDsnow, max_Dexc, max_O17exc, max_O17nlev)
     131!      logical bidouille_anti_divergence ! .TRUE. => xt relaxed to q to avoid errors accumulation leading to divergence
     132                                         ! Moved to wateriso2.h, rzad at execution
     133
     134  !=== QUANTITIES DEPENDING ON THE ISOTOPES FAMILIES (LENGTH: nbIso) ; SET UP FIRST TIME ONLY IN "isoCheck_init"
     135  TYPE(r), ALLOCATABLE, SAVE :: dmin(:),dmax(:),&!--- Thresholds for correct delta  for each isotope (len: isotopes families nb)
     136                                emin(:),emax(:)  !--- Thresholds for correct excess for each isotope (len: isotopes families nb)
     137  TYPE(i), ALLOCATABLE, SAVE :: iMajr(:)         !--- Indices of the major isotopes (to be compared with parent concentration)
     138  CHARACTER(LEN=6), ALLOCATABLE, SAVE ::       &
     139                       sIsoCheck(:)              !--- Activated checking routines keywords list ('npbem' ; 'x' if not activated)
     140
     141  !=== VALUES TO BE SET UP BEFORE CALLING ELEMENTAL FUNCTIONS
     142  REAL,    SAVE :: epsRel, epsAbs, epsPos        !--- NEEDED BY "isEq???".               Set before each checkEqual call
     143  REAL,    SAVE :: min_bnd, max_bnd              !--- NEEDED BY "isBounded".             Set before each checkBound call
     144  REAL,    SAVE :: tnat0                         !--- NEEDED BY "r2delta" and "delta2r". Set before each delta*     call
     145  REAL,             SAVE :: slope                !--- NEEDED BY "isoExcess".             Set using "set_isoParams"
     146  CHARACTER(LEN=3), SAVE :: lawType              !---  "  "
     147  INTEGER,          SAVE :: iso_ref              !---  "  "
    13148
    14149! variables de vérifications
    15150real errmax ! erreur maximale en absolu.
    16 parameter (errmax=1e-8)
     151parameter (errmax=abs_err)
    17152real errmax_sol ! erreur maximale en absolu.
    18153parameter (errmax_sol=5e-7)
    19154      real errmaxrel ! erreur maximale en relatif autorisée
    20155!      parameter (errmaxrel=1e10)
    21       parameter (errmaxrel=1e-3)
     156      parameter (errmaxrel=rel_err)
    22157      real borne ! valeur maximale que n'importe quelle variable peut
    23158                 ! atteindre, en val abs; utile pour vérif que pas NAN
    24       parameter (borne=1e19)
     159      parameter (borne=max_val)
    25160      real, save ::  deltalim ! deltalim est le maximum de deltaD qu'on
    26161                             ! autorise dans la vapeur, au-dela, en suppose que c'est abérrant.
    … …  
    57192!$OMP THREADPRIVATE(dexcess_min,dexcess_max)
    58193
    59       real faccond  ! dans le liquide, on autorise R(deltalim)*faccond.
    60       parameter (faccond=1.1)
     194!      real faccond  ! dans le liquide, on autorise R(deltalim)*faccond.
     195!      parameter (faccond=1.1)
    61196     
    62197!      logical bidouille_anti_divergence ! si true, alors on fait un
    … …  
    71206        !que on peut accepter que la remise à l'équilibre du sol avec
    72207        ! cette vapeur donne des deltaDevap aberrants.
    73         parameter (deltaDfaible=-300)
    74         real deltaDfaible_lax ! deltaD en dessous duquel la vapeur est tellement faible
    75         !que on peut accepter que la remise à l'équilibre du sol avec
    76         ! cette vapeur donne des deltaDevap aberrants.
    77         parameter (deltaDfaible_lax=-180)
     208        parameter (deltaDfaible=low_delD)
    78209
    79210        real faible_evap ! faible évaporation: on est plus laxiste
    80211        !pour les deltaD aberrant dans le cas de l'évap venant d'orchidee
    81         parameter (faible_evap=3.0)
     212        parameter (faible_evap=min_delDevap)
    82213
    83214        real Tmin_verif
    84         parameter (Tmin_verif=5.0) ! en K
     215        parameter (Tmin_verif=min_T) ! en K
    85216        real Tmax_verif
    86         parameter (Tmax_verif=1000.0) ! en K
     217        parameter (Tmax_verif=max_T) ! en K
    87218
    88219
    … …  
    92223CONTAINS
    93224
     225
     226LOGICAL FUNCTION checkiIso(iIso, sub, iFamily) RESULT(lerr)
     227  INTEGER,                    INTENT(IN) :: iIso
     228  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub
     229  INTEGER,          OPTIONAL, INTENT(IN) :: iFamily
     230  CHARACTER(LEN=maxlen) :: subn
     231  INTEGER :: isav
     232  LOGICAL :: ll
     233  subn = 'checkiIso'; IF(PRESENT(sub)) subn = TRIM(sub)//TRIM(subn)
     234  IF(PRESENT(iFamily)) THEN
     235     iSav = ixIso; lerr = isoSelect(iFamily); IF(lerr) RETURN
     236     lerr = iIso <= 0 .OR. iIso > isotope%niso
     237     ll = isoSelect(iSav)
     238  ELSE
     239     lerr = iIso <= 0 .OR. iIso > isotope%niso
     240  END IF
     241  CALL msg('wrong isotopic index iIso = '//TRIM(num2str(iIso))//' (must be >0, <= '//TRIM(num2str(iIso))//')', subn, lerr)
     242END FUNCTION checkiIso
     243
     244
     245LOGICAL FUNCTION isoCheck_init() RESULT(lerr)
     246  CHARACTER(LEN=maxlen) :: isoFamily, modname='isoCheck_init'
     247  INTEGER :: iIso
     248  LOGICAL :: isoCheck
     249
     250  CALL msg('Starting...', modname)
     251  ALLOCATE(dmin(nbIso), dmax(nbIso), emin(nbIso), emax(nbIso), sIsoCheck(nbIso), iMajr(nbIso))
     252  CALL getin_p('isoCheck', isoCheck, .FALSE.)
     253  sIsoCheck = 'xxxxxx'; IF(isoCheck) sIsoCheck = 'npmdeb'  !=== Single keyword for all isotopes families -> CAN BE GENERALIZED !
     254
     255  DO iIso = 1, nbIso
     256     isoFamily = isotope%parent
     257     SELECT CASE(isoFamily)
     258        !-----------------------------------------------------------------------------------------------------------------------
     259        CASE('H2O')
     260        !-----------------------------------------------------------------------------------------------------------------------
     261           iMajr(iIso)%i = [iso_eau]
     262           WRITE(*,*)'iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO = ', iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO
     263           iso_ref     = iso_eau
     264           CALL getin_p('min_delD',     min_delD,    -900.0)
     265           CALL getin_p('max_delD',     max_delD,     300.0)
     266           CALL getin_p('max_delDtrac', max_delDtrac, 300.0)
     267           max_delDsnow =  200.0 + max_delD
     268           WRITE(*,*)'max_delDsnow = ', max_delDsnow
     269           IF(iso_O17 > 0 .AND. iso_O18 > 0) THEN
     270              CALL getin_p('checkO17',    checkO17, .FALSE.)
     271              CALL getin_p('min_O17exc',  min_O17exc, -200.0)
     272              CALL getin_p('max_O17exc',  max_O17exc,  120.0)
     273              CALL getin_p('max_O17nlev', max_O17nlev, 1000)
     274           END IF
     275           IF(iso_HDO > 0 .AND. iso_O18 > 0) THEN
     276              CALL getin_p('min_Dexc', min_Dexc, -100.0)
     277              CALL getin_p('max_Dexc', max_Dexc, 1000.0)
     278           END IF
     279
     280           !=== HANDY VECTORS OF ACCEPTABLE LIMITS FOR delta
     281           ALLOCATE(dmin(iIso)%r(niso)); dmin(iIso)%r(:)=-max_val
     282           ALLOCATE(dmax(iIso)%r(niso)); dmax(iIso)%r(:)=+max_val
     283           IF(iso_HDO /= 0) THEN; dmin(iIso)%r(iso_HDO) = min_delD;     dmax(iIso)%r(iso_HDO) = max_delD;     END IF
     284           IF(iso_O17 /= 0) THEN; dmin(iIso)%r(iso_O17) =-max_val;      dmax(iIso)%r(iso_O17) = max_val;      END IF
     285           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
     286
     287           !=== HANDY VECTORS OF ACCEPTABLE LIMITS FOR excess
     288           ALLOCATE(emin(iIso)%r(niso)); emin(iIso)%r(:)=-max_val
     289           ALLOCATE(emax(iIso)%r(niso)); emax(iIso)%r(:)=+max_val
     290           IF(iso_O17 /= 0) THEN; emin(iIso)%r(iso_O17) = min_O17exc;   emax(iIso)%r(iso_O17) = max_O17exc;   END IF
     291           IF(iso_O18 /= 0) THEN; emin(iIso)%r(iso_O18) = min_Dexc;     emax(iIso)%r(iso_O18) = max_Dexc;     END IF
     292        !-----------------------------------------------------------------------------------------------------------------------
     293        CASE DEFAULT
     294           lerr = .TRUE.
     295           CALL msg('routine not yet set up for isotopes family "'//TRIM(isoFamily)//'"', modname, lerr)
     296        !-----------------------------------------------------------------------------------------------------------------------
     297     END SELECT
     298  END DO
     299
     300END FUNCTION isoCheck_init
     301
     302
     303!===============================================================================================================================
     304!=== ANCILLARY ROUTINES
     305!===============================================================================================================================
     306!=== CHECK ISOTOPIC INDEX "iIso" AND DETERMINE FEW USEFULL QUANTITIES, FOR ELEMENTAL FUNCTIONS AND FOR ERROR MESSAGES:
     307!===      lerr = set_isoParams(iIso, err_tag, sub[, mss][, phas][, iqIso][, iqPar])
     308!===
     309!=== ARGUMENTS:                                                                         Intent  Optional?  Default value
     310!===    * iIso     Index to be checked                                                    IN    MANDATORY
     311!===    * err_tag  Tag for error messages, appended at the beginning of "mss"             IN    OPTIONAL
     312!===    * sub      Calling sequence, appended with current routine name (":" separator)   IN    OPTIONAL
     313!===    * mss      Message in case absurd values are found                                OUT   OPTIONAL
     314!===    * phas     Phase number, needed for the computation of "iqIso" and "iqPar"        IN    OPTIONAL
     315!===    * iqIso    Index in "qx" of the "iIso"th isotope                                  OUT   OPTIONAL
     316!===    * iqPar    Index in "qx" of the "iIso"th isotope parent                           OUT   OPTIONAL
     317!===
     318!=== NOTES:          * SET CORRECT "slope", "iRef" AND "iso_ref" FOR ISOTOPIC EXCESS COMPUTATION.
     319!===                      SO FAR: H[2]HO-excess, H2[17]O-excess (with respect to H2[18]O)
     320!===                 * THE FOLLOWING DELTAs NEED TO BE COMPUTED:  delta-H[2]HO,  delta-H2[17]O,  delta-H2[18]O
     321!===============================================================================================================================
     322LOGICAL FUNCTION set_isoParams(iIso, err_tag, sub, mss, phas, iqIso, iqPar) RESULT(lerr)
     323  INTEGER,                         INTENT(IN)  :: iIso
     324  CHARACTER(LEN=*),                INTENT(IN)  :: err_tag
     325  CHARACTER(LEN=maxlen), OPTIONAL, INTENT(OUT) :: sub, mss
     326  CHARACTER(LEN=*),      OPTIONAL, INTENT(IN)  :: phas
     327  INTEGER,               OPTIONAL, INTENT(OUT) :: iqIso, iqPar
     328  INTEGER :: iq, iPha, ix
     329  LOGICAL :: lExc
     330  CHARACTER(LEN=maxlen) :: iname, subn
     331  subn = 'set_isoParams'; IF(PRESENT(sub)) subn = TRIM(sub)//TRIM(subn)
     332  lerr = checkiIso(iIso, subn); IF(lerr) RETURN
     333  IF(PRESENT(phas)) THEN
     334     iPha = INDEX(isoPhas, phas); lerr = iPha == 0
     335     CALL msg('isotopes family "'//TRIM(isotope%parent)//'" has no "'//TRIM(phas)//'" phase', subn, lerr)
     336  ELSE
     337     iPha = 0; lerr = PRESENT(iqIso) .OR. PRESENT(iqPar)
     338     CALL msg('missing phase, needed to compute iqIso or iqPar', subn, lerr)
     339  END IF
     340  IF(lerr) RETURN
     341  ix = INDEX(sub, ':')
     342  SELECT CASE(sub(ix+1:LEN_TRIM(sub)))
     343     CASE('delta' ); iname = vname(iIso, 'd', iPha)
     344        IF(iPha /= 0)   iq = iqIsoPha(iIso, iPha)
     345        IF(PRESENT(iqIso)) iqIso = iq
     346        IF(PRESENT(iqPar)) iqPar = tracers(iq)%iqParent
     347     CASE('excess'); iname = vname(iIso, 'e', iPha)
     348        IF     (iIso == iso_HDO) THEN
     349           slope = slope_HDO; lawType = progr_HDO; iso_ref = iso_O18
     350        ELSE IF(iIso == iso_O17) THEN
     351          slope = slope_O17; lawType = progr_O17; iso_ref = iso_O18
     352        ELSE
     353           CALL msg('missing parameters for "'//TRIM(iname)//'" (iIso = '//TRIM(num2str(iIso))//'")', subn)
     354           lerr = .TRUE.; RETURN
     355        END IF
     356  END SELECT
     357  IF(PRESENT(mss)) THEN
     358     mss = 'absurd '//TRIM(iname)//' values found'
     359     IF(err_tag /= '') mss = TRIM(err_tag)//':'//TRIM(mss)
     360  END IF
     361END FUNCTION set_isoParams
     362!===============================================================================================================================
     363!=== DETERMINE THE NAME OF THE QUANTITY COMPUTED DEPENDING ON "typ": "d"elta, "e"xcess OR "r"atio.
     364!===============================================================================================================================
     365CHARACTER(LEN=maxlen) FUNCTION vName(iIso, typ, iPha, iZon) RESULT(nam)
     366  INTEGER,           INTENT(IN) :: iIso
     367  CHARACTER(LEN=*),  INTENT(IN) :: typ
     368  INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon
     369  INTEGER :: ix, ip
     370  IF(iIso == 0) THEN
     371     IF(typ == 'd') nam = 'delta'
     372     IF(typ == 'e') nam = 'excess'
     373     IF(typ == 'r') nam = 'ratio'
     374  ELSE
     375     ip = 0;    IF(PRESENT(iPha)) ip = iPha
     376     ix = iIso; IF(PRESENT(iZon)) ix = itZonIso(iZon, iIso)
     377     IF(ip  /=  0 ) nam = tracers(iqIsoPha(ix, iPha))%name
     378     IF(ip  ==  0 ) nam = isoName(ix)
     379     IF(typ == 'd') nam = 'delta-'//TRIM(nam)
     380     IF(typ == 'e') nam = TRIM(nam)//'-excess'
     381     IF(typ == 'r') nam = TRIM(nam)//'-ratio'
     382  END IF
     383END FUNCTION vName
     384!===============================================================================================================================
     385!=== SAME AS vName, FOR SEVERAL INDICES "iIso"
     386!===============================================================================================================================
     387FUNCTION vNames(iIso, typs, iPha, iZon) RESULT(nam)
     388  INTEGER,           INTENT(IN) :: iIso(:)
     389  CHARACTER(LEN=*),  INTENT(IN) :: typs
     390  INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon
     391  CHARACTER(LEN=maxlen)         :: nam(SIZE(iIso))
     392  INTEGER :: it
     393  DO it = 1, SIZE(nam)
     394     nam(it) = vName(iIso(it), typs(it:it), iPha, iZon)
     395  END DO
     396END FUNCTION vNames
     397!===============================================================================================================================
     398!=== LIST OF TWO VARIABLES FOR DELTA CASE: ratio(iIso), delta(iIso)
     399!===============================================================================================================================
     400FUNCTION vNamDe(iIso, iPha, iZon) RESULT(nam)
     401  INTEGER,           INTENT(IN) :: iIso
     402  INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon
     403  CHARACTER(LEN=maxlen) :: nam(2)
     404  nam = vNames([iIso, iso_ref], 'rd', iPha, iZon)
     405END FUNCTION vNamDe
     406!===============================================================================================================================
     407!=== LIST OF THREE VARIABLES FOR EXCESS CASE: delta(iIso), delta(iParent), excess(iIso)
     408!===============================================================================================================================
     409FUNCTION vNamEx(iIso, iPha, iZon) RESULT(nam)
     410  INTEGER,           INTENT(IN) :: iIso
     411  INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon
     412  CHARACTER(LEN=maxlen) :: nam(3)
     413  nam = vNames([iIso, iso_ref, iIso], 'dde', iPha, iZon)
     414END FUNCTION vNamEx
     415!===============================================================================================================================
     416!=== GET THE VARIABLES NAMES LIST (DEFAULT VALUE OR FROM tracers(:)%name, isoName(:) OR OPTIONAL ARGUMENT nam(:))
     417!===============================================================================================================================
     418LOGICAL FUNCTION gettNam(nq, sub, tnam, nam) RESULT(lerr)
     419  INTEGER,                       INTENT(IN)  ::   nq(:)
     420  CHARACTER(LEN=*),              INTENT(IN)  ::  sub
     421  CHARACTER(LEN=*), ALLOCATABLE, INTENT(OUT) :: tnam(:)
     422  CHARACTER(LEN=*), OPTIONAL,    INTENT(IN)  ::  nam(:)
     423  INTEGER :: lq, ln, rk
     424  rk = SIZE(nq); lq = SIZE(nq, DIM=rk)
     425  IF(PRESENT(nam)) THEN; ln = SIZE(nam)
     426     lerr = ALL([1, lq] /= ln)
     427     CALL msg('SIZE(q,'//TRIM(num2str(rk))//')='//TRIM(num2str(lq))//' /= SIZE(nam)='//num2str(ln), sub, lerr); IF(lerr) RETURN
     428     tnam = nam
     429  ELSE
     430     tnam = ['q']
     431     IF(lq == ntiso) tnam = isoName(:)
     432     IF(lq == nqtot) tnam = tracers(:)%name
     433  END IF
     434END FUNCTION gettNam
     435!===============================================================================================================================
     436
     437
     438
     439
     440!===============================================================================================================================
     441!===        BASIC ELEMENTAL FUNCTIONS
     442!===   FUNCTION isNaN  (a)                         Not a Number (NaNs) detection: TRUE if |a|   > max_val
     443!===   FUNCTION isNeg  (a)                         Positiveness detection:        TRUE IF  a    < epsPos
     444!===   FUNCTION isEqAbs(a,b)                (1)    Absolute equality checking:    TRUE if |a-b| < epsAbs
     445!===   FUNCTION isEqRel(a,b)                (1)    Relative equality checking:    TRUE if |a-b|/max(|a|,|b|,1e-18) < epsRel
     446!===   FUNCTION isEqual(a,b)                (1)    Equality checking (both criteria)
     447!===   FUNCTION delta2ratio(ratioIso)       (2)    Delta  function from a   ratio value
     448!===   FUNCTION ratio2delta(ratioIso)       (2)    Ratio  function from a   delta value
     449!===   FUNCTION isoExcess(delIso, delRef)   (3)    Excess function from two delta values
     450!===
     451!===   NOTES:     (1): epsPos and epsAbs must be defined before calling
     452!===              (2): tnat0             must be defined before calling
     453!===              (3): slope and lawType must be defined before calling
     454!===
     455!===============================================================================================================================
     456ELEMENTAL LOGICAL FUNCTION isNaN(a) RESULT(lout)                     !--- IS "a" AN OUTLIER,   ie: a < -max_val or a > max_val ?
     457  REAL, INTENT(IN) :: a
     458  lout = a < -max_val .OR. max_val < a
     459END FUNCTION isNaN
     460!-------------------------------------------------------------------------------------------------------------------------------
     461ELEMENTAL LOGICAL FUNCTION isBounded(a) RESULT(lout)                 !--- IS "a" BOUNDED,      ie:    min_bnd < a < max_bnd ?
     462  REAL, INTENT(IN) :: a
     463  lout= min_bnd < a .AND. a < max_bnd
     464END FUNCTION isBounded
     465!-------------------------------------------------------------------------------------------------------------------------------
     466ELEMENTAL LOGICAL FUNCTION isNeg(a) RESULT(lout)                     !--- IS "a" NEGATIVE OR NOT, ie: a ≤ epsPos ?
     467  REAL, INTENT(IN) :: a
     468  lout = a <= epsPos
     469END FUNCTION isNeg
     470!-------------------------------------------------------------------------------------------------------------------------------
     471ELEMENTAL LOGICAL FUNCTION isPos(a) RESULT(lout)                     !--- IS "a" POSITIVE OR NOT, ie: a ≥ epsPos ?
     472  REAL, INTENT(IN) :: a
     473  lout = a >= epsPos
     474END FUNCTION isPos
     475!-------------------------------------------------------------------------------------------------------------------------------
     476ELEMENTAL LOGICAL FUNCTION isEqAbs(a,b) RESULT(lout)                 !--- IS "a" ABSOLUTELY EQ TO "b", ie: |a-b|<eps ?
     477  REAL, INTENT(IN) :: a, b
     478  lout = abs(a-b) < epsAbs
     479END FUNCTION isEqAbs
     480!-------------------------------------------------------------------------------------------------------------------------------
     481ELEMENTAL LOGICAL FUNCTION isEqRel(a,b) RESULT(lout)                 !--- IS "a" RELATIVELY EQ TO "b", ie: |a-b|/max(|a|,|b|,1e-18)<eps ?
     482  REAL, INTENT(IN) :: a, b
     483  lout = ABS(a-b)/MAX(ABS(a),ABS(b),1e-18) < epsRel
     484END FUNCTION isEqRel
     485!-------------------------------------------------------------------------------------------------------------------------------
     486ELEMENTAL LOGICAL FUNCTION isEqual(a,b) RESULT(lout)                 !--- IS "a" EQUAL TO "b", DUAL REL. AND ABS. CRITERIA
     487  REAL, INTENT(IN) :: a, b
     488  lout = isEqAbs(a,b); IF(lout) lout=isEqRel(a,b)
     489END FUNCTION isEqual
     490!-------------------------------------------------------------------------------------------------------------------------------
     491ELEMENTAL REAL FUNCTION ratio2delta(ratioIso) RESULT(delta)          !--- CONVERT AN ISOTOPIC RATIO INTO A DELTA
     492  REAL, INTENT(IN) :: ratioIso
     493  delta = (ratioIso / tnat0 - 1.0) * 1000.0                          !  /!\ tnat0 must be defined !!!
     494END FUNCTION ratio2delta
     495!-------------------------------------------------------------------------------------------------------------------------------
     496ELEMENTAL REAL FUNCTION delta2ratio(delta) RESULT(ratioIso)          !--- CONVERT A DELTA INTO AN ISOTOPIC RATIO
     497  REAL, INTENT(IN) :: delta
     498  ratioIso = (delta / 1000.0 + 1.0) * tnat0                          !  /!\ tnat0 must be defined !!!
     499END FUNCTION delta2ratio
     500!-------------------------------------------------------------------------------------------------------------------------------
     501ELEMENTAL REAL FUNCTION isoExcess(deltaIso, deltaRef) RESULT(excess) !--- COMPUTE EXCESS USING TWO ISOTOPIC RATIOS
     502  REAL, INTENT(IN) :: deltaIso, deltaRef                             !  /!\ slope and lawType must be defined !!!
     503  SELECT CASE(lawType)
     504    CASE('lin'); excess =              deltaIso        - slope *       deltaRef
     505    CASE('log'); excess = 1.e6*(LOG(1.+deltaIso/1000.) - slope *LOG(1.+deltaRef/1000.))
     506  END SELECT
     507END FUNCTION isoExcess
     508!-------------------------------------------------------------------------------------------------------------------------------
     509
     510
     511!===============================================================================================================================
     512!===   R = delta2r(iIso, delta[, lerr]):  convert an anomaly "delta" into a ratio "R"
     513!=== ARGUMENTS:                                                                     Intent  Optional?
     514!===   * iIso     Isotope index in "isoName"                                        INPUT   MANDATORY
     515!===   * delta    delta ratio for isotope of index "iIso"                           INPUT   MANDATORY
     516!===   * lerr     Error code                                                        OUTPUT  OPTIONAL
     517!===   * R        Isotopic ratio (output value, rank 1 to 3)                        OUTPUT  MANDATORY
     518!===============================================================================================================================
     519REAL FUNCTION delta2r_0D(iIso, delta, lerr) RESULT(R)
     520  INTEGER,           INTENT(IN)  :: iIso
     521  REAL,              INTENT(IN)  :: delta
     522  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
     523  LOGICAL :: ler
     524  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
     525  tnat0 = tnat(iIso); R = delta2ratio(delta)
     526END FUNCTION delta2r_0D
     527!-------------------------------------------------------------------------------------------------------------------------------
     528FUNCTION delta2r_1D(iIso, delta, lerr) RESULT(R)
     529  REAL,              ALLOCATABLE :: R(:)
     530  INTEGER,           INTENT(IN)  :: iIso
     531  REAL,              INTENT(IN)  :: delta(:)
     532  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
     533  LOGICAL :: ler
     534  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
     535  tnat0 = tnat(iIso); R = delta2ratio(delta)
     536END FUNCTION delta2r_1D
     537!-------------------------------------------------------------------------------------------------------------------------------
     538FUNCTION delta2r_2D(iIso, delta, lerr) RESULT(R)
     539  REAL,              ALLOCATABLE :: R(:,:)
     540  INTEGER,           INTENT(IN)  :: iIso
     541  REAL,              INTENT(IN)  :: delta(:,:)
     542  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
     543  LOGICAL :: ler
     544  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
     545  tnat0 = tnat(iIso); R = delta2ratio(delta)
     546END FUNCTION delta2r_2D
     547!===============================================================================================================================
     548!===   delta = r2delta(iIso, R[, lerr]):  convert a ratio "R" into an anomaly "delta"
     549!=== ARGUMENTS:                                                                     Intent  Optional?
     550!===   * iIso      Isotope index in "isoName"                                       INPUT   MANDATORY
     551!===   * R         Isotopic ratio for isotope of index "iIso"                       INPUT   MANDATORY
     552!===   * lerr      Error code                                                       OUTPUT  OPTIONAL
     553!===   * delta     delta ratio (output value, rank 1 to 3)                          OUTPUT  MANDATORY
     554!===============================================================================================================================
     555REAL FUNCTION r2delta_0D(iIso, R, lerr) RESULT(delta)
     556  INTEGER,           INTENT(IN)  :: iIso
     557  REAL,              INTENT(IN)  :: R
     558  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
     559  LOGICAL :: ler
     560  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
     561  tnat0 = tnat(iIso); delta = ratio2delta(R)
     562END FUNCTION r2delta_0D
     563!-------------------------------------------------------------------------------------------------------------------------------
     564FUNCTION r2delta_1D(iIso, R, lerr) RESULT(delta)
     565  REAL,              ALLOCATABLE :: delta(:)
     566  INTEGER,           INTENT(IN)  :: iIso
     567  REAL,              INTENT(IN)  :: R(:)
     568  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
     569  LOGICAL :: ler
     570  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
     571  tnat0 = tnat(iIso); delta = ratio2delta(R)
     572END FUNCTION r2delta_1D
     573!-------------------------------------------------------------------------------------------------------------------------------
     574FUNCTION r2delta_2D(iIso, R, lerr) RESULT(delta)
     575  REAL,              ALLOCATABLE :: delta(:,:)
     576  INTEGER,           INTENT(IN)  :: iIso
     577  REAL,              INTENT(IN)  :: R(:,:)
     578  LOGICAL, OPTIONAL, INTENT(OUT) :: lerr
     579  LOGICAL :: ler
     580  ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN
     581  tnat0 = tnat(iIso); delta = ratio2delta(R)
     582END FUNCTION r2delta_2D
     583!===============================================================================================================================
     584
     585
     586
     587!===============================================================================================================================
     588!===============================================================================================================================
     589!=== CHECK WHETHER INPUT FIELD "q" CONTAINS Not A Number (NaN) VALUES OR NOT
     590!===
     591!===   lerr = checkNaN(q[, err_tag][, subn][, nam][, nmax])
     592!===
     593!=== ARGUMENTS:                                                                     intent  optional?  default value
     594!===   * q         Concentration, scalar or array of rank 1 to 3    (checkNaN)      INPUT   MANDATORY
     595!===   * err_tag   Error message to display if NaNs are found                       INPUT   OPTIONAL   "NaNs found"
     596!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkNaN"
     597!===   * nam       Name(s) of the tracers (last index)                              INPUT   OPTIONAL   "q"
     598!===   * nmax      Maximum number of outliers values to compute for each tracer     INPUT   OPTIONAL   All values
     599!===============================================================================================================================
     600LOGICAL FUNCTION checkNaN_0D(q, err_tag, nam, subn) RESULT(lerr)
     601  REAL,                       INTENT(IN) :: q
     602  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn
     603  CHARACTER(LEN=maxlen)                  :: tag, sub
     604  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     605  tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     606  sub ='checkNaN';   IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     607  lerr = gettNam([1], sub, tnm, [nam]); IF(lerr) RETURN
     608  lerr = dispOutliers( [isNaN(q)], [q], [1], tag, tnm, sub)
     609END FUNCTION checkNaN_0D
     610!-------------------------------------------------------------------------------------------------------------------------------
     611LOGICAL FUNCTION checkNaN_1D(q, err_tag, nam, subn, nmax) RESULT(lerr)
     612  REAL,                       INTENT(IN) :: q(:)
     613  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     614  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     615  CHARACTER(LEN=maxlen)                  :: tag, sub
     616  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     617  tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     618  sub ='checkNaN';   IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     619  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     620  lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     621END FUNCTION checkNaN_1D
     622!-------------------------------------------------------------------------------------------------------------------------------
     623LOGICAL FUNCTION checkNaN_2D(q, err_tag, nam, subn, nmax) RESULT(lerr)
     624  REAL,                       INTENT(IN) :: q(:,:)
     625  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     626  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     627  CHARACTER(LEN=maxlen)                  :: tag, sub
     628  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     629  tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     630  sub ='checkNaN';   IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     631  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     632  lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     633END FUNCTION checkNaN_2D
     634!-------------------------------------------------------------------------------------------------------------------------------
     635LOGICAL FUNCTION checkNaN_3D(q, err_tag, nam, subn, nmax) RESULT(lerr)
     636  REAL,                       INTENT(IN) :: q(:,:,:)
     637  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     638  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     639  CHARACTER(LEN=maxlen)                  :: tag, sub
     640  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     641  tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     642  sub ='checkNaN';   IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     643  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     644  lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     645END FUNCTION checkNaN_3D
     646!===============================================================================================================================
     647
     648
     649!===============================================================================================================================
     650!===============================================================================================================================
     651!=== CHECK WHETHER INPUT FIELD "q" HAVE ONLY POSITIVE VALUES OR NOT
     652!===
     653!=== lerr = checkPos(q[, err_tag][, nam][, subn][, nmax][, tny])                    Check for negative values in "q" .
     654!===
     655!=== ARGUMENTS:                                                                     Intent  Optional?  Default value
     656!===   * q         Concentration, scalar or array of rank 1 to 3    (checkNaN)      INPUT   MANDATORY
     657!===   * err_tag   Error message to display if NaNs are found                       INPUT   OPTIONAL   "Negative values found"
     658!===   * nam       Name(s) of the tracers (last index)                              INPUT   OPTIONAL   "q"
     659!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkPositive"
     660!===   * nmax      Maximum number of outliers values to compute for each tracer     INPUT   OPTIONAL   All values
     661!===   * tny       Value (-) under which a "q" is considered to be negative         INPUT   OPTIONAL   small
     662!===============================================================================================================================
     663LOGICAL FUNCTION checkPos_0D(q, err_tag, nam, subn, tny) RESULT(lerr)
     664  REAL,                       INTENT(IN) :: q
     665  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn
     666  REAL,             OPTIONAL, INTENT(IN) :: tny
     667  CHARACTER(LEN=maxlen)                  :: tag, sub
     668  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     669  tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     670  sub = 'checkPos';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     671  epsPos = -small;              IF(PRESENT(tny )) epsPos = -tny
     672  lerr = gettNam(SHAPE(q), sub, tnm, [nam]); IF(lerr) RETURN
     673  lerr = dispOutliers( [isNeg(q)], [q], [1], tag, tnm, sub)
     674END FUNCTION checkPos_0D
     675!-------------------------------------------------------------------------------------------------------------------------------
     676LOGICAL FUNCTION checkPos_1D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr)
     677  REAL,                       INTENT(IN) :: q(:)
     678  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     679  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     680  REAL,             OPTIONAL, INTENT(IN) :: tny
     681  CHARACTER(LEN=maxlen)                  :: tag, sub
     682  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     683  tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     684  sub = 'checkPos';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     685  epsPos = -small;              IF(PRESENT(tny )) epsPos = -tny
     686  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     687  lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     688END FUNCTION checkPos_1D
     689!-------------------------------------------------------------------------------------------------------------------------------
     690LOGICAL FUNCTION checkPos_2D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr)
     691  REAL,                       INTENT(IN) :: q(:,:)
     692  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     693  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     694  REAL,             OPTIONAL, INTENT(IN) :: tny
     695  CHARACTER(LEN=maxlen)                  :: tag, sub
     696  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     697  tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     698  sub = 'checkPos';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     699  epsPos = -small;              IF(PRESENT(tny )) epsPos = -tny
     700  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     701  lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     702END FUNCTION checkPos_2D
     703!-------------------------------------------------------------------------------------------------------------------------------
     704LOGICAL FUNCTION checkPos_3D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr)
     705  REAL,                       INTENT(IN) :: q(:,:,:)
     706  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     707  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     708  REAL,             OPTIONAL, INTENT(IN) :: tny
     709  CHARACTER(LEN=maxlen)                  :: tag, sub
     710  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     711  tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     712  sub = 'checkPos';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     713  epsPos = -small;              IF(PRESENT(tny )) epsPos = -tny
     714  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     715  lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     716END FUNCTION checkPos_3D
     717!===============================================================================================================================
     718
     719
     720!===============================================================================================================================
     721!===============================================================================================================================
     722!=== CHECK WHETHER INPUT FIELD "q" HAVE VALUES INSIDE GIVEN BOUNDS OR NOT (IE "REASONNABLE" VALUES)
     723!===
     724!=== lerr = checkBnd(q[, err_tag][, nam][, subn][, nmax][, qmin][, qmax])
     725!===
     726!=== ARGUMENTS:                                                                Intent  Optional?  Default value
     727!===   * q       Tested field, scalar or array of rank 1 to 3                  INPUT   MANDATORY
     728!===   * err_tag Error message to display if out-of-bounds are found           INPUT   OPTIONAL   "Negative values found"
     729!===   * nam     Name(s) of the tracers (last index)                           INPUT   OPTIONAL   "q"
     730!===   * subn    Calling subroutine name                                       INPUT   OPTIONAL   "checkBnd"
     731!===   * nmax    Maximum number of outliers values to compute for each tracer  INPUT   OPTIONAL   All values
     732!===   * qmin    Lower bound. Correct values must be > lBnd                    INPUT   OPTIONAL    0
     733!===   * qmax    Upper bound. Correct values must be < uBnd                    INPUT   OPTIONAL   max_val
     734!===============================================================================================================================
     735LOGICAL FUNCTION checkBnd_0D(q, err_tag, nam, subn, qmin, qmax) RESULT(lerr)
     736  REAL,                       INTENT(IN) :: q
     737  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn
     738  REAL,             OPTIONAL, INTENT(IN) :: qmin, qmax
     739  CHARACTER(LEN=maxlen)                  :: tag, sub
     740  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     741  tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     742  sub = 'checkBnd';           IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     743  min_bnd = 0.;               IF(PRESENT(qmin)) min_bnd = qmin
     744  max_bnd = max_val;          IF(PRESENT(qmax)) max_bnd = qmax
     745  lerr = gettNam(SHAPE(q), sub, tnm, [nam]); IF(lerr) RETURN
     746  tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)'
     747  lerr = dispOutliers( [isBounded(q)], [q], [1], tag, tnm, sub)
     748END FUNCTION checkBnd_0D
     749!-------------------------------------------------------------------------------------------------------------------------------
     750LOGICAL FUNCTION checkBnd_1D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr)
     751  REAL,                       INTENT(IN) :: q(:)
     752  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     753  REAL,             OPTIONAL, INTENT(IN) :: qmin, qmax
     754  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     755  CHARACTER(LEN=maxlen)                  :: tag, sub
     756  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     757  tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     758  sub = 'checkBnd';           IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     759  min_bnd = 0.;               IF(PRESENT(qmin)) min_bnd = qmin
     760  max_bnd = max_val;          IF(PRESENT(qmax)) max_bnd = qmax
     761  tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)'
     762  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     763  lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     764END FUNCTION checkBnd_1D
     765!-------------------------------------------------------------------------------------------------------------------------------
     766LOGICAL FUNCTION checkBnd_2D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr)
     767  REAL,                       INTENT(IN) :: q(:,:)
     768  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     769  REAL,             OPTIONAL, INTENT(IN) :: qmin, qmax
     770  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     771  CHARACTER(LEN=maxlen)                  :: tag, sub
     772  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     773  tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     774  sub = 'checkBnd';           IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     775  min_bnd = 0.;               IF(PRESENT(qmin)) min_bnd = qmin
     776  max_bnd = max_val;          IF(PRESENT(qmax)) max_bnd = qmax
     777  tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)'
     778  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     779  lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     780END FUNCTION checkBnd_2D
     781!-------------------------------------------------------------------------------------------------------------------------------
     782LOGICAL FUNCTION checkBnd_3D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr)
     783  REAL,                       INTENT(IN) :: q(:,:,:)
     784  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn
     785  REAL,             OPTIONAL, INTENT(IN) :: qmin, qmax
     786  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     787  CHARACTER(LEN=maxlen)                  :: tag, sub
     788  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     789  tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     790  sub = 'checkBnd';           IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)
     791  min_bnd = 0.;               IF(PRESENT(qmin)) min_bnd = qmin
     792  max_bnd = max_val;          IF(PRESENT(qmax)) max_bnd = qmax
     793  tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)'
     794  lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN
     795  lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax)
     796END FUNCTION checkBnd_3D
     797!===============================================================================================================================
     798
     799
     800!===============================================================================================================================
     801!=== CHECK WHETHER FIELDS "a" AND "b" ARE EQUAL OR NOT (ABSOLUTELY, THE RELATIVELY)
     802!===
     803!=== lerr = checkEqu(a, b[, err_tag][, nam][, subn][, nmax][, absErr][, relErr])
     804!===
     805!=== ARGUMENTS:                                                                     Intent  Optional?  Default value
     806!===   * a, b      Fields,  rank 0 to 3                                             INPUT   MANDATORY
     807!===   * err_tag   Error message to display if fields are not matching              INPUT   OPTIONAL   "mismatching values"
     808!===   * nam       Name(s) of the tracers (last index)                              INPUT   OPTIONAL   "a","b" or from isoNames
     809!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkEqu"
     810!===   * nmax      Maximum number of outliers values to compute for each tracer     INPUT   OPTIONAL   All values
     811!===   * absErr    Maximum value for |q-r| to consider "q" and "q" as equal         INPUT   OPTIONAL    abs_err
     812!===   * relErr    Maximum value for |q-r|/max(|q|,|r|,1e-18)      " "              INPUT   OPTIONAL    rel_err
     813!===============================================================================================================================
     814LOGICAL FUNCTION checkEqu_0D(a, b, err_tag, nam, subn, absErr, relErr) RESULT(lerr)
     815  REAL,                       INTENT(IN) :: a, b
     816  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn
     817  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
     818  CHARACTER(LEN=maxlen)                  :: tag, sub
     819  CHARACTER(LEN=maxlen), ALLOCATABLE     :: tnm(:)
     820  REAL   :: aErr, rErr
     821  tag = 'mismatching value:'; IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     822  tnm = ['a  ','b  ','a-b'];  IF(PRESENT(nam    )) tnm = nam                   !--- Variables names
     823  sub = 'checkEqu';           IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)
     824  epsAbs = abs_err;           IF(PRESENT(absErr )) epsAbs = absErr             !--- Absolute error
     825  epsRel = rel_err;           IF(PRESENT(relErr )) epsRel = relErr             !--- Relative error
     826  lerr = dispOutliers([isEqual(a,b)], cat(a, b, a-b), [1], tag,tnm,sub)
     827END FUNCTION checkEqu_0D
     828!-------------------------------------------------------------------------------------------------------------------------------
     829LOGICAL FUNCTION checkEqu_1D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr)
     830  REAL,                       INTENT(IN) :: a(:), b(:)
     831  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn
     832  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     833  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
     834  CHARACTER(LEN=maxlen) :: tag,  tnm(3), sub
     835  REAL    :: aErr, rErr
     836  tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF
     837  tnm = ['a  ','b  ','a-b'];  IF(PRESENT(nam    )) tnm = nam                   !--- Variables names
     838  sub = 'checkEqu';           IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)
     839  epsAbs = abs_err;           IF(PRESENT(absErr )) epsAbs = absErr             !--- Absolute error
     840  epsRel = rel_err;           IF(PRESENT(relErr )) epsRel = relErr             !--- Relative error
     841  lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax)
     842END FUNCTION checkEqu_1D
     843!-------------------------------------------------------------------------------------------------------------------------------
     844LOGICAL FUNCTION checkEqu_2D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr)
     845  REAL,                       INTENT(IN) :: a(:,:), b(:,:)
     846  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn
     847  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     848  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
     849  CHARACTER(LEN=maxlen) :: tag,  tnm(3), sub=''
     850  REAL    :: aErr, rErr
     851  tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//TRIM(tag); END IF
     852  tnm = ['a  ','b  ','a-b'];  IF(PRESENT(nam    )) tnm = nam                   !--- Variables names
     853  sub = 'checkEqu';           IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)
     854  epsAbs = abs_err;           IF(PRESENT(absErr )) epsAbs = absErr             !--- Absolute error
     855  epsRel = rel_err;           IF(PRESENT(relErr )) epsRel = relErr             !--- Relative error
     856  lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax)
     857END FUNCTION checkEqu_2D
     858!-------------------------------------------------------------------------------------------------------------------------------
     859LOGICAL FUNCTION checkEqu_3D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr)
     860  REAL,                       INTENT(IN) :: a(:,:,:), b(:,:,:)
     861  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn
     862  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     863  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
     864  CHARACTER(LEN=maxlen) :: tag,  tnm(3), sub=''
     865  REAL    :: aErr, rErr
     866  tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//TRIM(tag); END IF
     867  tnm = ['a  ','b  ','a-b'];  IF(PRESENT(nam    )) tnm = nam                   !--- Variables names
     868  sub = 'checkEqu';           IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)
     869  epsAbs = abs_err;           IF(PRESENT(absErr )) epsAbs = absErr             !--- Absolute error
     870  epsRel = rel_err;           IF(PRESENT(relErr )) epsRel = relErr             !--- Relative error
     871  lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax)
     872END FUNCTION checkEqu_3D
     873!===============================================================================================================================
     874
     875
     876!===============================================================================================================================
     877!=== CHECK FOR MASS CONSERVATION, IE. WHETHER THE CONCENTRATION OF A TRACER IS EQUAL TO THE SUM OF ITS CHILDREN CONCENTRATIONS.
     878!=== EXCEPTION (PROBABLY TO BE SUPPRESSED): GENERATION 1 TRACERS ARE COMPARED TO THE MAJOR NEXT GEBNERATION ISOTOPES ONLY
     879!===
     880!=== lerr = checkMass([qt][, err_tag][, subn][, nmax][, absErr][, relErr])
     881!===
     882!=== ARGUMENTS:    Meaning                                                          Intent  Optional?  Default value
     883!===   * qt        Tracers+isotopes+tags stacked along last index, rank 1 to 3      INPUT   OPTIONAL   qx (RANK 3 CASE ONLY !)
     884!===   * err_tag   Error message to display if fields are not matching              INPUT   OPTIONAL   "mismatching values"
     885!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkEqu"
     886!===   * nmax      Maximum number of outliers values to compute for each tracer     INPUT   OPTIONAL   All values
     887!===   * absErr    Maximum value for |q-r| to consider "q" and "q" as equal.        INPUT   OPTIONAL    abs_err
     888!===   * relErr    Maximum value for |q-r|/max(|q|,|r|,1e-18)      " "              INPUT   OPTIONAL    rel_err
     889!===
     890!=== REMARKS:
     891!===   * The concentration of each tracer is compared to the sum of its childrens concentrations (they must match).
     892!===   * In the isotopes case, only a subset of childrens is used (usually only the major contributor).
     893!===   * "qt" last dimension size can either be:
     894!===      - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes (single unknown phase)
     895!===      - nqtot   ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx").
     896!===============================================================================================================================
     897LOGICAL FUNCTION checkMassQ_0D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr)
     898  REAL,                       INTENT(IN) :: qt(:)
     899  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn
     900  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     901  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
     902  CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3)
     903  REAL                  :: qs, qp
     904  INTEGER, ALLOCATABLE :: iqChld(:)
     905  INTEGER :: iIso, iPha,  nqChld, iq, it
     906  tag = '';           IF(PRESENT(err_tag)) tag = err_tag                       !--- Tag for error message
     907  sub  = 'checkMass'; IF(PRESENT(subn   )) sub = TRIM(sub)//':'//subn          !--- Calling subroutine chain
     908  epsAbs = abs_err;   IF(PRESENT(absErr )) epsAbs = absErr                     !--- Absolute error (for isEqual routine)
     909  epsRel = rel_err;   IF(PRESENT(relErr )) epsRel = relErr                     !--- Relative error (for isEqual routine)
     910  IF(SIZE(qt) == ntiso+1) THEN
     911     iqChld = iMajr(it)%i
     912     pNam = isotope%parent
     913     DO it = 0, niso
     914        IF(it /= 0) iqChld = itZonIso(:,it)
     915        IF(it /= 0) pNam = isoName(it)
     916        qp = qt(it+1)
     917        qs = SUM(qt(1+iqChld), DIM=1)                                          !--- Tracer compared to its major isotopes
     918        sNam = TRIM(strStack(isoName(iqChld),'+'))                             !--- Names of contributing childrens
     919        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
     920        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
     921        lerr = dispOutliers([isEqual(qp, qs)], cat(qp, qs, qp-qs), SHAPE(qt), tag, nam, sub, nmax); IF(lerr) RETURN
     922     END DO
     923  ELSE
     924     DO iq = 1, SIZE(qt)
     925        pNam   = tracers(iq)%name                                              !--- Current tracer name with phase
     926        nqChld = tracers(iq)%nqChildren                                        !--- Number of childrens (next generation only)
     927        IF(nqChld == 0) CYCLE                                                  !--- No descendants
     928        iIso = tracers(iq)%iso_iGroup                                          !--- Isotopes family index
     929        iPha = tracers(iq)%iso_iPhase                                          !--- Index of current phase in isoPhases(:)
     930        IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha)                   !--- Indices of major contributors in qt
     931        IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld)                  !--- Indices of every contributor  in qt
     932        qp = qt(iq)                                                            !--- Linearized (reprofiled) parent field
     933        qs = SUM(qt(iqChld))                                                   !--- Sum of contributing child(ren) field(s)
     934        sNam = TRIM(strStack(tracers(iqChld)%name ,'+'))                       !--- Names of contributing child(ren)
     935        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
     936        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
     937        lerr = dispOutliers([isEqual(qp, qs)], cat(qp, qs, qp-qs), [1], tag, nam, sub, nmax); IF(lerr) RETURN
     938     END DO
     939  END IF
     940END FUNCTION checkMassQ_0D
     941!-------------------------------------------------------------------------------------------------------------------------------
     942LOGICAL FUNCTION checkMassQ_1D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr)
     943  REAL,             TARGET,   INTENT(IN) :: qt(:,:)
     944  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn
     945  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     946  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
     947  CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3)
     948  REAL, DIMENSION(SIZE(qt(:,1))) :: qs, qp
     949  INTEGER, ALLOCATABLE :: iqChld(:)
     950  INTEGER :: iIso, iPha,  nqChld, iq, it
     951  tag = '';           IF(PRESENT(err_tag)) tag = err_tag                       !--- Tag for error message
     952  sub  = 'checkMass'; IF(PRESENT(subn   )) sub = TRIM(sub)//':'//subn          !--- Calling subroutine chain
     953  epsAbs = abs_err;   IF(PRESENT(absErr )) epsAbs = absErr                     !--- Absolute error (for isEqual routine)
     954  epsRel = rel_err;   IF(PRESENT(relErr )) epsRel = relErr                     !--- Relative error (for isEqual routine)
     955  IF(SIZE(qt, DIM=2) == ntiso+1) THEN
     956     iqChld = iMajr(it)%i
     957     pNam = isotope%parent
     958     DO it = 0, niso
     959        IF(it /= 0) iqChld = itZonIso(:,it)
     960        IF(it /= 0) pNam = isoName(it)
     961        qp = PACK(qt(:,it+1), MASK=.TRUE.)
     962        qs = SUM(qt(:,1+iqChld), DIM=2)                                        !--- Tracer compared to its major isotopes
     963        sNam = TRIM(strStack(isoName(iqChld),'+'))                             !--- Names of contributing childrens
     964        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
     965        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
     966        lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,1)), tag, nam, sub, nmax); IF(lerr) RETURN
     967     END DO
     968  ELSE
     969     DO iq = 1, SIZE(qt, DIM=2)
     970        pNam   = tracers(iq)%name                                              !--- Current tracer name with phase
     971        nqChld = tracers(iq)%nqChildren                                        !--- Number of childrens (next generation only)
     972        IF(nqChld == 0) CYCLE                                                  !--- No descendants
     973        iIso = tracers(iq)%iso_iGroup                                          !--- Isotopes family index
     974        iPha = tracers(iq)%iso_iPhase                                          !--- Index of current phase in isoPhases(:)
     975        IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha)                   !--- Indices of major contributors in qt
     976        IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld)                  !--- Indices of every contributor  in qt
     977        qp = PACK(qt(:,iq), MASK=.TRUE.)                                       !--- Linearized (reprofiled) parent field
     978        qs = PACK(SUM(qt(:,iqChld), DIM=2), MASK=.TRUE.)                       !--- Sum of contributing child(ren) field(s)
     979        sNam = TRIM(strStack(tracers(iqChld)%name ,'+'))                       !--- Names of contributing child(ren)
     980        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
     981        nam(1) = pnam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
     982       lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,1)), tag, nam, sub, nmax); IF(lerr) RETURN
     983     END DO
     984  END IF
     985END FUNCTION checkMassQ_1D
     986!-------------------------------------------------------------------------------------------------------------------------------
     987LOGICAL FUNCTION checkMassQ_2D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr)
     988  REAL,             TARGET,   INTENT(IN) :: qt(:,:,:)
     989  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn
     990  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     991  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
     992  CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3)
     993  REAL, DIMENSION(SIZE(qt(:,:,1))) :: qs, qp
     994  INTEGER, ALLOCATABLE :: iqChld(:)
     995  INTEGER :: iIso, iPha,  nqChld, iq, it
     996  tag = '';           IF(PRESENT(err_tag)) tag = err_tag                       !--- Tag for error message
     997  sub  = 'checkMass'; IF(PRESENT(subn   )) sub = TRIM(sub)//':'//subn          !--- Calling subroutine chain
     998  epsAbs = abs_err;   IF(PRESENT(absErr )) epsAbs = absErr                     !--- Absolute error (for isEqual routine)
     999  epsRel = rel_err;   IF(PRESENT(relErr )) epsRel = relErr                     !--- Relative error (for isEqual routine)
     1000  IF(SIZE(qt, DIM=3) == ntiso+1) THEN
     1001     iqChld = iMajr(it)%i
     1002     pNam = isotope%parent
     1003     DO it = 0, niso
     1004        IF(it /= 0) iqChld = itZonIso(:,it)
     1005        IF(it /= 0) pNam = isoName(it)
     1006        qp = PACK(qt(:,:,it+1), MASK=.TRUE.)
     1007        qs = PACK(SUM(qt(:,:,1+iqChld), DIM=3), MASK=.TRUE.)                   !--- Tracer compared to its major isotopes
     1008        sNam = TRIM(strStack(isoName(iqChld),'+'))                             !--- Names of contributing childrens
     1009        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
     1010        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
     1011        lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,:,1)), tag, nam, sub, nmax); IF(lerr) RETURN
     1012     END DO
     1013  ELSE
     1014     DO iq = 1, SIZE(qt, DIM=3)
     1015        pNam   = tracers(iq)%name                                              !--- Current tracer name with phase
     1016        nqChld = tracers(iq)%nqChildren                                        !--- Number of childrens (next generation only)
     1017        IF(nqChld == 0) CYCLE                                                  !--- No descendants
     1018        iIso = tracers(iq)%iso_iGroup                                          !--- Isotopes family index
     1019        iPha = tracers(iq)%iso_iPhase                                          !--- Index of current phase in isoPhases(:)
     1020        IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha)                   !--- Indices of major contributors in qt
     1021        IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld)                  !--- Indices of every contributor  in qt
     1022        qp = PACK(qt(:,:,iq), MASK=.TRUE.)                                     !--- Linearized (reprofiled) parent field
     1023        qs = PACK(SUM(qt(:,:,iqChld), DIM=3), MASK=.TRUE.)                     !--- Sum of contributing child(ren) field(s)
     1024        sNam = TRIM(strStack(tracers(iqChld)%name ,'+'))                       !--- Names of contributing child(ren)
     1025        tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam)
     1026        nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference'                   !--- Names for the table
     1027        lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,:,1)), tag, nam, sub, nmax); IF(lerr) RETURN
     1028     END DO
     1029  END IF
     1030END FUNCTION checkMassQ_2D
     1031!-------------------------------------------------------------------------------------------------------------------------------
     1032LOGICAL FUNCTION checkMassQx_2D(err_tag, subn, nmax, absErr, relErr) RESULT(lerr)
     1033  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn
     1034  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     1035  REAL,             OPTIONAL, INTENT(IN) :: absErr, relErr
     1036  lerr = checkMassQ_2D(qx, err_tag, subn, nmax, absErr, relErr)
     1037END FUNCTION checkMassQx_2D
     1038!===============================================================================================================================
     1039
     1040
     1041!===============================================================================================================================
     1042!===============================================================================================================================
     1043!=== COMPUTE THE ISOTOPIC ANOMALY "DELTA" AND CHECK THE VALUES
     1044!===
     1045!=== lerr = deltaR(rIso, iIso[, err_tag][, subn][, nmax][, delMax][, delIso]                [, lCheck])
     1046!=== lerr = deltaQ(qt,   iIso[, err_tag][, subn][, nmax][, delMax][, delIso][, Phas][, qmin][, lCheck])
     1047!=== lerr = deltaQ(      iIso[, err_tag][, subn][, nmax][, delMax][, delIso][, Phas][, qmin][, lCheck])
     1048!===
     1049!=== ARGUMENTS:    Meaning                                                          Intent  Optional?  Default      Present:
     1050!===   * rIso      Field of rank 0 to 3                                             INPUT   MANDATORY               CASE 1 ONLY
     1051!===   * qt        Stacked fields, rank 1 to 3 (last index for species)             INPUT   MANDATORY               CASE 2 ONLY
     1052!===   * iIso      Index of tested isotope in "isoName"                             INPUT   MANDATORY
     1053!===   * err_tag   Error message to display if fields are not matching              INPUT   OPTIONAL   "??????????"
     1054!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "delta[R]"
     1055!===   * nmax      Maximum number of lines to print for outliers                    INPUT   OPTIONAL   <all lines>  RANK>0 ONLY
     1056!===   * delMax    Maximum value for the isotopic delta                             INPUT   OPTIONAL    dmax(iIso)
     1057!===   * delIso    Isotopic anomaly delta for isotope iIso                          OUTPUT  OPTIONAL
     1058!===   * phas      Phase                                                            INPUT   OPTIONAL                CASE 2 ONLY
     1059!===   * qmin      Values lower than this threshold are not checked                 INPUT   OPTIONAL    small       CASE 2 ONLY
     1060!===   * lCheck    Trigger checking routines with tables for strage values          INPUT   OPTIONAL   .FALSE.
     1061!===
     1062!=== REMARKS:
     1063!===   * In the version 2 and 3 of the routine, values are checked only if q>=qmin
     1064!===   * In the version 2, the size of the last dimension of "qt" can either be:
     1065!===     - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes.
     1066!===     - nqtot   ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx").
     1067!===============================================================================================================================
     1068LOGICAL FUNCTION deltaR_0D(rIso, iIso, err_tag, subn, delMax, delIso, lCheck) RESULT(lerr)
     1069  REAL,                       INTENT(IN)  :: rIso
     1070  INTEGER,                    INTENT(IN)  :: iIso
     1071  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn
     1072  REAL,             OPTIONAL, INTENT(IN)  :: delMax
     1073  REAL,             OPTIONAL, INTENT(OUT) :: delIso
     1074  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1075  CHARACTER(LEN=maxlen) :: tag, sub, mss
     1076  REAL    :: dmn, dmx, dIso
     1077  LOGICAL :: lc
     1078  IF(niso == 0) RETURN
     1079  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
     1080  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1081  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1082  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
     1083  dmn = dMin(ixIso)%r(iIso)
     1084  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Check the index + define few vars
     1085  mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
     1086  tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                                 !=== Compute delta(isotope iIso in phase iPhas)
     1087  IF(PRESENT(delIso)) delIso = dIso
     1088  IF(.NOT.lc)   RETURN
     1089  lerr = dispOutliers([dIso<dmn .OR. dIso>dmx], cat(rIso, dIso), [1], mss, vNamDe(iIso), sub)
     1090END FUNCTION deltaR_0D
     1091!-------------------------------------------------------------------------------------------------------------------------------
     1092LOGICAL FUNCTION deltaR_1D(rIso, iIso, err_tag, subn, nmax, delMax, delIso, lCheck) RESULT(lerr)
     1093  REAL,                       INTENT(IN)  :: rIso(:)
     1094  INTEGER,                    INTENT(IN)  :: iIso
     1095  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn
     1096  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
     1097  REAL,             OPTIONAL, INTENT(IN)  :: delMax
     1098  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(rIso,1))
     1099  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1100  CHARACTER(LEN=maxlen) :: tag, sub, mss
     1101  REAL    :: dmn, dmx, dIso(SIZE(rIso,1))
     1102  LOGICAL :: lc
     1103  IF(niso == 0) RETURN
     1104  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
     1105  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1106  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1107  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
     1108  dmn = dMin(ixIso)%r(iIso)
     1109  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Check the index + define few vars
     1110  mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
     1111  tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                                 !=== Compute delta(isotope iIso in phase iPhas)
     1112  IF(PRESENT(delIso)) delIso = dIso
     1113  IF(.NOT.lc)   RETURN
     1114  lerr = dispOutliers(PACK(dIso<dmn .OR. dIso>dmx,.TRUE.), cat(rIso, dIso), SHAPE(dIso), mss, vNamDe(iIso), sub, nmax)
     1115END FUNCTION deltaR_1D
     1116!-------------------------------------------------------------------------------------------------------------------------------
     1117LOGICAL FUNCTION deltaR_2D(rIso, iIso, err_tag, subn, nmax, delMax, delIso, lCheck) RESULT(lerr)
     1118  REAL,                       INTENT(IN)  :: rIso(:,:)
     1119  INTEGER,                    INTENT(IN)  :: iIso
     1120  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn
     1121  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
     1122  REAL,             OPTIONAL, INTENT(IN)  :: delMax
     1123  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(rIso,1),SIZE(rIso,2))
     1124  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1125  CHARACTER(LEN=maxlen) :: tag, sub, mss
     1126  REAL    :: dmn, dmx, emx, dIso(SIZE(rIso,1),SIZE(rIso,2))
     1127  LOGICAL :: lc
     1128  IF(niso == 0) RETURN
     1129  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
     1130  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1131  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1132  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
     1133  dmn = dMin(ixIso)%r(iIso)
     1134  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Check the index + define few vars
     1135  mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
     1136  tnat0 = tnat(iIso); dIso = ratio2delta(rIso)
     1137  IF(PRESENT(delIso)) delIso = dIso
     1138  IF(.NOT.lc)   RETURN
     1139  lerr = dispOutliers( [dIso<dmn.OR.dIso>dmx], cat(PACK(rIso,.TRUE.),PACK(dIso,.TRUE.)), SHAPE(dIso),mss,vNamDe(iIso),sub,nmax)
     1140END FUNCTION deltaR_2D
     1141!===============================================================================================================================
     1142LOGICAL FUNCTION deltaQ_0D(qt, iIso, err_tag, subn, delMax, delIso, phas, qmin, lCheck) RESULT(lerr)
     1143  REAL, DIMENSION(:),         INTENT(IN)  :: qt
     1144  INTEGER,                    INTENT(IN)  :: iIso
     1145  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
     1146  REAL,             OPTIONAL, INTENT(IN)  :: delMax, qmin
     1147  REAL,             OPTIONAL, INTENT(OUT) :: delIso
     1148  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1149  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
     1150  REAL    :: dmn, dmx, qmn, rIso, dIso
     1151  LOGICAL :: lc, lq
     1152  INTEGER :: iqIso, iqPar, iZon, iPha, ii
     1153  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
     1154  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1155  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1156  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
     1157  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
     1158  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
     1159  lq = SIZE(qt, DIM=1) == niso+1
     1160
     1161  !=== CHECK PHASE
     1162  IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
     1163  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
     1164  IF(lerr)       RETURN
     1165  dmn = dMin(ixIso)%r(iIso)
     1166  m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
     1167
     1168  !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
     1169  ii = iIso
     1170  DO iZon = 0, nzone
     1171     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     1172     DO iPha = 1, nphas
     1173        p = isoPhas(iPha:iPha)
     1174        lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN    !--- Ckeck the index + define few vars
     1175        mss = TRIM(mss)//TRIM(m)
     1176        IF(     lq) rIso = qt(iIso+1)/qt(1)                                    !--- "qt(1+niso)": parent, then isotopes
     1177        IF(.NOT.lq) rIso = qt(iqIso)/qt(iqPar)                                 !--- Fields taken from a "qt" similar to "qx"
     1178        tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                           !=== Compute delta(isotope iIso in phase "p")
     1179        IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso       !--- Return  delta(iIso, phase=p)
     1180        IF(.NOT.lc) CYCLE
     1181        lerr = dispOutliers( [qt(iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], cat(rIso, dIso), [1], mss, vNamde(iIso), sub)
     1182        IF(lerr) RETURN
     1183     END DO
     1184     IF(.NOT.lc) EXIT
     1185  END DO
     1186END FUNCTION deltaQ_0D
     1187!-------------------------------------------------------------------------------------------------------------------------------
     1188LOGICAL FUNCTION deltaQ_1D(qt, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr)
     1189  REAL, DIMENSION(:,:),       INTENT(IN)  :: qt
     1190  INTEGER,                    INTENT(IN)  :: iIso
     1191  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
     1192  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
     1193  REAL,             OPTIONAL, INTENT(IN)  :: delMax, qmin
     1194  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(qt,1))
     1195  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1196  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
     1197  REAL    :: dmn, dmx, qmn
     1198  REAL, DIMENSION(SIZE(qt,1)) :: rIso, dIso
     1199  LOGICAL :: lc, lq
     1200  INTEGER :: iqIso, iqPar, iZon, iPha, ii
     1201  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
     1202  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1203  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1204  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
     1205  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
     1206  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
     1207  lq = SIZE(qt, DIM=2) == niso+1
     1208
     1209  !=== CHECK PHASE
     1210  IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
     1211  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
     1212  IF(lerr)       RETURN
     1213  dmn = dMin(ixIso)%r(iIso)
     1214  m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
     1215
     1216  !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
     1217  ii = iIso
     1218  DO iZon = 0, nzone
     1219     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     1220     DO iPha = 1, nphas
     1221        p = isoPhas(iPha:iPha)
     1222        lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN    !--- Ckeck the index + define few vars
     1223        mss = TRIM(mss)//TRIM(m)
     1224        IF(     lq) rIso = qt(:,iIso+1)/qt(:,1)                                !--- "qt(1+niso)": parent, then isotopes
     1225        IF(.NOT.lq) rIso = qt(:,iqIso)/qt(:,iqPar)                             !--- Fields taken from a "qt" similar to "qx"
     1226        tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                           !=== Compute delta(iIso, phase=p)
     1227        IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso       !--- Return  delta(iIso, phase=p)
     1228        IF(.NOT.lc) CYCLE
     1229        lerr = dispOutliers( [qt(:,iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], &
     1230                        cat(PACK(rIso,.TRUE.), PACK(dIso,.TRUE.)), SHAPE(dIso), mss, vNamde(iIso), sub, nmax)
     1231        IF(lerr) RETURN
     1232     END DO
     1233     IF(.NOT.lc) EXIT
     1234  END DO
     1235END FUNCTION deltaQ_1D
     1236!-------------------------------------------------------------------------------------------------------------------------------
     1237LOGICAL FUNCTION deltaQ_2D(qt, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr)
     1238  REAL, DIMENSION(:,:,:),     INTENT(IN)  :: qt
     1239  INTEGER,                    INTENT(IN)  :: iIso
     1240  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
     1241  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
     1242  REAL,             OPTIONAL, INTENT(IN)  :: delMax, qmin
     1243  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(qt,1),SIZE(qt,2))
     1244  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1245  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
     1246  REAL    :: dmn, dmx, qmn
     1247  REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)) :: rIso, dIso
     1248  LOGICAL :: lc, lq
     1249  INTEGER :: iqIso, iqPar, iZon, iPha, ii
     1250  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag                !--- Tag for error message
     1251  sub = 'delta';             IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1252  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1253  lc  = .FALSE.;             IF(PRESENT(lCheck )) lc = lCheck                  !--- Checking routines execution trigger
     1254  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
     1255  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
     1256  lq = SIZE(qt, DIM=3) == niso+1
     1257
     1258  !=== CHECK PHASE
     1259  IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
     1260  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
     1261  IF(lerr)       RETURN
     1262  dmn = dMin(ixIso)%r(iIso)
     1263  m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')'
     1264
     1265  !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
     1266  ii = iIso
     1267  DO iZon = 0, nzone
     1268     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     1269     DO iPha = 1, nphas
     1270        p = isoPhas(iPha:iPha)
     1271        lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN    !--- Ckeck the index + define few vars
     1272        mss = TRIM(mss)//' '//TRIM(m)
     1273        IF(     lq) rIso = qt(:,:,iIso+1)/qt(:,:,1)                            !--- "qt(1+niso)": parent, then isotopes
     1274        IF(.NOT.lq) rIso = qt(:,:,iqIso)/qt(:,:,iqPar)                         !--- Fields taken from a "qt" similar to "qx"
     1275        tnat0 = tnat(iIso); dIso = ratio2delta(rIso)                           !=== Compute delta(iIso, phase=p)
     1276        IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso       !--- Return  delta(iIso, phase=p)
     1277        IF(.NOT.lc) CYCLE
     1278        lerr = dispOutliers( [qt(:,:,iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], &
     1279                        cat(PACK(rIso,.TRUE.), PACK(dIso,.TRUE.)), SHAPE(dIso), mss, vNamDe(iIso), sub, nmax)
     1280        IF(lerr) RETURN
     1281     END DO
     1282     IF(.NOT.lc) EXIT
     1283  END DO
     1284END FUNCTION deltaQ_2D
     1285!===============================================================================================================================
     1286LOGICAL FUNCTION deltaQx_2D(iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr)
     1287  INTEGER,                    INTENT(IN)  :: iIso
     1288  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
     1289  INTEGER,          OPTIONAL, INTENT(IN)  :: nmax
     1290  REAL,             OPTIONAL, INTENT(IN)  :: delMax, qmin
     1291  REAL,             OPTIONAL, INTENT(OUT) :: delIso(SIZE(qx,1),SIZE(qx,2))
     1292  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1293  lerr = deltaQ_2D(qx, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck)
     1294END FUNCTION deltaQx_2D
     1295!===============================================================================================================================
     1296
     1297
     1298!===============================================================================================================================
     1299!===============================================================================================================================
     1300!=== COMPUTE THE ISOTOPIC EXCESS AND CHECK THE VALUE (D-excess and O17-excess so far)
     1301!===
     1302!=== lerr = excess(rIso,rRef,iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef]              [,lCheck])
     1303!=== lerr = excess(qt,       iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef][,phas][,qmin][,lCheck])
     1304!=== lerr = excess(          iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef][,phas][,qmin][,lCheck])
     1305!===
     1306!=== ARGUMENTS:                                                                     Intent  Optional?  Default value
     1307!===   * rIso,rRef Considered isotope + reference ratios: field of rank 0 to 3      INPUT   MANDATORY              (CASE 1 ONLY)
     1308!===   * qt        Stacked fields, rank 1 to 3 (last index for species)             INPUT   MANDATORY              (CASE 2 ONLY)
     1309!===   * iIso      Index of tested species                                          INPUT   MANDATORY
     1310!===   * err_tag   Error tag sent from the calling routine                          INPUT   OPTIONAL    none
     1311!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "excess"
     1312!===   * nmax      Maximum number of lines to print for outliers                    INPUT   OPTIONAL   <all lines> (RANK>0 ONLY)
     1313!===   * delMax    Maximum value for the isotopic delta                             INPUT   OPTIONAL
     1314!===   * excMax    Maximum value for the isotopic excess                            INPUT   OPTIONAL    dmax(iIso)
     1315!===   * delIso    Isotopic anomaly delta for isotope iIso                          OUTPUT  OPTIONAL
     1316!===   * excIso    Isotopic excess for isotope iIso                                 OUTPUT  OPTIONAL
     1317!===   * delRef    Isotopic anomaly delta for isotope iIso reference                OUTPUT  OPTIONAL
     1318!===   * phas      Phase name (one element of "isoPhas")                            INPUT   MANDATORY              (CASE 2 ONLY)
     1319!===   * qmin      Values lower than this threshold are not checked                 INPUT   OPTIONAL    small      (CASE 2 ONLY)
     1320!===
     1321!=== REMARKS:
     1322!===   * In the version 2 and 3 of the routine, values are checked only if q>=qmin
     1323!===   * In the version 2, the size of the last dimension of "qt" can either be:
     1324!===     - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes.
     1325!===     - nqtot   ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx").
     1326!===============================================================================================================================
     1327LOGICAL FUNCTION excessR_0D(rIso, rRef, iIso, err_tag, subn, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr)
     1328  REAL,                       INTENT(IN)  :: rIso, rRef
     1329  INTEGER,                    INTENT(IN)  :: iIso
     1330  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn
     1331  REAL,             OPTIONAL, INTENT(IN)  :: delMax, excMax
     1332  REAL,             OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
     1333  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1334  CHARACTER(LEN=maxlen) :: tag, sub, mss
     1335  REAL    :: deIso, deRef, exIso, dmx, drx, emn, emx
     1336  LOGICAL :: lc
     1337  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
     1338  sub = 'excess';            IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)  !--- Calling subroutine name
     1339  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1340  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
     1341  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
     1342  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Ckeck the index + define few vars
     1343  emn = eMin(ixIso)%r(iIso)
     1344  drx = dMax(ixIso)%r(iso_ref)
     1345  mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
     1346  lerr = deltaR_0D(rIso, iIso,    tag, sub, dmx, deIso, lCheck) &              !=== COMPUTE deltaIso AND deltaRef
     1347    .OR. deltaR_0D(rRef, iso_ref, tag, sub, drx, deRef, lCheck)
     1348  IF(lerr)    RETURN
     1349  exIso = isoExcess(deIso, deRef)                                              !=== COMPUTE THE EXCESS
     1350  IF(PRESENT(delIso)) delIso = deIso
     1351  IF(PRESENT(excIso)) excIso = exIso
     1352  IF(PRESENT(delRef)) delRef = deRef
     1353  IF(.NOT.lc) RETURN
     1354  lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), [1], mss, vNamEx(iIso), sub)
     1355END FUNCTION excessR_0D
     1356!-------------------------------------------------------------------------------------------------------------------------------
     1357LOGICAL FUNCTION excessR_1D(rIso, rRef, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr)
     1358  REAL, DIMENSION(:),                      INTENT(IN)  :: rIso, rRef
     1359  INTEGER,                                 INTENT(IN)  :: iIso
     1360  CHARACTER(LEN=*),              OPTIONAL, INTENT(IN)  :: err_tag, subn
     1361  INTEGER,                       OPTIONAL, INTENT(IN)  :: nmax
     1362  REAL,                          OPTIONAL, INTENT(IN)  :: delMax, excMax
     1363  REAL, DIMENSION(SIZE(rIso,1)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
     1364  LOGICAL,                       OPTIONAL, INTENT(IN)  :: lCheck
     1365  CHARACTER(LEN=maxlen) :: tag, sub, mss, m
     1366  REAL, DIMENSION(SIZE(rIso,1)) :: deIso, deRef, exIso
     1367  REAL    :: dmx, drx, emn, emx
     1368  LOGICAL :: lc
     1369  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
     1370  sub = 'excess';            IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)  !--- Calling subroutine name
     1371  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1372  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
     1373  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
     1374  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Ckeck the index + define few vars
     1375  emn = eMin(ixIso)%r(iIso)
     1376  drx = dMax(ixIso)%r(iso_ref)
     1377  mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
     1378  lerr = deltaR_1D(rIso, iIso,    tag, sub, nmax, dmx, deIso, lCheck) &        !=== COMPUTE deltaIso AND deltaRef
     1379    .OR. deltaR_1D(rRef, iso_ref, tag, sub, nmax, drx, deRef, lCheck)
     1380  IF(lerr)    RETURN
     1381  exIso = isoExcess(deIso, deRef)                                              !=== COMPUTE THE EXCESS
     1382  IF(PRESENT(delIso)) delIso = deIso
     1383  IF(PRESENT(excIso)) excIso = exIso
     1384  IF(PRESENT(delRef)) delRef = deRef
     1385  IF(.NOT.lc) RETURN
     1386  lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), SHAPE(exIso), mss, vNamEx(iIso), sub, nmax)
     1387END FUNCTION excessR_1D
     1388!-------------------------------------------------------------------------------------------------------------------------------
     1389LOGICAL FUNCTION excessR_2D(rIso, rRef, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr)
     1390  REAL, DIMENSION(:,:),                                 INTENT(IN)  :: rIso, rRef
     1391  INTEGER,                                              INTENT(IN)  :: iIso
     1392  CHARACTER(LEN=*),                           OPTIONAL, INTENT(IN)  :: err_tag, subn
     1393  INTEGER,                                    OPTIONAL, INTENT(IN)  :: nmax
     1394  REAL,                                       OPTIONAL, INTENT(IN)  :: delMax, excMax
     1395  REAL, DIMENSION(SIZE(rIso,1),SIZE(rIso,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
     1396  LOGICAL,                                    OPTIONAL, INTENT(IN)  :: lCheck
     1397  CHARACTER(LEN=maxlen) :: tag, sub, mss, m
     1398  REAL, DIMENSION(SIZE(rIso,1),SIZE(rIso,2)) :: deIso, deRef, exIso
     1399  REAL    :: dmx, drx, emn, emx
     1400  LOGICAL :: lc
     1401  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
     1402  sub = 'excess';            IF(PRESENT(subn   )) sub = TRIM(subn)//':'//TRIM(sub)  !--- Calling subroutine name
     1403  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1404  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
     1405  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
     1406  lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN                   !--- Ckeck the index + define few vars
     1407  emn = eMin(ixIso)%r(iIso)
     1408  drx = dMax(ixIso)%r(iso_ref)
     1409  mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
     1410  lerr = deltaR_2D(rIso, iIso,    tag, sub, nmax, dmx, deIso, lCheck) &        !=== COMPUTE deltaIso AND deltaRef
     1411    .OR. deltaR_2D(rRef, iso_ref, tag, sub, nmax, drx, deRef, lCheck)
     1412  IF(lerr)    RETURN
     1413  exIso = isoExcess(deIso, deRef)                                              !=== COMPUTE THE EXCESS
     1414  IF(PRESENT(delIso)) delIso = deIso
     1415  IF(PRESENT(excIso)) excIso = exIso
     1416  IF(PRESENT(delRef)) delRef = deRef
     1417  IF(.NOT.lc) RETURN
     1418  lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(PACK(deIso,.TRUE.), PACK(deRef,.TRUE.), PACK(exIso,.TRUE.)), &
     1419                     SHAPE(exIso), TRIM(mss)//' '//TRIM(m), vNamEx(iIso), sub, nmax)
     1420END FUNCTION excessR_2D
     1421!===============================================================================================================================
     1422LOGICAL FUNCTION excessQ_0D(qt, iIso, err_tag, subn, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck) RESULT(lerr)
     1423  REAL,                       INTENT(IN)  :: qt(:)
     1424  INTEGER,                    INTENT(IN)  :: iIso
     1425  CHARACTER(LEN=*), OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
     1426  REAL,             OPTIONAL, INTENT(IN)  :: delMax, excMax, qmin
     1427  REAL,             OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
     1428  LOGICAL,          OPTIONAL, INTENT(IN)  :: lCheck
     1429  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
     1430  REAL    :: deIso, exIso, deRef
     1431  REAL    :: dmx, drx, emn, emx, qmn
     1432  INTEGER :: iZon, iPha, ii
     1433  LOGICAL :: lc
     1434  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
     1435  sub = 'excess';            IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1436  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1437  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
     1438  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
     1439  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
     1440  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
     1441
     1442  !=== CHECK PHASE
     1443  lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
     1444  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
     1445  IF(lerr)       RETURN
     1446  emn = eMin(ixIso)%r(iIso)
     1447  drx = dMax(ixIso)%r(iso_ref)
     1448  m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
     1449
     1450  !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
     1451  ii = iIso
     1452  DO iZon = 0, nzone
     1453     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     1454     lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN                  !--- Check ii, set iso_ref + excess params
     1455     mss = TRIM(mss)//TRIM(m)
     1456     DO iPha = 1, nphas
     1457        p = isoPhas(iPha:iPha)                                                 !--- Current phase
     1458        lerr = deltaQ_0D(qt, ii,      tag, sub, dmx, deIso, p, qmn, lCheck) &  !=== COMPUTE deltaIso AND deltaRef
     1459          .OR. deltaQ_0D(qt, iso_ref, tag, sub, drx, deRef, p, qmn, lCheck)
     1460        exIso = isoExcess(deIso, deRef)                                        !=== COMPUTE THE EXCESS
     1461        IF(iZon == 0 .AND. p == pha) THEN
     1462           IF(PRESENT(delIso)) delIso = deIso                                  !--- Return delta(iIso,   phase=p)
     1463           IF(PRESENT(delRef)) delRef = deRef                                  !--- Return delta(iso_ref,phase=p)
     1464           IF(PRESENT(excIso)) excIso = exIso                                  !--- Return excess(iIso,  phase=p)
     1465        END IF
     1466        IF(.NOT.lc) CYCLE
     1467        lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), [1], mss, vNamEx(ii, iPha, iZon), sub)
     1468        IF(lerr) RETURN
     1469     END DO
     1470     IF(.NOT.lc) EXIT
     1471  END DO
     1472END FUNCTION excessQ_0D
     1473!-------------------------------------------------------------------------------------------------------------------------------
     1474LOGICAL FUNCTION excessQ_1D(qt, iIso, err_tag, subn, nmax, delMax,excMax,delIso,excIso,delRef, phas, qmin, lCheck) RESULT(lerr)
     1475  REAL,                                  INTENT(IN)  :: qt(:,:)
     1476  INTEGER,                               INTENT(IN)  :: iIso
     1477  CHARACTER(LEN=*),            OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
     1478  INTEGER,                     OPTIONAL, INTENT(IN)  :: nmax
     1479  REAL,                        OPTIONAL, INTENT(IN)  :: delMax, excMax, qmin
     1480  REAL, DIMENSION(SIZE(qt,1)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
     1481  LOGICAL,                     OPTIONAL, INTENT(IN)  :: lCheck
     1482  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
     1483  REAL, DIMENSION(SIZE(qt,1)) :: deIso, exIso, deRef
     1484  REAL    :: dmx, drx, emn, emx, qmn
     1485  INTEGER :: iZon, iPha, ii
     1486  LOGICAL :: lc
     1487  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
     1488  sub = 'excess';            IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1489  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1490  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
     1491  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
     1492  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
     1493  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
     1494
     1495  !=== CHECK PHASE
     1496  lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
     1497  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
     1498  IF(lerr)       RETURN
     1499  emn = eMin(ixIso)%r(iIso)
     1500  drx = dMax(ixIso)%r(iso_ref)
     1501  m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
     1502
     1503  !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
     1504  ii = iIso
     1505  DO iZon = 0, nzone
     1506     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     1507     lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN                  !--- Check ii, set iso_ref + excess params
     1508     mss = TRIM(mss)//TRIM(m)
     1509     DO iPha = 1, nphas
     1510        p = isoPhas(iPha:iPha)                                                 !--- Current phase
     1511        lerr = deltaQ_1D(qt, ii,      tag,sub,nmax, dmx, deIso,p,qmn,lCheck) & !=== COMPUTE deltaIso AND deltaRef
     1512          .OR. deltaQ_1D(qt, iso_ref, tag,sub,nmax, drx, deRef,p,qmn,lCheck)
     1513        exIso = isoExcess(deIso, deRef)                                        !=== COMPUTE THE EXCESS
     1514        IF(iZon == 0 .AND. p == pha) THEN
     1515           IF(PRESENT(delIso)) delIso = deIso                                  !--- Return delta(iIso,   phase=p)
     1516           IF(PRESENT(delRef)) delRef = deRef                                  !--- Return delta(iso_ref,phase=p)
     1517           IF(PRESENT(excIso)) excIso = exIso                                  !--- Return excess(iIso,  phase=p)
     1518        END IF
     1519        IF(.NOT.lc) CYCLE
     1520        lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), &
     1521                             SHAPE(exIso), mss, vNamEx(ii, iPha, iZon), sub, nmax)
     1522        IF(lerr) RETURN
     1523     END DO
     1524     IF(.NOT.lc) EXIT
     1525  END DO
     1526END FUNCTION excessQ_1D
     1527!-------------------------------------------------------------------------------------------------------------------------------
     1528LOGICAL FUNCTION excessQ_2D(qt, iIso, err_tag, subn, nmax, delMax,excMax,delIso,excIso,delRef, phas, qmin, lCheck) RESULT(lerr)
     1529  REAL,                                             INTENT(IN)  :: qt(:,:,:)
     1530  INTEGER,                                          INTENT(IN)  :: iIso
     1531  CHARACTER(LEN=*),                       OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
     1532  INTEGER,                                OPTIONAL, INTENT(IN)  :: nmax
     1533  REAL,                                   OPTIONAL, INTENT(IN)  :: delMax, excMax, qmin
     1534  REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
     1535  LOGICAL,                                OPTIONAL, INTENT(IN)  :: lCheck
     1536  CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m
     1537  REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)) :: deIso, exIso, deRef
     1538  REAL    :: dmx, drx, emn, emx, qmn
     1539  INTEGER :: iZon, iPha, ii
     1540  LOGICAL :: lc
     1541  tag = '';                  IF(PRESENT(err_tag)) tag = err_tag
     1542  sub = 'excess';            IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name
     1543  dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax                 !--- Maximum delta
     1544  emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax                 !--- Maximum excess
     1545  pha = 'g';                 IF(PRESENT(phas   )) pha = phas
     1546  qmn = small;               IF(PRESENT(qmin   )) qmn = qmin                   !--- q negligible if q<=dmn   (default: small)
     1547  lc  = .FALSE.;             IF(PRESENT(lCheck ))  lc = lCheck                 !--- Checking routines execution trigger
     1548
     1549  !=== CHECK PHASE
     1550  lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0
     1551  CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr)
     1552  IF(lerr)       RETURN
     1553  emn = eMin(ixIso)%r(iIso)
     1554  drx = dMax(ixIso)%r(iso_ref)
     1555  m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')'
     1556
     1557  !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone)
     1558  ii = iIso
     1559  DO iZon = 0, nzone
     1560     IF(iZon /= 0) ii = itZonIso(iZon, iIso)
     1561     lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN                  !--- Check ii, set iso_ref + excess params
     1562     mss = TRIM(mss)//TRIM(m)
     1563     DO iPha = 1, nphas
     1564        p = isoPhas(iPha:iPha)                                                 !--- Current phase
     1565        lerr = deltaQ_2D(qt, ii,      tag,sub,nmax, dmx, deIso,p,qmn,lCheck) & !=== COMPUTE deltaIso AND deltaRef
     1566          .OR. deltaQ_2D(qt, iso_ref, tag,sub,nmax, drx, deRef,p,qmn,lCheck)
     1567        exIso = isoExcess(deIso, deRef)                                        !=== COMPUTE THE EXCESS
     1568        IF(iZon == 0 .AND. p == pha) THEN
     1569           IF(PRESENT(delIso)) delIso = deIso                                  !--- Return delta(iIso,   phase=p)
     1570           IF(PRESENT(delRef)) delRef = deRef                                  !--- Return delta(iso_ref,phase=p)
     1571           IF(PRESENT(excIso)) excIso = exIso                                  !--- Return excess(iIso,  phase=p)
     1572        END IF
     1573        IF(.NOT.lc) CYCLE
     1574        lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(PACK(deIso,.TRUE.), PACK(deRef,.TRUE.), PACK(exIso,.TRUE.)), &
     1575                             SHAPE(exIso), mss, vNamEx(ii, iPha, iZon), sub, nmax)
     1576        IF(lerr) RETURN
     1577     END DO
     1578     IF(.NOT.lc) EXIT
     1579  END DO
     1580END FUNCTION excessQ_2D
     1581!===============================================================================================================================
     1582LOGICAL FUNCTION excessQx_2D(iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck) RESULT(lerr)
     1583  INTEGER,                                          INTENT(IN)  :: iIso
     1584  CHARACTER(LEN=*),                       OPTIONAL, INTENT(IN)  :: err_tag, subn, phas
     1585  INTEGER,                                OPTIONAL, INTENT(IN)  :: nmax
     1586  REAL,                                   OPTIONAL, INTENT(IN)  :: delMax, excMax, qmin
     1587  REAL, DIMENSION(SIZE(qx,1),SIZE(qx,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef
     1588  LOGICAL,                                OPTIONAL, INTENT(IN)  :: lCheck
     1589  lerr = excessQ_2D(qx, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck)
     1590END FUNCTION excessQx_2D
     1591!===============================================================================================================================
     1592
     1593
     1594!===============================================================================================================================
     1595!=== GENERAL PURPOSE CHECKING ROUTINE:   lerr = checkTrac(err_tag[, subn][, nmax][, sCheck])
     1596!===
     1597!=== ARGUMENTS:                                                                     Intent  Optional?  Default value
     1598!===   * err_tag   Error tag sent from the calling routine                          INPUT   OPTIONAL    none
     1599!===   * subn      Calling subroutine name                                          INPUT   OPTIONAL   "checkTrac"
     1600!===   * nmax      Maximum number of lines to print for outliers                    INPUT   OPTIONAL    32
     1601!===   * sCheck    Triggers for verifications                                       INPUT   OPTIONAL    sIsoChech
     1602!===
     1603!=== REMARKS:
     1604!===   Tunable thresholds available in the individual routines (delMax, excMax, qmin, tny, absErr, relErr) have not been kept
     1605!===   as optional arguments in this routine, because the adequate values are tracer or isotope-dependent.
     1606!===   For special usages, a specific version can be written, or individual routines with special thresholds can be called.
     1607!===============================================================================================================================
     1608LOGICAL FUNCTION checkTrac(err_tag, subn, nmax, sCheck) RESULT(lerr)
     1609  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, sCheck
     1610  INTEGER,          OPTIONAL, INTENT(IN) :: nmax
     1611  INTEGER :: nmx, i, iPha, iIso, iq
     1612  LOGICAL :: lNan, lPos, lMas, lDel, lExc, lBnd
     1613  CHARACTER(LEN=maxlen) :: tag, sub, chk
     1614  CHARACTER(LEN=maxlen), ALLOCATABLE :: tnam(:)
     1615
     1616!!!!! GERER TNAM
     1617
     1618  IF(niso == 0) RETURN
     1619  tag = '';                 IF(PRESENT(err_tag)) tag = err_tag
     1620  sub = 'checkTrac';        IF(PRESENT(subn   )) sub = TRIM(sub)//TRIM(subn)
     1621  nmx = 32;                 IF(PRESENT(nmax   )) nmx = nmax
     1622  chk =  sIsoCheck(ixIso) ; IF(PRESENT(sCheck )) chk = sCheck
     1623  SELECT CASE(isotope%parent)
     1624     CASE('H2O')
     1625        DO iPha = 1, isotope%nphas
     1626           DO iIso = 1, niso
     1627              iq = iqIsoPha(iIso,iPha); tnam = tracers(iq)%name!; qmin = tracers(iq)%qmin; qmax = tracers(iq)%qmax
     1628              IF(chk(1:1) == 'n') lerr = checkNaN (qx(:,:,iq),tag,tnam,sub,nmx); IF(lerr) RETURN
     1629              IF(chk(2:2) == 'p') lerr = checkPos (qx(:,:,iq),tag,tnam,sub,nmx); IF(lerr) RETURN ! tny
     1630              IF(chk(3:3) == 'm') lerr = checkMass(           tag,sub,     nmx); IF(lerr) RETURN ! absErr, relErr
     1631             !IF(chk(6:6) == 'b') lerr = checkBnd (qx(:,:,iq),tag,tnam,sub,nmx,qmin,qmax); IF(lerr) RETURN
     1632           END DO
     1633           iIso = iso_HDO; IF(chk(4:4) == 'd') lerr =  delta(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN
     1634           iIso = iso_O18; IF(chk(4:4) == 'd') lerr =  delta(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN
     1635           iIso = iso_HDO; IF(chk(5:5) == 'e') lerr = excess(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN
     1636           iIso = iso_O17; IF(chk(5:5) == 'e') lerr = excess(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN
     1637        END DO
     1638  END SELECT
     1639END FUNCTION checkTrac
     1640!===============================================================================================================================
     1641
     1642
     1643 
    941644        SUBROUTINE iso_verif_init()
    951645        use ioipsl_getin_p_mod, ONLY : getin_p
Note: See TracChangeset for help on using the changeset viewer.