Ignore:
Timestamp:
Apr 5, 2022, 3:44:30 PM (3 years ago)
Author:
dcugnet
Message:
  • New water names: H2Ov, H2Ol, H2Oi, H2Or -> H2O_g, H2O_l, H2O_s, H2O_r.
  • Corrections for the lOldCode=.FALSE., not activated yet.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d_common/infotrac.F90

    r4082 r4120  
    33MODULE infotrac
    44
    5    USE       strings_mod, ONLY: msg, find, strIdx,  strFind, strParse, dispTable, int2str,  reduceExpr,  &
     5   USE       strings_mod, ONLY: msg, find, strIdx,  strFind, strParse, dispTable, int2str,  reduceExpr, &
    66                          cat, fmsg, test, strTail, strHead, strStack, strReduce, bool2str, maxlen, testFile
    7    USE readTracFiles_mod, ONLY: trac_type, readTracersFiles, addPhase,  phases_sep,  nphases, ancestor,  &
    8                                 isot_type, readIsotopesFile, delPhase,   old_phases, getKey_init, tran0, &
    9                                 keys_type, initIsotopes,  indexUpdate, known_phases, getKey, setGeneration, &
    10                                 new2oldPhase
    11 
     7   USE readTracFiles_mod, ONLY: trac_type, readTracersFiles, addPhase, indexUpdate,  nphases, ancestor,  &
     8                                isot_type, old2newName,      delPhase,               getKey_init, tran0, &
     9                                keys_type, initIsotopes,     getPhase, known_phases, getKey, setGeneration
    1210   IMPLICIT NONE
    1311
    … …  
    2321   PUBLIC :: isotopes,  nbIso                              !--- Derived type, full isotopes families database + nb of families
    2422   PUBLIC :: isoSelect, ixIso                              !--- Isotopes family selection tool + selected family index
    25    PUBLIC :: min_qParent, min_qMass, min_ratio             !--- Min. values for various isotopic quantities
    2623   !=== FOR ISOTOPES: Specific to water
    2724   PUBLIC :: iH2O, tnat, alpha_ideal                       !--- H2O isotopes index, natural abundance, fractionning coeff.
     25   PUBLIC :: min_qParent, min_qMass, min_ratio             !--- Min. values for various isotopic quantities
    2826   !=== FOR ISOTOPES: Depending on the selected isotopes family
    2927   PUBLIC :: isotope, isoKeys                              !--- Selected isotopes database + associated keys (cf. getKey)
    3028   PUBLIC :: isoName, isoZone, isoPhas                     !--- Isotopes and tagging zones names, phases
    3129   PUBLIC :: niso,    nzone,   nphas,   ntiso              !---  " " numbers + isotopes & tagging tracers number
    32    PUBLIC :: iZonIso, iTraPha                              !--- 2D index tables to get "iq" index
     30   PUBLIC :: itZonIso, index_trac                          !--- idx "it" (in "isoName(1:niso)") = function(tagging idx, isotope idx)
     31   PUBLIC :: iqTraPha, iqiso                               !--- idx "iq" (in "qx") = function(isotope idx, phase idx) + aliases
    3332   PUBLIC :: isoCheck                                      !--- Run isotopes checking routines
    3433   !=== FOR BOTH TRACERS AND ISOTOPES
    3534   PUBLIC :: getKey                                        !--- Get a key from "tracers" or "isotope"
    3635
    37    PUBLIC :: ntraciso, ntraceurs_zone, iqiso
    38    PUBLIC :: ok_isotopes, ok_iso_verif, ok_isotrac, ok_init_iso, use_iso
    39    PUBLIC :: index_trac, iso_indnum, indnum_fn_num, niso_possibles
     36   !=== OLD QUANTITIES OR ALIASES FOR OLDER NAMES (TO BE REMOVED SOON)
     37   PUBLIC :: ntraciso, ntraceurs_zone
     38   PUBLIC :: ok_isotopes, ok_iso_verif, ok_isotrac, use_iso
     39   PUBLIC :: iso_num, iso_indnum, indnum_fn_num, niso_possibles
    4040   PUBLIC :: qperemin, masseqmin, ratiomin
    4141
    … …  
    9393!=== DERIVED TYPE EMBEDDING MOST OF THE ISOTOPES-RELATED QUANTITIES (LENGTH: nbIso, NUMBER OF ISOTOPES FAMILIES)
    9494!    Each entry is accessible using "%" sign.
    95 !  |-----------------+--------------------------------------------------+----------------+-----------------+
    96 !  |  entry | length | Meaning                                          | Former name    | Possible values |
    97 !  |-----------------+--------------------------------------------------+----------------+-----------------+
    98 !  | parent          | Parent tracer (isotopes family name)             |                |                 |
    99 !  | keys   | niso   | Isotopes keys/values pairs list + number         |                |                 |
    100 !  | trac   | ntiso  | Isotopes + tagging tracers list + number         |                |                 |
    101 !  | zone   | nzone  | Geographic tagging zones   list + number         |                |                 |
    102 !  | phase  | nphas  | Phases                     list + number         |                | [g][l][s], 1:3  |
    103 !  | niso            | Number of isotopes, excluding tagging tracers    |                |                 |
    104 !  | ntiso           | Number of isotopes, including tagging tracers    | ntraciso       |                 |
    105 !  | nzone           | Number of geographic tagging zones               | ntraceurs_zone |                 |
    106 !  | nphas           | Number of phases                                 |                |                 |
    107 !  | iTraPha         | Index in "trac(1:niso)" = f(name(1:ntiso)),phas) | iqiso          | 1:niso          |
    108 !  | iZonIso         | Index in "trac(1:ntiso)" = f(zone, name(1:niso)) | index_trac     | 1:nzone         |
    109 !  |-----------------+--------------------------------------------------+----------------+-----------------+
     95!  |-----------------+--------------------------------------------------+--------------------+-----------------+
     96!  |  entry | length | Meaning                                          |    Former name     | Possible values |
     97!  |-----------------+--------------------------------------------------+--------------------+-----------------+
     98!  | parent          | Parent tracer (isotopes family name)             |                    |                 |
     99!  | keys   | niso   | Isotopes keys/values pairs list + number         |                    |                 |
     100!  | trac   | ntiso  | Isotopes + tagging tracers list + number         | / | ntraciso       |                 |
     101!  | zone   | nzone  | Geographic tagging zones   list + number         | / | ntraceurs_zone |                 |
     102!  | phase  | nphas  | Phases                     list + number         |                    | [g][l][s], 1:3  |
     103!  | iqTraPha        | Index in "qx"           = f(name(1:ntiso)),phas) | iqiso              | 1:nqtot         |
     104!  | itZonIso        | Index in "trac(1:ntiso)"= f(zone, name(1:niso))  | index_trac         | 1:ntiso         |
     105!  +-----------------+--------------------------------------------------+--------------------+-----------------+
    110106
    111107   REAL, PARAMETER :: min_qParent = 1.e-30, min_qMass = 1.e-18, min_ratio = 1.e-16 ! MVals et CRisi
    … …  
    127123   TYPE(isot_type),         SAVE, POINTER :: isotope            !--- CURRENTLY SELECTED ISOTOPES FAMILY DESCRIPTOR
    128124   INTEGER,                 SAVE          :: ixIso, iH2O        !--- Index of the selected isotopes family and H2O family
    129    LOGICAL,                 SAVE          :: isoCheck           !--- Flag to trigger the checking routines
     125   LOGICAL,                 SAVE, POINTER :: isoCheck           !--- Flag to trigger the checking routines
    130126   TYPE(keys_type),         SAVE, POINTER :: isoKeys(:)         !--- ONE SET OF KEYS FOR EACH ISOTOPE (LISTED IN isoName)
    131127   CHARACTER(LEN=maxlen),   SAVE, POINTER :: isoName(:),   &    !--- ISOTOPES NAMES FOR THE CURRENTLY SELECTED FAMILY
    132128                                             isoZone(:),   &    !--- TAGGING ZONES  FOR THE CURRENTLY SELECTED FAMILY
    133129                                             isoPhas            !--- USED PHASES    FOR THE CURRENTLY SELECTED FAMILY
    134    INTEGER, TARGET,         SAVE          ::  niso, nzone, &    !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
    135                                              nphas, ntiso       !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
    136    INTEGER,                 SAVE, POINTER :: iZonIso(:,:)       !--- INDEX IN "isoTrac" AS f(tagging zone, isotope)
    137    INTEGER,                 SAVE, POINTER :: iTraPha(:,:)       !--- INDEX IN "isoTrac" AS f(isotopic tracer, phase)
    138    INTEGER, ALLOCATABLE,    SAVE ::  index_trac(:,:) ! numero ixt en fn izone, indnum entre 1 et niso
    139    INTEGER, ALLOCATABLE,    SAVE ::  iqiso(:,:)      ! donne indice iq en fn de (ixt,phase)
    140 
    141    !--- Aliases for older names
    142    INTEGER, POINTER, SAVE :: ntraciso, ntraceurs_zone
    143    REAL,             SAVE :: qperemin, masseqmin, ratiomin
    144 
    145 ! CRisi: cas particulier des isotopes
    146    INTEGER, PARAMETER :: niso_possibles = 5
    147    LOGICAL, SAVE      :: ok_isotopes, ok_iso_verif, ok_isotrac, ok_init_iso
     130   INTEGER,                 SAVE, POINTER ::  niso, nzone, &    !--- NUMBER OF ISOTOPES, TAGGING ZONES AND PHASES
     131                                             nphas, ntiso, &    !--- NUMBER OF PHASES AND ISOTOPES + ISOTOPIC TAGGING TRACERS
     132                                            itZonIso(:,:), &    !--- INDEX IN "isoTrac" AS f(tagging zone idx,  isotope idx)
     133                                            iqTraPha(:,:)       !--- INDEX IN "qx"      AS f(isotopic tracer idx, phase idx)
     134
     135   !--- Aliases for older names + quantities to be removed soon
     136   INTEGER,                 SAVE, POINTER ::  index_trac(:,:)   ! numero ixt en fn izone, indnum entre 1 et niso
     137   INTEGER,                 SAVE, POINTER ::  iqiso(:,:)        ! donne indice iq en fn de (ixt,phase)
     138   INTEGER,                 SAVE, POINTER :: ntraciso, ntraceurs_zone
     139   REAL,    SAVE :: qperemin, masseqmin, ratiomin
     140   INTEGER, SAVE :: niso_possibles
     141   LOGICAL, SAVE :: ok_isotopes, ok_iso_verif, ok_isotrac
    148142   LOGICAL, SAVE, ALLOCATABLE ::       use_iso(:)
    149    INTEGER, SAVE, ALLOCATABLE ::    iso_indnum(:)     !--- Gives 1<=idx<=niso_possibles as function(1<=iq <=nqtot)
    150    INTEGER, SAVE, ALLOCATABLE :: indnum_fn_num(:)     !--- Gives 1<=idx<=niso           as function(1<=idx<=niso_possibles)
     143   INTEGER, SAVE, ALLOCATABLE ::       iso_num(:)               !--- idx in [1,niso_possibles] = f(1<=iq <=nqtot)
     144   INTEGER, SAVE, ALLOCATABLE ::    iso_indnum(:)               !--- idx in [1,niso]           = f(1<=iq <=nqtot)
     145   INTEGER, SAVE, ALLOCATABLE :: indnum_fn_num(:)               !--- idx in [1,niso]           = f(1<=idx<=niso_possibles)
    151146
    152147   !=== VARIABLES FOR ISOTOPES INITIALIZATION AND FOR INCA
    153    REAL,               SAVE, ALLOCATABLE ::     tnat(:), &     !--- Natural relative abundance of water isotope        (niso)
    154                                          alpha_ideal(:)         !--- Ideal fractionning coefficient (for initial state) (niso)
    155    INTEGER,            SAVE, ALLOCATABLE :: conv_flg(:), &     !--- Convection     activation ; needed for INCA        (nbtr)
    156                                              pbl_flg(:)         !--- Boundary layer activation ; needed for INCA        (nbtr)
    157    CHARACTER(LEN=8),   SAVE, ALLOCATABLE ::   solsym(:)         !--- Names from INCA                                    (nbtr)
     148   REAL,                SAVE, ALLOCATABLE ::     tnat(:), &     !--- Natural relative abundance of water isotope        (niso)
     149                                          alpha_ideal(:)        !--- Ideal fractionning coefficient (for initial state) (niso)
     150   INTEGER,             SAVE, ALLOCATABLE :: conv_flg(:), &     !--- Convection     activation ; needed for INCA        (nbtr)
     151                                              pbl_flg(:)        !--- Boundary layer activation ; needed for INCA        (nbtr)
     152   CHARACTER(LEN=8),    SAVE, ALLOCATABLE ::   solsym(:)        !--- Names from INCA                                    (nbtr)
    158153   LOGICAL, PARAMETER :: lOldCode = .TRUE.
    159154
    … …  
    175170!   05/94: F.Forget      Modif special traceur
    176171!   02/02: M-A Filiberti Lecture de traceur.def
    177 !   01/22: D. Cugnet     Nouveaux tracer.def et tracer_*.def + encapsulation (types tr et iso)
     172!   01/22: D. Cugnet     Nouveaux tracer.def et tracer_*.def + encapsulation (types trac_type et isot_type)
    178173!
    179174!   Objet:
    … …  
    212207   TYPE(isot_type), POINTER             :: iso
    213208
    214    CHARACTER(LEN=maxlen), ALLOCATABLE :: tnom_0(:), tnom_transp(:)        !--- Tracer short name + transporting fluid name
     209   CHARACTER(LEN=maxlen), ALLOCATABLE :: tnom_0(:), tnom_transp(:)   !--- Tracer short name + transporting fluid name
    215210   CHARACTER(LEN=maxlen)              :: tchaine
    216211   INTEGER :: ierr
    217    LOGICAL :: lINCA
    218212
    219213   CHARACTER(LEN=*), PARAMETER :: modname="infotrac_init"
    … …  
    238232      !--- MESSAGE ABOUT THE CHOSEN CONFIGURATION
    239233      msg1 = 'For type_trac = "'//TRIM(str(it))//'":'
    240       SELECT CASE(type_trac)
     234      SELECT CASE(str(it))
    241235         CASE('inca'); CALL msg(TRIM(msg1)//' coupling with INCA chemistry model, config_inca='//config_inca, modname)
    242236         CASE('inco'); CALL msg(TRIM(msg1)//' coupling jointly with INCA and CO2 cycle',  modname)
    … …  
    254248      !--- COHERENCE TEST BETWEEN "type_trac" AND PREPROCESSING KEYS
    255249      SELECT CASE(str(it))
    256          CASE('inca','inco')
     250         CASE('inca', 'inco')
    257251#ifndef INCA
    258252            CALL abort_gcm(modname, 'You must add cpp key INCA and compile with INCA code', 1)
    … …  
    283277
    284278!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    285    IF(lOldCode) THEN
     279   IF(lOldCode) THEN         !--- "type_trac" is a single keyword => no need to loop on its parsed version "str(:)"
    286280!------------------------------------------------------------------------------------------------------------------------------
    287281   !--- Determine nqtrue and (INCA only) nqo, nbtr
    … …  
    289283   IF(ierr /= 0) CALL abort_gcm(modname, 'file "traceur.def" not found !', 1)
    290284   CALL msg('File "traceur.def" successfully opened.', modname)
    291    lINCA = ANY(['inca','inco'] == type_trac)
    292 
    293    IF(lINCA) THEN
     285
     286   IF(ANY(['inca','inco'] == type_trac)) THEN
    294287#ifdef INCA
    295288      READ(90,*) nqo
    … …  
    299292      nqtrue = nbtr + nqo
    300293      IF(ALL([2,3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo='//TRIM(int2str(nqo)), 1)
    301       CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
    302       CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
    303       CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
    304       CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
    305       CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
    306294      ALLOCATE(hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
    307295      ALLOCATE(vadv_inca(nqINCA),  pbl_flg_inca(nqINCA))
    308296      CALL init_transport(hadv_inca, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
    309       ! DC passive CO2 tracer is at position 1: H2O was removed ; nqCO2/=0 in "inco" case only
    310       ALLOCATE(conv_flg(nbtr),pbl_flg(nbtr),solsym(nbtr))
    311       conv_flg = [(  1,        ic=1, nqCO2),conv_flg_inca]
    312        pbl_flg = [(  1,        ic=1, nqCO2), pbl_flg_inca]
    313        solsym  = [('CO2     ', ic=1, nqCO2), solsym_inca]
    314       DEALLOCATE(conv_flg_inca, pbl_flg_inca)
    315297#endif
    316298   ELSE
    … …  
    337319      END IF
    338320#endif
    339       CALL msg('237: iq='//TRIM(int2str(iq)), modname)
    340321      READ(90,'(I2,X,I2,X,A)',IOSTAT=ierr) hadv(iq),vadv(iq),tchaine
    341       WRITE(msg1,'("hadv(",i0,"), vadv(",i0,") = ",i0,", ",i0)')iq, iq, hadv(iq), vadv(iq)
    342       CALL msg(TRIM(msg1), modname)
    343       CALL msg('tchaine = "'//TRIM(tchaine)//'"', modname)
    344       CALL msg('infotrac 238: IOstatus='//TRIM(int2str(ierr)), modname)
    345322      IF(ierr/=0) CALL abort_gcm('infotrac_init', 'Pb dans la lecture de traceur.def', 1)
    346323      jq = INDEX(tchaine(1:LEN_TRIM(tchaine)),' ')
    347324      CALL msg("Ancienne version de traceur.def: traceurs d'air uniquement", modname, iq==1 .AND. jq==0)
    348325      CALL msg("Nouvelle version de traceur.def",                            modname, iq==1 .AND. jq/=0)
     326      CALL msg('iq, hadv, vadv, tchaine ='//TRIM(strStack(int2str([iq, hadv(iq), vadv(iq)])))//', '//TRIM(tchaine), modname)
    349327      IF(jq /= 0) THEN                                               !--- Space in the string chain => new format
    350328         tnom_0     (iq) = tchaine(1:jq-1)
    … …  
    354332         tnom_transp(iq) = 'air'
    355333      END IF
    356       CALL msg(     'tnom_0(iq)=<'//TRIM(tnom_0(iq))     //'>', modname)
    357       CALL msg('tnom_transp(iq)=<'//TRIM(tnom_transp(iq))//'>', modname)
    358    END DO
    359 #ifdef INCA
    360    DEALLOCATE(solsym_inca)
    361 #endif
    362 
     334   END DO
    363335   CLOSE(90)
    364336
    365337#ifndef INCA
    366    CALL msg('Valeur de traceur.def :', modname)
    367    CALL msg('nombre total de traceurs '//TRIM(int2str(nqtrue)), modname)
    368    DO iq = 1, nqtrue
    369       CALL msg(strStack([int2str(hadv(iq)), int2str(vadv(iq)), tnom_0(iq), tnom_transp(iq)]), modname)
    370    END DO
    371338   IF(planet_type /= 'earth') nqo = 0                                !--- Same number of tracers in dynamics and physics
    372339   IF(planet_type == 'earth') nqo = COUNT(delPhase(tnom_0) == 'H2O') !--- for all planets except for Earth
    373    nbtr = nqtrue - nqo               
    374    ALLOCATE(conv_flg(nbtr),pbl_flg(nbtr),solsym(nbtr))
    375    conv_flg(1:nbtr) = 1                                     !--- Convection activated for all tracers
    376    pbl_flg(1:nbtr) = 1                                     !--- Boundary layer activated for all tracers
    377 #endif
     340   nbtr = nqtrue - nqo
     341#endif
     342
     343   CALL msg('RAW CONTENT OF "traceur.def" FILE:', modname)
     344   IF(dispTable('iiss', ['hadv  ', 'vadv  ', 'name  ', 'parent'], cat(tnom_0, tnom_transp), cat(hadv, vadv))) &
     345      CALL abort_gcm(modname, "problem with the tracers table content", 1)
    378346
    379347   !--- SET FIELDS %name, %parent, %phase, %component
    380    tracers(:)%name      = tnom_0
    381    tracers(:)%parent    = tnom_transp
    382    tracers(:)%phase     = 'g'
     348   tracers(:)%name      = old2newName(tnom_0)
     349   tracers(:)%parent    = old2newName(tnom_transp)
     350   tracers(:)%phase     = [( getPhase(tracers(iq)%name), iq=1, nqtrue )]
    383351   tracers(:)%component = type_trac
    384 
    385 
    386352   DO iq = 1, nqtrue
    387       ip = strIdx([(addPhase('H2O',old_phases(ix:ix),''), ix=1, nphases)], strHead(tracers(iq)%name,'_',.TRUE.))
    388       IF(ip == 0) CYCLE
    389       tracers(iq)%phase = known_phases(ip:ip)
    390       tracers(iq)%component = 'lmdz'
    391    END DO
    392    IF(lINCA) tracers(1+nqo:nqCO2+nqo)%component = 'co2i'
     353      IF(addPhase('H2O',tracers(iq)%phase) == tracers(iq)%name) tracers(iq)%component = 'lmdz'
     354   END DO
     355   IF(ANY(['inca','inco'] == type_trac)) tracers(1+nqo:nqCO2+nqo)%component = 'co2i'
    393356   CALL setGeneration(tracers)                                       !--- SET FIELDS %iGeneration, %gen0Name
    394 ! manque "type"
     357   WHERE(tracers(:)%iGeneration == 2) tracers(:)%type = 'tag'        !--- DEFAULT VALUE: "tracer"
     358
     359   !--- FINALIZE
     360   DEALLOCATE(tnom_0, tnom_transp)
     361#ifdef INCA
     362   DEALLOCATE(hadv_inca, vadv_inca, conv_flag_inca, pbl_flag_inca, solsym_inca)
     363#endif
    395364
    396365!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    … …  
    400369   IF(fType == 0) CALL abort_gcm(modname, 'Missing tracers file: "traceur.def", "tracer.def" or "tracer_<keyword>.def file.',1)
    401370   !---------------------------------------------------------------------------------------------------------------------------
    402    IF(fType == 1) THEN                                               !=== FOUND AN OLD STYLE "traceur.def"
     371   IF(fType == 1 .AND. ANY(['inca','inco'] == type_trac)) THEN       !=== FOUND OLD STYLE INCA "traceur.def" (single type_trac)
    403372   !---------------------------------------------------------------------------------------------------------------------------
    404373#ifdef INCA
    … …  
    409378      nqtrue = nbtr + nqo                                            !--- Total number of "true" tracers
    410379      IF(ALL([2,3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo='//TRIM(int2str(nqo)), 1)
    411       CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
    412       CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
    413       CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
    414       CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
    415       CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
    416       ALLOCATE(hadv(nqtrue), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
    417       ALLOCATE(vadv(nqtrue), vadv_inca(nqINCA),  pbl_flg_inca(nqINCA))
     380      ALLOCATE(hadv(nqtrue), conv_flg(nbtr), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
     381      ALLOCATE(vadv(nqtrue),  pbl_flg(nbtr), vadv_inca(nqINCA),  pbl_flg_inca(nqINCA), solsym(nbtr))
    418382      CALL init_transport(hadv_inca, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
    419       ! DC passive CO2 tracer is at position 1: H2O was removed ; nqCO2/=0 in "inco" case only
    420      
    421       conv_flg = [(  1       , k=1, nqCO2), conv_flg_inca]
    422       pbl_flg  = [(  1       , k=1, nqCO2), pbl_flg_inca]
    423       solsym   = [('CO2     ', k=1, nqCO2), solsym_inca]
    424       DEALLOCATE(conv_flg_inca, pbl_flg_inca, solsym_inca)
     383      !--- Passive CO2 tracer is at position 1 because: H2O has been removed ; nqCO2/=0 in "inco" case only
     384      conv_flg(1:nbtr) = [(1,          k=1, nqCO2), conv_flg_inca]
     385       pbl_flg(1:nbtr) = [(1,          k=1, nqCO2),  pbl_flg_inca]
     386       solsym (1:nbtr) = [('CO2     ', k=1, nqCO2),   solsym_inca]
    425387      ALLOCATE(ttr(nqtrue))
    426388      ttr(1:nqo+nqCO2)                    = tracers
    … …  
    433395      lerr = getKey('hadv', had, ky=ttr(:)%keys); hadv(:) = [had, hadv_inca]
    434396      lerr = getKey('vadv', vad, ky=ttr(:)%keys); vadv(:) = [vad, vadv_inca]
    435       DEALLOCATE(had, hadv_inca, vad, vadv_inca)
    436397      CALL MOVE_ALLOC(FROM=ttr, TO=tracers)
    437398      CALL setGeneration(tracers)                                    !--- SET FIELDS %iGeneration, %gen0Name
    438 #else
    439       nqo    = COUNT(delPhase(tracers(:)%name) == 'H2O')             !--- Number of water phases
     399      DEALLOCATE(had, hadv_inca, vad, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
     400#endif
     401   !---------------------------------------------------------------------------------------------------------------------------
     402   ELSE                                                              !=== FOUND NEW STYLE TRACERS CONFIGURATION FILE(S)
     403   !---------------------------------------------------------------------------------------------------------------------------
     404      nqo    =        COUNT(delPhase(tracers(:)%name)     == 'H2O' &
     405                               .AND. tracers(:)%component == 'lmdz') !--- Number of water phases
    440406      nqtrue = SIZE(tracers)                                         !--- Total number of "true" tracers
    441       nbtr   = nqtrue - nqo                                          !--- Number of tracers passed to phytrac
     407      nbtr   = nqtrue-COUNT(delPhase(tracers(:)%gen0Name) == 'H2O' &
     408                               .AND. tracers(:)%component == 'lmdz') !--- Number of tracers passed to phytrac
    442409      lerr = getKey('hadv', hadv, ky=tracers(:)%keys)
    443410      lerr = getKey('vadv', vadv, ky=tracers(:)%keys)
    444       ALLOCATE(solsym(nbtr))
    445       conv_flg(1:nbtr)=1  !--- Convection activated for all tracers
    446       pbl_flg(1:nbtr)=1   !--- Boundary layer activated for all tracers
    447 #endif
    448    !---------------------------------------------------------------------------------------------------------------------------
    449    ELSE                                                              !=== FOUND NEW STYLE TRACERS CONFIGURATION FILE(S)
    450    !---------------------------------------------------------------------------------------------------------------------------
    451       nqo    = COUNT(delPhase(tracers(:)%name) == 'H2O')             !--- Number of water phases
    452       nqtrue = SIZE(tracers)                                         !--- Total number of "true" tracers
    453       nbtr   = nqtrue - nqo                                          !--- Number of tracers passed to phytrac
    454       lerr = getKey('hadv', hadv, ky=tracers(:)%keys)
    455       lerr = getKey('vadv', vadv, ky=tracers(:)%keys)
    456       ALLOCATE(solsym(nbtr))
    457       conv_flg(1:nbtr)=1  !--- Convection activated for all tracers
    458        pbl_flg(1:nbtr)=1  !--- Boundary layer activated for all tracers
     411      ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
     412      conv_flg(1:nbtr) = [(1, it=1, nbtr)]                           !--- Convection activated for all tracers
     413       pbl_flg(1:nbtr) = [(1, it=1, nbtr)]                           !--- Boundary layer activated for all tracers
    459414   !---------------------------------------------------------------------------------------------------------------------------
    460415   END IF
    … …  
    488443      CALL msg('The total number of tracers needed is '//TRIM(int2str(nqtot)))
    489444   END IF
    490    CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
    491    CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
    492    CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
    493    CALL msg('nqtot  = '//TRIM(int2str(nqtot)),  modname)
    494445
    495446!==============================================================================================================================
    … …  
    527478      t1%iadv       = iad
    528479      t1%isAdvected = iad >= 0
    529       t1%isInPhysics= .not. (delPhase(t1%gen0Name) == 'H2O' .and. t1%component=='lmdz')  !=== TO BE COMPLETED WITH OTHER EXCEPTIONS: CO2i, SURSATURATED CLOUDS...
     480      t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' &
     481                          .OR. t1%component /= 'lmdz' !=== OTHER EXCEPTIONS TO BE ADDED: CO2i, SURSATURATED WATER CLOUD...
    530482      ttr(iq)       = t1
    531483
    532484      !--- DEFINE THE HIGHER ORDER TRACERS, IF ANY
    533485      nm = 0
    534       IF(iad == 20) nm = 3                                             !--- 2nd order scheme
    535       IF(iad == 30) nm = 9                                             !--- 3rd order scheme
    536       IF(nm == 0) CYCLE                                                !--- No higher moments
     486      IF(iad == 20) nm = 3                                           !--- 2nd order scheme
     487      IF(iad == 30) nm = 9                                           !--- 3rd order scheme
     488      IF(nm == 0) CYCLE                                              !--- No higher moments
    537489      ttr(jq+1:jq+nm)             = t1
    538490      ttr(jq+1:jq+nm)%name        = [(TRIM(t1%name)    //'-'//TRIM(suff(im)), im=1, nm) ]
    … …  
    549501   CALL indexUpdate(tracers)
    550502
    551    CALL msg('Information stored in infotrac :', modname)
    552    CALL msg('iadv  name  long_name :', modname)
    553 
    554503   !=== TEST ADVECTION SCHEME
    555504   DO iq=1,nqtot ; t1 => tracers(iq); iad = t1%iadv
    … …  
    568517
    569518      !--- ONLY VALID SCHEME NUMBER FOR WATER VAPOUR:            iadv = 14
    570       ll = t1%name /= addPhase('H2O','g'); IF(lOldCode) ll = t1%name /= 'H2Ov'
     519      ll = t1%name /= addPhase('H2O','g')
    571520      IF(fmsg('WARNING ! iadv=14 is valid for water vapour only. Setting iadv=10 for "'//TRIM(t1%name)//'".', &
    572521         modname, iad == 14 .AND. ll))                 t1%iadv = 10
    … …  
    578527   CALL infotrac_isoinit                    !--- SET FIELDS %type, %iso_iName, %iso_iZone, %iso_iPhase
    579528   CALL getKey_init(tracers, isotopes)
    580    IF(isoSelect('H2O')) RETURN                                    !--- Select water isotopes ; finished if no water isotopes
    581    iH2O = ixIso                                                   !--- Keep track of water family index
     529   IF(isoSelect('H2O')) RETURN              !--- Select water isotopes ; finished if no water isotopes
     530   iH2O = ixIso                             !--- Keep track of water family index
    582531
    583532   !--- Remove the isotopic tracers from the tracers list passed to phytrac
    584533   nbtr    = nbtr -nqo*   ntiso             !--- ISOTOPIC TAGGING TRACERS ARE NOT PASSED TO THE PHYSICS
    585534   nqtottr = nqtot-nqo*(1+ntiso)            !--- NO H2O-FAMILY    TRACER  IS      PASSED TO THE PHYSICS
    586    CALL msg('702: nbtr, ntiso='//strStack(int2str([nbtr, ntiso])), modname)
    587    CALL msg('704: nqtottr, nqtot, nqo = '//strStack(int2str([nqtottr, nqtot, nqo])), modname)
    588    ! Rq: nqtottr n'est pas forcement egal a nbtr dans le cas ou nmom/=0
    589    IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers(:)%name) == 'H2O' .AND. tracers(:)%component=='lmdz') /= nqtottr) &
    590       CALL abort_gcm('infotrac_init', 'pb dans le calcul de nqtottr', 1)
    591 
    592    !--- Finalize :
    593    DEALLOCATE(tnom_0, tnom_transp)
     535
     536   ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), solsym(nbtr))
     537#ifndef INCA
     538   conv_flg(1:nbtr) = 1                                              !--- Convection activated for all tracers
     539    pbl_flg(1:nbtr) = 1                                              !--- Boundary layer activated for all tracers
     540#else
     541   !--- Passive CO2 tracer is at position 1 because: H2O has been removed ; nqCO2/=0 in "inco" case only
     542   conv_flg(1:nbtr) = [(  1,        ic=1, nqCO2),conv_flg_inca]
     543    pbl_flg(1:nbtr) = [(  1,        ic=1, nqCO2), pbl_flg_inca]
     544     solsym(1:nbtr) = [('CO2     ', ic=1, nqCO2),  solsym_inca]
     545#endif
    594546
    595547ELSE
    596548
    597549   CALL initIsotopes(tracers, isotopes)
    598    nbIso = SIZE(isotopes); IF(nbIso==0) RETURN                    !--- No isotopes: finished.
    599 
    600    !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE SPECIFIC TO WATER ISOTOPES
    601    !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat, alpha_ideal)
    602    CALL getKey_init(tracers, isotopes)
    603    IF(isoSelect('H2O')) RETURN                                    !--- Select water isotopes ; finished if no water isotopes
    604    iH2O = ixIso                                                   !--- Keep track of water family index
    605    IF(getKey('tnat' , tnat,        isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "tnat"', 1)
    606    IF(getKey('alpha', alpha_ideal, isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "alpha_ideal"', 1)
    607 
    608    !=== ENSURE THE MEMBERS OF AN ISOTOPES FAMILY ARE PRESENT IN THE SAME PHASES
    609    DO ix = 1, nbIso
    610       iso => isotopes(ix)
    611       !--- Check whether each isotope and tagging isotopic tracer is present in the same number of phases
    612       DO it = 1, iso%ntiso
    613          np = SUM([(COUNT(tracers(:)%name == addPhase(iso%trac(it), iso%phase(ip:ip))), ip=1, iso%nphas)])
    614          IF(np == iso%nphas) CYCLE
    615          WRITE(msg1,'("Found ",i0," phases for ",s," instead of ",i0)')np, iso%trac(it), iso%nphas
    616          CALL abort_gcm(modname, msg1, 1)
     550   nbIso = SIZE(isotopes)
     551   nqtottr = nqtot - COUNT(tracers%gen0Name == 'H2O' .AND. tracers%component == 'lmdz')
     552   IF(nbIso/=0) THEN                        !--- ISOTOPES FOUND
     553
     554      !=== READ PHYSICAL PARAMETERS FROM "isotopes_params.def" FILE SPECIFIC TO WATER ISOTOPES
     555      !    DONE HERE, AND NOT ONLY IN "infotrac_phy", BECAUSE SOME PHYSICAL PARAMS ARE NEEDED FOR RESTARTS (tnat, alpha_ideal)
     556      CALL getKey_init(tracers, isotopes)
     557      IF(isoSelect('H2O')) RETURN           !--- Select water isotopes ; finished if no water isotopes
     558      iH2O = ixIso                          !--- Keep track of water family index
     559      IF(getKey('tnat' , tnat,        isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "tnat"', 1)
     560      IF(getKey('alpha', alpha_ideal, isoName(1:niso))) CALL abort_gcm(modname, 'can''t read "alpha_ideal"', 1)
     561
     562      !=== MAKE SURE THE MEMBERS OF AN ISOTOPES FAMILY ARE PRESENT IN THE SAME PHASES
     563      DO ix = 1, nbIso
     564         iso => isotopes(ix)
     565         !--- Check whether each isotope and tagging isotopic tracer is present in the same number of phases
     566         DO it = 1, iso%ntiso
     567            np = SUM([(COUNT(tracers(:)%name == addPhase(iso%trac(it), iso%phase(ip:ip))), ip=1, iso%nphas)])
     568            IF(np == iso%nphas) CYCLE
     569            WRITE(msg1,'("Found ",i0," phases for ",a," instead of ",i0)')np, TRIM(iso%trac(it)), iso%nphas
     570            CALL abort_gcm(modname, msg1, 1)
     571         END DO
     572         DO it = 1, iso%niso
     573            nz = SUM([(COUNT(iso%trac == TRIM(iso%trac(it))//'_'//iso%zone(iz)), iz=1, iso%nzone)])
     574            IF(nz == iso%nzone) CYCLE
     575            WRITE(msg1,'("Found ",i0," tagging zones for ",a," instead of ",i0)')nz, TRIM(iso%trac(it)), iso%nzone
     576            CALL abort_gcm(modname, msg1, 1)
     577         END DO
    617578      END DO
    618       DO it = 1, iso%niso
    619          nz = SUM([(COUNT(iso%trac == iso%trac(it)//'_'//iso%zone(iz)), iz=1, iso%nzone)])
    620          IF(nz == iso%nzone) CYCLE
    621          WRITE(msg1,'("Found ",i0," tagging zones for ",s," instead of ",i0)')nz, iso%trac(it), iso%nzone
    622          CALL abort_gcm(modname, msg1, 1)
    623       END DO
    624    END DO
    625    nqtottr = COUNT(tracers%iso_iName == 0)
     579   END IF
    626580
    627581END IF
    628582
    629    !=== DISPLAY THE RESULTING LIST
    630    t => tracers
    631    CALL msg('Information stored in infotrac :')
    632    IF(dispTable('isssssssssiiiiiiiii', &
    633       ['iq      ', 'name    ', 'longN.  ', 'gen0N.  ', 'parent  ', 'type    ', 'phase   ', 'compon. ', 'isAdv.  ', 'isPhy.  '&
    634       ,'iadv    ', 'iGen.   ', 'iqPar.  ', 'nqDes.  ', 'nqChil. ', 'iso_iG. ', 'iso_iN. ', 'iso_iZ. ', 'iso_iP. '],          &
    635       cat(t%name,  t%longName,  t%gen0Name,  t%parent,  t%type,  t%phase, &
    636           t%component, bool2str(t%isAdvected), bool2str(t%isInPhysics)),  &
    637       cat([(iq, iq=1, nqtot)],  t%iadv,  t%iGeneration, t%iqParent, t%nqDescen, &
    638          t%nqChilds, t%iso_iGroup, t%iso_iName, t%iso_iZone, t%iso_iPhase))) &
    639       CALL abort_gcm(modname, "problem with the tracers table content", 1)
     583   !--- Note: nqtottr can differ from nbtr when nmom/=0
     584!   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
     585!      CALL abort_gcm('infotrac_init', 'pb dans le calcul de nqtottr', 1)
    640586
    641587   !--- Some aliases to be removed later
    642    ntraciso       => isotope%ntiso
    643    ntraceurs_zone => isotope%nzone
     588   ntraciso       => ntiso
     589   ntraceurs_zone => nzone
    644590   qperemin       =  min_qParent
    645591   masseqmin      =  min_qMass
    646592   ratiomin       =  min_ratio
     593   iqiso          => iqTraPha
     594   index_trac     => itZonIso
     595
     596   !=== DISPLAY THE RESULTS
     597   CALL msg('nqo    = '//TRIM(int2str(nqo)),    modname)
     598   CALL msg('nbtr   = '//TRIM(int2str(nbtr)),   modname)
     599   CALL msg('nqtrue = '//TRIM(int2str(nqtrue)), modname)
     600   CALL msg('nqtot  = '//TRIM(int2str(nqtot)),  modname)
     601   CALL msg('niso   = '//TRIM(int2str(niso)),   modname)
     602   CALL msg('ntiso  = '//TRIM(int2str(ntiso)),  modname)
     603#ifdef INCA
     604   CALL msg('nqCO2  = '//TRIM(int2str(nqCO2)),  modname)
     605   CALL msg('nqINCA = '//TRIM(int2str(nqINCA)), modname)
     606#endif
     607   t => tracers
     608   CALL msg('Information stored in infotrac :', modname)
     609   IF(dispTable('isssssssssiiiiiiiii', &
     610      ['iq    ', 'name  ', 'lName ', 'gen0N ', 'parent', 'type  ', 'phase ', 'compon', 'isAdv ', 'isPhy ', &
     611       'iadv  ', 'iGen  ', 'iqPar ', 'nqDes ', 'nqChld', 'iGroup', 'iName ', 'iZone ', 'iPhase'],          &
     612      cat(t%name, t%longName, t%gen0Name, t%parent, t%type, t%phase, t%component, bool2str(t%isAdvected),  &
     613                                                                                  bool2str(t%isInPhysics)),&
     614      cat([(iq, iq=1, nqtot)], t%iadv, t%iGeneration, t%iqParent, t%nqDescen, t%nqChilds, t%iso_iGroup,    &
     615                                                    t%iso_iName, t%iso_iZone, t%iso_iPhase), sub=modname)) &
     616      CALL abort_gcm(modname, "problem with the tracers table content", 1)
     617   IF(niso > 0) THEN
     618      CALL msg('Where, for isotopes family "'//TRIM(isotope%parent)//'":', modname)
     619      CALL msg('  isoKeys = '//strStack(isoKeys%name), modname)
     620      CALL msg('  isoName = '//strStack(isoName),      modname)
     621      CALL msg('  isoZone = '//strStack(isoZone),      modname)
     622      CALL msg('  isoPhas = '//TRIM(isoPhas),          modname)
     623   ELSE
     624      CALL msg('No isotopes identified.', modname)
     625   END IF
    647626   CALL msg('end', modname)
    648627
    … …  
    654633   !--- Purpose: Set fields %iqParent, %nqChilds, %iGeneration, %iqDescen, %nqDescen (old method)
    655634   USE strings_mod, ONLY: strIdx
    656    INTEGER               :: iq, ipere, ifils
     635   INTEGER               :: iq, jq, ipere, ifils
    657636   INTEGER, ALLOCATABLE  :: iqfils(:,:)
    658637   CHARACTER(LEN=maxlen) :: msg1, modname='infotrac_init'
    … …  
    661640   !=== SET FIELDS %iqParent, %nqChilds
    662641   ALLOCATE(iqfils(nqtot,nqtot)); iqfils(:,:) = 0
     642   tracers(:)%nqChilds = 0
     643   tracers(:)%iqParent = 0
    663644
    664645   DO iq = 1, nqtot
    … …  
    684665   CALL msg('iqChilds = '//strStack(int2str(PACK(iqfils,MASK=.TRUE.))),modname)
    685666
    686 
    687667   !=== SET FIELDS %iGeneration, %iqDescen, %nqDescen
    688668   tracers(:)%iGeneration = 0
     669   tracers(:)%nqDescen = 0
    689670   DO iq = 1, nqtot
    690671      ifils = iq
    … …  
    703684   END DO
    704685
     686   !=== SET FIELD %gen0Name
     687   DO iq = 1, nqtot
     688      jq=iq; DO WHILE(tracers(jq)%iGeneration > 0); jq=tracers(jq)%iqParent; END DO
     689      tracers(iq)%gen0Name = tracers(jq)%name
     690   END DO
     691
    705692   CALL msg('nqDescen = '//TRIM(strStack(int2str(tracers(:)%nqDescen))), modname)
    706693   CALL msg('nqDescen_tot = ' //TRIM(int2str(SUM(tracers(:)%nqDescen))), modname)
    707    CALL msg('iqChilds = '//strStack(int2str(PACK(iqfils, MASK=.TRUE.))), modname)
    708 
     694   CALL msg('iqDescen = '//strStack(int2str(PACK(iqfils, MASK=.TRUE.))), modname)
    709695
    710696END SUBROUTINE infotrac_setHeredity
    … …  
    719705   USE ioipsl_getincom
    720706#endif
     707   USE readTracFiles_mod, ONLY: tnom_iso => newH2OIso
    721708   IMPLICIT NONE
    722    CHARACTER(LEN=3)      :: tnom_iso(niso_possibles)
    723    INTEGER, ALLOCATABLE  :: nb_iso(:,:), nb_traciso(:,:), nb_isoind(:)
    724    INTEGER               :: ii, ip, iq, it, iz, ixt, nzone_prec
     709   INTEGER, ALLOCATABLE  :: nb_iso(:), nb_tiso(:), nb_zone(:), ix(:)
     710   INTEGER               :: ii, ip, iq, it, iz, ixt
    725711   TYPE(isot_type), POINTER :: i
    726    TYPE(trac_type), POINTER :: t(:)
    727    CHARACTER(LEN=maxlen)    :: tnom_trac
     712   TYPE(trac_type), POINTER :: t(:), t1
     713   CHARACTER(LEN=maxlen)    :: tnom_trac, modname, t0
    728714   CHARACTER(LEN=maxlen), ALLOCATABLE :: str(:)
    729715   LOGICAL, DIMENSION(:), ALLOCATABLE :: mask
     716   REAL, ALLOCATABLE :: tnat0(:), alpha_ideal0(:)
    730717   INCLUDE "iniprint.h"
    731718
    732    tnom_iso = ['eau', 'HDO', 'O18', 'O17', 'HTO']
    733    ALLOCATE(nb_iso       (niso_possibles,nqo))
    734    ALLOCATE(nb_traciso   (niso_possibles,nqo))
    735    ALLOCATE(use_iso      (niso_possibles))
    736    ALLOCATE(indnum_fn_num(niso_possibles))
    737    ALLOCATE(iso_indnum(nqtot))
    738    ALLOCATE(nb_isoind(nqo))
    739 
    740    iso_indnum   (:) = 0
    741    use_iso      (:) = .FALSE.
    742    indnum_fn_num(:) = 0
    743    nb_iso     (:,:) = 0 
    744    nb_traciso (:,:) = 0
    745    nb_isoind    (:) = 0
    746 
    747    DO iq=1, nqtot
    748       IF(delPhase(tracers(iq)%name) == 'H2O' .OR. .NOT.tracers(iq)%isAdvected) CYCLE
    749 outer:DO ip = 1, nqo
    750          DO ixt= 1,niso_possibles
    751             tnom_trac = 'H2O'//old_phases(ip:ip)//'_'//TRIM(tnom_iso(ixt))
    752             IF (tracers(iq)%name == tnom_trac) THEN
    753                nb_iso(ixt,ip)         = nb_iso(ixt,ip)+1
    754                nb_isoind (ip)         = nb_isoind (ip)+1
    755                tracers(iq)%type       = 'tracer'
    756                tracers(iq)%iso_iGroup = 1
    757                tracers(iq)%iso_iName  = ixt
    758                iso_indnum(iq)         = nb_isoind(ip)
    759                indnum_fn_num(ixt)     = iso_indnum(iq)
    760                tracers(iq)%iso_iPhase = ip
    761                EXIT outer
    762             ELSE IF(tracers(iq)%iqParent> 0) THEN
    763                IF(tracers(tracers(iq)%iqParent)%name == tnom_trac) THEN
    764                   nb_traciso(ixt,ip)     = nb_traciso(ixt,ip)+1
    765                   iso_indnum(iq)         = indnum_fn_num(ixt)
    766                   tracers(iq)%type       = 'tag'
    767                   tracers(iq)%iso_iGroup = 1
    768                   tracers(iq)%iso_iName  = ixt
    769                   tracers(iq)%iso_iZone  = nb_traciso(ixt,ip)
    770                   tracers(iq)%iso_iPhase = ip
    771                   EXIT outer
    772                END IF
    773             END IF
    774          END DO
    775       END DO outer
    776    END DO
    777 
    778    niso = 0; nzone_prec = nb_traciso(1,1)
    779    DO ixt = 1, niso_possibles
    780       IF(nb_iso(ixt,1) == 0) CYCLE
    781       IF(nb_iso(ixt,1) /= 1) CALL abort_gcm('infotrac_init', 'Isotopes are not well defined in traceur.def', 1)
    782 
    783       ! on verifie que toutes les phases ont le meme nombre d'isotopes
    784       IF(ANY(nb_iso(ixt,:) /= 1)) CALL abort_gcm('infotrac_init', 'Phases must have same number of isotopes', 1)
    785 
    786       niso = niso+1
    787       use_iso(ixt) = .TRUE.
    788       nzone = nb_traciso(ixt,1)
    789 
    790       ! on verifie que toutes les phases ont le meme nombre de traceurs d'isotopes
    791       IF(ANY(nb_traciso(ixt,2:nqo) /= nzone)) CALL abort_gcm('infotrac_init','Phases must have same number of tracers',1)
    792 
    793       ! on verifie que tous les isotopes ont le meme nombre de traceurs d'isotopes
    794       IF(nzone /= nzone_prec) CALL abort_gcm('infotrac_init','Isotope tracers are not well defined in traceur.def',1)
    795       nzone_prec = nzone
    796    END DO
    797 
    798    ! dimensions et flags isotopiques:
    799    ntiso = niso*(nzone+1)
    800    ok_isotopes = niso  > 0
    801    ok_isotrac  = nzone > 0
    802  
    803    IF(ok_isotopes) THEN
    804       ok_iso_verif = .FALSE.; CALL getin('ok_iso_verif', ok_iso_verif)
    805       ok_init_iso  = .FALSE.; CALL getin('ok_init_iso',  ok_init_iso)
    806    END IF
    807       tnat        = [1.0, 155.76e-6, 2005.2e-6, 0.004/100., 0.0]
    808       alpha_ideal = [1.0, 1.01,      1.006,     1.003,      1.0]
    809 !   END IF
    810 
    811    ! remplissage du tableau iqiso(ntiso,phase)
    812    ALLOCATE(iqiso(ntiso,nqo))   
    813    iqiso(:,:)=0     
    814    DO iq = 1, nqtot
    815       IF(tracers(iq)%iso_iName <= 0) CYCLE
    816       ixt = iso_indnum(iq) + tracers(iq)%iso_iZone*niso
    817       iqiso(ixt, tracers(iq)%iso_iPhase) = iq
    818    END DO
    819 
    820    ! remplissage du tableau index_trac(nzone,niso)
    821    ALLOCATE(index_trac(nzone, niso)) 
    822    IF(ok_isotrac) then
    823       DO ii = 1, niso; index_trac(:, ii) = ii + niso*[(iz, iz=1, nzone)]; END DO
    824    ELSE
    825       index_trac(:,:)=0.0
    826    END IF
    827 
     719   modname = 'infotrac_isoinit'
     720   tnat0       = [ 1.0 ,  155.76e-6, 2005.2e-6, 0.004/100., 0.0 ]    !--- Same length as tnom_iso
     721   alpha_ideal0= [ 1.0 ,  1.01,      1.006,     1.003,      1.0 ]    !--- Same length as tnom_iso
    828722   ALLOCATE(isotopes(1))                                             !--- Only water
    829    nbIso = 1
     723   nbIso = SIZE(isotopes)
    830724   t => tracers
    831725   i => isotopes(1)
    832726   i%parent = 'H2O'
    833727
    834    !--- Isotopes names list (embedded in the "keys" field)
    835    i%niso  = niso
    836    ALLOCATE(i%keys(i%niso))
    837    mask = t%type=='tracer' .AND. delPhase(t%gen0Name)=='H2O' .AND. t%phase == 'g' .AND. t%iGeneration==1
    838    str = strTail(PACK(delPhase(t%name), MASK=mask), '_')
    839    CALL strReduce(str)
    840    i%keys(:)%name = str
    841 
    842    !--- Full isotopes list, with isotopes tagging tracers (if any) following the previous list
    843    i%ntiso = ntiso
    844    ALLOCATE(i%trac(i%ntiso))
    845    mask = t%type=='tag'    .AND. delPhase(t%gen0Name)=='H2O' .AND. t%phase == 'g' .AND. t%iGeneration==2
     728   !--- Effective isotopes names list (embedded in the "keys" field)
     729   mask = t%type=='tracer' .AND. t%gen0Name==addPhase('H2O', 'g') .AND. t%iGeneration==1
    846730   str = PACK(delPhase(t%name), MASK=mask)
    847731   CALL strReduce(str)
     732   i%niso = SIZE(str)
     733   ALLOCATE(i%keys(i%niso))
     734   i%keys(:)%name = str
     735
     736   !--- Check whether found isotopes are known
     737   mask = [(ALL(tnom_iso /= str(ii)), ii=1, i%niso)]
     738   IF(ANY(mask)) CALL abort_gcm(modname, 'The following isotopes are unknown: '//strStack(PACK(str, MASK=mask)), 1)
     739
     740   !--- Full isotopes list, with isotopes tagging tracers (if any) following the previous list
     741   mask = t%type=='tag'    .AND. t%gen0Name==addPhase('H2O', 'g') .AND. t%iGeneration==2
     742   str = PACK(delPhase(t%name), MASK=mask)
     743   i%ntiso = i%niso + SIZE(str)
     744   ALLOCATE(i%trac(i%ntiso))
    848745   i%trac(:) = [i%keys(:)%name, str]
    849746
    850    !--- Tagging zones names list
    851    i%nzone = nzone
     747   !--- Effective tagging zones names list
    852748   i%zone = strTail(str, '_', .TRUE.)
     749   CALL strReduce(i%zone)
     750   i%nzone = SIZE(i%zone)
     751   IF(i%ntiso /= i%niso*(i%nzone+1)) CALL abort_gcm(modname, 'Error in "ntiso" calculation', 1)
    853752
    854753   !--- Effective phases list
    855    i%nphas = nqo
    856754   i%phase = ''
    857    DO ip=1,nphases; IF(strIdx(t%name, addPhase('H2O',old_phases(ip:ip),''))/=0) i%phase=TRIM(i%phase)//known_phases(ip:ip); END DO
    858 
    859    !--- Table: index in "qx" of an isotope, knowing its indices "it","ip" in "isotope%iName,%iPhase"
    860    i%iTraPha = RESHAPE([((strIdx(t%name, TRIM(addPhase('H2O', new2oldPhase(i%phase(ip:ip)), ''))//'_'//TRIM(i%trac(it))), &
    861                           it=1,i%ntiso), ip=1,i%nphas)], [i%ntiso,i%nphas])
     755   DO ip=1,nphases; IF(strIdx(t%name, addPhase('H2O', ip))/=0) i%phase=TRIM(i%phase)//known_phases(ip:ip); END DO
     756   i%nphas = LEN_TRIM(i%phase)
     757
     758   !--- Indexes related to isotopes
     759   DO iq = 1, nqtot
     760      t1 => tracers(iq)
     761      t0 = t1%gen0Name
     762      IF(t1%iGeneration==0 .OR. .NOT.t1%isAdvected .OR. delPhase(t0)/='H2O') CYCLE
     763      t1%iso_iGroup = 1
     764      t1%iso_iPhase = INDEX(i%phase, getPhase(t0))
     765      t1%iso_iZone = strIdx(i%zone,  strTail(t1%name, '_'))
     766      IF(t1%iso_iZone /= 0) t1%iso_iName = strIdx(i%keys(:)%name, delPhase(t1%parent))
     767      IF(t1%iso_iZone == 0) t1%iso_iName = strIdx(i%keys(:)%name, delPhase(t1%name  ))
     768   END DO
     769
     770   niso_possibles = SIZE(tnom_iso)
     771!   ix = strIdx(tnom_iso, i%trac)
     772!   tnat        = tnat0       (PACK(ix, MASK=ix/=0))
     773!   alpha_ideal = alpha_ideal0(PACK(ix, MASK=ix/=0))
     774   tnat        = tnat0
     775   alpha_ideal = alpha_ideal0
     776
     777   !--- Tests
     778   nb_iso  = [(COUNT(t%iso_iPhase == ip .AND. t%iGeneration == 1), ip=1, i%nphas)]
     779   nb_tiso = [(COUNT(t%iso_iPhase == ip .AND. t%iGeneration == 2), ip=1, i%nphas)]
     780   nb_zone = [(COUNT(t%iso_iZone  == iz), iz=1, i%nzone)]
     781   IF(ANY(nb_iso (:) /= nb_iso (1))) CALL abort_gcm(modname, 'Phases must have same number of isotopes', 1)
     782   IF(ANY(nb_tiso(:) /= nb_tiso(1))) CALL abort_gcm(modname, 'Phases must have same number of tagging tracers', 1)
     783   IF(ANY(nb_zone(:) /= nb_zone(1))) CALL abort_gcm(modname, 'Isotopes must have the same number of tagging tracers', 1)
     784
     785   !--- Isotopic checking routines activation flag
     786   i%check = .FALSE.; IF(i%niso > 0) CALL getin('ok_iso_verif', i%check)
     787
     788   !--- Table: index in "qx(:)" of an isotope, knowing its indices "it","ip" in "isotope%iName,%iPhase"
     789   i%iqTraPha = RESHAPE([((strIdx(t%name, TRIM(addPhase(i%trac(it),ip,i%phase))),it=1,i%ntiso),ip=1,i%nphas)],[i%ntiso,i%nphas])
    862790
    863791   !--- Table: index in "isotope%tracs(:)%name" of an isotopic tagging tracer, knowing its indices "iz","ip" in "isotope%iZone,%iName"
    864    i%iZonIso = RESHAPE([((strIdx(i%trac,TRIM(i%trac(it))//'_'//TRIM(i%zone(iz))),iz=1,i%nzone),it=1,i%niso )],[i%nzone,i%niso ])
    865    DO it=1,ntiso
    866       WRITE(lunout,'(a,i0,a)')TRIM('infotrac_init')//': iqiso  (',it,',:) = '//strStack(int2str(iqiso(it,:)))
    867       WRITE(lunout,'(a,i0,a)')TRIM('infotrac_init')//': iTraPha(',it,',:) = '//strStack(int2str(i%iTraPha(it,:)))
    868    END DO
    869    DO iz=1,nzone
    870       WRITE(lunout,'(a,i0,a)')TRIM('infotrac_init')//': index_trac(',iz,',:) = '//strStack(int2str(index_trac(iz,:)))
    871       WRITE(lunout,'(a,i0,a)')TRIM('infotrac_init')//': iZonIso   (',iz,',:) = '//strStack(int2str(i%iZonIso(iz,:)))
    872    END DO
    873 
    874    ! Finalize :
     792   i%itZonIso = RESHAPE([((strIdx(i%trac,TRIM(i%trac(it))//'_'//TRIM(i%zone(iz))),iz=1,i%nzone),it=1,i%niso )],[i%nzone,i%niso])
     793
     794   DO it=1,i%ntiso; CALL msg('iqTraPha('//TRIM(int2str(it))//',:) = '//strStack(int2str(i%iqTraPha(it,:))), modname); END DO
     795   DO iz=1,i%nzone; CALL msg('itZonIso('//TRIM(int2str(iz))//',:) = '//strStack(int2str(i%itZonIso(iz,:))), modname); END DO
     796
     797   !--- Isotopic quantities (to be removed soon)
     798   ok_isotopes   = i%niso  > 0
     799   ok_isotrac    = i%nzone > 0
     800   ok_iso_verif  = i%check
     801   niso_possibles = SIZE(tnom_iso)
     802   iso_num       = [(strIdx(tnom_iso(:),    strHead(delPhase(tracers(iq)%name), '_')), iq=1, nqtot)]
     803   iso_indnum    = [(strIdx(i%keys(:)%name, strHead(delPhase(tracers(iq)%name), '_')), iq=1, nqtot)]
     804   indnum_fn_num = [(strIdx(i%keys(:)%name, tnom_iso(ixt)), ixt=1, niso_possibles)]
     805   use_iso       = indnum_fn_num /= 0            !--- .TRUE. for the effectively used isotopes of the possible isotopes list
     806
     807   !--- Finalize :
    875808   DEALLOCATE(nb_iso)
    876809
    … …  
    903836   lv = .FALSE.; IF(PRESENT(lVerbose)) lv = lVerbose
    904837   lerr = .FALSE.
    905    IF(iIso == ixIso) RETURN                                      !--- Nothing to do if the index is already OK
     838   IF(iIso == ixIso) RETURN                                          !--- Nothing to do if the index is already OK
    906839   lerr = iIso<=0 .OR. iIso>nbIso
    907840   CALL msg('Inconsistent isotopes family index '//TRIM(int2str(iIso))//': should be > 0 and <= '//TRIM(int2str(nbIso))//'"',&
    908841            ll=lerr .AND. lV)
    909842   IF(lerr) RETURN
    910    ixIso = iIso                                                  !--- Update currently selected family index
    911    isotope => isotopes(ixIso)                                    !--- Select corresponding component
    912    isoKeys => isotope%keys;    niso     = isotope%niso
    913    isoName => isotope%trac;    ntiso    = isotope%ntiso
    914    isoZone => isotope%zone;    nzone    = isotope%nzone
    915    isoPhas => isotope%phase;   nphas    = isotope%nphas
    916    iZonIso => isotope%iZonIso; isoCheck = isotope%check
    917    iTraPha => isotope%iTraPha
     843   ixIso = iIso                                                      !--- Update currently selected family index
     844   isotope  => isotopes(ixIso)                                       !--- Select corresponding component
     845   isoKeys  => isotope%keys;     niso     => isotope%niso
     846   isoName  => isotope%trac;     ntiso    => isotope%ntiso
     847   isoZone  => isotope%zone;     nzone    => isotope%nzone
     848   isoPhas  => isotope%phase;    nphas    => isotope%nphas
     849   itZonIso => isotope%itZonIso; isoCheck => isotope%check
     850   iqTraPha => isotope%iqTraPha
    918851END FUNCTION isoSelectByIndex
    919852!==============================================================================================================================
Note: See TracChangeset for help on using the changeset viewer.