source: LMDZ6/trunk/libf/phylmdiso/isotopes_verif_mod.F90 @ 5802

Last change on this file since 5802 was 5774, checked in by dcugnet, 5 months ago

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

  • Property svn:executable set to *
File size: 182.3 KB
Line 
1
2#ifdef ISOVERIF
3! $Id: $
4
5MODULE isotopes_verif_mod
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              !---  "  "
148
149! variables de vérifications
150real errmax ! erreur maximale en absolu.
151parameter (errmax=abs_err)
152real errmax_sol ! erreur maximale en absolu.
153parameter (errmax_sol=5e-7)
154      real errmaxrel ! erreur maximale en relatif autorisée
155!      parameter (errmaxrel=1e10)
156      parameter (errmaxrel=rel_err)
157      real borne ! valeur maximale que n'importe quelle variable peut
158                 ! atteindre, en val abs; utile pour vérif que pas NAN
159      parameter (borne=max_val)
160      real, save ::  deltalim ! deltalim est le maximum de deltaD qu'on
161                             ! autorise dans la vapeur, au-dela, en suppose que c'est abérrant.
162                             ! dans le liquide, on autorise deltalim*faccond.
163!$OMP THREADPRIVATE(deltalim)
164!      parameter (deltalim=1e10)
165!      parameter (deltalim=300.0)
166       ! maintenant défini dans iso.def
167
168       real, save :: deltalimtrac ! max de deltaD dans les traceurs, défini dans iso.def
169!$OMP THREADPRIVATE(deltalimtrac)
170
171      real, save ::  deltalim_snow ! deltalim est le maximum de deltaD qu'on
172                             ! autorise dans la neige, au-dela, en suppose que c'est abérrant.
173!$OMP THREADPRIVATE(deltalim_snow)
174!      parameter (deltalim_snow=500.0) ! initalisé dans iso_init
175
176        real, save ::  deltaDmin
177!$OMP THREADPRIVATE(deltaDmin)
178!       parameter (deltaDmin=-900.0)
179    ! maintentant, défini dans iso.def
180
181      real, save ::  o17excess_bas,o17excess_haut ! bornes inf et sup de l'O17-excess
182!      parameter(o17excess_bas=-200.0,o17excess_haut=120)
183!$OMP THREADPRIVATE(o17excess_bas,o17excess_haut)
184      integer nlevmaxO17
185!$OMP THREADPRIVATE(nlevmaxO17)
186
187      logical, save ::  O17_verif
188!$OMP THREADPRIVATE(O17_verif)
189!      parameter (O17_verif=.true.)
190
191        real, save ::  dexcess_min,dexcess_max
192!$OMP THREADPRIVATE(dexcess_min,dexcess_max)
193
194!      real faccond  ! dans le liquide, on autorise R(deltalim)*faccond.
195!      parameter (faccond=1.1)
196     
197!      logical bidouille_anti_divergence ! si true, alors on fait un
198!                                        ! rappel régulier des xt4 vers les q pour
199!                                        !éviter accumulation d'érreurs et  divergences
200!      parameter (bidouille_anti_divergence=.true.)
201!      parameter (bidouille_anti_divergence=.false.)
202    ! bidouille_anti_divergence a été déplacé dans wateriso2.h et est lue à l'execution
203
204       
205        real deltaDfaible ! deltaD en dessous duquel la vapeur est tellement faible
206        !que on peut accepter que la remise à l'équilibre du sol avec
207        ! cette vapeur donne des deltaDevap aberrants.
208        parameter (deltaDfaible=low_delD)
209
210        real faible_evap ! faible évaporation: on est plus laxiste
211        !pour les deltaD aberrant dans le cas de l'évap venant d'orchidee
212        parameter (faible_evap=min_delDevap)
213
214        real Tmin_verif
215        parameter (Tmin_verif=min_T) ! en K
216        real Tmax_verif
217        parameter (Tmax_verif=max_T) ! en K
218
219
220! subroutines de vérifications génériques, à ne pas modifier
221
222 
223CONTAINS
224
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 
1644        SUBROUTINE iso_verif_init()
1645        use ioipsl_getin_p_mod, ONLY : getin_p
1646        use isotopes_mod, ONLY: iso_O17, iso_O18, iso_HDO
1647        implicit none
1648
1649        write(*,*) 'iso_verif_init 99: entree'
1650        deltalim=300.0
1651        deltalim_snow=500.0
1652        deltaDmin=-900.0
1653        deltalimtrac=2000.0
1654        O17_verif=.false.
1655        o17excess_bas=-200.0
1656        o17excess_haut=120.0
1657        dexcess_min=-100.0
1658        dexcess_max=1000.0
1659
1660call getin_p('deltalim',deltalim)
1661deltalim_snow=deltalim+200.0
1662call getin_p('deltaDmin',deltaDmin)
1663call getin_p('deltalimtrac',deltalimtrac)
1664
1665write(*,*) 'iso_O17,iso_O18,iso_HDO=',iso_O17,iso_O18,iso_HDO
1666if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
1667  call getin_p('O17_verif',O17_verif) 
1668  call getin_p('o17excess_bas',o17excess_bas) 
1669  call getin_p('o17excess_haut',o17excess_haut) 
1670  call getin_p('nlevmaxO17',nlevmaxO17) 
1671endif
1672if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
1673  call getin_p('dexcess_min',dexcess_min) 
1674  call getin_p('dexcess_max',dexcess_max) 
1675endif
1676
1677write(*,*) 'deltalim=',deltalim
1678write(*,*) 'deltaDmin=',deltaDmin     
1679write(*,*) 'deltalimtrac=',deltalimtrac 
1680write(*,*) 'O17_verif=',O17_verif
1681write(*,*) 'o17excess_bas=',o17excess_bas
1682write(*,*) 'o17excess_haut=',o17excess_haut 
1683write(*,*) 'dexcess_min=',dexcess_min
1684write(*,*) 'dexcess_max=',dexcess_max
1685
1686        end SUBROUTINE iso_verif_init
1687
1688      subroutine iso_verif_egalite(a,b,err_msg)
1689        implicit none
1690        ! compare a et b. Si pas egal à errmax près, on affiche message
1691        ! d'erreur et stoppe
1692
1693        ! input:
1694        real a, b
1695        character*(*) err_msg ! message d''erreur à afficher
1696
1697        ! local
1698        !integer iso_verif_egalite_choix_nostop
1699
1700
1701        if (iso_verif_egalite_choix_nostop(a,b,err_msg, &
1702     &           errmax,errmaxrel).eq.1) then
1703                stop
1704        endif
1705       
1706#ifdef ISOVERIF
1707#else
1708        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
1709        stop
1710#endif           
1711
1712        return
1713        end subroutine iso_verif_egalite
1714
1715        !*****************
1716
1717        function iso_verif_egalite_nostop(a,b,err_msg)
1718        implicit none
1719        ! compare a et b. Si pas egal à errmax près, on affiche message
1720        ! d'erreur et stoppe
1721
1722        ! input:
1723        real a, b
1724        character*(*) err_msg ! message d''erreur à afficher
1725        !ouptut
1726        integer iso_verif_egalite_nostop
1727        ! local
1728        !integer iso_verif_egalite_choix_nostop
1729
1730
1731        iso_verif_egalite_nostop=iso_verif_egalite_choix_nostop &
1732     &          (a,b,err_msg,errmax,errmaxrel)
1733
1734#ifdef ISOVERIF
1735#else
1736        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
1737        stop
1738#endif                 
1739
1740        return
1741        end function iso_verif_egalite_nostop
1742
1743
1744
1745        !************************************
1746
1747        subroutine iso_verif_aberrant(R,err_msg)
1748        use isotopes_mod, ONLY: ridicule, iso_HDO
1749        implicit none
1750        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
1751        ! d'erreur et stoppe
1752
1753        ! input:
1754        real R
1755        character*(*) err_msg ! message d''erreur à afficher
1756        !real deltaD
1757        !integer iso_verif_aberrant_choix_nostop
1758
1759
1760        if (iso_verif_aberrant_choix_nostop(R,1.0,ridicule, &
1761     &           deltalim,err_msg).eq.1) then
1762             stop
1763        endif
1764
1765#ifdef ISOVERIF
1766        if (.not.(iso_HDO.gt.0)) then   
1767         write(*,*) 'iso86: err_msg=',err_msg,': pourquoi verif?'
1768          stop
1769        endif
1770#else       
1771        write(*,*) 'iso 90: err_msg=',err_msg,': pourquoi verif?'
1772        stop
1773#endif                   
1774               
1775        end subroutine iso_verif_aberrant
1776
1777        subroutine iso_verif_aberrant_encadre(R,err_msg)
1778        use isotopes_mod, ONLY: ridicule, iso_HDO
1779        implicit none
1780        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
1781        ! d'erreur et stoppe
1782
1783        ! input:
1784        real R
1785        character*(*) err_msg ! message d''erreur à afficher
1786        !real deltaD
1787
1788        !integer iso_verif_aberrant_enc_choix_nostop
1789
1790
1791        if (iso_verif_aberrant_enc_choix_nostop(R,1.0,ridicule, &
1792     &           deltalim,err_msg).eq.1) then
1793             write(*,*) 'deltaD=',deltaD(R)
1794             call abort_physic('isotopes_verif_mod > iso_verif_aberrant_encadre',err_msg,1)
1795             !stop             
1796        endif
1797
1798#ifdef ISOVERIF
1799        if (.not.(iso_HDO.gt.0)) then   
1800         write(*,*) 'iso86: err_msg=',err_msg,': pourquoi verif?'
1801          stop
1802        endif
1803#else       
1804        write(*,*) 'iso 90: err_msg=',err_msg,': pourquoi verif?'
1805        stop
1806#endif                   
1807       
1808        end subroutine iso_verif_aberrant_encadre
1809
1810        !************************************
1811
1812        subroutine iso_verif_aberrant_choix(xt,q,qmin,deltaDmax,err_msg)
1813        use isotopes_mod, ONLY: iso_HDO
1814        implicit none
1815        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
1816        ! d'erreur et stoppe
1817
1818        ! input:
1819        real xt,q,qmin,deltaDmax
1820        character*(*) err_msg ! message d''erreur à afficher
1821        !real deltaD
1822
1823        ! locals
1824        !integer iso_verif_aberrant_choix_nostop
1825
1826
1827        if (iso_verif_aberrant_choix_nostop(xt,q,qmin, &
1828     &           deltaDmax,err_msg).eq.1) then
1829             stop
1830        endif
1831
1832#ifdef ISOVERIF
1833        if (.not.(iso_HDO.gt.0)) then   
1834         write(*,*) 'iso122: err_msg=',err_msg,': pourquoi verif?'
1835          stop
1836        endif
1837#else
1838        write(*,*) 'iso126: err_msg=',err_msg,': pourquoi verif?'
1839        stop
1840#endif                   
1841       
1842        end subroutine iso_verif_aberrant_choix
1843
1844        !************************************
1845
1846        function iso_verif_aberrant_nostop(R,err_msg)
1847        use isotopes_mod, ONLY: ridicule,iso_HDO
1848        implicit none
1849        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
1850        ! d'erreur et stoppe
1851
1852        ! input:
1853        real R
1854        character*(*) err_msg ! message d''erreur à afficher
1855        integer iso_verif_aberrant_nostop ! output: 1 si erreur, 0 sinon
1856        !real deltaD
1857
1858        ! locals
1859        !integer iso_verif_aberrant_choix_nostop
1860
1861        iso_verif_aberrant_nostop=iso_verif_aberrant_choix_nostop &
1862     &           (R,1.0,ridicule,deltalim,err_msg)
1863
1864#ifdef ISOVERIF
1865        if (.not.(iso_HDO.gt.0)) then   
1866         write(*,*) 'iso156: err_msg=',err_msg,': pourquoi verif?'
1867          stop
1868        endif
1869#else
1870        write(*,*) 'iso160: err_msg=',err_msg,': pourquoi verif?'
1871        stop
1872#endif           
1873       
1874        return
1875        end function iso_verif_aberrant_nostop
1876
1877        function iso_verif_aberrant_enc_nostop(R,err_msg)
1878        use isotopes_mod, ONLY: ridicule,iso_HDO
1879        implicit none
1880        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
1881        ! d'erreur et stoppe
1882
1883        ! input:
1884        real R
1885        character*(*) err_msg ! message d''erreur à afficher
1886        integer iso_verif_aberrant_enc_nostop ! output: 1 si erreur, 0 sinon
1887        !real deltaD
1888
1889        ! locals
1890        !integer iso_verif_aberrant_enc_choix_nostop
1891
1892        iso_verif_aberrant_enc_nostop= &
1893     &           iso_verif_aberrant_enc_choix_nostop &
1894     &           (R,1.0,ridicule,deltalim,err_msg)
1895
1896#ifdef ISOVERIF
1897        if (.not.(iso_HDO.gt.0)) then   
1898         write(*,*) 'iso156: err_msg=',err_msg,': pourquoi verif?'
1899          stop
1900        endif
1901#else
1902        write(*,*) 'iso160: err_msg=',err_msg,': pourquoi verif?'
1903        stop
1904#endif                   
1905        return
1906        end function iso_verif_aberrant_enc_nostop
1907
1908        !************************************
1909
1910        function iso_verif_aberrant_choix_nostop(xt,q,   &
1911     &            qmin,deltaDmax,err_msg)
1912
1913        use isotopes_mod, ONLY: iso_HDO
1914        implicit none
1915        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
1916        ! d'erreur et stoppe
1917
1918        ! input:
1919        real xt,q,qmin,deltaDmax
1920        character*(*) err_msg ! message d''erreur à afficher
1921        ! output
1922        integer iso_verif_aberrant_choix_nostop
1923        ! locals
1924        !real deltaD
1925        !integer iso_verif_noNaN_nostop       
1926
1927
1928        iso_verif_aberrant_choix_nostop=0
1929
1930#ifdef ISOVERIF       
1931        if (iso_verif_noNaN_nostop(q,err_msg).eq.1) then
1932            write(*,*) 'q=',q
1933            iso_verif_aberrant_choix_nostop=1
1934        endif     
1935        if (iso_verif_noNaN_nostop(xt,err_msg).eq.1) then
1936            write(*,*) 'xt=',xt
1937            iso_verif_aberrant_choix_nostop=1
1938        endif
1939#endif
1940
1941        if (q.gt.qmin) then
1942        if ((deltaD(xt/q).gt.deltaDmax).or. &
1943     &          (deltaD(xt/q).lt.-borne).or. &
1944     &          (xt.lt.-borne).or. &
1945     &          (xt.gt.borne)) then                 
1946            write(*,*) 'erreur detectee par '// &
1947     &             'iso_verif_aberrant_choix_nostop:'
1948            write(*,*) err_msg
1949            write(*,*) 'q,deltaD=',q,deltaD(xt/q)
1950            write(*,*) 'deltaDmax=',deltaDmax
1951            iso_verif_aberrant_choix_nostop=1
1952            if (abs(xt-q)/q.lt.errmax) then
1953                write(*,*) 'attention, n''a-t-on pas confondu'// &
1954     &           ' iso_HDO et iso_eau dans l''appel à la verif?'
1955            endif
1956        endif
1957        endif
1958
1959#ifdef ISOVERIF
1960        if (.not.(iso_HDO.gt.0)) then   
1961         write(*,*) 'iso205: err_msg=',err_msg,': pourquoi verif?'
1962          stop
1963        endif
1964#else
1965        write(*,*) 'iso 209: err_msg=',err_msg,': pourquoi verif?'
1966        stop
1967#endif                   
1968       
1969        return
1970        end function iso_verif_aberrant_choix_nostop
1971
1972        function iso_verif_aberrant_enc_choix_nostop(xt,q,   &
1973     &            qmin,deltaDmax,err_msg)
1974        use isotopes_mod, ONLY: iso_HDO
1975        implicit none
1976        ! si le rapprot iso R est plus grand que deltaDlim, on affiche message
1977        ! d'erreur et stoppe
1978
1979        ! input:
1980        real xt,q,qmin,deltaDmax
1981        character*(*) err_msg ! message d''erreur à afficher
1982        ! output
1983        integer iso_verif_aberrant_enc_choix_nostop
1984        ! locals
1985        !real deltaD
1986
1987        iso_verif_aberrant_enc_choix_nostop=0
1988        if (q.gt.qmin) then
1989        if ((deltaD(xt/q).gt.deltaDmax).or. &
1990     &          (deltaD(xt/q).lt.deltaDmin)) then                 
1991            write(*,*) 'erreur detectee par '// &
1992     &             'iso_verif_aberrant_choix_nostop:'
1993            write(*,*) err_msg
1994            write(*,*) 'q,deltaD=',q,deltaD(xt/q)
1995            iso_verif_aberrant_enc_choix_nostop=1
1996            if (abs(xt-q)/q.lt.errmax) then
1997                write(*,*) 'attention, n''a-t-on pas confondu'// &
1998     &           ' iso_HDO et iso_eau dans l''appel à la verif?'
1999            endif
2000        endif
2001        endif
2002
2003#ifdef ISOVERIF
2004        if (.not.(iso_HDO.gt.0)) then   
2005         write(*,*) 'iso205: err_msg=',err_msg,': pourquoi verif?'
2006          stop
2007        endif
2008#else
2009        write(*,*) 'iso 209: err_msg=',err_msg,': pourquoi verif?'
2010        stop
2011#endif                   
2012       
2013        return
2014        end function iso_verif_aberrant_enc_choix_nostop
2015
2016        !*******************
2017
2018        subroutine iso_verif_aberrant_o17(R17,R18,err_msg)
2019        implicit none
2020        ! si l'O17-excess est aberrant, on affiche un message
2021        ! d'erreur et stoppe
2022
2023        ! input:
2024        real R17,R18
2025        character*(*) err_msg ! message d''erreur à afficher
2026        !real o17excess
2027
2028        ! locals
2029        !integer iso_verif_aberrant_o17_nostop
2030
2031!        write(*,*) 'O17_verif=',O17_verif
2032        if (O17_verif) then
2033            if (iso_verif_aberrant_o17_nostop(R17,R18,err_msg) &
2034     &           .eq.1) then
2035                stop
2036            endif
2037        endif !if (O17_verif) then
2038
2039#ifdef ISOVERIF
2040#else
2041        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
2042        stop
2043#endif                   
2044       
2045        return
2046        end subroutine iso_verif_aberrant_o17
2047
2048        !*******************
2049
2050        function iso_verif_aberrant_o17_nostop(R17,R18,err_msg)
2051        USE isotopes_mod, ONLY: tnat,iso_O17,iso_O18
2052        implicit none
2053        ! si l'O17-excess est aberrant, on affiche un message
2054        ! d'erreur et renvoit 1
2055
2056        ! input:
2057        real R17,R18
2058        character*(*) err_msg ! message d''erreur à afficher
2059        !local
2060        !real o17excess
2061        ! output
2062        integer iso_verif_aberrant_o17_nostop
2063
2064        if (O17_verif) then
2065
2066        iso_verif_aberrant_o17_nostop=0
2067        if ((o17excess(R17,R18).gt.o17excess_haut).or. &
2068     &            (o17excess(R17,R18).lt.o17excess_bas)) then 
2069            write(*,*) 'erreur detectee par iso_verif_aberrant_O17:'
2070            write(*,*) err_msg
2071            write(*,*) 'o17excess=',o17excess(R17,R18)
2072            write(*,*) 'deltaO17=',(R17/tnat(iso_o17)-1.0)*1000.0
2073            write(*,*) 'deltaO18=',(R18/tnat(iso_O18)-1.0)*1000.0
2074            ! attention, vérifier que la ligne suivante est bien activée
2075            iso_verif_aberrant_o17_nostop=1
2076        endif
2077
2078        endif  !if (O17_verif) then
2079
2080#ifdef ISOVERIF
2081#else
2082        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
2083        stop
2084#endif                   
2085       
2086        return
2087        end function iso_verif_aberrant_o17_nostop
2088
2089
2090        !************************************
2091
2092        subroutine iso_verif_noNaN(x,err_msg)
2093        implicit none
2094        ! si x est NaN, on affiche message
2095        ! d'erreur et stoppe
2096
2097        ! input:
2098        real x
2099        character*(*) err_msg ! message d''erreur à afficher
2100
2101        ! locals
2102        !integer iso_verif_noNAN_nostop
2103
2104        if (iso_verif_noNAN_nostop(x,err_msg).eq.1) then
2105            stop
2106        endif
2107#ifdef ISOVERIF
2108#else
2109        write(*,*) 'err_msg iso443=',err_msg,': pourquoi verif?'
2110        stop
2111#endif           
2112       
2113        end subroutine iso_verif_noNaN
2114
2115                !************************************
2116
2117        function iso_verif_noNaN_nostop(x,err_msg)
2118        implicit none
2119        ! si x est NaN, on affiche message
2120        ! d'erreur et return 1 si erreur
2121
2122        ! input:
2123        real x
2124        character*(*) err_msg ! message d''erreur à afficher
2125
2126        ! output
2127        integer iso_verif_noNAN_nostop
2128
2129        if ((x.gt.-borne).and.(x.lt.borne)) then
2130                iso_verif_noNAN_nostop=0
2131        else
2132            write(*,*) 'erreur detectee par iso_verif_nonNAN:'
2133            write(*,*) err_msg
2134            write(*,*) 'x=',x
2135            iso_verif_noNAN_nostop=1
2136        endif     
2137
2138#ifdef ISOVERIF
2139#else
2140        write(*,*) 'err_msg iso 482=',err_msg,': pourquoi verif?'
2141        stop
2142#endif           
2143
2144        return
2145        end function iso_verif_noNaN_nostop
2146
2147        subroutine iso_verif_noNaN_vect2D(x,err_msg,ni,n,m)
2148        implicit none
2149        ! si x est NaN, on affiche message
2150        ! d'erreur et return 1 si erreur
2151
2152        ! input:
2153          integer n,m,ni
2154        real x(ni,n,m)
2155        character*(*) err_msg ! message d''erreur à afficher
2156
2157        ! output
2158
2159        ! locals       
2160        integer i,j,ixt
2161
2162      do i=1,n
2163       do j=1,m
2164        do ixt=1,ni
2165         if ((x(ixt,i,j).gt.-borne).and. &
2166     &            (x(ixt,i,j).lt.borne)) then
2167         else !if ((x(ixt,i,j).gt.-borne).and.
2168            write(*,*) 'erreur detectee par iso_verif_nonNAN:'
2169            write(*,*) err_msg
2170            write(*,*) 'x,ixt,i,j=',x(ixt,i,j),ixt,i,j
2171            stop
2172         endif  !if ((x(ixt,i,j).gt.-borne).and.
2173        enddo !do ixt=1,ni
2174       enddo !do j=1,m
2175      enddo !do i=1,n     
2176
2177#ifdef ISOVERIF
2178#else
2179        write(*,*) 'err_msg iso525=',err_msg,': pourquoi verif?'
2180        stop
2181#endif           
2182
2183        end subroutine iso_verif_noNaN_vect2D
2184
2185
2186        subroutine iso_verif_noNaN_vect1D(x,err_msg,ni,n)
2187        implicit none
2188        ! si x est NaN, on affiche message
2189        ! d'erreur et return 1 si erreur
2190
2191        ! input:
2192          integer n,ni
2193        real x(ni,n)
2194        character*(*) err_msg ! message d''erreur à afficher
2195
2196        ! output
2197
2198        ! locals       
2199        integer i,ixt
2200
2201      do i=1,n
2202        do ixt=1,ni
2203         if ((x(ixt,i).gt.-borne).and. &
2204     &            (x(ixt,i).lt.borne)) then
2205         else !if ((x(ixt,i,j).gt.-borne).and.
2206            write(*,*) 'erreur detectee par iso_verif_nonNAN:'
2207            write(*,*) err_msg
2208            write(*,*) 'x,ixt,i=',x(ixt,i),ixt,i
2209            stop
2210         endif  !if ((x(ixt,i,j).gt.-borne).and.
2211        enddo !do ixt=1,ni
2212      enddo !do i=1,n     
2213
2214#ifdef ISOVERIF
2215#else
2216        write(*,*) 'err_msg iso525=',err_msg,': pourquoi verif?'
2217        stop
2218#endif           
2219
2220        end subroutine iso_verif_noNaN_vect1D
2221
2222
2223
2224        !************************
2225        subroutine iso_verif_egalite_choix(a,b,err_msg,erabs,errel)
2226        implicit none
2227        ! compare a et b. Si pas egal, on affiche message
2228        ! d'erreur et stoppe
2229        ! pour egalite, on verifie erreur absolue et arreur relative
2230
2231        ! input:
2232        real a, b
2233        real erabs,errel !erreur absolue et relative
2234        character*(*) err_msg ! message d''erreur à afficher
2235
2236        ! locals
2237        !integer iso_verif_egalite_choix_nostop
2238
2239        if (iso_verif_egalite_choix_nostop( &
2240     &            a,b,err_msg,erabs,errel).eq.1) then
2241           stop
2242        endif
2243
2244#ifdef ISOVERIF
2245#else
2246        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
2247        stop
2248#endif           
2249
2250        end subroutine iso_verif_egalite_choix
2251
2252!************************
2253        function iso_verif_egalite_choix_nostop &
2254     &           (a,b,err_msg,erabs,errel)
2255        implicit none
2256        ! compare a et b. Si pas egal, on affiche message
2257        ! d'erreur et stoppe
2258        ! pour egalite, on verifie erreur absolue et arreur relative
2259
2260        ! input:
2261        real a, b
2262        real erabs,errel !erreur absolue et relative
2263        character*(*) err_msg ! message d''erreur à afficher
2264
2265        ! output
2266        integer iso_verif_egalite_choix_nostop
2267        ! locals
2268        !integer iso_verif_noNaN_nostop
2269
2270        iso_verif_egalite_choix_nostop=0
2271
2272#ifdef ISOVERIF
2273        if (iso_verif_noNaN_nostop(a,err_msg).eq.1) then
2274            write(*,*) 'a=',a
2275            iso_verif_egalite_choix_nostop=1
2276        endif     
2277        if (iso_verif_noNaN_nostop(b,err_msg).eq.1) then
2278            write(*,*) 'b=',b
2279            iso_verif_egalite_choix_nostop=1
2280        endif     
2281#endif
2282
2283        if (abs(a-b).gt.erabs) then
2284          if (abs((a-b)/max(max(abs(b),abs(a)),1e-18)) &
2285     &            .gt.errel) then
2286            write(*,*) 'erreur detectee par iso_verif_egalite:'
2287            write(*,*) err_msg
2288            write(*,*) 'a=',a
2289            write(*,*) 'b=',b
2290            write(*,*) 'erabs,errel=',erabs,errel
2291            iso_verif_egalite_choix_nostop=1
2292!            stop
2293          endif
2294        endif
2295
2296#ifdef ISOVERIF
2297#else
2298        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
2299        stop
2300#endif           
2301       
2302        return
2303        end function iso_verif_egalite_choix_nostop       
2304
2305
2306
2307        !******************
2308        subroutine iso_verif_positif(x,err_msg)
2309        use isotopes_mod, ONLY: ridicule
2310        implicit none
2311        ! si x<0, on plante.
2312        ! si très limite, on le met à 0.
2313
2314        ! input:
2315        real x
2316        character*(*) err_msg ! message d''erreur à afficher
2317
2318        ! locals
2319        !integer iso_verif_positif_choix_nostop
2320
2321        if (iso_verif_positif_choix_nostop(x,ridicule,err_msg) &
2322     &           .eq.1) then
2323           stop
2324        endif
2325     
2326#ifdef ISOVERIF
2327#else
2328        write(*,*) 'iso_verif 637: err_msg=',err_msg, &
2329     &           ': pourquoi verif?'
2330        stop
2331#endif
2332       
2333        end subroutine iso_verif_positif
2334
2335        !******************
2336        subroutine iso_verif_positif_vect(x,n,err_msg)
2337        use isotopes_mod, ONLY: ridicule
2338        implicit none
2339        ! si x<0, on plante.
2340
2341        ! input:
2342        integer n
2343        real x(n)
2344        character*(*) err_msg ! message d''erreur à afficher
2345
2346        ! locals
2347        integer i
2348        integer iso_verif_positif_nostop
2349        integer ifaux
2350
2351        iso_verif_positif_nostop=0
2352        do i=1,n
2353          if (x(i).lt.-ridicule) then
2354            iso_verif_positif_nostop=1
2355            ifaux=i
2356          endif
2357        enddo
2358        if (iso_verif_positif_nostop.eq.1) then
2359          write(*,*) 'erreur detectee par iso_verif_positif_vect:'
2360          write(*,*) err_msg
2361          write(*,*) 'i,x=',ifaux,x(ifaux)
2362          stop
2363        endif   
2364       
2365        end subroutine iso_verif_positif_vect
2366
2367        subroutine iso_verif_positif_choix_vect(x,n,ridic,err_msg)
2368        implicit none
2369        ! si x<0, on plante.
2370
2371        ! input:
2372        integer n
2373        real x(n)
2374        real ridic
2375        character*(*) err_msg ! message d''erreur à afficher
2376        ! locals
2377        integer i
2378        integer iso_verif_positif_nostop
2379        integer ifaux
2380
2381        iso_verif_positif_nostop=0
2382        do i=1,n
2383          if (x(i).lt.-ridic) then
2384                iso_verif_positif_nostop=1
2385                ifaux=i
2386          endif
2387        enddo
2388        if (iso_verif_positif_nostop.eq.1) then                       
2389         write(*,*) 'erreur detectee par iso_verif_positif_choix_vect:'
2390         write(*,*) err_msg
2391         write(*,*) 'i,x=',ifaux,x(ifaux)
2392         stop
2393        endif   
2394       
2395        end subroutine iso_verif_positif_choix_vect
2396
2397
2398!******************
2399        subroutine iso_verif_positif_strict(x,err_msg)
2400        implicit none
2401        ! si x<0, on plante.
2402        ! si très limite, on le met à 0.
2403
2404        ! input:
2405        real x
2406        character*(*) err_msg ! message d''erreur à afficher
2407
2408        ! locals
2409        !integer iso_verif_positif_strict_nostop
2410
2411        if (iso_verif_positif_strict_nostop(x,err_msg) &
2412     &           .eq.1) then
2413           stop
2414        endif           
2415       
2416        end subroutine iso_verif_positif_strict
2417
2418!******************
2419
2420        function iso_verif_positif_strict_nostop(x,err_msg)
2421        implicit none
2422        ! si x<0, on plante.
2423        ! si très limite, on le met à 0.
2424
2425        ! input:
2426        real x
2427        character*(*) err_msg ! message d''erreur à afficher*
2428
2429        ! output
2430        integer iso_verif_positif_strict_nostop
2431
2432        if (x.gt.0.0) then
2433            iso_verif_positif_strict_nostop=0
2434        else     
2435            write(*,*) 'erreur detectee par iso_verif_positif_strict:'
2436            write(*,*) err_msg
2437            write(*,*) 'x=',x 
2438            iso_verif_positif_strict_nostop=1   
2439!            stop 
2440        endif   
2441       
2442        return
2443        end function iso_verif_positif_strict_nostop
2444
2445!******************
2446
2447        subroutine iso_verif_positif_choix(x,ridic,err_msg)
2448        implicit none
2449        ! si x<0, on plante.
2450        ! si très limite, on le met à 0.
2451
2452        ! input:
2453        real x
2454        real ridic
2455        character*(*) err_msg ! message d''erreur à afficher
2456
2457        ! locals
2458        !integer iso_verif_positif_choix_nostop
2459
2460        if (iso_verif_positif_choix_nostop(x,ridic,err_msg) &
2461     &           .eq.1) then
2462           stop
2463        endif
2464     
2465#ifdef ISOVERIF
2466#else
2467        write(*,*) 'iso_verif 801: err_msg=',err_msg, &
2468     &           ': pourquoi verif?'
2469        stop
2470#endif           
2471       
2472        end subroutine iso_verif_positif_choix       
2473
2474        !******************
2475        function iso_verif_positif_nostop(x,err_msg)
2476        use isotopes_mod, ONLY: ridicule
2477        implicit none
2478        ! si x<0, on plante.
2479        ! si très limite, on le met à 0.
2480
2481        ! input:
2482        real x
2483        character*(*) err_msg ! message d''erreur à afficher
2484
2485        ! output
2486        integer iso_verif_positif_nostop
2487
2488        ! locals
2489        !integer iso_verif_positif_choix_nostop
2490
2491        iso_verif_positif_nostop=iso_verif_positif_choix_nostop &
2492     &          (x,ridicule,err_msg)
2493
2494#ifdef ISOVERIF
2495#else
2496        write(*,*) 'iso_verif 837: err_msg=',err_msg, &
2497     &           ': pourquoi verif?'
2498        stop
2499
2500#endif         
2501       
2502        return
2503        end function iso_verif_positif_nostop
2504
2505        !******************
2506        function iso_verif_positif_choix_nostop(x,ridic,err_msg)
2507        implicit none
2508        ! si x<0, on plante.
2509        ! si très limite, on le met à 0.
2510
2511        ! input:
2512        real x
2513        real ridic
2514        character*(*) err_msg ! message d''erreur à afficher
2515
2516        ! output
2517        integer iso_verif_positif_choix_nostop
2518
2519        if (x.lt.-ridic) then
2520            write(*,*) 'erreur detectee par iso_verif_positif_nostop:'
2521            write(*,*) err_msg
2522            write(*,*) 'x=',x
2523            iso_verif_positif_choix_nostop=1
2524        else
2525!            x=max(x,0.0)
2526            iso_verif_positif_choix_nostop=0
2527        endif
2528
2529#ifdef ISOVERIF
2530#else
2531        write(*,*) 'iso_verif 877: err_msg=',err_msg, &
2532     &           ': pourquoi verif?'
2533        stop
2534#endif
2535       
2536        return
2537        end function iso_verif_positif_choix_nostop
2538
2539
2540        !**************
2541
2542       
2543        subroutine iso_verif_O18_aberrant(Rd,Ro,err_msg)
2544        implicit none
2545
2546        ! vérifie que:
2547        ! 1) deltaD et deltaO18 dans bonne gamme
2548        ! 2) dexcess dans bonne gamme
2549        ! input:
2550        real Rd,Ro
2551        character*(*) err_msg ! message d''erreur à afficher
2552
2553        ! local
2554        !integer iso_verif_O18_aberrant_nostop
2555
2556        if (iso_verif_O18_aberrant_nostop(Rd,Ro,err_msg).eq.1) then
2557            stop
2558        endif
2559
2560        end subroutine iso_verif_O18_aberrant
2561       
2562        function iso_verif_O18_aberrant_nostop(Rd,Ro,err_msg)
2563        USE isotopes_mod, ONLY: tnat, iso_HDO, iso_O18
2564        implicit none
2565
2566        ! vérifie que:
2567        ! 1) deltaD et deltaO18 dans bonne gamme
2568        ! 2) dexcess dans bonne gamme
2569
2570        ! input:
2571        real Rd,Ro
2572        character*(*) err_msg ! message d''erreur à afficher
2573
2574        ! outputs
2575        integer iso_verif_O18_aberrant_nostop
2576
2577        !locals
2578        real deltaD,deltao,dexcess
2579
2580        deltaD=(Rd/tnat(iso_hdo)-1)*1000
2581        deltao=(Ro/tnat(iso_O18)-1)*1000
2582        dexcess=deltaD-8*deltao
2583
2584        iso_verif_O18_aberrant_nostop=0
2585        if ((deltaD.lt.deltaDmin).or.(deltao.lt.deltaDmin/2.0).or. &
2586     &        (deltaD.gt.deltalim).or.(deltao.gt.deltalim/8.0).or. &
2587     &        ((deltaD.gt.-500.0).and.((dexcess.lt.dexcess_min) &
2588     &        .or.(dexcess.gt.dexcess_max)))) then
2589            write(*,*) 'erreur detectee par iso_verif_O18_aberrant:'
2590            write(*,*) err_msg
2591            write(*,*) 'delta180=',deltao
2592            write(*,*) 'deltaD=',deltaD
2593            write(*,*) 'Dexcess=',dexcess
2594            write(*,*) 'tnat=',tnat
2595!            stop
2596            iso_verif_O18_aberrant_nostop=1
2597          endif
2598
2599#ifdef ISOVERIF
2600#else
2601        write(*,*) 'err_msg=',err_msg,': pourquoi verif?'
2602        stop
2603#endif                   
2604
2605        return
2606        end function iso_verif_O18_aberrant_nostop
2607
2608
2609        ! **********
2610        function deltaD(R)
2611        USE isotopes_mod, ONLY: tnat,iso_HDO
2612        implicit none
2613        real R,deltaD
2614
2615       
2616        if (iso_HDO.gt.0) then
2617           deltaD=(R/tnat(iso_HDO)-1)*1000.0
2618        else
2619            write(*,*) 'iso_verif_egalite>deltaD 260: iso_HDO.gt.0=', &
2620     &           iso_HDO.gt.0
2621        endif
2622        return
2623        end function deltaD
2624
2625        ! **********
2626        function deltaO(R)
2627        USE isotopes_mod, ONLY: tnat,iso_O18
2628        implicit none
2629        real R,deltaO
2630       
2631        if (iso_O18.gt.0) then
2632           deltaO=(R/tnat(iso_O18)-1)*1000.0
2633        else
2634            write(*,*) 'iso_verif_egalite>deltaO18 260: iso_O18.gt.0=', &
2635     &           iso_O18.gt.0
2636        endif
2637        return
2638        end function deltaO
2639
2640        ! **********
2641        function dexcess(RD,RO)
2642        USE isotopes_mod, ONLY: tnat,iso_O18,iso_HDO
2643        implicit none
2644        real RD,RO,deltaD,deltaO,dexcess
2645       
2646        if ((iso_O18.gt.0).and.(iso_HDO.gt.0)) then
2647           deltaD=(RD/tnat(iso_HDO)-1)*1000.0
2648           deltaO=(RO/tnat(iso_O18)-1)*1000.0
2649           dexcess=deltaD-8*deltaO
2650        else
2651            write(*,*) 'iso_verif_egalite 1109: iso_O18,iso_HDO=',iso_O18,iso_HDO
2652        endif
2653        return
2654        end function dexcess
2655
2656
2657        ! **********
2658        function delta_all(R,ixt)
2659        USE isotopes_mod, ONLY: tnat
2660        implicit none
2661        real R,delta_all
2662        integer ixt
2663       
2664        delta_all=(R/tnat(ixt)-1)*1000.0
2665        return
2666        end function delta_all
2667
2668        ! **********
2669        function delta_to_R(delta,ixt)
2670        USE isotopes_mod, ONLY: tnat
2671        implicit none
2672        real delta,delta_to_R
2673        integer ixt
2674       
2675        delta_to_R=(delta/1000.0+1.0)*tnat(ixt)
2676        return
2677        end function delta_to_R
2678
2679         ! **********
2680        function o17excess(R17,R18)
2681        USE isotopes_mod, ONLY: tnat,iso_O18,iso_O17
2682        implicit none
2683        real R17,R18,o17excess
2684       
2685        if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then
2686           
2687           o17excess=1e6*(log(R17/tnat(iso_o17)) &
2688     &           -0.528*log(R18/tnat(iso_O18)))
2689!           write(*,*) 'o17excess=',o17excess
2690        else
2691            write(*,*) 'iso_verif_egalite>deltaD 260: iso_O17.gt.0,18=', &
2692     &           iso_O17.gt.0,iso_O18.gt.0
2693        endif
2694        return
2695        end function o17excess
2696
2697        !       ****************
2698
2699          subroutine iso_verif_egalite_vect2D( &
2700     &           xt,q,err_msg,ni,n,m)
2701       
2702        USE isotopes_mod, ONLY: iso_eau
2703          implicit none
2704
2705          ! inputs
2706          integer n,m,ni
2707          real q(n,m)
2708          real xt(ni,n,m)
2709          character*(*) err_msg
2710
2711        ! locals
2712        integer iso_verif_egalite_nostop_loc
2713        integer i,j,ixt
2714        integer ifaux,jfaux
2715
2716        !write(*,*) 'iso_verif_egalite_vect2D 1099 tmp: q(2,1),xt(iso_eau,2,1)=',q(2,1),xt(iso_eau,2,1)
2717        !write(*,*) 'ni,n,m=',ni,n,m,errmax,errmaxrel
2718        if (iso_eau.gt.0) then
2719        iso_verif_egalite_nostop_loc=0
2720        do i=1,n
2721         do j=1,m
2722          if (abs(q(i,j)-xt(iso_eau,i,j)).gt.errmax) then
2723           if (abs((q(i,j)-xt(iso_eau,i,j))/ &
2724     &           max(max(abs(q(i,j)),abs(xt(iso_eau,i,j))),1e-18)) &
2725     &           .gt.errmaxrel) then
2726              iso_verif_egalite_nostop_loc=1
2727              ifaux=i
2728              jfaux=j
2729           endif
2730          endif
2731         enddo !do j=1,m
2732        enddo !do i=1,n
2733
2734        if (iso_verif_egalite_nostop_loc.eq.1) then
2735          write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:'
2736          write(*,*) err_msg
2737          write(*,*) 'i,j=',ifaux,jfaux
2738          write(*,*) 'xt,q=',xt(iso_eau,ifaux,jfaux),q(ifaux,jfaux)
2739          stop
2740        endif
2741        endif
2742       
2743#ifdef ISOVERIF
2744        call iso_verif_noNaN_vect2D(xt,err_msg,ni,n,m)
2745#endif         
2746
2747        return
2748        end subroutine iso_verif_egalite_vect2D
2749
2750        subroutine iso_verif_egalite_vect1D( &
2751     &           xt,q,err_msg,ni,n)
2752
2753        USE isotopes_mod, ONLY: iso_eau
2754        implicit none
2755
2756        ! inputs
2757        integer n,ni
2758        real q(n)
2759        real xt(ni,n)
2760        character*(*) err_msg
2761
2762        ! locals
2763        integer iso_verif_egalite_nostop_loc
2764        integer i
2765        integer ifaux
2766
2767        if (iso_eau.gt.0) then
2768        iso_verif_egalite_nostop_loc=0
2769        do i=1,n
2770          if (abs(q(i)-xt(iso_eau,i)).gt.errmax) then
2771           if (abs((q(i)-xt(iso_eau,i))/ &
2772     &           max(max(abs(q(i)),abs(xt(iso_eau,i))),1e-18)) &
2773     &           .gt.errmaxrel) then
2774              iso_verif_egalite_nostop_loc=1
2775              ifaux=i
2776           endif !if (abs((q(i)-xt(iso_eau,i))/
2777          endif !if (abs(q(i)-xt(iso_eau,i)).gt.errmax) then
2778        enddo !do i=1,n
2779
2780        if (iso_verif_egalite_nostop_loc.eq.1) then
2781          write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:'
2782          write(*,*) err_msg
2783          write(*,*) 'i=',ifaux
2784          write(*,*) 'xt,q=',xt(iso_eau,ifaux),q(ifaux)
2785          stop
2786        endif  !if (iso_verif_egalite_nostop.eq.1) then
2787        endif !if (iso_eau.gt.0) then
2788
2789        end subroutine iso_verif_egalite_vect1D       
2790
2791        subroutine iso_verif_egalite_std_vect( &
2792     &           a,b,err_msg,n,m,errmax,errmaxrel)
2793
2794          implicit none
2795
2796          ! inputs
2797          integer n,m,ni
2798          real a(n,m)
2799          real b(n,m)
2800          character*(*) err_msg
2801          real errmax,errmaxrel
2802
2803        ! locals
2804        integer iso_verif_egalite_nostop_loc
2805        integer i,j
2806        integer ifaux,jfaux
2807
2808        iso_verif_egalite_nostop_loc=0
2809        do i=1,n
2810         do j=1,m
2811          if (abs(a(i,j)-b(i,j)).gt.errmax) then
2812           if (abs((a(i,j)-b(i,j))/ &
2813     &           max(max(abs(a(i,j)),abs(b(i,j))),1e-18)) &
2814     &           .gt.errmaxrel) then
2815              iso_verif_egalite_nostop_loc=1
2816              ifaux=i
2817              jfaux=j
2818           endif
2819          endif
2820         enddo !do j=1,m
2821        enddo !do i=1,n
2822
2823        if (iso_verif_egalite_nostop_loc.eq.1) then
2824          write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:'
2825          write(*,*) err_msg
2826          write(*,*) 'i,j=',ifaux,jfaux
2827          write(*,*) 'a,b=',a(ifaux,jfaux),b(ifaux,jfaux)
2828          stop
2829        endif
2830
2831        return
2832        end subroutine iso_verif_egalite_std_vect
2833
2834        subroutine iso_verif_aberrant_vect2D( &
2835     &           xt,q,err_msg,ni,n,m)
2836        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
2837          implicit none
2838
2839          ! inputs   
2840          integer n,m,ni
2841          real q(n,m)
2842          real xt(ni,n,m)
2843          character*(*) err_msg
2844
2845        ! locals
2846        integer iso_verif_aberrant_nostop_loc
2847        integer i,j
2848        integer ifaux,jfaux
2849        !real deltaD
2850
2851        if (iso_HDO.gt.0) then
2852        iso_verif_aberrant_nostop_loc=0
2853        do i=1,n
2854         do j=1,m
2855          if (q(i,j).gt.ridicule) then
2856            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
2857     &                   .gt.deltalim).or. &
2858     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
2859     &                   .lt.-borne)) then   
2860              iso_verif_aberrant_nostop_loc=1
2861              ifaux=i
2862              jfaux=j
2863           endif
2864          endif
2865         enddo !do j=1,m
2866        enddo !do i=1,n
2867
2868        if (iso_verif_aberrant_nostop_loc.eq.1) then
2869          write(*,*) 'erreur detectee par iso_verif_aberrant_vect2D:'
2870          write(*,*) err_msg
2871          write(*,*) 'i,j=',ifaux,jfaux
2872          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
2873     &           /q(ifaux,jfaux))
2874          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
2875          stop
2876        endif 
2877        endif !if (iso_HDO.gt.0) then
2878
2879        end subroutine iso_verif_aberrant_vect2D       
2880
2881        subroutine iso_verif_aberrant_enc_vect2D( &
2882     &           xt,q,err_msg,ni,n,m)
2883
2884        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
2885          implicit none
2886
2887          ! inputs   
2888          integer n,m,ni
2889          real q(n,m)
2890          real xt(ni,n,m)
2891          character*(*) err_msg
2892
2893        ! locals
2894        integer iso_verif_aberrant_nostop_loc
2895        integer i,j
2896        integer ifaux,jfaux
2897        !real deltaD
2898
2899        if (iso_HDO.gt.0) then
2900        iso_verif_aberrant_nostop_loc=0
2901        do i=1,n
2902         do j=1,m
2903          if (q(i,j).gt.ridicule) then
2904            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
2905     &                   .gt.deltalim).or. &
2906     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
2907     &                   .lt.deltaDmin).or. &
2908     &           (xt(iso_HDO,i,j).lt.-borne).or. &
2909     &           (xt(iso_HDO,i,j).gt.borne)) then     
2910              iso_verif_aberrant_nostop_loc=1
2911              ifaux=i
2912              jfaux=j
2913           endif
2914          endif
2915         enddo !do j=1,m
2916        enddo !do i=1,n
2917
2918        if (iso_verif_aberrant_nostop_loc.eq.1) then
2919          write(*,*) 'erreur detectee par ', &
2920     &           'iso_verif_aberrant_enc_vect2D:'
2921          write(*,*) err_msg
2922          write(*,*) 'i,j=',ifaux,jfaux
2923          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
2924     &           /q(ifaux,jfaux))
2925          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
2926          write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux)
2927          call abort_physic('isotopes_verif_mod','iso_verif_aberrant_enc_vect2D',1)
2928        endif 
2929        endif !if (iso_HDO.gt.0) then
2930
2931        end subroutine iso_verif_aberrant_enc_vect2D       
2932
2933       
2934        subroutine iso_verif_aberrant_enc_vect2D_ns( &
2935     &           xt,q,err_msg,ni,n,m)
2936
2937        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
2938          implicit none
2939
2940          ! inputs   
2941          integer n,m,ni
2942          real q(n,m)
2943          real xt(ni,n,m)
2944          character*(*) err_msg
2945
2946        ! locals
2947        integer iso_verif_aberrant_nostop_loc
2948        integer i,j
2949        integer ifaux,jfaux
2950        !real deltaD
2951
2952        if (iso_HDO.gt.0) then
2953        iso_verif_aberrant_nostop_loc=0
2954        do i=1,n
2955         do j=1,m
2956          if (q(i,j).gt.ridicule) then
2957            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
2958     &                   .gt.deltalim).or. &
2959     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
2960     &                   .lt.deltaDmin)) then   
2961              iso_verif_aberrant_nostop_loc=1
2962              ifaux=i
2963              jfaux=j
2964           endif
2965          endif
2966         enddo !do j=1,m
2967        enddo !do i=1,n
2968
2969        if (iso_verif_aberrant_nostop_loc.eq.1) then
2970          write(*,*) 'erreur detectee par ', &
2971     &           'iso_verif_aberrant_vect2D_ns:'
2972          write(*,*) err_msg
2973          write(*,*) 'i,j=',ifaux,jfaux
2974          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
2975     &           /q(ifaux,jfaux))
2976          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
2977!          stop
2978        endif 
2979        endif !if (iso_HDO.gt.0) then
2980
2981        end subroutine iso_verif_aberrant_enc_vect2D_ns       
2982
2983
2984         subroutine iso_verif_aberrant_vect2Dch( &
2985     &           xt,q,err_msg,ni,n,m,deltaDmax)
2986
2987        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
2988          implicit none
2989
2990
2991          ! inputs   
2992          integer n,m,ni
2993          real q(n,m)
2994          real xt(ni,n,m)
2995          character*(*) err_msg
2996          real deltaDmax
2997
2998        ! locals
2999        integer iso_verif_aberrant_nostop_loc
3000        integer i,j
3001        integer ifaux,jfaux
3002        !real deltaD
3003
3004        if (iso_HDO.gt.0) then
3005        iso_verif_aberrant_nostop_loc=0
3006        do i=1,n
3007         do j=1,m
3008          if (q(i,j).gt.ridicule) then
3009            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
3010     &                   .gt.deltaDmax).or. &
3011     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
3012     &                   .lt.-borne)) then   
3013              iso_verif_aberrant_nostop_loc=1
3014              ifaux=i
3015              jfaux=j
3016           endif
3017          endif
3018         enddo !do j=1,m
3019        enddo !do i=1,n
3020
3021        if (iso_verif_aberrant_nostop_loc.eq.1) then
3022          write(*,*) 'erreur detectee par iso_verif_aberrant_vect2D:'
3023          write(*,*) err_msg
3024          write(*,*) 'i,j=',ifaux,jfaux
3025          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
3026     &           /q(ifaux,jfaux))
3027          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
3028          stop
3029        endif 
3030        endif !if (iso_HDO.gt.0) then
3031
3032        end subroutine iso_verif_aberrant_vect2Dch     
3033
3034        subroutine iso_verif_O18_aberrant_enc_vect2D( &
3035     &           xt,q,err_msg,ni,n,m)
3036
3037        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO,iso_O18
3038          implicit none
3039
3040          ! inputs   
3041          integer n,m,ni
3042          real q(n,m)
3043          real xt(ni,n,m)
3044          character*(*) err_msg
3045
3046        ! locals
3047        integer iso_verif_aberrant_nostop_loc
3048        integer i,j
3049        integer ifaux,jfaux
3050        real deltaDloc,deltaoloc,dexcessloc
3051
3052        if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then
3053        iso_verif_aberrant_nostop_loc=0
3054        do i=1,n
3055         do j=1,m
3056          if (q(i,j).gt.ridicule) then
3057
3058            deltaDloc=(xt(iso_HDO,i,j)/q(i,j)/tnat(iso_hdo)-1)*1000
3059            deltaoloc=(xt(iso_O18,i,j)/q(i,j)/tnat(iso_O18)-1)*1000
3060            dexcessloc=deltaDloc-8*deltaoloc
3061            if ((deltaDloc.lt.deltaDmin).or.(deltaoloc.lt.deltaDmin/2.0).or. &
3062     &        (deltaDloc.gt.deltalim).or.(deltaoloc.gt.deltalim/8.0).or. &
3063     &        ((deltaDloc.gt.-500.0).and.((dexcessloc.lt.dexcess_min) &
3064     &        .or.(dexcessloc.gt.dexcess_max)))) then
3065              iso_verif_aberrant_nostop_loc=1
3066              ifaux=i
3067              jfaux=j
3068              write(*,*) 'deltaD,deltao,dexcess=',deltaDloc,deltaoloc,dexcessloc
3069           endif
3070          endif
3071         enddo !do j=1,m
3072        enddo !do i=1,n
3073
3074        if (iso_verif_aberrant_nostop_loc.eq.1) then
3075          write(*,*) 'erreur detectee par ', &
3076     &           'iso_verif_aberrant_enc_vect2D:'
3077          write(*,*) err_msg
3078          write(*,*) 'i,j=',ifaux,jfaux
3079          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
3080          write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux)
3081          call abort_physic('isotopes_verif_mod','iso_verif_aberrant_enc_vect2D',1)
3082        endif 
3083        endif !if (iso_HDO.gt.0) then
3084
3085        end subroutine iso_verif_O18_aberrant_enc_vect2D   
3086
3087
3088        subroutine select_dim23_from4D(n1,n2,n3,n4, &
3089     &          var,vec,i1,i4)
3090        implicit none
3091
3092        ! inputs
3093        integer n1,n2,n3,n4
3094        real var(n1,n2,n3,n4)
3095        integer i1,i4
3096        ! outputs
3097        real vec(n2,n3)
3098        ! locals
3099        integer i2,i3
3100
3101        do i2=1,n2
3102         do i3=1,n3
3103          vec(i2,i3)=var(i1,i2,i3,i4)
3104         enddo
3105        enddo
3106
3107        return
3108        end subroutine select_dim23_from4D
3109
3110       
3111        subroutine select_dim4_from4D(ntime,nlev,nlat,nlon, &
3112     &          var,vec,itime,ilev,ilat)
3113        implicit none
3114
3115        ! inputs
3116        integer ntime,nlev,nlat,nlon
3117        real var(ntime,nlev,nlat,nlon)
3118        integer itime,ilev,ilat
3119        ! outputs
3120        real vec(nlon)
3121        ! locals
3122        integer ilon
3123
3124        do ilon=1,nlon
3125          vec(ilon)=var(itime,ilev,ilat,ilon)
3126        enddo
3127
3128        return
3129        end subroutine select_dim4_from4D
3130
3131        subroutine select_dim5_from5D(n1,n2,n3,n4,n5, &
3132     &          var,vec,i1,i2,i3,i4)
3133        implicit none
3134
3135        ! inputs
3136        integer n1,n2,n3,n4,n5
3137        real var(n1,n2,n3,n4,n5)
3138        integer i1,i2,i3,i4
3139        ! outputs
3140        real vec(n5)
3141        ! locals
3142        integer i5
3143
3144        do i5=1,n5
3145          vec(i5)=var(i1,i2,i3,i4,i5)
3146        enddo
3147
3148        end subroutine select_dim5_from5D
3149
3150       
3151        subroutine select_dim3_from3D(ntime,nlat,nlon, &
3152     &          var,vec,itime,ilat)
3153        implicit none
3154
3155        ! inputs
3156        integer ntime,nlat,nlon
3157        real var(ntime,nlat,nlon)
3158        integer itime,ilat
3159        ! outputs
3160        real vec(nlon)
3161        ! locals
3162        integer ilon
3163
3164        do ilon=1,nlon
3165          vec(ilon)=var(itime,ilat,ilon)
3166        enddo
3167
3168        end subroutine select_dim3_from3D
3169
3170       
3171        subroutine select_dim23_from3D(n1,n2,n3, &
3172     &          var,vec,i1)
3173        implicit none
3174
3175        ! inputs
3176        integer n1,n2,n3
3177        real var(n1,n2,n3)
3178        integer i1
3179        ! outputs
3180        real vec(n2,n3)
3181        ! locals
3182        integer i2,i3
3183
3184        do i2=1,n2
3185         do i3=1,n3
3186          vec(i2,i3)=var(i1,i2,i3)
3187         enddo
3188        enddo
3189
3190        end subroutine select_dim23_from3D
3191
3192        subroutine putinto_dim23_from4D(n1,n2,n3,n4, &
3193     &          var,vec,i1,i4)
3194        implicit none
3195
3196        ! inputs
3197        integer n1,n2,n3,n4
3198        real vec(n2,n3)
3199        integer i1,i4
3200        ! inout
3201        real var(n1,n2,n3,n4)
3202        ! locals
3203        integer i2,i3
3204
3205       do i2=1,n2
3206        do i3=1,n3
3207          var(i1,i2,i3,i4)=vec(i2,i3)
3208         enddo
3209        enddo
3210
3211        end subroutine putinto_dim23_from4D
3212
3213       
3214        subroutine putinto_dim12_from4D(n1,n2,n3,n4, &
3215     &          var,vec,i3,i4)
3216        implicit none
3217
3218        ! inputs
3219        integer n1,n2,n3,n4
3220        real vec(n1,n2)
3221        integer i3,i4
3222        ! inout
3223        real var(n1,n2,n3,n4)
3224        ! locals
3225        integer i1,i2
3226
3227       do i1=1,n1
3228        do i2=1,n2
3229          var(i1,i2,i3,i4)=vec(i1,i2)
3230         enddo
3231        enddo
3232
3233        end subroutine putinto_dim12_from4D
3234       
3235        subroutine putinto_dim23_from3D(n1,n2,n3, &
3236     &          var,vec,i1)
3237        implicit none
3238
3239        ! inputs
3240        integer n1,n2,n3
3241        real vec(n2,n3)
3242        integer i1
3243        ! inout
3244        real var(n1,n2,n3)
3245        ! locals
3246        integer i2,i3
3247
3248       do i2=1,n2
3249        do i3=1,n3
3250          var(i1,i2,i3)=vec(i2,i3)
3251         enddo
3252        enddo
3253
3254        end subroutine putinto_dim23_from3D
3255
3256       
3257
3258        subroutine iso_verif_noNaN_par2D(x,err_msg,ni,n,m,ib,ie)
3259        implicit none
3260        ! si x est NaN, on affiche message
3261        ! d'erreur et return 1 si erreur
3262
3263        ! input:
3264          integer n,m,ni,ib,ie
3265        real x(ni,n,m)
3266        character*(*) err_msg ! message d''erreur à afficher
3267
3268        ! output
3269
3270        ! locals       
3271        integer i,j,ixt
3272
3273      do i=ib,ie
3274       do j=1,m
3275        do ixt=1,ni
3276         if ((x(ixt,i,j).gt.-borne).and. &
3277     &            (x(ixt,i,j).lt.borne)) then
3278         else !if ((x(ixt,i,j).gt.-borne).and.
3279            write(*,*) 'erreur detectee par iso_verif_nonNAN:'
3280            write(*,*) err_msg
3281            write(*,*) 'x,ixt,i,j=',x(ixt,i,j),ixt,i,j
3282            stop
3283         endif  !if ((x(ixt,i,j).gt.-borne).and.
3284        enddo !do ixt=1,ni
3285       enddo !do j=1,m
3286      enddo !do i=1,n     
3287
3288#ifdef ISOVERIF
3289#else
3290        write(*,*) 'err_msg iso1772=',err_msg,': pourquoi verif?'
3291        stop
3292#endif           
3293
3294        return
3295        end subroutine iso_verif_noNaN_par2D
3296
3297       
3298        subroutine iso_verif_aberrant_enc_par2D( &
3299     &           xt,q,err_msg,ni,n,m,ib,ie)
3300
3301        use isotopes_mod, ONLY: ridicule,tnat,iso_HDO
3302          implicit none
3303
3304          ! inputs   
3305          integer n,m,ni,ib,ie
3306          real q(n,m)
3307          real xt(ni,n,m)
3308          character*(*) err_msg
3309
3310        ! locals
3311        integer iso_verif_aberrant_nostop_loc
3312        integer i,j
3313        integer ifaux,jfaux
3314        !real deltaD
3315
3316        if (iso_HDO.gt.0) then
3317        iso_verif_aberrant_nostop_loc=0
3318        do i=ib,ie
3319         do j=1,m
3320          if (q(i,j).gt.ridicule) then
3321            if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
3322     &                   .gt.deltalim).or. &
3323     &          ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 &
3324     &                   .lt.deltaDmin)) then   
3325              iso_verif_aberrant_nostop_loc=1
3326              ifaux=i
3327              jfaux=j
3328           endif
3329          endif
3330         enddo !do j=1,m
3331        enddo !do i=1,n
3332
3333        if (iso_verif_aberrant_nostop_loc.eq.1) then
3334          write(*,*) 'erreur detectee par iso_verif_aberrant_par2D:'
3335          write(*,*) err_msg
3336          write(*,*) 'i,j=',ifaux,jfaux
3337          write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) &
3338     &           /q(ifaux,jfaux))
3339          write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux)
3340          write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux)
3341          stop
3342        endif 
3343        endif !if (iso_HDO.gt.0) then
3344
3345        end subroutine iso_verif_aberrant_enc_par2D       
3346
3347       
3348          subroutine iso_verif_egalite_par2D( &
3349     &           xt,q,err_msg,ni,n,m,ib,ie)
3350       
3351        USE isotopes_mod, ONLY: iso_eau
3352          implicit none
3353
3354          ! inputs
3355          integer n,m,ni,ib,ie
3356          real q(n,m)
3357          real xt(ni,n,m)
3358          character*(*) err_msg
3359
3360        ! locals
3361        integer iso_verif_egalite_nostop_loc
3362        integer i,j
3363        integer ifaux,jfaux
3364
3365        if (iso_eau.gt.0) then
3366        iso_verif_egalite_nostop_loc=0
3367        do i=ib,ie
3368         do j=1,m
3369          if (abs(q(i,j)-xt(iso_eau,i,j)).gt.errmax) then
3370           if (abs((q(i,j)-xt(iso_eau,i,j))/ &
3371     &           max(max(abs(q(i,j)),abs(xt(iso_eau,i,j))),1e-18)) &
3372     &           .gt.errmaxrel) then
3373              iso_verif_egalite_nostop_loc=1
3374              ifaux=i
3375              jfaux=j
3376           endif
3377          endif
3378         enddo !do j=1,m
3379        enddo !do i=1,n
3380
3381        if (iso_verif_egalite_nostop_loc.eq.1) then
3382          write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:'
3383          write(*,*) err_msg
3384          write(*,*) 'i,j=',ifaux,jfaux
3385          write(*,*) 'xt,q=',xt(iso_eau,ifaux,jfaux),q(ifaux,jfaux)
3386          stop
3387        endif
3388        endif
3389
3390        end subroutine iso_verif_egalite_par2D
3391
3392#ifdef ISOTRAC
3393
3394      function iso_verif_traceur_choix_nostop(x,err_msg, &
3395     &       errmax,errmaxrel,ridicule_trac,deltalimtrac)
3396        use isotopes_mod, ONLY: iso_HDO
3397        implicit none
3398        ! vérifier des choses sur les traceurs
3399        ! * toutes les zones donne t l'istope total
3400        ! * pas de deltaD aberrant
3401       
3402       ! inputs
3403       real x(ntraciso)
3404       character*(*) err_msg ! message d''erreur à afficher
3405       real errmax,errmaxrel,ridicule_trac,deltalimtrac
3406
3407       ! output
3408       integer iso_verif_traceur_choix_nostop
3409
3410       ! locals
3411       !integer iso_verif_traceur_noNaN_nostop
3412       !integer iso_verif_tracm_choix_nostop
3413       !integer iso_verif_tracdD_choix_nostop
3414       !integer iso_verif_tracpos_choix_nostop
3415
3416        iso_verif_traceur_choix_nostop=0 
3417
3418        ! verif noNaN
3419        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
3420             iso_verif_traceur_choix_nostop=1
3421        endif
3422       
3423        ! verif masse
3424        if (iso_verif_tracm_choix_nostop(x,err_msg, &
3425     &           errmax,errmaxrel).eq.1) then
3426             iso_verif_traceur_choix_nostop=1
3427        endif             
3428
3429        ! verif deltaD
3430        if (iso_HDO.gt.0) then
3431        if (iso_verif_tracdD_choix_nostop(x,err_msg, &
3432     &           ridicule_trac,deltalimtrac).eq.1) then
3433             iso_verif_traceur_choix_nostop=1
3434        endif 
3435        endif !if (iso_HDO.gt.0) then 
3436
3437        ! verif pas aberramment negatif
3438        if (iso_verif_tracpos_choix_nostop(x,err_msg, &
3439     &           1e-5).eq.1) then
3440             iso_verif_traceur_choix_nostop=1
3441        endif
3442
3443        end function iso_verif_traceur_choix_nostop
3444
3445        function iso_verif_tracnps_choix_nostop(x,err_msg, &
3446     &       errmax,errmaxrel,ridicule_trac,deltalimtrac)
3447        USE isotopes_mod, ONLY: iso_HDO
3448        implicit none
3449        ! vérifier des choses sur les traceurs
3450        ! * toutes les zones donne t l'istope total
3451        ! * pas de deltaD aberrant
3452        ! on ne vérfie pas la positivité
3453       
3454       ! inputs
3455       real x(ntraciso)
3456       character*(*) err_msg ! message d''erreur à afficher
3457       real errmax,errmaxrel,ridicule_trac,deltalimtrac
3458
3459       ! output
3460       integer iso_verif_tracnps_choix_nostop
3461
3462       ! locals
3463       !integer iso_verif_traceur_noNaN_nostop
3464       !integer iso_verif_tracm_choix_nostop
3465       !integer iso_verif_tracdD_choix_nostop
3466
3467        iso_verif_tracnps_choix_nostop=0 
3468
3469        ! verif noNaN
3470        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
3471             iso_verif_tracnps_choix_nostop=1
3472        endif
3473       
3474        ! verif masse
3475        if (iso_verif_tracm_choix_nostop(x,err_msg, &
3476     &           errmax,errmaxrel).eq.1) then
3477             iso_verif_tracnps_choix_nostop=1
3478        endif             
3479
3480        ! verif deltaD
3481        if (iso_HDO.gt.0) then
3482        if (iso_verif_tracdD_choix_nostop(x,err_msg, &
3483     &           ridicule_trac,deltalimtrac).eq.1) then
3484             iso_verif_tracnps_choix_nostop=1
3485        endif   
3486        endif ! if (iso_HDO.gt.0) then 
3487
3488        return
3489        end function iso_verif_tracnps_choix_nostop
3490
3491        function iso_verif_tracpos_choix_nostop(x,err_msg,seuil)
3492        use isotopes_mod, only: isoName
3493        implicit none
3494
3495        ! inputs
3496       real x(ntraciso)
3497       character*(*) err_msg ! message d''erreur à afficher
3498       real seuil
3499
3500       ! output
3501       integer iso_verif_tracpos_choix_nostop
3502
3503       ! locals
3504       integer lnblnk
3505       integer iiso,ixt
3506       !integer iso_verif_positif_choix_nostop
3507
3508       iso_verif_tracpos_choix_nostop=0
3509
3510       do ixt=niso+1,ntraciso
3511          if (iso_verif_positif_choix_nostop(x(ixt),seuil,err_msg// &
3512     &           ', verif positif, iso'//TRIM(isoName(ixt))).eq.1) then
3513            iso_verif_tracpos_choix_nostop=1
3514          endif
3515        enddo
3516
3517        end function iso_verif_tracpos_choix_nostop
3518
3519
3520        function iso_verif_traceur_noNaN_nostop(x,err_msg)
3521        use isotopes_mod, only: isoName
3522        implicit none
3523
3524        ! on vérifie juste que pas NaN
3525        ! inputs
3526       real x(ntraciso)
3527       character*(*) err_msg ! message d''erreur à afficher
3528
3529       ! output
3530       integer iso_verif_traceur_noNaN_nostop
3531
3532       ! locals
3533       integer lnblnk
3534       integer iiso,ixt
3535       !integer iso_verif_nonaN_nostop
3536
3537       iso_verif_traceur_noNaN_nostop=0
3538
3539        do ixt=niso+1,ntraciso
3540!          write(*,*) 'iso_verif_traceurs 154: iiso,ixt=',iiso,ixt
3541          if (iso_verif_noNaN_nostop(x(ixt),err_msg// &
3542     &           ', verif trac no NaN, iso'//TRIM(isoName(ixt))) &
3543     &           .eq.1) then
3544            iso_verif_traceur_noNaN_nostop=1
3545          endif
3546        enddo
3547
3548        end function iso_verif_traceur_noNaN_nostop
3549
3550        function iso_verif_tracm_choix_nostop(x,err_msg, &
3551     &           errmaxin,errmaxrelin)
3552
3553        use isotopes_mod, ONLY: ridicule,isoName
3554        ! on vérifie juste bilan de masse
3555        implicit none
3556       
3557        ! inputs
3558        real x(ntraciso)
3559        character*(*) err_msg ! message d''erreur à afficher
3560        real errmaxin,errmaxrelin
3561
3562        ! output
3563        integer iso_verif_tracm_choix_nostop
3564
3565       ! locals
3566       !integer iso_verif_egalite_choix_nostop
3567       integer iiso,izone,ixt
3568       real xtractot
3569
3570       iso_verif_tracm_choix_nostop=0
3571
3572        do iiso=1,niso
3573
3574          xtractot=0.0
3575          do izone=1,nzone 
3576            ixt=itZonIso(izone,iiso)
3577            xtractot=xtractot+x(ixt)
3578          enddo
3579
3580          if (iso_verif_egalite_choix_nostop(xtractot,x(iiso), &
3581     &        err_msg//', verif trac egalite1, iso '// &
3582     &        TRIM(isoName(iiso)), &
3583     &        errmaxin,errmaxrelin).eq.1) then
3584            write(*,*) 'iso_verif_traceur 202: x=',x
3585!            write(*,*) 'xtractot=',xtractot
3586            do izone=1,nzone 
3587              ixt=itZonIso(izone,iiso)
3588              write(*,*) 'izone,iiso,ixt=',izone,iiso,ixt
3589            enddo
3590            iso_verif_tracm_choix_nostop=1
3591          endif
3592
3593          ! verif ajoutée le 19 fev 2011
3594          if ((abs(xtractot).lt.ridicule**2).and. &
3595     &           (abs(x(iiso)).gt.ridicule)) then
3596            write(*,*) err_msg,', verif masse traceurs, iso ', &
3597     &          TRIM(isoName(iiso))
3598            write(*,*) 'iso_verif_traceur 209: x=',x
3599!            iso_verif_tracm_choix_nostop=1
3600          endif
3601
3602        enddo !do iiso=1,ntraceurs_iso 
3603
3604        end function iso_verif_tracm_choix_nostop
3605
3606        function iso_verif_tracdD_choix_nostop(x,err_msg, &
3607     &           ridicule_trac,deltalimtrac)
3608
3609        USE isotopes_mod, ONLY: iso_eau, iso_HDO
3610        use isotrac_mod, only: strtrac
3611        ! on vérifie juste deltaD
3612        implicit none
3613               
3614        ! inputs
3615        real x(ntraciso)
3616        character*(*) err_msg ! message d''erreur à afficher
3617        real ridicule_trac,deltalimtrac
3618
3619        ! output
3620        integer iso_verif_tracdD_choix_nostop       
3621
3622       ! locals
3623       integer izone,ieau,ixt
3624       !integer iso_verif_aberrant_choix_nostop
3625
3626        iso_verif_tracdD_choix_nostop=0
3627
3628        if ((iso_eau.gt.0).and.(iso_HDO.gt.0)) then
3629        do izone=1,nzone
3630             ieau=itZonIso(izone,iso_eau)
3631             ixt=itZonIso(izone,iso_HDO)
3632
3633             if (iso_verif_aberrant_choix_nostop(x(ixt),x(ieau), &
3634     &           ridicule_trac,deltalimtrac,err_msg// &
3635     &           ', verif trac no aberrant zone '//strtrac(izone)) &
3636     &           .eq.1) then
3637               iso_verif_tracdD_choix_nostop=1
3638             endif
3639!             if (x(ieau).gt.ridicule) then
3640!              call iso_verif_aberrant(x(ixt)/x(ieau),
3641!     :           err_msg//', verif trac no aberrant zone '
3642!     :           //strtrac(izone))
3643!             endif
3644        enddo !do izone=1,nzone
3645       endif ! if ((iso_eau.gt.0).and.(iso_HDO.gt.0)) then
3646
3647       end function iso_verif_tracdD_choix_nostop
3648
3649INTEGER FUNCTION iso_verif_tag17_q_deltaD_chns(x,err_msg) RESULT(res)
3650  USE isotopes_mod, ONLY: iso_HDO, iso_eau, ridicule
3651  USE isotrac_mod,  ONLY: nzone_temp, option_traceurs
3652  IMPLICIT NONE
3653  REAL,             INTENT(IN) :: x(ntraciso)
3654  CHARACTER(LEN=*), INTENT(IN) :: err_msg
3655  INTEGER :: ieau, ixt, ieau1
3656  res = 0
3657  IF(ALL([17,18]/=option_traceurs)) RETURN
3658  !--- Check whether * deltaD(highest tagging layer) < 200 permil
3659  !                  * q <
3660  ieau=itZonIso(nzone_temp,iso_eau)
3661  ixt=itZonIso(nzone_temp,iso_HDO)
3662  IF(x(ieau)>ridicule) THEN
3663    IF(iso_verif_positif_nostop(-200.0-deltaD(x(ixt)/x(ieau)), err_msg//': deltaDt05 trop fort')==1) THEN
3664      res=1; write(*,*) 'x=',x
3665    END IF
3666  END IF
3667  IF(iso_verif_positif_nostop(2.0e-3-x(ieau),err_msg//': qt05 trop fort')==1) THEN
3668    res=1; write(*,*) 'x=',x
3669  END IF
3670  !--- Check whether q is small ; then, qt01 < 10%
3671  IF(x(iso_eau)<2.0e-3) THEN
3672    ieau1= itZonIso(1,iso_eau)
3673    IF(iso_verif_positif_nostop(0.1-(x(ieau1)/x(iso_eau)),err_msg//': qt01 trop abondant')==1) THEN
3674      res=1; write(*,*) 'x=',x
3675    END IF
3676  END IF
3677END FUNCTION iso_verif_tag17_q_deltaD_chns
3678
3679SUBROUTINE iso_verif_trac17_q_deltaD(x,err_msg)
3680  USE isotrac_mod,  ONLY: nzone_temp, option_traceurs
3681  IMPLICIT NONE
3682  REAL,             INTENT(IN) :: x(ntraciso)
3683  CHARACTER(LEN=*), INTENT(IN) :: err_msg
3684  IF(ALL([17,18]/=option_traceurs)) RETURN
3685  IF(nzone_temp>=5) THEN
3686    IF(iso_verif_tag17_q_deltaD_chns(x,err_msg)==1) STOP
3687  END IF
3688END SUBROUTINE iso_verif_trac17_q_deltaD
3689
3690      subroutine iso_verif_traceur(x,err_msg)
3691        use isotrac_mod, only: ridicule_trac
3692        implicit none
3693        ! vérifier des choses sur les traceurs
3694        ! * toutes les zones donne t l'istope total
3695        ! * pas de deltaD aberrant
3696
3697        ! on prend les valeurs pas défaut pour
3698        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
3699       
3700       ! inputs
3701       real x(ntraciso)
3702       character*(*) err_msg ! message d''erreur à afficher
3703
3704       ! locals
3705       !integer iso_verif_traceur_choix_nostop 
3706
3707        if (iso_verif_traceur_choix_nostop(x,err_msg, &
3708     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
3709     &       .eq.1) then
3710                stop
3711        endif
3712
3713        end subroutine iso_verif_traceur
3714
3715       
3716      subroutine iso_verif_traceur_retourne3D(x,n1,n2,n3, &
3717     &           i1,i2,i3,err_msg)
3718        use isotrac_mod, only: ridicule_trac
3719
3720        implicit none
3721        ! vérifier des choses sur les traceurs
3722        ! * toutes les zones donne t l'istope total
3723        ! * pas de deltaD aberrant
3724
3725        ! on prend les valeurs pas défaut pour
3726        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
3727       
3728       ! inputs
3729       integer n1,n2,n3
3730       real x(n1,n2,n3,ntraciso)
3731       character*(*) err_msg ! message d''erreur à afficher
3732       integer i1,i2,i3
3733
3734       ! locals
3735       !integer iso_verif_traceur_choix_nostop 
3736       real xiso(ntraciso)
3737
3738        call select_dim4_from4D(n1,n2,n3,ntraciso, &
3739     &      x,xiso,i1,i2,i3)
3740        if (iso_verif_traceur_choix_nostop(xiso,err_msg, &
3741     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
3742     &       .eq.1) then
3743                stop
3744        endif
3745
3746        end subroutine iso_verif_traceur_retourne3D
3747
3748        subroutine iso_verif_traceur_retourne4D(x,n1,n2,n3,n4, &
3749     &           i1,i2,i3,i4,err_msg)
3750        use isotrac_mod, only: ridicule_trac
3751
3752        implicit none
3753        ! vérifier des choses sur les traceurs
3754        ! * toutes les zones donne t l'istope total
3755        ! * pas de deltaD aberrant
3756
3757        ! on prend les valeurs pas défaut pour
3758        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
3759       
3760       ! inputs
3761       integer n1,n2,n3,n4
3762       real x(n1,n2,n3,n4,ntraciso)
3763       character*(*) err_msg ! message d''erreur à afficher
3764       integer i1,i2,i3,i4
3765
3766       ! locals
3767       !integer iso_verif_traceur_choix_nostop 
3768       real xiso(ntraciso)
3769
3770        call select_dim5_from5D(n1,n2,n3,n4,ntraciso, &
3771     &      x,xiso,i1,i2,i3,i4)
3772        if (iso_verif_traceur_choix_nostop(xiso,err_msg, &
3773     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
3774     &       .eq.1) then
3775                stop
3776        endif
3777
3778        end subroutine iso_verif_traceur_retourne4D
3779
3780       
3781      subroutine iso_verif_traceur_retourne2D(x,n1,n2, &
3782     &           i1,i2,err_msg)
3783        use isotrac_mod, only: ridicule_trac
3784        implicit none
3785        ! vérifier des choses sur les traceurs
3786        ! * toutes les zones donne t l'istope total
3787        ! * pas de deltaD aberrant
3788
3789        ! on prend les valeurs pas défaut pour
3790        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
3791       
3792       ! inputs
3793       integer n1,n2
3794       real x(n1,n2,ntraciso)
3795       character*(*) err_msg ! message d''erreur à afficher
3796       integer i1,i2
3797
3798       ! locals
3799       !integer iso_verif_traceur_choix_nostop 
3800       real xiso(ntraciso)
3801
3802        call select_dim3_from3D(n1,n2,ntraciso, &
3803     &      x,xiso,i1,i2)
3804        if (iso_verif_traceur_choix_nostop(xiso,err_msg, &
3805     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
3806     &       .eq.1) then
3807                stop
3808        endif
3809
3810        end subroutine iso_verif_traceur_retourne2D
3811
3812        subroutine iso_verif_traceur_vect(x,n,m,err_msg)
3813        USE isotopes_mod, ONLY: iso_HDO
3814        implicit none
3815        ! vérifier des choses sur les traceurs
3816        ! * toutes les zones donne t l'istope total
3817        ! * pas de deltaD aberrant
3818
3819        ! on prend les valeurs pas défaut pour
3820        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
3821       
3822       ! inputs
3823       integer n,m
3824       real x(ntraciso,n,m)
3825       character*(*) err_msg ! message d''erreur à afficher
3826
3827       ! locals
3828       logical iso_verif_traceur_nostop
3829       integer i,j,ixt,iiso,izone,ieau
3830       integer ifaux,jfaux,ixtfaux
3831       
3832       call iso_verif_traceur_noNaN_vect(x,n,m,err_msg)       
3833
3834        ! verif masse: iso_verif_tracm_choix_nostop
3835        call iso_verif_trac_masse_vect(x,n,m,err_msg,errmax,errmaxrel)
3836
3837        ! verif deltaD: iso_verif_tracdD_choix_nostop   
3838        if (iso_HDO.gt.0) then 
3839        call iso_verif_tracdd_vect(x,n,m,err_msg)     
3840        endif !if (iso_HDO.gt.0) then       
3841
3842        ! verif pas aberramment negatif: iso_verif_tracpos_choix_nostop
3843        call iso_verif_tracpos_vect(x,n,m,err_msg,1e-5) 
3844       
3845        end subroutine iso_verif_traceur_vect
3846
3847        subroutine iso_verif_tracnps_vect(x,n,m,err_msg)
3848        USE isotopes_mod, ONLY: iso_HDO
3849        implicit none
3850        ! vérifier des choses sur les traceurs
3851        ! * toutes les zones donne t l'istope total
3852        ! * pas de deltaD aberrant
3853
3854        ! on prend les valeurs pas défaut pour
3855        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
3856       
3857       ! inputs
3858       integer n,m
3859       real x(ntraciso,n,m)
3860       character*(*) err_msg ! message d''erreur à afficher
3861
3862       ! locals
3863       logical iso_verif_traceur_nostop
3864       integer i,j,ixt,iiso,izone,ieau
3865       integer ifaux,jfaux,ixtfaux
3866       
3867       call iso_verif_traceur_noNaN_vect(x,n,m,err_msg)       
3868
3869        ! verif masse: iso_verif_tracm_choix_nostop
3870        call iso_verif_trac_masse_vect(x,n,m,err_msg,errmax,errmaxrel)
3871
3872        ! verif deltaD: iso_verif_tracdD_choix_nostop   
3873        if (iso_HDO.gt.0) then 
3874        call iso_verif_tracdd_vect(x,n,m,err_msg)     
3875        endif !if (iso_HDO.gt.0) then       
3876       
3877        end subroutine iso_verif_tracnps_vect
3878
3879
3880        subroutine iso_verif_traceur_noNaN_vect(x,n,m,err_msg)
3881        implicit none
3882       
3883       ! inputs
3884       integer n,m
3885       real x(ntraciso,n,m)
3886       character*(*) err_msg ! message d''erreur à afficher
3887
3888       ! locals
3889       logical iso_verif_traceur_nostop
3890       integer i,j,ixt,iiso
3891       integer ifaux,jfaux,ixtfaux
3892
3893
3894       iso_verif_traceur_nostop=.false.
3895        ! verif noNaN: iso_verif_traceur_noNaN_nostop       
3896        do j=1,m
3897        do i=1,n
3898        do ixt=niso+1,ntraciso
3899          if ((x(ixt,i,j).gt.-borne).and.(x(ixt,i,j).lt.borne)) then
3900          else !if ((x.gt.-borne).and.(x.lt.borne)) then
3901              iso_verif_traceur_nostop=.true.
3902              ifaux=i
3903              jfaux=i
3904          endif !if ((x.gt.-borne).and.(x.lt.borne)) then
3905        enddo !do ixt=niso+1,ntraciso
3906        enddo ! do i=1,n
3907        enddo ! do j=1,m
3908       
3909
3910        if (iso_verif_traceur_nostop) then
3911            write(*,*) 'erreur detectée par iso_verif_nonNAN ', &
3912     &           'dans iso_verif_traceur_vect'
3913            write(*,*) ''
3914            write(*,*) err_msg
3915            write(*,*) 'x=',x(:,ifaux,jfaux)
3916            stop
3917        endif
3918
3919        end subroutine iso_verif_traceur_noNaN_vect
3920
3921       
3922        subroutine iso_verif_trac_masse_vect(x,n,m,err_msg, &
3923     &            errmax,errmaxrel)
3924        use isotopes_mod, only: isoName
3925        implicit none
3926       
3927        ! inputs
3928       integer n,m
3929       real x(ntraciso,n,m)
3930       character*(*) err_msg ! message d''erreur à afficher
3931       real errmax,errmaxrel
3932
3933       ! locals
3934       logical iso_verif_traceur_nostop
3935       integer i,j,ixt,iiso,izone
3936       integer ifaux,jfaux,ixtfaux
3937       real xtractot(n,m)
3938       real xiiso(n,m)
3939
3940        do iiso=1,niso       
3941        do j=1,m
3942         do i=1,n       
3943          xtractot(i,j)=0.0
3944          xiiso(i,j)=x(iiso,i,j)
3945          do izone=1,nzone
3946            ixt=itZonIso(izone,iiso)
3947            xtractot(i,j)=xtractot(i,j)+x(ixt,i,j)           
3948          enddo !do izone=1,nzone
3949         enddo !do i=1,n
3950        enddo !do j=1,m
3951       
3952
3953        call iso_verif_egalite_std_vect( &
3954     &           xtractot,xiiso, &
3955     &           err_msg//', verif trac egalite2, iso ' &
3956     &           //TRIM(isoName(iiso)), &
3957     &           n,m,errmax,errmaxrel)
3958        enddo !do iiso=1,niso
3959
3960        end subroutine iso_verif_trac_masse_vect
3961
3962        subroutine iso_verif_tracdd_vect(x,n,m,err_msg)
3963        use isotopes_mod, only: iso_HDO,iso_eau
3964        use isotrac_mod, only: strtrac
3965        implicit none
3966       
3967        ! inputs
3968       integer n,m
3969       real x(ntraciso,n,m)
3970       character*(*) err_msg ! message d''erreur à afficher
3971
3972       ! locals
3973       integer i,j,iiso,izone,ieau,ixt
3974       real xiiso(niso,n,m)
3975       real xeau(n,m)
3976       integer lnblnk
3977
3978       if (iso_HDO.gt.0) then
3979        do izone=1,nzone
3980          ieau=itZonIso(izone,iso_eau)
3981          do iiso=1,niso
3982           ixt=itZonIso(izone,iiso)
3983           do j=1,m
3984            do i=1,n
3985             xiiso(iiso,i,j)=x(ixt,i,j)
3986            enddo !do i=1,n
3987           enddo !do j=1,m
3988          enddo !do iiso=1,niso
3989           
3990          do j=1,m
3991           do i=1,n
3992            xeau(i,j)=x(ieau,i,j)
3993           enddo !do i=1,n
3994          enddo !do j=1,m
3995         
3996          call iso_verif_aberrant_vect2Dch( &
3997     &           xiiso,xeau,err_msg//strtrac(izone),niso,n,m, &
3998     &           deltalimtrac)
3999         enddo !do izone=1,nzone
4000        endif !if (iso_HDO.gt.0) then
4001
4002        end subroutine iso_verif_tracdd_vect
4003
4004        subroutine iso_verif_tracpos_vect(x,n,m,err_msg,seuil)
4005        implicit none
4006
4007       ! inputs
4008       integer n,m
4009       real x(ntraciso,n,m)
4010       character*(*) err_msg ! message d''erreur à afficher
4011       real seuil
4012
4013       ! locals
4014       integer i,j,ixt
4015       logical iso_verif_traceur_nostop
4016       integer ifaux,jfaux,ixtfaux
4017
4018        iso_verif_traceur_nostop=.false.       
4019        do j=1,m
4020        do i=1,n
4021        do ixt=niso+1,ntraciso
4022          if (x(ixt,i,j).lt.-seuil) then
4023              ifaux=i
4024              jfaux=j
4025              ixtfaux=ixt
4026              iso_verif_traceur_nostop=.true.
4027          endif
4028        enddo !do ixt=niso+1,ntraciso
4029        enddo !do i=1,n
4030        enddo !do j=1,m       
4031
4032        if (iso_verif_traceur_nostop) then
4033            write(*,*) 'erreur detectée par verif positif ', &
4034     &           'dans iso_verif_traceur_vect'
4035            write(*,*) ''
4036            write(*,*) err_msg
4037            write(*,*) 'x=',x(:,ifaux,jfaux)
4038            stop
4039        endif
4040
4041        end subroutine iso_verif_tracpos_vect
4042
4043
4044
4045        subroutine iso_verif_tracnps(x,err_msg)
4046        use isotrac_mod, only: ridicule_trac
4047
4048        implicit none
4049        ! vérifier des choses sur les traceurs
4050        ! * toutes les zones donne t l'istope total
4051        ! * pas de deltaD aberrant
4052
4053        ! on prend les valeurs pas défaut pour
4054        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
4055       
4056       ! inputs
4057       real x(ntraciso)
4058       character*(*) err_msg ! message d''erreur à afficher
4059
4060       ! locals
4061       !integer iso_verif_tracnps_choix_nostop 
4062
4063        if (iso_verif_tracnps_choix_nostop(x,err_msg, &
4064     &       errmax,errmaxrel,ridicule_trac,deltalimtrac) &
4065     &       .eq.1) then
4066                stop
4067        endif
4068
4069        end subroutine iso_verif_tracnps
4070
4071        subroutine iso_verif_tracpos_choix(x,err_msg,seuil)
4072        implicit none
4073        ! vérifier des choses sur les traceurs
4074        ! * toutes les zones donne t l'istope total
4075        ! * pas de deltaD aberrant
4076
4077        ! on prend les valeurs pas défaut pour
4078        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
4079       
4080       ! inputs
4081       real x(ntraciso)
4082       character*(*) err_msg ! message d''erreur à afficher
4083       real seuil
4084
4085       ! locals
4086       !integer iso_verif_tracpos_choix_nostop 
4087
4088        if (iso_verif_tracpos_choix_nostop(x,err_msg,seuil) &
4089     &       .eq.1) then
4090                stop
4091        endif
4092
4093        end subroutine iso_verif_tracpos_choix
4094
4095        subroutine iso_verif_traceur_choix(x,err_msg, &
4096     &       errmax,errmaxrel,ridicule_trac_loc,deltalimtrac)
4097        implicit none
4098        ! vérifier des choses sur les traceurs
4099        ! * toutes les zones donne t l'istope total
4100        ! * pas de deltaD aberrant
4101       
4102       ! inputs
4103       real x(ntraciso)
4104       character*(*) err_msg ! message d''erreur à afficher
4105       real errmax,errmaxrel,ridicule_trac_loc,deltalimtrac
4106
4107       ! locals
4108       !integer iso_verif_traceur_choix_nostop 
4109
4110        if (iso_verif_traceur_choix_nostop(x,err_msg, &
4111     &       errmax,errmaxrel,ridicule_trac_loc,deltalimtrac) &
4112     &       .eq.1) then
4113                stop
4114        endif
4115
4116        end subroutine iso_verif_traceur_choix
4117
4118        function iso_verif_traceur_nostop(x,err_msg)
4119        use isotrac_mod, only: ridicule_trac
4120        !use isotopes_verif, only: errmax,errmaxrel,deltalimtrac
4121        implicit none
4122        ! vérifier des choses sur les traceurs
4123        ! * toutes les zones donne t l'istope total
4124        ! * pas de deltaD aberrant
4125
4126        ! on prend les valeurs pas défaut pour
4127        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
4128       
4129       ! inputs
4130       real x(ntraciso)
4131       character*(*) err_msg ! message d''erreur à afficher
4132
4133       ! output
4134       integer iso_verif_traceur_nostop
4135
4136       ! locals
4137       !integer iso_verif_traceur_choix_nostop 
4138
4139        iso_verif_traceur_nostop= &
4140     &       iso_verif_traceur_choix_nostop(x,err_msg, &
4141     &       errmax,errmaxrel,ridicule_trac,deltalimtrac)
4142
4143        end function iso_verif_traceur_nostop
4144
4145
4146      subroutine iso_verif_traceur_justmass(x,err_msg)
4147        implicit none
4148        ! on vérifie que noNaN et masse
4149
4150        ! on prend les valeurs pas défaut pour
4151        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
4152       
4153       ! inputs
4154       real x(ntraciso)
4155       character*(*) err_msg ! message d''erreur à afficher
4156
4157       ! locals
4158       !integer iso_verif_traceur_noNaN_nostop
4159       !integer iso_verif_tracm_choix_nostop
4160
4161        ! verif noNaN
4162        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
4163             stop
4164        endif
4165       
4166        ! verif masse
4167        if (iso_verif_tracm_choix_nostop(x,err_msg, &
4168     &           errmax,errmaxrel).eq.1) then
4169             stop
4170        endif   
4171       
4172        end subroutine iso_verif_traceur_justmass
4173
4174        function iso_verif_traceur_jm_nostop(x,err_msg)
4175        implicit none
4176        ! on vérifie que noNaN et masse
4177
4178        ! on prend les valeurs pas défaut pour
4179        ! errmax,errmaxrel,ridicule_trac,deltalimtrac
4180       
4181       ! inputs
4182       real x(ntraciso)
4183       character*(*) err_msg ! message d''erreur à afficher
4184
4185       ! output
4186       integer iso_verif_traceur_jm_nostop
4187
4188       ! locals
4189!       integer iso_verif_traceur_noNaN_nostop
4190       !integer iso_verif_tracm_choix_nostop
4191
4192        iso_verif_traceur_jm_nostop=0
4193!        ! verif noNaN
4194!        if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then
4195!             iso_verif_traceur_jm_nostop=1
4196!        endif
4197       
4198        ! verif masse
4199        if (iso_verif_tracm_choix_nostop(x,err_msg, &
4200     &           errmax,errmaxrel).eq.1) then
4201             iso_verif_traceur_jm_nostop=1
4202        endif   
4203       
4204        end function iso_verif_traceur_jm_nostop
4205
4206        subroutine iso_verif_tag17_q_deltaD_vect(x,n,m,err_msg)
4207        USE isotopes_mod, ONLY: tnat,iso_eau, ridicule,iso_HDO
4208        use isotrac_mod, only: option_traceurs,nzone_temp
4209        implicit none
4210
4211        ! inputs
4212        integer n,m
4213        real x(ntraciso,n,m)
4214        character*(*) err_msg
4215
4216        ! locals
4217        !integer iso_verif_positif_nostop
4218        !real deltaD
4219        integer ieau,ixt,ieau1
4220        integer i,k
4221
4222        if ((option_traceurs.eq.17).or. &
4223     &           (option_traceurs.eq.18)) then
4224        ! verifier que deltaD du tag de la couche la plus haute <
4225        ! 200 permil, et vérifier que son q est inférieur à
4226        ieau=itZonIso(nzone_temp,iso_eau)
4227        ixt=itZonIso(nzone_temp,iso_HDO)
4228        ieau1=itZonIso(1,iso_eau)
4229        do i=1,n
4230         do k=1,m
4231           if (x(ieau,i,k).gt.ridicule) then
4232             if ((x(ixt,i,k)/x(ieau,i,k)/tnat(iso_HDO)-1)*1000 &
4233     &            .gt.-200.0) then
4234                write(*,*) err_msg,', vect:deltaDt05 trop fort'
4235                write(*,*) 'i,k=',i,k
4236                write(*,*) 'x(:,i,k)=',x(:,i,k)
4237                stop
4238             endif !if (iso_verif_positif_nostop(-200.0-deltaD(x(ixt),x(ieau)),
4239           endif !if (x(ieau).gt.ridicule) then
4240           if (x(ieau,i,k).gt.2.0e-3) then
4241                write(*,*) err_msg,', vect:qt05 trop fort'
4242                write(*,*) 'i,k=',i,k
4243                write(*,*) 'x(:,i,k)=',x(:,i,k)
4244                stop
4245           endif !if (iso_verif_positif_nostop(1.0e-3-x(ieau),
4246           if (x(iso_eau,i,k).lt.2.0e-3) then
4247                if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then
4248                   write(*,*) err_msg,', vect: qt01 trop abondant'
4249                   write(*,*) 'i,k=',i,k
4250                   write(*,*) 'ieau1,iso_eau,x(ieau1,i,k),', &
4251     &                 'x(iso_eau,i,k)=',ieau1,iso_eau, &
4252     &                  x(ieau1,i,k),x(iso_eau,i,k) 
4253                   write(*,*) 'x(:,i,k)=',x(:,i,k)
4254                   stop
4255                endif !if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then
4256            endif
4257          enddo !do k=1,m
4258        enddo !do i=1,n
4259
4260        endif !if (option_traceurs.eq.17) then
4261
4262        end subroutine iso_verif_tag17_q_deltaD_vect
4263
4264
4265        subroutine iso_verif_tag17_q_deltaD_vect_ret3D(x,n,m,nq,err_msg)
4266        USE isotopes_mod, ONLY: tnat,iso_eau,iso_HDO,ridicule
4267        use isotrac_mod, only: option_traceurs,nzone_temp
4268        implicit none
4269
4270        ! inputs
4271        integer n,m,nq
4272        real x(n,m,nq,ntraciso)
4273        character*(*) err_msg
4274
4275        ! locals
4276        !integer iso_verif_positif_nostop
4277        !real deltaD
4278        integer ieau,ixt,ieau1
4279        integer i,k,iq
4280
4281        if ((option_traceurs.eq.17).or. &
4282     &           (option_traceurs.eq.18)) then
4283        ! verifier que deltaD du tag de la couche la plus haute <
4284        ! 200 permil, et vérifier que son q est inférieur à
4285        ieau=itZonIso(nzone_temp,iso_eau)
4286        ixt=itZonIso(nzone_temp,iso_HDO)
4287        ieau1=itZonIso(1,iso_eau)
4288        do iq=1,nq
4289        do i=1,n
4290         do k=1,m
4291           if (x(i,k,iq,ieau).gt.ridicule) then
4292             if ((x(i,k,iq,ixt)/x(i,k,iq,ieau)/tnat(iso_HDO)-1)*1000 &
4293     &            .gt.-200.0) then
4294                write(*,*) err_msg,', vect:deltaDt05 trop fort'
4295                write(*,*) 'i,k=',i,k
4296                write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:)
4297                stop
4298             endif !if (iso_verif_positif_nostop(-200.0-deltaD(x(ixt),x(ieau)),
4299           endif !if (x(ieau).gt.ridicule) then
4300           if (x(i,k,iq,ieau).gt.2.0e-3) then
4301                write(*,*) err_msg,', vect:qt05 trop fort'
4302                write(*,*) 'i,k=',i,k
4303                write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:)
4304                stop
4305           endif !if (iso_verif_positif_nostop(1.0e-3-x(ieau),
4306           if (x(i,k,iq,iso_eau).lt.2.0e-3) then
4307                if (x(i,k,iq,ieau1)/x(i,k,iq,iso_eau).gt.0.1) then
4308                   write(*,*) err_msg,', vect: qt01 trop abondant'
4309                   write(*,*) 'i,k=',i,k
4310                   write(*,*) 'ieau1,iso_eau,x(i,k,iq,ieau1),', &
4311     &                 'x(i,k,iq,ieau)=',ieau1,iso_eau, &
4312     &                  x(i,k,iq,ieau1),x(i,k,iq,iso_eau) 
4313                   write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:)
4314                   stop
4315                endif !if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then
4316            endif
4317          enddo !do k=1,m
4318        enddo !do i=1,n
4319        enddo ! do iq=1,nq
4320
4321        endif !if (option_traceurs.eq.17) then
4322
4323        end subroutine iso_verif_tag17_q_deltaD_vect_ret3D
4324
4325
4326#endif
4327! endif ISOTRAC
4328
4329END MODULE isotopes_verif_mod
4330
4331#endif         
4332! endif ISOVERIF
4333
Note: See TracBrowser for help on using the repository browser.