#ifdef ISOVERIF ! $Id: $ MODULE isotopes_verif_mod USE isotopes_mod, ONLY: iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO, small => ridicule, tnat USE infotrac_phy, ONLY: isoName, isoPhas, isoZone, iqIsoPha, tracers, nqtot, isoSelect, ntraciso => ntiso, & nbIso, niso, ntiso, nphas, nzone, itZonIso, isotope, ixIso, isoFamilies, delPhase USE strings_mod, ONLY: cat, dispOutliers, lunout, num2str, msg, strStack, strIdx, maxlen USE phys_local_var_mod, ONLY: qx => qx_seri USE ioipsl_getin_p_mod, ONLY : getin_p IMPLICIT NONE PRIVATE PUBLIC :: delta2r, delta, checkNaN, checkEqu, checkMass !, deltaR PUBLIC :: r2delta, excess, checkPos, checkBnd, checkTrac, abs_err, rel_err PUBLIC :: errmax, errmaxrel, deltalim, faccond, o17_verif, tmin_verif, faible_evap, deltaDfaible, deltalim_snow, tmax_verif, nlevmaxo17, & errmax_sol PUBLIC :: iso_verif_init PUBLIC :: iso_verif_aberrant, iso_verif_aberrant_nostop PUBLIC :: iso_verif_aberrant_choix, iso_verif_aberrant_choix_nostop PUBLIC :: iso_verif_aberrant_vect2D PUBLIC :: iso_verif_aberrant_vect2Dch PUBLIC :: iso_verif_aberrant_encadre, iso_verif_aberrant_enc_nostop PUBLIC :: iso_verif_aberrant_enc_vect2D, iso_verif_aberrant_enc_vect2D_ns PUBLIC :: iso_verif_aberrant_enc_par2D PUBLIC :: iso_verif_aberrant_enc_choix_nostop PUBLIC :: iso_verif_aberrant_o17, iso_verif_aberrant_o17_nostop PUBLIC :: iso_verif_O18_aberrant, iso_verif_O18_aberrant_nostop PUBLIC :: iso_verif_O18_aberrant_enc_vect2D PUBLIC :: iso_verif_egalite, iso_verif_egalite_nostop PUBLIC :: iso_verif_egalite_choix, iso_verif_egalite_choix_nostop PUBLIC :: iso_verif_egalite_par2D PUBLIC :: iso_verif_egalite_vect1D, iso_verif_egalite_vect2D PUBLIC :: iso_verif_egalite_std_vect PUBLIC :: iso_verif_noNaN, iso_verif_noNaN_nostop PUBLIC :: iso_verif_noNaN_par2D PUBLIC :: iso_verif_noNaN_vect1D PUBLIC :: iso_verif_noNaN_vect2D PUBLIC :: iso_verif_positif, iso_verif_positif_nostop PUBLIC :: iso_verif_positif_vect PUBLIC :: iso_verif_positif_choix, iso_verif_positif_choix_nostop PUBLIC :: iso_verif_positif_choix_vect PUBLIC :: iso_verif_positif_strict, iso_verif_positif_strict_nostop #ifdef ISOTRAC PUBLIC :: iso_verif_tag17_q_deltaD_vect PUBLIC :: iso_verif_tag17_q_deltaD_vect_ret3D PUBLIC :: iso_verif_tag17_q_deltaD_chns PUBLIC :: iso_verif_trac17_q_deltaD PUBLIC :: iso_verif_tracdd_vect PUBLIC :: iso_verif_tracdD_choix_nostop PUBLIC :: iso_verif_traceur, iso_verif_traceur_nostop PUBLIC :: iso_verif_traceur_vect PUBLIC :: iso_verif_traceur_choix, iso_verif_traceur_choix_nostop PUBLIC :: iso_verif_tracnps PUBLIC :: iso_verif_tracnps_vect PUBLIC :: iso_verif_tracnps_choix_nostop PUBLIC :: iso_verif_tracpos_vect PUBLIC :: iso_verif_tracpos_choix, iso_verif_tracpos_choix_nostop PUBLIC :: iso_verif_trac_masse_vect PUBLIC :: iso_verif_traceur_justmass, iso_verif_traceur_jm_nostop PUBLIC :: iso_verif_traceur_noNaN_nostop PUBLIC :: iso_verif_traceur_noNaN_vect PUBLIC :: iso_verif_tracm_choix_nostop #endif PUBLIC :: deltaD, deltaO, dexcess, delta_all, delta_to_R, o17excess ! SUPPRIMEE: deltaDfaible_lax=-180.0 !=============================================================================================================================== INTERFACE delta2r; MODULE PROCEDURE delta2r_0D, delta2r_1D, delta2r_2D; END INTERFACE delta2r INTERFACE r2delta; MODULE PROCEDURE r2delta_0D, r2delta_1D, r2delta_2D; END INTERFACE r2delta INTERFACE checkNaN; MODULE PROCEDURE checkNaN_0D, checkNaN_1D, checkNaN_2D, checkNaN_3D; END INTERFACE checkNaN INTERFACE checkPos; MODULE PROCEDURE checkPos_0D, checkPos_1D, checkPos_2D, checkPos_3D; END INTERFACE checkPos INTERFACE checkBnd; MODULE PROCEDURE checkBnd_0D, checkBnd_1D, checkBnd_2D, checkBnd_3D; END INTERFACE checkBnd INTERFACE checkEqu; MODULE PROCEDURE checkEqu_0D, checkEqu_1D, checkEqu_2D, checkEqu_3D; END INTERFACE checkEqu INTERFACE checkMass; MODULE PROCEDURE checkMassQ_0D, checkMassQ_1D, checkMassQ_2D, checkMassQx_2D; END INTERFACE checkMass INTERFACE delta; MODULE PROCEDURE deltaQ_0D, deltaQ_1D, deltaQ_2D, deltaQx_2D; END INTERFACE delta INTERFACE deltaR; MODULE PROCEDURE deltaR_0D, deltaR_1D, deltaR_2D; END INTERFACE deltaR INTERFACE excess; MODULE PROCEDURE excessQ_0D, excessQ_1D, excessQ_2D, & excessR_0D, excessR_1D, excessR_2D, excessQx_2D; END INTERFACE excess !=============================================================================================================================== TYPE :: r; REAL, ALLOCATABLE :: r(:); END TYPE r TYPE :: i; INTEGER, ALLOCATABLE :: i(:); END TYPE i !=============================================================================================================================== ! DESCRIPTION PREVIOUS NAME REAL, PARAMETER :: abs_err = 1.e-8, & !--- Max. absolute error allowed errmax rel_err = 1.e-3, & !--- Max. relative error allowed errmax_rel max_val = 1.e19, & !--- Max. value allowed (NaN if greater) borne faccond = 1.1, & !--- In liquids, R(max_delD)*faccond is allowed min_T = 5.0, & !--- Min. realistic temperature value (in K) Tmin_verif max_T = 1000.0, & !--- Max. realistic temperature value (in K) Tmax_verif low_delD = -300.0, & !--- Max quantity of vapour making outliers from deltaDfaible ! balance with ground acceptable min_delDevap = 3.0, & !--- Less demanding with ORCHIDEE evap outliers faible_evap slope_HDO = 8.0, & !--- H[2]HO excess law: slope slope_O17 = 0.528 !--- H2[17]O excess law: slope CHARACTER(LEN=3), PARAMETER :: & progr_HDO = 'lin', & !--- H[2]HO excess law: progression type progr_O17 = 'log' !--- H2[17]O excess law: progression type REAL, SAVE :: min_delD, & !--- Min. deltaD allowed (defined in iso.def) deltaDmin max_delD, & !--- Max. deltaD for vapour (outlier if greater) deltalim max_delDtrac, & !--- Max. deltaD for tracers (defined in iso.def) max_delDsnow, & !--- Max. deltaD for snow (outlier if greater) min_Dexc, & !--- Min. value of D-excess (excess of H2[18]O dexcess_min max_Dexc, & !--- Max. value " " w/resp. to H[2]HO) dexcess_max min_O17exc, & !--- Min. value of O17-excess (excess of H2[17]O O17excess_bas max_O17exc !--- Max. value " " w/resp. to H2[18]O) O17excess_haut LOGICAL, SAVE :: checkO17 !--- Flag to trigger delta(H2[17]O) checking O17_verif INTEGER, SAVE :: max_O17nlev !--- nlevmaxO17 !$OMP THREADPRIVATE(min_delD, max_delDtrac, min_Dexc, min_O17exc, checkO17) !$OMP THREADPRIVATE(max_delD, max_delDsnow, max_Dexc, max_O17exc, max_O17nlev) ! logical bidouille_anti_divergence ! .TRUE. => xt relaxed to q to avoid errors accumulation leading to divergence ! Moved to wateriso2.h, rzad at execution !=== QUANTITIES DEPENDING ON THE ISOTOPES FAMILIES (LENGTH: nbIso) ; SET UP FIRST TIME ONLY IN "isoCheck_init" TYPE(r), ALLOCATABLE, SAVE :: dmin(:),dmax(:),&!--- Thresholds for correct delta for each isotope (len: isotopes families nb) emin(:),emax(:) !--- Thresholds for correct excess for each isotope (len: isotopes families nb) TYPE(i), ALLOCATABLE, SAVE :: iMajr(:) !--- Indices of the major isotopes (to be compared with parent concentration) CHARACTER(LEN=6), ALLOCATABLE, SAVE :: & sIsoCheck(:) !--- Activated checking routines keywords list ('npbem' ; 'x' if not activated) !=== VALUES TO BE SET UP BEFORE CALLING ELEMENTAL FUNCTIONS REAL, SAVE :: epsRel, epsAbs, epsPos !--- NEEDED BY "isEq???". Set before each checkEqual call REAL, SAVE :: min_bnd, max_bnd !--- NEEDED BY "isBounded". Set before each checkBound call REAL, SAVE :: tnat0 !--- NEEDED BY "r2delta" and "delta2r". Set before each delta* call REAL, SAVE :: slope !--- NEEDED BY "isoExcess". Set using "set_isoParams" CHARACTER(LEN=3), SAVE :: lawType !--- " " INTEGER, SAVE :: iso_ref !--- " " ! variables de vérifications real errmax ! erreur maximale en absolu. parameter (errmax=abs_err) real errmax_sol ! erreur maximale en absolu. parameter (errmax_sol=5e-7) real errmaxrel ! erreur maximale en relatif autorisée ! parameter (errmaxrel=1e10) parameter (errmaxrel=rel_err) real borne ! valeur maximale que n'importe quelle variable peut ! atteindre, en val abs; utile pour vérif que pas NAN parameter (borne=max_val) real, save :: deltalim ! deltalim est le maximum de deltaD qu'on ! autorise dans la vapeur, au-dela, en suppose que c'est abérrant. ! dans le liquide, on autorise deltalim*faccond. !$OMP THREADPRIVATE(deltalim) ! parameter (deltalim=1e10) ! parameter (deltalim=300.0) ! maintenant défini dans iso.def real, save :: deltalimtrac ! max de deltaD dans les traceurs, défini dans iso.def !$OMP THREADPRIVATE(deltalimtrac) real, save :: deltalim_snow ! deltalim est le maximum de deltaD qu'on ! autorise dans la neige, au-dela, en suppose que c'est abérrant. !$OMP THREADPRIVATE(deltalim_snow) ! parameter (deltalim_snow=500.0) ! initalisé dans iso_init real, save :: deltaDmin !$OMP THREADPRIVATE(deltaDmin) ! parameter (deltaDmin=-900.0) ! maintentant, défini dans iso.def real, save :: o17excess_bas,o17excess_haut ! bornes inf et sup de l'O17-excess ! parameter(o17excess_bas=-200.0,o17excess_haut=120) !$OMP THREADPRIVATE(o17excess_bas,o17excess_haut) integer nlevmaxO17 !$OMP THREADPRIVATE(nlevmaxO17) logical, save :: O17_verif !$OMP THREADPRIVATE(O17_verif) ! parameter (O17_verif=.true.) real, save :: dexcess_min,dexcess_max !$OMP THREADPRIVATE(dexcess_min,dexcess_max) ! real faccond ! dans le liquide, on autorise R(deltalim)*faccond. ! parameter (faccond=1.1) ! logical bidouille_anti_divergence ! si true, alors on fait un ! ! rappel régulier des xt4 vers les q pour ! !éviter accumulation d'érreurs et divergences ! parameter (bidouille_anti_divergence=.true.) ! parameter (bidouille_anti_divergence=.false.) ! bidouille_anti_divergence a été déplacé dans wateriso2.h et est lue à l'execution real deltaDfaible ! deltaD en dessous duquel la vapeur est tellement faible !que on peut accepter que la remise à l'équilibre du sol avec ! cette vapeur donne des deltaDevap aberrants. parameter (deltaDfaible=low_delD) real faible_evap ! faible évaporation: on est plus laxiste !pour les deltaD aberrant dans le cas de l'évap venant d'orchidee parameter (faible_evap=min_delDevap) real Tmin_verif parameter (Tmin_verif=min_T) ! en K real Tmax_verif parameter (Tmax_verif=max_T) ! en K ! subroutines de vérifications génériques, à ne pas modifier CONTAINS LOGICAL FUNCTION checkiIso(iIso, sub, iFamily) RESULT(lerr) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: sub INTEGER, OPTIONAL, INTENT(IN) :: iFamily CHARACTER(LEN=maxlen) :: subn INTEGER :: isav LOGICAL :: ll subn = 'checkiIso'; IF(PRESENT(sub)) subn = TRIM(sub)//TRIM(subn) IF(PRESENT(iFamily)) THEN iSav = ixIso; lerr = isoSelect(iFamily); IF(lerr) RETURN lerr = iIso <= 0 .OR. iIso > isotope%niso ll = isoSelect(iSav) ELSE lerr = iIso <= 0 .OR. iIso > isotope%niso END IF CALL msg('wrong isotopic index iIso = '//TRIM(num2str(iIso))//' (must be >0, <= '//TRIM(num2str(iIso))//')', subn, lerr) END FUNCTION checkiIso LOGICAL FUNCTION isoCheck_init() RESULT(lerr) CHARACTER(LEN=maxlen) :: isoFamily, modname='isoCheck_init' INTEGER :: iIso LOGICAL :: isoCheck CALL msg('Starting...', modname) ALLOCATE(dmin(nbIso), dmax(nbIso), emin(nbIso), emax(nbIso), sIsoCheck(nbIso), iMajr(nbIso)) CALL getin_p('isoCheck', isoCheck, .FALSE.) sIsoCheck = 'xxxxxx'; IF(isoCheck) sIsoCheck = 'npmdeb' !=== Single keyword for all isotopes families -> CAN BE GENERALIZED ! DO iIso = 1, nbIso isoFamily = isotope%parent SELECT CASE(isoFamily) !----------------------------------------------------------------------------------------------------------------------- CASE('H2O') !----------------------------------------------------------------------------------------------------------------------- iMajr(iIso)%i = [iso_eau] WRITE(*,*)'iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO = ', iso_eau, iso_O17, iso_O18, iso_HDO, iso_HTO iso_ref = iso_eau CALL getin_p('min_delD', min_delD, -900.0) CALL getin_p('max_delD', max_delD, 300.0) CALL getin_p('max_delDtrac', max_delDtrac, 300.0) max_delDsnow = 200.0 + max_delD WRITE(*,*)'max_delDsnow = ', max_delDsnow IF(iso_O17 > 0 .AND. iso_O18 > 0) THEN CALL getin_p('checkO17', checkO17, .FALSE.) CALL getin_p('min_O17exc', min_O17exc, -200.0) CALL getin_p('max_O17exc', max_O17exc, 120.0) CALL getin_p('max_O17nlev', max_O17nlev, 1000) END IF IF(iso_HDO > 0 .AND. iso_O18 > 0) THEN CALL getin_p('min_Dexc', min_Dexc, -100.0) CALL getin_p('max_Dexc', max_Dexc, 1000.0) END IF !=== HANDY VECTORS OF ACCEPTABLE LIMITS FOR delta ALLOCATE(dmin(iIso)%r(niso)); dmin(iIso)%r(:)=-max_val ALLOCATE(dmax(iIso)%r(niso)); dmax(iIso)%r(:)=+max_val IF(iso_HDO /= 0) THEN; dmin(iIso)%r(iso_HDO) = min_delD; dmax(iIso)%r(iso_HDO) = max_delD; END IF IF(iso_O17 /= 0) THEN; dmin(iIso)%r(iso_O17) =-max_val; dmax(iIso)%r(iso_O17) = max_val; END IF IF(iso_O18 /= 0) THEN; dmin(iIso)%r(iso_O18) = min_delD/2.0; dmax(iIso)%r(iso_O18) = max_delD/8.0; END IF !=== HANDY VECTORS OF ACCEPTABLE LIMITS FOR excess ALLOCATE(emin(iIso)%r(niso)); emin(iIso)%r(:)=-max_val ALLOCATE(emax(iIso)%r(niso)); emax(iIso)%r(:)=+max_val IF(iso_O17 /= 0) THEN; emin(iIso)%r(iso_O17) = min_O17exc; emax(iIso)%r(iso_O17) = max_O17exc; END IF IF(iso_O18 /= 0) THEN; emin(iIso)%r(iso_O18) = min_Dexc; emax(iIso)%r(iso_O18) = max_Dexc; END IF !----------------------------------------------------------------------------------------------------------------------- CASE DEFAULT lerr = .TRUE. CALL msg('routine not yet set up for isotopes family "'//TRIM(isoFamily)//'"', modname, lerr) !----------------------------------------------------------------------------------------------------------------------- END SELECT END DO END FUNCTION isoCheck_init !=============================================================================================================================== !=== ANCILLARY ROUTINES !=============================================================================================================================== !=== CHECK ISOTOPIC INDEX "iIso" AND DETERMINE FEW USEFULL QUANTITIES, FOR ELEMENTAL FUNCTIONS AND FOR ERROR MESSAGES: !=== lerr = set_isoParams(iIso, err_tag, sub[, mss][, phas][, iqIso][, iqPar]) !=== !=== ARGUMENTS: Intent Optional? Default value !=== * iIso Index to be checked IN MANDATORY !=== * err_tag Tag for error messages, appended at the beginning of "mss" IN OPTIONAL !=== * sub Calling sequence, appended with current routine name (":" separator) IN OPTIONAL !=== * mss Message in case absurd values are found OUT OPTIONAL !=== * phas Phase number, needed for the computation of "iqIso" and "iqPar" IN OPTIONAL !=== * iqIso Index in "qx" of the "iIso"th isotope OUT OPTIONAL !=== * iqPar Index in "qx" of the "iIso"th isotope parent OUT OPTIONAL !=== !=== NOTES: * SET CORRECT "slope", "iRef" AND "iso_ref" FOR ISOTOPIC EXCESS COMPUTATION. !=== SO FAR: H[2]HO-excess, H2[17]O-excess (with respect to H2[18]O) !=== * THE FOLLOWING DELTAs NEED TO BE COMPUTED: delta-H[2]HO, delta-H2[17]O, delta-H2[18]O !=============================================================================================================================== LOGICAL FUNCTION set_isoParams(iIso, err_tag, sub, mss, phas, iqIso, iqPar) RESULT(lerr) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), INTENT(IN) :: err_tag CHARACTER(LEN=maxlen), OPTIONAL, INTENT(OUT) :: sub, mss CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: phas INTEGER, OPTIONAL, INTENT(OUT) :: iqIso, iqPar INTEGER :: iq, iPha, ix LOGICAL :: lExc CHARACTER(LEN=maxlen) :: iname, subn subn = 'set_isoParams'; IF(PRESENT(sub)) subn = TRIM(sub)//TRIM(subn) lerr = checkiIso(iIso, subn); IF(lerr) RETURN IF(PRESENT(phas)) THEN iPha = INDEX(isoPhas, phas); lerr = iPha == 0 CALL msg('isotopes family "'//TRIM(isotope%parent)//'" has no "'//TRIM(phas)//'" phase', subn, lerr) ELSE iPha = 0; lerr = PRESENT(iqIso) .OR. PRESENT(iqPar) CALL msg('missing phase, needed to compute iqIso or iqPar', subn, lerr) END IF IF(lerr) RETURN ix = INDEX(sub, ':') SELECT CASE(sub(ix+1:LEN_TRIM(sub))) CASE('delta' ); iname = vname(iIso, 'd', iPha) IF(iPha /= 0) iq = iqIsoPha(iIso, iPha) IF(PRESENT(iqIso)) iqIso = iq IF(PRESENT(iqPar)) iqPar = tracers(iq)%iqParent CASE('excess'); iname = vname(iIso, 'e', iPha) IF (iIso == iso_HDO) THEN slope = slope_HDO; lawType = progr_HDO; iso_ref = iso_O18 ELSE IF(iIso == iso_O17) THEN slope = slope_O17; lawType = progr_O17; iso_ref = iso_O18 ELSE CALL msg('missing parameters for "'//TRIM(iname)//'" (iIso = '//TRIM(num2str(iIso))//'")', subn) lerr = .TRUE.; RETURN END IF END SELECT IF(PRESENT(mss)) THEN mss = 'absurd '//TRIM(iname)//' values found' IF(err_tag /= '') mss = TRIM(err_tag)//':'//TRIM(mss) END IF END FUNCTION set_isoParams !=============================================================================================================================== !=== DETERMINE THE NAME OF THE QUANTITY COMPUTED DEPENDING ON "typ": "d"elta, "e"xcess OR "r"atio. !=============================================================================================================================== CHARACTER(LEN=maxlen) FUNCTION vName(iIso, typ, iPha, iZon) RESULT(nam) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), INTENT(IN) :: typ INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon INTEGER :: ix, ip IF(iIso == 0) THEN IF(typ == 'd') nam = 'delta' IF(typ == 'e') nam = 'excess' IF(typ == 'r') nam = 'ratio' ELSE ip = 0; IF(PRESENT(iPha)) ip = iPha ix = iIso; IF(PRESENT(iZon)) ix = itZonIso(iZon, iIso) IF(ip /= 0 ) nam = tracers(iqIsoPha(ix, iPha))%name IF(ip == 0 ) nam = isoName(ix) IF(typ == 'd') nam = 'delta-'//TRIM(nam) IF(typ == 'e') nam = TRIM(nam)//'-excess' IF(typ == 'r') nam = TRIM(nam)//'-ratio' END IF END FUNCTION vName !=============================================================================================================================== !=== SAME AS vName, FOR SEVERAL INDICES "iIso" !=============================================================================================================================== FUNCTION vNames(iIso, typs, iPha, iZon) RESULT(nam) INTEGER, INTENT(IN) :: iIso(:) CHARACTER(LEN=*), INTENT(IN) :: typs INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon CHARACTER(LEN=maxlen) :: nam(SIZE(iIso)) INTEGER :: it DO it = 1, SIZE(nam) nam(it) = vName(iIso(it), typs(it:it), iPha, iZon) END DO END FUNCTION vNames !=============================================================================================================================== !=== LIST OF TWO VARIABLES FOR DELTA CASE: ratio(iIso), delta(iIso) !=============================================================================================================================== FUNCTION vNamDe(iIso, iPha, iZon) RESULT(nam) INTEGER, INTENT(IN) :: iIso INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon CHARACTER(LEN=maxlen) :: nam(2) nam = vNames([iIso, iso_ref], 'rd', iPha, iZon) END FUNCTION vNamDe !=============================================================================================================================== !=== LIST OF THREE VARIABLES FOR EXCESS CASE: delta(iIso), delta(iParent), excess(iIso) !=============================================================================================================================== FUNCTION vNamEx(iIso, iPha, iZon) RESULT(nam) INTEGER, INTENT(IN) :: iIso INTEGER, OPTIONAL, INTENT(IN) :: iPha, iZon CHARACTER(LEN=maxlen) :: nam(3) nam = vNames([iIso, iso_ref, iIso], 'dde', iPha, iZon) END FUNCTION vNamEx !=============================================================================================================================== !=== GET THE VARIABLES NAMES LIST (DEFAULT VALUE OR FROM tracers(:)%name, isoName(:) OR OPTIONAL ARGUMENT nam(:)) !=============================================================================================================================== LOGICAL FUNCTION gettNam(nq, sub, tnam, nam) RESULT(lerr) INTEGER, INTENT(IN) :: nq(:) CHARACTER(LEN=*), INTENT(IN) :: sub CHARACTER(LEN=*), ALLOCATABLE, INTENT(OUT) :: tnam(:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: nam(:) INTEGER :: lq, ln, rk rk = SIZE(nq); lq = SIZE(nq, DIM=rk) IF(PRESENT(nam)) THEN; ln = SIZE(nam) lerr = ALL([1, lq] /= ln) CALL msg('SIZE(q,'//TRIM(num2str(rk))//')='//TRIM(num2str(lq))//' /= SIZE(nam)='//num2str(ln), sub, lerr); IF(lerr) RETURN tnam = nam ELSE tnam = ['q'] IF(lq == ntiso) tnam = isoName(:) IF(lq == nqtot) tnam = tracers(:)%name END IF END FUNCTION gettNam !=============================================================================================================================== !=============================================================================================================================== !=== BASIC ELEMENTAL FUNCTIONS !=== FUNCTION isNaN (a) Not a Number (NaNs) detection: TRUE if |a| > max_val !=== FUNCTION isNeg (a) Positiveness detection: TRUE IF a < epsPos !=== FUNCTION isEqAbs(a,b) (1) Absolute equality checking: TRUE if |a-b| < epsAbs !=== FUNCTION isEqRel(a,b) (1) Relative equality checking: TRUE if |a-b|/max(|a|,|b|,1e-18) < epsRel !=== FUNCTION isEqual(a,b) (1) Equality checking (both criteria) !=== FUNCTION delta2ratio(ratioIso) (2) Delta function from a ratio value !=== FUNCTION ratio2delta(ratioIso) (2) Ratio function from a delta value !=== FUNCTION isoExcess(delIso, delRef) (3) Excess function from two delta values !=== !=== NOTES: (1): epsPos and epsAbs must be defined before calling !=== (2): tnat0 must be defined before calling !=== (3): slope and lawType must be defined before calling !=== !=============================================================================================================================== ELEMENTAL LOGICAL FUNCTION isNaN(a) RESULT(lout) !--- IS "a" AN OUTLIER, ie: a < -max_val or a > max_val ? REAL, INTENT(IN) :: a lout = a < -max_val .OR. max_val < a END FUNCTION isNaN !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL LOGICAL FUNCTION isBounded(a) RESULT(lout) !--- IS "a" BOUNDED, ie: min_bnd < a < max_bnd ? REAL, INTENT(IN) :: a lout= min_bnd < a .AND. a < max_bnd END FUNCTION isBounded !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL LOGICAL FUNCTION isNeg(a) RESULT(lout) !--- IS "a" NEGATIVE OR NOT, ie: a ≤ epsPos ? REAL, INTENT(IN) :: a lout = a <= epsPos END FUNCTION isNeg !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL LOGICAL FUNCTION isPos(a) RESULT(lout) !--- IS "a" POSITIVE OR NOT, ie: a ≥ epsPos ? REAL, INTENT(IN) :: a lout = a >= epsPos END FUNCTION isPos !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL LOGICAL FUNCTION isEqAbs(a,b) RESULT(lout) !--- IS "a" ABSOLUTELY EQ TO "b", ie: |a-b|<eps ? REAL, INTENT(IN) :: a, b lout = abs(a-b) < epsAbs END FUNCTION isEqAbs !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL LOGICAL FUNCTION isEqRel(a,b) RESULT(lout) !--- IS "a" RELATIVELY EQ TO "b", ie: |a-b|/max(|a|,|b|,1e-18)<eps ? REAL, INTENT(IN) :: a, b lout = ABS(a-b)/MAX(ABS(a),ABS(b),1e-18) < epsRel END FUNCTION isEqRel !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL LOGICAL FUNCTION isEqual(a,b) RESULT(lout) !--- IS "a" EQUAL TO "b", DUAL REL. AND ABS. CRITERIA REAL, INTENT(IN) :: a, b lout = isEqAbs(a,b); IF(lout) lout=isEqRel(a,b) END FUNCTION isEqual !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL REAL FUNCTION ratio2delta(ratioIso) RESULT(delta) !--- CONVERT AN ISOTOPIC RATIO INTO A DELTA REAL, INTENT(IN) :: ratioIso delta = (ratioIso / tnat0 - 1.0) * 1000.0 ! /!\ tnat0 must be defined !!! END FUNCTION ratio2delta !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL REAL FUNCTION delta2ratio(delta) RESULT(ratioIso) !--- CONVERT A DELTA INTO AN ISOTOPIC RATIO REAL, INTENT(IN) :: delta ratioIso = (delta / 1000.0 + 1.0) * tnat0 ! /!\ tnat0 must be defined !!! END FUNCTION delta2ratio !------------------------------------------------------------------------------------------------------------------------------- ELEMENTAL REAL FUNCTION isoExcess(deltaIso, deltaRef) RESULT(excess) !--- COMPUTE EXCESS USING TWO ISOTOPIC RATIOS REAL, INTENT(IN) :: deltaIso, deltaRef ! /!\ slope and lawType must be defined !!! SELECT CASE(lawType) CASE('lin'); excess = deltaIso - slope * deltaRef CASE('log'); excess = 1.e6*(LOG(1.+deltaIso/1000.) - slope *LOG(1.+deltaRef/1000.)) END SELECT END FUNCTION isoExcess !------------------------------------------------------------------------------------------------------------------------------- !=============================================================================================================================== !=== R = delta2r(iIso, delta[, lerr]): convert an anomaly "delta" into a ratio "R" !=== ARGUMENTS: Intent Optional? !=== * iIso Isotope index in "isoName" INPUT MANDATORY !=== * delta delta ratio for isotope of index "iIso" INPUT MANDATORY !=== * lerr Error code OUTPUT OPTIONAL !=== * R Isotopic ratio (output value, rank 1 to 3) OUTPUT MANDATORY !=============================================================================================================================== REAL FUNCTION delta2r_0D(iIso, delta, lerr) RESULT(R) INTEGER, INTENT(IN) :: iIso REAL, INTENT(IN) :: delta LOGICAL, OPTIONAL, INTENT(OUT) :: lerr LOGICAL :: ler ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN tnat0 = tnat(iIso); R = delta2ratio(delta) END FUNCTION delta2r_0D !------------------------------------------------------------------------------------------------------------------------------- FUNCTION delta2r_1D(iIso, delta, lerr) RESULT(R) REAL, ALLOCATABLE :: R(:) INTEGER, INTENT(IN) :: iIso REAL, INTENT(IN) :: delta(:) LOGICAL, OPTIONAL, INTENT(OUT) :: lerr LOGICAL :: ler ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN tnat0 = tnat(iIso); R = delta2ratio(delta) END FUNCTION delta2r_1D !------------------------------------------------------------------------------------------------------------------------------- FUNCTION delta2r_2D(iIso, delta, lerr) RESULT(R) REAL, ALLOCATABLE :: R(:,:) INTEGER, INTENT(IN) :: iIso REAL, INTENT(IN) :: delta(:,:) LOGICAL, OPTIONAL, INTENT(OUT) :: lerr LOGICAL :: ler ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN tnat0 = tnat(iIso); R = delta2ratio(delta) END FUNCTION delta2r_2D !=============================================================================================================================== !=== delta = r2delta(iIso, R[, lerr]): convert a ratio "R" into an anomaly "delta" !=== ARGUMENTS: Intent Optional? !=== * iIso Isotope index in "isoName" INPUT MANDATORY !=== * R Isotopic ratio for isotope of index "iIso" INPUT MANDATORY !=== * lerr Error code OUTPUT OPTIONAL !=== * delta delta ratio (output value, rank 1 to 3) OUTPUT MANDATORY !=============================================================================================================================== REAL FUNCTION r2delta_0D(iIso, R, lerr) RESULT(delta) INTEGER, INTENT(IN) :: iIso REAL, INTENT(IN) :: R LOGICAL, OPTIONAL, INTENT(OUT) :: lerr LOGICAL :: ler ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN tnat0 = tnat(iIso); delta = ratio2delta(R) END FUNCTION r2delta_0D !------------------------------------------------------------------------------------------------------------------------------- FUNCTION r2delta_1D(iIso, R, lerr) RESULT(delta) REAL, ALLOCATABLE :: delta(:) INTEGER, INTENT(IN) :: iIso REAL, INTENT(IN) :: R(:) LOGICAL, OPTIONAL, INTENT(OUT) :: lerr LOGICAL :: ler ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN tnat0 = tnat(iIso); delta = ratio2delta(R) END FUNCTION r2delta_1D !------------------------------------------------------------------------------------------------------------------------------- FUNCTION r2delta_2D(iIso, R, lerr) RESULT(delta) REAL, ALLOCATABLE :: delta(:,:) INTEGER, INTENT(IN) :: iIso REAL, INTENT(IN) :: R(:,:) LOGICAL, OPTIONAL, INTENT(OUT) :: lerr LOGICAL :: ler ler = checkiIso(iIso, 'delta2R'); IF(PRESENT(lerr)) lerr = ler; IF(ler) RETURN tnat0 = tnat(iIso); delta = ratio2delta(R) END FUNCTION r2delta_2D !=============================================================================================================================== !=============================================================================================================================== !=============================================================================================================================== !=== CHECK WHETHER INPUT FIELD "q" CONTAINS Not A Number (NaN) VALUES OR NOT !=== !=== lerr = checkNaN(q[, err_tag][, subn][, nam][, nmax]) !=== !=== ARGUMENTS: intent optional? default value !=== * q Concentration, scalar or array of rank 1 to 3 (checkNaN) INPUT MANDATORY !=== * err_tag Error message to display if NaNs are found INPUT OPTIONAL "NaNs found" !=== * subn Calling subroutine name INPUT OPTIONAL "checkNaN" !=== * nam Name(s) of the tracers (last index) INPUT OPTIONAL "q" !=== * nmax Maximum number of outliers values to compute for each tracer INPUT OPTIONAL All values !=============================================================================================================================== LOGICAL FUNCTION checkNaN_0D(q, err_tag, nam, subn) RESULT(lerr) REAL, INTENT(IN) :: q CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub ='checkNaN'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) lerr = gettNam([1], sub, tnm, [nam]); IF(lerr) RETURN lerr = dispOutliers( [isNaN(q)], [q], [1], tag, tnm, sub) END FUNCTION checkNaN_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkNaN_1D(q, err_tag, nam, subn, nmax) RESULT(lerr) REAL, INTENT(IN) :: q(:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub ='checkNaN'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkNaN_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkNaN_2D(q, err_tag, nam, subn, nmax) RESULT(lerr) REAL, INTENT(IN) :: q(:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub ='checkNaN'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkNaN_2D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkNaN_3D(q, err_tag, nam, subn, nmax) RESULT(lerr) REAL, INTENT(IN) :: q(:,:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='NaNs found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub ='checkNaN'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isNaN(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkNaN_3D !=============================================================================================================================== !=============================================================================================================================== !=============================================================================================================================== !=== CHECK WHETHER INPUT FIELD "q" HAVE ONLY POSITIVE VALUES OR NOT !=== !=== lerr = checkPos(q[, err_tag][, nam][, subn][, nmax][, tny]) Check for negative values in "q" . !=== !=== ARGUMENTS: Intent Optional? Default value !=== * q Concentration, scalar or array of rank 1 to 3 (checkNaN) INPUT MANDATORY !=== * err_tag Error message to display if NaNs are found INPUT OPTIONAL "Negative values found" !=== * nam Name(s) of the tracers (last index) INPUT OPTIONAL "q" !=== * subn Calling subroutine name INPUT OPTIONAL "checkPositive" !=== * nmax Maximum number of outliers values to compute for each tracer INPUT OPTIONAL All values !=== * tny Value (-) under which a "q" is considered to be negative INPUT OPTIONAL small !=============================================================================================================================== LOGICAL FUNCTION checkPos_0D(q, err_tag, nam, subn, tny) RESULT(lerr) REAL, INTENT(IN) :: q CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn REAL, OPTIONAL, INTENT(IN) :: tny CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub = 'checkPos'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) epsPos = -small; IF(PRESENT(tny )) epsPos = -tny lerr = gettNam(SHAPE(q), sub, tnm, [nam]); IF(lerr) RETURN lerr = dispOutliers( [isNeg(q)], [q], [1], tag, tnm, sub) END FUNCTION checkPos_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkPos_1D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr) REAL, INTENT(IN) :: q(:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: tny CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub = 'checkPos'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) epsPos = -small; IF(PRESENT(tny )) epsPos = -tny lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkPos_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkPos_2D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr) REAL, INTENT(IN) :: q(:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: tny CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub = 'checkPos'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) epsPos = -small; IF(PRESENT(tny )) epsPos = -tny lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkPos_2D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkPos_3D(q, err_tag, nam, subn, nmax, tny) RESULT(lerr) REAL, INTENT(IN) :: q(:,:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: tny CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='Negative values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub = 'checkPos'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) epsPos = -small; IF(PRESENT(tny )) epsPos = -tny lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isNeg(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkPos_3D !=============================================================================================================================== !=============================================================================================================================== !=============================================================================================================================== !=== CHECK WHETHER INPUT FIELD "q" HAVE VALUES INSIDE GIVEN BOUNDS OR NOT (IE "REASONNABLE" VALUES) !=== !=== lerr = checkBnd(q[, err_tag][, nam][, subn][, nmax][, qmin][, qmax]) !=== !=== ARGUMENTS: Intent Optional? Default value !=== * q Tested field, scalar or array of rank 1 to 3 INPUT MANDATORY !=== * err_tag Error message to display if out-of-bounds are found INPUT OPTIONAL "Negative values found" !=== * nam Name(s) of the tracers (last index) INPUT OPTIONAL "q" !=== * subn Calling subroutine name INPUT OPTIONAL "checkBnd" !=== * nmax Maximum number of outliers values to compute for each tracer INPUT OPTIONAL All values !=== * qmin Lower bound. Correct values must be > lBnd INPUT OPTIONAL 0 !=== * qmax Upper bound. Correct values must be < uBnd INPUT OPTIONAL max_val !=============================================================================================================================== LOGICAL FUNCTION checkBnd_0D(q, err_tag, nam, subn, qmin, qmax) RESULT(lerr) REAL, INTENT(IN) :: q CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam, subn REAL, OPTIONAL, INTENT(IN) :: qmin, qmax CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub = 'checkBnd'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) min_bnd = 0.; IF(PRESENT(qmin)) min_bnd = qmin max_bnd = max_val; IF(PRESENT(qmax)) max_bnd = qmax lerr = gettNam(SHAPE(q), sub, tnm, [nam]); IF(lerr) RETURN tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)' lerr = dispOutliers( [isBounded(q)], [q], [1], tag, tnm, sub) END FUNCTION checkBnd_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkBnd_1D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr) REAL, INTENT(IN) :: q(:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn REAL, OPTIONAL, INTENT(IN) :: qmin, qmax INTEGER, OPTIONAL, INTENT(IN) :: nmax CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub = 'checkBnd'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) min_bnd = 0.; IF(PRESENT(qmin)) min_bnd = qmin max_bnd = max_val; IF(PRESENT(qmax)) max_bnd = qmax tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)' lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkBnd_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkBnd_2D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr) REAL, INTENT(IN) :: q(:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn REAL, OPTIONAL, INTENT(IN) :: qmin, qmax INTEGER, OPTIONAL, INTENT(IN) :: nmax CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub = 'checkBnd'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) min_bnd = 0.; IF(PRESENT(qmin)) min_bnd = qmin max_bnd = max_val; IF(PRESENT(qmax)) max_bnd = qmax tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)' lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkBnd_2D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkBnd_3D(q, err_tag, nam, subn, nmax, qmin, qmax) RESULT(lerr) REAL, INTENT(IN) :: q(:,:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(:), subn REAL, OPTIONAL, INTENT(IN) :: qmin, qmax INTEGER, OPTIONAL, INTENT(IN) :: nmax CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) tag ='Absurd values found'; IF(PRESENT(err_tag)) THEN; IF(TRIM(err_tag)/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF sub = 'checkBnd'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub) min_bnd = 0.; IF(PRESENT(qmin)) min_bnd = qmin max_bnd = max_val; IF(PRESENT(qmax)) max_bnd = qmax tag = TRIM(tag)//' (must be in ]'//TRIM(num2str(min_bnd))//', '//TRIM(num2str(max_bnd))//'[)' lerr = gettNam(SHAPE(q), sub, tnm, nam); IF(lerr) RETURN lerr = dispOutliers(PACK(isBounded(q),.TRUE.), PACK(q,.TRUE.), SHAPE(q), tag, tnm, sub, nmax) END FUNCTION checkBnd_3D !=============================================================================================================================== !=============================================================================================================================== !=== CHECK WHETHER FIELDS "a" AND "b" ARE EQUAL OR NOT (ABSOLUTELY, THE RELATIVELY) !=== !=== lerr = checkEqu(a, b[, err_tag][, nam][, subn][, nmax][, absErr][, relErr]) !=== !=== ARGUMENTS: Intent Optional? Default value !=== * a, b Fields, rank 0 to 3 INPUT MANDATORY !=== * err_tag Error message to display if fields are not matching INPUT OPTIONAL "mismatching values" !=== * nam Name(s) of the tracers (last index) INPUT OPTIONAL "a","b" or from isoNames !=== * subn Calling subroutine name INPUT OPTIONAL "checkEqu" !=== * nmax Maximum number of outliers values to compute for each tracer INPUT OPTIONAL All values !=== * absErr Maximum value for |q-r| to consider "q" and "q" as equal INPUT OPTIONAL abs_err !=== * relErr Maximum value for |q-r|/max(|q|,|r|,1e-18) " " INPUT OPTIONAL rel_err !=============================================================================================================================== LOGICAL FUNCTION checkEqu_0D(a, b, err_tag, nam, subn, absErr, relErr) RESULT(lerr) REAL, INTENT(IN) :: a, b CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn REAL, OPTIONAL, INTENT(IN) :: absErr, relErr CHARACTER(LEN=maxlen) :: tag, sub CHARACTER(LEN=maxlen), ALLOCATABLE :: tnm(:) REAL :: aErr, rErr tag = 'mismatching value:'; IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF tnm = ['a ','b ','a-b']; IF(PRESENT(nam )) tnm = nam !--- Variables names sub = 'checkEqu'; IF(PRESENT(subn )) sub = TRIM(subn)//':'//TRIM(sub) epsAbs = abs_err; IF(PRESENT(absErr )) epsAbs = absErr !--- Absolute error epsRel = rel_err; IF(PRESENT(relErr )) epsRel = relErr !--- Relative error lerr = dispOutliers([isEqual(a,b)], cat(a, b, a-b), [1], tag,tnm,sub) END FUNCTION checkEqu_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkEqu_1D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr) REAL, INTENT(IN) :: a(:), b(:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: absErr, relErr CHARACTER(LEN=maxlen) :: tag, tnm(3), sub REAL :: aErr, rErr tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//': '//TRIM(tag); END IF tnm = ['a ','b ','a-b']; IF(PRESENT(nam )) tnm = nam !--- Variables names sub = 'checkEqu'; IF(PRESENT(subn )) sub = TRIM(subn)//':'//TRIM(sub) epsAbs = abs_err; IF(PRESENT(absErr )) epsAbs = absErr !--- Absolute error epsRel = rel_err; IF(PRESENT(relErr )) epsRel = relErr !--- Relative error lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax) END FUNCTION checkEqu_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkEqu_2D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr) REAL, INTENT(IN) :: a(:,:), b(:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: absErr, relErr CHARACTER(LEN=maxlen) :: tag, tnm(3), sub='' REAL :: aErr, rErr tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//TRIM(tag); END IF tnm = ['a ','b ','a-b']; IF(PRESENT(nam )) tnm = nam !--- Variables names sub = 'checkEqu'; IF(PRESENT(subn )) sub = TRIM(subn)//':'//TRIM(sub) epsAbs = abs_err; IF(PRESENT(absErr )) epsAbs = absErr !--- Absolute error epsRel = rel_err; IF(PRESENT(relErr )) epsRel = relErr !--- Relative error lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax) END FUNCTION checkEqu_2D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkEqu_3D(a, b, err_tag, nam, subn, nmax, absErr, relErr) RESULT(lerr) REAL, INTENT(IN) :: a(:,:,:), b(:,:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, nam(3), subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: absErr, relErr CHARACTER(LEN=maxlen) :: tag, tnm(3), sub='' REAL :: aErr, rErr tag = 'mismatching values:';IF(PRESENT(err_tag)) THEN; IF(err_tag/='') tag = TRIM(err_tag)//TRIM(tag); END IF tnm = ['a ','b ','a-b']; IF(PRESENT(nam )) tnm = nam !--- Variables names sub = 'checkEqu'; IF(PRESENT(subn )) sub = TRIM(subn)//':'//TRIM(sub) epsAbs = abs_err; IF(PRESENT(absErr )) epsAbs = absErr !--- Absolute error epsRel = rel_err; IF(PRESENT(relErr )) epsRel = relErr !--- Relative error lerr = dispOutliers(PACK(isEqual(a,b),.TRUE.), cat(PACK(a,.TRUE.),PACK(b,.TRUE.),PACK(a-b,.TRUE.)),SHAPE(a), tag,tnm,sub,nmax) END FUNCTION checkEqu_3D !=============================================================================================================================== !=============================================================================================================================== !=== CHECK FOR MASS CONSERVATION, IE. WHETHER THE CONCENTRATION OF A TRACER IS EQUAL TO THE SUM OF ITS CHILDREN CONCENTRATIONS. !=== EXCEPTION (PROBABLY TO BE SUPPRESSED): GENERATION 1 TRACERS ARE COMPARED TO THE MAJOR NEXT GEBNERATION ISOTOPES ONLY !=== !=== lerr = checkMass([qt][, err_tag][, subn][, nmax][, absErr][, relErr]) !=== !=== ARGUMENTS: Meaning Intent Optional? Default value !=== * qt Tracers+isotopes+tags stacked along last index, rank 1 to 3 INPUT OPTIONAL qx (RANK 3 CASE ONLY !) !=== * err_tag Error message to display if fields are not matching INPUT OPTIONAL "mismatching values" !=== * subn Calling subroutine name INPUT OPTIONAL "checkEqu" !=== * nmax Maximum number of outliers values to compute for each tracer INPUT OPTIONAL All values !=== * absErr Maximum value for |q-r| to consider "q" and "q" as equal. INPUT OPTIONAL abs_err !=== * relErr Maximum value for |q-r|/max(|q|,|r|,1e-18) " " INPUT OPTIONAL rel_err !=== !=== REMARKS: !=== * The concentration of each tracer is compared to the sum of its childrens concentrations (they must match). !=== * In the isotopes case, only a subset of childrens is used (usually only the major contributor). !=== * "qt" last dimension size can either be: !=== - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes (single unknown phase) !=== - nqtot ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx"). !=============================================================================================================================== LOGICAL FUNCTION checkMassQ_0D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr) REAL, INTENT(IN) :: qt(:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: absErr, relErr CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3) REAL :: qs, qp INTEGER, ALLOCATABLE :: iqChld(:) INTEGER :: iIso, iPha, nqChld, iq, it tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'checkMass'; IF(PRESENT(subn )) sub = TRIM(sub)//':'//subn !--- Calling subroutine chain epsAbs = abs_err; IF(PRESENT(absErr )) epsAbs = absErr !--- Absolute error (for isEqual routine) epsRel = rel_err; IF(PRESENT(relErr )) epsRel = relErr !--- Relative error (for isEqual routine) IF(SIZE(qt) == ntiso+1) THEN iqChld = iMajr(it)%i pNam = isotope%parent DO it = 0, niso IF(it /= 0) iqChld = itZonIso(:,it) IF(it /= 0) pNam = isoName(it) qp = qt(it+1) qs = SUM(qt(1+iqChld), DIM=1) !--- Tracer compared to its major isotopes sNam = TRIM(strStack(isoName(iqChld),'+')) !--- Names of contributing childrens tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam) nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference' !--- Names for the table lerr = dispOutliers([isEqual(qp, qs)], cat(qp, qs, qp-qs), SHAPE(qt), tag, nam, sub, nmax); IF(lerr) RETURN END DO ELSE DO iq = 1, SIZE(qt) pNam = tracers(iq)%name !--- Current tracer name with phase nqChld = tracers(iq)%nqChildren !--- Number of childrens (next generation only) IF(nqChld == 0) CYCLE !--- No descendants iIso = tracers(iq)%iso_iGroup !--- Isotopes family index iPha = tracers(iq)%iso_iPhase !--- Index of current phase in isoPhases(:) IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha) !--- Indices of major contributors in qt IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld) !--- Indices of every contributor in qt qp = qt(iq) !--- Linearized (reprofiled) parent field qs = SUM(qt(iqChld)) !--- Sum of contributing child(ren) field(s) sNam = TRIM(strStack(tracers(iqChld)%name ,'+')) !--- Names of contributing child(ren) tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam) nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference' !--- Names for the table lerr = dispOutliers([isEqual(qp, qs)], cat(qp, qs, qp-qs), [1], tag, nam, sub, nmax); IF(lerr) RETURN END DO END IF END FUNCTION checkMassQ_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkMassQ_1D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr) REAL, TARGET, INTENT(IN) :: qt(:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: absErr, relErr CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3) REAL, DIMENSION(SIZE(qt(:,1))) :: qs, qp INTEGER, ALLOCATABLE :: iqChld(:) INTEGER :: iIso, iPha, nqChld, iq, it tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'checkMass'; IF(PRESENT(subn )) sub = TRIM(sub)//':'//subn !--- Calling subroutine chain epsAbs = abs_err; IF(PRESENT(absErr )) epsAbs = absErr !--- Absolute error (for isEqual routine) epsRel = rel_err; IF(PRESENT(relErr )) epsRel = relErr !--- Relative error (for isEqual routine) IF(SIZE(qt, DIM=2) == ntiso+1) THEN iqChld = iMajr(it)%i pNam = isotope%parent DO it = 0, niso IF(it /= 0) iqChld = itZonIso(:,it) IF(it /= 0) pNam = isoName(it) qp = PACK(qt(:,it+1), MASK=.TRUE.) qs = SUM(qt(:,1+iqChld), DIM=2) !--- Tracer compared to its major isotopes sNam = TRIM(strStack(isoName(iqChld),'+')) !--- Names of contributing childrens tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam) nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference' !--- Names for the table lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,1)), tag, nam, sub, nmax); IF(lerr) RETURN END DO ELSE DO iq = 1, SIZE(qt, DIM=2) pNam = tracers(iq)%name !--- Current tracer name with phase nqChld = tracers(iq)%nqChildren !--- Number of childrens (next generation only) IF(nqChld == 0) CYCLE !--- No descendants iIso = tracers(iq)%iso_iGroup !--- Isotopes family index iPha = tracers(iq)%iso_iPhase !--- Index of current phase in isoPhases(:) IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha) !--- Indices of major contributors in qt IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld) !--- Indices of every contributor in qt qp = PACK(qt(:,iq), MASK=.TRUE.) !--- Linearized (reprofiled) parent field qs = PACK(SUM(qt(:,iqChld), DIM=2), MASK=.TRUE.) !--- Sum of contributing child(ren) field(s) sNam = TRIM(strStack(tracers(iqChld)%name ,'+')) !--- Names of contributing child(ren) tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam) nam(1) = pnam; nam(2) = 'SUM'; nam(3) = 'difference' !--- Names for the table lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,1)), tag, nam, sub, nmax); IF(lerr) RETURN END DO END IF END FUNCTION checkMassQ_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkMassQ_2D(qt, err_tag, subn, nmax, absErr, relErr) RESULT(lerr) REAL, TARGET, INTENT(IN) :: qt(:,:,:) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: absErr, relErr CHARACTER(LEN=maxlen) :: tag, sub, sNam, pNam, nam(3) REAL, DIMENSION(SIZE(qt(:,:,1))) :: qs, qp INTEGER, ALLOCATABLE :: iqChld(:) INTEGER :: iIso, iPha, nqChld, iq, it tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'checkMass'; IF(PRESENT(subn )) sub = TRIM(sub)//':'//subn !--- Calling subroutine chain epsAbs = abs_err; IF(PRESENT(absErr )) epsAbs = absErr !--- Absolute error (for isEqual routine) epsRel = rel_err; IF(PRESENT(relErr )) epsRel = relErr !--- Relative error (for isEqual routine) IF(SIZE(qt, DIM=3) == ntiso+1) THEN iqChld = iMajr(it)%i pNam = isotope%parent DO it = 0, niso IF(it /= 0) iqChld = itZonIso(:,it) IF(it /= 0) pNam = isoName(it) qp = PACK(qt(:,:,it+1), MASK=.TRUE.) qs = PACK(SUM(qt(:,:,1+iqChld), DIM=3), MASK=.TRUE.) !--- Tracer compared to its major isotopes sNam = TRIM(strStack(isoName(iqChld),'+')) !--- Names of contributing childrens tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam) nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference' !--- Names for the table lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,:,1)), tag, nam, sub, nmax); IF(lerr) RETURN END DO ELSE DO iq = 1, SIZE(qt, DIM=3) pNam = tracers(iq)%name !--- Current tracer name with phase nqChld = tracers(iq)%nqChildren !--- Number of childrens (next generation only) IF(nqChld == 0) CYCLE !--- No descendants iIso = tracers(iq)%iso_iGroup !--- Isotopes family index iPha = tracers(iq)%iso_iPhase !--- Index of current phase in isoPhases(:) IF(iIso /= 0) iqChld = iqIsoPha(iMajr(iIso)%i, iPha) !--- Indices of major contributors in qt IF(iIso == 0) iqChld = tracers(iq)%iqDescen(1:nqChld) !--- Indices of every contributor in qt qp = PACK(qt(:,:,iq), MASK=.TRUE.) !--- Linearized (reprofiled) parent field qs = PACK(SUM(qt(:,:,iqChld), DIM=3), MASK=.TRUE.) !--- Sum of contributing child(ren) field(s) sNam = TRIM(strStack(tracers(iqChld)%name ,'+')) !--- Names of contributing child(ren) tag = TRIM(tag)//' checking difference between '//TRIM(pNam)//' and '//TRIM(sNam) nam(1) = pNam; nam(2) = 'SUM'; nam(3) = 'difference' !--- Names for the table lerr = dispOutliers(isEqual(qp, qs), cat(qp, qs, qp-qs), SHAPE(qt(:,:,1)), tag, nam, sub, nmax); IF(lerr) RETURN END DO END IF END FUNCTION checkMassQ_2D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION checkMassQx_2D(err_tag, subn, nmax, absErr, relErr) RESULT(lerr) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: absErr, relErr lerr = checkMassQ_2D(qx, err_tag, subn, nmax, absErr, relErr) END FUNCTION checkMassQx_2D !=============================================================================================================================== !=============================================================================================================================== !=============================================================================================================================== !=== COMPUTE THE ISOTOPIC ANOMALY "DELTA" AND CHECK THE VALUES !=== !=== lerr = deltaR(rIso, iIso[, err_tag][, subn][, nmax][, delMax][, delIso] [, lCheck]) !=== lerr = deltaQ(qt, iIso[, err_tag][, subn][, nmax][, delMax][, delIso][, Phas][, qmin][, lCheck]) !=== lerr = deltaQ( iIso[, err_tag][, subn][, nmax][, delMax][, delIso][, Phas][, qmin][, lCheck]) !=== !=== ARGUMENTS: Meaning Intent Optional? Default Present: !=== * rIso Field of rank 0 to 3 INPUT MANDATORY CASE 1 ONLY !=== * qt Stacked fields, rank 1 to 3 (last index for species) INPUT MANDATORY CASE 2 ONLY !=== * iIso Index of tested isotope in "isoName" INPUT MANDATORY !=== * err_tag Error message to display if fields are not matching INPUT OPTIONAL "??????????" !=== * subn Calling subroutine name INPUT OPTIONAL "delta[R]" !=== * nmax Maximum number of lines to print for outliers INPUT OPTIONAL <all lines> RANK>0 ONLY !=== * delMax Maximum value for the isotopic delta INPUT OPTIONAL dmax(iIso) !=== * delIso Isotopic anomaly delta for isotope iIso OUTPUT OPTIONAL !=== * phas Phase INPUT OPTIONAL CASE 2 ONLY !=== * qmin Values lower than this threshold are not checked INPUT OPTIONAL small CASE 2 ONLY !=== * lCheck Trigger checking routines with tables for strage values INPUT OPTIONAL .FALSE. !=== !=== REMARKS: !=== * In the version 2 and 3 of the routine, values are checked only if q>=qmin !=== * In the version 2, the size of the last dimension of "qt" can either be: !=== - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes. !=== - nqtot ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx"). !=============================================================================================================================== LOGICAL FUNCTION deltaR_0D(rIso, iIso, err_tag, subn, delMax, delIso, lCheck) RESULT(lerr) REAL, INTENT(IN) :: rIso INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn REAL, OPTIONAL, INTENT(IN) :: delMax REAL, OPTIONAL, INTENT(OUT) :: delIso LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss REAL :: dmn, dmx, dIso LOGICAL :: lc IF(niso == 0) RETURN tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'delta'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger dmn = dMin(ixIso)%r(iIso) lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN !--- Check the index + define few vars mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')' tnat0 = tnat(iIso); dIso = ratio2delta(rIso) !=== Compute delta(isotope iIso in phase iPhas) IF(PRESENT(delIso)) delIso = dIso IF(.NOT.lc) RETURN lerr = dispOutliers([dIso<dmn .OR. dIso>dmx], cat(rIso, dIso), [1], mss, vNamDe(iIso), sub) END FUNCTION deltaR_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION deltaR_1D(rIso, iIso, err_tag, subn, nmax, delMax, delIso, lCheck) RESULT(lerr) REAL, INTENT(IN) :: rIso(:) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax REAL, OPTIONAL, INTENT(OUT) :: delIso(SIZE(rIso,1)) LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss REAL :: dmn, dmx, dIso(SIZE(rIso,1)) LOGICAL :: lc IF(niso == 0) RETURN tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'delta'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger dmn = dMin(ixIso)%r(iIso) lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN !--- Check the index + define few vars mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')' tnat0 = tnat(iIso); dIso = ratio2delta(rIso) !=== Compute delta(isotope iIso in phase iPhas) IF(PRESENT(delIso)) delIso = dIso IF(.NOT.lc) RETURN lerr = dispOutliers(PACK(dIso<dmn .OR. dIso>dmx,.TRUE.), cat(rIso, dIso), SHAPE(dIso), mss, vNamDe(iIso), sub, nmax) END FUNCTION deltaR_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION deltaR_2D(rIso, iIso, err_tag, subn, nmax, delMax, delIso, lCheck) RESULT(lerr) REAL, INTENT(IN) :: rIso(:,:) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax REAL, OPTIONAL, INTENT(OUT) :: delIso(SIZE(rIso,1),SIZE(rIso,2)) LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss REAL :: dmn, dmx, emx, dIso(SIZE(rIso,1),SIZE(rIso,2)) LOGICAL :: lc IF(niso == 0) RETURN tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'delta'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger dmn = dMin(ixIso)%r(iIso) lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN !--- Check the index + define few vars mss = TRIM(mss)//' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')' tnat0 = tnat(iIso); dIso = ratio2delta(rIso) IF(PRESENT(delIso)) delIso = dIso IF(.NOT.lc) RETURN lerr = dispOutliers( [dIso<dmn.OR.dIso>dmx], cat(PACK(rIso,.TRUE.),PACK(dIso,.TRUE.)), SHAPE(dIso),mss,vNamDe(iIso),sub,nmax) END FUNCTION deltaR_2D !=============================================================================================================================== LOGICAL FUNCTION deltaQ_0D(qt, iIso, err_tag, subn, delMax, delIso, phas, qmin, lCheck) RESULT(lerr) REAL, DIMENSION(:), INTENT(IN) :: qt INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, phas REAL, OPTIONAL, INTENT(IN) :: delMax, qmin REAL, OPTIONAL, INTENT(OUT) :: delIso LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m REAL :: dmn, dmx, qmn, rIso, dIso LOGICAL :: lc, lq INTEGER :: iqIso, iqPar, iZon, iPha, ii tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'delta'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger pha = 'g'; IF(PRESENT(phas )) pha = phas qmn = small; IF(PRESENT(qmin )) qmn = qmin !--- q negligible if q<=dmn (default: small) lq = SIZE(qt, DIM=1) == niso+1 !=== CHECK PHASE IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0 CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr) IF(lerr) RETURN dmn = dMin(ixIso)%r(iIso) m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')' !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone) ii = iIso DO iZon = 0, nzone IF(iZon /= 0) ii = itZonIso(iZon, iIso) DO iPha = 1, nphas p = isoPhas(iPha:iPha) lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN !--- Ckeck the index + define few vars mss = TRIM(mss)//TRIM(m) IF( lq) rIso = qt(iIso+1)/qt(1) !--- "qt(1+niso)": parent, then isotopes IF(.NOT.lq) rIso = qt(iqIso)/qt(iqPar) !--- Fields taken from a "qt" similar to "qx" tnat0 = tnat(iIso); dIso = ratio2delta(rIso) !=== Compute delta(isotope iIso in phase "p") IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso !--- Return delta(iIso, phase=p) IF(.NOT.lc) CYCLE lerr = dispOutliers( [qt(iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], cat(rIso, dIso), [1], mss, vNamde(iIso), sub) IF(lerr) RETURN END DO IF(.NOT.lc) EXIT END DO END FUNCTION deltaQ_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION deltaQ_1D(qt, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr) REAL, DIMENSION(:,:), INTENT(IN) :: qt INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, phas INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax, qmin REAL, OPTIONAL, INTENT(OUT) :: delIso(SIZE(qt,1)) LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m REAL :: dmn, dmx, qmn REAL, DIMENSION(SIZE(qt,1)) :: rIso, dIso LOGICAL :: lc, lq INTEGER :: iqIso, iqPar, iZon, iPha, ii tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'delta'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger pha = 'g'; IF(PRESENT(phas )) pha = phas qmn = small; IF(PRESENT(qmin )) qmn = qmin !--- q negligible if q<=dmn (default: small) lq = SIZE(qt, DIM=2) == niso+1 !=== CHECK PHASE IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0 CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr) IF(lerr) RETURN dmn = dMin(ixIso)%r(iIso) m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')' !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone) ii = iIso DO iZon = 0, nzone IF(iZon /= 0) ii = itZonIso(iZon, iIso) DO iPha = 1, nphas p = isoPhas(iPha:iPha) lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN !--- Ckeck the index + define few vars mss = TRIM(mss)//TRIM(m) IF( lq) rIso = qt(:,iIso+1)/qt(:,1) !--- "qt(1+niso)": parent, then isotopes IF(.NOT.lq) rIso = qt(:,iqIso)/qt(:,iqPar) !--- Fields taken from a "qt" similar to "qx" tnat0 = tnat(iIso); dIso = ratio2delta(rIso) !=== Compute delta(iIso, phase=p) IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso !--- Return delta(iIso, phase=p) IF(.NOT.lc) CYCLE lerr = dispOutliers( [qt(:,iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], & cat(PACK(rIso,.TRUE.), PACK(dIso,.TRUE.)), SHAPE(dIso), mss, vNamde(iIso), sub, nmax) IF(lerr) RETURN END DO IF(.NOT.lc) EXIT END DO END FUNCTION deltaQ_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION deltaQ_2D(qt, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr) REAL, DIMENSION(:,:,:), INTENT(IN) :: qt INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, phas INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax, qmin REAL, OPTIONAL, INTENT(OUT) :: delIso(SIZE(qt,1),SIZE(qt,2)) LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m REAL :: dmn, dmx, qmn REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)) :: rIso, dIso LOGICAL :: lc, lq INTEGER :: iqIso, iqPar, iZon, iPha, ii tag = ''; IF(PRESENT(err_tag)) tag = err_tag !--- Tag for error message sub = 'delta'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger pha = 'g'; IF(PRESENT(phas )) pha = phas qmn = small; IF(PRESENT(qmin )) qmn = qmin !--- q negligible if q<=dmn (default: small) lq = SIZE(qt, DIM=3) == niso+1 !=== CHECK PHASE IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0 CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr) IF(lerr) RETURN dmn = dMin(ixIso)%r(iIso) m = ' (dMin='//TRIM(num2str(dmn))//', dMax = '//TRIM(num2str(dmx))//')' !=== COMPUTE AND CHECK THE DELTA VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone) ii = iIso DO iZon = 0, nzone IF(iZon /= 0) ii = itZonIso(iZon, iIso) DO iPha = 1, nphas p = isoPhas(iPha:iPha) lerr = set_isoParams(ii,tag,sub,mss,p,iqIso,iqPar); IF(lerr) RETURN !--- Ckeck the index + define few vars mss = TRIM(mss)//' '//TRIM(m) IF( lq) rIso = qt(:,:,iIso+1)/qt(:,:,1) !--- "qt(1+niso)": parent, then isotopes IF(.NOT.lq) rIso = qt(:,:,iqIso)/qt(:,:,iqPar) !--- Fields taken from a "qt" similar to "qx" tnat0 = tnat(iIso); dIso = ratio2delta(rIso) !=== Compute delta(iIso, phase=p) IF(iZon == 0 .AND. p == pha .AND. PRESENT(delIso)) delIso = dIso !--- Return delta(iIso, phase=p) IF(.NOT.lc) CYCLE lerr = dispOutliers( [qt(:,:,iqPar)>=qmn .AND. (dIso<dmn .OR. dIso>dmx)], & cat(PACK(rIso,.TRUE.), PACK(dIso,.TRUE.)), SHAPE(dIso), mss, vNamDe(iIso), sub, nmax) IF(lerr) RETURN END DO IF(.NOT.lc) EXIT END DO END FUNCTION deltaQ_2D !=============================================================================================================================== LOGICAL FUNCTION deltaQx_2D(iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) RESULT(lerr) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, phas INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax, qmin REAL, OPTIONAL, INTENT(OUT) :: delIso(SIZE(qx,1),SIZE(qx,2)) LOGICAL, OPTIONAL, INTENT(IN) :: lCheck lerr = deltaQ_2D(qx, iIso, err_tag, subn, nmax, delMax, delIso, phas, qmin, lCheck) END FUNCTION deltaQx_2D !=============================================================================================================================== !=============================================================================================================================== !=============================================================================================================================== !=== COMPUTE THE ISOTOPIC EXCESS AND CHECK THE VALUE (D-excess and O17-excess so far) !=== !=== lerr = excess(rIso,rRef,iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef] [,lCheck]) !=== lerr = excess(qt, iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef][,phas][,qmin][,lCheck]) !=== lerr = excess( iIso[,err_tag][,subn][,nmax][,delMax][,excMax][,delIso][,excIso][,delRef][,phas][,qmin][,lCheck]) !=== !=== ARGUMENTS: Intent Optional? Default value !=== * rIso,rRef Considered isotope + reference ratios: field of rank 0 to 3 INPUT MANDATORY (CASE 1 ONLY) !=== * qt Stacked fields, rank 1 to 3 (last index for species) INPUT MANDATORY (CASE 2 ONLY) !=== * iIso Index of tested species INPUT MANDATORY !=== * err_tag Error tag sent from the calling routine INPUT OPTIONAL none !=== * subn Calling subroutine name INPUT OPTIONAL "excess" !=== * nmax Maximum number of lines to print for outliers INPUT OPTIONAL <all lines> (RANK>0 ONLY) !=== * delMax Maximum value for the isotopic delta INPUT OPTIONAL !=== * excMax Maximum value for the isotopic excess INPUT OPTIONAL dmax(iIso) !=== * delIso Isotopic anomaly delta for isotope iIso OUTPUT OPTIONAL !=== * excIso Isotopic excess for isotope iIso OUTPUT OPTIONAL !=== * delRef Isotopic anomaly delta for isotope iIso reference OUTPUT OPTIONAL !=== * phas Phase name (one element of "isoPhas") INPUT MANDATORY (CASE 2 ONLY) !=== * qmin Values lower than this threshold are not checked INPUT OPTIONAL small (CASE 2 ONLY) !=== !=== REMARKS: !=== * In the version 2 and 3 of the routine, values are checked only if q>=qmin !=== * In the version 2, the size of the last dimension of "qt" can either be: !=== - ntiso+1 ; then the parent tracer is at position 1, followed by the isotopes. !=== - nqtot ; then the tracers are described by the desciptor "tracers(:)" (same caracteristics as "qx"). !=============================================================================================================================== LOGICAL FUNCTION excessR_0D(rIso, rRef, iIso, err_tag, subn, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr) REAL, INTENT(IN) :: rIso, rRef INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn REAL, OPTIONAL, INTENT(IN) :: delMax, excMax REAL, OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss REAL :: deIso, deRef, exIso, dmx, drx, emn, emx LOGICAL :: lc tag = ''; IF(PRESENT(err_tag)) tag = err_tag sub = 'excess'; IF(PRESENT(subn )) sub = TRIM(subn)//':'//TRIM(sub) !--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax !--- Maximum excess lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN !--- Ckeck the index + define few vars emn = eMin(ixIso)%r(iIso) drx = dMax(ixIso)%r(iso_ref) mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')' lerr = deltaR_0D(rIso, iIso, tag, sub, dmx, deIso, lCheck) & !=== COMPUTE deltaIso AND deltaRef .OR. deltaR_0D(rRef, iso_ref, tag, sub, drx, deRef, lCheck) IF(lerr) RETURN exIso = isoExcess(deIso, deRef) !=== COMPUTE THE EXCESS IF(PRESENT(delIso)) delIso = deIso IF(PRESENT(excIso)) excIso = exIso IF(PRESENT(delRef)) delRef = deRef IF(.NOT.lc) RETURN lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), [1], mss, vNamEx(iIso), sub) END FUNCTION excessR_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION excessR_1D(rIso, rRef, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr) REAL, DIMENSION(:), INTENT(IN) :: rIso, rRef INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax, excMax REAL, DIMENSION(SIZE(rIso,1)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss, m REAL, DIMENSION(SIZE(rIso,1)) :: deIso, deRef, exIso REAL :: dmx, drx, emn, emx LOGICAL :: lc tag = ''; IF(PRESENT(err_tag)) tag = err_tag sub = 'excess'; IF(PRESENT(subn )) sub = TRIM(subn)//':'//TRIM(sub) !--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax !--- Maximum excess lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN !--- Ckeck the index + define few vars emn = eMin(ixIso)%r(iIso) drx = dMax(ixIso)%r(iso_ref) mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')' lerr = deltaR_1D(rIso, iIso, tag, sub, nmax, dmx, deIso, lCheck) & !=== COMPUTE deltaIso AND deltaRef .OR. deltaR_1D(rRef, iso_ref, tag, sub, nmax, drx, deRef, lCheck) IF(lerr) RETURN exIso = isoExcess(deIso, deRef) !=== COMPUTE THE EXCESS IF(PRESENT(delIso)) delIso = deIso IF(PRESENT(excIso)) excIso = exIso IF(PRESENT(delRef)) delRef = deRef IF(.NOT.lc) RETURN lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), SHAPE(exIso), mss, vNamEx(iIso), sub, nmax) END FUNCTION excessR_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION excessR_2D(rIso, rRef, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, lCheck) RESULT(lerr) REAL, DIMENSION(:,:), INTENT(IN) :: rIso, rRef INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax, excMax REAL, DIMENSION(SIZE(rIso,1),SIZE(rIso,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss, m REAL, DIMENSION(SIZE(rIso,1),SIZE(rIso,2)) :: deIso, deRef, exIso REAL :: dmx, drx, emn, emx LOGICAL :: lc tag = ''; IF(PRESENT(err_tag)) tag = err_tag sub = 'excess'; IF(PRESENT(subn )) sub = TRIM(subn)//':'//TRIM(sub) !--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax !--- Maximum excess lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger lerr = set_isoParams(iIso, tag, sub, mss); IF(lerr) RETURN !--- Ckeck the index + define few vars emn = eMin(ixIso)%r(iIso) drx = dMax(ixIso)%r(iso_ref) mss = TRIM(mss)//'(eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')' lerr = deltaR_2D(rIso, iIso, tag, sub, nmax, dmx, deIso, lCheck) & !=== COMPUTE deltaIso AND deltaRef .OR. deltaR_2D(rRef, iso_ref, tag, sub, nmax, drx, deRef, lCheck) IF(lerr) RETURN exIso = isoExcess(deIso, deRef) !=== COMPUTE THE EXCESS IF(PRESENT(delIso)) delIso = deIso IF(PRESENT(excIso)) excIso = exIso IF(PRESENT(delRef)) delRef = deRef IF(.NOT.lc) RETURN lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(PACK(deIso,.TRUE.), PACK(deRef,.TRUE.), PACK(exIso,.TRUE.)), & SHAPE(exIso), TRIM(mss)//' '//TRIM(m), vNamEx(iIso), sub, nmax) END FUNCTION excessR_2D !=============================================================================================================================== LOGICAL FUNCTION excessQ_0D(qt, iIso, err_tag, subn, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck) RESULT(lerr) REAL, INTENT(IN) :: qt(:) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, phas REAL, OPTIONAL, INTENT(IN) :: delMax, excMax, qmin REAL, OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m REAL :: deIso, exIso, deRef REAL :: dmx, drx, emn, emx, qmn INTEGER :: iZon, iPha, ii LOGICAL :: lc tag = ''; IF(PRESENT(err_tag)) tag = err_tag sub = 'excess'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax !--- Maximum excess pha = 'g'; IF(PRESENT(phas )) pha = phas qmn = small; IF(PRESENT(qmin )) qmn = qmin !--- q negligible if q<=dmn (default: small) lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger !=== CHECK PHASE lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0 CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr) IF(lerr) RETURN emn = eMin(ixIso)%r(iIso) drx = dMax(ixIso)%r(iso_ref) m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')' !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone) ii = iIso DO iZon = 0, nzone IF(iZon /= 0) ii = itZonIso(iZon, iIso) lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN !--- Check ii, set iso_ref + excess params mss = TRIM(mss)//TRIM(m) DO iPha = 1, nphas p = isoPhas(iPha:iPha) !--- Current phase lerr = deltaQ_0D(qt, ii, tag, sub, dmx, deIso, p, qmn, lCheck) & !=== COMPUTE deltaIso AND deltaRef .OR. deltaQ_0D(qt, iso_ref, tag, sub, drx, deRef, p, qmn, lCheck) exIso = isoExcess(deIso, deRef) !=== COMPUTE THE EXCESS IF(iZon == 0 .AND. p == pha) THEN IF(PRESENT(delIso)) delIso = deIso !--- Return delta(iIso, phase=p) IF(PRESENT(delRef)) delRef = deRef !--- Return delta(iso_ref,phase=p) IF(PRESENT(excIso)) excIso = exIso !--- Return excess(iIso, phase=p) END IF IF(.NOT.lc) CYCLE lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), [1], mss, vNamEx(ii, iPha, iZon), sub) IF(lerr) RETURN END DO IF(.NOT.lc) EXIT END DO END FUNCTION excessQ_0D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION excessQ_1D(qt, iIso, err_tag, subn, nmax, delMax,excMax,delIso,excIso,delRef, phas, qmin, lCheck) RESULT(lerr) REAL, INTENT(IN) :: qt(:,:) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, phas INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax, excMax, qmin REAL, DIMENSION(SIZE(qt,1)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m REAL, DIMENSION(SIZE(qt,1)) :: deIso, exIso, deRef REAL :: dmx, drx, emn, emx, qmn INTEGER :: iZon, iPha, ii LOGICAL :: lc tag = ''; IF(PRESENT(err_tag)) tag = err_tag sub = 'excess'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax !--- Maximum excess pha = 'g'; IF(PRESENT(phas )) pha = phas qmn = small; IF(PRESENT(qmin )) qmn = qmin !--- q negligible if q<=dmn (default: small) lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger !=== CHECK PHASE lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0 CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr) IF(lerr) RETURN emn = eMin(ixIso)%r(iIso) drx = dMax(ixIso)%r(iso_ref) m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')' !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone) ii = iIso DO iZon = 0, nzone IF(iZon /= 0) ii = itZonIso(iZon, iIso) lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN !--- Check ii, set iso_ref + excess params mss = TRIM(mss)//TRIM(m) DO iPha = 1, nphas p = isoPhas(iPha:iPha) !--- Current phase lerr = deltaQ_1D(qt, ii, tag,sub,nmax, dmx, deIso,p,qmn,lCheck) & !=== COMPUTE deltaIso AND deltaRef .OR. deltaQ_1D(qt, iso_ref, tag,sub,nmax, drx, deRef,p,qmn,lCheck) exIso = isoExcess(deIso, deRef) !=== COMPUTE THE EXCESS IF(iZon == 0 .AND. p == pha) THEN IF(PRESENT(delIso)) delIso = deIso !--- Return delta(iIso, phase=p) IF(PRESENT(delRef)) delRef = deRef !--- Return delta(iso_ref,phase=p) IF(PRESENT(excIso)) excIso = exIso !--- Return excess(iIso, phase=p) END IF IF(.NOT.lc) CYCLE lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(deIso, deRef, exIso), & SHAPE(exIso), mss, vNamEx(ii, iPha, iZon), sub, nmax) IF(lerr) RETURN END DO IF(.NOT.lc) EXIT END DO END FUNCTION excessQ_1D !------------------------------------------------------------------------------------------------------------------------------- LOGICAL FUNCTION excessQ_2D(qt, iIso, err_tag, subn, nmax, delMax,excMax,delIso,excIso,delRef, phas, qmin, lCheck) RESULT(lerr) REAL, INTENT(IN) :: qt(:,:,:) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, phas INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax, excMax, qmin REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef LOGICAL, OPTIONAL, INTENT(IN) :: lCheck CHARACTER(LEN=maxlen) :: tag, sub, mss, p, pha, m REAL, DIMENSION(SIZE(qt,1),SIZE(qt,2)) :: deIso, exIso, deRef REAL :: dmx, drx, emn, emx, qmn INTEGER :: iZon, iPha, ii LOGICAL :: lc tag = ''; IF(PRESENT(err_tag)) tag = err_tag sub = 'excess'; IF(PRESENT(subn)) sub = TRIM(subn)//':'//TRIM(sub)!--- Calling subroutine name dmx = dMax(ixIso)%r(iIso); IF(PRESENT(delMax )) dmx = delMax !--- Maximum delta emx = eMax(ixIso)%r(iIso); IF(PRESENT(excMax )) emx = excMax !--- Maximum excess pha = 'g'; IF(PRESENT(phas )) pha = phas qmn = small; IF(PRESENT(qmin )) qmn = qmin !--- q negligible if q<=dmn (default: small) lc = .FALSE.; IF(PRESENT(lCheck )) lc = lCheck !--- Checking routines execution trigger !=== CHECK PHASE lerr = .FALSE.; IF(PRESENT(delIso)) lerr = INDEX(isoPhas, pha) == 0 CALL msg('missing phase "'//TRIM(pha)//'" for isotopes family "'//TRIM(isotope%parent)//'"', sub, lerr) IF(lerr) RETURN emn = eMin(ixIso)%r(iIso) drx = dMax(ixIso)%r(iso_ref) m = ' (eMin='//TRIM(num2str(emn))//', eMax = '//TRIM(num2str(emx))//')' !=== COMPUTE AND CHECK THE EXCESS VALUES OF THE ISOTOPE iIso (iZon==0) AND ITS TAGGING TRACERS (iZon = 1:nzone) ii = iIso DO iZon = 0, nzone IF(iZon /= 0) ii = itZonIso(iZon, iIso) lerr = set_isoParams(ii, tag, sub, mss); IF(lerr) RETURN !--- Check ii, set iso_ref + excess params mss = TRIM(mss)//TRIM(m) DO iPha = 1, nphas p = isoPhas(iPha:iPha) !--- Current phase lerr = deltaQ_2D(qt, ii, tag,sub,nmax, dmx, deIso,p,qmn,lCheck) & !=== COMPUTE deltaIso AND deltaRef .OR. deltaQ_2D(qt, iso_ref, tag,sub,nmax, drx, deRef,p,qmn,lCheck) exIso = isoExcess(deIso, deRef) !=== COMPUTE THE EXCESS IF(iZon == 0 .AND. p == pha) THEN IF(PRESENT(delIso)) delIso = deIso !--- Return delta(iIso, phase=p) IF(PRESENT(delRef)) delRef = deRef !--- Return delta(iso_ref,phase=p) IF(PRESENT(excIso)) excIso = exIso !--- Return excess(iIso, phase=p) END IF IF(.NOT.lc) CYCLE lerr = dispOutliers( [exIso<emn .OR. exIso>emx], cat(PACK(deIso,.TRUE.), PACK(deRef,.TRUE.), PACK(exIso,.TRUE.)), & SHAPE(exIso), mss, vNamEx(ii, iPha, iZon), sub, nmax) IF(lerr) RETURN END DO IF(.NOT.lc) EXIT END DO END FUNCTION excessQ_2D !=============================================================================================================================== LOGICAL FUNCTION excessQx_2D(iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck) RESULT(lerr) INTEGER, INTENT(IN) :: iIso CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, phas INTEGER, OPTIONAL, INTENT(IN) :: nmax REAL, OPTIONAL, INTENT(IN) :: delMax, excMax, qmin REAL, DIMENSION(SIZE(qx,1),SIZE(qx,2)), OPTIONAL, INTENT(OUT) :: delIso, excIso, delRef LOGICAL, OPTIONAL, INTENT(IN) :: lCheck lerr = excessQ_2D(qx, iIso, err_tag, subn, nmax, delMax, excMax, delIso, excIso, delRef, phas, qmin, lCheck) END FUNCTION excessQx_2D !=============================================================================================================================== !=============================================================================================================================== !=== GENERAL PURPOSE CHECKING ROUTINE: lerr = checkTrac(err_tag[, subn][, nmax][, sCheck]) !=== !=== ARGUMENTS: Intent Optional? Default value !=== * err_tag Error tag sent from the calling routine INPUT OPTIONAL none !=== * subn Calling subroutine name INPUT OPTIONAL "checkTrac" !=== * nmax Maximum number of lines to print for outliers INPUT OPTIONAL 32 !=== * sCheck Triggers for verifications INPUT OPTIONAL sIsoChech !=== !=== REMARKS: !=== Tunable thresholds available in the individual routines (delMax, excMax, qmin, tny, absErr, relErr) have not been kept !=== as optional arguments in this routine, because the adequate values are tracer or isotope-dependent. !=== For special usages, a specific version can be written, or individual routines with special thresholds can be called. !=============================================================================================================================== LOGICAL FUNCTION checkTrac(err_tag, subn, nmax, sCheck) RESULT(lerr) CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: err_tag, subn, sCheck INTEGER, OPTIONAL, INTENT(IN) :: nmax INTEGER :: nmx, i, iPha, iIso, iq LOGICAL :: lNan, lPos, lMas, lDel, lExc, lBnd CHARACTER(LEN=maxlen) :: tag, sub, chk CHARACTER(LEN=maxlen), ALLOCATABLE :: tnam(:) !!!!! GERER TNAM IF(niso == 0) RETURN tag = ''; IF(PRESENT(err_tag)) tag = err_tag sub = 'checkTrac'; IF(PRESENT(subn )) sub = TRIM(sub)//TRIM(subn) nmx = 32; IF(PRESENT(nmax )) nmx = nmax chk = sIsoCheck(ixIso) ; IF(PRESENT(sCheck )) chk = sCheck SELECT CASE(isotope%parent) CASE('H2O') DO iPha = 1, isotope%nphas DO iIso = 1, niso iq = iqIsoPha(iIso,iPha); tnam = tracers(iq)%name!; qmin = tracers(iq)%qmin; qmax = tracers(iq)%qmax IF(chk(1:1) == 'n') lerr = checkNaN (qx(:,:,iq),tag,tnam,sub,nmx); IF(lerr) RETURN IF(chk(2:2) == 'p') lerr = checkPos (qx(:,:,iq),tag,tnam,sub,nmx); IF(lerr) RETURN ! tny IF(chk(3:3) == 'm') lerr = checkMass( tag,sub, nmx); IF(lerr) RETURN ! absErr, relErr !IF(chk(6:6) == 'b') lerr = checkBnd (qx(:,:,iq),tag,tnam,sub,nmx,qmin,qmax); IF(lerr) RETURN END DO iIso = iso_HDO; IF(chk(4:4) == 'd') lerr = delta(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN iIso = iso_O18; IF(chk(4:4) == 'd') lerr = delta(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN iIso = iso_HDO; IF(chk(5:5) == 'e') lerr = excess(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN iIso = iso_O17; IF(chk(5:5) == 'e') lerr = excess(iIso, tag, sub, nmx, lCheck=iIso/=0); IF(lerr) RETURN END DO END SELECT END FUNCTION checkTrac !=============================================================================================================================== SUBROUTINE iso_verif_init() use ioipsl_getin_p_mod, ONLY : getin_p use isotopes_mod, ONLY: iso_O17, iso_O18, iso_HDO implicit none write(*,*) 'iso_verif_init 99: entree' deltalim=300.0 deltalim_snow=500.0 deltaDmin=-900.0 deltalimtrac=2000.0 O17_verif=.false. o17excess_bas=-200.0 o17excess_haut=120.0 dexcess_min=-100.0 dexcess_max=1000.0 call getin_p('deltalim',deltalim) deltalim_snow=deltalim+200.0 call getin_p('deltaDmin',deltaDmin) call getin_p('deltalimtrac',deltalimtrac) write(*,*) 'iso_O17,iso_O18,iso_HDO=',iso_O17,iso_O18,iso_HDO if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then call getin_p('O17_verif',O17_verif) call getin_p('o17excess_bas',o17excess_bas) call getin_p('o17excess_haut',o17excess_haut) call getin_p('nlevmaxO17',nlevmaxO17) endif if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then call getin_p('dexcess_min',dexcess_min) call getin_p('dexcess_max',dexcess_max) endif write(*,*) 'deltalim=',deltalim write(*,*) 'deltaDmin=',deltaDmin write(*,*) 'deltalimtrac=',deltalimtrac write(*,*) 'O17_verif=',O17_verif write(*,*) 'o17excess_bas=',o17excess_bas write(*,*) 'o17excess_haut=',o17excess_haut write(*,*) 'dexcess_min=',dexcess_min write(*,*) 'dexcess_max=',dexcess_max end SUBROUTINE iso_verif_init subroutine iso_verif_egalite(a,b,err_msg) implicit none ! compare a et b. Si pas egal à errmax près, on affiche message ! d'erreur et stoppe ! input: real a, b character*(*) err_msg ! message d''erreur à afficher ! local !integer iso_verif_egalite_choix_nostop if (iso_verif_egalite_choix_nostop(a,b,err_msg, & & errmax,errmaxrel).eq.1) then stop endif #ifdef ISOVERIF #else write(*,*) 'err_msg=',err_msg,': pourquoi verif?' stop #endif return end subroutine iso_verif_egalite !***************** function iso_verif_egalite_nostop(a,b,err_msg) implicit none ! compare a et b. Si pas egal à errmax près, on affiche message ! d'erreur et stoppe ! input: real a, b character*(*) err_msg ! message d''erreur à afficher !ouptut integer iso_verif_egalite_nostop ! local !integer iso_verif_egalite_choix_nostop iso_verif_egalite_nostop=iso_verif_egalite_choix_nostop & & (a,b,err_msg,errmax,errmaxrel) #ifdef ISOVERIF #else write(*,*) 'err_msg=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_egalite_nostop !************************************ subroutine iso_verif_aberrant(R,err_msg) use isotopes_mod, ONLY: ridicule, iso_HDO implicit none ! si le rapprot iso R est plus grand que deltaDlim, on affiche message ! d'erreur et stoppe ! input: real R character*(*) err_msg ! message d''erreur à afficher !real deltaD !integer iso_verif_aberrant_choix_nostop if (iso_verif_aberrant_choix_nostop(R,1.0,ridicule, & & deltalim,err_msg).eq.1) then stop endif #ifdef ISOVERIF if (.not.(iso_HDO.gt.0)) then write(*,*) 'iso86: err_msg=',err_msg,': pourquoi verif?' stop endif #else write(*,*) 'iso 90: err_msg=',err_msg,': pourquoi verif?' stop #endif end subroutine iso_verif_aberrant subroutine iso_verif_aberrant_encadre(R,err_msg) use isotopes_mod, ONLY: ridicule, iso_HDO implicit none ! si le rapprot iso R est plus grand que deltaDlim, on affiche message ! d'erreur et stoppe ! input: real R character*(*) err_msg ! message d''erreur à afficher !real deltaD !integer iso_verif_aberrant_enc_choix_nostop if (iso_verif_aberrant_enc_choix_nostop(R,1.0,ridicule, & & deltalim,err_msg).eq.1) then write(*,*) 'deltaD=',deltaD(R) call abort_physic('isotopes_verif_mod > iso_verif_aberrant_encadre',err_msg,1) !stop endif #ifdef ISOVERIF if (.not.(iso_HDO.gt.0)) then write(*,*) 'iso86: err_msg=',err_msg,': pourquoi verif?' stop endif #else write(*,*) 'iso 90: err_msg=',err_msg,': pourquoi verif?' stop #endif end subroutine iso_verif_aberrant_encadre !************************************ subroutine iso_verif_aberrant_choix(xt,q,qmin,deltaDmax,err_msg) use isotopes_mod, ONLY: iso_HDO implicit none ! si le rapprot iso R est plus grand que deltaDlim, on affiche message ! d'erreur et stoppe ! input: real xt,q,qmin,deltaDmax character*(*) err_msg ! message d''erreur à afficher !real deltaD ! locals !integer iso_verif_aberrant_choix_nostop if (iso_verif_aberrant_choix_nostop(xt,q,qmin, & & deltaDmax,err_msg).eq.1) then stop endif #ifdef ISOVERIF if (.not.(iso_HDO.gt.0)) then write(*,*) 'iso122: err_msg=',err_msg,': pourquoi verif?' stop endif #else write(*,*) 'iso126: err_msg=',err_msg,': pourquoi verif?' stop #endif end subroutine iso_verif_aberrant_choix !************************************ function iso_verif_aberrant_nostop(R,err_msg) use isotopes_mod, ONLY: ridicule,iso_HDO implicit none ! si le rapprot iso R est plus grand que deltaDlim, on affiche message ! d'erreur et stoppe ! input: real R character*(*) err_msg ! message d''erreur à afficher integer iso_verif_aberrant_nostop ! output: 1 si erreur, 0 sinon !real deltaD ! locals !integer iso_verif_aberrant_choix_nostop iso_verif_aberrant_nostop=iso_verif_aberrant_choix_nostop & & (R,1.0,ridicule,deltalim,err_msg) #ifdef ISOVERIF if (.not.(iso_HDO.gt.0)) then write(*,*) 'iso156: err_msg=',err_msg,': pourquoi verif?' stop endif #else write(*,*) 'iso160: err_msg=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_aberrant_nostop function iso_verif_aberrant_enc_nostop(R,err_msg) use isotopes_mod, ONLY: ridicule,iso_HDO implicit none ! si le rapprot iso R est plus grand que deltaDlim, on affiche message ! d'erreur et stoppe ! input: real R character*(*) err_msg ! message d''erreur à afficher integer iso_verif_aberrant_enc_nostop ! output: 1 si erreur, 0 sinon !real deltaD ! locals !integer iso_verif_aberrant_enc_choix_nostop iso_verif_aberrant_enc_nostop= & & iso_verif_aberrant_enc_choix_nostop & & (R,1.0,ridicule,deltalim,err_msg) #ifdef ISOVERIF if (.not.(iso_HDO.gt.0)) then write(*,*) 'iso156: err_msg=',err_msg,': pourquoi verif?' stop endif #else write(*,*) 'iso160: err_msg=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_aberrant_enc_nostop !************************************ function iso_verif_aberrant_choix_nostop(xt,q, & & qmin,deltaDmax,err_msg) use isotopes_mod, ONLY: iso_HDO implicit none ! si le rapprot iso R est plus grand que deltaDlim, on affiche message ! d'erreur et stoppe ! input: real xt,q,qmin,deltaDmax character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_aberrant_choix_nostop ! locals !real deltaD !integer iso_verif_noNaN_nostop iso_verif_aberrant_choix_nostop=0 #ifdef ISOVERIF if (iso_verif_noNaN_nostop(q,err_msg).eq.1) then write(*,*) 'q=',q iso_verif_aberrant_choix_nostop=1 endif if (iso_verif_noNaN_nostop(xt,err_msg).eq.1) then write(*,*) 'xt=',xt iso_verif_aberrant_choix_nostop=1 endif #endif if (q.gt.qmin) then if ((deltaD(xt/q).gt.deltaDmax).or. & & (deltaD(xt/q).lt.-borne).or. & & (xt.lt.-borne).or. & & (xt.gt.borne)) then write(*,*) 'erreur detectee par '// & & 'iso_verif_aberrant_choix_nostop:' write(*,*) err_msg write(*,*) 'q,deltaD=',q,deltaD(xt/q) write(*,*) 'deltaDmax=',deltaDmax iso_verif_aberrant_choix_nostop=1 if (abs(xt-q)/q.lt.errmax) then write(*,*) 'attention, n''a-t-on pas confondu'// & & ' iso_HDO et iso_eau dans l''appel à la verif?' endif endif endif #ifdef ISOVERIF if (.not.(iso_HDO.gt.0)) then write(*,*) 'iso205: err_msg=',err_msg,': pourquoi verif?' stop endif #else write(*,*) 'iso 209: err_msg=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_aberrant_choix_nostop function iso_verif_aberrant_enc_choix_nostop(xt,q, & & qmin,deltaDmax,err_msg) use isotopes_mod, ONLY: iso_HDO implicit none ! si le rapprot iso R est plus grand que deltaDlim, on affiche message ! d'erreur et stoppe ! input: real xt,q,qmin,deltaDmax character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_aberrant_enc_choix_nostop ! locals !real deltaD iso_verif_aberrant_enc_choix_nostop=0 if (q.gt.qmin) then if ((deltaD(xt/q).gt.deltaDmax).or. & & (deltaD(xt/q).lt.deltaDmin)) then write(*,*) 'erreur detectee par '// & & 'iso_verif_aberrant_choix_nostop:' write(*,*) err_msg write(*,*) 'q,deltaD=',q,deltaD(xt/q) iso_verif_aberrant_enc_choix_nostop=1 if (abs(xt-q)/q.lt.errmax) then write(*,*) 'attention, n''a-t-on pas confondu'// & & ' iso_HDO et iso_eau dans l''appel à la verif?' endif endif endif #ifdef ISOVERIF if (.not.(iso_HDO.gt.0)) then write(*,*) 'iso205: err_msg=',err_msg,': pourquoi verif?' stop endif #else write(*,*) 'iso 209: err_msg=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_aberrant_enc_choix_nostop !******************* subroutine iso_verif_aberrant_o17(R17,R18,err_msg) implicit none ! si l'O17-excess est aberrant, on affiche un message ! d'erreur et stoppe ! input: real R17,R18 character*(*) err_msg ! message d''erreur à afficher !real o17excess ! locals !integer iso_verif_aberrant_o17_nostop ! write(*,*) 'O17_verif=',O17_verif if (O17_verif) then if (iso_verif_aberrant_o17_nostop(R17,R18,err_msg) & & .eq.1) then stop endif endif !if (O17_verif) then #ifdef ISOVERIF #else write(*,*) 'err_msg=',err_msg,': pourquoi verif?' stop #endif return end subroutine iso_verif_aberrant_o17 !******************* function iso_verif_aberrant_o17_nostop(R17,R18,err_msg) USE isotopes_mod, ONLY: tnat,iso_O17,iso_O18 implicit none ! si l'O17-excess est aberrant, on affiche un message ! d'erreur et renvoit 1 ! input: real R17,R18 character*(*) err_msg ! message d''erreur à afficher !local !real o17excess ! output integer iso_verif_aberrant_o17_nostop if (O17_verif) then iso_verif_aberrant_o17_nostop=0 if ((o17excess(R17,R18).gt.o17excess_haut).or. & & (o17excess(R17,R18).lt.o17excess_bas)) then write(*,*) 'erreur detectee par iso_verif_aberrant_O17:' write(*,*) err_msg write(*,*) 'o17excess=',o17excess(R17,R18) write(*,*) 'deltaO17=',(R17/tnat(iso_o17)-1.0)*1000.0 write(*,*) 'deltaO18=',(R18/tnat(iso_O18)-1.0)*1000.0 ! attention, vérifier que la ligne suivante est bien activée iso_verif_aberrant_o17_nostop=1 endif endif !if (O17_verif) then #ifdef ISOVERIF #else write(*,*) 'err_msg=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_aberrant_o17_nostop !************************************ subroutine iso_verif_noNaN(x,err_msg) implicit none ! si x est NaN, on affiche message ! d'erreur et stoppe ! input: real x character*(*) err_msg ! message d''erreur à afficher ! locals !integer iso_verif_noNAN_nostop if (iso_verif_noNAN_nostop(x,err_msg).eq.1) then stop endif #ifdef ISOVERIF #else write(*,*) 'err_msg iso443=',err_msg,': pourquoi verif?' stop #endif end subroutine iso_verif_noNaN !************************************ function iso_verif_noNaN_nostop(x,err_msg) implicit none ! si x est NaN, on affiche message ! d'erreur et return 1 si erreur ! input: real x character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_noNAN_nostop if ((x.gt.-borne).and.(x.lt.borne)) then iso_verif_noNAN_nostop=0 else write(*,*) 'erreur detectee par iso_verif_nonNAN:' write(*,*) err_msg write(*,*) 'x=',x iso_verif_noNAN_nostop=1 endif #ifdef ISOVERIF #else write(*,*) 'err_msg iso 482=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_noNaN_nostop subroutine iso_verif_noNaN_vect2D(x,err_msg,ni,n,m) implicit none ! si x est NaN, on affiche message ! d'erreur et return 1 si erreur ! input: integer n,m,ni real x(ni,n,m) character*(*) err_msg ! message d''erreur à afficher ! output ! locals integer i,j,ixt do i=1,n do j=1,m do ixt=1,ni if ((x(ixt,i,j).gt.-borne).and. & & (x(ixt,i,j).lt.borne)) then else !if ((x(ixt,i,j).gt.-borne).and. write(*,*) 'erreur detectee par iso_verif_nonNAN:' write(*,*) err_msg write(*,*) 'x,ixt,i,j=',x(ixt,i,j),ixt,i,j stop endif !if ((x(ixt,i,j).gt.-borne).and. enddo !do ixt=1,ni enddo !do j=1,m enddo !do i=1,n #ifdef ISOVERIF #else write(*,*) 'err_msg iso525=',err_msg,': pourquoi verif?' stop #endif end subroutine iso_verif_noNaN_vect2D subroutine iso_verif_noNaN_vect1D(x,err_msg,ni,n) implicit none ! si x est NaN, on affiche message ! d'erreur et return 1 si erreur ! input: integer n,ni real x(ni,n) character*(*) err_msg ! message d''erreur à afficher ! output ! locals integer i,ixt do i=1,n do ixt=1,ni if ((x(ixt,i).gt.-borne).and. & & (x(ixt,i).lt.borne)) then else !if ((x(ixt,i,j).gt.-borne).and. write(*,*) 'erreur detectee par iso_verif_nonNAN:' write(*,*) err_msg write(*,*) 'x,ixt,i=',x(ixt,i),ixt,i stop endif !if ((x(ixt,i,j).gt.-borne).and. enddo !do ixt=1,ni enddo !do i=1,n #ifdef ISOVERIF #else write(*,*) 'err_msg iso525=',err_msg,': pourquoi verif?' stop #endif end subroutine iso_verif_noNaN_vect1D !************************ subroutine iso_verif_egalite_choix(a,b,err_msg,erabs,errel) implicit none ! compare a et b. Si pas egal, on affiche message ! d'erreur et stoppe ! pour egalite, on verifie erreur absolue et arreur relative ! input: real a, b real erabs,errel !erreur absolue et relative character*(*) err_msg ! message d''erreur à afficher ! locals !integer iso_verif_egalite_choix_nostop if (iso_verif_egalite_choix_nostop( & & a,b,err_msg,erabs,errel).eq.1) then stop endif #ifdef ISOVERIF #else write(*,*) 'err_msg=',err_msg,': pourquoi verif?' stop #endif end subroutine iso_verif_egalite_choix !************************ function iso_verif_egalite_choix_nostop & & (a,b,err_msg,erabs,errel) implicit none ! compare a et b. Si pas egal, on affiche message ! d'erreur et stoppe ! pour egalite, on verifie erreur absolue et arreur relative ! input: real a, b real erabs,errel !erreur absolue et relative character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_egalite_choix_nostop ! locals !integer iso_verif_noNaN_nostop iso_verif_egalite_choix_nostop=0 #ifdef ISOVERIF if (iso_verif_noNaN_nostop(a,err_msg).eq.1) then write(*,*) 'a=',a iso_verif_egalite_choix_nostop=1 endif if (iso_verif_noNaN_nostop(b,err_msg).eq.1) then write(*,*) 'b=',b iso_verif_egalite_choix_nostop=1 endif #endif if (abs(a-b).gt.erabs) then if (abs((a-b)/max(max(abs(b),abs(a)),1e-18)) & & .gt.errel) then write(*,*) 'erreur detectee par iso_verif_egalite:' write(*,*) err_msg write(*,*) 'a=',a write(*,*) 'b=',b write(*,*) 'erabs,errel=',erabs,errel iso_verif_egalite_choix_nostop=1 ! stop endif endif #ifdef ISOVERIF #else write(*,*) 'err_msg=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_egalite_choix_nostop !****************** subroutine iso_verif_positif(x,err_msg) use isotopes_mod, ONLY: ridicule implicit none ! si x<0, on plante. ! si très limite, on le met à 0. ! input: real x character*(*) err_msg ! message d''erreur à afficher ! locals !integer iso_verif_positif_choix_nostop if (iso_verif_positif_choix_nostop(x,ridicule,err_msg) & & .eq.1) then stop endif #ifdef ISOVERIF #else write(*,*) 'iso_verif 637: err_msg=',err_msg, & & ': pourquoi verif?' stop #endif end subroutine iso_verif_positif !****************** subroutine iso_verif_positif_vect(x,n,err_msg) use isotopes_mod, ONLY: ridicule implicit none ! si x<0, on plante. ! input: integer n real x(n) character*(*) err_msg ! message d''erreur à afficher ! locals integer i integer iso_verif_positif_nostop integer ifaux iso_verif_positif_nostop=0 do i=1,n if (x(i).lt.-ridicule) then iso_verif_positif_nostop=1 ifaux=i endif enddo if (iso_verif_positif_nostop.eq.1) then write(*,*) 'erreur detectee par iso_verif_positif_vect:' write(*,*) err_msg write(*,*) 'i,x=',ifaux,x(ifaux) stop endif end subroutine iso_verif_positif_vect subroutine iso_verif_positif_choix_vect(x,n,ridic,err_msg) implicit none ! si x<0, on plante. ! input: integer n real x(n) real ridic character*(*) err_msg ! message d''erreur à afficher ! locals integer i integer iso_verif_positif_nostop integer ifaux iso_verif_positif_nostop=0 do i=1,n if (x(i).lt.-ridic) then iso_verif_positif_nostop=1 ifaux=i endif enddo if (iso_verif_positif_nostop.eq.1) then write(*,*) 'erreur detectee par iso_verif_positif_choix_vect:' write(*,*) err_msg write(*,*) 'i,x=',ifaux,x(ifaux) stop endif end subroutine iso_verif_positif_choix_vect !****************** subroutine iso_verif_positif_strict(x,err_msg) implicit none ! si x<0, on plante. ! si très limite, on le met à 0. ! input: real x character*(*) err_msg ! message d''erreur à afficher ! locals !integer iso_verif_positif_strict_nostop if (iso_verif_positif_strict_nostop(x,err_msg) & & .eq.1) then stop endif end subroutine iso_verif_positif_strict !****************** function iso_verif_positif_strict_nostop(x,err_msg) implicit none ! si x<0, on plante. ! si très limite, on le met à 0. ! input: real x character*(*) err_msg ! message d''erreur à afficher* ! output integer iso_verif_positif_strict_nostop if (x.gt.0.0) then iso_verif_positif_strict_nostop=0 else write(*,*) 'erreur detectee par iso_verif_positif_strict:' write(*,*) err_msg write(*,*) 'x=',x iso_verif_positif_strict_nostop=1 ! stop endif return end function iso_verif_positif_strict_nostop !****************** subroutine iso_verif_positif_choix(x,ridic,err_msg) implicit none ! si x<0, on plante. ! si très limite, on le met à 0. ! input: real x real ridic character*(*) err_msg ! message d''erreur à afficher ! locals !integer iso_verif_positif_choix_nostop if (iso_verif_positif_choix_nostop(x,ridic,err_msg) & & .eq.1) then stop endif #ifdef ISOVERIF #else write(*,*) 'iso_verif 801: err_msg=',err_msg, & & ': pourquoi verif?' stop #endif end subroutine iso_verif_positif_choix !****************** function iso_verif_positif_nostop(x,err_msg) use isotopes_mod, ONLY: ridicule implicit none ! si x<0, on plante. ! si très limite, on le met à 0. ! input: real x character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_positif_nostop ! locals !integer iso_verif_positif_choix_nostop iso_verif_positif_nostop=iso_verif_positif_choix_nostop & & (x,ridicule,err_msg) #ifdef ISOVERIF #else write(*,*) 'iso_verif 837: err_msg=',err_msg, & & ': pourquoi verif?' stop #endif return end function iso_verif_positif_nostop !****************** function iso_verif_positif_choix_nostop(x,ridic,err_msg) implicit none ! si x<0, on plante. ! si très limite, on le met à 0. ! input: real x real ridic character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_positif_choix_nostop if (x.lt.-ridic) then write(*,*) 'erreur detectee par iso_verif_positif_nostop:' write(*,*) err_msg write(*,*) 'x=',x iso_verif_positif_choix_nostop=1 else ! x=max(x,0.0) iso_verif_positif_choix_nostop=0 endif #ifdef ISOVERIF #else write(*,*) 'iso_verif 877: err_msg=',err_msg, & & ': pourquoi verif?' stop #endif return end function iso_verif_positif_choix_nostop !************** subroutine iso_verif_O18_aberrant(Rd,Ro,err_msg) implicit none ! vérifie que: ! 1) deltaD et deltaO18 dans bonne gamme ! 2) dexcess dans bonne gamme ! input: real Rd,Ro character*(*) err_msg ! message d''erreur à afficher ! local !integer iso_verif_O18_aberrant_nostop if (iso_verif_O18_aberrant_nostop(Rd,Ro,err_msg).eq.1) then stop endif end subroutine iso_verif_O18_aberrant function iso_verif_O18_aberrant_nostop(Rd,Ro,err_msg) USE isotopes_mod, ONLY: tnat, iso_HDO, iso_O18 implicit none ! vérifie que: ! 1) deltaD et deltaO18 dans bonne gamme ! 2) dexcess dans bonne gamme ! input: real Rd,Ro character*(*) err_msg ! message d''erreur à afficher ! outputs integer iso_verif_O18_aberrant_nostop !locals real deltaD,deltao,dexcess deltaD=(Rd/tnat(iso_hdo)-1)*1000 deltao=(Ro/tnat(iso_O18)-1)*1000 dexcess=deltaD-8*deltao iso_verif_O18_aberrant_nostop=0 if ((deltaD.lt.deltaDmin).or.(deltao.lt.deltaDmin/2.0).or. & & (deltaD.gt.deltalim).or.(deltao.gt.deltalim/8.0).or. & & ((deltaD.gt.-500.0).and.((dexcess.lt.dexcess_min) & & .or.(dexcess.gt.dexcess_max)))) then write(*,*) 'erreur detectee par iso_verif_O18_aberrant:' write(*,*) err_msg write(*,*) 'delta180=',deltao write(*,*) 'deltaD=',deltaD write(*,*) 'Dexcess=',dexcess write(*,*) 'tnat=',tnat ! stop iso_verif_O18_aberrant_nostop=1 endif #ifdef ISOVERIF #else write(*,*) 'err_msg=',err_msg,': pourquoi verif?' stop #endif return end function iso_verif_O18_aberrant_nostop ! ********** function deltaD(R) USE isotopes_mod, ONLY: tnat,iso_HDO implicit none real R,deltaD if (iso_HDO.gt.0) then deltaD=(R/tnat(iso_HDO)-1)*1000.0 else write(*,*) 'iso_verif_egalite>deltaD 260: iso_HDO.gt.0=', & & iso_HDO.gt.0 endif return end function deltaD ! ********** function deltaO(R) USE isotopes_mod, ONLY: tnat,iso_O18 implicit none real R,deltaO if (iso_O18.gt.0) then deltaO=(R/tnat(iso_O18)-1)*1000.0 else write(*,*) 'iso_verif_egalite>deltaO18 260: iso_O18.gt.0=', & & iso_O18.gt.0 endif return end function deltaO ! ********** function dexcess(RD,RO) USE isotopes_mod, ONLY: tnat,iso_O18,iso_HDO implicit none real RD,RO,deltaD,deltaO,dexcess if ((iso_O18.gt.0).and.(iso_HDO.gt.0)) then deltaD=(RD/tnat(iso_HDO)-1)*1000.0 deltaO=(RO/tnat(iso_O18)-1)*1000.0 dexcess=deltaD-8*deltaO else write(*,*) 'iso_verif_egalite 1109: iso_O18,iso_HDO=',iso_O18,iso_HDO endif return end function dexcess ! ********** function delta_all(R,ixt) USE isotopes_mod, ONLY: tnat implicit none real R,delta_all integer ixt delta_all=(R/tnat(ixt)-1)*1000.0 return end function delta_all ! ********** function delta_to_R(delta,ixt) USE isotopes_mod, ONLY: tnat implicit none real delta,delta_to_R integer ixt delta_to_R=(delta/1000.0+1.0)*tnat(ixt) return end function delta_to_R ! ********** function o17excess(R17,R18) USE isotopes_mod, ONLY: tnat,iso_O18,iso_O17 implicit none real R17,R18,o17excess if ((iso_O17.gt.0).and.(iso_O18.gt.0)) then o17excess=1e6*(log(R17/tnat(iso_o17)) & & -0.528*log(R18/tnat(iso_O18))) ! write(*,*) 'o17excess=',o17excess else write(*,*) 'iso_verif_egalite>deltaD 260: iso_O17.gt.0,18=', & & iso_O17.gt.0,iso_O18.gt.0 endif return end function o17excess ! **************** subroutine iso_verif_egalite_vect2D( & & xt,q,err_msg,ni,n,m) USE isotopes_mod, ONLY: iso_eau implicit none ! inputs integer n,m,ni real q(n,m) real xt(ni,n,m) character*(*) err_msg ! locals integer iso_verif_egalite_nostop_loc integer i,j,ixt integer ifaux,jfaux !write(*,*) 'iso_verif_egalite_vect2D 1099 tmp: q(2,1),xt(iso_eau,2,1)=',q(2,1),xt(iso_eau,2,1) !write(*,*) 'ni,n,m=',ni,n,m,errmax,errmaxrel if (iso_eau.gt.0) then iso_verif_egalite_nostop_loc=0 do i=1,n do j=1,m if (abs(q(i,j)-xt(iso_eau,i,j)).gt.errmax) then if (abs((q(i,j)-xt(iso_eau,i,j))/ & & max(max(abs(q(i,j)),abs(xt(iso_eau,i,j))),1e-18)) & & .gt.errmaxrel) then iso_verif_egalite_nostop_loc=1 ifaux=i jfaux=j endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_egalite_nostop_loc.eq.1) then write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'xt,q=',xt(iso_eau,ifaux,jfaux),q(ifaux,jfaux) stop endif endif #ifdef ISOVERIF call iso_verif_noNaN_vect2D(xt,err_msg,ni,n,m) #endif return end subroutine iso_verif_egalite_vect2D subroutine iso_verif_egalite_vect1D( & & xt,q,err_msg,ni,n) USE isotopes_mod, ONLY: iso_eau implicit none ! inputs integer n,ni real q(n) real xt(ni,n) character*(*) err_msg ! locals integer iso_verif_egalite_nostop_loc integer i integer ifaux if (iso_eau.gt.0) then iso_verif_egalite_nostop_loc=0 do i=1,n if (abs(q(i)-xt(iso_eau,i)).gt.errmax) then if (abs((q(i)-xt(iso_eau,i))/ & & max(max(abs(q(i)),abs(xt(iso_eau,i))),1e-18)) & & .gt.errmaxrel) then iso_verif_egalite_nostop_loc=1 ifaux=i endif !if (abs((q(i)-xt(iso_eau,i))/ endif !if (abs(q(i)-xt(iso_eau,i)).gt.errmax) then enddo !do i=1,n if (iso_verif_egalite_nostop_loc.eq.1) then write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:' write(*,*) err_msg write(*,*) 'i=',ifaux write(*,*) 'xt,q=',xt(iso_eau,ifaux),q(ifaux) stop endif !if (iso_verif_egalite_nostop.eq.1) then endif !if (iso_eau.gt.0) then end subroutine iso_verif_egalite_vect1D subroutine iso_verif_egalite_std_vect( & & a,b,err_msg,n,m,errmax,errmaxrel) implicit none ! inputs integer n,m,ni real a(n,m) real b(n,m) character*(*) err_msg real errmax,errmaxrel ! locals integer iso_verif_egalite_nostop_loc integer i,j integer ifaux,jfaux iso_verif_egalite_nostop_loc=0 do i=1,n do j=1,m if (abs(a(i,j)-b(i,j)).gt.errmax) then if (abs((a(i,j)-b(i,j))/ & & max(max(abs(a(i,j)),abs(b(i,j))),1e-18)) & & .gt.errmaxrel) then iso_verif_egalite_nostop_loc=1 ifaux=i jfaux=j endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_egalite_nostop_loc.eq.1) then write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'a,b=',a(ifaux,jfaux),b(ifaux,jfaux) stop endif return end subroutine iso_verif_egalite_std_vect subroutine iso_verif_aberrant_vect2D( & & xt,q,err_msg,ni,n,m) use isotopes_mod, ONLY: ridicule,tnat,iso_HDO implicit none ! inputs integer n,m,ni real q(n,m) real xt(ni,n,m) character*(*) err_msg ! locals integer iso_verif_aberrant_nostop_loc integer i,j integer ifaux,jfaux !real deltaD if (iso_HDO.gt.0) then iso_verif_aberrant_nostop_loc=0 do i=1,n do j=1,m if (q(i,j).gt.ridicule) then if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .gt.deltalim).or. & & ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .lt.-borne)) then iso_verif_aberrant_nostop_loc=1 ifaux=i jfaux=j endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_aberrant_nostop_loc.eq.1) then write(*,*) 'erreur detectee par iso_verif_aberrant_vect2D:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) & & /q(ifaux,jfaux)) write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux) stop endif endif !if (iso_HDO.gt.0) then end subroutine iso_verif_aberrant_vect2D subroutine iso_verif_aberrant_enc_vect2D( & & xt,q,err_msg,ni,n,m) use isotopes_mod, ONLY: ridicule,tnat,iso_HDO implicit none ! inputs integer n,m,ni real q(n,m) real xt(ni,n,m) character*(*) err_msg ! locals integer iso_verif_aberrant_nostop_loc integer i,j integer ifaux,jfaux !real deltaD if (iso_HDO.gt.0) then iso_verif_aberrant_nostop_loc=0 do i=1,n do j=1,m if (q(i,j).gt.ridicule) then if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .gt.deltalim).or. & & ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .lt.deltaDmin).or. & & (xt(iso_HDO,i,j).lt.-borne).or. & & (xt(iso_HDO,i,j).gt.borne)) then iso_verif_aberrant_nostop_loc=1 ifaux=i jfaux=j endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_aberrant_nostop_loc.eq.1) then write(*,*) 'erreur detectee par ', & & 'iso_verif_aberrant_enc_vect2D:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) & & /q(ifaux,jfaux)) write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux) write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux) call abort_physic('isotopes_verif_mod','iso_verif_aberrant_enc_vect2D',1) endif endif !if (iso_HDO.gt.0) then end subroutine iso_verif_aberrant_enc_vect2D subroutine iso_verif_aberrant_enc_vect2D_ns( & & xt,q,err_msg,ni,n,m) use isotopes_mod, ONLY: ridicule,tnat,iso_HDO implicit none ! inputs integer n,m,ni real q(n,m) real xt(ni,n,m) character*(*) err_msg ! locals integer iso_verif_aberrant_nostop_loc integer i,j integer ifaux,jfaux !real deltaD if (iso_HDO.gt.0) then iso_verif_aberrant_nostop_loc=0 do i=1,n do j=1,m if (q(i,j).gt.ridicule) then if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .gt.deltalim).or. & & ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .lt.deltaDmin)) then iso_verif_aberrant_nostop_loc=1 ifaux=i jfaux=j endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_aberrant_nostop_loc.eq.1) then write(*,*) 'erreur detectee par ', & & 'iso_verif_aberrant_vect2D_ns:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) & & /q(ifaux,jfaux)) write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux) ! stop endif endif !if (iso_HDO.gt.0) then end subroutine iso_verif_aberrant_enc_vect2D_ns subroutine iso_verif_aberrant_vect2Dch( & & xt,q,err_msg,ni,n,m,deltaDmax) use isotopes_mod, ONLY: ridicule,tnat,iso_HDO implicit none ! inputs integer n,m,ni real q(n,m) real xt(ni,n,m) character*(*) err_msg real deltaDmax ! locals integer iso_verif_aberrant_nostop_loc integer i,j integer ifaux,jfaux !real deltaD if (iso_HDO.gt.0) then iso_verif_aberrant_nostop_loc=0 do i=1,n do j=1,m if (q(i,j).gt.ridicule) then if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .gt.deltaDmax).or. & & ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .lt.-borne)) then iso_verif_aberrant_nostop_loc=1 ifaux=i jfaux=j endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_aberrant_nostop_loc.eq.1) then write(*,*) 'erreur detectee par iso_verif_aberrant_vect2D:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) & & /q(ifaux,jfaux)) write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux) stop endif endif !if (iso_HDO.gt.0) then end subroutine iso_verif_aberrant_vect2Dch subroutine iso_verif_O18_aberrant_enc_vect2D( & & xt,q,err_msg,ni,n,m) use isotopes_mod, ONLY: ridicule,tnat,iso_HDO,iso_O18 implicit none ! inputs integer n,m,ni real q(n,m) real xt(ni,n,m) character*(*) err_msg ! locals integer iso_verif_aberrant_nostop_loc integer i,j integer ifaux,jfaux real deltaDloc,deltaoloc,dexcessloc if ((iso_HDO.gt.0).and.(iso_O18.gt.0)) then iso_verif_aberrant_nostop_loc=0 do i=1,n do j=1,m if (q(i,j).gt.ridicule) then deltaDloc=(xt(iso_HDO,i,j)/q(i,j)/tnat(iso_hdo)-1)*1000 deltaoloc=(xt(iso_O18,i,j)/q(i,j)/tnat(iso_O18)-1)*1000 dexcessloc=deltaDloc-8*deltaoloc if ((deltaDloc.lt.deltaDmin).or.(deltaoloc.lt.deltaDmin/2.0).or. & & (deltaDloc.gt.deltalim).or.(deltaoloc.gt.deltalim/8.0).or. & & ((deltaDloc.gt.-500.0).and.((dexcessloc.lt.dexcess_min) & & .or.(dexcessloc.gt.dexcess_max)))) then iso_verif_aberrant_nostop_loc=1 ifaux=i jfaux=j write(*,*) 'deltaD,deltao,dexcess=',deltaDloc,deltaoloc,dexcessloc endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_aberrant_nostop_loc.eq.1) then write(*,*) 'erreur detectee par ', & & 'iso_verif_aberrant_enc_vect2D:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux) write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux) call abort_physic('isotopes_verif_mod','iso_verif_aberrant_enc_vect2D',1) endif endif !if (iso_HDO.gt.0) then end subroutine iso_verif_O18_aberrant_enc_vect2D subroutine select_dim23_from4D(n1,n2,n3,n4, & & var,vec,i1,i4) implicit none ! inputs integer n1,n2,n3,n4 real var(n1,n2,n3,n4) integer i1,i4 ! outputs real vec(n2,n3) ! locals integer i2,i3 do i2=1,n2 do i3=1,n3 vec(i2,i3)=var(i1,i2,i3,i4) enddo enddo return end subroutine select_dim23_from4D subroutine select_dim4_from4D(ntime,nlev,nlat,nlon, & & var,vec,itime,ilev,ilat) implicit none ! inputs integer ntime,nlev,nlat,nlon real var(ntime,nlev,nlat,nlon) integer itime,ilev,ilat ! outputs real vec(nlon) ! locals integer ilon do ilon=1,nlon vec(ilon)=var(itime,ilev,ilat,ilon) enddo return end subroutine select_dim4_from4D subroutine select_dim5_from5D(n1,n2,n3,n4,n5, & & var,vec,i1,i2,i3,i4) implicit none ! inputs integer n1,n2,n3,n4,n5 real var(n1,n2,n3,n4,n5) integer i1,i2,i3,i4 ! outputs real vec(n5) ! locals integer i5 do i5=1,n5 vec(i5)=var(i1,i2,i3,i4,i5) enddo end subroutine select_dim5_from5D subroutine select_dim3_from3D(ntime,nlat,nlon, & & var,vec,itime,ilat) implicit none ! inputs integer ntime,nlat,nlon real var(ntime,nlat,nlon) integer itime,ilat ! outputs real vec(nlon) ! locals integer ilon do ilon=1,nlon vec(ilon)=var(itime,ilat,ilon) enddo end subroutine select_dim3_from3D subroutine select_dim23_from3D(n1,n2,n3, & & var,vec,i1) implicit none ! inputs integer n1,n2,n3 real var(n1,n2,n3) integer i1 ! outputs real vec(n2,n3) ! locals integer i2,i3 do i2=1,n2 do i3=1,n3 vec(i2,i3)=var(i1,i2,i3) enddo enddo end subroutine select_dim23_from3D subroutine putinto_dim23_from4D(n1,n2,n3,n4, & & var,vec,i1,i4) implicit none ! inputs integer n1,n2,n3,n4 real vec(n2,n3) integer i1,i4 ! inout real var(n1,n2,n3,n4) ! locals integer i2,i3 do i2=1,n2 do i3=1,n3 var(i1,i2,i3,i4)=vec(i2,i3) enddo enddo end subroutine putinto_dim23_from4D subroutine putinto_dim12_from4D(n1,n2,n3,n4, & & var,vec,i3,i4) implicit none ! inputs integer n1,n2,n3,n4 real vec(n1,n2) integer i3,i4 ! inout real var(n1,n2,n3,n4) ! locals integer i1,i2 do i1=1,n1 do i2=1,n2 var(i1,i2,i3,i4)=vec(i1,i2) enddo enddo end subroutine putinto_dim12_from4D subroutine putinto_dim23_from3D(n1,n2,n3, & & var,vec,i1) implicit none ! inputs integer n1,n2,n3 real vec(n2,n3) integer i1 ! inout real var(n1,n2,n3) ! locals integer i2,i3 do i2=1,n2 do i3=1,n3 var(i1,i2,i3)=vec(i2,i3) enddo enddo end subroutine putinto_dim23_from3D subroutine iso_verif_noNaN_par2D(x,err_msg,ni,n,m,ib,ie) implicit none ! si x est NaN, on affiche message ! d'erreur et return 1 si erreur ! input: integer n,m,ni,ib,ie real x(ni,n,m) character*(*) err_msg ! message d''erreur à afficher ! output ! locals integer i,j,ixt do i=ib,ie do j=1,m do ixt=1,ni if ((x(ixt,i,j).gt.-borne).and. & & (x(ixt,i,j).lt.borne)) then else !if ((x(ixt,i,j).gt.-borne).and. write(*,*) 'erreur detectee par iso_verif_nonNAN:' write(*,*) err_msg write(*,*) 'x,ixt,i,j=',x(ixt,i,j),ixt,i,j stop endif !if ((x(ixt,i,j).gt.-borne).and. enddo !do ixt=1,ni enddo !do j=1,m enddo !do i=1,n #ifdef ISOVERIF #else write(*,*) 'err_msg iso1772=',err_msg,': pourquoi verif?' stop #endif return end subroutine iso_verif_noNaN_par2D subroutine iso_verif_aberrant_enc_par2D( & & xt,q,err_msg,ni,n,m,ib,ie) use isotopes_mod, ONLY: ridicule,tnat,iso_HDO implicit none ! inputs integer n,m,ni,ib,ie real q(n,m) real xt(ni,n,m) character*(*) err_msg ! locals integer iso_verif_aberrant_nostop_loc integer i,j integer ifaux,jfaux !real deltaD if (iso_HDO.gt.0) then iso_verif_aberrant_nostop_loc=0 do i=ib,ie do j=1,m if (q(i,j).gt.ridicule) then if (((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .gt.deltalim).or. & & ((xt(iso_HDO,i,j)/q(i,j)/tnat(iso_HDO)-1)*1000.0 & & .lt.deltaDmin)) then iso_verif_aberrant_nostop_loc=1 ifaux=i jfaux=j endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_aberrant_nostop_loc.eq.1) then write(*,*) 'erreur detectee par iso_verif_aberrant_par2D:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'deltaD=',deltaD(xt(iso_HDO,ifaux,jfaux) & & /q(ifaux,jfaux)) write(*,*) 'xt(:,ifaux,jfaux)=',xt(:,ifaux,jfaux) write(*,*) 'q(ifaux,jfaux)=',q(ifaux,jfaux) stop endif endif !if (iso_HDO.gt.0) then end subroutine iso_verif_aberrant_enc_par2D subroutine iso_verif_egalite_par2D( & & xt,q,err_msg,ni,n,m,ib,ie) USE isotopes_mod, ONLY: iso_eau implicit none ! inputs integer n,m,ni,ib,ie real q(n,m) real xt(ni,n,m) character*(*) err_msg ! locals integer iso_verif_egalite_nostop_loc integer i,j integer ifaux,jfaux if (iso_eau.gt.0) then iso_verif_egalite_nostop_loc=0 do i=ib,ie do j=1,m if (abs(q(i,j)-xt(iso_eau,i,j)).gt.errmax) then if (abs((q(i,j)-xt(iso_eau,i,j))/ & & max(max(abs(q(i,j)),abs(xt(iso_eau,i,j))),1e-18)) & & .gt.errmaxrel) then iso_verif_egalite_nostop_loc=1 ifaux=i jfaux=j endif endif enddo !do j=1,m enddo !do i=1,n if (iso_verif_egalite_nostop_loc.eq.1) then write(*,*) 'erreur detectee par iso_verif_egalite_vect2D:' write(*,*) err_msg write(*,*) 'i,j=',ifaux,jfaux write(*,*) 'xt,q=',xt(iso_eau,ifaux,jfaux),q(ifaux,jfaux) stop endif endif end subroutine iso_verif_egalite_par2D #ifdef ISOTRAC function iso_verif_traceur_choix_nostop(x,err_msg, & & errmax,errmaxrel,ridicule_trac,deltalimtrac) use isotopes_mod, ONLY: iso_HDO implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher real errmax,errmaxrel,ridicule_trac,deltalimtrac ! output integer iso_verif_traceur_choix_nostop ! locals !integer iso_verif_traceur_noNaN_nostop !integer iso_verif_tracm_choix_nostop !integer iso_verif_tracdD_choix_nostop !integer iso_verif_tracpos_choix_nostop iso_verif_traceur_choix_nostop=0 ! verif noNaN if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then iso_verif_traceur_choix_nostop=1 endif ! verif masse if (iso_verif_tracm_choix_nostop(x,err_msg, & & errmax,errmaxrel).eq.1) then iso_verif_traceur_choix_nostop=1 endif ! verif deltaD if (iso_HDO.gt.0) then if (iso_verif_tracdD_choix_nostop(x,err_msg, & & ridicule_trac,deltalimtrac).eq.1) then iso_verif_traceur_choix_nostop=1 endif endif !if (iso_HDO.gt.0) then ! verif pas aberramment negatif if (iso_verif_tracpos_choix_nostop(x,err_msg, & & 1e-5).eq.1) then iso_verif_traceur_choix_nostop=1 endif end function iso_verif_traceur_choix_nostop function iso_verif_tracnps_choix_nostop(x,err_msg, & & errmax,errmaxrel,ridicule_trac,deltalimtrac) USE isotopes_mod, ONLY: iso_HDO implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on ne vérfie pas la positivité ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher real errmax,errmaxrel,ridicule_trac,deltalimtrac ! output integer iso_verif_tracnps_choix_nostop ! locals !integer iso_verif_traceur_noNaN_nostop !integer iso_verif_tracm_choix_nostop !integer iso_verif_tracdD_choix_nostop iso_verif_tracnps_choix_nostop=0 ! verif noNaN if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then iso_verif_tracnps_choix_nostop=1 endif ! verif masse if (iso_verif_tracm_choix_nostop(x,err_msg, & & errmax,errmaxrel).eq.1) then iso_verif_tracnps_choix_nostop=1 endif ! verif deltaD if (iso_HDO.gt.0) then if (iso_verif_tracdD_choix_nostop(x,err_msg, & & ridicule_trac,deltalimtrac).eq.1) then iso_verif_tracnps_choix_nostop=1 endif endif ! if (iso_HDO.gt.0) then return end function iso_verif_tracnps_choix_nostop function iso_verif_tracpos_choix_nostop(x,err_msg,seuil) use isotopes_mod, only: isoName implicit none ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher real seuil ! output integer iso_verif_tracpos_choix_nostop ! locals integer lnblnk integer iiso,ixt !integer iso_verif_positif_choix_nostop iso_verif_tracpos_choix_nostop=0 do ixt=niso+1,ntraciso if (iso_verif_positif_choix_nostop(x(ixt),seuil,err_msg// & & ', verif positif, iso'//TRIM(isoName(ixt))).eq.1) then iso_verif_tracpos_choix_nostop=1 endif enddo end function iso_verif_tracpos_choix_nostop function iso_verif_traceur_noNaN_nostop(x,err_msg) use isotopes_mod, only: isoName implicit none ! on vérifie juste que pas NaN ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_traceur_noNaN_nostop ! locals integer lnblnk integer iiso,ixt !integer iso_verif_nonaN_nostop iso_verif_traceur_noNaN_nostop=0 do ixt=niso+1,ntraciso ! write(*,*) 'iso_verif_traceurs 154: iiso,ixt=',iiso,ixt if (iso_verif_noNaN_nostop(x(ixt),err_msg// & & ', verif trac no NaN, iso'//TRIM(isoName(ixt))) & & .eq.1) then iso_verif_traceur_noNaN_nostop=1 endif enddo end function iso_verif_traceur_noNaN_nostop function iso_verif_tracm_choix_nostop(x,err_msg, & & errmaxin,errmaxrelin) use isotopes_mod, ONLY: ridicule,isoName ! on vérifie juste bilan de masse implicit none ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher real errmaxin,errmaxrelin ! output integer iso_verif_tracm_choix_nostop ! locals !integer iso_verif_egalite_choix_nostop integer iiso,izone,ixt real xtractot iso_verif_tracm_choix_nostop=0 do iiso=1,niso xtractot=0.0 do izone=1,nzone ixt=itZonIso(izone,iiso) xtractot=xtractot+x(ixt) enddo if (iso_verif_egalite_choix_nostop(xtractot,x(iiso), & & err_msg//', verif trac egalite1, iso '// & & TRIM(isoName(iiso)), & & errmaxin,errmaxrelin).eq.1) then write(*,*) 'iso_verif_traceur 202: x=',x ! write(*,*) 'xtractot=',xtractot do izone=1,nzone ixt=itZonIso(izone,iiso) write(*,*) 'izone,iiso,ixt=',izone,iiso,ixt enddo iso_verif_tracm_choix_nostop=1 endif ! verif ajoutée le 19 fev 2011 if ((abs(xtractot).lt.ridicule**2).and. & & (abs(x(iiso)).gt.ridicule)) then write(*,*) err_msg,', verif masse traceurs, iso ', & & TRIM(isoName(iiso)) write(*,*) 'iso_verif_traceur 209: x=',x ! iso_verif_tracm_choix_nostop=1 endif enddo !do iiso=1,ntraceurs_iso end function iso_verif_tracm_choix_nostop function iso_verif_tracdD_choix_nostop(x,err_msg, & & ridicule_trac,deltalimtrac) USE isotopes_mod, ONLY: iso_eau, iso_HDO use isotrac_mod, only: strtrac ! on vérifie juste deltaD implicit none ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher real ridicule_trac,deltalimtrac ! output integer iso_verif_tracdD_choix_nostop ! locals integer izone,ieau,ixt !integer iso_verif_aberrant_choix_nostop iso_verif_tracdD_choix_nostop=0 if ((iso_eau.gt.0).and.(iso_HDO.gt.0)) then do izone=1,nzone ieau=itZonIso(izone,iso_eau) ixt=itZonIso(izone,iso_HDO) if (iso_verif_aberrant_choix_nostop(x(ixt),x(ieau), & & ridicule_trac,deltalimtrac,err_msg// & & ', verif trac no aberrant zone '//strtrac(izone)) & & .eq.1) then iso_verif_tracdD_choix_nostop=1 endif ! if (x(ieau).gt.ridicule) then ! call iso_verif_aberrant(x(ixt)/x(ieau), ! : err_msg//', verif trac no aberrant zone ' ! : //strtrac(izone)) ! endif enddo !do izone=1,nzone endif ! if ((iso_eau.gt.0).and.(iso_HDO.gt.0)) then end function iso_verif_tracdD_choix_nostop INTEGER FUNCTION iso_verif_tag17_q_deltaD_chns(x,err_msg) RESULT(res) USE isotopes_mod, ONLY: iso_HDO, iso_eau, ridicule USE isotrac_mod, ONLY: nzone_temp, option_traceurs IMPLICIT NONE REAL, INTENT(IN) :: x(ntraciso) CHARACTER(LEN=*), INTENT(IN) :: err_msg INTEGER :: ieau, ixt, ieau1 res = 0 IF(ALL([17,18]/=option_traceurs)) RETURN !--- Check whether * deltaD(highest tagging layer) < 200 permil ! * q < ieau=itZonIso(nzone_temp,iso_eau) ixt=itZonIso(nzone_temp,iso_HDO) IF(x(ieau)>ridicule) THEN IF(iso_verif_positif_nostop(-200.0-deltaD(x(ixt)/x(ieau)), err_msg//': deltaDt05 trop fort')==1) THEN res=1; write(*,*) 'x=',x END IF END IF IF(iso_verif_positif_nostop(2.0e-3-x(ieau),err_msg//': qt05 trop fort')==1) THEN res=1; write(*,*) 'x=',x END IF !--- Check whether q is small ; then, qt01 < 10% IF(x(iso_eau)<2.0e-3) THEN ieau1= itZonIso(1,iso_eau) IF(iso_verif_positif_nostop(0.1-(x(ieau1)/x(iso_eau)),err_msg//': qt01 trop abondant')==1) THEN res=1; write(*,*) 'x=',x END IF END IF END FUNCTION iso_verif_tag17_q_deltaD_chns SUBROUTINE iso_verif_trac17_q_deltaD(x,err_msg) USE isotrac_mod, ONLY: nzone_temp, option_traceurs IMPLICIT NONE REAL, INTENT(IN) :: x(ntraciso) CHARACTER(LEN=*), INTENT(IN) :: err_msg IF(ALL([17,18]/=option_traceurs)) RETURN IF(nzone_temp>=5) THEN IF(iso_verif_tag17_q_deltaD_chns(x,err_msg)==1) STOP END IF END SUBROUTINE iso_verif_trac17_q_deltaD subroutine iso_verif_traceur(x,err_msg) use isotrac_mod, only: ridicule_trac implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher ! locals !integer iso_verif_traceur_choix_nostop if (iso_verif_traceur_choix_nostop(x,err_msg, & & errmax,errmaxrel,ridicule_trac,deltalimtrac) & & .eq.1) then stop endif end subroutine iso_verif_traceur subroutine iso_verif_traceur_retourne3D(x,n1,n2,n3, & & i1,i2,i3,err_msg) use isotrac_mod, only: ridicule_trac implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs integer n1,n2,n3 real x(n1,n2,n3,ntraciso) character*(*) err_msg ! message d''erreur à afficher integer i1,i2,i3 ! locals !integer iso_verif_traceur_choix_nostop real xiso(ntraciso) call select_dim4_from4D(n1,n2,n3,ntraciso, & & x,xiso,i1,i2,i3) if (iso_verif_traceur_choix_nostop(xiso,err_msg, & & errmax,errmaxrel,ridicule_trac,deltalimtrac) & & .eq.1) then stop endif end subroutine iso_verif_traceur_retourne3D subroutine iso_verif_traceur_retourne4D(x,n1,n2,n3,n4, & & i1,i2,i3,i4,err_msg) use isotrac_mod, only: ridicule_trac implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs integer n1,n2,n3,n4 real x(n1,n2,n3,n4,ntraciso) character*(*) err_msg ! message d''erreur à afficher integer i1,i2,i3,i4 ! locals !integer iso_verif_traceur_choix_nostop real xiso(ntraciso) call select_dim5_from5D(n1,n2,n3,n4,ntraciso, & & x,xiso,i1,i2,i3,i4) if (iso_verif_traceur_choix_nostop(xiso,err_msg, & & errmax,errmaxrel,ridicule_trac,deltalimtrac) & & .eq.1) then stop endif end subroutine iso_verif_traceur_retourne4D subroutine iso_verif_traceur_retourne2D(x,n1,n2, & & i1,i2,err_msg) use isotrac_mod, only: ridicule_trac implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs integer n1,n2 real x(n1,n2,ntraciso) character*(*) err_msg ! message d''erreur à afficher integer i1,i2 ! locals !integer iso_verif_traceur_choix_nostop real xiso(ntraciso) call select_dim3_from3D(n1,n2,ntraciso, & & x,xiso,i1,i2) if (iso_verif_traceur_choix_nostop(xiso,err_msg, & & errmax,errmaxrel,ridicule_trac,deltalimtrac) & & .eq.1) then stop endif end subroutine iso_verif_traceur_retourne2D subroutine iso_verif_traceur_vect(x,n,m,err_msg) USE isotopes_mod, ONLY: iso_HDO implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs integer n,m real x(ntraciso,n,m) character*(*) err_msg ! message d''erreur à afficher ! locals logical iso_verif_traceur_nostop integer i,j,ixt,iiso,izone,ieau integer ifaux,jfaux,ixtfaux call iso_verif_traceur_noNaN_vect(x,n,m,err_msg) ! verif masse: iso_verif_tracm_choix_nostop call iso_verif_trac_masse_vect(x,n,m,err_msg,errmax,errmaxrel) ! verif deltaD: iso_verif_tracdD_choix_nostop if (iso_HDO.gt.0) then call iso_verif_tracdd_vect(x,n,m,err_msg) endif !if (iso_HDO.gt.0) then ! verif pas aberramment negatif: iso_verif_tracpos_choix_nostop call iso_verif_tracpos_vect(x,n,m,err_msg,1e-5) end subroutine iso_verif_traceur_vect subroutine iso_verif_tracnps_vect(x,n,m,err_msg) USE isotopes_mod, ONLY: iso_HDO implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs integer n,m real x(ntraciso,n,m) character*(*) err_msg ! message d''erreur à afficher ! locals logical iso_verif_traceur_nostop integer i,j,ixt,iiso,izone,ieau integer ifaux,jfaux,ixtfaux call iso_verif_traceur_noNaN_vect(x,n,m,err_msg) ! verif masse: iso_verif_tracm_choix_nostop call iso_verif_trac_masse_vect(x,n,m,err_msg,errmax,errmaxrel) ! verif deltaD: iso_verif_tracdD_choix_nostop if (iso_HDO.gt.0) then call iso_verif_tracdd_vect(x,n,m,err_msg) endif !if (iso_HDO.gt.0) then end subroutine iso_verif_tracnps_vect subroutine iso_verif_traceur_noNaN_vect(x,n,m,err_msg) implicit none ! inputs integer n,m real x(ntraciso,n,m) character*(*) err_msg ! message d''erreur à afficher ! locals logical iso_verif_traceur_nostop integer i,j,ixt,iiso integer ifaux,jfaux,ixtfaux iso_verif_traceur_nostop=.false. ! verif noNaN: iso_verif_traceur_noNaN_nostop do j=1,m do i=1,n do ixt=niso+1,ntraciso if ((x(ixt,i,j).gt.-borne).and.(x(ixt,i,j).lt.borne)) then else !if ((x.gt.-borne).and.(x.lt.borne)) then iso_verif_traceur_nostop=.true. ifaux=i jfaux=i endif !if ((x.gt.-borne).and.(x.lt.borne)) then enddo !do ixt=niso+1,ntraciso enddo ! do i=1,n enddo ! do j=1,m if (iso_verif_traceur_nostop) then write(*,*) 'erreur detectée par iso_verif_nonNAN ', & & 'dans iso_verif_traceur_vect' write(*,*) '' write(*,*) err_msg write(*,*) 'x=',x(:,ifaux,jfaux) stop endif end subroutine iso_verif_traceur_noNaN_vect subroutine iso_verif_trac_masse_vect(x,n,m,err_msg, & & errmax,errmaxrel) use isotopes_mod, only: isoName implicit none ! inputs integer n,m real x(ntraciso,n,m) character*(*) err_msg ! message d''erreur à afficher real errmax,errmaxrel ! locals logical iso_verif_traceur_nostop integer i,j,ixt,iiso,izone integer ifaux,jfaux,ixtfaux real xtractot(n,m) real xiiso(n,m) do iiso=1,niso do j=1,m do i=1,n xtractot(i,j)=0.0 xiiso(i,j)=x(iiso,i,j) do izone=1,nzone ixt=itZonIso(izone,iiso) xtractot(i,j)=xtractot(i,j)+x(ixt,i,j) enddo !do izone=1,nzone enddo !do i=1,n enddo !do j=1,m call iso_verif_egalite_std_vect( & & xtractot,xiiso, & & err_msg//', verif trac egalite2, iso ' & & //TRIM(isoName(iiso)), & & n,m,errmax,errmaxrel) enddo !do iiso=1,niso end subroutine iso_verif_trac_masse_vect subroutine iso_verif_tracdd_vect(x,n,m,err_msg) use isotopes_mod, only: iso_HDO,iso_eau use isotrac_mod, only: strtrac implicit none ! inputs integer n,m real x(ntraciso,n,m) character*(*) err_msg ! message d''erreur à afficher ! locals integer i,j,iiso,izone,ieau,ixt real xiiso(niso,n,m) real xeau(n,m) integer lnblnk if (iso_HDO.gt.0) then do izone=1,nzone ieau=itZonIso(izone,iso_eau) do iiso=1,niso ixt=itZonIso(izone,iiso) do j=1,m do i=1,n xiiso(iiso,i,j)=x(ixt,i,j) enddo !do i=1,n enddo !do j=1,m enddo !do iiso=1,niso do j=1,m do i=1,n xeau(i,j)=x(ieau,i,j) enddo !do i=1,n enddo !do j=1,m call iso_verif_aberrant_vect2Dch( & & xiiso,xeau,err_msg//strtrac(izone),niso,n,m, & & deltalimtrac) enddo !do izone=1,nzone endif !if (iso_HDO.gt.0) then end subroutine iso_verif_tracdd_vect subroutine iso_verif_tracpos_vect(x,n,m,err_msg,seuil) implicit none ! inputs integer n,m real x(ntraciso,n,m) character*(*) err_msg ! message d''erreur à afficher real seuil ! locals integer i,j,ixt logical iso_verif_traceur_nostop integer ifaux,jfaux,ixtfaux iso_verif_traceur_nostop=.false. do j=1,m do i=1,n do ixt=niso+1,ntraciso if (x(ixt,i,j).lt.-seuil) then ifaux=i jfaux=j ixtfaux=ixt iso_verif_traceur_nostop=.true. endif enddo !do ixt=niso+1,ntraciso enddo !do i=1,n enddo !do j=1,m if (iso_verif_traceur_nostop) then write(*,*) 'erreur detectée par verif positif ', & & 'dans iso_verif_traceur_vect' write(*,*) '' write(*,*) err_msg write(*,*) 'x=',x(:,ifaux,jfaux) stop endif end subroutine iso_verif_tracpos_vect subroutine iso_verif_tracnps(x,err_msg) use isotrac_mod, only: ridicule_trac implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher ! locals !integer iso_verif_tracnps_choix_nostop if (iso_verif_tracnps_choix_nostop(x,err_msg, & & errmax,errmaxrel,ridicule_trac,deltalimtrac) & & .eq.1) then stop endif end subroutine iso_verif_tracnps subroutine iso_verif_tracpos_choix(x,err_msg,seuil) implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher real seuil ! locals !integer iso_verif_tracpos_choix_nostop if (iso_verif_tracpos_choix_nostop(x,err_msg,seuil) & & .eq.1) then stop endif end subroutine iso_verif_tracpos_choix subroutine iso_verif_traceur_choix(x,err_msg, & & errmax,errmaxrel,ridicule_trac_loc,deltalimtrac) implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher real errmax,errmaxrel,ridicule_trac_loc,deltalimtrac ! locals !integer iso_verif_traceur_choix_nostop if (iso_verif_traceur_choix_nostop(x,err_msg, & & errmax,errmaxrel,ridicule_trac_loc,deltalimtrac) & & .eq.1) then stop endif end subroutine iso_verif_traceur_choix function iso_verif_traceur_nostop(x,err_msg) use isotrac_mod, only: ridicule_trac !use isotopes_verif, only: errmax,errmaxrel,deltalimtrac implicit none ! vérifier des choses sur les traceurs ! * toutes les zones donne t l'istope total ! * pas de deltaD aberrant ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_traceur_nostop ! locals !integer iso_verif_traceur_choix_nostop iso_verif_traceur_nostop= & & iso_verif_traceur_choix_nostop(x,err_msg, & & errmax,errmaxrel,ridicule_trac,deltalimtrac) end function iso_verif_traceur_nostop subroutine iso_verif_traceur_justmass(x,err_msg) implicit none ! on vérifie que noNaN et masse ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher ! locals !integer iso_verif_traceur_noNaN_nostop !integer iso_verif_tracm_choix_nostop ! verif noNaN if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then stop endif ! verif masse if (iso_verif_tracm_choix_nostop(x,err_msg, & & errmax,errmaxrel).eq.1) then stop endif end subroutine iso_verif_traceur_justmass function iso_verif_traceur_jm_nostop(x,err_msg) implicit none ! on vérifie que noNaN et masse ! on prend les valeurs pas défaut pour ! errmax,errmaxrel,ridicule_trac,deltalimtrac ! inputs real x(ntraciso) character*(*) err_msg ! message d''erreur à afficher ! output integer iso_verif_traceur_jm_nostop ! locals ! integer iso_verif_traceur_noNaN_nostop !integer iso_verif_tracm_choix_nostop iso_verif_traceur_jm_nostop=0 ! ! verif noNaN ! if (iso_verif_traceur_noNaN_nostop(x,err_msg).eq.1) then ! iso_verif_traceur_jm_nostop=1 ! endif ! verif masse if (iso_verif_tracm_choix_nostop(x,err_msg, & & errmax,errmaxrel).eq.1) then iso_verif_traceur_jm_nostop=1 endif end function iso_verif_traceur_jm_nostop subroutine iso_verif_tag17_q_deltaD_vect(x,n,m,err_msg) USE isotopes_mod, ONLY: tnat,iso_eau, ridicule,iso_HDO use isotrac_mod, only: option_traceurs,nzone_temp implicit none ! inputs integer n,m real x(ntraciso,n,m) character*(*) err_msg ! locals !integer iso_verif_positif_nostop !real deltaD integer ieau,ixt,ieau1 integer i,k if ((option_traceurs.eq.17).or. & & (option_traceurs.eq.18)) then ! verifier que deltaD du tag de la couche la plus haute < ! 200 permil, et vérifier que son q est inférieur à ieau=itZonIso(nzone_temp,iso_eau) ixt=itZonIso(nzone_temp,iso_HDO) ieau1=itZonIso(1,iso_eau) do i=1,n do k=1,m if (x(ieau,i,k).gt.ridicule) then if ((x(ixt,i,k)/x(ieau,i,k)/tnat(iso_HDO)-1)*1000 & & .gt.-200.0) then write(*,*) err_msg,', vect:deltaDt05 trop fort' write(*,*) 'i,k=',i,k write(*,*) 'x(:,i,k)=',x(:,i,k) stop endif !if (iso_verif_positif_nostop(-200.0-deltaD(x(ixt),x(ieau)), endif !if (x(ieau).gt.ridicule) then if (x(ieau,i,k).gt.2.0e-3) then write(*,*) err_msg,', vect:qt05 trop fort' write(*,*) 'i,k=',i,k write(*,*) 'x(:,i,k)=',x(:,i,k) stop endif !if (iso_verif_positif_nostop(1.0e-3-x(ieau), if (x(iso_eau,i,k).lt.2.0e-3) then if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then write(*,*) err_msg,', vect: qt01 trop abondant' write(*,*) 'i,k=',i,k write(*,*) 'ieau1,iso_eau,x(ieau1,i,k),', & & 'x(iso_eau,i,k)=',ieau1,iso_eau, & & x(ieau1,i,k),x(iso_eau,i,k) write(*,*) 'x(:,i,k)=',x(:,i,k) stop endif !if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then endif enddo !do k=1,m enddo !do i=1,n endif !if (option_traceurs.eq.17) then end subroutine iso_verif_tag17_q_deltaD_vect subroutine iso_verif_tag17_q_deltaD_vect_ret3D(x,n,m,nq,err_msg) USE isotopes_mod, ONLY: tnat,iso_eau,iso_HDO,ridicule use isotrac_mod, only: option_traceurs,nzone_temp implicit none ! inputs integer n,m,nq real x(n,m,nq,ntraciso) character*(*) err_msg ! locals !integer iso_verif_positif_nostop !real deltaD integer ieau,ixt,ieau1 integer i,k,iq if ((option_traceurs.eq.17).or. & & (option_traceurs.eq.18)) then ! verifier que deltaD du tag de la couche la plus haute < ! 200 permil, et vérifier que son q est inférieur à ieau=itZonIso(nzone_temp,iso_eau) ixt=itZonIso(nzone_temp,iso_HDO) ieau1=itZonIso(1,iso_eau) do iq=1,nq do i=1,n do k=1,m if (x(i,k,iq,ieau).gt.ridicule) then if ((x(i,k,iq,ixt)/x(i,k,iq,ieau)/tnat(iso_HDO)-1)*1000 & & .gt.-200.0) then write(*,*) err_msg,', vect:deltaDt05 trop fort' write(*,*) 'i,k=',i,k write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:) stop endif !if (iso_verif_positif_nostop(-200.0-deltaD(x(ixt),x(ieau)), endif !if (x(ieau).gt.ridicule) then if (x(i,k,iq,ieau).gt.2.0e-3) then write(*,*) err_msg,', vect:qt05 trop fort' write(*,*) 'i,k=',i,k write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:) stop endif !if (iso_verif_positif_nostop(1.0e-3-x(ieau), if (x(i,k,iq,iso_eau).lt.2.0e-3) then if (x(i,k,iq,ieau1)/x(i,k,iq,iso_eau).gt.0.1) then write(*,*) err_msg,', vect: qt01 trop abondant' write(*,*) 'i,k=',i,k write(*,*) 'ieau1,iso_eau,x(i,k,iq,ieau1),', & & 'x(i,k,iq,ieau)=',ieau1,iso_eau, & & x(i,k,iq,ieau1),x(i,k,iq,iso_eau) write(*,*) 'x(i,k,iq,:)=',x(i,k,iq,:) stop endif !if (x(ieau1,i,k)/x(iso_eau,i,k).gt.0.1) then endif enddo !do k=1,m enddo !do i=1,n enddo ! do iq=1,nq endif !if (option_traceurs.eq.17) then end subroutine iso_verif_tag17_q_deltaD_vect_ret3D #endif ! endif ISOTRAC END MODULE isotopes_verif_mod #endif ! endif ISOVERIF